矩阵分解SVD原理

常用的经典矩阵分解算法:
矩阵分解SVD原理
  • 经典算法PCA、SVD
  • 主题模型算法LDA
  • 概率矩阵分解PMF,由深度学习大牛Ruslan Salakhutdinov所写,主要应用于推荐系统中,在大规模的稀疏不平衡性Netflix数据集上取得较好的效果;
  • 非负矩阵分解(Non-negative Matrix Factorization)NMF,由Lee和 Seung 在《Nature》上提出,应用于文本聚类;

1. 特征值、特征向量

如果一个向量xxn×nn\times n矩阵A的特征向量,那么可表示成以下形式:Ax=λxAx=\lambda x其中,A是一个n×nn \times n的实对称矩阵,xx是一个n维特征向量,λ\lambda是矩阵A的特征值,xx是矩阵A的特征值λ\lambda所对应的特征向量。

思考:为什么一个矩阵乘以一个向量的效果与一个实数乘以相同向量的效果是一样的呢?

事实上,矩阵是线性空间里变换的描述。矩阵A与向量相乘,本质上对向量xx进行一次线性变换(旋转或者拉伸变换),而该转换的效果等价于常数λ\lambda乘以向量xx(拉伸)的效果。所以,当求解矩阵的特征值与对应的特征向量时,就是为了求矩阵A能使得哪些向量只发生拉伸变换,而拉伸的程度用特征值λ\lambda来度量。

线性变换

矩阵是线性空间里变换的描述。一个矩阵乘以一个向量,实质上是对向量做线性变换。对于一个对称矩阵M:M=[3001]M=\begin{bmatrix} 3 & 0 \\ 0 & 1 \end{bmatrix}
对应的线性变换是下面的形式:

矩阵分解SVD原理

因为这个矩阵M乘以一个向量(x,y)的结果是:[3001][xy]=[3xy]\begin{bmatrix} 3 & 0 \\ 0 & 1 \end{bmatrix}\begin{bmatrix} x \\ y \end{bmatrix}=\begin{bmatrix} 3x \\ y \end{bmatrix}上面的矩阵是对称的,所以这个变换是一个对x,y轴的方向一个拉伸变换(每一个对角线上的元素将会对一个维度进行拉伸变换,当值大于1时,是拉长变换,当值小于1时是缩短变换),当矩阵不是对称的时候,假如说矩阵是下面的样子: [3001]\begin{bmatrix} 3 & 0 \\ 0 & 1 \end{bmatrix}
它所描述的变换是下面的样子:

矩阵分解SVD原理

这其实是在平面上对一个轴进行的拉伸变换(如蓝色的箭头所示),在图中,蓝色的箭头是一个最主要的变化方向(变化方向可能有不止一个),如果我们想要描述好一个变换,那我们就描述好这个变换主要的变化方向就好了。

2. 特征分解

对于一个矩阵A,将其特征分解,得到矩阵A的n个特征值λ1λ2...λn\lambda_1 \leq \lambda_2 \leq ...\leq \lambda_n以及对应的特征向量{w1,w2,...,wn}\left\{w_1,w_2,...,w_n\right\},那么特征分解可表示为:A=WΣW1A=W\Sigma W^{-1}其中,WW是矩阵A的特征向量所张成的n×nn\times n维矩阵,Σ\Sigma是n个特征值为主对角线的n×nn\times n维对角矩阵。
通常情况下,将得到的一组特征向量进行Schmidt正交化单位化,即wi2=1||w_i||_2=1,那么这组特征向量为标准正交基,满足WTW=IW^TW=I,即WT=W1W^{T}=W^{-1},也就是说W为酉矩阵,表达式为:A=WΣWTA=W\Sigma W^T

特征分解的式子,分解得到的Σ矩阵是一个对角阵,里面的特征值是由大到小排列的,这些特征值所对应的特征向量就是描述这个矩阵变化方向(从主要的变化到次要的变化排列)。

当矩阵是高维的情况下,那么这个矩阵就是高维空间下的一个线性变换,这个变换有很多的变换方向,我们通过特征值分解得到的前N个特征向量,那么就对应了这个矩阵最主要的N个变化方向。我们利用这前N个变化方向,就可以近似这个矩阵(变换)。也就是之前说的:提取这个矩阵最重要的特征

总结一下,矩阵特征分解可以得到特征值与特征向量,特征值表示的是这个特征到底有多重要,而特征向量表示这个特征是什么

特征分解的局限性

特征分解有一些局限性,比如变换矩阵必须是方阵,即nnn * n的矩阵。而在实际应用场景中,大部分不是这种矩阵。举个最简单的例子,关系型数据库中的某一张表的数据存储结构就类似于一个二维矩阵,假设这个表有m行,有n个字段,那么这个表数据矩阵的规模就是mnm∗n。很明显,在绝大部分情况下,m与n并不相等。如果对这个矩阵要进行特征提取,特征值分解的方法显然就行不通了。那么这个时候就轮到SVD登场。

3. 奇异值分解(SVD)

SVD是一种适用于任意矩阵的分解方法,不限于方阵。对于一个m×nm\times n矩阵A,那么矩阵A的奇异值分解为:A=UΣVTA=U\Sigma V^T其中,U是一个m×mm\times m矩阵,矩阵U中的正交向量为左奇异向量;V是一个n×nn\times n矩阵,矩阵V中的正交向量为右奇异向量,U和V都是酉矩阵,满足UTU=I,VTV=IU^TU=I,V^TV=IΣ\Sigma是一个m×nm\times n矩阵,除了对角线元素以外都为0,对角线上的元素为奇异值。下图表示SVD的分解过程:

矩阵分解SVD原理

思考:任意矩阵可通过SVD进行分解,那么如何求解UUVVΣ\Sigma呢?

如果用矩阵A乘以A的转置得到一个m×mm \times m的方阵AATAA^T,方阵进行特征分解,得到特征值以对应的特征向量:(AAT)ui=λiui(AA^T)u_i=\lambda_iu_i方阵AATAA^T分解得到m个特征值和对应的特征向量uu,将所有特征向量张成一个矩阵U,就是SVD分解公式的U矩阵。
同理,如果用矩阵A的转置乘以A得到一个n×nn \times n的方阵ATAA^TA,方阵进行特征分解,得到特征值以对应的特征向量:(ATA)vi=λivi(A^TA)v_i=\lambda_iv_i方阵AATAA^T分解得到n个特征值和对应的特征向量vv,将所有特征向量张成一个矩阵V,就是SVD分解公式的V矩阵。

思考:ATAA^TA的特征向量组成的矩阵是SVD中的V矩阵,而AATAA^T的特征向量组成的矩阵是SVD中的U矩阵,这个是怎么证明的?
证明:
A=UΣVTA=U\Sigma V^TAT=VΣTUTA^T=V\Sigma^T U^TATA=VΣTUTUΣVT=VΣ2VT\Rightarrow A^TA=V\Sigma^TU^TU\Sigma V^T=V\Sigma^2V^T
上式使用了UTU=I,ΣTΣ=Σ2U^TU=I,\Sigma^T\Sigma=\Sigma^2。不难看出ATAA^TA的特征向量组成的矩阵就是SVD中的V矩阵。同理,AATAA^T的特征向量组成的矩阵就是SVD中的U矩阵。

接下来求解的是奇异值,其解法有两种:

  1. 通过上式证明可发现,ATAA^TA的特征值矩阵是奇异值矩阵的平方,也就是说特征值和奇异值满足如下关系:σi=λi\sigma_i=\sqrt{\lambda_i}
  2. A=UΣVTAV=UΣVTVAV=UΣAvi=σuiσi=AviuiA=U\Sigma V^T\Rightarrow AV=U\Sigma V^TV\Rightarrow AV=U\Sigma \Rightarrow Av_i=\sigma u_i \Rightarrow \sigma_i=\frac{Av_i}{u_i}

思考:任意矩阵可进行SVD分解,那么问题又来了,一个mnm*n的矩阵A,你把它分解成mmm*m的矩阵U、mnm*n的矩阵Σ和nnn*n的矩阵。这三个矩阵中任何一个的维度似乎一点也不比A的维度小,而且还要做两次矩阵的乘法,这不是没事找事干嘛!把简单的事情搞复杂了么!并且我们知道矩阵乘法的时间复杂度为。O那奇异值分解到底要怎么做呢?

回答:在奇异值分解矩阵中Σ里面的奇异值按从大到小的顺序排列,奇异值从大到小的顺序减小的特别快。在很多情况下,前10%甚至1%的奇异值的和就占了全部的奇异值之和的99%以上。也就是说,剩下的90%甚至99%的奇异值几乎没有什么作用。因此,我们可以用前面r个大的奇异值来近似描述矩阵,于是奇异值分解公式可以写成如下:AmnUmrΣrrVrnTA_{m*n}\approx U_{m*r}\Sigma_{r*r}V^T_{r*n}

矩阵分解SVD原理

其中,r<<m,r<<nr<<m,r<<n。将一个矩阵A分解为三个小矩阵。如果r越大,与原来的矩阵相似度越大,但存储和计算成本也会越大。因此,使用SVD时,需要根据不同的业务场景、资源情况来合理选择r的大小。本质上是在计算精度与空间时间成本之前做折中。

4. SVD举例

用一个简单的例子来详细化矩阵的奇异值分解过程。假如一个矩阵A为:A=(011110)A=\begin{pmatrix}0 & 1 \\ 1 & 1 \\ 1 & 0 \end{pmatrix}先计算ATAA^TAAATAA^TATA=(011110)(011110)=(2112)A^TA=\begin{pmatrix}0 & 1 & 1 \\ 1 & 1 & 0 \end{pmatrix}\begin{pmatrix}0 & 1 \\ 1 & 1 \\ 1 & 0 \end{pmatrix}=\begin{pmatrix}2 & 1 \\ 1 & 2 \end{pmatrix}AAT=(011110)(011110)=(110121011)AA^T=\begin{pmatrix}0 & 1 \\ 1 & 1 \\ 1 & 0 \end{pmatrix}\begin{pmatrix}0 & 1 & 1 \\ 1 & 1 & 0 \end{pmatrix}=\begin{pmatrix}1 & 1 & 0 \\ 1 & 2 & 1 \\ 0 & 1 & 1 \end{pmatrix}
然后求解ATAA^TA的特征值及对应的特征向量:λ1=3;v1=(1/21/2);λ2=1;v2=(1/21/2)\lambda_1=3;v_1=\begin{pmatrix}1/\sqrt2 \\ 1/\sqrt2\end{pmatrix};\lambda_2=1;v_2=\begin{pmatrix}-1/\sqrt2 \\ 1/\sqrt2\end{pmatrix}同理求解AATAA^T的特征值及对应的特征向量:λ1=3;u1=(1/62/61/6);λ2=1;u2=(1/201/2);λ3=0;u3=(1/31/31/3)\lambda_1=3;u_1=\begin{pmatrix}1/\sqrt6 \\ 2/\sqrt6 \\ 1/\sqrt6 \end{pmatrix};\lambda_2=1;u_2=\begin{pmatrix}1/\sqrt2 \\ 0 \\ -1/\sqrt2 \end{pmatrix};\lambda_3=0;u_3=\begin{pmatrix}1/\sqrt3 \\ -1/\sqrt3 \\ 1/\sqrt3 \end{pmatrix}通过σi=λi\sigma_i=\sqrt{\lambda_i}求解奇异值为3\sqrt3和1
最终矩阵A的奇异值分解为:A=UΣVT=(1/61/21/32/601/31/61/21/3)(300100)(1/21/21/21/2)A=U\Sigma V^T=\begin{pmatrix}1/\sqrt6 & 1/\sqrt2 & 1/\sqrt3 \\ 2/\sqrt6 & 0 & -1/\sqrt3 \\ 1/\sqrt6 & -1/\sqrt2 & 1/\sqrt3 \end{pmatrix}\begin{pmatrix}\sqrt3 & 0 \\ 0 & 1 \\ 0 & 0 \end{pmatrix}\begin{pmatrix}1/\sqrt2 & 1/\sqrt2 \\ -1/\sqrt2 & 1/\sqrt2 \end{pmatrix}

5. SVD分解的应用

  1. 降维
    矩阵A的特征有n维,经过SVD分解之后,完全可以用前r个非零奇异值对应的奇异向量表示矩阵A的主要特征。这样,就起到了蒋伟德作用。
  2. 压缩
    经过SVD分解之后,表示原来的矩阵A,只需要存U,Σ,VU,\Sigma ,V三个较小的矩阵即可。而这三个小矩阵的规模加起来也远远小于原始矩阵A。这样,就达到压缩的作用。
  3. PCA
    PCA降维需要找到样本协方差矩阵XTXX^TX的最大d个特征向量,然后用这些特征向量张成的矩阵来做低维投影降维。这个过程中,需要先求出协方差矩阵XTXX^TX,但是,当样本数和特征数很多的时候,计算量是相当大的。
    SVD可以应用于PCA降维。注意到SVD可得到协方差矩阵XTXX^TX最大的d个特征向量张成的矩阵,但SVD有一个好处是先不求协方差矩阵XTXX^TX,也能求解出右奇异矩阵V。也就是说,PCA算法可以不用做特征分解,而是做SVD来完成。这个方法在样本量很大的时候很有效。实际上,scikit-learn的PCA算法的背后真正实现就是SVD,而不是暴力特征求解。
    另一方面,PCA仅仅使用SVD的其中一个奇异矩阵,如右奇异矩阵,没有使用左奇异矩阵。假设样本是m×nm\times n矩阵X,如果使用SVD得到矩阵XXTXX^T最大的d个特征向量张成的m×dm\times d维矩阵U,则进行如下处理:Xd×n=Ud×mTXm×nX^{'}_{d\times n}=U^T_{d\times m}X_{m\times n}得到一个d×nd\times n的矩阵XX^{'},与原来的m×nm\times n维矩阵X相比,行数由m降到d,可见对行数进行了压缩。换句话说,左奇异矩阵用于行数的压缩;相对的,右奇异矩阵用于列数(特征维度)的压缩,也就是PCA降维。

参考文献