矩阵分解¶

矩阵分解将复杂的矩阵分解成更简单的因子，用于求解方程组、计算逆矩阵和压缩数据。本文涵盖高斯消元、LU分解、QR分解、Cholesky分解、特征分解和奇异值分解（SVD）——这些是主成分分析（PCA）、推荐系统和机器学习中数值稳定性背后的算法。

矩阵分解（或因子分解）将一个矩阵分解成更易于处理的简单部分。可以把它类比为整数的因子分解：\(12 = 3 \times 4\) 比单独的 12 更容易理解。
我们分解矩阵是为了更快地求解方程组、稳定地计算逆矩阵、求特征值、压缩数据以及理解变换的几何意义。
最基本的技术是高斯消元法（行化简）。其思想很简单：给定一个方程组 \(A\mathbf{x} = \mathbf{b}\)，使用三种允许的操作简化 \(A\)，直到答案显而易见。
允许的操作是：交换两行，将某一行乘以一个非零标量，或将一行的倍数加到另一行上。
例如，为了消去主元下方的第一列，从下面的行中减去第一行的倍数：

\[ \begin{bmatrix} 2 & 1 & 5 \\ 4 & 3 & 7 \\ 6 & 5 & 9 \end{bmatrix} \xrightarrow{R_2 - 2R_1} \begin{bmatrix} 2 & 1 & 5 \\ 0 & 1 & -3 \\ 6 & 5 & 9 \end{bmatrix} \xrightarrow{R_3 - 3R_1} \begin{bmatrix} 2 & 1 & 5 \\ 0 & 1 & -3 \\ 0 & 2 & -6 \end{bmatrix} \]

目标是行阶梯形：每个主元（每行第一个非零条目）下方全为零，且每个主元位于它上面主元的右侧。矩阵变成阶梯形状。

高斯消元：行操作产生三角形形式，然后自底向上求解

进一步化为简化行阶梯形，使每个主元为 1 且是其列中唯一的非零条目。每个矩阵都有唯一的简化行阶梯形。
一旦变成三角形形式，我们通过回代求解：最下面一行直接给出最后一个变量，然后向上回代。
这是所有其他分解所构建的基础，分解的目标就是将矩阵化简为三角形形式，这样我们就可以回代并求解变量。
LU分解 通过将方阵分解为 \(A = LU\)（或包含行交换的 \(A = PLU\)）来形式化高斯消元，其中 \(L\) 是下三角矩阵，\(U\) 是上三角矩阵。

LU分解：将一个困难的矩阵拆分成两个简单的三角矩阵

求解 \(A\mathbf{x} = \mathbf{b}\)：首先通过前代（自上而下）求解 \(L\mathbf{y} = \mathbf{b}\)，然后通过回代（自下而上）求解 \(U\mathbf{x} = \mathbf{y}\)。用两次三角求解代替一次困难的一般求解。
相对于原始高斯消元的优势在于可重用性。一旦你得到了 \(L\) 和 \(U\)，就可以针对许多不同的 \(\mathbf{b}\) 向量求解，而无需重新进行分解。
如果你需要用 1000 个不同的右端项求解同一个方程组（在仿真中很常见），那么只需分解一次，然后重复使用。
当矩阵是对称且正定的（如协方差矩阵）时，我们可以做得更好。
Cholesky分解 将其分解为 \(A = LL^T\)，其中 \(L\) 是下三角矩阵。例如：

\[ \begin{bmatrix} 4 & 2 \\ 2 & 5 \end{bmatrix} = \begin{bmatrix} 2 & 0 \\ 1 & 2 \end{bmatrix} \begin{bmatrix} 2 & 1 \\ 0 & 2 \end{bmatrix} \]

这大约比 LU 快两倍，并且保证数值稳定。可以把它看作矩阵的“平方根”。
如果分解失败（平方根下出现负值），则矩阵不是正定的。因此 Cholesky 也兼作正定性的检验。
方阵 \(A\) 的特征向量是那些特殊的向量方向，变换仅沿这些方向拉伸或收缩它们，而不旋转。特征值是缩放因子：

\[A\mathbf{x} = \lambda\mathbf{x}\]

特征向量保持在同一条直线上（仅缩放），普通向量被旋转

大多数向量被矩阵相乘后方向会改变。但特征向量是特殊的：输出方向与输入方向相同，仅缩放 \(\lambda\) 倍。如果 \(\lambda = 2\)，特征向量长度加倍。如果 \(\lambda = -1\)，方向翻转。如果 \(\lambda = 0\)，则被压缩为零。
例如，对于

\[ A = \begin{bmatrix} 3 & 1 \\ 0 & 2 \end{bmatrix} \]

向量 \([1, 0]^T\) 是特征值为 \(\lambda = 3\) 的特征向量，因为 \(A[1, 0]^T = [3, 0]^T = 3[1, 0]^T\)。

求特征值需要解特征多项式 \(\det(A - \lambda I) = 0\)，其根即为特征值。然后将每个 \(\lambda\) 代回 \((A - \lambda I)\mathbf{x} = \mathbf{0}\) 求出对应的特征向量。
关键性质：
- \(A\) 的迹等于其特征值之和。
- \(A\) 的行列式等于其特征值之积。
- 对称矩阵的特征向量互相垂直，特征值为实数。
- 正定矩阵的所有特征值都为正。
- 协方差矩阵（我们将在统计学中遇到）总是半正定的。
通过特征多项式计算特征值对于大型矩阵是不现实的。取而代之的是使用迭代方法：
- 幂迭代：反复乘以 \(A\) 并归一化。收敛到主特征向量（最大特征值）。简单但只能找到一个特征对。
- QR算法：是主力方法。反复使用 QR 分解并重组，直到矩阵收敛到三角形形式，从而在对角线上揭示所有特征值。
- 反迭代：找到最接近给定目标值的特征向量。当你大致知道想要哪个特征值时很有用。
- 对于大型稀疏矩阵，Arnoldi 和 Lanczos 迭代利用稀疏性提高效率。
如果方阵有一组完整的线性无关的特征向量，它可以被对角化：\(A = PDP^{-1}\)，其中 \(D\) 是以特征值为对角元素的对角矩阵，\(P\) 的列是特征向量。
为什么这很有用？对角矩阵很容易处理。需要 \(A^{100}\)？无需将 \(A\) 自乘 100 次，只需计算 \(PD^{100}P^{-1}\)，而对角矩阵的幂只需要独立地对每个条目求幂。这将昂贵的操作变成了廉价的操作。
特征基 是由特征向量构成的基。在这个基下，矩阵变成对角形式，变换就是沿每个特征向量方向独立缩放。这就像是找到了变换的自然坐标系。
QR分解 将任意矩阵 \(A\) 分解为 \(A = QR\)，其中 \(Q\) 是正交矩阵（其列是标准正交的），\(R\) 是上三角矩阵。可以把它看作将“方向”信息（\(Q\)）与“缩放和混合”信息（\(R\)）分离。
Gram-Schmidt过程 逐列构建 \(Q\)。取 \(A\) 的第一列并归一化。取第二列，减去其在第一列上的投影（使其垂直），然后归一化。对每一列重复。结果是一个标准正交向量集。
QR分解是用于求特征值的QR算法背后的引擎。它也直接用于求解最小二乘问题：当 \(A\mathbf{x} = \mathbf{b}\) 没有精确解（方程数多于未知数）时，QR 找到最佳的近似解。
奇异值分解（SVD） 是最通用且可以说是最重要的分解。每个矩阵（任何形状、任何秩）都有 SVD：\(A = U\Sigma V^T\)
- \(V^T\)（\(n \times n\)，正交）：旋转输入
- \(\Sigma\)（\(m \times n\)，对角）：沿正交轴缩放（奇异值，非负，降序排列）
- \(U\)（\(m \times m\)，正交）：旋转输出

SVD：任何变换 = 旋转，然后缩放，再旋转

从几何上讲，SVD 表明每一个线性变换，无论多么复杂，都只是先旋转，然后沿轴拉伸，再旋转。一个圆变成了椭圆。
奇异值（\(\sigma_1 \geq \sigma_2 \geq \ldots\)）揭示了每个方向的“重要性”。大的奇异值对应于最重要的方向。\(A\) 的秩等于非零奇异值的个数。
低秩近似：只保留最大的 \(k\) 个奇异值，将其余置零，就得到了 \(A\) 的最优秩 \(k\) 近似。这就是图像压缩的原理：一个 \(1000 \times 1000\) 的图像可能只需要 \(k = 50\) 个奇异值就能看起来几乎一样，从而将其压缩 20 倍。
SVD 也提供了伪逆：\(A^+ = V\Sigma^+U^T\)，其中 \(\Sigma^+\) 对非零奇异值取倒数。
特征分解只适用于方阵，而 SVD 适用于任何矩阵。这是它的关键优势。
主成分分析（PCA） 使用特征分解（或 SVD）进行降维。
想象一个数据集，每个样本有 100 个特征（向量维度为 100，堆叠成一个矩阵）。其中许多特征是相关的、冗余的。
PCA 找到数据实际变化的方向，让你只保留重要的部分。

PCA 找到数据中方差最大的方向

第一主成分（PC1）是方差最大的方向。
第二主成分（PC2）捕捉剩余方差中最大的部分，并且垂直于第一主成分。
如果大部分方差只沿着少数几个方向分布，你就可以将数据投影到这些维度上，丢弃其余的维度，而损失极小。
步骤：
- 标准化数据（减去均值，除以标准差），使所有特征贡献均等
- 计算协方差矩阵
- 求其特征值和特征向量
- 选择特征值最大的 \(k\) 个特征向量（这些就是主成分）
- 将数据投影到这些主成分上
标准化至关重要：没有它，以公里为单位的特征将主导以厘米为单位的特征，无论其实际重要性如何。
在实践中，PCA 用于可视化（将高维数据投影到二维或三维）、降噪（丢弃主要是噪声的低方差方向）以及通过减少输入特征数量来加速机器学习模型。
核PCA 将 PCA 扩展到非线性关系。它通过核函数将数据映射到更高维的空间，在那里结构变为线性，然后应用标准 PCA 并投影回来。
Schur分解 将方阵分解为 \(A = QTQ^\ast\)，其中 \(Q\) 是酉矩阵，\(T\) 是上三角矩阵。每个方阵都有 Schur 分解，即使它不能被对角化。
非负矩阵分解（NMF） 将一个矩阵分解为两个非负矩阵：\(A \approx WH\)，其中 \(W\) 和 \(H\) 中的所有条目都 \(\geq 0\)。与可能产生负条目的 SVD 不同，NMF 只做加法，从不做减法。这使得分解部分具有可解释性：在主题建模中，\(W\) 给出每个文档的主题权重，\(H\) 给出每个主题的词权重，所有值均为非负，这与我们对文档包含“每个主题多少”的认知一致。
谱定理 指出对称（或 Hermitian）矩阵总是可以用正交（或酉）矩阵对角化。它们的特征值总是实数，特征向量总是正交。这是 PCA 背后的理论基础。

编程任务（使用 CoLab 或 notebook）¶

计算对称矩阵的特征值和特征向量。验证特征向量互相垂直，并从特征分解重建矩阵。

import jax.numpy as jnp

A = jnp.array([[4.0, 2.0],
               [2.0, 3.0]])

eigenvalues, eigenvectors = jnp.linalg.eigh(A)
print(f"特征值: {eigenvalues}")
print(f"特征向量正交: {jnp.dot(eigenvectors[:,0], eigenvectors[:,1]):.6f}")

# 重建：A = P D P^T
D = jnp.diag(eigenvalues)
A_reconstructed = eigenvectors @ D @ eigenvectors.T
print(f"重建匹配: {jnp.allclose(A, A_reconstructed)}")

实现幂迭代求最大特征值，以及反迭代求最小特征值。与 jnp.linalg.eigh 比较。然后尝试自己实现 QR 算法。

import jax.numpy as jnp

A = jnp.array([[4.0, 2.0],
               [2.0, 3.0]])

# 幂迭代：求最大的特征值
v = jnp.array([1.0, 0.0])
for _ in range(20):
    v = A @ v
    v = v / jnp.linalg.norm(v)
print(f"最大特征值:  {v @ A @ v:.4f}")

# 反迭代：乘以 A^{-1} 而不是 A，求最小的特征值
v = jnp.array([1.0, 0.0])
for _ in range(20):
    v = jnp.linalg.solve(A, v)
    v = v / jnp.linalg.norm(v)
print(f"最小特征值: {1.0 / (v @ jnp.linalg.solve(A, v)):.4f}")

print(f"jnp.linalg.eigh:    {jnp.linalg.eigh(A)[0]}")

计算矩阵的 SVD，然后仅使用前 k 个奇异值重建矩阵，观察近似质量随 k 的变化。

import jax.numpy as jnp

A = jnp.array([[1.0, 2.0, 3.0],
               [4.0, 5.0, 6.0],
               [7.0, 8.0, 9.0]])

U, S, Vt = jnp.linalg.svd(A)

for k in [1, 2, 3]:
    approx = U[:, :k] @ jnp.diag(S[:k]) @ Vt[:k, :]
    error = jnp.linalg.norm(A - approx)
    print(f"k={k}, 重建误差: {error:.4f}")