矩阵分解技术解析

矩阵分解技术解析

矩阵分解是线性代数中最核心的技术之一，它将复杂的矩阵运算转化为更简单、更易处理的形式。这项技术不仅在理论数学中占据重要地位，更是现代机器学习、数据科学和工程计算的基础工具。

矩阵分解的本质与意义

矩阵分解的核心思想是将一个复杂的矩阵 $\mathbf{A}$ 表示为若干个具有特殊性质的矩阵的乘积。这种分解揭示了原始矩阵的内在结构，使得许多困难的计算变得简单高效。就像将一个复合函数分解为基本函数的组合，矩阵分解使我们能够通过理解各个组成部分来理解整体。在实际应用中，矩阵分解的价值体现在多个方面。首先，它能够大幅降低计算复杂度——例如求解线性方程组 $\mathbf{Ax}=\mathbf{b}$，直接求解需要 $O(n^3)$ 的计算量，而通过 LU 分解后只需要 $O(n^2)$ 的前向和后向替换。其次，矩阵分解能够提高数值稳定性，避免直接计算中的数值误差累积。最重要的是，不同的分解方法揭示了矩阵的不同性质，为各种应用提供了理论基础。

奇异值分解（SVD）：最强大的分解工具

奇异值分解可以说是所有矩阵分解技术中最重要、最通用的一种。对于任意 $m\times n$ 矩阵 $\mathbf{A}$，无论是方阵还是矩形矩阵，实数还是复数矩阵，都存在奇异值分解：

$$\mathbf{A}={\mathbf{U\Sigma V}}^T$$

其中 $\mathbf{U}$ 是 $m\times m$ 的正交矩阵，包含左奇异向量；$\mathbf{\Sigma }$ 是 $m\times n$ 的对角矩阵，对角线上的元素称为奇异值，按照 ${\sigma }_1\ge {\sigma }_2\ge \dots \ge {\sigma }_r\ge 0$ 的顺序排列；${\mathbf{V}}^T$ 是 $n\times n$ 正交矩阵的转置，包含右奇异向量。

SVD 的几何直观

从几何角度理解，SVD 将任何线性变换分解为三个基本操作的组合。想象一个二维平面上的单位圆，经过矩阵 $\mathbf{A}$ 的变换后会变成一个椭圆。SVD 告诉我们这个变换可以分解为：首先通过 ${\mathbf{V}}^T$ 进行旋转或反射，将坐标系对齐到特定方向；然后通过 $\mathbf{\Sigma }$ 沿着坐标轴进行拉伸或压缩，拉伸的比例就是奇异值；最后通过 $\mathbf{U}$ 再次旋转或反射到最终位置。这种几何理解对于理解主成分分析（PCA）特别有帮助。在高维数据空间中，SVD 能够找到数据变化最大的方向（主成分），这些方向就是右奇异向量，而对应的奇异值表示了在该方向上的变化程度。

计算 SVD 的算法

计算 SVD 的标准方法是先计算 ${\mathbf{A}}^T\mathbf{A}$ 的特征值和特征向量得到 $\mathbf{V}$，再计算 ${\mathbf{AA}}^T$ 的特征值和特征向量得到 $\mathbf{U}$，奇异值则是 ${\mathbf{A}}^T\mathbf{A}$ 特征值的平方根。然而在实际计算中，这种方法会导致数值精度损失，因为 ${\mathbf{A}}^T\mathbf{A}$ 的条件数是 $\mathbf{A}$ 条件数的平方。现代算法采用更稳定的方法，如 Golub-Kahan 双对角化算法。该算法首先将矩阵化为双对角形式，然后使用 QR 迭代或分治算法计算奇异值。整个过程的时间复杂度为 $O(\min (m{n}^2,m^{2}n))$，对于方阵就是 $O(n^3)$。

特征值分解：揭示矩阵的本质特性

特征值分解是理解方阵性质的基本工具。对于 $n\times n$ 方阵 $\mathbf{A}$，如果存在 $n$ 个线性无关的特征向量，则可以分解为：

$$\mathbf{A}={\mathbf{V\Lambda V}}^{-1}$$

其中 $\mathbf{V}$ 的列是特征向量，$\mathbf{\Lambda }$ 是对角矩阵，对角线元素是对应的特征值。对于实对称矩阵，特征向量是正交的，分解简化为 $\mathbf{A}={\mathbf{Q\Lambda Q}}^T$。

特征值的物理意义

特征值和特征向量有着深刻的物理意义。特征向量代表了变换的不变方向——沿着特征向量方向的向量在变换后仍然保持同一方向，只是长度发生了变化，变化的倍数就是对应的特征值。在振动分析中，特征向量对应振动模态，特征值对应振动频率的平方；在量子力学中，特征向量是量子态，特征值是可观测量的测量值。其中矩阵的许多重要性质都可以通过特征值来表达。矩阵的迹（对角线元素之和）等于所有特征值之和：$\text{tr}(\mathbf{A})={\sum }_i{\lambda }_i$；行列式等于所有特征值之积：$\det (\mathbf{A})={\prod }_i{\lambda }_i$；矩阵的秩等于非零特征值的个数。这些关系在理论分析和实际计算中都非常有用。

计算特征值的数值方法

理论上，特征值是特征多项式 $\det (\mathbf{A}-\lambda \mathbf{I})=0$ 的根，但直接求解高次方程在数值上是不稳定的。实际计算中使用迭代方法，其中最重要的是 QR 算法。QR 算法的基本思想是通过相似变换将矩阵化为上三角形式（Schur 形式），对角线元素就是特征值。算法步骤是：将 $\mathbf{A}$ 进行 QR 分解得到 $\mathbf{A}={\mathbf{Q}}_1{\mathbf{R}}_1$，然后计算 ${\mathbf{A}}_1={\mathbf{R}}_1{\mathbf{Q}}_1$，重复这个过程。在适当的条件下，序列 {Ak}\{\mathbf{A}_k\}{Ak​} 收敛到上三角矩阵。为了加速收敛，实际实现中会先将矩阵化为 Hessenberg 形式，并使用位移策略。对于大规模稀疏矩阵，Lanczos 算法和 Arnoldi 算法更加高效。这些方法基于 Krylov 子空间，只计算部分特征值和特征向量，特别适合只需要最大或最小几个特征值的情况。

LU 分解：求解线性方程组的主力

LU 分解将方阵分解为下三角矩阵 $\mathbf{L}$ 和上三角矩阵 $\mathbf{U}$ 的乘积：

$$\mathbf{PA}=\mathbf{LU}$$

其中 $\mathbf{P}$ 是置换矩阵，用于数值稳定性。LU 分解本质上是高斯消元法的矩阵形式，$\mathbf{L}$ 记录了消元过程中的乘数，$\mathbf{U}$ 是消元后的上三角矩阵。

LU 分解的计算过程

Doolittle 算法是计算 LU 分解的标准方法。算法逐列处理矩阵，对于第 $k$ 列，首先选择主元（部分选主元策略选择该列绝对值最大的元素），然后计算 $\mathbf{L}$ 的第 $k$ 列和 $\mathbf{U}$ 的第 $k$ 行。具体计算公式为：

对于 $\mathbf{U}$ 的元素：
$$u_{ij}=a_{ij}-\sum _{k=1}^{i-1}{l}_{ik}{u}_{kj}$$

对于 $\mathbf{L}$ 的元素：
$$l_{ij}=\frac{a_{ij}-{\sum }_{k=1}^{j-1}{l}_{ik}{u}_{kj}}{u_{jj}}$$

整个算法的时间复杂度为 $O(2n^3/3)$，约为矩阵乘法的 $2/3$。

数值稳定性与选主元策略

不使用选主元的 LU 分解可能在数值上不稳定。考虑矩阵：
$$\left[\begin{array}{cc}ϵ& 1\\ 1& 1\end{array}\right]$$

当 $ϵ$ 很小时，直接 LU 分解会产生很大的数值误差。部分选主元通过交换行来确保 $|l_{ij}|\le 1$，大大提高了数值稳定性。理论上，增长因子（衡量数值稳定性的指标）最坏情况下可达 $2^{n-1}$，但实际问题中很少出现如此极端的情况。

求解线性方程组

有了 LU 分解，求解 $\mathbf{Ax}=\mathbf{b}$ 变成两步：首先求解 $\mathbf{Ly}=\mathbf{Pb}$（前向替换），然后求解 $\mathbf{Ux}=\mathbf{y}$（后向替换）。每步只需要 $O(n^2)$ 的运算，远比直接求逆矩阵高效。当需要求解多个右端向量时，只需进行一次 LU 分解，然后对每个右端向量进行前向和后向替换，这在许多工程应用中非常有用。

QR 分解：数值稳定性的典范

QR 分解将矩阵分解为正交矩阵 $\mathbf{Q}$ 和上三角矩阵 $\mathbf{R}$ 的乘积：

$$\mathbf{A}=\mathbf{QR}$$

其中 ${\mathbf{Q}}^T\mathbf{Q}=\mathbf{I}$，保证了良好的数值性质。QR 分解在最小二乘问题、特征值计算等领域有着重要应用。

Gram-Schmidt 正交化与数值稳定性

经典的 Gram-Schmidt 过程通过逐个正交化列向量来构造 $\mathbf{Q}$：对于第 $k$ 个列向量，减去它在前 $k-1$ 个正交向量上的投影，然后归一化。然而，由于舍入误差的累积，经典 Gram-Schmidt 在数值上是不稳定的。改进的 Gram-Schmidt 算法通过改变计算顺序提高了稳定性，但最稳定的方法是 Householder 变换。Householder 变换使用反射矩阵 $\mathbf{H}=\mathbf{I}-2{\mathbf{vv}}^T$ 来消去向量的特定分量。每次变换将一列的子对角元素化为零，$n-1$ 次变换后得到上三角矩阵 $\mathbf{R}$。

QR 分解在最小二乘问题中的应用

对于超定方程组 $\mathbf{Ax}\approx \mathbf{b}$（$m>n$），最小二乘解满足正规方程 ${\mathbf{A}}^T\mathbf{Ax}={\mathbf{A}}^T\mathbf{b}$。直接求解正规方程会使条件数平方，降低数值稳定性。使用 QR 分解，问题转化为：

$$\Vert \mathbf{Ax}-\mathbf{b}{\Vert }^2=\Vert \mathbf{QRx}-\mathbf{b}{\Vert }^2=\Vert \mathbf{Rx}-{\mathbf{Q}}^T\mathbf{b}{\Vert }^2$$

由于 $\mathbf{Q}$ 是正交矩阵，不改变向量的范数。因此只需求解上三角方程组 $\mathbf{Rx}={\mathbf{Q}}^T\mathbf{b}$ 的前 $n$ 个方程即可。这种方法避免了形成 ${\mathbf{A}}^T\mathbf{A}$，保持了原问题的条件数，在数值上更加稳定。

QR 算法计算特征值

QR 算法是计算矩阵全部特征值的标准方法。基本思想是通过一系列相似变换将矩阵化为上三角形式。算法的每一步进行 QR 分解 ${\mathbf{A}}_k={\mathbf{Q}}_k{\mathbf{R}}_k$，然后形成 ${\mathbf{A}}_{k+1}={\mathbf{R}}_k{\mathbf{Q}}_k$。由于 ${\mathbf{A}}_{k+1}={\mathbf{Q}}_k^T{\mathbf{A}}_k{\mathbf{Q}}_k$，保持了特征值不变。在适当条件下，{Ak}\{\mathbf{A}_k\}{Ak​} 收敛到上三角矩阵（实 Schur 形式），对角线元素就是特征值。实际实现中，先将矩阵化为 Hessenberg 形式可以大大减少计算量，每步 QR 分解只需要 $O(n^2)$ 而不是 $O(n^3)$。

Cholesky 分解：对称正定矩阵的专属方法

对于对称正定矩阵，Cholesky 分解提供了最高效的分解方法：

$$\mathbf{A}={\mathbf{LL}}^T$$

其中 $\mathbf{L}$ 是下三角矩阵，对角线元素为正。Cholesky 分解可以看作对称正定矩阵的"平方根"。

计算 Cholesky 分解

Cholesky 分解的计算相对简单。对角线元素通过下式计算，确保了实数和正值：
$$l_{ii}=\sqrt{a_{ii}-\sum _{k=1}^{i-1}{l}_{ik}^2}$$

非对角线元素通过下式计算：
$$l_{ij}=\frac{a_{ij}-{\sum }_{k=1}^{j-1}{l}_{ik}{l}_{jk}}{l_{jj}}$$

算法只需要处理下三角部分，时间复杂度为 $O(n^3/3)$，约为 LU 分解的一半。空间上也只需要存储下三角部分，节省了一半的内存。

正定性的几何意义

对称正定矩阵定义了一个椭球体：${\mathbf{x}}^T\mathbf{Ax}=1$。Cholesky 分解 $\mathbf{A}={\mathbf{LL}}^T$ 将这个椭球体的描述分解为两步变换：首先通过 ${\mathbf{L}}^T$ 将标准球面变换到新的坐标系，然后通过 $\mathbf{L}$ 进行相同的变换。这保证了椭球体的所有主轴长度都是正的，对应于正定性的要求。在概率论中，多元正态分布的协方差矩阵总是对称正定的。Cholesky 分解在生成多元正态随机数时特别有用：如果 $\mathbf{z}$ 是标准正态随机向量，那么 $\mathbf{x}=\mathbf{Lz}+\mathbf{\mu }$ 服从均值为 $\mathbf{\mu }$、协方差为 $\mathbf{A}={\mathbf{LL}}^T$ 的多元正态分布。

在优化算法中的应用

在数值优化中，许多算法需要求解形如 $\mathbf{Hx}=-\mathbf{g}$ 的方程，其中 $\mathbf{H}$ 是 Hessian 矩阵。对于凸优化问题，Hessian 矩阵是对称正定的，使用 Cholesky 分解最为高效。牛顿法的每次迭代需要求解这样的线性方程组。通过 Cholesky 分解，不仅计算效率高，还能检测 Hessian 是否正定——如果分解过程中出现负的对角元素，说明当前点不是严格凸的，需要修正算法。

其他重要的矩阵分解

Schur 分解

Schur 分解是特征值分解的推广，对任意方阵都存在：

$$\mathbf{A}={\mathbf{QTQ}}^H$$

其中 $\mathbf{Q}$ 是酉矩阵，$\mathbf{T}$ 是上三角矩阵，对角线元素是特征值。与 Jordan 标准形相比，Schur 分解在数值上更稳定，是许多特征值算法的基础。

实 Schur 分解特别重要，它保证了实矩阵的分解结果也是实的，其中 $\mathbf{T}$ 是准上三角矩阵（对角线上可能有 $2\times 2$ 的块对应复共轭特征值对）。

极分解

极分解将矩阵分解为酉矩阵和半正定矩阵的乘积：

$$\mathbf{A}=\mathbf{UP}$$

其中 $\mathbf{U}$ 是酉矩阵（代表旋转），$\mathbf{P}$ 是半正定 Hermitian 矩阵（代表拉伸）。这种分解在计算机图形学中用于分离刚体变换和形变，在力学中用于分析应变张量。

极分解与 SVD 密切相关：如果 $\mathbf{A}={\mathbf{U\Sigma V}}^T$ 是 SVD，那么极分解为 $\mathbf{A}=({\mathbf{UV}}^T)({\mathbf{V\Sigma V}}^T)$。

矩阵分解方法的选择策略

选择合适的矩阵分解方法需要考虑多个因素。首先是问题的数学性质：矩阵是否对称？是否正定？是否稀疏？条件数如何？其次是计算需求：需要全部特征值还是部分？是一次性求解还是多次求解？对精度和速度的要求如何？

对于求解线性方程组 $\mathbf{Ax}=\mathbf{b}$，如果 $\mathbf{A}$ 是对称正定的，Cholesky 分解是最佳选择，速度最快且数值稳定。对于一般的方阵，LU 分解配合部分选主元是标准方法。如果矩阵接近奇异或条件数很大，应该使用 QR 分解或 SVD 以获得更好的数值稳定性。

对于最小二乘问题 $\min \Vert \mathbf{Ax}-\mathbf{b}{\Vert }_2$，QR 分解是首选方法，它在保持数值稳定性的同时有较好的计算效率。当矩阵严重病态或需要正则化时，SVD 提供了最稳定的解决方案，可以方便地实现 Tikhonov 正则化。

在特征值问题中，对于中小规模的密集矩阵，QR 算法计算全部特征值；对于大规模稀疏矩阵，Lanczos 或 Arnoldi 方法计算部分特征值更高效。对称矩阵的特征值计算通常比一般矩阵更快更稳定。

数值稳定性的深入理解

数值稳定性是选择和实现矩阵分解算法的关键考虑因素。条件数 $\kappa (\mathbf{A})=\Vert \mathbf{A}\Vert \Vert {\mathbf{A}}^{-1}\Vert$ 衡量了矩阵对扰动的敏感程度。当 $\kappa (\mathbf{A})\gg 1$ 时，矩阵是病态的，小的输入误差会导致大的输出误差。

SVD 直接提供了条件数的计算：$\kappa (\mathbf{A})={\sigma }_{\max }/{\sigma }_{\min }$。这也解释了为什么 SVD 在处理病态问题时特别有用——它能够识别并处理接近零的奇异值，避免数值问题。

向后稳定性是评价算法的重要标准。一个向后稳定的算法产生的结果等价于对略微扰动的输入进行精确计算。Householder QR 分解是向后稳定的典范，而经典 Gram-Schmidt 则不是。这解释了为什么实际应用中优先选择 Householder 方法。

在实现中，选主元策略对提高稳定性至关重要。部分选主元确保 LU 分解中的乘数不会过大，完全选主元虽然理论上更稳定但计算成本较高。对于对称矩阵，对角选主元可以保持对称性同时提高稳定性。

附录：核心定理的详细推导

A. 奇异值分解的存在性与构造

定理 A.1（SVD 存在性定理）：对于任意 $m\times n$ 的复矩阵 $\mathbf{A}$，存在 $m\times m$ 酉矩阵 $\mathbf{U}$，$n\times n$ 酉矩阵 $\mathbf{V}$，以及 $m\times n$ 对角矩阵 $\mathbf{\Sigma }$，使得：
$$\mathbf{A}={\mathbf{U\Sigma V}}^{\ast }$$

证明：

首先考虑 $n\times n$ Hermitian 半正定矩阵 ${\mathbf{A}}^{\ast }\mathbf{A}$。由谱定理，存在酉矩阵 $\mathbf{V}$ 和非负实对角矩阵 $\mathbf{D}$，使得：
$${\mathbf{A}}^{\ast }\mathbf{A}={\mathbf{VDV}}^{\ast }$$

设 $\mathbf{D}=\text{diag}({\lambda }_1,{\lambda }_2,\dots ,{\lambda }_n)$，其中 ${\lambda }_1\ge {\lambda }_2\ge \dots \ge {\lambda }_r>0={\lambda }_{r+1}=\dots ={\lambda }_n$，$r=\text{rank}(\mathbf{A})$。

定义奇异值 ${\sigma }_i=\sqrt{{\lambda }_i}$ 对于 $i=1,\dots ,r$。构造 $m\times n$ 矩阵：
$$\mathbf{\Sigma }=\left[\begin{array}{cc}\text{diag}({\sigma }_1,\dots ,{\sigma }_r)& 0\\ 0& 0\end{array}\right]$$

令 ${\mathbf{v}}_1,\dots ,{\mathbf{v}}_n$ 为 $\mathbf{V}$ 的列向量（${\mathbf{A}}^{\ast }\mathbf{A}$ 的特征向量）。对于 $i=1,\dots ,r$，定义：
$${\mathbf{u}}_i=\frac{1}{{\sigma }_i}{\mathbf{Av}}_i$$

需要验证 {u1,…,ur}\{\mathbf{u}_1, \ldots, \mathbf{u}_r\}{u1​,…,ur​} 是标准正交的。对于 $i,j\le r$：
$${\mathbf{u}}_i^{\ast }{\mathbf{u}}_j=\frac{1}{{\sigma }_i{\sigma }_j}{\mathbf{v}}_i^{\ast }{\mathbf{A}}^{\ast }{\mathbf{Av}}_j=\frac{1}{{\sigma }_i{\sigma }_j}{\mathbf{v}}_i^{\ast }({\lambda }_j{\mathbf{v}}_j)=\frac{{\lambda }_j}{{\sigma }_i{\sigma }_j}{\mathbf{v}}_i^{\ast }{\mathbf{v}}_j=\frac{{\sigma }_j^2}{{\sigma }_i{\sigma }_j}{\delta }_{ij}={\delta }_{ij}$$

若 $r<m$，通过 Gram-Schmidt 过程扩展 {u1,…,ur}\{\mathbf{u}_1, \ldots, \mathbf{u}_r\}{u1​,…,ur​} 到 ${\mathbb{C}}^m$ 的完整标准正交基 {u1,…,um}\{\mathbf{u}_1, \ldots, \mathbf{u}_m\}{u1​,…,um​}。

构造 $\mathbf{U}=[{\mathbf{u}}_1|\dots |{\mathbf{u}}_m]$。现在验证 $\mathbf{A}={\mathbf{U\Sigma V}}^{\ast }$：

对于 $i\le r$：
$${\mathbf{Av}}_i={\sigma }_i{\mathbf{u}}_i=(\mathbf{U\Sigma }{)}_i$$

对于 $i>r$：
$${\mathbf{Av}}_i=0=(\mathbf{U\Sigma }{)}_i$$

因此 $\mathbf{AV}=\mathbf{U\Sigma }$，即 $\mathbf{A}={\mathbf{U\Sigma V}}^{\ast }$。

B. Eckart-Young-Mirsky 定理：最优低秩近似

定理 B.1：设 $\mathbf{A}={\mathbf{U\Sigma V}}^T$ 是 $m\times n$ 矩阵的 SVD，其中奇异值为 ${\sigma }_1\ge {\sigma }_2\ge \dots \ge {\sigma }_{\min (m,n)}$。对于 $k<\text{rank}(\mathbf{A})$，定义：
$${\mathbf{A}}_k=\sum _{i=1}^k{\sigma }_i{\mathbf{u}}_i{\mathbf{v}}_i^T$$

则对于任意秩不超过 $k$ 的矩阵 $\mathbf{B}$：
$$\Vert \mathbf{A}-{\mathbf{A}}_k{\Vert }_F\le \Vert \mathbf{A}-\mathbf{B}{\Vert }_F$$

其中 $\Vert \cdot {\Vert }_F$ 是 Frobenius 范数。

证明：

首先注意到：
$$\Vert \mathbf{A}-{\mathbf{A}}_k{\Vert }_F^2={\Vert \sum _{i=k+1}^r{\sigma }_i{\mathbf{u}}_i{\mathbf{v}}_i^T\Vert }_F^2=\sum _{i=k+1}^r{\sigma }_i^2$$

这是因为 {uiviT}\{\mathbf{u}_i\mathbf{v}_i^T\}{ui​viT​} 构成正交基。

设 $\mathbf{B}$ 是任意秩不超过 $k$ 的矩阵。存在维数至少为 $n-k$ 的子空间 $\mathcal{V}\subseteq {\mathbb{R}}^n$ 使得 $\mathbf{Bx}=0$ 对所有 $\mathbf{x}\in \mathcal{V}$。

考虑由 {v1,…,vk+1}\{\mathbf{v}_1, \ldots, \mathbf{v}_{k+1}\}{v1​,…,vk+1​} 张成的 $(k+1)$ 维子空间 $\mathcal{W}$。由于 $\dim (\mathcal{V})+\dim (\mathcal{W})=(n-k)+(k+1)=n+1>n$，必存在非零向量 $\mathbf{x}\in \mathcal{V}\cap \mathcal{W}$。

不失一般性，设 $\Vert \mathbf{x}\Vert =1$。由于 $\mathbf{x}\in \mathcal{W}$，可以写成：
$$\mathbf{x}=\sum _{i=1}^{k+1}{\alpha }_i{\mathbf{v}}_i,\sum _{i=1}^{k+1}{\alpha }_i^2=1$$

因此：
$$\mathbf{Ax}=\sum _{i=1}^{k+1}{\alpha }_i{\sigma }_i{\mathbf{u}}_i$$

由于 $\mathbf{x}\in \mathcal{V}$，有 $\mathbf{Bx}=0$，所以：
$$\Vert (\mathbf{A}-\mathbf{B})\mathbf{x}{\Vert }^2=\Vert \mathbf{Ax}{\Vert }^2=\sum _{i=1}^{k+1}{\alpha }_i^2{\sigma }_i^2\ge {\sigma }_{k+1}^2\sum _{i=1}^{k+1}{\alpha }_i^2={\sigma }_{k+1}^2$$

由 Frobenius 范数的变分特征：
$$\Vert \mathbf{A}-\mathbf{B}{\Vert }_F^2\ge {\sigma }_{k+1}^2+{\sigma }_{k+2}^2+\dots +{\sigma }_r^2=\Vert \mathbf{A}-{\mathbf{A}}_k{\Vert }_F^2$$

这完成了证明。

C. QR 分解的 Householder 构造

定理 C.1：对于任意 $m\times n$ 矩阵 $\mathbf{A}$（$m\ge n$），存在正交矩阵 $\mathbf{Q}$ 和上三角矩阵 $\mathbf{R}$ 使得 $\mathbf{A}=\mathbf{QR}$。

构造性证明：

Householder 反射器定义为：
$$\mathbf{H}=\mathbf{I}-2\frac{{\mathbf{vv}}^T}{{\mathbf{v}}^T\mathbf{v}}$$

其中 $\mathbf{v}\ne 0$ 是反射向量。$\mathbf{H}$ 是对称正交矩阵，满足 ${\mathbf{H}}^2=\mathbf{I}$。

给定向量 $\mathbf{x}\in {\mathbb{R}}^m$，要找到 $\mathbf{H}$ 使得 $\mathbf{Hx}=\alpha {\mathbf{e}}_1$，其中 ${\mathbf{e}}_1$ 是第一个标准基向量，$\alpha =\pm \Vert \mathbf{x}\Vert$。

设 $\mathbf{v}=\mathbf{x}-\alpha {\mathbf{e}}_1$。则：
$$\mathbf{Hx}=\mathbf{x}-2\frac{{\mathbf{v}}^T\mathbf{x}}{{\mathbf{v}}^T\mathbf{v}}\mathbf{v}$$

计算：
$${\mathbf{v}}^T\mathbf{x}=(\mathbf{x}-\alpha {\mathbf{e}}_1{)}^T\mathbf{x}=\Vert \mathbf{x}{\Vert }^2-\alpha x_1$$
$${\mathbf{v}}^T\mathbf{v}=\Vert \mathbf{x}{\Vert }^2-2\alpha x_1+{\alpha }^2=2({\alpha }^2-\alpha x_1)$$

因此：
$$\frac{{\mathbf{v}}^T\mathbf{x}}{{\mathbf{v}}^T\mathbf{v}}=\frac{{\alpha }^2-\alpha x_1}{2({\alpha }^2-\alpha x_1)}=\frac{1}{2}$$

所以：
$$\mathbf{Hx}=\mathbf{x}-\mathbf{v}=\alpha {\mathbf{e}}_1$$

为了数值稳定性，选择 $\alpha =-\text{sign}(x_1)\Vert \mathbf{x}\Vert$ 以避免消去。

QR 分解算法：

对矩阵 ${\mathbf{A}}^{(0)}=\mathbf{A}$ 进行迭代。在第 $k$ 步，设：
$${\mathbf{A}}^{(k-1)}=\left[\begin{array}{cc}{\mathbf{R}}_{11}^{(k-1)}& {\mathbf{R}}_{12}^{(k-1)}\\ 0& {\mathbf{A}}_{22}^{(k-1)}\end{array}\right]$$

其中 ${\mathbf{R}}_{11}^{(k-1)}$ 是 $(k-1)\times (k-1)$ 上三角矩阵。

构造 Householder 矩阵 ${\mathbf{H}}_k$ 将 ${\mathbf{A}}_{22}^{(k-1)}$ 的第一列化为 $\alpha {\mathbf{e}}_1$：
$${\mathbf{Q}}_k=\left[\begin{array}{cc}{\mathbf{I}}_{k-1}& 0\\ 0& {\mathbf{H}}_k\end{array}\right]$$

令 ${\mathbf{A}}^{(k)}={\mathbf{Q}}_k{\mathbf{A}}^{(k-1)}$。

经过 $n$ 步后，${\mathbf{A}}^{(n)}=\mathbf{R}$ 是上三角矩阵，且：
$$\mathbf{R}={\mathbf{Q}}_n{\mathbf{Q}}_{n-1}\cdots {\mathbf{Q}}_1\mathbf{A}$$

因此：
$$\mathbf{A}={\mathbf{Q}}_1^T{\mathbf{Q}}_2^T\cdots {\mathbf{Q}}_n^T\mathbf{R}=\mathbf{QR}$$

其中 $\mathbf{Q}={\mathbf{Q}}_1^T{\mathbf{Q}}_2^T\cdots {\mathbf{Q}}_n^T$ 是正交矩阵。

D. Cholesky 分解的唯一性与稳定性

定理 D.1：若 $\mathbf{A}$ 是 $n\times n$ 对称正定矩阵，则存在唯一的下三角矩阵 $\mathbf{L}$（对角元素为正）使得 $\mathbf{A}={\mathbf{LL}}^T$。

证明：

存在性：使用归纳法。

对于 $n=1$，$\mathbf{A}=[a]$ 其中 $a>0$，取 $\mathbf{L}=[\sqrt{a}]$。

假设定理对 $(n-1)\times (n-1)$ 矩阵成立。将 $\mathbf{A}$ 分块为：
$$\mathbf{A}=\left[\begin{array}{cc}{a}_{11}& {\mathbf{a}}^T\\ \mathbf{a}& {\mathbf{A}}_{22}\end{array}\right]$$

其中 $a_{11}>0$（因为 $\mathbf{A}$ 正定），$\mathbf{a}\in {\mathbb{R}}^{n-1}$，${\mathbf{A}}_{22}$ 是 $(n-1)\times (n-1)$ 矩阵。

设：
$$\mathbf{L}=\left[\begin{array}{cc}{l}_{11}& 0^T\\ \mathbf{l}& {\mathbf{L}}_{22}\end{array}\right]$$

从 $\mathbf{A}={\mathbf{LL}}^T$ 得：
$$a_{11}=l_{11}^2\Rightarrow l_{11}=\sqrt{a_{11}}$$
$$\mathbf{a}=l_{11}\mathbf{l}\Rightarrow \mathbf{l}=\frac{1}{l_{11}}\mathbf{a}$$
$${\mathbf{A}}_{22}={\mathbf{ll}}^T+{\mathbf{L}}_{22}{\mathbf{L}}_{22}^T$$

因此：
$${\mathbf{L}}_{22}{\mathbf{L}}_{22}^T={\mathbf{A}}_{22}-{\mathbf{ll}}^T={\mathbf{A}}_{22}-\frac{1}{a_{11}}{\mathbf{aa}}^T$$

需要证明 $\mathbf{S}={\mathbf{A}}_{22}-\frac{1}{a_{11}}{\mathbf{aa}}^T$ 是正定的。对于任意 $\mathbf{x}\ne 0$：
$${\mathbf{x}}^T\mathbf{Sx}={\mathbf{x}}^T{\mathbf{A}}_{22}\mathbf{x}-\frac{1}{a_{11}}({\mathbf{a}}^T\mathbf{x}{)}^2$$

考虑向量 $\mathbf{y}=\left[\begin{array}{c}-\frac{{\mathbf{a}}^T\mathbf{x}}{a_{11}}\\ \mathbf{x}\end{array}\right]$。则：
$${\mathbf{y}}^T\mathbf{Ay}=a_{11}{\left(-\frac{{\mathbf{a}}^T\mathbf{x}}{a_{11}}\right)}^2-2\frac{{\mathbf{a}}^T\mathbf{x}}{a_{11}}{\mathbf{a}}^T\mathbf{x}+{\mathbf{x}}^T{\mathbf{A}}_{22}\mathbf{x}$$
$$=\frac{({\mathbf{a}}^T\mathbf{x}{)}^2}{a_{11}}-2\frac{({\mathbf{a}}^T\mathbf{x}{)}^2}{a_{11}}+{\mathbf{x}}^T{\mathbf{A}}_{22}\mathbf{x}={\mathbf{x}}^T\mathbf{Sx}$$

由于 $\mathbf{A}$ 正定且 $\mathbf{y}\ne 0$，故 ${\mathbf{x}}^T\mathbf{Sx}={\mathbf{y}}^T\mathbf{Ay}>0$。

因此 $\mathbf{S}$ 是正定的，由归纳假设，存在唯一的 ${\mathbf{L}}_{22}$ 使得 $\mathbf{S}={\mathbf{L}}_{22}{\mathbf{L}}_{22}^T$。

唯一性：假设 $\mathbf{A}={\mathbf{L}}_1{\mathbf{L}}_1^T={\mathbf{L}}_2{\mathbf{L}}_2^T$，其中 ${\mathbf{L}}_1,{\mathbf{L}}_2$ 都是对角元素为正的下三角矩阵。

令 $\mathbf{M}={\mathbf{L}}_1{\mathbf{L}}_2^{-1}$。则：
$${\mathbf{MM}}^T={\mathbf{L}}_1{\mathbf{L}}_2^{-1}({\mathbf{L}}_2^{-1}{)}^T{\mathbf{L}}_1^T={\mathbf{L}}_1{\mathbf{L}}_1^T({\mathbf{L}}_2^T{)}^{-1}{\mathbf{L}}_2^{-1}=\mathbf{A}{\mathbf{A}}^{-1}=\mathbf{I}$$

所以 $\mathbf{M}$ 是正交矩阵。但 $\mathbf{M}$ 同时是下三角矩阵（两个下三角矩阵之积），且对角元素为正（$m_{ii}=l_{1,ii}/l_{2,ii}>0$）。

正交下三角矩阵必为对角矩阵，且对角元素为 $\pm 1$。由于对角元素为正，故 $\mathbf{M}=\mathbf{I}$，即 ${\mathbf{L}}_1={\mathbf{L}}_2$。

E. 条件数与扰动分析

定理 E.1（线性方程组的扰动界）：设 $\mathbf{Ax}=\mathbf{b}$，$\mathbf{A}(\mathbf{x}+\delta \mathbf{x})=\mathbf{b}+\delta \mathbf{b}$。则：
$$\frac{\Vert \delta \mathbf{x}\Vert }{\Vert \mathbf{x}\Vert }\le \kappa (\mathbf{A})\frac{\Vert \delta \mathbf{b}\Vert }{\Vert \mathbf{b}\Vert }$$

其中 $\kappa (\mathbf{A})=\Vert \mathbf{A}\Vert \Vert {\mathbf{A}}^{-1}\Vert$ 是条件数。

证明：

从 $\mathbf{A}\delta \mathbf{x}=\delta \mathbf{b}$ 得：
$$\delta \mathbf{x}={\mathbf{A}}^{-1}\delta \mathbf{b}$$

因此：
$$\Vert \delta \mathbf{x}\Vert \le \Vert {\mathbf{A}}^{-1}\Vert \Vert \delta \mathbf{b}\Vert$$

从 $\mathbf{Ax}=\mathbf{b}$ 得：
$$\Vert \mathbf{b}\Vert \le \Vert \mathbf{A}\Vert \Vert \mathbf{x}\Vert$$

即：
$$\frac{1}{\Vert \mathbf{x}\Vert }\le \frac{\Vert \mathbf{A}\Vert }{\Vert \mathbf{b}\Vert }$$

结合两个不等式：
$$\frac{\Vert \delta \mathbf{x}\Vert }{\Vert \mathbf{x}\Vert }\le \Vert {\mathbf{A}}^{-1}\Vert \Vert \delta \mathbf{b}\Vert \cdot \frac{\Vert \mathbf{A}\Vert }{\Vert \mathbf{b}\Vert }=\kappa (\mathbf{A})\frac{\Vert \delta \mathbf{b}\Vert }{\Vert \mathbf{b}\Vert }$$

这个界是紧的：存在 $\mathbf{b}$ 和 $\delta \mathbf{b}$ 使等号成立。

推论 E.1：若同时考虑 $\mathbf{A}$ 和 $\mathbf{b}$ 的扰动：
$$(\mathbf{A}+\delta \mathbf{A})(\mathbf{x}+\delta \mathbf{x})=\mathbf{b}+\delta \mathbf{b}$$

且 $\Vert \delta \mathbf{A}\Vert <\Vert {\mathbf{A}}^{-1}{\Vert }^{-1}$，则：
$$\frac{\Vert \delta \mathbf{x}\Vert }{\Vert \mathbf{x}\Vert }\le \frac{\kappa (\mathbf{A})}{1-\kappa (\mathbf{A})\frac{\Vert \delta \mathbf{A}\Vert }{\Vert \mathbf{A}\Vert }}\left(\frac{\Vert \delta \mathbf{A}\Vert }{\Vert \mathbf{A}\Vert }+\frac{\Vert \delta \mathbf{b}\Vert }{\Vert \mathbf{b}\Vert }\right)$$

这表明条件数大的矩阵对扰动极其敏感，这是数值计算中必须考虑的关键因素。