end B=B(:,1:n); B=B.'; V=qr(B); V1=V(:,1:r); U(:,1:r)=A*V1*(inv(D(1:r,1:r))); U(:,r+1:m)=null(A*A'); end
2.5 运行与数据分析
以教材上的 A=[1 0;0 1;1 0]为例来验证上述求矩阵的奇异值分解 程序的正确性。在 matlab 运行结果如下： >> A=[1 0;0 1;1 0]; >> [U1,D1,V1] = SVDecom(A) U1 = 0.7071 0 0.7071 D1 = 1.4142 0 0 V1 = 1 0 0 1 0 1.0000 0 0 1.0000 0 0.7071 0 -0.7071
s11 1 1 即有 U1=AV1 .其中 =
s2 1
sr 1
第四步： 解方程组 AAHy = 0， 对基础解系单位正交化可以求得 γr+1， γr+2，…，γm，令 U =（γ1 , γ2 , … , γr , γr+1 , γr+2 , … , γm）.
2 矩阵的奇异值分解
2.1 原理
设 A∈Cm×n，s1，s2，…，sr 是 A 的非零奇异值，则存在 m 阶酉矩 阵 U∈Cm×n 及 n 阶酉矩阵 V，m× n 矩阵 D，
s1 0 0 0 0 D= 0 0 sr 0 0 0 0 0 0 = 0 0 0 0
使得 A=UDVH 这就是矩阵 A 的奇异值分解.
2.2 算法
第一步：求出 AHA 的特征值 1 ≥ 2 ≥…≥ r ＞0= r 1 =…= n ，确定非 零奇异值 si = i ，i=1,2 …, r. 第二步：分别求出矩阵 AHA 的对应于特征值 i 的特征向量并将其 单位正交化，得到标准正交向量组 α1 , α2 , … , αn 令 V=（α1 , α2 , … , αn）=（V1 , V2） ，V1=（α1 , α2 , … , αr） ，V2= （αr+1 ,αr+2 , …, αn） 第三步：若 U=（γ1 , γ2 , … , γr , γr+1 , γr+2 , … , γm）=（U1 , U2） ，其 中 U1=（γ1 , γ2 , … , γr） ， U2=（γr+1 , γr+2 , … , γm） ， 则因（Aα1 , Aα2 , … , Aαr）=（s1γ1 , s2γ2 , … , srγratlab 自带求解矩阵奇异值分解函数： [U，S，V] = svd(A)其 中 U 就是所求的 U 矩阵，S 是所求的对角阵，V 就是所求的酉矩阵
V. 在 matlab 中运行下述指令： >> A=[1 0;0 1;1 0]; >> [U0,D0,V0] = svd(A) U0 = -0.7071 0 -0.7071 D0 = 1.4142 0 0 V0 = -1 0 0 1 0 1.0000 0 0 1.0000 0 -0.7071 0 0.7071
2.3 程序流程图
输入矩阵A
AHA的特征值及对应特征 向量
AHA的特征值由大到小排 列并排列及对应特征向量
单位化得V
计算U1 ，U2
计算D
结束
2.4 MATLAB 程序
function [U,D,V] = SVDecom(A) [m,n]=size(A); U=zeros(m);V=zeros(n); r=rank(A); D=zeros(m,n); [B,C]=eig(A'*A); x=diag(C); B=[B.',x]; B=sortrows(B,-(n+1)); for i=1:r D(i,i)=sqrt(B(i,n+1));