奇异值2.m

SVD的应用

m=input('please iuput m number:')
n=input('please iuput n number:')
A=rand(m,n)% 随机矩阵A
r=rank(A);%求秩
[U,S,V]=svd(A);%SVD分解
[m,n]=size(A);%矩阵规模
S1=S(:,1:r-1);
S2=zeros(m,n-r+1);
S1=[S1 S2];%令原矩阵S最后一个奇异值为0的新矩阵S1
A1=U*S1*V';%由S1，U,V组成的新矩阵A1
[U,S,V]=svd(A1);
r1=rank(A1);
A2=orth(A1);%将矩阵A1正交规范化，A2的列与A1的列具有相同的空间
            %A2的列向量是正交向量,且满足：A2'*A2=eye(rank(A1))。
[R,jb]=rref(A2);%%求A1的最大线性无关组。jb是一个向量，其含义为：r = length(jb)为A的秩；
                %A(:, jb)为A的列向量基；jb中元素表示基向量所在的列
U1=A2(:,jb);%U1是U的一个最大线性无关组，因为U=(U1,U2)
P1=U1*U1';%空间A1的有效投影矩阵
sum=0;
for i=1:n
    b(:,i)=A(:,i);
    a2=(eye(m)-P1)*b(:,i);%向量xi在A1上的正交投影
    sum=norm(a2)+sum;
end 
sum%A中每个向量xi在A1上的正交投影之和

a=100;%求1000次随机情况
t=0;
j=2;
while(a~=0)
    B=rand(m,n);% 随机矩阵B
    [m,n]=size(B);
    B1=orth(B);
    [R1,jb1]=rref(B1);
    B1=B(:,jb1);
    P2=B1*B1';%空间B的有效投影矩阵
    sum2(1)=0;
    for i=1:n
        b(:,i)=A(:,i);
        a2=(eye(m)-P2)*b(:,i);%向量xi在A1上的正交投影矩阵
        sum2(j)=norm(a2)+sum2(j-1)       
    end 
    x=sum2(j)%A中每个向量xi在A1上的正交投影之和    
    if(sum<sum2(j))
            t=t+1;
    end
    j=j+1;
    a=a-1;    
end
t
sum
u=2:j-2;
m=min(sum2(u))