矩阵的奇异值分解基本概念

1. 背景知识

矩阵的奇异值分解，可以简单的理解为特征值分解问题的推广，对于一个m×n（m≥n）的矩阵A，可以具有如下形式的分解：

A=U×S×V'

其中U和V分别是m×m和n×n的正交阵，S则是一个m×n的对角阵（当然在n+1~m行，全是零）。

S就是所谓的奇异值（Singular Value）矩阵，U和V分别被称为左/右奇异向量（Singular Vector）。

对奇异值S的定义是：一个m×n矩阵A的奇异值就是n×n对称阵（A'×A）的非零特征值的正平方根。

通常，为了便于组织和描述，会将S的奇异值按从大到小的顺序排列，例如MATLAB中使用函数svd进行奇异值分解：

>> A=[1 0 1;0 1 0;0 1 1;0 1 0;1 1 0];
>> [U,S,V]=svd(A)

U =

   -0.3651    0.8165    0.0000    0.1184   -0.4313
   -0.3651   -0.4082   -0.0000   -0.5635   -0.6185
   -0.5477   -0.0000    0.7071   -0.1184    0.4313
   -0.3651   -0.4082   -0.0000    0.8003   -0.2441
   -0.5477    0.0000   -0.7071   -0.1184    0.4313


S =

    2.2361         0         0
         0    1.4142         0
         0         0    1.0000
         0         0         0
         0         0         0


V =

   -0.4082    0.5774   -0.7071
   -0.8165   -0.5774   -0.0000
   -0.4082    0.5774    0.7071

关于左/右奇异向量的确定

从定义可以知道，对奇异值向量可以很容易确定，但是对于左右奇异向量，在定义中似乎并没有明确的方法，不过还是可以通过固定的过程进行计算来确定。

首先，做一点点简短的矩阵计算：

A'×A=V×D×V'

左边是一个对称阵，一定可以分解为对角阵D和一个正交阵V，在这种情况下，如果特征值按降序排列，则D是唯一的，而且，此时D的对角线上的元素值就是奇异值的平方，将其开方，就可以构成奇异值矩阵。

但是，由于A'×A可能存在的多重特征值，所以特征向量的选取有可能并不是唯一的，但可以肯定，V一定是存在的。

确定左奇异矩阵，可以使用如下方式：

通过计算A×A'所对应的一组标准正交特征向量，将其构造成m×m的正交矩阵。就可以实现A=U×S×V'。

计算A×A'的特征值，其实和A'×A的特征值一样，只是一个的有(m-n)个特征值为0，这样的话，也会导致U的选择也不唯一，但同样的，U肯定存在。

将公式总结一下，当A=U×S×V'时：

A×A'=(U×S×V)×(U×S×V)'=U×S×V×V'×S'×U'=U×D×U'

A'×A=(U×S×V)'×(U×S×V)=V'×S'×U'×U×S×V=V×D×V'

注意两个公式中，上面的的D是m×m的，下面的D是n×n的。

奇异值分解的验算

可以通过MATLAB做个实验验算一下。

% usv.m
% A=U×S×V' method test
% you can use matlab function svd 
A=[1 0 1;0 1 0;0 1 1;0 1 0;1 1 0];
B=A'*A;
[V,D1]=eig(B);
[n1,~]=size(B);
D1=sqrt(D1);
%sorted D1 and corresponding V
for p=1:n1-1
	for q=p+1:n1
		if D1(p,p)<D1(q,q)
			t=D1(p,p);
			D1(p,p)=D1(q,q);
			D1(q,q)=t;
			V(:,[p,q])=V(:,[q,p]);
		end
	end
end
C=A*A';
[U,D2]=eig(C);
[n2,~]=size(C);
D2=sqrt(D2);
%sorted D2 and corresponding U
for p=1:n2-1
	for q=p+1:n2
		if D2(p,p)<D2(q,q)
			t=D2(p,p);
			D2(p,p)=D2(q,q);
			D2(q,q)=t;
			U(:,[p,q])=U(:,[q,p]);
		end
	end
end
S=D2(:,1:n1);
% Ax should be equal to A (of course within tolerance)
Ax=U*S*V'; 
diff=Ax-A;
disp(norm(diff,inf));

运行结果为：

>> usv
   2.1094e-15

也可以通过MATLAB官方函数svd校验：

>> [U0,S0,V0]=svd(A)

U0 =

   -0.3651    0.8165    0.0000    0.1184   -0.4313
   -0.3651   -0.4082   -0.0000   -0.5635   -0.6185
   -0.5477   -0.0000    0.7071   -0.1184    0.4313
   -0.3651   -0.4082   -0.0000    0.8003   -0.2441
   -0.5477    0.0000   -0.7071   -0.1184    0.4313


S0 =

    2.2361         0         0
         0    1.4142         0
         0         0    1.0000
         0         0         0
         0         0         0


V0 =

   -0.4082    0.5774   -0.7071
   -0.8165   -0.5774   -0.0000
   -0.4082    0.5774    0.7071

>> norm(U0*S*V0'-A,inf)

ans =

   4.4409e-16

可以看到，官方算法具有更小的误差。同样可以对比一下左/右奇异矩阵：

>> V

V =

    0.4082    0.5774    0.7071
    0.8165   -0.5774    0.0000
    0.4082    0.5774   -0.7071

>> V0

V0 =

   -0.4082    0.5774   -0.7071
   -0.8165   -0.5774   -0.0000
   -0.4082    0.5774    0.7071

>> U

U =

    0.3651    0.8165    0.0000   -0.3955    0.2087
    0.3651   -0.4082   -0.0000   -0.0656    0.8341
    0.5477    0.0000   -0.7071    0.3955   -0.2087
    0.3651   -0.4082    0.0000   -0.7255   -0.4167
    0.5477   -0.0000    0.7071    0.3955   -0.2087

>> U0

U0 =

   -0.3651    0.8165    0.0000    0.1184   -0.4313
   -0.3651   -0.4082   -0.0000   -0.5635   -0.6185
   -0.5477   -0.0000    0.7071   -0.1184    0.4313
   -0.3651   -0.4082   -0.0000    0.8003   -0.2441
   -0.5477    0.0000   -0.7071   -0.1184    0.4313

可以看到，显然V有些列向量的符号不相同，此外，对于U的最后两列，则是完全不同，说明U和V的选取可能并不唯一。