线性方程组:最古老且仍在跑神经网络的数学

线性方程组:最古老且仍在跑神经网络的数学

本节摘要:求解 Ax = b 是数学里最古老、却仍在跑你神经网络的问题。每次训练线性回归、每次最小二乘拟合、每个 y = Wx + b 层、每个高斯过程的矩阵分解、每次为马氏距离求协方差逆,都是在解线性系统。本节从几何直觉讲起(每个方程定义一张超平面,解是所有超平面的交点;行图景是方程组、列图景是「A 的列的何种线性组合得到 b」,后者更根本——b 在列空间则有解,否则找最近点即最小二乘);给出高斯消元(带部分主元)与回代;推导三大分解——LU(因子化摊销 O(n³) 成本,后续每个 b 只 O(n²))、QR(正交三角化,最小二乘更稳)、Cholesky(对称正定,比 LU 快两倍、省一半存储);连上最小二乘的正规方程 AᵀAx = Aᵀb(就是线性回归闭式解,加 λI 即岭回归);用条件数 κ = σ_max/σ_min 诊断病态(损失约 log₁₀(κ) 位精度);最后给共轭梯度等迭代法与「何时用何法」总表。

对应原课程:Phase 01 · Lesson 17 · linear-systems(原英文 phases/01-math-foundations/17-linear-systems/docs/en.md)。前置:第 1、2、3 节(线性代数)。

学习目标

阅读完本节,你应当能够:

1. 用带部分主元的高斯消元与回代求解 Ax = b。
2. 用 LU、QR、Cholesky 分解因子化矩阵,并说明各自适用场景。
3. 推导最小二乘的正规方程,并连上线性回归与岭回归。
4. 用条件数诊断病态系统,用正则化稳定之。

一、问题与直觉

Ax = b 无处不在:A 是已知系数矩阵,b 是已知输出向量,x 是你想求的未知向量。线性回归里 A 是数据矩阵、b 是目标向量、x 是权重向量,整个模型归结为「找 x 使 Ax 尽量接近 b」。本节从零搭建解这个方程的所有主要方法,让你明白为何有的快、有的稳、有的只适合方阵、有的能处理超定,以及为何矩阵的条件数决定你的答案是否还有意义。

1.1 Ax = b 的几何含义

线性方程组的几何解读:每个方程定义一张超平面,解是所有超平面的交点。

2x + y = 5          2D 里两条直线
x − y  = 1          交于 x=2, y=1
​

三种情况:① 唯一解——A 可逆;② 无解——系统不相容;③ 无穷多解——A 有零空间(列线性相关)。大多数 ML 问题落在「无精确解」类,因为方程(数据点)比未知量(参数)多——这就是最小二乘登场之处。

1.2 列图景 vs 行图景

读 Ax = b 有两种方式:

• 行图景:A 的每一行定义一个方程(一张超平面),解是它们全体的交。
• 列图景:A 的每一列是一个向量,问题变成「A 各列的何种线性组合产生 b?」

A = |2  1|,b = |5|
    |1 -1|    |1|
列图景:找 x1,x2 使 x1·[2,1] + x2·[1,-1] = [5,1]
        2·[2,1] + 1·[1,-1] = [5,1] ✓
​

列图景更根本:b 在 A 的列空间则有解;b 不在则找列空间中最近点,那就是最小二乘解。

1.3 高斯消元

高斯消元把 Ax = b 变成上三角 Ux = c,再回代求解。

1. 对每列 k(主元列):
   a. 找列 k 中第 k 行及以下的最大元(部分主元)
   b. 把那行换到第 k 行
   c. 对 k 以下每行 i:算乘子 m = A[i][k]/A[k][k];第 i 行减 m 倍第 k 行
2. 回代:从最后一个方程向上解
​

成本 O(n³),1000×1000 系统约 10 亿次浮点运算。快,但若要解同一 A 的多个 b,还有更好办法。

1.4 部分主元:为何要紧

不加主元,高斯消元会失败或产出垃圾:主元为零则除零,主元太小则放大舍入误差。部分主元总选最大可用主元以最小化误差放大。在有限精度浮点下,无主元版会丢有效位——选最大主元是数值安全的根本。

1.5 LU 分解

LU 把 A 因子化成下三角 L 与上三角 U:A = LU。L 存高斯消元的乘子,U 是消元结果。为何要因子化而非直接消元?因为一旦有了 L 与 U,对任何新 b 解 Ax = b 只需 O(n²):LUx = b,令 y=Ux,先 Ly = b(前代,O(n²)),再 Ux = y(回代,O(n²))。O(n³) 成本在因子化时付一次,之后每次求解 O(n²)。若同一 A 解 1000 个不同 b,LU 省功 1000/3 倍。带部分主元得 PA = LU,P 是记录换行的排列矩阵。

1.6 QR 分解

QR 把 A 因子化成正交矩阵 Q 与上三角 R:A = QR。正交矩阵 QᵀQ = I,列是标准正交向量,乘 Q 保持长度与角度。解 Ax = b:QRx = b → Rx = Qᵀb(乘 Qᵀ,无需求逆)→ 回代。QR 解最小二乘比 LU 数值更稳。Gram-Schmidt 过程逐列建 Q:每步去掉沿所有已有 q 的分量,只留新的正交方向。

1.7 Cholesky 分解

当 A 对称(A = Aᵀ)且正定(所有特征值为正),可因子化 A = LLᵀ,L 是下三角。Cholesky 比 LU 快两倍、省一半存储,只对对称正定矩阵有效,但这类矩阵常出现:协方差矩阵(对称半正定,加正则即正定);高斯过程的核矩阵(对称正定);凸函数在极小值的 Hessian(对称正定);AᵀA 恒对称半正定。高斯过程里用 Cholesky 分解核矩阵 K,再解 Kα = y 得预测均值;Cholesky 因子还给边缘似然的对数行列式 log det(K) = 2·Σ log(diag(L))。

1.8 最小二乘:Ax = b 无精确解时

A 是 m×n 且 m>n(方程比未知量多)时系统超定,无精确解,改为最小化平方误差:min ‖Ax − b‖²。极小点满足正规方程 AᵀAx = Aᵀb。推导:展开 ‖Ax-b‖² = (Ax-b)ᵀ(Ax-b) = xᵀAᵀAx - 2xᵀAᵀb + bᵀb,对 x 求梯度置零得 2AᵀAx - 2Aᵀb = 0。

1.9 正规方程 = 线性回归

联系是精确的。线性回归里数据矩阵 X 每行一样本、每列一特征,目标向量 y 每样本一项,权重 w 满足 XᵀXw = Xᵀy,w = (XᵀX)⁻¹Xᵀy——这是线性回归闭式解,每次 sklearn.linear_model.LinearRegression.fit() 都在算它(或等价的 QR/SVD)。给矩阵加 λI 正则项即岭回归:(XᵀX + λI)w = Xᵀy。正则使矩阵更好条件(更易精确求逆)、把权重向零收缩防过拟合;XᵀX + λI 在 λ>0 时恒对称正定,故可用 Cholesky 解。

1.10 伪逆(Moore-Penrose)

伪逆 A⁺ 把矩阵求逆推广到非方阵与奇异矩阵。对任意 A:x = A⁺b,其中 A⁺ = VΣ⁺Uᵀ(经 SVD),Σ⁺ 是把每个非零奇异值取倒数再转置。伪逆给最小范数最小二乘解:有唯一解则 A⁺b 给之;无解则给最小二乘解;无穷多解则给 ‖x‖ 最小者。NumPy 的 np.linalg.lstsq 与 np.linalg.pinv 内部都用 SVD。

1.11 条件数

条件数量化解对小输入扰动的敏感度:κ(A) = ‖A‖·‖A⁻¹‖ = σ_max/σ_min。

良态(κ~1):b 小变 → x 小变        病态(κ~10¹⁵):b 小变 → x 巨变
|2 0|,κ=2/1=2,安全解              |1   1         |,κ~10¹⁵,解是垃圾
|0 1|                              |1   1+10⁻¹⁵   |
​

经验:κ<100 安全;κ~10^k 损失约 k 位精度;κ~10¹⁶(float64)解无意义,矩阵事实上奇异。ML 里特征近共线时病态,正则化(加 λI)把条件数从 σ_max/σ_min 改善到 (σ_max+λ)/(σ_min+λ)。

1.12 迭代法:共轭梯度

极大稀疏系统(百万未知量)直接法(LU、Cholesky)太贵。迭代法通过多轮改进猜测逼近解。共轭梯度(CG) 解 A 对称正定的 Ax = b,精确算术下最多 n 步得精确解,但 A 特征值聚集时通常快得多。用于:大规模优化(Newton-CG)、解 PDE 离散化、核矩阵太大无法因子化、作其他迭代解器的预条件。收敛率取决于条件数,条件越好收敛越快——这又是正则化有用的理由。

1.13 全景:何时用何法

1.14 与 ML 的联系

每个方法都在生产 ML 里出现:

• 线性回归:闭式解解正规方程 XᵀXw = Xᵀy,用 Cholesky(n 小)、QR(数值要紧)或 SVD(可能秩亏)。
• 岭回归:加 λI,正则系统 (XᵀX+λI)w = Xᵀy 在 λ>0 时恒对称正定,可用 Cholesky。
• 高斯过程:预测均值需解 Kα = y,Cholesky 分解 K 是标准做法;对数边缘似然用 log det(K) = 2·Σ log(diag(L))。
• 神经网络初始化:正交初始化用 QR 分解造列标准正交的权重矩阵,防深网信号塌缩。
• 预条件:大规模优化器用不完全 Cholesky 或不完全 LU 作共轭梯度的预条件子。
• 特征工程:XᵀX 的条件数告诉你特征是否共线,κ 大就丢特征或加正则。

二、从零实现

完整源码见 phases/01-math-foundations/17-linear-systems/code/linear_systems.py。

2.1 带部分主元的高斯消元

def gaussian_elimination(A, b):
    n = len(b)
    Ab = np.hstack([A.astype(float), b.reshape(-1, 1).astype(float)])
    for k in range(n):
        max_row = k + np.argmax(np.abs(Ab[k:, k]))      # 部分主元
        Ab[[k, max_row]] = Ab[[max_row, k]]
        if abs(Ab[k, k]) < 1e-12:
            raise ValueError(f"主元 {k} 处矩阵奇异或近奇异")
        for i in range(k + 1, n):
            m = Ab[i, k] / Ab[k, k]
            Ab[i, k:] -= m * Ab[k, k:]                   # 消元
    x = np.zeros(n)
    for i in range(n - 1, -1, -1):                       # 回代
        x[i] = (Ab[i, -1] - Ab[i, i+1:n] @ x[i+1:n]) / Ab[i, i]
    return x
​

2.2 LU 分解

def lu_decompose(A):
    n = A.shape[0]
    L = np.eye(n); U = A.astype(float).copy(); P = np.eye(n)
    for k in range(n):
        max_row = k + np.argmax(np.abs(U[k:, k]))
        if max_row != k:
            U[[k, max_row]] = U[[max_row, k]]; P[[k, max_row]] = P[[max_row, k]]
            if k > 0: L[[k, max_row], :k] = L[[max_row, k], :k]
        for i in range(k + 1, n):
            L[i, k] = U[i, k] / U[k, k]
            U[i, k:] -= L[i, k] * U[k, k:]
    return P, L, U
​

2.3 Cholesky 分解

def cholesky(A):
    n = A.shape[0]; L = np.zeros_like(A, dtype=float)
    for i in range(n):
        for j in range(i + 1):
            s = A[i, j] - L[i, :j] @ L[j, :j]
            if i == j:
                if s <= 0: raise ValueError("矩阵非正定")
                L[i, j] = np.sqrt(s)
            else:
                L[i, j] = s / L[j, j]
    return L
​

2.4 正规方程最小二乘与岭回归

def least_squares_normal(A, b):
    return gaussian_elimination(A.T @ A, A.T @ b)        # 注意:平方条件数,见第 11、13 节

def ridge_regression(A, b, lam):
    n = A.shape[1]
    L = cholesky(A.T @ A + lam * np.eye(n))              # 对称正定,用 Cholesky
    y = forward_substitute(L, A.T @ b)
    return back_substitute(L.T, y)
​

2.5 条件数

def condition_number(A):
    U, S, Vt = np.linalg.svd(A)
    return S[0] / S[-1]
​

三、框架对比

把零件拼起来,在真实数据上跑线性回归与岭回归:

np.random.seed(42)
X_raw = np.random.randn(100, 3)
w_true = np.array([2.0, -1.0, 0.5])
y = X_raw @ w_true + np.random.randn(100) * 0.1
X = np.column_stack([np.ones(100), X_raw])               # 加截距列

w_ols = least_squares_normal(X, y)
w_np = np.linalg.lstsq(X, y, rcond=None)[0]              # NumPy(内部 SVD,更稳)
print(f"最大差异: {np.max(np.abs(w_ols - w_np)):.2e}")

w_ridge = ridge_regression(X, y, lam=1.0)
from sklearn.linear_model import Ridge
ridge_sk = Ridge(alpha=1.0, fit_intercept=False).fit(X, y)
​

💡 正式场合优先用 np.linalg.lstsq 或 np.linalg.solve 而非手写正规方程——正规方程 AᵀA 平方条件数(见第 11、13 节),高条件数矩阵会丢精度;lstsq 内部走 SVD,数值更稳。

四、可复用产物

• code/linear_systems.py:高斯消元、LU、Cholesky、最小二乘、岭回归的从零实现。
• 一个验证 demo:正规方程与 sklearn 的 LinearRegression 给出相同权重。

源码见 phases/01-math-foundations/17-linear-systems/code/。

五、练习

1. (Easy) 用你的高斯消元、LU 解器、np.linalg.solve 三种方法解 [[1,2,3],[4,5,6],[7,8,10]] x = [6,15,27],验证三者给同一答案(浮点容差内)。
2. (Medium) 生成 50×5 随机矩阵 X 与目标 y,用正规方程、QR(np.linalg.qr)、SVD(np.linalg.svd)、np.linalg.lstsq 四种方法求 w,比较;测 XᵀX 条件数,解释它如何影响你信哪种方法。
3. (Medium) 造近奇异矩阵(第二列 ≈ 第一列 + 1e-10·噪声),算其条件数,有/无正则(加 0.01·I)分别解 Ax = b,比较解与残差,解释正则为何有用。
4. (Hard) 实现 100×100 随机对称正定矩阵的共轭梯度,数收敛到 1e-8 容差需几步,与理论上限 n 步比较。
5. (Hard) 对大小 10、50、200、500 的对称正定矩阵,计时你的 Cholesky vs LU vs np.linalg.solve,画图,验证 Cholesky 约比 LU 快 2 倍。

本节要点回顾

1. Ax = b 是 ML 的底层运算——线性回归、最小二乘、Wx+b 层、高斯过程、马氏距离逆协方差,全是它在跑。
2. 列图景更根本:b 在列空间则有解,不在则找最近点(最小二乘)。
3. 高斯消元 + 部分主元:O(n³),总选最大主元最小化误差放大。
4. LU 摊销成本:因子化 O(n³) 付一次,之后每个 b 只 O(n²)(前代 + 回代)。
5. QR 解最小二乘更稳:正交三角化,Gram-Schmidt 逐列去冗余建 Q。
6. Cholesky 对称正定专用:比 LU 快两倍、省一半存储,协方差、核矩阵、AᵀA、岭回归的 (XᵀX+λI) 都用它。
7. 正规方程 AᵀAx = Aᵀb = 线性回归闭式解,加 λI 即岭回归。
8. 伪逆 A⁺ = VΣ⁺Uᵀ(经 SVD)给最小范数最小二乘解,np.linalg.lstsq/pinv 内部都用 SVD。
9. 条件数 κ = σ_max/σ_min 量敏感度,损失约 log₁₀(κ) 位精度;κ~10¹⁶ 解无意义。
10. 共轭梯度 解大稀疏对称正定系统,最多 n 步,条件越好收敛越快。

下一节,我们看优化何时有保证——凸优化:凸问题只有一个谷,任何下山算法都到全局最优;牛顿法用 Hessian 二阶信息;Lagrange 乘子与 KKT 把约束变无约束;以及为何非凸的神经网络 SGD 仍能找到好解。