四、矩阵的奇异值分解

Singular Value decomposition是一种矩阵分解方法，它是科学工程计算和数值代数的最有用和最有效的工具之一。

奇异值分解的由来

1873年Beltrami从双线性函数\(f(x,y) = x^{T}Ay,A \in R^{n \times n}\)出发，引入线性变换\(x = U\xi,y = V\eta\)，这样双线性函数变为\(f(x,y) = \xi^{T}S\eta,S = U^{T}AV\)，Beltrami观测到，如果选择约束\(U\)和\(V\)为正交矩阵，则有可能使矩阵\(S\)的对角线以外的元素全部为零，即矩阵

\[S = \Sigma = diag(\sigma_{1},\sigma_{2},\ldots,\sigma_{n})\]

为对角阵，然后可得\(A = U\Sigma V^{T}\).

从几何的观点来看\(SVD\)，这里用到一个基本事实： 任意一个\(m \times n\)矩阵将\(n\)维空间的单位球面映射为\(m\)维空间的一个超球面。

从高维空间看，设\(A \in R^{n \times n},rank(A) = n\),\(S\)是\(R^{n}\)中的单位球面，\(S\)经\(A\)映射后，其像\(AS\)在\(R^{n}\)中为一椭球面，它的\(n\)个主半轴的长度记成\(\sigma_{1},\sigma_{2},\ldots,\sigma_{n}\)，习惯上假设奇异值以降序编号，\(\sigma_{1} \geq \sigma_{2} \geq \ldots \geq \sigma_{n} > 0\).

其次我们定义\(AS\)的\(n\)个主半轴方向的单位向量\(u_{1},u_{2},\ldots,u_{n}\)为左奇异向量，其编号对应与奇异值，向量\(\sigma_{i}u_{i}\)就是\(AS\)的第\(i\)大的主半轴。

最后我们定义\(S\)中的单位向量\(v_{1},v_{2},\ldots,v_{n} \in S\)，它是\(AS\)的\(n\)个主半轴的原像，称为\(A\)的右奇异向量，编号使得\(Av_{j} = \sigma_{j}u_{j}\).

约化SVD

\(AV = \overline{U}\overline{\Sigma}\)

其中\(\overline{\Sigma}\)是\(n \times n\)矩阵，它有正的实元素（假设\(A\)是列满秩\(n\)），\(\overline{U}\)是一个有正交列的\(m \times n\)矩阵。\(V\)是一个具有正交列的\(n \times n\)矩阵，因此\(V\)是酉矩阵，可以右乘其逆\(V^{H}\)，这样得到

\[A = \overline{U}\overline{\Sigma}V^{H}\]

这个\(A\)的矩阵分解称为\(A\)的约化奇异值分解（reduced singular value decomposition）或约化SVD。

完全SVD

定理4.1（实矩阵的奇异值分解）令\(A \in R^{m \times n}\)，则存在正交矩阵\(U \in R^{m \times m}\)和\(V \in R^{n \times n}\)使得

\[A = U\Sigma V^{T}\]

其中\(\Sigma = \begin{bmatrix} \Sigma_{1} & 0 \\ 0 & 0 \\ \end{bmatrix}_{m \times n}\), 且\(\Sigma_{1} = diag(\sigma_{1},\sigma_{2},\ldots,\sigma_{r})\)，对角线元素按照从大到小的次序排列，即

\[\sigma_{1} \geq \sigma_{2} \geq \ldots \geq \sigma_{r} > 0,r = rank(A)\]

数值\(\sigma_{1},\sigma_{2},\ldots,\sigma_{r}\)连同\(\sigma_{r + 1} = \sigma_{r + 2} = \ldots = \sigma_{\min(m,n)} = 0\)称为\(A\)的奇异值。

定理设\(A \in R^{m \times n}\)，则总存在\(x \in R^{m},y \in R^{n}\)，且\(\parallel x \parallel_{2} = 1, \parallel y \parallel_{2} = 1\)，使得\(Ay = \sigma x,\sigma = \parallel A \parallel_{2}\).

证明矩阵的SVD分解定理：

令\(x \in R^{n},y \in R^{n}\) 分别是\(n\)维、\(m\)维单位向量，满足\(Ax = \sigma y\)，其中\(\sigma = \parallel A \parallel_{2}\).

则总存在\(V_{2} = R^{n \times (n - 1)}\)和\(U_{2} \in R^{m \times (m - 1)}\)使得\(V = \lbrack x:V_{2}\rbrack \in R^{n \times n},U = \lbrack y:U_{2}\rbrack \in R^{m \times m}\)为正交矩阵，这样就有\(U^{T}AV = \begin{bmatrix} \sigma & w^{T} \\ 0 & B \\ \end{bmatrix} \equiv A_{1},w \in R^{n - 1},B \in R^{(m - 1) \times (n - 1)}\)，由于\(\parallel A_{1}\begin{bmatrix} \sigma \\ w \\ \end{bmatrix} \parallel_{2}^{2} = \parallel \begin{bmatrix} \sigma^{2} + w^{T}w \\ Bw \\ \end{bmatrix} \parallel_{2}^{2} \geq (\sigma^{2} + w^{T}w)^{2}\)，我们有\(\parallel A_{1} \parallel_{2}^{2} \geq (\sigma^{2} + w^{T}w)^{2}\)（这里用到\(\parallel A_{1} \parallel_{2}^{2} \geq \frac{\parallel A_{1}\begin{bmatrix} \sigma \\ w \\ \end{bmatrix} \parallel_{2}^{2}}{\parallel \begin{bmatrix} \sigma \\ w \\ \end{bmatrix} \parallel_{2}^{2}}\)）。

而\(\parallel A \parallel_{2}^{2} = \parallel A_{1} \parallel_{2}^{2}\)(这里用到了矩阵2-范数为正交不变范数)。

所以\(w^{T}w = 0,w = 0\)，然后对\(B\)做同样的证明。

定理4.2(复数矩阵的奇异值分解) 令\(A \in R^{m \times n}\)，则存在酉矩阵\(U \in R^{m \times m}\)和\(V \in R^{n \times n}\)使得

\[A = U\Sigma V^{H}\]

\[\sigma_{1} \geq \sigma_{2} \geq \ldots \geq \sigma_{r} > 0,r = rank(A)\]

数值\(\sigma_{1},\sigma_{2},\ldots,\sigma_{r}\)连同\(\sigma_{r + 1} = \sigma_{r + 2} = \ldots = \sigma_{\min(m,n)} = 0\)称为\(A\)的奇异值。

SVD的若干性质

\(n \times n\)矩阵\(V\)为酉矩阵，用\(V\)右乘\(A = U\Sigma V^{H}\)就得到\(AV = U\Sigma\)，这个向量的形式就

\[Av_{i} = \left\{ \begin{matrix} \sigma_{i}u_{i} & i = 1,2,\ldots,r \\ 0 & i = r + 1,r + 2,\ldots,\min(m,n) \\ \end{matrix} \right.\ \]

因此，\(V\)的列向量\(v_{i}\)称为矩阵\(A\)的右奇异向量(right singular veector)，\(V\)称为\(A\)的右奇异向量矩阵。
\(m \times m\)矩阵\(U\)为酉矩阵，用\(U^{H}\)左乘式\(A = U\Sigma V^{H}\)就得到\(U^{H}A = \Sigma V\)，这个列向量的形式就是

\[u_{i}^{H}A = \left\{ \begin{matrix} \sigma_{i}v_{i}^{H} & i = 1,2,\ldots,r \\ 0 & i = r + 1,r + 2,\ldots,\min(m,n) \\ \end{matrix} \right.\ \]

因此，\(U\)的列向量\(u_{i}\)称为矩阵\(A\)的左奇异向量(left singular vector)，矩阵\(U\)称为\(A\)的左奇异向量矩阵。
矩阵\(A\)的奇异值分解式\(A = U\Sigma V^{H}\)可以改写成向量表达式

\[A = \Sigma_{i = 1}^{r}\sigma_{i}u_{i}v_{i}^{H}\]

就是\(r\)个秩1矩阵之和。
当矩阵\(A\)的秩\(r = rank(A) < \min\{ m,n\}\)时，由于奇异值\(0 = \sigma_{r + 1} = \ldots = \sigma_{\min\{ m,n\}}\)，故奇异值分解可以简化为

\[A = U_{r}\Sigma_{r}V_{r}^{H}\]

其中\(U_{r} = \lbrack u_{1},u_{2},\ldots,u_{r}\rbrack,V_{r} = \lbrack v_{1},v_{2},\ldots,v_{r}\rbrack,\Sigma_{r} = diag(\sigma_{1},\sigma_{2},\ldots,\sigma_{r})\)则上述简化分解式称为\(A\)的截断奇异值分解(truncated SVD)或者紧凑奇异值分解。
用\(u_{i}^{H}\)左乘（1）中的式子，由于\(u_{i}^{h}u_{i} = 1\)，因此得到

\[u_{i}^{H}Av_{i} = \sigma_{i},\quad i = 1,2,\ldots,\min\{ m,n\}\]

容易得到\(AA^{H} = U\Sigma^{2}U^{H}\), 这说明一个重要的事实，就是\(m \times n\)的奇异值\(\sigma_{i}\)就是矩阵乘积\(AA^{H}\)的特征值的开方的算数平方根。
如果矩阵\(A_{m \times n}\)的秩为\(r\)，那么对应于\(A\)的奇异值分解

\[A = U\Sigma V^{H}\]

，就有下列结果\(m \times m\)酉矩阵\(U\)的前\(r\)列组成矩阵\(A\)的列空间\(R(A) = \{ y \in C^{m}|Ax = y,\forall x \in C^{n}\}\)的标准正交基底；

\(U\)的后\(m - r\)列向量组成了\(A^{H}\)的零空间\(N(A^{H}) = \{ x \in C^{m}|A^{H}x = 0\}\)的标准正交基底；

\(n \times n\)酉矩阵\(V\)的前\(r\)列组成的矩阵\(A\)的行空间\(R(A^{H})\)的标准正交基底；

\(V\)的后\(n - r\)列组成矩阵\(A\)的零空间\(N(A)\)的标准正交基底。

酉矩阵\(U\)和\(V\)的列向量分别给出了相关矩阵\(R(A),R(A^{H}),N(A),N(A^{H})\)的标准正交基底。

其实当\(A\)是一个\(m \times n\)的矩阵时，\(A\)的非零奇异值就是\(A^{H}A\)或者\(AA^{H}\)的非零特征值的开方，所不同的是零特征值的重数不同。而\(U\)的\(m\)个列向量是\(AA^{H}\)的特征向量；\(V\)的\(n\)个列向量是\(A^{H}A\)的特征向量。

这样如果\(A\)是一个\(n \times n\)的正方矩阵，只要\(A\)有一个奇异值为零，那么\(A\)一定是奇异矩阵；事实上如果\(A\)的奇异值按照从大到小的次序排列，那么\(A\)的最小奇异值就反映了矩阵\(A\)接近奇异矩阵的距离，如果其最小奇异值为零，那么\(A\)就是一个奇异矩阵；如果最小奇异值不为零但很小，就说明\(A\)尽管不是奇异矩阵，但它接近于奇异矩阵。

SVD与条件数等数值的关系

定义设非奇异矩阵\(A \in R^{m \times m}\)，则定义\(\parallel A \parallel \parallel A^{- 1} \parallel\)为\(A\)的条件数，记为\(cond(A) = \parallel A \parallel \parallel A^{- 1} \parallel\)。

显然我们这里仅定义了非奇异矩阵的条件数，以后引进了广义逆矩阵的概念可以对一般的\(m \times n\)矩阵定义条件数。

这里需要指出两点：

矩阵的条件数是有矩阵\(A\)和\(A^{- 1}\)的范数来确定的，所用范数不同条件数也不同，如果使用矩阵的\(2 -\)范数，则\(A\)的条件数是\(cond_{2}(A) = \frac{\sigma_{1}}{\sigma_{n}}\)，其中\(\sigma_{1}\)是\(A\)的最大奇异值，\(\sigma_{n}\)是\(A\)的最小奇异值。
矩阵的条件数是矩阵\(A\)的固有性质，它对解的精度会产生重大的影响。

奇异值与矩阵的范数、行列式、特征值都有密切的关系：

矩阵\(A\)的\(2 -\)范数（谱范数）等于\(A\)的最大奇异值，即

\[\parallel A \parallel_{2} = \sigma_{1} = \sigma_{\max}\]

这是因为\(A\)的\(2 -\)范数\(\parallel A \parallel_{2} = \sqrt{\lambda_{\max}(A^{H}A)} = \sigma_{1} = \sigma_{\max}\).如果\(A\)是实对称矩阵，则\(A^{T}A = A^{2}\)，所以

\[\parallel A \parallel_{2} = \sqrt{\lambda_{\max}(A^{2})} = |\lambda_{\max}(A)|\]

由于\(A\)的\(F -\)范数\(\parallel A \parallel_{F} = \sqrt{\Sigma_{i = 1}^{m}\Sigma_{j = 1}^{n}|a_{ij}|^{2}}\)是酉不变范数，即\(\parallel U^{H}AV \parallel_{F} = \parallel A \parallel_{F}\)，因此就有

\[\parallel A \parallel_{F} = \sqrt{\Sigma_{i = 1}^{m}\Sigma_{j = 1}^{n}|a_{ij}|^{2}} = \parallel U^{H}AV \parallel_{F} = \parallel \Sigma \parallel_{F} = \sqrt{\sigma_{1}^{2} + \sigma_{2}^{2} + \ldots + \sigma_{r}^{2}}\]

，从这里可以看出\(\parallel A \parallel_{2} \leq \parallel A \parallel_{F}\).
奇异值与行列式的关系：如果\(A\)是一个\(n \times n\)的正方矩阵，由于酉矩阵的行列式的绝对值等于1，所以

\[|\det(A)| = |\det(U\Sigma V^{H})| = |\det(U)| \times |\det(V^{H})| \times |\det(\Sigma)| = |\det(\Sigma)| = \sigma_{1}\sigma_{2}\ldots\sigma_{n}\]

奇异值与特征值的关系：讲到特征值只能是方阵，设\(A\)是\(n \times n\)的实对称矩阵，且其特征值为\(\lambda_{1},\lambda_{2},\ldots,\lambda_{n}(|\lambda_{1}| \geq |\lambda_{2}| \geq \ldots \geq |\lambda_{n}|)\)，奇异值为\(\sigma_{1},\sigma_{2},\ldots,\sigma_{n}(|\sigma_{1}| \geq |\sigma_{2}| \geq \ldots \geq |\sigma_{n}|\) ，则\(\sigma_{i} = |\lambda_{i}|,i = 1,2,\ldots,n\)). 如果\(A\)为\(m \times n\)矩阵，如果\(p = \min\{ m,n\}\)，且\(\sigma_{1},\sigma_{2},\ldots,\sigma_{p}\)是\(A\)的奇异值，则

\[tr(A^{H}A) = \Sigma_{i = 1}^{n}a_{i}^{H}a_{i} = \Sigma_{i}^{n}\lambda_{i}(A^{H}A) = \Sigma_{i}^{p}\sigma_{i}^{2}\]

\(A\)的谱范数\(\parallel A \parallel_{2} = \sigma_{\max}\)，矩阵\(A\)的谱条件数

\[cond_{2}(A) = \parallel A \parallel_{2} \parallel A^{- 1} \parallel_{2} = \frac{\sigma_{\max}}{\sigma_{\min}}\]

矩阵的低秩逼近

假如\(A\)是一个\(m \times n\)矩阵，不妨设它的\(rank(A) = p\)，低秩逼近问题就是要在所有的\(m \times n\)矩阵类中，确定一个rank为 \(k < p\)的\(m \times n\)矩阵\(A_{k}\)，使其与\(A\)在一定的度量意义下与\(A\)最接近。

定理4.3（The Eckhart-Young Theorem）

设\(A \in R^{m \times n}\),\(rank(A) = r\),\(A\)的奇异值分解为\(A = \Sigma_{i = 1}^{r}\sigma_{i}u_{i}v_{i}^{T}\), 令\(k < r\)，则取低秩逼近矩阵为\(A_{k} = \Sigma_{i = 1}^{k}\sigma_{i}u_{i}v_{i}^{T}\), 逼近质量为\(\min_{rank(B) = k} \parallel A - B \parallel_{2} = \parallel A - A_{k} \parallel_{2} = \sigma_{k + 1}\)

前面讨论的是\(2 -\)范数，对于\(F -\)范数，我们有

定理4.4对任意满足\(0 \leq k \leq r = rank(A)\),上述给出的\(A_{k}\)也满足

\[\parallel A - A_{k} \parallel_{F} = \min_{B \in R^{m \times n},rank(b) = k} \parallel A - B \parallel_{F} = \sqrt{\sigma_{k + 1}^{2} + \ldots + \sigma_{r}^{2}}\]

SVD的应用-图像压缩（1）

图像压缩之所以可行，是因为图像数据是高度相关的，表现在：

大多数图像内相邻像素之间有较大的相关性，存在很大的冗余度，这是空间冗余度。
序列图像前后帧之间有较大的相关性，这是时间冗余度。
用相同的码长表示不同出现概率的符号也会造成比特数的浪费，这是符号冗余度。
允许图像编码有一定的失真也是图像可压缩的一个重要原因。

图像压缩编码：

根据编码过程中是否存在信息的损耗将编码分成有损压缩编码和无损压缩编码。
根据编码原理，又可以将图像编码氛围熵编码、预测编码、变换编码和混合编码等。

奇异值分解的一个重要作用是可以降维。

如果\(A\)表示为\(n\)个\(m\)维向量即\(A \in R^{m \times n}\)，则可以用\(m + n + 1\)个\(r\)维的向量来表示，即\(U_{r} \in R^{m \times r}\)，\(\Sigma_{r} = diag(\sigma_{1},\sigma_{2},\ldots,\sigma_{r})\)，\(V_{r} \in R^{n \times r}\)。如果\(A\)的秩\(r\)远远小于\(m\)和\(n\)，则通过奇异值分解可以大大降低\(A\)的维数。

用奇异值分解来压缩图像的基本思想是对图像矩阵进行奇异值分解，选取部分奇异值和对应的左、右奇异向量来重构图像矩阵。

有损压缩\(k < r\)，只使用\(k(2n + 1)\)个数据来代替\(n^{2}\)个图像数据，这\(k(2n + 1)\)个数据分别是矩阵\(A\)的前\(k\)个奇异值，\(n \times n\)左奇异向量矩阵\(U\)的前\(k\)列向量\(U_{k} = \lbrack u_{1},u_{2},\ldots,u_{k}\rbrack\)，和\(n \times n\)右奇异向量矩阵\(V\)的前\(k\)列向量\(V_{k} = \lbrack v_{1},v_{2},\ldots,v_{k}\rbrack\)，因此图像的压缩比应该是

\[\rho = \frac{n^{2}}{k(2n + 1)}\]

一般情况下，被选择的奇异值的个数\(k\)应该满足条件

\[k(2n + 1) < < n^{2}\]

重构原图像矩阵:

\[A_{k} = \sigma_{1}u_{1}v_{1}^{T} + \sigma_{2}u_{2}v_{2}^{T} + \ldots + \sigma_{k}u_{k}v_{k}^{T}\]

重构图像\(A_{k}\)与原图像\(A\)的误差为

\[\parallel A - A_{k} \parallel_{F}^{2} = \sigma_{k + 1}^{2} + \sigma_{k + 2}^{2} + \ldots + \sigma_{r}^{2}\]

第\(i\)个奇异值\(\sigma_{i}\)对整个图像的贡献率可以定义为

\[\mathcal{E}_{i} = \frac{\sigma_{i}^{2}}{\Sigma_{j = 1}^{r}\sigma_{j}^{2}}\]

对压缩图像的质量评价主要有：均方差(MSE)和信噪比(SNR).

数字水印

数字水印模型由以下3部分组成：

水印信号（watermark）
水印嵌入算法（insertion algorithm）
水印验证、提取和检测算法（verification detection or extraction algorithm）

这个模型有两个具体过程：编码过程（encoding process）和解码过程（decoding peocess）。

(1)嵌入过程，假设\(A\)是载体图像，\(W\)是要嵌入的水印，\(\alpha\)是水印强度参数，按照以下方式构造水印图像\(A_{w}\)：

\[\left\{ \begin{matrix} A \Rightarrow USV^{T} \\ L \Leftarrow S + \alpha W \\ L \Rightarrow U_{1}S_{1}V_{1}^{T} \\ A_{w} \Leftarrow US_{1}V^{T} \\ \end{matrix} \right.\ \]

(2)提取过程，假设\(P\)为待检测图像，嵌入过程中的\(U_{1},V_{1},\alpha,S\)是保留的参数，提取步骤如下：

\[\left\{ \begin{matrix} P \Rightarrow U_{p}S_{p}V_{p}^{T} \\ F \Leftarrow U_{1}S_{p}V_{1}^{T} \\ W_{E} \Leftarrow (F - S)/\alpha \\ \end{matrix} \right.\ \]

数字水印算法的误差分析：

若\(A = \lbrack a_{ij}\rbrack \in R^{m \times n}\)，则有\(\parallel A \parallel_{2} = \sqrt{\lambda_{\max}(A^{T}A)} = \sigma_{\max}\);
若\(U \in R^{m \times m},V \in R^{n \times n}\)，是正交矩阵，并且\(A = \lbrack a_{ij}\rbrack \in R^{m \times n}\),则有\(\parallel UAV \parallel_{2} = \parallel A \parallel_{2}\);(这时矩阵\(2 -\)范数的正交不变性)
令\(A \in R^{m \times n},\delta A\)是\(A\)的误差矩阵，并且记\(A_{w} = A + \delta A\).若\(A,A_{w}\)的奇异值从大到小分别排序为\(S_{i}(A),S_{i}(A_{w})\)，则有

\[|S_{i}(A) - S_{i}(A_{w})| < \parallel \delta A \parallel_{2},i = 1,2,\ldots,n\]

结合水印加密与解密的公式，我们可以得到结论如下：若\(A,A_{w}\)及\(W\)定义如前，则有

\[|S_{i}(A) - S_{i}(A_{w})| \leq \alpha \parallel W \parallel_{2},i = 1,2,\ldots,n\]

因为图像的误差仅来自其奇异值的不同，因此可以根据不同水印、不同系数\(\alpha\)的选取，可以对图像误差等情况进行控制，这也是用\(SVD\)进行数字水印算法所具有的优势。

基于SVD矩阵分解的数字水印算法具有误差容易估计、水印叠加位置容易确定、水印叠加能量和容量容易控制的优点，由于方法中提取水印的时候需要使用嵌入水印时的正交矩阵，而这两个正交矩阵代表的正是有少许修改的原始水印的子空间，由此对原始未嵌入水印的图像和随机选取图像误检测为含水印图像的可能性比较大。这是这种方法的主要不足，有待进一步改进。

可惜流光

4-矩阵的奇异值分解