矩阵分解

矩阵分解

矩阵分解把复杂的矩阵拆成更简单的因子，用于求解方程组、计算逆和数据压缩。本文件讲解高斯消元、LU、QR、Cholesky、特征值分解与 SVD——它们是 PCA、推荐系统以及 ML 中数值稳定性背后的算法。

• 矩阵分解（或因子分解）把矩阵拆成更容易处理的几块。就像把一个数分解因数：12 = 3 \times 4 比单独的 12 更好分析。

• 我们分解矩阵是为了：更快地求解方程组、稳定地计算逆、找特征值、压缩数据，以及理解变换的几何。

• 最根本的技术是高斯消元（Gaussian elimination）（行简化）。思路很简单：给定方程组 A\mathbf{x} = \mathbf{b}，用三种允许的操作把 A 化简，直到答案一目了然。

• 这三种操作是：交换两行、把某一行乘以一个非零标量、把一行的若干倍加到另一行上。

• 例如，为了消去主元下方第一列的元素，把第 1 行的若干倍从下面的行中减去：

\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}

• 目标是行阶梯形（row echelon form，REF）：每个主元（每行第一个非零项）下方全为零，且每个主元都在它上方那个主元的右侧。矩阵变成阶梯形。

• 进一步化为简化行阶梯形（reduced row echelon form，RREF），让每个主元都等于 1，并且是它所在列唯一的非零项。每个矩阵都有唯一的 RREF。

• 一旦化为三角形式，我们就用**回代（back substitution）**求解：最底下一行直接给出最后一个变量，然后向上推进。

• 这是所有其他分解所建立的根基。分解的目标就是把矩阵化为三角形式，这样我们就能回代并解出变量。

• **LU 分解（LU decomposition）**把高斯消元形式化，把方阵分解为 A = LU（或带行交换的 A = PLU），其中 L 是下三角，U 是上三角。

• 求解 A\mathbf{x} = \mathbf{b} 时：先用前代换解 L\mathbf{y} = \mathbf{b}（自顶向下），再用回代解 U\mathbf{x} = \mathbf{y}（自底向上）。两次简单的三角求解，取代了一次困难的通用求解。

• 相比原始高斯消元，它的优势在于复用。一旦你有了 L 和 U，就可以对许多不同的 \mathbf{b} 向量求解，而不必重新做分解。

• 如果你需要对 1000 个不同的右端项求解同一个方程组（在仿真中很常见），你只需分解一次然后反复复用。

• 当矩阵对称且正定时（比如协方差矩阵），我们还能做得更好。

• **Cholesky 分解（Cholesky decomposition）**把它分解为 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 的**特征向量（eigenvectors）**是那些变换后只被拉伸或压缩、而不旋转的特殊方向。**特征值（eigenvalue）**就是那个缩放因子：

• 大多数向量在乘以一个矩阵后会改变方向。但特征向量很特别：输出与输入指向同一方向，只是被 \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。

• 要找特征值，解特征多项式（characteristic polynomial） \det(A - \lambda I) = 0。它的根就是特征值。然后把每个 \lambda 代回 (A - \lambda I)\mathbf{x} = \mathbf{0} 求出对应的特征向量。

• 关键性质：

  • A 的迹等于其特征值之和。
  • A 的行列式等于其特征值之积。
  • 对称矩阵的特征向量两两垂直，特征值全是实数。
  • 正定矩阵的特征值全为正。
  • 协方差矩阵（我们在统计学中会遇到）总是半正定的。

• 对于大矩阵，用特征多项式算特征值是不现实的，我们改用迭代方法：

  • 幂迭代（power iteration）：反复乘以 A 并归一化。收敛到主特征向量（最大特征值）。简单，但只能找到一个特征对。

  • QR 算法（QR algorithm）：主力方法。反复用 QR 分解来分解再重组，直到矩阵收敛为三角形式，所有特征值就出现在对角线上。

  • 逆迭代（inverse iteration）：找到最接近某个目标值的特征向量。当你大致知道想要哪个特征值时很有用。

  • 对大型稀疏矩阵，Arnoldi 和 Lanczos 迭代利用稀疏性来提高效率。

• 如果方阵有一组完整的线性无关特征向量，它就可以被对角化（diagonalised）：A = PDP^{-1}，其中 D 是由特征值组成的对角矩阵，P 的各列是特征向量。

• 这有什么用？对角矩阵处理起来轻而易举。需要 A^{100}？与其把 A 自乘 100 次，不如计算 PD^{100}P^{-1}——把对角矩阵自乘只需独立地把每个元素自乘。这把一个昂贵的运算变成了廉价的。

• **特征基（eigenbasis）**是完全由特征向量构成的基。在这个基下，矩阵变成对角阵，变换就是沿每个特征向量方向的独立缩放。这就像为这个变换找到了最自然的坐标系。

• **QR 分解（QR decomposition）**把任意矩阵 A 分解为 A = QR，其中 Q 是正交矩阵（各列标准正交），R 是上三角。可以把它看成把「方向」信息（Q）和「缩放与混合」信息（R）分离开来。

• Gram-Schmidt 过程逐列构造 Q。取 A 的第一列并归一化。取第二列，减去它在第一列上的投影（让它与第一列垂直），再归一化。对每一列重复。结果就是一组标准正交的向量。

• QR 分解是求特征值的 QR 算法背后的引擎。它也直接用于求解最小二乘问题：当 A\mathbf{x} = \mathbf{b} 没有精确解（方程数多于未知数）时，QR 给出最佳的近似解。

• SVD（奇异值分解，Singular Value Decomposition）是最通用、也可以说是最重要的分解。每个矩阵（任意形状、任意秩）都有 SVD：A = U\Sigma V^T

  • V^T（n \times n，正交）：旋转输入
  • \Sigma（m \times n，对角）：沿正交轴缩放（奇异值，非负，按降序排列）
  • U（m \times m，正交）：旋转输出

• 从几何上看，SVD 告诉我们：每个线性变换无论多复杂，都只是一次旋转、沿各轴的拉伸、再一次旋转。一个圆变成一个椭圆。

• 奇异值（\sigma_1 \geq \sigma_2 \geq \ldots）揭示了每个方向的「重要性」。大的奇异值对应最重要的方向。A 的秩等于非零奇异值的个数。

• 低秩近似（low-rank approximation）：只保留最大的 k 个奇异值、把其余置零，就得到 A 最好的秩-k 近似。图像压缩就是这样工作的：一张 1000 \times 1000 的图像可能只需 k = 50 个奇异值就看起来几乎一样，把它压缩了 20 倍。

• SVD 还给出伪逆：A^+ = V\Sigma^+U^T，其中 \Sigma^+ 把非零奇异值取倒数。

• 特征值分解只适用于方阵，而 SVD 适用于任意矩阵。这是它的关键优势。

• PCA（主成分分析，Principal Component Analysis）用特征值分解（或 SVD）来降维。

• 想象一个数据集，每个样本有 100 个特征（一个 100 维向量叠成一个矩阵）。这些特征中很多是相关且冗余的。

• PCA 找出数据真正变化的方向，让你只保留重要的部分。

• 第一主成分（PC1）是方差最大的方向。

• 第二主成分（PC2）捕捉剩余部分中最大的方差，并且与第一个垂直。

• 如果大部分方差都集中在少数几个方向上，你就可以把数据投影到这些维度，丢弃其余部分，而几乎不损失信息。

• 步骤：

  • 标准化数据（减去均值、除以标准差），让所有特征的贡献均等
  • 计算协方差矩阵
  • 求它的特征值和特征向量
  • 选出特征值最大的 k 个特征向量（它们就是主成分）
  • 把数据投影到这些成分上

• 标准化至关重要：没有它，以千米为单位的特征会压倒以厘米为单位的特征，不管实际重要性如何。

• 在实际中，PCA 用于可视化（把高维数据投影到二维或三维）、降噪（丢弃主要是噪声的低方差方向），以及通过减少输入特征数来加速 ML 模型。

• **核 PCA（Kernel PCA）**把 PCA 推广到非线性关系。它通过一个核函数把数据映射到更高维的空间，在那里结构变成线性的，然后应用标准 PCA 并投影回来。

• **Schur 分解（Schur decomposition）**把方阵分解为 A = QTQ^\ast，其中 Q 是酉矩阵，T 是上三角。每个方阵都有 Schur 分解，即使它不能被对角化。

• **非负矩阵分解（Non-negative Matrix Factorisation，NMF）**把矩阵分解成两个非负矩阵：A \approx WH，其中 W 和 H 的所有元素都 \geq 0。与可能产生负值的 SVD 不同，NMF 只做加法、不做减法。这让各部分变得可解释：在主题建模中，W 给出每篇文档的主题权重，H 给出每个主题的词权重，全部非负——正好契合我们思考「一篇文档包含多少每个主题」的方式。

• **谱定理（spectral theorem）**表明，对称（或 Hermitian）矩阵总能用正交（或酉）矩阵对角化。它们的特征值总是实数，特征向量总是正交。这是 PCA 背后的理论基础。

编程练习（使用 CoLab 或 notebook）

1. 计算一个对称矩阵的特征值和特征向量。验证特征向量两两垂直，并用特征值分解重建这个矩阵。

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: {eigenvalues}")
print(f"Eigenvectors orthogonal: {jnp.dot(eigenvectors[:,0], eigenvectors[:,1]):.6f}")

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

2. 实现幂迭代来求最大特征值，用逆迭代求最小特征值，并与 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"Largest eigenvalue:  {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"Smallest eigenvalue: {1.0 / (v @ jnp.linalg.solve(A, v)):.4f}")

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

3. 计算一个矩阵的 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}, reconstruction error: {error:.4f}")
​