【记录读论文时遇到的一些算法5】—— SVD分解

1.特征值、特征向量

如果一个向量 $x$ 是 $n\times n$ 矩阵 $\mathbf{A}$ 的特征向量，那么可表示成以下形式：

$$\mathbf{A}x=\lambda x$$

其中，$\mathbf{A}$ 是一个 $n\times n$ 的实对称矩阵， $x$ 是一个 $n$ 维特征向量，$\lambda$ 是矩阵$\mathbf{A}$ 的特征值，$x$ 是矩阵 $\mathbf{A}$ 的特征值 $\lambda$ 所对应的特征向量。

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

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

线性变换

矩阵是线性空间里变换的描述。一个矩阵乘以一个向量，实质上是对向量做线性变换。对于一个对称矩阵$\mathbf{M}$：

$$\mathbf{M}=\left[\begin{array}{ll}3& 0\\ 0& 1\end{array}\right]$$
对应的线性变换是下面的形式:

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

它所描述的变换是下面的样子：

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

2. 特征分解

对于一个矩阵 $\mathbf{A}$，将其特征分解，得到矩阵 $\mathbf{A}$ 的 $n$ 个特征值 ${\lambda }_1\le {\lambda }_2\le \dots \le {\lambda }_n$ 以及对应的特征向量 { w 1 , w 2 , … , w n } \left\{w_{1}, w_{2}, \ldots, w_{n}\right\} {w1​,w2​,…,wn​} ，那么特征分解可 表示为:
$$A=W\Sigma W^{-1}$$
其中， $W$ 是矩阵$\mathbf{A}$的特征向量所构成的 $n\times n$ 维矩阵， $\Sigma$ 是 $n$ 个特征值为主对角线的 $n\times n$ 维对角矩阵。
通常情况下，将得到的一组特征向量进行Schmidt正交化单位化，即 ${\Vert w_i\Vert }_2=1$ ，那么这组特征向量为标准正交基，满足 $W^{T}W=I$ ， 即 $W^T=W^{-1}$ ，也就是说W为酉矩阵，表达式为:
$$A=W\Sigma W^T$$
特征分解的式子，分解得到的 $\Sigma$ 矩阵是一个对角阵，里面的特征值是由大到小排列的，这些特征值所对应的特征向量就是描述这个矩阵变 化方向 (从主要的变化到次要的变化排列)。
当矩阵是高维的情况下，那么这个矩阵就是高维空间下的一个线性变换，这个变换有很多的变换方向，我们通过特征值分解得到的前 $\mathrm{N}$ 个 特征向量，那么就对应了这个矩阵最主要的 $\mathrm{N}$ 个变化方向。我们利用这前 $\mathrm{N}$ 个变化方向，就可以近似这个矩阵（变换）。也就是之前说 的：提取这个矩阵最重要的特征。

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

特征分解的局限性

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

3. 奇异值分解（SVD）

SVD是一种适用于任意矩阵(表述不确切，有特殊情况)的分解方法，不限于方阵。对于一个 $m\times n$ 矩阵 $\mathbf{A}$，那么矩阵 $\mathbf{A}$ 的奇异值分解为:
$$\mathbf{A}=U\Sigma V^T$$
其中， $U$ 是一个 $m\times m$ 矩阵，矩阵 $U$ 中的正交向量为左奇异向量; $\mathrm{V}$ 是一个 $n\times n$ 矩阵，矩阵 $V$ 中的正交向量为右奇异向量， $U$ 和 $V$ 都是酉 矩阵，满足 $U^{T}U=I,V^{T}V=I;\Sigma$ 是一个 $m\times n$ 矩阵，除了对角线元素以外都为 0 ，对角线上的元素为奇异值。下图表示SVD的分解 过程:

思考：任意矩阵可通过SVD进行分解，那么如何求解 $U、V$ 和 $\Sigma$ 呢?

如果用矩阵$\mathbf{A}$乘以$\mathbf{A}$的转置得到一个 $m\times m$ 的方阵 $\mathbf{A}{\mathbf{A}}^T$ ，方阵进行特征分解，得到特征值以对应的特征向量:

$$\left(\mathbf{A}{\mathbf{A}}^T\right)u_i={\lambda }_i{u}_i$$

方阵 $\mathbf{A}{\mathbf{A}}^T$ 分解得到 $m$ 个特征值和对应的特征向量 $u$ ，将所有特征向量张成一个矩阵 $U$ ，就是SVD分解公式的 $U$ 矩阵。

同理，如果用矩阵$\mathbf{A}$的转置乘以$\mathbf{A}$得到一个 $n\times n$ 的方阵 $A^{T}A$ ，方阵进行特征分解，得到特征值以对应的特征向量:
$$\left({\mathbf{A}}^T\mathbf{A}\right)v_i={\lambda }_i{v}_i$$
方阵 $A{A}^T$ 分解得到 $n$ 个特征值和对应的特征向量 $v$ ，将所有特征向量张成一个矩阵 $V$ ，就是SVD分解公式的 $V$ 矩阵。

思考: ${\mathbf{A}}^T\mathbf{A}$ 的特征向量组成的矩阵是SVD中的$V$矩阵，而 $\mathbf{A}{\mathbf{A}}^T$ 的特征向量组成的矩阵是SVD中的$U$矩阵，这个是怎么证明的?

证明:
$$\begin{array}{c}\mathbf{A}=U\Sigma V^T\\ {\mathbf{A}}^T=V{\Sigma }^T{U}^T\\ \Rightarrow {\mathbf{A}}^T\mathbf{A}=V{\Sigma }^T{U}^{T}U\Sigma V^T=V{\Sigma }^2{V}^T\end{array}$$
上式使用了 $U^{T}U=I,{\Sigma }^T\Sigma ={\Sigma }^2$ 。不难看出 ${\mathbf{A}}^T\mathbf{A}$ 的特征向量组成的矩阵就是SVD中的 $V$ 矩阵。同 理， $\mathbf{A}{\mathbf{A}}^T$ 的特征向量组成的矩阵就是SVD中的 $U$ 矩阵。

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

1. 通过上式证明可发现， ${\mathbf{A}}^T\mathbf{A}$ 的特征值矩阵是奇异值矩阵的平方，也就是说特征值和奇异值满足 如下关系：
   $${\sigma }_i=\sqrt{{\lambda }_i}$$
2. $\mathbf{A}=U\Sigma V^T\Rightarrow \mathbf{A}V=U\Sigma V^{T}V\Rightarrow \mathbf{A}V=U\Sigma \Rightarrow \mathbf{A}{v}_i=\sigma u_i\Rightarrow {\sigma }_i=\frac{\mathbf{A}{v}_i}{u_i}$

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

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

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

4. SVD分解意义

SVD 分解可以看成先旋转，然后进行分别缩放，然后再旋转的过程。

5. SVD举例

$$A=\left(\begin{array}{ll}0& 1\\ 1& 1\\ 1& 0\end{array}\right)$$
先计算 $A^{T}A$ 和 $A{A}^T$
$$\begin{array}{c}{A}^{T}A=\left(\begin{array}{lll}0& 1& 1\\ 1& 1& 0\end{array}\right)\left(\begin{array}{ll}0& 1\\ 1& 1\\ 1& 0\end{array}\right)=\left(\begin{array}{ll}2& 1\\ 1& 2\end{array}\right)\\ A{A}^T=\left(\begin{array}{ll}0& 1\\ 1& 1\\ 1& 0\end{array}\right)\left(\begin{array}{lll}0& 1& 1\\ 1& 1& 0\end{array}\right)=\left(\begin{array}{lll}1& 1& 0\\ 1& 2& 1\\ 0& 1& 1\end{array}\right)\end{array}$$
然后求解 $A^{T}A$ 的特征值及对应的特征向量:
$${\lambda }_1=3;v_1=\left(\begin{array}{l}1/\sqrt{2}\\ 1/\sqrt{2}\end{array}\right);{\lambda }_2=1;v_2=\left(\begin{array}{c}-1/\sqrt{2}\\ 1/\sqrt{2}\end{array}\right)$$
同理求解 $A{A}^T$ 的特征值及对应的特征向量:
$${\lambda }_1=3;u_1=\left(\begin{array}{c}1/\sqrt{6}\\ 2/\sqrt{6}\\ 1/\sqrt{6}\end{array}\right);{\lambda }_2=1;u_2=\left(\begin{array}{c}1/\sqrt{2}\\ 0\\ -1/\sqrt{2}\end{array}\right);{\lambda }_3=0;u_3=\left(\begin{array}{c}1/\sqrt{3}\\ -1/\sqrt{3}\\ 1/\sqrt{3}\end{array}\right)$$
通过 ${\sigma }_i=\sqrt{{\lambda }_i}$ 求解奇异值为 $\sqrt{3}$ 和1
最终矩阵A的奇异值分解为:
$$A=U\Sigma V^T=\left(\begin{array}{ccc}1/\sqrt{6}& 1/\sqrt{2}& 1/\sqrt{3}\\ 2/\sqrt{6}& 0& -1/\sqrt{3}\\ 1/\sqrt{6}& -1/\sqrt{2}& 1/\sqrt{3}\end{array}\right)\left(\begin{array}{cc}\sqrt{3}& 0\\ 0& 1\\ 0& 0\end{array}\right)\left(\begin{array}{cc}1/\sqrt{2}& 1/\sqrt{2}\\ -1/\sqrt{2}& 1/\sqrt{2}\end{array}\right)$$

5. SVD分解的应用

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

6.为什么Ax=0的解为最小奇异值对应的向量？

工程中很多问题会归结为求超定方程 $\mathbf{A}\mathbf{x}=0，\mathbf{A}$ 是 $m\times n$的矩阵，且 $m>n$ 。如 SLAM中三角化地图点，PnP等一些问题都是求解这个方程。
很显然，这个方程有一个解，但这不是我们想要的，我们实际想求非零解。

为了求非零解，我们对 $\mathbf{x}$ 加上一个约束 $\Vert \mathbf{x}{\Vert }^2=1$ 。也就是限制 $\mathbf{x}$ 的长度为 1 。并构建成一 个带约束的最小二乘问题:
$$\hat{\mathbf{x}}=\arg \min \Vert \mathbf{A}\mathbf{x}{\Vert }^2\text{, subject to }\Vert \mathbf{x}{\Vert }^2=1\text{\hspace{1em}(1)}$$
这是一个带约束的最小二乘问题，我们把拉格朗日搬出来:
$$\begin{array}{rl}L(\mathbf{x},\lambda )& =\Vert \mathbf{A}\mathbf{x}{\Vert }^2+\lambda \left(1-\Vert \mathbf{x}{\Vert }^2\right)\\ & ={\mathbf{x}}^T{\mathbf{A}}^T\mathbf{A}\mathbf{x}+\lambda \left(1-{\mathbf{x}}^T\mathbf{x}\right)\end{array}\text{\hspace{1em}(2)}$$
为了求极值，我们分别对 $\mathbf{x}$ 和 $\lambda$ 求偏导数，令为 0 :
$$\begin{array}{rl}& \frac{\partial L(\mathbf{x},\lambda )}{\partial \mathbf{x}}=2{\mathbf{A}}^T\mathbf{A}\mathbf{x}-2\lambda \mathbf{x}=0\\ & \frac{\partial L(\mathbf{x},\lambda )}{\partial \lambda }=1-{\mathbf{x}}^T\mathbf{x}=0\end{array}\text{\hspace{1em}(3)}$$
把(3)式整理一下:
$$\begin{array}{r}\left({\mathbf{A}}^T\mathbf{A}-\lambda \mathbf{I}\right)\mathbf{x}=0\\ {\mathbf{A}}^T\mathbf{A}\mathbf{x}=\lambda \mathbf{x}\end{array}\text{\hspace{1em}(4)}$$
可以看出 $\lambda$ 和 $\mathbf{x}$ 分别是 ${\mathbf{A}}^T\mathbf{A}$ 的特征值和特征向量。也就是说(1)式的解，就是这些特征向量 中的一个。
问题来了，那么多的特征向量，应该选择哪个作为解呢? 我们展开 $\Vert \mathbf{A}\mathbf{x}{\Vert }^2$ 看一下:
$$\Vert \mathbf{A}\mathbf{x}{\Vert }^2={\mathbf{x}}^T{\mathbf{A}}^T\mathbf{A}\mathbf{x}={\mathbf{x}}^T\lambda \mathbf{x}=\lambda {\mathbf{x}}^T\mathbf{x}=\lambda \text{\hspace{1em}(5)}$$

（5）式的推导，用到了式(4)和 $\Vert \mathbf{x}{\Vert }^2=1$ 。

也就是说，我们想要$|\mathbf{A}\mathbf{x}{\Vert }^2$最小，就需要$\lambda$最小。

那方程 (1) 的非零解就是 ${\mathbf{A}}^T\mathbf{A}$ 最小特征值 $\lambda$ 对应的特征向量。

例子：

假设$\mathbf{A}$为8*9的矩阵，则SVD分解结果为:
$$\mathbf{A}=UD{V}^T$$

U：左奇异向量，为8×8的正交矩阵
V：右奇异向量，为9×9的正交矩阵

D：一个8*9的对角矩阵，除了对角线元素均为0，对角线元素称为奇异值，一般来说奇异值是按照从大到小的顺序降序排列。因为每一个奇异值都是一个残差项，因此最后一个奇异值最小，其含义是最优的残差。因此其对用的奇异值向量就是最优解。

参考文献

[1] https://blog.csdn.net/wangshuailpp/article/details/80209863 视觉SLAM常见的QR分解SVD分解等矩阵分解方式求解满秩和亏秩最小二乘问题（最全的方法分析总结）
[2] https://blog.csdn.net/code__online/article/details/90705000 矩阵分解SVD原理
[3] https://blog.csdn.net/jdy_lyy/article/details/104828112 矩阵（一）：SVD分解
[4] https://blog.csdn.net/jdy_lyy/article/details/117171934 矩阵（二）：为什么Ax=0的解为最小奇异值对应的向量？