奇异值分解 -- singular value decomposition (SVD)

摘要

在本文中我们首先介绍什么是奇异值分解(singular value decomposition, SVD)。我们会给出 SVD 的公式，以及其在Python中的实现。我们也会讨论 SVD 与 主成分分析 (principal component analysis, PCA) 之间的关系。

什么是奇异值分解

SVD 是对矩阵的一种分解方法。对任意的矩阵，我们都可以对其进行 SVD 分解。具体来说，对于任意一个 $m\times n$ 的矩阵 $X$，我们可以将 $X$ 分解为
$$X=U\Sigma V^T\text{\hspace{1em}(1)}$$
$U\in {\mathbb{R}}^{m\times m},\Sigma \in {\mathbb{R}}^{m\times n},V\in {\mathbb{R}}^{n\times n}$ 。

其中 $U^{T}U=I_{m\times m},V^{T}V=I_{n\times n}$， $\Sigma$ 是一个对角矩阵，
$$\Sigma =\left(\begin{array}{cccccc}{\sigma }_1& & & & & \\ & \ddots & & & \text{0}& \\ & & {\sigma }_r& & & \\ & & & 0& & \\ & \text{0}& & & \ddots & \\ & & & & & 0\end{array}\right)$$
这里 ${\sigma }_1\ge {\sigma }_2\cdots \ge {\sigma }_r\ge 0$。$r$ 是矩阵 $X^{T}X$ 的秩 (rank)。 并且 ${\sigma }_i$ 为矩阵 $X^{T}X$ 的特征值 ${\lambda }_i$ 的开方。即 ${\lambda }_i={\sigma }_i^2$。

而 $V$ 是由矩阵 $X^{T}X$ 的特征向量为列组成的矩阵。即 $V=(v_1,v_2,\cdots ,v_n)$，每一个 $v_i$ 是一个 $m\times 1$ 的列向量， $X^{T}X{v}_i={\lambda }_i{v}_i$ 。这里我们的 $v_i$ 是归一化的，即每个 $v_i$ 的模为1。所以 $V$ 的列向量就是 ${\mathbb{R}}^n$ 空间的一组归一化的正交基向量 (orthonormal basis)。
$$v_i^T{v}_j={\delta }_{ij}=\{\begin{array}{l}1,i=j\\ 0,i\ne j\end{array}$$。
值得注意的是，如果 $r<n$，即如果我们只有 $r$ 个非零的特征值，我们须要补齐 $n-r$ 个单位向量，来与原来的 $r$ 个单位向量（特征向量）组成 $n$ 个向量 [1]。

我们定义 $u_i=\frac{1}{{\sigma }_i}X{v}_i$，那么每个 $u_i$ 均为 $n\times 1$ 的列向量。因为 $X^{T}X{v}_i={\lambda }_i{v}_i\ne 0$，故 $X{v}_i\ne 0$。
并且，我们有
$$u_i^T{u}_j=\frac{1}{{\sigma }_i{\sigma }_j}{v}_i^T{X}^{T}X{v}_j=\frac{1}{{\sigma }_i{\sigma }_j}{v}_i^T{\lambda }_j{v}_j={\delta }_{ij}$$。

我们就验证了 { u 1 ,   u 2 ,   ⋯   ,   u r } \{ u_1, \, u_2, \, \cdots, \, u_r \} {u1​,u2​,⋯,ur​} 是 $r$ 维空间的一组基。同样的，如果 $r<m$，我们可以将这个 $r$ 维空间的基向量扩展成 $m$ 维空间的基向量。记为 { u 1 ,   u 2 , ⋯   ,   u m } \{u_1, \, u_2, \cdots, \, u_m \} {u1​,u2​,⋯,um​}，那么以 $u_1,u_2,\cdots ,u_m$ 为列向量组成的矩阵就为 $U$。

这里须要指出的是，矩阵 $X^{T}X$ 的秩与矩阵 $X$ (或 $X^T$)的秩相同，即 $rank(X^{T}X)=rank(X)=rank(X^T)$。因为 $rank(X)\le \min (m,n)$，所以 $rank(X^{T}X)\le \min (m,n)$。从而我们无须担心由于 $r$ 过大而无法补全基向量的情况。关于 $rank(X^{T}X)=rank(X)$ 的证明可见附录。

由于 $rank(X^{T}X)=rank(X)$，如果 $rank(X^{T}X)=r<n$，那么当我们用 $r$ 个归一化的 $X^{T}X$ 的特征向量构造 $V$ 之后，还要补齐 $n-r$ 个模为 1 的基向量。对于这 $n-r$ 个补齐的基向量 $v_k,r+1\le k\le n$，我们有 $X{v}_k=0,r+1\le k\le n$。这是因为 $v_1,v_2,\cdots ,v_r$ 都不在 $X$ 的零空间中，并且 $v_1,v_2,\cdots ,v_r$ 形成了一个 $r$ 维空间。而 $v_{r+1},\cdots ,v_n$ 这 $n-r$ 个基向量均与 $v_1,v_2,\cdots ,v_r$ 垂直，所以 $v_{r+1},\cdots ,v_n$ 都在 $X$ 的零空间中。

接着上面的分析，如果 $r<n$，我们补齐的基向量 $v_{k+1},\cdots ,v_n$ 均在 $X$ 的零空间中，我们便不能再取 $X{v}_j$ 当作 $u_j$ （如果 $U$ 也须要补齐的话）。那这个时候如何选取 $u_j$ 向量呢？我们要做的很简单，这时候只须要根据已有的 $u_j,1\le j\le r\le m$，来构造出 $m$ 维空间中剩下的 $m-r$ 个基向量，来当作 $U$ 的剩下的 $m-r$ 个列向量即可。

左右奇异向量

根据上面的定义，我们称 $V$ 的列向量 $v_i$ 是 $X$ 的右奇异向量 （right singular vectors）；称 $U$ 的列向量是 $X$ 的左奇异向量（left singular vectors）。 我们有 $v_i$ 是矩阵 $X^{T}X$ 的特征向量；$u_i$ 是矩阵 $X{X}^T$ 的特征向量。

并且，如果 $X^{T}X{v}_i={\lambda }_i{v}_i,{\lambda }_i>0$，那么 $u_i=\frac{1}{\sqrt{{\lambda }_i}}X{v}_i$ 便是 $X{X}^T$ 的特征向量，且有 $X{X}^T{u}_i={\lambda }_i{u}_i$。同样的，如果 $u_i$ 是矩阵 $X{X}^T$ 的特征向量，$X{X}^T{u}_i={\lambda }_i{u}_i,{\lambda }_i>0$，那么 $v_i=\frac{1}{\sqrt{{\lambda }_i}}{X}^T{u}_i$ 便是矩阵 $X^{T}X$ 的特征向量，且有 $X^{T}X{v}_i={\lambda }_i{v}_i$。

通过上述论述，我们同时也证明了矩阵 $X^{T}X$ 和 矩阵 $X{X}^T$ 的非零特征根是一样的。

程序实现

在 Python 中，我们可以用 numpy 中 linalg.svd 来对矩阵进行 SVD 分解。

import numpy as np
from numpy.linalg import svd
from sklearn.decomposition import TruncatedSVD

X = np.array([[1, 2, 3], [2, 4, 2]])
svd(X, full_matrices=True)

(array([[-0.59233648, -0.8056907 ],
[-0.8056907 , 0.59233648]]),
array([5.98022195, 1.49564213]),
array([[-3.68501017e-01, -7.37002033e-01, -5.66599509e-01],
[ 2.53391004e-01, 5.06782008e-01, -8.23993323e-01],
[-8.94427191e-01, 4.47213595e-01, -2.22044605e-16]]))

不确定性

根据之前的分析， 如果 $r<n$，我们在确定 $V$ 的时候须要对 $n$ 维空间的基向量进行补齐。而这种补齐并不是唯一的。举例来说，在四维空间中，如果我们已经有了两个标准基向量
$$u_1=\left(\begin{array}{c}1\\ 0\\ 0\\ 0\end{array}\right),u_2=\left(\begin{array}{c}0\\ 1\\ 0\\ 0\end{array}\right)$$，
那么我们再补齐的两个归一化基向量可以选择
$$u_3=\left(\begin{array}{c}0\\ 0\\ 1\\ 0\end{array}\right),u_4=\left(\begin{array}{c}0\\ 0\\ 0\\ 1\end{array}\right)$$，
也可以选择
$$u_3^{\prime }=\left(\begin{array}{c}0\\ 0\\ \frac{\sqrt{2}}{2}\\ \frac{\sqrt{2}}{2}\end{array}\right),u_4^{\prime }=\left(\begin{array}{c}0\\ 0\\ \frac{\sqrt{2}}{2}\\ -\frac{\sqrt{2}}{2}\end{array}\right)$$。

那么对于 SVD，如果 $r<n$，分解也是不确定的。
比如说对于 $X=\left(\begin{array}{ccc}1& 1& 1\\ 2& 2& 2\end{array}\right)$，根据 python 的结果我们有

X = np.array([[1, 1, 1], [2, 2, 2]])
svd(X, full_matrices=True)

(array([[-0.4472136 , -0.89442719],
[-0.89442719, 0.4472136 ]]),
array([3.87298335e+00, 1.40433339e-16]),
array([[-0.57735027, -0.57735027, -0.57735027],
[-0.81649658, 0.40824829, 0.40824829],
[ 0. , -0.70710678, 0.70710678]]))

即 $U=\left(\begin{array}{cc}-0.4472& -0.8944\\ -0.8944& 0.4472\end{array}\right)$，
$\Sigma =\left(\begin{array}{ccc}\sqrt{15}& 0& 0\\ 0& 0& 0\end{array}\right)$，
$V=\left(\begin{array}{ccc}-0.5774& -0.8164& 0\\ -0.5774& 0.4082& -0.7071\\ -0.5774& 0.4082& 0.7071\end{array}\right)$。

我们也可以选择 $V^{\prime }=\left(\begin{array}{ccc}-0.5774& \frac{1+\sqrt{5}}{4}& -0.1102\\ -0.5774& \frac{1-\sqrt{5}}{4}& 0.7558\\ -0.5774& -0.5& -0.6455\end{array}\right)$。
可以验证，这时 $U\Sigma V^{\prime T}$ 依然等于 $X$ 。

与 PCA 的关系

对于熟悉 PCA [3] 的读者，可以发现 SVD 与 PCA本质是一样的。我们求数据集 $X$ 的主成分，实际是对 $X$ 进行 SVD 分解，求解 $X^{T}X$ 的特征向量 （这里假设 $X$ 的列向量的均值为0）。

附录

假设 $X\in {\mathbb{R}}^{m\times n}$，要证明 $rank(X^{T}X)=rank(X)$，我们可以从零空间的角度 [2] 去考虑。如果我们能证明矩阵 $X^{T}X$ 的零空间和矩阵 $X$ 的零空间是相同的，那么我们就证明了 $rank(X^{T}X)=rank(X)$。这是因为矩阵 $X^{T}X$ 和矩阵 $X$ 都有 $n$ 列，所以
$$rank(X)+dim(null(X))=n$$ $$rank(X^{T}X)+dim(null(X^{T}X))=n$$ 。
这里 $dim(null(X))$ 表示矩阵 $X$ 的零空间的维度。如果 $null(X)=null(X^{T}X)$，即矩阵 $X$ 与矩阵 $X^{T}X$ 的零空间相同，那么它们的维度自然也相同。下面我们就证明矩阵 $X$ 与矩阵 $X^{T}X$ 的零空间相同。

如果 $Xw=0$，那么自然有 $X^{T}Xw=0$。也就是说每一个 $null(X)$ 的向量都在 $null(X^{T}X)$ 之中。所以 $null(X)\subset null(X^{T}X)$。
另一方面，如果 $X^{T}Xw=0$，那么 $w^T{X}^{T}Xw=0$。所以 $||Xw|{|}^2=0$，即 $Xw=0$。从而 $null(X^{T}X)\subset null(X)$。于是 $null(X)=null(X^{T}X)$。
$\square$

参考文献

[1] A Tutorial on Principal Component Analysis, Jonathon Shlens, arXiv, 2014
[2] 机器学习：矩阵的秩和矩阵的四个子空间, Matrix_11, CSDN blog
[3] 主成分分析（principal component analysis, PCA）公式, kdaHugh, CSDN blog