奇异值分解（SVD）

Jiaxxxxxx

已于 2023-04-30 22:33:54 修改

阅读量8.4k

点赞数 15

分类专栏：数学基础机器学习文章标签：线性代数算法机器学习

于 2021-08-29 20:16:41 首次发布

本文链接：https://blog.csdn.net/yyy_77/article/details/119981659

版权

机器学习同时被 2 个专栏收录

8 篇文章 1 订阅

订阅专栏

数学基础

7 篇文章 5 订阅

订阅专栏

奇异值分解SVD

0 什么是奇异值分解
1 特征值分解和几何意义
2 SVD分解
3 SVD代码实现
- 3.1 基于opencv实现SVD
- 3.2 基于eigen实现SVD
4 应用
5 最后，简单说下SVD和PCA之间的关系

0 什么是奇异值分解

奇异值分解（Singular-Value Decomposition，简称SVD）是一种矩阵分解算法，用于将矩阵简化为其组成部分，以使某些后续矩阵计算更简单。它和特征分解不同，SVD并不要求要分解的矩阵为方阵。那为什么将矩阵进行分解呢？先看以下几个例子：

例子1.以下是一个所有元素值都相等的矩阵

$A=\begin{pmatrix} 1 &1 &1 &1 &1 &1 \\ 1 &1 &1 &1 &1 & 1\\ 1 &1 &1 &1 &1 &1 \\ 1 &1 &1 &1 &1 &1 \\ 1 &1 &1 &1 &1 &1 \\ 1 &1 &1 &1 &1 &1 \end{pmatrix}$
如果要保存这个矩阵，通常可以压缩这个矩阵，压缩的方法有很多，此时如果用向量积的方式实现，如下所示：
$A=\begin{pmatrix} 1\\ 1\\ 1\\ 1\\ 1\\ 1 \end{pmatrix}\begin{pmatrix} 1 &1 &1 &1 &1 &1 \end{pmatrix}$
此时只需保存12个元素，而不是36个元素
例子2.意大利国旗
在这里插入图片描述
意大利国旗可以写成以下矩阵形式：
$A=\begin{pmatrix} g &g &w &w &r &r \\ g &g &w &w &r &r \\ g &g &w &w &r &r \\ g &g &w &w &r &r \\ g &g &w &w &r &r \\ g &g &w &w &r &r \end{pmatrix}$
同样的，也能用一种非常好的方式来压缩它，如下：
$A=\begin{pmatrix} 1\\ 1\\ 1\\ 1\\ 1\\ 1 \end{pmatrix}\begin{pmatrix} g & g & w &w &r &r \end{pmatrix}$

例子3:现有如下矩阵A
$A=\begin{pmatrix} 1 &0 \\ 1 &1 \end{pmatrix}$
分解这样的矩阵可以考虑用以上的方法，显然，不能只用一对行列向量完成分解。可以考虑将矩阵分解为两项的和（差），如下：

$A=\begin{pmatrix} 1\\ 1 \end{pmatrix}\begin{pmatrix} 1 &1 \end{pmatrix}-\begin{pmatrix} 1\\ 0 \end{pmatrix}\begin{pmatrix} 0 &1 \end{pmatrix}$
从以上的例子可以看出，SVD中数据压缩上起着很重要的作用。比如一个典型的机器学习问题可能有着几百个或者更多的变量，而许多机器学习算法如果出现超过几十个就会崩溃，这使得奇异值分解在ML中对于变量减少是必不可少的。除此以外，
那么奇异值分解实际上在做什么？
SVD将矩阵分解为一组向量v（可以看成多个列向量）、一组向量u（可以看成多个行向量）和一个对角矩阵（可以看成多个标量），即将矩阵分解为多项式。
如果有一个秩r=2的矩阵，我们可以将矩阵分解为
$u_{1}v_{1}^{T}+u_{2}v_{2}^{T}$
如果有一个秩r=1的矩阵，结果如下
$u_{1}v_{1}^{T}$
稍微复杂一点的分解是添加标量sigma存储在对角矩阵中，对于秩r=2的矩阵，分解结果如下
$\sigma _{1}u_{1}v_{1}^{T}+\sigma _{2}u_{2}v_{2}^{T}$
显然，sigma值可以认为是加权系数。

1 特征值分解和几何意义

再理解SVD之前，先了解矩阵的特征值分解。
对于某一个矩阵A可以看作是通过乘法作用于向量x以产生新向量Ax的变换，如果A是m x p矩阵，B是p x n矩阵，则矩阵乘积C=AB(即m x n矩阵)定义为：
$C_{ij} = \sum_{k=1}^{p}a_{ik}b_{kj}$
例如，二维空间的旋转矩阵定义如下，即将某一向量围绕原点旋转一定的角度。
$A=\begin{bmatrix} \cos\theta & -\sin \theta \\ \sin \theta & \cos \theta \end{bmatrix}$
二维空间中的拉伸矩阵定义如下，即将某一向量沿x轴拉伸一个常数因子k。
$B=\begin{bmatrix} k &0 \\ 0&1 \end{bmatrix}$
类似的，在y轴有一个拉伸，矩阵定义如下。
$C=\begin{bmatrix} 1 &0 \\ 0&k \end{bmatrix}$
如果有一个向量：
$x=\begin{bmatrix} 1\\ 0 \end{bmatrix}$
那么y=Ax是在向量x上旋转一定角度得到的向量，而Bx是在x方向上进行拉伸得到的向量。如下图所示：
在这里插入图片描述
如果是其他类型的变换矩阵呢？假设
$A=\begin{bmatrix} 3 & 2 \\ 0 &2 \end{bmatrix}$
如果我们将该矩阵作用于一堆向量上而不是仅一个向量，假设我们有一个圆，其中包含距离原点一个单位的所有向量，则该向量一般形式为：
$x=\begin{bmatrix} x_{i} \\ y_{i} \end{bmatrix}$
将矩阵A作用于向量x上，则最初的向量x形成的是一个圆，矩阵A作用后将其变成了一个椭圆，如下图所示。
请添加图片描述
假设圆中的样本向量 x1 和 x2 分别变换为 t1 和 t2，则有：
$t_{1} = A x_{1} \\ t_{2} = A x_{2}$
向量是一个既有大小又有方向的量。上式例子中，矩阵A对x向量的影响是旋转和拉伸的组合，其中对x1的影响是既改变了大小又改变了方向，对x2只有幅度上的变化，因此改变矢量大小而不改变其方向的唯一方法是将其乘以标量。对于上式x2这样的向量，乘以A的效果等同于与标量相乘，如下所示。
$t_{2} = A x_{2} = \lambda x_{2}$
以上式子不适用于x中的所有向量。实际上，对于每个矩阵A，只有部分向量具有这个性质，这个特殊向量称为A的特征向量，对应的标量lambda称为该特征向量的特征值。因此，n x n的矩阵A的特征向量被定义为非零向量u，使得：
$\lambda u$
任何与特征向量u具有相同方向的向量也具有对应特征值的特征向量。

举个例子，有一个矩阵:
$B=\begin{bmatrix} -1 &1 \\ 0&-2 \end{bmatrix}$
其特征值分别为-1和-2，对应的特征向量为：
$u_{1} = \begin{bmatrix} 1 \\ 0 \end{bmatrix}$
$u_{1} = \begin{bmatrix} -1 \\ 1 \end{bmatrix}$
这意味着矩阵B作用于所有可能向量时，它不会改变这两个向量的方向，而只会拉伸它们，所以对于特征向量，矩阵乘法变成了一个简单的标量乘法。

已知矩阵C并计算其特征向量：
$C=\begin{bmatrix} 3 & 2 \\ 0 &2 \end{bmatrix}$
得到两个特征向量分别为：
$u_{1} = \begin{bmatrix} 1 \\ 0 \end{bmatrix}$
$u_{2} = \begin{bmatrix} -0.8944\\ 0.4472 \end{bmatrix}$
对应的特征值为：
$\lambda _{1} =3\\ \lambda _{2} =2$
请添加图片描述
对于矩阵D并计算其特征向量：
$D=\begin{bmatrix} 3 & 1 \\ 1 &2 \end{bmatrix}$
其特征向量分别为：
$u_{1} = \begin{bmatrix} 0.8507\\ 0.5257 \end{bmatrix}$

$u_{2} = \begin{bmatrix} -0.5257\\ 0.8507 \end{bmatrix}$
对应的特征值为：
$\lambda _{1} =3.618\\ \lambda _{2} =1.382$
请添加图片描述
观察可发现，特征向量是沿着椭圆的长轴和短轴，那为什么矩阵C的特征向量没有这个性质呢？那是因为矩阵D是一个对称矩阵，对称矩阵是等于其转置矩阵，该矩阵的主对角线上的元素是任意的，但对于其他元素，第i行第j列上的每一个元素等于第j行第i列的元素，以下是一个对称矩阵的例子：
$\begin{bmatrix} 5&0 &4 &3 \\ 0 & 7 & 9 & 2\\ 4 & 9 & 6 & 1\\ 3 & 2 & 1&3 \end{bmatrix}$
对称矩阵始终是方阵，它是用过沿其特征向量进行拉伸或缩放来变换向量，因此则有：
$Du_{1}=\lambda u_{1}\\ Du_{2}=\lambda u_{2}$
请添加图片描述
对称矩阵的性质：
以上实例可以看出：对称矩阵D的特征向量相互垂直并形成正交向量。
对称矩阵的一个重要性质是n x n对称矩阵有n个线性独立且正交的特征向量，并且它有n个与这些特征向量对应的实特征值，这些特征值不一定彼此不相等。对称矩阵另一个重要特性是它们是正交可对角化的。

特征分解：
对称矩阵是可以正交对角化的，意味着如果我们有一个对称矩阵A，可以将其分解为：
$A=PDP^{T}$
其中，D是由A的n个特征值组成的n x n对角矩阵，P也是n x n矩阵，P的列是A的n个线性独立特征向量，分别对应于D中的那些特征值，如果u1, u2, u3 …, un是A的特征向量， λ1, λ2, …, λn分别是它们对应的特征值，那么A可以写成：
$\begin{bmatrix} u_{1} & u_{2} & \dots &u_{n} \end{bmatrix}\begin{bmatrix} \lambda _{1} & 0 & \dots& 0\\ 0 & \lambda _{2} & \dots &0 \\ \vdots & \vdots& \vdots &\vdots \\ 0 & 0 & \dots &\lambda _{n} \end{bmatrix}\begin{bmatrix} u_{1}^{T} \\ u_{2}^{T} \\ \vdots \\ u_{n}^{T} \end{bmatrix}$
以上, A的这种因式分解成为特征分解。

特征分解的几何解释：
为了更好的理解特征分解方程，首先需要对其进行化简，如果假设每一个特征向量ui是一个n x 1的列向量。
$u_{i}= \begin{bmatrix} u_{i1} \\ u_{i2} \\ \vdots \\ u_{in} \end{bmatrix}$
转置后：
$u_{i}^{T} = \begin{bmatrix} u_{i1} & u_{i2} & \dots &u_{in} \end{bmatrix}$
相乘后，得到n x n矩阵：
$u_{i}u_{i}^{T} = \begin{bmatrix} u_{i1} \\ u_{i2} \\ \vdots \\ u_{in} \end{bmatrix} \begin{bmatrix} u_{i1} & u_{i2} & \dots &u_{in} \end{bmatrix}=\\ \begin{bmatrix} u_{i1}u_{i1} & u_{i1}u_{i2} & \dots & u_{i1}u_{in}\\ u_{i2}u_{i1} & u_{i2}u_{i2} & \dots & u_{i2}u_{in}\\ \vdots & \vdots & \vdots &\vdots \\ u_{in}u_{i1} & u_{in}u_{i2}& \dots &u_{in}u_{in} \end{bmatrix}$
首先，我们计算 DP^T 以简化特征分解方程：
$\begin{bmatrix} \lambda _{1} &0 &\dots &0 \\ 0 &\lambda _{2} &\dots & 0\\ \vdots & \vdots & \vdots & \vdots\\ 0 & 0 &\dots &\lambda _{n} \end{bmatrix}\begin{bmatrix} u_{1}^{T} \\ u_{2}^{T} \\ \vdots \\ u_{n}^{T} \end{bmatrix}=\begin{bmatrix} \lambda _{1}u_{1}^{T} \\ \lambda _{2}u_{2}^{T} \\ \vdots \\ \lambda _{n}u_{n}^{T} \end{bmatrix}$
现在特征分解方程变为：
$A=\begin{bmatrix} u_{1} & u_{2} & \dots &u_{n} \end{bmatrix}\begin{bmatrix} \lambda _{1}u_{1}^{T} \\ \lambda _{2}u_{2}^{T} \\ \vdots \\ \lambda _{n}u_{n}^{T} \end{bmatrix}\\ =\lambda _{1}u_{1}u_{1}^{T} + \lambda _{1}u_{1}u_{1}^{T} + \dots +\lambda _{n}u_{n}u_{n}^{T}$
因此nxn矩阵A可以分解成n个形状相同（n x n）的矩阵，每个矩阵都有对应的特征值λi。每个矩阵称为投影矩阵：
$u_{i}u_{i}^{T}$
假设我们有一个向量 x 和一个单位向量 v，v 和 x 的内积等于 vx=v^T x 给出了 x 在 v 上的标量投影（这是 x 到 v 的向量投影的长度 v)，如果我们再次将它乘以 v，它会给出一个向量，称为 x 到 v 上的正交投影。如下图所示：
请添加图片描述
因此将 ui ui^T 乘以 x，我们得到 x 到 ui 的正交投影。
现在我们计算一下之前提到的矩阵D的投影矩阵：
$D=\begin{bmatrix} 3 & 1 \\ 1 &2 \end{bmatrix}$
得到特征分解方程中的第一项：
$D_{1}= \lambda _{1}u_{1}u_{1}^{T} =3.618\begin{bmatrix} 0.8507 \\ 0.5257 \end{bmatrix}\begin{bmatrix} 0.8507 &0.5257 \end{bmatrix}\\ =\begin{bmatrix} 2.618 &1.618 \\ 1.618 &1 \end{bmatrix}$
可见，它是一个对称矩阵。事实上，特征分解方程中的所有投影矩阵都是对称的，是因为每一项的第m行第n列的元素与第n行第m列的元素具有相同的值。
因此，特征分解方程是具有 n 个特征向量的对称 n×n 矩阵。特征向量与原始矩阵 D相同，即 u1, u2, … un， ui 对应的特征值为 λi（与 D 相同），但其他所有特征值均为零。对称矩阵转换向量是将沿其特征向量拉伸或收缩向量，拉伸或收缩的量与对应的特征值成正比，所以这个矩阵将沿着 ui 拉伸一个向量，但由于其他特征值为零，它会在这些方向上将其缩小为零。

假设我们将对称矩D作用于任意向量x，现在特征分解方程变为：
$Dx=\lambda _{1}u_{1}u_{1}^{T}x + \lambda _{1}u_{1}u_{1}^{T} x+ \dots +\lambda _{n}u_{n}u_{n}^{T}x$
因此在特征分解方程的每一项中，它是x到ui到正交投影，然后再乘以λi，由于λi是标量，它乘以一个向量只会改变向量到大小，而不改变向量到方向。如下图所示。

请添加图片描述
所以，特征分解在数学上解释了我们之前在图中看到的对称矩阵的一个重要性质，对称矩阵通过沿其特征向量拉伸或者收缩向量来变换向量，并且沿每个特征向量的拉伸或收缩量与对应的特征值呈正比。此外，特征分解可以将一个n x n对称矩阵分解为n个具有相同形状（n x n）的矩阵乘以一个特征值。投影矩阵仅将 x 投影到每个 ui 上，但特征值缩放矢量投影的长度 (ui ui^ T x)。特征值越大，得到的向量 (λi ui ui^ T x) 的长度越大，赋予其对应矩阵 (ui ui^ T) 的权重越大。因此，我们可以通过对具有最高特征值的项求和来近似我们的原始对称矩阵 D。例如，如果我们假设特征值 λi 已按降序排序：
$D=\lambda _{1}u_{1}u_{1}^{T} + \lambda _{1}u_{1}u_{1}^{T} + \dots +\lambda _{n}u_{n}u_{n}^{T}\\ \lambda _{1}\ge \lambda _{2}\ge \dots \ge \lambda _{k}\ge \dots \ge \lambda _{n}$
那么我们取特征分解方程中的前k项来对原始矩阵有一个很好的近似：
$\approx D_{k} =\lambda _{1}u_{1}u_{1}^{T} + \lambda _{1}u_{1}u_{1}^{T} + \dots +\lambda _{k}u_{k}u_{k}^{T}$

因此在原始矩阵 D 中，我们遗漏的其他 (n-k) 个非常小且接近于零的特征值，则近似矩阵与原始矩阵非常相似，并且我们有一个很好的近似。如下图所示，将矩阵C作用于向量x（左图）上得到Cx（中图），而特征分解方程的第一项（右图）就已经很接近Cx（中图）了。
请添加图片描述

那么为什么特征分解方程是有效的？为什么它需要一个对称矩阵呢？
源于对称矩阵的重要性质。假设 x 是一个 n×1 列向量。如果 A 是一个 n×n 对称矩阵，那么它有 n 个线性独立且正交的特征向量，可以作为一个新的基，并且基于基的定义，任何向量 x 都可以唯一地写为 D 的特征向量的线性组合。如果是非对称矩阵，则这里的特征向量是线性独立的，而非正交，且它们没有显示变换后该矩阵的正确拉伸方向。

**因此，特征分解方法是非常有用的，但仅适用于对称矩阵，对称矩阵一定是一个方阵，如果有一个不是方阵的矩阵或者一个方阵但是是非对称矩阵，那么就不能使用特征分解方法来用其他矩阵逼近它，SVD分解能克服这个问题。

2 SVD分解

2.1 奇异值

在讨论SVD之前，我们应该找到一种方法来计算非对称矩阵的拉伸方向。
假设A是一个 m x n 矩阵，它不是一个对称矩阵。那么可以证明矩阵
$A^{T} A$
是一个n x n对称矩阵。由线性代数相关知识可知：
$\left (A^{T}A \right ) ^{T} =A^{T}\left ( A^{T} \right ) ^{T} =A^{T}A$
所以A^T A 等于它的转置，是一个对称矩阵。我们想计算非对称矩阵的拉伸方向。但是我们如何在数学上定义拉伸方向？

假设现在列向量有 3 个元素。最初，我们有一个球体，其中包含距离原点一个单位的所有向量，如下图所示。如果我们将这些向量称为 x，则 ||x||=1。现在，如果我们将它们乘以一个 3×3 对称矩阵，Ax 就变成了一个 3-d 椭圆。第一拉伸方向可以定义为该椭圆中长度最大的向量的方向（下图中的 Av1）。实际上Av1是||Ax||的最大值在所有单位向量 x 上。这个向量是向量 v1 通过 A 的变换。
请添加图片描述
拉伸的第二个方向是沿矢量 Av2， Av2 是 ||Ax|| 的最大值在 x 中垂直于 v1 的所有向量上，所以在 x 的所有向量中，我们最大化 ||Ax||在这个约束下，x 垂直于 v1，最后v3 是垂直于 v1 和 v2 的向量，并且在这些约束下给出了 Ax 的最大长度，Av3 的方向决定了拉伸的第三方向。所以一般在n维空间中，第i个拉伸方向是向量Avi的最大长度方向，并且垂直于之前的(i-1)个拉伸方向。

因此，对于一个m x n的矩阵A，A^ T A是一个对称矩阵，它有 n 个实特征值和 n 个线性独立和正交特征向量，这些特征向量可以构成它可以变换的 n 个元素向量的基础（在 R^n 空间中）。

2.2 Singular Value Decomposition (SVD)

令 A 为 m×n 矩阵，秩 A = r，所以 A 的非零奇异值的个数是 r。由于它们是正数并按降序标记，我们可以将它们写为
$\sigma _{1} \ge \sigma _{2} \ge \dots \ge \sigma _{n}$
其中
$\sigma _{r+1} = \sigma _{r+2} = \dots = \sigma _{n} =0$

我们知道每个奇异值 σi 是 λi 的平方根（A^T A的特征值），对应一个同阶的特征向量 vi。现在我们可以将 A 的奇异值分解写为：
$A=U\Sigma V^{T}$
其中 V 是一个 n×n 矩阵，其列是 vi。则有：
$V=\begin{bmatrix} v_{1} & v_{2} &\dots &v_{n} \end{bmatrix}$
我们称一组正交且归一化向量为正交集合，所以集合 {vi} 是一个正交集合，V是正交矩阵。

Σ 是一个 m×n 对角矩阵，其形式为：
请添加图片描述
所以我们首先制作一个 r × r 对角矩阵，其中对角元素为 σ1, σ2, …, σr，然后我们用零填充它，使其成为一个 m × n 矩阵。
我们还知道集合 {Av1, Av2, …, Avr} 是 Col A 的正交基，并且 σi = ||Avi||。所以我们可以通过将 Avi 向量除以它们的长度来归一化它们：
$u_{i}=\frac{Av_{i}}{\parallel A v_{i}\parallel }=\frac{A v_{i}}{\sigma _{i} } （1< i< r）$
现在我们有一个集合 {u1, u2, …, ur}，它是 r 维的 Ax 的标准正交基。我们知道A是一个m×n的矩阵，A的秩最多可以是m（当A的所有列都是线性独立的时候）。由于我们需要 U 的 m×m 矩阵，我们将 (mr) 个向量添加到 ui 集合中，使其成为 m 维空间 R^m 的归一化基。
所以现在我们有一个正交基{u1, u2, … ,um}。这些向量将是 U 的列，它是一个正交 m×m 矩阵。
$U=\begin{bmatrix} u_{1} & u_{2}& \dots &u_{m} \end{bmatrix}$
综上，矩阵A最终分解为
$A=\begin{bmatrix} u_{1} & u_{2}& \dots &u_{m} \end{bmatrix}\begin{bmatrix} \sigma _{1} &0 &\dots &0 &0 &\dots &0 \\ 0 &\sigma _{2} &\dots &0 &0 &\dots &0 \\ \vdots &\vdots &\ddots &\vdots &\vdots &\ddots &\vdots \\ 0 &0 &\dots &\sigma _{r} &0 &\dots &0 \\ 0 &0 &\dots &0 &0 &\dots &0 \\ \vdots &\vdots &\ddots &\vdots &\vdots &\ddots &\vdots \\ 0 &0 &\dots &0 &0 &\dots &0 \end{bmatrix}\begin{bmatrix} v_{1}^{T} \\ v_{2}^{T} \\ \vdots \\ v_{n}^{T} \end{bmatrix}$

简化为：
$A=\begin{bmatrix} u_{1} & u_{2}& \dots &u_{m} \end{bmatrix}\begin{bmatrix} \sigma _{1}v_{1}^{T} \\ \sigma _{2}v_{2}^{T} \\ \vdots \\ \sigma _{r}v_{r}^{T} \\ 0 \\ \vdots \\ 0 \end{bmatrix}\\ =\sigma _{1}u_{1}v_{1}^{T} +\sigma _{2}u_{2}v_{2}^{T}+\dots +\sigma _{r}u_{r}v_{r}^{T}$
所以SVD方程将矩阵A分解成r个形状相同（m x n）的矩阵。
所以上式中的每一项 ai = σivi ^ Tx 是 Ax 在 ui 上的标量投影，如果乘以 ui，结果是一个向量，它是 Ax 在 ui 上的正交投影，奇异值 σi 沿 ui 缩放该向量的长度。所以在特征分解方程中，每个 ui ui^T 是一个投影矩阵，它将给出 x 到 ui 的正交投影。

2.3 summary

从以上可以看出，任何矩阵都可以分解为向量u和向量v，现在需要思考的是矩阵是如何分解的，SVD的思想如下：
$A=U\Sigma V^{T}$

U是一个为mxm的矩阵；
V是一个为nxn的矩阵，和矩阵U一样也为酉矩阵；
sigma是一个为mxn的矩阵，除了主对角线上的元素以外全为0，主对角线上的每一个元素都称为奇异值。

在这里插入图片描述

通过以下公式求出SVD分解后的三个矩阵：
首先，将A的转置和A做矩阵乘法运算，得到一个nxn的方阵A^{T}A，此时可以进行特征分解，得到特征值和特征向量满足下式：
$\left ( A^{T}A \right )v_{i}=\lambda _{i}v_{i}$
因此得到了该nxn方阵的n个特征值和对应的n个特征向量v，如果将所有的特征向量张成一个nxn的矩阵V，就是SVD公式里面的V矩阵，一般称V中的每个特征向量叫做A的右奇异向量。

其次，将A和A的转置做矩阵乘法，得到一个mxm的方阵AA^{T}，同样，对它进行特征分解，得到特征只和特征向量满足下式：
$\left ( AA^{T} \right )v_{i}=\lambda _{i}v_{i}$
因此得到了该mxm方阵的m个特征值和对应的m个特征向量u，同样，将所有的特征向量张成一个mxm的矩阵U，就是SVD公式里的U矩阵，一般称U中的每个特征向量叫做A的左奇异值。

对于奇异值矩阵，因为该矩阵除了对角线上是奇异值其他位置都是0，因此只需要求出奇异值即可。
通过以下公式可以求出每一个奇异值，进而求出奇异值矩阵。
$A=U\Sigma V^{T}\Rightarrow AV=U\Sigma V^{T}V\Rightarrow AV=U\Sigma \Rightarrow Av_{i}=\sigma _{i}u_{i}\Rightarrow \sigma _{i}=Av_{i}/u_{i}$

用图像进行解释说明，矩阵A可以认为将图像进行矢量化并存储为列，矩阵U可以认为与A中的列对应的相似列。与特征向量类似，我们将这些“向量面”成为特征面，用特征面来表示人脸，这是图像处理和计算机视觉中用于人脸分析的基本方法。
在这里插入图片描述
矩阵U和V为实对称矩阵，即将U和UT相乘、V和VT相乘，得到的是单位矩阵。

因此，

3 SVD代码实现

3.1 基于opencv实现SVD

#include <iostream>

#include "eigen3/Eigen/Core"
#include "eigen3/Eigen/SVD"

#include <opencv2/opencv.hpp>
#include <opencv2/core/eigen.hpp>

cv::Mat load(const std::string & p_filePath)
{
  cv::Mat mat;
  cv::FileStorage fs;
  if (!fs.open(p_filePath, cv::FileStorage::READ))
    {
      std::cout << "Can't open matrix file" << std::endl;
      return cv::Mat();
    }
  fs["matrix"] >> mat;
  fs.release();
  return mat;
}

int main (void)
{
  cv::Mat src = load("svd_src.xml"); // input (rows: 6 ; cols: 515 808 ; type: CV_64FC1)
  cv::Mat u; // output
  

  // Eigen implementation
  Eigen::MatrixXd m;
  cv::cv2eigen(src, m);

  std::cout << "[start] Eigen SVD decomposition" << std::endl;
  Eigen::JacobiSVD<Eigen::MatrixXd> svd(m, Eigen::ComputeThinU);
  std::cout << "[end] Eigen SVD decomposition" << std::endl;
    
  cv::eigen2cv(svd.matrixU(), u);


  // OpenCV implementation
  cv::Mat w, vt;
  std::cout << "[start] OpenCV SVD decomposition" << std::endl;
  cv::SVD::compute(src, w, u, vt);
  std::cout << "[end] OpenCV SVD decomposition" << std::endl;

  return 0;
}

3.2 基于eigen实现SVD

#include <Eigen/Dense>
#include <iostream>
int main()
{
Eigen::MatrixXf A = Eigen::MatrixXf::Random(3, 2);
std::cout << "A: " << A << std::endl;
Eigen::JacobiSVD<Eigen::MatrixXf> svd(m, Eigen::ComputeFullV | Eigen::ComputeFullU);
std::cout << "singular_values: " << svd.singularValues() << std::endl;
std::cout << "U: " << svd.matrixU() << std::endl;
std::cout << "V: " << svd.matrixV() << std::endl;
return 0;
}

4 应用

输入一张人脸图像

faces = X[0,:].reshape((64,64))

通过奇异值进行分解

# we decompose the image of a face
U, S, VT = np.linalg.svd(faces)
S = np.diag(S)
faces = faces.astype('float')
#U, VT = truncate_u_v(S, U, VT)
draw_svd(faces, U, S, VT, 'gray')

输出为：
在这里插入图片描述
可以求出前30个特征值

diag_elem = np.diag(S)
plt.stem(diag_elem[1:])

在这里插入图片描述
当然，这里可以取前n个特征值进行重建，剩下的特征值舍弃，显然，取得的特征值越多，图像的轮廓越清晰。因为越靠前的特征值和特征向量，权重越大，越靠后，权重越小。

ps：如果矩阵A为mxn，用奇异值分解，如果A是方阵，用特征值分解即可。

5 最后，简单说下SVD和PCA之间的关系

PCA（Principal Component Analysis）和SVD（Singular Value Decomposition）是常用的用于数据降维和特征提取的线性代数算法。

他们的相似点是，都利用了矩阵分解的思想，旨在将一个高维矩阵用一组低维度的线性关系来描述。

这里简单的说明一下矩阵分解的方法，矩阵分解有很多种方法：

对于方阵，可以进行特征分解；
对于非方阵，不能进行特征分解，但可以做SVD分解。

那么PCA和SVD的区别在于分解的方式不同：

PCA求解的关键在于求解协方差矩阵 $C = (1/m) X X^T$ 的特征值分解，然后对其排序，选取前k大的特征值对应的特征向量作为基矢，从而得到投影矩阵。
SVD分解是直接对数据矩阵进行分解，即分解为三个矩阵的乘积， $A=U\Sigma V^T$ ，其中U和V是正交矩阵， $\Sigma$ 是除对角线外其他元素都为0的矩阵。其中 $\Sigma$ 的对角线元素就是奇异值，是一个降序排列的数列。将其前k个奇异值及其对应的左右奇异向量作为基矢，就可以得到投影矩阵。