SVD分解及其应用

• SVD起源
  • 对角化概述
• SVD
• SVD应用
  • 图像压缩2
  • 数据去噪
  • LSA
  • 推荐系统
• 注意
• 参考资料

这里写图片描述

SVD可谓线性代数的登峰造极者。 其本质就是找到将任何一个矩阵对角化分解的两组标准正交的基底，同时对应的奇异值反映了对应基底变换的性质，为0表示对应的维度缺少信息，越大表明对应的维度容纳的信息方差越大。

SVD起源

认识一个问题总要追根究底，为什么要有SVD这个东西呢？ 要了解这一点，必须知道矩阵对角化是个奇妙的东西，以及并不是所有的矩阵都可以对角化。同时引出了今天的重点——如何让每个矩阵都可以对角化？

对角化概述

矩阵分析中，我们都想要好的矩阵，好的矩阵的一大特点就是可以对角化。

对角化的对象矩阵有两类：

1. 方矩阵的对角化
2. 长方形矩阵的对角化

对角化的方法也有两类：

1. 输入和输出空间的基完全一样，对应的特征值特征向量分解A=SΛS−1A = S\Lambda S^{-1}。但是这种情况下SS基中的向量不一定是正交的。当AA是对称矩阵的话SS中的基可以是标准正交的。同时基也不是一定存在的，只有有足够的特征向量，比如n×nn \times n的矩阵对角化的充要条件是有nn个不相关的特征向量。
2. 输入空间和输出空间的基不一样，这就对应了SVDSVD，也就是AV=USAV=US，V,UV,U分别是输入和输出空间的基，这种情况下的对角化总是存在的。

所以，综合对角化对象矩阵的形状以及对角化的方法，有以下结论：

1. 如果矩阵是nn阶方阵，可以尝试同一组基下面的对角化，也就是特征值特征向量分解。这种情况下对角化存在当且仅当存在nn个线性无关的特征向量。如果不存在的话说明不能找到同一组基使矩阵对角化。同时如果矩阵是对称矩阵的话，那么特征值肯定存在，肯定存在标准正交的特征向量。
2. 其他矩阵可以采用不同的两组基底实现SVDSVD对角化。
3. 对于nn阶方阵来说，eigeig分解与svdsvd分解相等当且仅当矩阵是半正定矩阵，也就是方阵在同一组标准正交基上对角化并且特定维度上的方向不变。
4. 特征值和奇异值分别表示对角化解耦后对应的基底的长度，从线性变换的角度上是对不同的基的延伸程度，从方差的角度上来说是方差的大小信息的多少。
5. 特征值或奇异值如果等于0，说明矩阵存在某一个维度上的信息缺失。因此可以得到如果矩阵AmnA_{mn}的秩为rr，那么它肯定有rr个不等于0的特征值和奇异值。

对角化的优点是（以特征值分解举例）：

1. 可以进行对角化分解，A=SΛS−1A= S \Lambda S^{-1}
2. 矩阵的kk次方Ak=SΛkS−1A^k =S \Lambda^k S^{-1}
3. 从对角化的矩阵中可以知道矩阵是不是缺失了某些维度的信息（特征值或者奇异值等于0），如果存在0的话那么矩阵不可逆（因为损失信息回不到以前了）。
4. 如果基底是标准正交基，那么从特征值或者奇异值的绝对值上可以找到哪个维度上的方差最大，利用这个思路可以实现数据压缩。

那么，具体如何将一个矩阵分解成对角矩阵和标准正交矩阵的乘积？

SVD

## 公式描述

AmnAmnAmn=UmrSrrVnrT=UmmSmnVnnT=u1σ1v1T+u1σ2v2T+⋯+urσ1vrT

\begin{split} A_{mn} &= U_{mr}S_{rr}{V_{nr}}^T \\ A_{mn} &= U_{mm}S_{mn}{V_{nn}}^T \\ A_{mn} &= \mathrm{u}_1 \sigma_1 {\mathrm{v}_1}^T+\mathrm{u}_1 \sigma_2 {\mathrm{v}_2}^T+\cdots +\mathrm{u}_r \sigma_1 {\mathrm{v}_r}^T \end{split} 上面的分解计算了奇异值σ\sigma不为0的情况，中间的分解考虑了奇异值σ\sigma为0的情况，最后的分解拆成了rr个列向量与行向量的乘积。 其中，

UVeigenValue(U)i=eigenVector(AAT)=eigenVector(ATA)=eigenValue(V)i=(diag(Sr)i)2i≤r

\begin{split} U &= eigenVector(AA^T)\\ V &= eigenVector(A^TA) \\ eigenValue(U)_i &= eigenValue(V)_i = ( diag(S_r)_i )^2 \quad i \le r \end{split} ## 几何描述

图片名称

SVD的几何意义是对于特定的矩阵AmnA_{mn}，寻找行空间RnR^n中的一组标准正交基，通过线性变换AA得到列空间RmR^m中的一组标准正交基。σi\sigma_i也可以理解成变换到AVAV空间后的模Avi=σiuiAv_i = \sigma_i u_i。

同时，如果矩阵的秩为rr，那么行空间、零空间、列空间、左零空间的示意与转换如上图所示。

采用matlab中的eigshow函数，可以得到类似的解释：

通过选择一组合适的正交基VV，使得AV=USAV=US也是正交的。然后分别以这两组正交且单位化的V,UV,U为基，SS中包含了他们的比例系数，构建了对角化的矩阵SS，实现了对角化解耦的线性变换。

SVD应用

##图像压缩 我们的目标就是这位美女——蒋勤勤，将对这幅图片实现压缩。

图片名称

首先，先来看看原来的RGB图片以及RGB分量的灰度图片：

图片名称

SVD之后，先来看看SVD的奇异值大小的分布情况和累计分布比率：

接着，看一看选择不同数量的奇异值的结果。

这里写图片描述

最后，看一看图像的压缩比，原图像是686×482686 \times 482，选择最大的100个奇异值已经能够得到相当好的结果了。这时候奇异值的累积比例为89.6%89.6\%，压缩比是

686×482(686+482)×100=2.83

\frac{686 \times 482}{ (686+482) \times 100} = 2.83

具体代码如下：

% 读入图像RGB数据，并从Uint8转换成double类型方便之后的处理。
p = imread('/Users/yangguangyao/Desktop/test/p.jpg');
pr = p(:,:,1);
pg = p(:,:,2);
pb = p(:,:,3);
pr = double(pr);
pg = double(pg);
pb = double(pb);

% 可视化图像
figure()
subplot(2,2,1);
imshow(p)
title('原来的RGB图像')
subplot(2,2,2);
imshow(mat2gray(pr))
title('R分量的灰度图像')
subplot(2,2,3);
imshow(mat2gray(pg))
title('G分量的灰度图像')
subplot(2,2,4);
imshow(mat2gray(pb))
title('B分量的灰度图像')

% SVD分解
[Ur,Sr,Vr] = svd(pr);
[Ug,Sg,Vg] = svd(pg);
[Ub,Sb,Vb] = svd(pb);

% 分析SVD，计划选取1 3 5 10 30 50 100 150
svdD = diag(Sr);
cumsumD = cumsum(svdD);
plot(svdD,'LineWidth',2)
plot(cumsumD,'LineWidth',2)

% 分解后按照singular value从大到小选择
fr = @(n) Ur(:,1:n)*Sr(1:n,1:n)*Vr(:,1:n)';
fg = @(n) Ug(:,1:n)*Sg(1:n,1:n)*Vg(:,1:n)';
fb = @(n) Ub(:,1:n)*Sb(1:n,1:n)*Vb(:,1:n)';

param = [1,3,5,10,30,50,100,150];
figure()
for i = 1:8
    subplot(2,4,i);
    n = param(i);
    pnew(:,:,1) = fr(n);
    pnew(:,:,2) = fg(n);
    pnew(:,:,3) = fb(n);
    pnew = uint8(pnew);
    imshow(pnew)
    title(strcat(num2str(param(i)),'个奇异值'))
end

图像压缩2

这里写图片描述

上图是一个15×2515 \times 25的图像，其本质上是15×2515 \times 25的矩阵，白色的元素代表相应位置上是1，黑色代表相应位置上是0。其一共有三种模式：

这里写图片描述

这个矩阵的特征值有3，对应了三种的模式，选择最大的三个奇异值进行SVD后的结果是：

这里写图片描述

代码如下：

% 构建目标图像矩阵X
x1 = ones(25,1);
x2 = [ones(5,1);zeros(15,1);ones(5,1);];
x3 = [ones(5,1);zeros(3,1);ones(9,1);zeros(3,1);ones(5,1);];
X = [repmat(x1,1,2),repmat(x2,1,3),repmat(x3,1,5),repmat(x2,1,3),repmat(x1,1,2)];

% 查看矩阵的秩
r = rank(X);

% 进行SVD
[U,S,V] = svd(X);
Xnew = U(:,1:r)*S(1:r,1:r)*V(:,1:r)';

% 分析画图
figure()
subplot(1,2,1)
imshow(X)
title('原始的图像')
subplot(1,2,2)
imshow(Xnew)
title('选择秩数量3的奇异值分解')

数据去噪

SVD总能找到标准化正交基后方差最大的维度，因此可以用它进行降维去噪等等。

下面分别用SVD和linear regression去拟合直线，结果如下，看来效果还不错哦。

这里写图片描述

代码如下：

% 模拟线性数据，加上一定的高斯噪声
X = 1:10;
Xb = ones(1,10);
Y = 2*X + random('Normal',0,1,1,10);

% 进行SVD分解并选择原输入空间的一个奇异值比较大的基，实现了数据降维
M = [X;Y];
[U,S,V] = svd(M);

U1 = U(:,1);
u1 = U1(1);
u2 = U1(2);
k = (u2/u1);

% 进行线性回归
w = pinv([X;Xb])'*Y';

% 分别是SVD和线性回归拟合的数据
Y1 = k*X;
Y2 = w(1)*X+w(2);

% 画图并比较
figure()
% 注释蛮方便的函数ezplot('y-2*x-1')
% refline(u2/u1,0) 
hold on
plot(X,Y,'ko','LineWidth',2)
plot(X,Y1,'r-','LineWidth',1)
plot(X,Y2,'b-','LineWidth',1)
legend('数据点','SVD拟合','线性回归拟合')

LSA

LSA（Latent Semantic Analysis也叫作Latent Semantic Indexing）分析文档发现潜在的概念意义。

推荐系统

注意

1. USVTUSV^T分解和SΛS−1S\Lambda S^{-1}分解等价当且仅当矩阵是半正定（对称矩阵，有大于等于0的特征值）。
2. 介绍下几种矩阵分解
   • A=USVTA=USV^T，U,VU,V都是标准正交的，S是对角化的上面按照由大到小的顺序存放着奇异值。
   • A=QHRA=QHR，QQ是标准正交的，RR对角线上是1，HH是对角化的上面存放着高度hih_i。
   • A=LDUA=LDU，L,UL,U的对角线上是1，主元存放在DD中。

A的行列式的绝对值等于奇异值的乘积，等于主元乘积的绝对值，等于高度hih_i乘积的绝对值，等于特征值的乘积。

参考资料

1. 关于SVD很棒的包含几何解释的资料
2. 关于LSA的很棒的文章