SVD 与低秩压缩

它到底有多省？原本存 \(A\) 要 \(m\times n\) 个数；存 \(A_k\) 只需存前 \(k\) 个 \(u_i\)（共 \(mk\) 个数）、\(k\) 个 \(\sigma_i\)、\(k\) 个 \(v_i\)（共 \(nk\) 个数），合计约 \(k(m+n)\)。当 \(k\) 远小于 \(m,n\) 时，省下的存储是数量级的。

Eckart–Young 定理：截断是"最优"的

你可能担心：随手砍掉尾巴，会不会有更聪明的低秩矩阵比 \(A_k\) 更接近 \(A\)？Eckart–Young 定理给了一个让人安心的回答——没有。在所有秩不超过 \(k\) 的矩阵里，\(A_k\) 是最优逼近，而且误差有简洁的闭式：

这里 \(\|\cdot\|_F\) 是 Frobenius 范数（把矩阵所有元素平方求和再开根，相当于把矩阵拉直成长向量后的欧氏长度）。这条定理是低秩压缩的理论基石：它保证"按奇异值大小截断"不是某种近似启发式，而是数学上可证明的最优策略。误差就等于"被你扔掉的那些奇异值的能量"（开根号）——正是上面定义的那个能量。

import numpy as np
np.set_printoptions(precision=4, suppress=True)

A = np.array([[3., 2., 2.],
              [2., 3., -2.]])
U, s, Vt = np.linalg.svd(A, full_matrices=False)

# 秩-1 近似 = 第一项外积
A1 = s[0] * np.outer(U[:, 0], Vt[0])
print("A1 (秩-1 近似) =\n", A1)

err = np.linalg.norm(A - A1)              # Frobenius 范数
print("秩-1 逼近误差 =", round(err, 4))     # 3.0
print("Eckart-Young 预言 = sqrt(σ2^2) =", round(s[1], 4))  # 3.0 —— 完全一致

7. 条件数：矩阵"健不健康"的体检指标

奇异值里最大和最小的那两个，组成一个判断矩阵数值健康度的关键量——条件数 condition number：

几何上它度量"椭圆有多扁"：\(\kappa=1\) 是完美的圆（各方向增益相同，最健康）；\(\kappa\) 很大说明椭圆被压得又细又长，某些方向几乎被压没了。这样的矩阵叫病态 ill-conditioned。

为什么要在意？因为解方程 \(Ax=b\) 或做矩阵求逆时，条件数刻画了相对误差的最坏放大上界：严格地说，\(\dfrac{\|\Delta x\|}{\|x\|}\le \kappa\cdot\dfrac{\|\Delta b\|}{\|b\|}\)。也就是说，若 \(b\) 有千分之一的相对扰动，解 \(x\) 的相对误差最坏可达 \(\kappa\times\) 千分之一。注意这是"相对对相对"的最坏情形上界，不是"任意绝对扰动一律被放大 \(\kappa\) 倍"。\(\kappa\) 上万的矩阵，计算结果的有效数字会被吃掉好几位。

import numpy as np
B = np.array([[2., 0.], [0., 2.]])
C = np.array([[1., 1.], [1., 1.0001]])
print("健康矩阵 B 的条件数:", np.linalg.cond(B))   # 1.0，完美的圆
print("病态矩阵 C 的条件数:", round(np.linalg.cond(C), 1))  # ~40002
print("C 的奇异值:", np.linalg.svd(C, compute_uv=False))

8. 伪逆：当矩阵不能求逆时

普通的逆 \(A^{-1}\) 只对"方且满秩"的矩阵存在。可现实里我们经常想"反解"一个长方形或病态的方程组——比如最小二乘拟合。SVD 给出了一个万能替代品：伪逆（Moore–Penrose pseudo-inverse），记作 \(A^+\)。

做法朴素得可爱：既然 \(A=U\Sigma V^\top\)，那就把三个因子各自"反着来"，并把 \(\Sigma\) 里每个非零奇异值取倒数（零保持为零）：

这里要补一句形状（第 3 节那种维度对账别丢）：\(\Sigma\) 是 \(m\times n\)，而 \(\Sigma^+\) 的形状是 \(n\times m\)——是把 \(\Sigma\) 先转置成 \(n\times m\)，再对非零对角元取倒数。这样 \(V(n\times n)\cdot\Sigma^+(n\times m)\cdot U^\top(m\times m)\) 才正好拼成 \(n\times m\) 的 \(A^+\)（恰好是 \(A\) 形状的"转置"，符合"反解"的直觉）。

当 \(A\) 可逆时，\(A^+\) 恰好就是 \(A^{-1}\)；当 \(A\) 是长方形或不满秩时，\(A^+\) 给出最小二乘意义下的"最佳反解"。这也解释了第 7 节的病态问题：伪逆本身在数学上是良定义的，麻烦出在极小的奇异值取倒数后变成一个巨大的数——当 \(b\) 含噪声时，它会把 \(b\) 在最小奇异方向上的噪声分量放大 \(1/\sigma_{\min}\) 倍，于是解被噪声主导、极不稳定（这不是浮点溢出，而是"噪声放大"）。所以实践中常把太小的 \(\sigma_i\) 直接当 0 截掉（这叫"截断伪逆"，本质又回到了低秩思想）。

import numpy as np
np.set_printoptions(precision=4, suppress=True)
A = np.array([[3., 2., 2.],
              [2., 3., -2.]])     # 2x3，无法求普通逆
Aplus = np.linalg.pinv(A)        # 伪逆，形状 3x2
print("A+ 形状:", Aplus.shape)
print("A @ A+ =\n", A @ Aplus)   # 应得到 2x2 单位阵 I

9. 和 ML 的联系：低秩思想无处不在

import numpy as np
# LoRA 省了多少参数？
d, r = 4096, 8
full = d * d                 # 完整 ΔW 的参数量
lora = 2 * d * r             # B(d×r) + A(r×d)
print("完整微调参数:", full)            # 16,777,216
print("LoRA 参数  :", lora)            # 65,536
print("占比:", round(lora / full * 100, 3), "%")   # 0.391 %

10. 综合代码：把小矩阵当"图像"压缩

下面这段把前面所有概念串起来：合成一个 \(8\times8\) 的"低秩+噪声"小图，看奇异值如何衰减、累计能量如何快速逼近 1、以及不同 \(k\) 的重建误差是否精确符合 Eckart–Young。

import numpy as np
np.set_printoptions(precision=3, suppress=True)
rng = np.random.default_rng(0)

# 合成一张 8x8 的"图像"：一个秩-1 主结构 + 少量噪声
u1 = np.array([1, 1, 1, 1, 2, 2, 2, 2.])
v1 = np.array([1, 1, 1, 1, 1, 1, 2, 2.])
M = np.outer(u1, v1) + 0.3 * rng.standard_normal((8, 8))

U, s, Vt = np.linalg.svd(M, full_matrices=False)
print("奇异值:", s)                       # 第一个远大于其余 —— 主结构

# 累计能量：前 k 个奇异值占了多少"信息"（能量 = σ_i^2）
energy = np.cumsum(s**2) / np.sum(s**2)
print("累计能量比:", energy.round(4))      # 秩-1 就已 ~99%

# 各 k 的秩-k 重建误差 vs Eckart-Young 闭式预言
for k in [1, 2, 3]:
    Mk = (U[:, :k] * s[:k]) @ Vt[:k]      # = Σ_{i


要点：本课一图流

分解：任意 \(A=U\Sigma V^\top\)，几何上 = 旋转 \(V^\top\) → 轴向缩放 \(\Sigma\) → 旋转 \(U\)（"旋转"泛指正交变换，可能含反射）。
构件：\(\sigma_i\) 是各方向增益；\(U\)（\(m\times m\)）、\(V\)（\(n\times n\)）的列是左/右奇异向量；\(\sigma_i^2\) 是 \(A^\top A\) / \(AA^\top\) 的特征值。
压缩：\(A=\sum\sigma_i u_i v_i^\top\)，截断到前 \(k\) 项即最优低秩逼近（Eckart–Young），误差 \(=\sqrt{\sum_{i>k}\sigma_i^2}\)；这里"能量"指 \(\sigma_i^2\)，总能量 \(=\|A\|_F^2\)。
健康：条件数 \(\kappa=\sigma_{\max}/\sigma_{\min}\)（相对误差的最坏放大上界）；伪逆 \(A^+=V\Sigma^+U^\top\)（\(\Sigma^+\) 为 \(n\times m\)）是万能"反解"。
ML：图像压缩 / PCA / 推荐 / 白化都是它，LoRA 是"事前低秩约束"版的同一思想。



调一调，观察现象
下面三个小任务都基于本课已经验证过的代码，改一个数、跑一下，看奇异值/能量/条件数怎么随之变化——把"低秩近似"和"病态"从文字变成你亲眼跑出来的数字。


任务 1：改截断秩 k，看重建误差
改什么：对本课的 \(A\)（奇异值 5 和 3）把 k 从 1 调到 2。

预期现象：\(k=1\) 时重建误差 \(=3.0000\)，恰好等于被丢弃的 \(\sigma_2\)；\(k=2\) 时误差降到 \(0.0000\)（两个非零奇异值全留下，精确还原）。

为什么：Eckart–Young 说截断误差 \(=\sqrt{\sum_{i>k}\sigma_i^2}\)。\(k=1\) 丢掉 \(\sigma_2=3\)，误差就是 3；留满秩则无可丢，误差为 0。
import numpy as np
A = np.array([[3.,2.,2.],[2.,3.,-2.]])
U, s, Vt = np.linalg.svd(A, full_matrices=False)
for k in [1, 2]:
    Ak = (U[:, :k] * s[:k]) @ Vt[:k]
    err = np.linalg.norm(A - Ak)
    pred = np.sqrt((s[k:]**2).sum())
    print(f"k={k}: 实测误差={err:.4f}, 理论={pred:.4f}")




任务 2：让两行越来越像，看条件数飙升
改什么：把第 7 节病态矩阵里第二行的小扰动 eps 从 1.0 一路调到 0.001。

预期现象：条件数依次约为 6.85 → 42 → 4002，越调越大；同时最小奇异值从 0.382 一路缩到 0.0005，几乎贴到 0。

为什么：eps 越小两行越接近线性相关，矩阵越逼近秩 1，于是 \(\sigma_{\min}\to 0\)，而 \(\kappa=\sigma_{\max}/\sigma_{\min}\) 随之爆炸——这正是"病态"的来源。
import numpy as np
for eps in [1.0, 0.1, 0.001]:
    C = np.array([[1., 1.], [1., 1. + eps]])
    s = np.linalg.svd(C, compute_uv=False)
    print(f"eps={eps}: 条件数={np.linalg.cond(C):.2f}, 奇异值={s.round(4)}")




任务 3：加大噪声，看能量从"集中"变"分散"
改什么：对一个秩-1 主结构 + 噪声的 \(8\times8\) 矩阵，把噪声强度 noise 从 0.05 调到 1.5。

预期现象：秩-1 占的累计能量比从约 0.9997 一路掉到约 0.76——噪声越大，能量越分散到其余奇异值上。

为什么："能量"就是 \(\sigma_i^2\)。低噪时主结构几乎独占第一个奇异值，秩-1 就够；噪声把能量摊到所有方向，第一个方向的占比自然下降——这也是"自然数据才低秩、纯噪声不低秩"的直观体现。
import numpy as np
rng = np.random.default_rng(0)
u1 = np.array([1,1,1,1,2,2,2,2.]); v1 = np.array([1,1,1,1,1,1,2,2.])
for noise in [0.05, 0.3, 1.5]:
    M = np.outer(u1, v1) + noise * rng.standard_normal((8, 8))
    s = np.linalg.svd(M, compute_uv=False)
    e1 = (s[0]**2) / (s**2).sum()
    print(f"noise={noise}: 秩-1 累计能量比={e1:.4f}")

动手练习


手算奇异值。对 \(A=\begin{bmatrix}2&0\\0&-3\end{bmatrix}\)，先用 \(A^\top A\) 手算两个奇异值（提示：对角阵很好算，注意奇异值非负），再用代码验证。
import numpy as np
A = np.array([[2., 0.], [0., -3.]])
s = np.linalg.svd(A, compute_uv=False)
print("奇异值:", s)   # 和你手算的对一下？(应为 3, 2)
print("特征值:", np.linalg.eigvals(A))  # 对比：特征值可以是负的

验证 Av=σu。对任意你喜欢的 \(3\times2\) 矩阵做 SVD，取出 \(v_1, u_1, \sigma_1\)，验证 \(Av_1\) 是否等于 \(\sigma_1 u_1\)。
import numpy as np
A = np.array([[1., 2.], [3., 4.], [5., 6.]])
U, s, Vt = np.linalg.svd(A, full_matrices=False)
v1 = Vt[0]          # 第一个右奇异向量
# TODO: 计算 A @ v1，并与 s[0]*U[:,0] 比较，打印两者的差的范数

能量曲线。合成一个 \(20\times20\) 的随机矩阵和一个"秩-3 + 噪声"矩阵，分别打印它们的累计能量比。对比：随机矩阵的能量是不是衰减得慢得多？（体会"自然数据才低秩"。）

截断伪逆。给第 7 节的病态矩阵 \(C\) 做 SVD，手动把最小的奇异值置 0 后再按 \(V\Sigma^+U^\top\) 拼伪逆，比较它与 np.linalg.pinv(C) 在解一个带噪声 \(b\) 时谁更稳。
import numpy as np
C = np.array([[1., 1.], [1., 1.0001]])
U, s, Vt = np.linalg.svd(C)
# TODO: 复制一份 s，把 s[-1] 设为很小或 0，重建截断伪逆并打印对比

（挑战）小图压缩比。合成一个 \(32\times32\)、秩约为 5 的"图像"，找出"累计能量达到 99%"所需的最小 \(k\)，并算出相对原始存储 \(32\times32\) 的压缩比 \(k(32+32)/(32\times32)\)。


掌握自检


给定一个 \(m\times n\) 矩阵，我能立刻说出 \(U,\Sigma,V\) 各自的形状，并解释 \(\Sigma\) 为什么是长方形的。
我能用"圆 → 旋转 → 椭圆 → 旋转"向别人讲清 SVD 的几何含义，并指出椭圆半轴就是奇异值。
我能说出 \(V\)、\(U\)、\(\sigma_i\) 分别和 \(A^\top A\)、\(AA^\top\) 的特征向量/特征值的对应关系，并解释为什么 SVD 对任意矩阵都存在。
我能写出秩-k 截断公式，说出 Eckart–Young 定理的内容，并知道逼近误差等于被丢弃奇异值的能量 \(\sqrt{\sum_{i>k}\sigma_i^2}\)（能量 = \(\sigma_i^2\)）。
我能解释条件数 \(\kappa=\sigma_{\max}/\sigma_{\min}\) 的几何与数值含义（相对误差的最坏放大上界），看到 \(\kappa\) 很大就警觉"病态"。
我能用一句话讲清 LoRA 为什么省参数，并指出它和"秩-k 截断"的关系（事前约束 vs 事后压缩）。


到这里，模块 1 的线性代数部分就全部讲完了——从向量、矩阵变换、内积投影、特征分解，一路走到今天的 SVD，你已经握齐了理解神经网络"前向计算"所需的全部代数工具。下一课我们换挡进入微积分：导数与梯度，去看模型是怎么学的——也就是它如何顺着"下坡方向"一步步调整自己的参数向量。
可以先放过的点
第 4 节里"奇异值为 0 的列怎么补齐"、四个基本子空间（零空间、左零空间……）的细节，现在只要记住"SVD 自动把方向分成通得过和通不过两堆"就够，等到后面专门的线代专题再回来抠。
Eckart–Young 定理的证明不必现在啃下来——你先把它当成"按奇异值大小截断就是最优"这条可放心使用的结论，会用就行。
第 8 节伪逆里 \(\Sigma^+\) 是 \(n\times m\)、要先转置再取倒数这类维度对账，第一遍读糊了没关系，真正用到最小二乘求解时（以后的回归课）再回头逐维数清楚。
条件数那条不等式 \(\frac{\|\Delta x\|}{\|x\|}\le\kappa\frac{\|\Delta b\|}{\|b\|}\) 的严格推导可以暂缓，现在记住"\(\kappa\) 很大 = 病态 = 要警觉"这个直觉即可。