目录
设矩阵A的大小m*n（m>n），A = UD
首先求出的特征值及特征值
,对应的正交单位特征向量
。将
的特征值从大到小按顺序排列，对应的正交特征向量也要按特征值的顺序排列。排列好的正交的特征向量即为V = [V1,V2]，其中V1为非零特征值的特征列向量组成的矩阵，V2为零值的特征列向量组成矩阵。V的大小为n*n。
D的形式为：
大小为m*n。其中
为一r*r的矩阵，其主对角线的元素为
的奇异值，对角线的元素按从大到小排列。即
=diag {
}。
U1 = ，它是一个m*r的形的次酉矩阵。U2为
的零特征值对应的正交单位列向量组成的矩阵。U = [U1,U2],是一个m*m的矩阵。
function [U,D,V]=MySVD(A)
[m,n] = size(A);
B=A;
if m<n %仅适用于m>n的情况 如果m<n 现将矩阵转置再分解
A = A';
end
[v,D] = eig(A'*A);%D特征值
V = flip(v,2);
si = size(D);
n=0;
D(D<1e-10)=0;%很关键的一个限制，将小于某一个数的特征值设为0 否则无法判断非零特征值的数量
for i=1:si(1)
if D(i,i)>0
n=n+1;
end
end
si = size(V);
V1 = zeros(si(1),n);
% for i=1:n
% V1(:,i) = V(:,i); %非零特征值向量
% end
V1(:,1:n) = V(:,1:n); %非零特征值向量
si = size(D); %分离零和非零特征向量
if n==si(1)%无零特征值
V = V1;
else
V2 = V(:,n+1:si(1));
V = [V1,V2];
end
h = zeros(n,n);%特征值对角矩阵
for i=1:n
h(i,i) = D(si(1)-i+1,si(1)-i+1)^0.5; %奇异值矩阵
end
h = inv(h);
U1 = A*V1*h;
[v1,D1] = eig(A*A');%D特征值 V特征向量
si = size(D1);
m = zeros(si(1),1);
for i=1:si(1)
m(i) = D1(i,i);
end
m(m<1e-10)=0; %很关键的一个限制，将小于某一个数的特征值设为0，否则无法判断非零特征值的数量
%求U2
% m = eig(A*A'); %求A*A'特征值个数
si = size(m);
m = flip(m); %从大到小排列
n=0;%奇异值个数
for i=1:si(1)%计算零特征向量
if m(i)==0
n = n+1;
end
end
si = size(v1);
U2 = zeros(si(1),n);
U2 = v1(:,1:n);
U = [U1,U2];
D = zeros(size(A));
si = size(m);
for i=1:si(1)-n
D(i,i) = m(i)^0.5;
end
[m,n] = size(B);
if m<n %A = UDV' A' = VD'U'
u = zeros(size(U));
u = U;
U = V;
V = u;
D = D';
end使用MATLAB计算矩阵奇异值分解时，会遇到计算过程中MATLAB 的输出为零，但是查看变量时却不为零，是为一个远小于零的数。这个数极小，但是在变量中MATLAB却不认为等于零，会影响判断非零特征值和零特征值的个数，因此在程序中特定一个限定值，将MATLAB的特征值计算结果与限定值比较，小于限定值则认为特征值为零。产生这种结果的原因应该是浮点运算的精度问题