矩阵分析与应用-6.2~6.3-奇异值分解-Section2
前言
本文学习过程来源是《矩阵分析与应用-张贤达》一书. 可以通过 z-lib 下载.
一、奇异值的性质汇总
1. 奇异值服从的等式关系
矩阵 \(A_{m \times n}\) 和其复共轭转置矩阵 \(A^{\mathrm{H}}\) 具有相同的奇异值.
矩阵 \(A_{m \times n}\) 的非零奇异值是 \(AA^{\mathrm{H}}\) 或者 \(A^{\mathrm{H}A}\) 的非零特征值的正平方根.
\(\sigma > 0\) 是矩阵 \(A_{m \times n}\) 的单奇异值, 当且仅当 \(\sigma^2\) 是 \(AA^{\mathrm{H}}\) 或 \(A^{\mathrm{H}}A\) 的单特征值.
若 \(p = \min\{m,n\}\), 且 \(\sigma_1,\sigma_2,\cdots,\sigma_p\) 是矩阵 \(A_{m \times n}\) 的奇异值, 则
\[ {\mathrm{tr}}(A^{\mathrm{H}}A) = \sum_{i=1}^{p} \sigma_i^2 \]
矩阵行列式的绝对值等于矩阵奇异值之乘积, 即 \(|{\mathrm{det}}(A)|=\sigma_1\sigma_2\cdots\sigma_n\).
矩阵 \(A\) 的谱范数等于 \(A\) 的最大奇异值, 即 \(\lVert A \rVert_{\mathrm{spec}}=\sigma_{\mathrm{max}}\)
若 \(m \ge n\), 对于矩阵 \(A_{m \times n}\), 有
\[ \begin{aligned} \sigma_{\mathrm{min}}(A) &= \min \left \{ \left ( \frac{x^{\mathrm{H}}A^{\mathrm{H}}Ax}{x^{\mathrm{H}}x} \right )^{1/2} : x \neq 0 \right \} \\ &= \min \left \{ (x^{\mathrm{H}}A^{\mathrm{H}}Ax)^{1/2}:x^{\mathrm{H}}x=1,x \in C^n \right \} \end{aligned} \]
- 若 \(m \ge n\), 对于矩阵 \(A_{m \times n}\), 有
\[ \begin{aligned} \sigma_{\mathrm{max}}(A) &= \max \left \{ \left ( \frac{x^{\mathrm{H}}A^{\mathrm{H}}Ax}{x^{\mathrm{H}}x} \right )^{1/2} : x \neq 0 \right \} \\ &= \max \left \{ (x^{\mathrm{H}}A^{\mathrm{H}}Ax)^{1/2}:x^{\mathrm{H}}x=1,x \in C^n \right \} \end{aligned} \]
- 若 \(m \times m\) 矩阵 \(A\) 非奇异, 则
\[ \frac{1}{\sigma_{\min}(A)} = \max \left \{ \left ( \frac{x^{\mathrm{H}}(A^{-1})^{\mathrm{H}}A^{-1}x}{x^{\mathrm{H}}x} \right )^{1/2}:x \neq 0, x \in C^n \right \} \]
- 若 \(A = U \begin{bmatrix} \Sigma_1&O \\O&O\end{bmatrix}V^{\mathrm{H}}\) 是 \(m \times n\) 矩阵 \(A\) 的奇异值分解, 则 \(A\) 的 \(\mathrm{Moore-Penrose}\) 逆矩阵
\[ A^{\dagger} = V \begin{bmatrix} \Sigma_1^{-1} & O \\ O & O \end{bmatrix} U^{\mathrm{H}} \]
- 若 \(\sigma_1,\sigma_2,\cdots,\sigma_p\) 是 \(m \times n\) 矩阵 \(A\) 的非奇异值 (其中, \(p=\min\{m,n\}\)), 则矩阵 \(\begin{bmatrix} O&A \\A^{\mathrm{H}}&O\end{bmatrix}\) 具有 \(2p\) 个非奇异值 \(\sigma_1,\cdots,\sigma_p,-\sigma_1,\cdots,-\sigma_p\) 和 \(|m - n|\) 个零奇异值.
2. 奇异值服从的不等式关系
若 \(A\) 和 \(B\) 是 \(m \times n\) 矩阵, 对于 \(1 \le i,j \le p, i+j \le p+1(p =\min\{m,n\})\), 有
\[ \sigma_{i+j-1}(A+B) \le \sigma_i(A) + \sigma_j(B) \]
特别地, 当 \(j=1\) 时, \(\sigma_i(A+B) \le \sigma_i(A) + \sigma_i(B), i = 1,2,\cdots,p\) 成立.
对矩阵 \(A_{m \times n}, B_{m \times n}\), 有
\[ \sigma_{\max}(A+B) \le \sigma_{\max}(A) + \sigma_{\max}(B) \]
- 若 \(A\) 和 \(B\) 是 \(m \times n\) 矩阵, 则
\[ \sum_{j=1}^{p}[\sigma_j(A+B) - \sigma_j(A)]^2 \le \lVert B \rVert_{\mathrm{F}}^2, \ p = \min\{m,n\} \]
- 若 \(A_{m \times m}=[a_1,a_2,\cdots,a_m]\) 的奇异值 \(\sigma_1(A) \ge \sigma_2(A) \ge \cdots \ge \sigma_m(A)\), 则
\[ \sum_{j=1}^{k}[\sigma_{m-k+j}(A)]^2 \le \sum_{j=1}^{k}a_j^{\mathrm{H}}a_j \le \sum_{j=1}^{k}[\sigma_j(A)]^2, \quad k=1,2,\cdots,m \]
- 若 \(p=\min\{m,n\}\), 且 \(A_{m \times n}\) 和 \(B_{m \times n}\) 的奇异值排列为 \(\sigma_1(A) \ge \sigma_2(A) \ge \cdots \ge \sigma_p(A), \sigma_1(B) \ge \sigma_2(B) \ge \cdots \ge \sigma_p(B)\) 和 \(\sigma_1(A+B) \ge \sigma_2(A+B) \ge \cdots \ge \sigma_p(A+B)\), 则
\[ \sigma_{i+j-1}(AB^{\mathrm{H}}) \le \sigma_i(A)\sigma_j(B), 1 \le i,j \le p,i+j \le p+1 \]
设 \(m \times (n-1)\) 矩阵 \(B\) 是删去 \(m \times n\) 矩阵 \(A\) 任意一列得到的矩阵, 并且它们的奇异值都按照非降顺序排列, 则
\[ \sigma_1(A) \ge \sigma_1(B) \ge \sigma_2(A) \ge \sigma_2(B) \ge \cdots \ge \sigma_h(A) \ge \sigma_h(B) \ge 0 \]
式中, \(h = \min\{m, n-1\}\)
设 \((m-1) \times n\) 矩阵 \(B\) 是删去 \(m \times n\) 矩阵 \(A\) 任意一行得到的矩阵, 并且它们的奇异值都按照非降顺序排列, 则
\[ \sigma_1(A) \ge \sigma_1(B) \ge \sigma_2(A) \ge \sigma_2(B) \ge \cdots \ge \sigma_h(A) \ge \sigma_h(B) \ge 0 \]
式中, \(h = \min\{m, n-1\}\)
矩阵 \(A_{m \times n}\) 的最大奇异值满足不等式
\[ \sigma_{\max}(A) \ge [\frac{1}{n}\mathrm{tr}(A^{\mathrm{H}}A)]^{1/2} \]
二、秩亏缺最小二乘解
1. 低秩矩阵逼近
在奇异值分析的应用中, 常常需要用一个低秩的矩阵逼近一个含噪声或扰动的矩阵. 低秩矩阵好理解, 但是这和现实世界中又有什么关系呢?
举个非常简单的例子, 当拍摄一张大海的照片, 看起来就好像一个像素经过不断复杂形成的. 想象最小复制单元是一行或一列, 那么类推到矩阵中是不是就是说明组成图像的所有向量其实是线性相关的, 然后矩阵的秩就为 1, 这就叫做低秩矩阵. 但是当图中出现轮船, 海鸟, 岛屿等事物, 图片的复杂程度就会上升, 也就不能够叫做低秩矩阵.
在信号处理中, 用低秩矩阵逼近信号矩阵, 就是为了去除噪声或扰动的信号.
下面定理给出了逼近质量的评价.
定理 1 : 令 \(A \in R^{m \times n}\) 的奇异值分解由 \(A = \sum_{i=1}^{p}\sigma_iu_iv_i^{\mathrm{T}}\) 给出, 其中, \(p=\mathrm{rank}(A)\)
若 \(k < p\), 并且 \(A_k=\sum_{i=1}^{k}\sigma_iu_iv_i^{\mathrm{T}}\), 则逼近质量可分别使用谱范数和 \(\mathrm{Frobenius}\) 范数度量.
\[ \min_{\mathrm{rank}(B)=k}\lVert A-B \rVert_{\mathrm{spec}} = \lVert A-A_k \rVert_{\mathrm{spec}} = \sigma_{k+1} \tag{1} \]
\[ \min_{\mathrm{rank}(B)=k}\lVert A-B \rVert_{\mathrm{F}} = \lVert A-A_k \rVert_{\mathrm{F}} = \sqrt{\sum_{i=k+1}^{q}\sigma_i^2} \tag{2} \]
式子中, \(q=\min\{m,n\}\).
2. 有效秩
在需要计算秩 \(r\) 的估计值 \(\hat{r}\) 的方法, 在信号处理和系统理论中, 常将该估计值称为 "有效秩"
有效秩确定有以下两种常用方法
归一化奇异值方法
计算归一化奇异值
\[ \bar{\sigma_i} = \frac{\hat{\sigma_i}}{\hat{\sigma_1}} \]
选择满足准则
\[ \bar{\sigma_i} \ge \epsilon \]
的最大整数作为有效秩的估计值 \(\hat{r}\).
范数比方法
令 \(m \times n\) 矩阵 \(A_k\) 是原 \(m \times n\) 矩阵 \(A\) 的秩 \(k\) 近似, 定义该近似矩阵与原矩阵的 \(\mathrm{Frobenius}\) 范数比为
\[ v(k) = \frac{\lVert A_k \rVert_{\mathrm{F}}}{\lVert A \rVert_{\mathrm{F}}} = \frac{\sqrt{\sigma_1^2 + \sigma_2^2 + \cdots + \sigma_k^2}}{\sigma_1^2 + \sigma_2^2 + \cdots + \sigma_h^2}, \quad h = \min\{m,n\} \]
并选择满足
\[ v(k) \ge \alpha \]
的最大整数作为有效秩估计值 \(\hat{r}\). 其中 \(\alpha\) 是接近于 1 的阈值.
3. \(\mathrm{SVD}\) 的子集选择方法
算法 1 : 子集选择算法
Step 1 : 计算 \(A\) 的 \(\mathrm{SVD}\), 并确定 \(A\) 的有效秩 \(\hat{r}\)
Step 2 : 计算置换矩阵 \(P\), 使得在 \(AP=[B_1,B_2]\) 中的矩阵 \(B_1 \in C^{m \times \hat{r}}\) 的列是 "足够线性无关的"
Step 3 : 将 \(Ax=0\) 的 \(\mathrm{LS}\) 问题变换为求 \(AP\begin{bmatrix}z\\0\end{bmatrix}=b\) 的 \(\mathrm{LS}\) 解 \(z \in C^{\hat{r}}\), 即求 \(\lVert B_1z = b \rVert_2\) 的极小化变量 \(z\).
算法 2 : 低秩 \(\mathrm{LS}\) 算法, 给定 \(A \in R^{m \times n}, b \in R^m\), 计算
\[ x = \mathrm{arg}\min_z \Bigg \lVert AP\begin{bmatrix}z\\0\end{bmatrix}-b \Bigg \rVert_2 \]
Step 1 : 计算 \(\mathrm{SVD}\), 确定有效秩 \(\hat{r}\), 根据 \(\hat{r}\) 将 \(V\) 分块为
\[ V = \begin{bmatrix} V_{11}& V_{12}\\ V_{21}& V_{22} \end{bmatrix} \]
其中, \(V_{11} \in R^{\hat{r} \times \hat{r}}\), 然后存储 \(V_{11}\) 和 \(V_{21}\)
Step 2 : 利用列主元 \(\mathrm{QR}\) 算法计算 \(Q^\mathrm{T}[V_{11}^\mathrm{T},V_{21}^\mathrm{T}]P=[R_{11},R_{12}]\), 然后计算 \(AP=[B_1,B_2]\), 其中 \(B_1 \in R^{m \times \hat{r}}\)
Step 3 : 计算 \(z = (B_1^\mathrm{T}B_1)^{-1}B_1^\mathrm{T}b\)
附 列主元 \(\mathrm{QR}\) 算法
给定矩阵 \(A \in R^{m \times n}\), 其中, \(m \ge n\). 下面的算法计算 \(r = \mathrm{rank}(A)\) 和列主元 \(\mathrm{QR}\) 分解
\[ Q^\mathrm{T}AP=\begin{bmatrix} R_{11}&R_{12} \\ O&O \end{bmatrix} \]
其中, \(R_{11} \in R^{r \times r}\) 是上三角的非奇异矩阵, 作为输出结果, \(A\) 的上三角部分存放 \(R\) 的上三角部分, 置换矩阵 \(P\) 用整数向量 \(piv\) 编码 (若 \(piv\) 的第 \(j\) 个元素等于整数 \(m\), 则 \(P\) 的第 \(j\) 列仅第 \(m\) 个元素为 1, 而其他元素皆为零)
例: 假定
\[ \begin{bmatrix} 3& 4& -1\\ 7& 4& -3\\ 2& 5& 3\\ -1& 4& 5 \end{bmatrix}, \quad b = \begin{bmatrix} 1\\ 1\\ 1\\ 1 \end{bmatrix} \]
得到 \(\mathrm{rank}(A)=2\) 和
\[ x_{\mathrm{LS}} = \begin{bmatrix} 0.0815\\ 0.1545\\ 0.0730 \end{bmatrix} \]
应用算法 低秩 \(\mathrm{LS}\) 方法得到
\[ P = \begin{bmatrix} 0& 1& 0\\ 1& 0& 0\\ 0& 0& 1 \end{bmatrix}, \quad x = \begin{bmatrix} 0.0845\\ 0.2275\\ 0.0000 \end{bmatrix} \]
三、奇异值分解的 \(\mathrm{QR}\) 分解算法
通常是两个阶段:
第一阶段为矩阵的二重对角化, 通过 \(\mathrm{Householder}\) 变换将矩阵 \(A_{m \times n}\) 变换为二重对角矩阵 (除对角线及其上面一条对角线的元素外, 其他元素全为零).
第二阶段利用 \(\mathrm{QR}\) 分解, 保存二重对角矩阵的形式不变, 利用正交变换使上一条对角线的元素逐渐减小, 使矩阵接近对角矩阵. 此时, 需要具体构造正交矩阵 \(V\), 但 \(U\) 可以不予保留.
和之前在线性代数中学习的人工进行 \(SVD\) 分解不同, 这里的算法更适合于计算机.
有点难度, 需要慢慢消化.
(...... 未完待续)
四、奇异值分解的精确计算
初始矩阵 \(A\) 的元素通过观察或计算得到, 存在一定的误差. 为了使得奇异值和奇异向量的计算结果是精确的, 提出了右边 \(\mathrm{Jacobi}\) 旋转.
(...... 未完待续)