奇异值（Singular value decomposition SVD)分解

本文摘自两篇博客，感谢博主分享

一、原文地址：http://blog.csdn.net/wangzhiqing3/article/details/7446444

SVD分解

SVD分解是LSA的数学基础，本文是我的LSA学习笔记的一部分，之所以单独拿出来，是因为SVD可以说是LSA的基础，要理解LSA必须了解SVD，因此将LSA笔记的SVD一节单独作为一篇文章。本节讨论SVD分解相关数学问题，一个分为3个部分，第一部分讨论线性代数中的一些基础知识，第二部分讨论SVD矩阵分解，第三部分讨论低阶近似。本节讨论的矩阵都是实数矩阵。

基础知识

数学知识不清楚可以参考：http://student.zjzk.cn/course_ware/web-gcsx/main.htm

1. 矩阵的秩：矩阵的秩是矩阵中线性无关的行或列的个数 参考：http://student.zjzk.cn/course_ware/web-gcsx/gcsx/chapter3/chapter3.1.htm

2. 对角矩阵：对角矩阵是除对角线外所有元素都为零的方阵

3. 单位矩阵：如果对角矩阵中所有对角线上的元素都为1，该矩阵称为单位矩阵

4. 特征值：对一个M x M矩阵C和向量X，如果存在λ使得下式成立



 

则称λ为矩阵C的特征值，X称为矩阵的特征向量。非零特征值的个数小于等于矩阵的秩。

5. 特征值和矩阵的关系：考虑以下矩阵





该矩阵特征值λ1 = 30,λ2 = 20,λ3 = 1。对应的特征向量





假设VT=(2,4,6) 计算S x VT









有上面计算结果可以看出，矩阵与向量相乘的结果与特征值，特征向量有关。观察三个特征值λ1 = 30,λ2 = 20,λ3 = 1，λ3值最小，对计算结果的影响也最小，如果忽略λ3，那么运算结果就相当于从(60,80,6)转变为(60,80,0)，这两个向量十分相近。这也表示了数值小的特征值对矩阵-向量相乘的结果贡献小，影响小。这也是后面谈到的低阶近似的数学基础。

矩阵分解

1. 方阵的分解

1) 设S是M x M方阵，则存在以下矩阵分解





其中U 的列为S的特征向量，

为对角矩阵，其中对角线上的值为S的特征值，按从大到小排列：





2) 设S是M x M 方阵，并且是对称矩阵，有M个特征向量。则存在以下分解





其中Q的列为矩阵S的单位正交特征向量，

仍表示对角矩阵，其中对角线上的值为S的特征值，按从大到小排列。最后，QT=Q-1，因为正交矩阵的逆等于其转置。

2. 奇异值分解

上面讨论了方阵的分解，但是在LSA中，我们是要对Term-Document矩阵进行分解，很显然这个矩阵不是方阵。这时需要奇异值分解对Term-Document进行分解。奇异值分解的推理使用到了上面所讲的方阵的分解。

假设C是M x N矩阵，U是M x M矩阵，其中U的列为CCT的正交特征向量，V为N x N矩阵，其中V的列为CTC的正交特征向量，再假设r为C矩阵的秩，则存在奇异值分解：





其中CCT和CTC的特征值相同，为



Σ为M X N，其中



，其余位置数值为0，

的值按大小降序排列。以下是Σ的完整数学定义：





σi称为矩阵C的奇异值。

用C乘以其转置矩阵CT得：





上式正是在上节中讨论过的对称矩阵的分解。

奇异值分解的图形表示：





从图中可以看到Σ虽然为M x N矩阵，但从第N+1行到M行全为零，因此可以表示成N x N矩阵，又由于右式为矩阵相乘，因此U可以表示为M x N矩阵，VT可以表示为N x N矩阵

3. 低阶近似

LSA潜在语义分析中，低阶近似是为了使用低维的矩阵来表示一个高维的矩阵，并使两者之差尽可能的小。本节主要讨论低阶近似和F-范数。

给定一个M x N矩阵C(其秩为r)和正整数k，我们希望找到一个M x N矩阵Ck，其秩不大于K。设X为C与Ck之间的差，X=C – Ck，X的F-范数为





当k远小于r时，称Ck为C的低阶近似，其中X也就是两矩阵之差的F范数要尽可能的小。

SVD可以被用与求低阶近似问题，步骤如下：

1. 给定一个矩阵C，对其奇异值分解：



2. 构造

，它是将

的第k+1行至M行设为零，也就是把

的最小的r-k个(the
r-k smallest)奇异值设为零。

3. 计算Ck：



回忆在基础知识一节里曾经讲过，特征值数值的大小对矩阵-向量相乘影响的大小成正比，而奇异值和特征值也是正比关系，因此这里选取数值最小的r-k个特征值设为零合乎情理，即我们所希望的C-Ck尽可能的小。完整的证明可以在Introduction to Information Retrieval[2]中找到。

我们现在也清楚了LSA的基本思路：LSA希望通过降低传统向量空间的维度来去除空间中的“噪音”，而降维可以通过SVD实现，因此首先对Term-Document矩阵进行SVD分解，然后降维并构造语义空间。

二、下面这篇文字深入浅出解释了奇异值分解的意义：

原文地址：http://blog.csdn.net/xiahouzuoxin/article/details/41118351

特征值与特征向量的几何意义

矩阵的乘法是什么，别只告诉我只是“前一个矩阵的行乘以后一个矩阵的列”，还会一点的可能还会说“前一个矩阵的列数等于后一个矩阵的行数才能相乘”，然而，这里却会和你说——那都是表象。

矩阵乘法真正的含义是变换，我们学《线性代数》一开始就学行变换列变换，那才是线代的核心——别会了点猫腻就忘了本——对，矩阵乘法 

就是线性变换，若以其中一个向量A为中心，则B的作用主要是使A发生如下变化：

伸缩

clf;
A = [0, 1, 1, 0, 0;...
1, 1, 0, 0, 1];  % 原空间
B = [3 0; 0 2];      % 线性变换矩阵

plot(A(1,:),A(2,:), '-*');hold on
grid on;axis([0 3 0 3]); gtext('变换前');

Y = B * A;

plot(Y(1,:),Y(2,:), '-r*');
grid on;axis([0 3 0 3]); gtext('变换后');

1

从上图可知，y方向进行了2倍的拉伸，x方向进行了3倍的拉伸，这就是B=[3 0; 0 2]的功劳,3和2就是伸缩比例。请注意，这时B除了对角线元素为各个维度的倍数外，非正对角线元素都为0，因为下面将要看到，对角线元素非0则将会发生切变及旋转的效果。

切变

clf;
A = [0, 1, 1, 0, 0;...
1, 1, 0, 0, 1];  % 原空间
B1 = [1 0; 1 1];       % 线性变换矩阵
B2 = [1 0; -1 1];       % 线性变换矩阵
B3 = [1 1; 0 1];       % 线性变换矩阵
B4 = [1 -1; 0 1];       % 线性变换矩阵

Y1 = B1 * A;
Y2 = B2 * A;
Y3 = B3 * A;
Y4 = B4 * A;

subplot(2,2,1);
plot(A(1,:),A(2,:), '-*'); hold on;plot(Y1(1,:),Y1(2,:), '-r*');
grid on;axis([-1 3 -1 3]);
subplot(2,2,2);
plot(A(1,:),A(2,:), '-*'); hold on;plot(Y2(1,:),Y2(2,:), '-r*');
grid on;axis([-1 3 -1 3]);
subplot(2,2,3);
plot(A(1,:),A(2,:), '-*'); hold on;plot(Y3(1,:),Y3(2,:), '-r*');
grid on;axis([-1 3 -1 3]);
subplot(2,2,4);
plot(A(1,:),A(2,:), '-*'); hold on;plot(Y4(1,:),Y4(2,:), '-r*');
grid on;axis([-1 3 -1 3]);

2

旋转

所有的变换其实都可以通过上面的伸缩和切变变换的到，如果合理地对变换矩阵B取值，能得到图形旋转的效果，如下，

clf;
A = [0, 1, 1, 0, 0;...
1, 1, 0, 0, 1];  % 原空间
theta = pi/6;
B = [cos(theta) sin(theta); -sin(theta) cos(theta)];
Y = B * A;
figure;
plot(A(1,:),A(2,:), '-*'); hold on;plot(Y(1,:),Y(2,:), '-r*');
grid on;axis([-1 3 -1 3]);

3

好，关于矩阵乘就这些了。那么，我们接着就进入主题了，对特定的向量，经过一种方阵变换，经过该变换后，向量的方向不变（或只是反向），而只是进行伸缩变化（伸缩值可以是负值，相当于向量的方向反向）？这个时候我们不妨将书上对特征向量的定义对照一遍：

数学教材定义： 设A是n阶方阵，如果存在 λ 和n维非零向量X，使 

 ，则 λ 称为方阵A的一个特征值，X为方阵A对应于或属于特征值 λ 的一个特征向量。

上面特定的向量不就是特征向量吗？ λ 不就是那个伸缩的倍数吗？因此，特征向量的代数上含义是：将矩阵乘法转换为数乘操作；特征向量的几何含义是：特征向量通过方阵A变换只进行伸缩，而保持特征向量的方向不变。特征值表示的是这个特征到底有多重要，类似于权重，而特征向量在几何上就是一个点，从原点到该点的方向表示向量的方向。

特征向量有一个重要的性质：同一特征值的任意多个特征向量的线性组合仍然是A属于同一特征值的特征向量。关于特征值，网上有一段关于“特征值是震动的谱”的解释：

戏说在朝代宋的时候，我国就与发现矩阵特征值理论的机会擦肩而过。话说没有出息的秦少游在往池塘里扔了一颗小石头后，刚得到一句“投石冲开水底天”的泡妞诗对之后，就猴急猴急地去洞房了，全然没有想到水波中隐含着矩阵的特征值及特征向量的科学大道理。大概地说，水面附近的任一点水珠在原处上下振动（实际上在做近似圆周运动），并没有随着波浪向外圈移动，同时这些上下振动的水珠的幅度在渐渐变小，直至趋于平静。在由某块有着特定质量和形状的石头被以某种角度和速度投入某个面积和深度特定的水池中所决定的某个矩阵中，纹波荡漾中水珠的渐变过程中其特征值起着决定性的作用，它决定着水珠振动的频率和幅度减弱的衰退率。

在理解关于振动的特征值和特征向量的过程中，需要加入复向量和复矩阵的概念，因为在实际应用中，实向量和实矩阵是干不了多少事的。机械振动和电振动有频谱，振动的某个频率具有某个幅度；那么矩阵也有矩阵的谱，矩阵的谱就是矩阵特征值的概念，是矩阵所固有的特性，所有的特征值形成了矩阵的一个频谱，每个特征值是矩阵的一个“谐振频点”。

美国数学家斯特让（G..Strang）在其经典教材《线性代数及其应用》中这样介绍了特征值作为频率的物理意义，他说：

大概最简单的例子（我从不相信其真实性，虽然据说1831年有一桥梁毁于此因）是一对士兵通过桥梁的例子。传统上，他们要停止齐步前进而要散步通过。这个理由是因为他们可能以等于桥的特征值之一的频率齐步行进，从而将发生共振。就像孩子的秋千那样，你一旦注意到一个秋千的频率，和此频率相配，你就使频率荡得更高。一个工程师总是试图使他的桥梁或他的火箭的自然频率远离风的频率或液体燃料的频率；而在另一种极端情况，一个证券经纪人则尽毕生精力于努力到达市场的自然频率线。特征值是几乎任何一个动力系统的最重要的特征。

其实，这个矩阵之所以能形成“频率的谱”，就是因为矩阵在特征向量所指的方向上具有对向量产生恒定的变换作用：增强（或减弱）特征向量的作用。进一步的，如果矩阵持续地叠代作用于向量，那么特征向量的就会凸现出来。

更多关于特征向量及特征值的实际例子参见Wikipedia： http://zh.wikipedia.org/wiki/特征向量 。

特征值分解

设A有n个特征值及特征向量，则：



将上面的写到一起成矩阵形式：



若(x1,x2,...,xn)可逆，则左右两边都求逆，则方阵A可直接通过特征值和特征向量进行唯一的表示，令

Q=(x1,x2,...,xn)

Σ = diag(λ1, λ2, ..., λn)

则 

 ，该表达式称为方阵的特征值分解，这样方阵A就被特征值和特征向量唯一表示。

一个变换方阵的所有特征向量组成了这个变换矩阵的一组基。所谓基，可以理解为坐标系的轴。我们平常用到的大多是直角坐标系，在线性代数中可以把这个坐标系扭曲、拉伸、旋转，称为基变换。我们可以按需求去设定基，但是基的轴之间必须是线性无关的，也就是保证坐标系的不同轴不要指向同一个方向或可以被别的轴组合而成，否则的话原来的空间就“撑”不起来了。从线性空间的角度看，在一个定义了内积的线性空间里，对一个N阶对称方阵进行特征分解，就是产生了该空间的N个标准正交基，然后把矩阵投影到这N个基上。N个特征向量就是N个标准正交基，而特征值的模则代表矩阵在每个基上的投影长度。特征值越大，说明矩阵在对应的特征向量上的方差越大，功率越大，信息量越多。不过，特征值分解也有很多的局限，比如说变换的矩阵必须是方阵。

在机器学习特征提取中，意思就是最大特征值对应的特征向量方向上包含最多的信息量，如果某几个特征值很小，说明这几个方向信息量很小，可以用来降维，也就是删除小特征值对应方向的数据，只保留大特征值方向对应的数据，这样做以后数据量减小，但有用信息量变化不大，PCA降维就是基于这种思路。

Matlab中通过eig函数就可求得特征值和特征向量矩阵。

>> B = [ 3     -2      -.9    2*eps
-2      4       1    -eps
-eps/4  eps/2  -1     0
-.5    -.5      .1    1   ]
B =
3.0000   -2.0000   -0.9000    0.0000
-2.0000    4.0000    1.0000   -0.0000
-0.0000    0.0000   -1.0000         0
-0.5000   -0.5000    0.1000    1.0000

>> [V D] = eig(B)
V =
0.6153   -0.4176   -0.0000   -0.1437
-0.7881   -0.3261   -0.0000    0.1264
-0.0000   -0.0000   -0.0000   -0.9196
0.0189    0.8481    1.0000    0.3432
D =
5.5616         0         0         0
0    1.4384         0         0
0         0    1.0000         0
0         0         0   -1.0000

D对角线的元素即为特征值（表示了伸缩的比例），D就是特征值分解公式中的Q，V的每一列与D没列对应，表示对应的特征向量，即特征值分解中的Σ。

奇异值分解

特征值分解是一个提取矩阵特征很不错的方法，但是它只适用于方阵。而在现实的世界中，我们看到的大部分矩阵都不是方阵，比如说有M个学生，每个学生有N科成绩，这样形成的一个M * N的矩阵就可能不是方阵，我们怎样才能像描述特征值一样描述这样一般矩阵呢的重要特征呢？奇异值分解就是用来干这个事的，奇异值分解是一个能适用于任意的矩阵的一种分解的方法。我们有必要先说说特征值和奇异值之间的关系。

对于特征值分解公式， ATA 是方阵，我们求 ATA 的特征值，即 

 ，此时求得的特征值就对应奇异值的平方，求得的特征向量v称为右奇异向量，另外还可以得到：





所求的ui就是左奇异向量， σi 就是奇异值。已有人对SVD的几何机理做了清晰的分析，非常受用，就不重复造轮子，下文为转载自http://blog.sciencenet.cn/blog-696950-699432.html 。



SVD分解

SVD之所以很有效，是因为：在很多情况下，前10%甚至1%的奇异值的和就占了全部的奇异值之和的99%以上了。在这里，我们用图像简单的实践一下SVD的大妙处，下面是matlab对图像进行SVD分解的例子，

<code class="sourceCode matlab" style="margin: 0px; padding: 0.5em; line-height: normal; font-family: Monaco, 'Courier New', 'DejaVu Sans Mono', 'Bitstream Vera Sans Mono', monospace; color: rgb(51, 51, 51); border: none; display: block; background-color: rgb(248, 248, 255); background-position: initial initial; background-repeat: initial initial;">I = imread(<span class="st" style="margin: 0px; padding: 0px; color: rgb(64, 112, 160);">'lena_gray.bmp'</span>);                  <span class="co" style="margin: 0px; padding: 0px; color: rgb(96, 160, 176); font-style: italic;">% 512x512的Lena图像</span>
im = double(I);
[s,v,d]=svd(im);                              <span class="co" style="margin: 0px; padding: 0px; color: rgb(96, 160, 176); font-style: italic;">% svd分解，svd分解后特征值v对角线按从大到小排列，因此可以选择特征值大的进行恢复</span>
recv1=s(:,<span class="fl" style="margin: 0px; padding: 0px; color: rgb(64, 160, 112);">1</span>:<span class="fl" style="margin: 0px; padding: 0px; color: rgb(64, 160, 112);">20</span>)*v(<span class="fl" style="margin: 0px; padding: 0px; color: rgb(64, 160, 112);">1</span>:<span class="fl" style="margin: 0px; padding: 0px; color: rgb(64, 160, 112);">20</span>,<span class="fl" style="margin: 0px; padding: 0px; color: rgb(64, 160, 112);">1</span>:<span class="fl" style="margin: 0px; padding: 0px; color: rgb(64, 160, 112);">50</span>)*d(:,<span class="fl" style="margin: 0px; padding: 0px; color: rgb(64, 160, 112);">1</span>:<span class="fl" style="margin: 0px; padding: 0px; color: rgb(64, 160, 112);">50</span>)';      <span class="co" style="margin: 0px; padding: 0px; color: rgb(96, 160, 176); font-style: italic;">% svd取最高的100个特征值进行恢复</span>
recv2=s(:,<span class="fl" style="margin: 0px; padding: 0px; color: rgb(64, 160, 112);">1</span>:<span class="fl" style="margin: 0px; padding: 0px; color: rgb(64, 160, 112);">50</span>)*v(<span class="fl" style="margin: 0px; padding: 0px; color: rgb(64, 160, 112);">1</span>:<span class="fl" style="margin: 0px; padding: 0px; color: rgb(64, 160, 112);">50</span>,<span class="fl" style="margin: 0px; padding: 0px; color: rgb(64, 160, 112);">1</span>:<span class="fl" style="margin: 0px; padding: 0px; color: rgb(64, 160, 112);">100</span>)*d(:,<span class="fl" style="margin: 0px; padding: 0px; color: rgb(64, 160, 112);">1</span>:<span class="fl" style="margin: 0px; padding: 0px; color: rgb(64, 160, 112);">100</span>)';    <span class="co" style="margin: 0px; padding: 0px; color: rgb(96, 160, 176); font-style: italic;">% svd取最高的100个特征值进行恢复</span>
recv3=s(:,<span class="fl" style="margin: 0px; padding: 0px; color: rgb(64, 160, 112);">1</span>:<span class="fl" style="margin: 0px; padding: 0px; color: rgb(64, 160, 112);">200</span>)*v(<span class="fl" style="margin: 0px; padding: 0px; color: rgb(64, 160, 112);">1</span>:<span class="fl" style="margin: 0px; padding: 0px; color: rgb(64, 160, 112);">200</span>,<span class="fl" style="margin: 0px; padding: 0px; color: rgb(64, 160, 112);">1</span>:<span class="fl" style="margin: 0px; padding: 0px; color: rgb(64, 160, 112);">200</span>)*d(:,<span class="fl" style="margin: 0px; padding: 0px; color: rgb(64, 160, 112);">1</span>:<span class="fl" style="margin: 0px; padding: 0px; color: rgb(64, 160, 112);">200</span>)';  <span class="co" style="margin: 0px; padding: 0px; color: rgb(96, 160, 176); font-style: italic;">% svd取最高的100个特征值进行恢复</span>

subplot(<span class="fl" style="margin: 0px; padding: 0px; color: rgb(64, 160, 112);">221</span>);imshow(I);
title(<span class="st" style="margin: 0px; padding: 0px; color: rgb(64, 112, 160);">'原图'</span>);
subplot(<span class="fl" style="margin: 0px; padding: 0px; color: rgb(64, 160, 112);">222</span>);imshow(mat2gray(recv1));
title(<span class="st" style="margin: 0px; padding: 0px; color: rgb(64, 112, 160);">'恢复：左奇异20、右奇异50'</span>);
subplot(<span class="fl" style="margin: 0px; padding: 0px; color: rgb(64, 160, 112);">223</span>);imshow(mat2gray(recv2));
title(<span class="st" style="margin: 0px; padding: 0px; color: rgb(64, 112, 160);">'恢复：左奇异50、右奇异100'</span>);
subplot(<span class="fl" style="margin: 0px; padding: 0px; color: rgb(64, 160, 112);">224</span>);imshow(mat2gray(recv3));
title(<span class="st" style="margin: 0px; padding: 0px; color: rgb(64, 112, 160);">'恢复：左奇异200、右奇异200'</span>);</code>

图注：SVD二维图像压缩恢复

如果按左下角的方式压缩原图，则存储量变为：50x50+100x100+50=12500，而存储原图像的存储量为512x512=262144，则压缩比为262144/12500=20.97，这里没有考虑存储数据类型的差异。

SVD分解相对于特征值分解的优势就是：
分解的矩阵可以是任意矩阵
在恢复信号的时候左右奇异值可以选择不同的维度

另外值得注意的一点：不论是奇异值分解还是特征值分解，分解出来的特征向量都是正交的。

奇异值分解与PCA

关于奇异值与PCA的关系， http://www.cnblogs.com/LeftNotEasy/archive/2011/01/19/svd-and-applications.html 给了很好的解释，也直接整理过来，感谢原作者：



图注：SVD与PCA

PCA就是一种用于对数据进行降维的方法（降维肯定会丢失数据，只不过能在减少大量存储量的同时损失尽可能少），参见之前matlab对图像进行SVD分解的例子，更容易理解：实现了SVD就实现了PCA，PCA仅是SVD的包装。

PCA的应用很广，主要用在机器学习中对特征进行降维，还能用于去噪，下面两图是PCA降维和PCA去噪的例子（图片来自邹博PPT：北京9月秋季班·机器学习初步）



图注：PCA降维

降维说白了就是将信息通过投影到更低得多维度，这样必然会带来信息的丢失，但就如上图，这种信息的丢失却有时对分类没有影响，反而能降低识别算法的维度，提高速度，缓解所谓的维度灾难。



PCA去噪的前提是噪声的特征值会比信号的特征值小，即信噪比高的情况，否则PCA去噪会产生逆效果——把信号去掉了而噪声没去掉。

SVD其它

SVD还有其它很多方面的应用，通过查找资料，这里先做个简单的罗列，有机会再一个个研究：
求伪逆。我们知道，矩阵求逆要求矩阵必须是方阵，SVD可以用来求任意矩阵的逆运算，求任意矩阵的逆矩阵称为求伪逆
最小二乘法求解。凭着对《矩阵论》的零星的记忆，SVD算法就是因为能求伪逆所以用来求解最小二乘法。
基于SVD的文本分类。首先接触是从吴军老师的《数学之美》一书上看到的，大致是：通过TF/IDF(term frequency/inverse document frequency)构建“一百万篇文章和五十万词的关联性”的矩阵 A1000000x500000 ，然后对A矩阵使用SVD分解之后，存储量减小，而且左奇异值矩阵和右奇异值矩阵以及奇异值矩阵的物理含义将非常明晰，http://blog.csdn.net/wangzhiqing3/article/details/7431276 有简单介绍，更多参见吴军老师的《数学之美》

另外，开源视觉库OpenCV中也提供SVD分解的算法。

参考

强大的矩阵奇异值分解(SVD)及其应用
http://blog.sciencenet.cn/blog-696950-699432.html ，关于SVD部分的讲解即出自此处。
http://www.cise.ufl.edu/~srajaman/Rajamanickam_S.pdf
July培训班PPT：“北京9月秋季班：贝叶斯、PCA、SVD”.
戴华编著《矩阵论》 教材，科学出版社. 有时候明白一点之后再回头参考数学教材要比网上大把的抓好得多，豁然开朗。