clear all
clc
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
%S=S-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=10000;%求a次隨機情況
t=0;
j=2;
sum2(1)=0;
while(a~=0)
    B=rand(m,n-1);% 隨機矩陣B
    [m,g]=size(B);
    %r2=rank(B);
    B1=orth(B);
    [R1,jb1]=rref(B1);
    B1=B(:,jb1);
    P2=B1*B1';%空間B的有效投影矩陣    
    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(j)=sum2(j);%A中每個向量xi在A1上的正交投影之和    
    if(sum<sum2(j))
            t=t+1;
    end
    sum2(j)=0;
    j=j+1;
    a=a-1;        
end
t%若a=t，說明a次隨機情況中sum是最小的
sum%A中每個向量xi在A1上的正交投影之和
u=2:j-2;
w=min(x(u))%A中每個向量xi在B上的正交投影之和的最小值