目录
1 计算原理
1.1求解V
1.2求解D
1.3求解U
2 MATLAB程序
2.1注意
1 计算原理
设矩阵A的大小m*n（m>n），A = UD
1.1求解V
首先求出的特征值及特征值,对应的正交单位特征向量。将的特征值从大到小按顺序排列，对应的正交特征向量也要按特征值的顺序排列。排列好的正交的特征向量即为V = [V1,V2]，其中V1为非零特征值的特征列向量组成的矩阵，V2为零值的特征列向量组成矩阵。V的大小为n*n。
1.2求解D
D的形式为：
大小为m*n。其中为一r*r的矩阵，其主对角线的元素为的奇异值，对角线的元素按从大到小排列。即=diag {}。
1.3求解U
U1 =，它是一个m*r的形的次酉矩阵。U2为的零特征值对应的正交单位列向量组成的矩阵。U = [U1,U2],是一个m*m的矩阵。
2 MATLAB程序
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)>0n=n+1; endendsi = size(V);V1 = zeros(si(1),n);% for i=1:n% V1(:,i) = V(:,i); %非零特征值向量% endV1(:,1:n) = V(:,1:n); %非零特征值向量si = size(D); %分离零和非零特征向量if n==si(1)%无零特征值V = V1;else V2 = V(:,n+1:si(1));V = [V1,V2];endh = zeros(n,n);%特征值对角矩阵for i=1:nh(i,i) = D(si(1)-i+1,si(1)-i+1)^0.5; %奇异值矩阵endh = 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);endm(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)==0n = n+1;endendsi = 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)-nD(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
2.1注意
使用MATLAB计算矩阵奇异值分解时，会遇到计算过程中MATLAB的输出为零，但是查看变量时却不为零，是为一个远小于零的数。这个数极小，但是在变量中MATLAB却不认为等于零，会影响判断非零特征值和零特征值的个数，因此在程序中特定一个限定值，将MATLAB的特征值计算结果与限定值比较，小于限定值则认为特征值为零。产生这种结果的原因应该是浮点运算的精度问题
