机器学习中的矩阵方法04:SVD 分解

 

前面我们讲了 QR 分解有一些优良的特性,但是 QR 分解仅仅是对矩阵的行进行操作(左乘一个酉矩阵),可以得到列空间。这一小节的 SVD 分解则是将行与列同等看待,既左乘酉矩阵,又右乘酉矩阵,可以得出更有意思的信息。奇异值分解( SVD, Singular Value Decomposition ) 在计算矩阵的伪逆( pseudoinverse ),最小二乘法最优解,矩阵近似,确定矩阵的列向量空间,秩以及线性系统的解集空间都有应用。

 

1. SVD 的形式

对于一个任意的 m×n 的矩阵 A,SVD 将一个矩阵分解为三个特殊矩阵的乘积, 国外还做了一个视频——SVD 之歌

其中, U 和 是酉矩阵, 是对角线矩阵。注意酉矩阵只是坐标的转换,数据本身分布的形状并没有改变,而对角矩阵,则是对数据进行了拉伸或者压缩。由于 m >= n,又可以写成下面这种 thin SVD 形式:

2. SVD 的几何解释

考虑 A 是 2×2 的简单情况。我们知道,一个几何形状左乘一个矩阵 A 实际上就是将该形状进行旋转、对称、拉伸变换等错切变换, A 就是所谓的 shear 矩阵。比如可以将一个圆通过左乘 A 得到一个旋转后的椭圆。

举个例子,假设 A 对平面上的一个圆进行变换:

如果只看矩阵 A, 我们几乎无法直观看到圆是如何变换的。变换前的圆是:

C, M, Y, K 分别表示第一、二、三、四象限。左乘 A 矩阵后得到:

我们的问题是,旋转了多少度?伸缩的方向是多少?最大伸缩比例是多少?

利用这个在线工具对 A 进行 SVD 分解:

明显,我们看到 A 变换实际上是先对圆顺时针旋转 45°(可以看做是坐标轴逆时针宣战了 45°,主成分方向),再关于 x 轴对称(第一行乘以 -1), 即左乘 V^T:

然后在 x 方向拉伸 3 倍(左乘S):

最后再顺时针旋转 45°,再关于 x 轴对称(第一行乘以 -1, 两次对称操作抵消了), 即左乘 U:

上述绘图过程的 python 代码如下:

 1 # -*- coding=utf-8 
 2 #!/usr/bin/python2.7
 3 
 4 from pylab import *
 5 
 6 
 7 def plotCircle(before, M=matrix([[1, 0], [0, 1]])):
 8     '''before: 变换前的矩阵
 9     M: 变换矩阵,默认为单位矩阵
10     返回变换之后的矩阵
11     '''
12     eclMat = M * before  # M 变换
13     eclX = array(eclMat[0]).reshape(-1)
14     eclY = array(eclMat[1]).reshape(-1)
15     axis('equal')
16     axis([-3, 3, -3, 3])
17     grid(True)
18     plot(eclX[:25], eclY[:25], 'c', linewidth=3)
19     plot(eclX[25:50], eclY[25:50], 'm', linewidth=3)
20     plot(eclX[50:75], eclY[50:75], 'y', linewidth=3)
21     plot(eclX[75:100], eclY[75:100], 'k', linewidth=3)
22     show()
23     return eclMat
24 
25 
26 ang = linspace(0, 2*pi, 100)
27 
28 x = cos(ang)
29 y = sin(ang)
30 cirMat = matrix([x, y]) # 2 × 100 的圆圈矩阵
31 
32 # 画最开始的图形——圆
33 plotCircle(cirMat)
34 
35 # 画变换之后的椭圆
36 M = matrix([[2, 1], [1, 2]]) # 2 × 2 的变换矩阵
37 clMat = plotCircle(cirMat, M)
38 
39 # 将 M 矩阵进行 svd 分解
40 U, s, V = np.linalg.svd(M, full_matrices=True)
41 S = np.diag(s)
42 
43 # SVD 矩阵对圆的变换
44 plotCircle(cirMat)
45 Tran1 = plotCircle(cirMat, V)
46 Tran2 = plotCircle(Tran1, S)
47 Tran3 = plotCircle(Tran2, U)

 

3. SVD 向量空间

假设矩阵 A 的秩是 r, 那么对角线矩阵的秩也是 r (乘以酉矩阵不会改变矩阵的秩), 我们假设:

  1. 那么, Ax = 0 的 null space 也就是解空间是什么呢?答案如下:


    证明非常简单,直接带入 SVD 分解三个矩阵的右边,计算的时候正交的向量相乘都等于0, 不为 0 的恰好都被对角线矩阵的 0 元素归零,等号成立。如果 A 是列满秩的,显然符合条件的只剩下零向量了。同样的,你可以知道 A' 的 null space。

     

  2. 矩阵 A 的线性子空间是啥?设想有一个很高的矩阵 A, 它的列向量很有可能不是正交的。对于任意一个坐标 x,矩阵 A 的线性子空间可以定义为:

    R(A) = {y | y = Ax, x 是任意的坐标}

    那么, u1, u2,..., ur 是 R(A) 的正交基。
    这个想法也非常直观,无论来了一个什么向量(或者叫坐标), 经过 V^T 和 对角线矩阵变换之后还是一个坐标,这个坐标就是 U 矩阵列向量线性组合的系数。

     

4. SVD 的计算

这就是正定矩阵的对角化。计算过程如下:

  1. 计算 A' (A 的转置)和 A'A
  2. 计算 A'A 的特征值,将特征值按照递减的顺序排列, 求均方根,得到 A 的奇异值
  3. 由奇异值可以构建出对角线矩阵 S, 同时求出 S 逆,以备后面的计算
  4. 有上述排好序的特征值可以求出对应的特征向量,以特征向量为列得到矩阵 V, 转置后得到 V'
  5. U = AVS逆,求出 U,完毕。

Lanczos algorithm 是一种迭代的计算方法,没来得及细看。

感觉 SVD 计算水很深,要用到的时候再看,现在暂不深入了。

 

5. PCA

主成分分析经常用于减少数据集的维数,同时保持数据集中的对方差贡献最大的特征。假设上面的椭圆中(二维空间,两个坐标值)均匀分布了非常多的点,如何在一维空间(一个坐标值)里面就能最大程度将这些点区分开来呢?这时候就要用到 PCA:

可以看到,上图中将所有的点投影到 +45° 直线上,将二维空间映射到一维空间,可以最大程度地区分开这些点,即投影后的样本分布方差最大,这个方向就是 v1 向量的方向。就上图来说,假设 100×2 的矩阵 X 表示样本点在平面上的坐标,将这些点投影到 v1 上保留最大的样本方差:

这样,就将 100 个 2 维空间的样本点压缩为 100 个 1 维空间的样本点,这里的列压缩实际上是对特征进行了压缩, 而不是简单地丢弃。PCA 是不是只能在样本点个数大于特征个数的时候才能压缩呢?例如,现在 3 个 100 维特征的样本用 100×3 的矩阵表示,用 SVD 分解后可以得到三个矩阵的乘积,做如下变换:

同样也最大限度保留了主成分。

总结起来,一个矩阵 A, 如果想对行进行压缩并保留主成分,那么左乘 u1', 如果相对列进行压缩并保留主成分(让我联想起了稀疏表示),那么右乘 v1。

当然,上面是简单的保留第一个主成分,PCA的全部工作简单点说,就是对原始的空间中顺序地找一组相互正交的坐标轴,第一个轴是使得方差最大的,第二个轴是在与第一个轴正交的平面中使得方差最大的,第三个轴是在与第1、2个轴正交的平面中方差最大的,这样假设在N维空间中,我们可以找到N个这样的坐标轴,我们取前r个去近似这个空间,这样就从一个N维的空间压缩到r维的空间了,但是我们选择的r个坐标轴能够使得空间的压缩使得数据的损失最小。

 

6. 最小二乘问题

想法和 QR 分解的办法类似,主要是利用酉矩阵变换的长度不变性:

得到最优解是:

掰掰~

posted on 2013-07-26 21:13  daniel-D  阅读(9301)  评论(1编辑  收藏  举报

导航