inv、det、rank、eig 线性代数四件套：求逆、行列式、秩、特征值各一行代码

一、主线确立：为什么"一行代码"是最危险的抽象

线性代数库的设计哲学是"把复杂留给自己，把简洁留给用户"。NumPy 的 linalg 模块把 LAPACK 的数千行 Fortran 代码封装成一行调用，这既是工程胜利，也是认知陷阱。当一位工程师写下 A_inv = np.linalg.inv(A) 时，他通常并不知道：这个调用内部走的是 LU 分解加三角回代，浮点运算量约为 2n³/3 次，而误差放大倍数由矩阵条件数 κ(A) 决定。

本文评述：这四个函数之所以值得单独成文，不是因为它们难，而是因为它们太容易被误用。误用的根源在于，API 的简洁性掩盖了三个关键信息——算法复杂度、数值稳定性、以及"是否存在更好的替代路径"。一个负责任的工程师需要建立一条判断链：我到底需要矩阵本身，还是只需要它作用在某个向量上的结果？

1.1 条件数：贯穿全文的度量尺

条件数定义为 κ(A) = ‖A‖·‖A⁻¹‖，在 2-范数下等于最大奇异值与最小奇异值之比 σ_max/σ_min。它衡量的是：输入端的相对扰动，会被放大多少倍传递到输出端。对于求解 Ax = b，相对误差满足经典界（Higham, 2002）：

‖Δx‖/‖x‖ ≲ κ(A) · (‖ΔA‖/‖A‖ + ‖Δb‖/‖b‖)

这条不等式是全文的分析主线。它告诉我们：任何线性代数运算的精度预算，本质上都是条件数的预算。当 κ(A) ≈ 10¹⁶ 时，双精度（约 16 位有效数字）下求解结果可能一位有效数字都不剩。这不是库的 bug，而是数学的硬约束。

一个实用的经验法则：双精度下，若 κ(A) > 10¹²，求逆和解方程的结果基本不可信；若 κ(A) > 10⁸，需要警惕；κ(A) < 10⁴ 才算"良态"。这个阈值不是绝对的，取决于你对结果精度的实际需求。

1.2 四件套的复杂度与稳定性速查

运算	典型算法	浮点复杂度	稳定性特征
inv	LU + 三角回代	≈ 2n³/3	误差随 κ(A) 放大
det	LU 对角元乘积	≈ 2n³/3	易上溢/下溢
rank	SVD 奇异值计数	≈ 4n³（含 SVD）	依赖阈值选择
eig	QR 迭代（Hessenberg 化）	≈ 10n³ 量级	非对称矩阵特征值可剧烈敏感

注：复杂度系数为数量级估计，具体常数因实现（LAPACK 例程、BLAS 后端、分块策略）而异。表中数据为算法分析中的经典结论，非实测数据。

二、inv：求逆的代价与替代路径

2.1 一行代码做了什么

import numpy as np
A = np.array([[4.0, 7.0], [2.0, 6.0]])
A_inv = np.linalg.inv(A)   # 一行求逆

这行代码在底层调用 LAPACK 的 getrf（LU 分解，带部分主元）和 getri（基于 LU 的逆矩阵计算）。对于 n×n 稠密矩阵，标准实现约需 2n³/3 次浮点乘加。但真正的问题不在速度，而在精度。

本文评述：求逆是线性代数中最被滥用的操作。教科书里写 x = A⁻¹b，但工程实现中几乎永远不该显式构造 A⁻¹。原因有三：其一，显式求逆的运算量是直接解方程的约 3 倍；其二，求逆后再乘 b 引入两次舍入误差，而直接解只引入一次；其三，当 A 稀疏时，A⁻¹ 通常稠密，存储和计算双重爆炸。

2.2 正确姿势：解方程而非求逆

# 反例：显式求逆
x_bad = np.linalg.inv(A) @ b

# 正例：直接求解（内部走 LU + 回代）
x_good = np.linalg.solve(A, b)

# 多右端项：一次分解，多次回代
X = np.linalg.solve(A, B)   # B 的每一列是一个右端项

当需要反复求解同一系数矩阵、不同右端项时，更高效的做法是显式做一次 LU 分解，然后复用：

from scipy.linalg import lu_factor, lu_solve
lu, piv = lu_factor(A)          # 分解一次，O(n^3)
x1 = lu_solve((lu, piv), b1)    # 回代，O(n^2)
x2 = lu_solve((lu, piv), b2)    # 复用分解

SciPy 官方文档明确指出，lu_factor + lu_solve 的组合在重复求解场景下显著优于反复调用 solve，因为后者每次都重新分解。

2.3 什么时候真的需要逆矩阵

确实存在需要显式逆的场景，但比想象中少：

• 协方差矩阵的逆用于马氏距离、高斯判别分析。此时若矩阵接近奇异，应改用伪逆或加正则项 (A + λI)⁻¹。
• 迭代算法中的预条件子，如牛顿法更新 x ← x − J⁻¹f，实践中用 solve 而非 inv。
• 需要逆矩阵的对角元（如方差估计），可用 inv 但更推荐 scipy.linalg.inv 配合 check_finite=False 提升性能。

本文评述：判断标准很简单——如果你接下来要做的是 A_inv @ something，那你几乎一定不需要 inv。这个"something"应该直接进入 solve。只有当逆矩阵本身作为最终产物被存储、传输或进一步做元素级运算时，显式求逆才有意义。

2.4 伪逆：奇异与病态矩阵的退路

当 A 奇异或接近奇异时，inv 会抛出 LinAlgError 或返回数值垃圾。此时 Moore-Penrose 伪逆 A⁺ 是标准退路，它通过 SVD 构造：A⁺ = V Σ⁺ Uᵀ，其中 Σ⁺ 把非零奇异值取倒数、零奇异值保持为零。

A_pinv = np.linalg.pinv(A, rcond=1e-15)   # rcond 控制截断阈值

NumPy 的 pinv 默认 rcond=1e-15，即把小于 σ_max · 1e-15 的奇异值视为零。这个阈值直接决定了伪逆的"有效秩"，是工程调参的关键旋钮。

三、det：行列式的数值陷阱与对数化改造

3.1 行列式的数学意义与计算路径

行列式 det(A) 的教科书定义是 Leibniz 展开式，涉及 n! 项求和——这是指数级复杂度，仅适合理论推导。实际计算走 LU 分解：det(A) = det(P) · ∏ᵢ uᵢᵢ，其中 uᵢᵢ 是上三角矩阵的对角元，det(P) 是置换矩阵的符号（±1）。

d = np.linalg.det(A)   # 内部：LU 分解 + 对角元乘积

本文评述：行列式在数值计算中是一个"理论优美、实践危险"的量。它的危险来自两个方面：其一，数值范围极端——一个 100×100 的矩阵，若每个元素量级为 2，行列式量级可达 2¹⁰⁰ ≈ 10³⁰，远超双精度上界 1.8×10³⁰⁸ 的安全区；其二，它对矩阵元素的扰动极其敏感，与条件数同阶放大。

3.2 上溢与下溢：真实可复现的失败

import numpy as np
np.random.seed(0)
A = np.random.rand(200, 200) * 2
print(np.linalg.det(A))   # 可能返回 0.0 或 inf，取决于量级

上述代码为模拟示例。对于 200 维、元素在 [0,2) 的随机矩阵，行列式的对数期望约为 n·E[log λ]，实际数值极易超出双精度表示范围，导致 det 返回 0.0（下溢）或 inf（上溢）。这不是 bug，是浮点表示的物理极限。

3.3 对数行列式：工程标准解法

几乎所有涉及行列式的工程场景（高斯过程、变分推断、概率图模型）真正需要的是 log det(A)，而非 det(A) 本身。SciPy 提供了直接接口：

from scipy.linalg import lu_factor
import numpy as np

def logdet(A):
    lu, piv = lu_factor(A)
    diag = np.diag(lu)
    # 置换符号：piv 中发生交换的次数决定符号
    sign = np.prod(np.sign(diag))
    return np.sum(np.log(np.abs(diag))), sign

log_d, s = logdet(A)

NumPy 也提供了 np.linalg.slogdet，直接返回符号与对数绝对值：

sign, logabsdet = np.linalg.slogdet(A)
# det(A) = sign * exp(logabsdet)

本文评述：slogdet 是 NumPy 中最被低估的函数之一。它把"可能溢出"的乘法转化为"稳定的"加法，代价只是多返回一个符号。在高斯过程的对数边际似然计算中，用 slogdet 替代 det 是避免数值崩溃的第一步。

3.4 行列式与可逆性的关系：一个常见误解

"det(A) = 0 意味着矩阵不可逆"——这句话在数学上正确，在数值上几乎无用。一个条件数为 10¹⁵ 的矩阵，其行列式可能是一个正常的非零小数，但求逆结果毫无意义。反之，一个行列式因下溢显示为 0 的矩阵，可能实际条件数并不大。

本文评述：判断可逆性应该看条件数或奇异值，而不是行列式。行列式是"体积缩放因子"，它混合了所有方向的信息；而可逆性取决于是否存在接近零的奇异值。前者是一个标量，后者是一个谱分布——信息量不对等。

四、rank：秩判定的阈值哲学

4.1 数值秩的定义

数学秩是精确概念：线性无关列的最大数目。数值秩则依赖阈值：给定容差 tol，数值秩是大于 tol 的奇异值个数。NumPy 的实现：

r = np.linalg.matrix_rank(A)          # 默认阈值
r = np.linalg.matrix_rank(A, tol=1e-8) # 自定义阈值

NumPy 文档说明，默认阈值 tol = S_max · max(M, N) · eps，其中 S_max 是最大奇异值，eps 是浮点精度（双精度约 2.22×10⁻¹⁶）。这个默认值来自 SVD 的数值误差分析。

4.2 阈值选择的工程含义

阈值策略	适用场景	风险
默认（机器精度）	纯数学验证、精确输入	对噪声数据过于宽松
相对阈值（如 1e-6）	含噪测量数据	可能低估真实秩
基于间隙（gap）	奇异值有明显断层	断层不明显时失效

本文评述：秩判定本质上是一个聚类问题——把奇异值分成"显著"和"可忽略"两类。阈值就是聚类边界。当奇异值谱平滑衰减、没有明显间隙时，任何阈值都是人为选择。此时更稳健的做法是报告完整的奇异值谱，或者用"有效秩"（如累积能量占比 99%）替代硬阈值。

4.3 低秩近似：比判秩更有用的操作

工程中真正需要的往往不是"秩是多少"，而是"用多少秩能近似到什么程度"。Eckart-Young-Mirsky 定理给出了答案：截断 SVD 是 Frobenius 范数和谱范数下的最优低秩近似。

U, S, Vt = np.linalg.svd(A, full_matrices=False)
k = 10
A_k = U[:, :k] @ np.diag(S[:k]) @ Vt[:k, :]   # 最优秩-k 近似
err = np.linalg.norm(A - A_k, 'fro')           # 近似误差

这个误差有闭式表达：‖A − A_k‖_F = sqrt(Σ_{i>k} σᵢ²)。也就是说，你可以在不构造 A_k 的情况下，直接从奇异值谱读出任意 k 的近似误差。

五、eig：特征值分解的敏感性与非对称困境

5.1 一行代码的两种命运

w, v = np.linalg.eig(A)     # 一般矩阵，可能返回复数
w, v = np.linalg.eigh(A)    # 对称/Hermitian，返回实数，更快更稳

这两个函数的选择是特征值计算中最重要的决策。eigh 利用对称性，把问题化为三对角矩阵的特征值问题，用更稳定的算法（如 MRRR 或分治法），保证实特征值和正交特征向量。eig 则要处理一般矩阵，走 Hessenberg 化 + QR 迭代，复杂度和不稳定性都更高。

本文评述：如果你的矩阵在数学上是对称的，但数值上因为浮点误差有微小不对称（例如 A = BᵀB 计算后），应该先做 (A + Aᵀ)/2 对称化，再用 eigh。直接用 eig 会得到复数特征值，且虚部大小反映的是数值噪声而非真实信息。

5.2 非对称矩阵的特征值敏感性

对称矩阵的特征值对扰动是 Lipschitz 连续的：扰动 ‖ΔA‖ = ε 导致特征值变化不超过 ε。但非对称矩阵没有这个保证。经典反例是 Wilkinson 矩阵和以下 2×2 例子：

A = np.array([[0, 1], [0, 0]])       # 特征值均为 0
A_pert = np.array([[0, 1], [1e-10, 0]]) # 扰动 1e-10
# 特征值变为 ±1e-5，放大了 1e5 倍

这个例子中，矩阵元素的扰动是 10⁻¹⁰，但特征值从 0 变为 ±10⁻⁵，放大了五个数量级。原因在于该矩阵是 defective 的（不可对角化），特征值对扰动的敏感度是 ε^(1/m) 量级，m 是 Jordan 块大小。

本文评述：这解释了为什么非对称矩阵的特征值分解在工程中要慎用。如果你的矩阵来自真实测量，且没有对称性保证，特征值可能只是"看起来收敛"的数值幻象。此时更稳健的替代是 Schur 分解（scipy.linalg.schur）或奇异值分解——后者对任意矩阵都有良好的扰动理论。

5.3 广义特征值问题

工程中经常遇到 Ax = λBx 形式，如结构动力学、广义 Rayleigh 商。SciPy 提供：

from scipy.linalg import eigh
w, v = eigh(A, B)   # 求解 A x = λ B x，A、B 对称，B 正定

注意：eigh(A, B) 要求 B 正定。若 B 只是半正定，需要改用 QZ 分解（scipy.linalg.eig(A, B)），但后者不保证实特征值。

5.4 特征值 vs 奇异值：别搞混

这是初学者最常见的混淆点。特征值定义在方阵上，可能为负、可能为复数；奇异值定义在任意矩阵上，恒为非负实数。对于对称正定矩阵，两者相等；对于一般矩阵，σᵢ = sqrt(λᵢ(AᵀA))。

属性	特征值	奇异值
适用矩阵	方阵	任意 m×n
取值范围	复数域	非负实数
扰动稳定性	非对称时可能极差	始终良好
典型算法	QR 迭代	Golub-Kahan 双对角化

本文评述：当你不确定该用哪个时，优先选 SVD。SVD 对任意矩阵存在，数值稳定，且能同时给出秩、条件数、低秩近似和最小二乘解。代价是计算量约为特征值分解的 2 倍，但这个代价通常值得。

六、四件套的联动：从分解到应用

6.1 最小二乘：为什么用 QR 而不是正规方程

求解超定系统 Ax ≈ b 的标准方法是正规方程 AᵀAx = Aᵀb，但这会把条件数平方：κ(AᵀA) = κ(A)²。若 κ(A) = 10⁸，正规方程的条件数就是 10¹⁶，双精度下直接失效。

# 反例：正规方程
x_bad = np.linalg.solve(A.T @ A, A.T @ b)

# 正例：QR 分解（NumPy 的 lstsq 内部走 SVD 或 QR）
x_good, residuals, rank, sv = np.linalg.lstsq(A, b, rcond=None)

本文评述：np.linalg.lstsq 是四件套之外最该掌握的函数。它一次性返回解、残差、秩和奇异值，是诊断线性系统健康状况的"体检报告"。看到 lstsq 返回的 rank 小于列数，就该警惕共线性问题。

6.2 PCA：SVD 与协方差特征值分解的等价性

主成分分析有两种实现路径：对中心化数据矩阵 X 做 SVD，或对协方差矩阵 C = XᵀX/(n−1) 做特征值分解。数学上等价，数值上差异显著。

Xc = X - X.mean(axis=0)

# 路径 A：SVD（推荐）
U, S, Vt = np.linalg.svd(Xc, full_matrices=False)
components = Vt
explained_var = S**2 / np.sum(S**2)

# 路径 B：协方差特征值分解
C = Xc.T @ Xc / (Xc.shape[0] - 1)
w, v = np.linalg.eigh(C)   # 注意 eigh 返回升序，需反转

本文评述：路径 B 的条件数是路径 A 的平方，且需要显式构造 p×p 协方差矩阵（p 为特征维度）。当 p 很大时，路径 B 的内存和计算成本都更高。scikit-learn 的 PCA 实现默认走 SVD 路径，这是有充分数值理由的工程决策。

6.3 线性方程组诊断清单

面对一个线性系统，建议按以下顺序检查：

1. 检查条件数：np.linalg.cond(A)，若大于 10¹²，考虑正则化。
2. 检查秩：np.linalg.matrix_rank(A)，若小于 n，系统欠定或奇异。
3. 选择求解器：良态方阵用 solve；超定用 lstsq；病态用 pinv 或加 Tikhonov 正则。
4. 验证残差：np.linalg.norm(A @ x - b)，与 κ(A)·eps·‖b‖ 比较。

七、工程实践：框架选型、性能与踩坑清单

7.1 四大框架的 API 对比

框架	求逆	行列式	秩	特征值
NumPy	linalg.inv	linalg.det / slogdet	linalg.matrix_rank	linalg.eig / eigh
SciPy	linalg.inv	linalg.det / slogdet	linalg.matrix_rank	linalg.eig / eigh / eigvalsh
PyTorch	torch.linalg.inv	torch.linalg.det / slogdet	torch.linalg.matrix_rank	torch.linalg.eig / eigh
JAX	jnp.linalg.inv	jnp.linalg.det / slogdet	jnp.linalg.matrix_rank	jnp.linalg.eig / eigh

本文评述：PyTorch 和 JAX 的 linalg 命名空间是对 NumPy 的刻意对齐，但有两个关键差异：其一，它们支持自动微分，det、eig 等操作可反传梯度；其二，它们支持批处理，输入可以是 (batch, n, n) 张量。在深度学习场景中，批处理能力往往比单次运算速度更重要。

7.2 性能基准的真实数据来源

关于 NumPy/SciPy 线性代数性能，可参考以下公开资源：

• NumPy 官方性能基准仓库：github.com/numpy/numpy/tree/main/benchmarks
• SciPy 性能文档：docs.scipy.org/doc/scipy/tutorial/linalg.html
• LAPACK 官方用户指南：netlib.org/lapack/lug/

本文评述：性能数据高度依赖硬件（CPU 型号、BLAS 后端、缓存大小）和矩阵规模。任何"inv 比 solve 慢 3 倍"的断言都只是数量级参考，实际项目应在自己的目标硬件上做基准测试，而不是照搬网上的数字。

7.3 踩坑清单

八、前沿演进：稀疏、GPU、自动微分与随机数值线性代数

8.1 稀疏矩阵：绝不做稠密分解

当矩阵规模超过 10⁴ 且稀疏度高于 99% 时，稠密算法（LAPACK）的内存和计算都不可接受。SciPy 的 sparse.linalg 提供迭代解法：

from scipy.sparse import csr_matrix
from scipy.sparse.linalg import spsolve, eigsh, svds

A = csr_matrix(A_dense)          # 转稀疏格式
x = spsolve(A, b)                 # 稀疏直接解（SuperLU）
w, v = eigsh(A, k=10, which='LM') # 前 10 个最大特征值
U, S, Vt = svds(A, k=10)          # 前 10 个奇异值

本文评述：稀疏场景下，"求全部特征值"通常是错误需求。工程中真正关心的是极端特征值（最大/最小）或谱的局部信息。Lanczos（对称）和 Arnoldi（非对称）迭代算法正是为此设计，复杂度从 O(n³) 降到 O(k · nnz)，其中 nnz 是非零元个数。

8.2 GPU 加速：cuSOLVER 与批处理

NVIDIA 的 cuSOLVER 库为 GPU 提供了 LAPACK 对应实现。在 PyTorch 中，torch.linalg 会自动调度到 cuSOLVER：

import torch
A = torch.randn(64, 512, 512, device='cuda')  # 批大小 64
A_inv = torch.linalg.inv(A)                    # 批量求逆
w = torch.linalg.eigvalsh(A)                   # 批量特征值

本文评述：GPU 线性代数的优势在于批处理吞吐量，而非单矩阵延迟。对于单个 100×100 矩阵，CPU 可能更快（因为数据传输开销）。只有当批大小足够大（通常 >32）时，GPU 才显现优势。这是工程选型中容易被忽视的一点。

8.3 自动微分：可微线性代数的挑战

在深度学习框架中，det、inv、eig 的反向传播公式并不平凡。例如 det 的梯度是 ∂det(A)/∂A = det(A) · A⁻ᵀ，当 det(A) 接近零时梯度爆炸或消失。

本文评述：可微线性代数的核心矛盾是——前向稳定不等于反向稳定。一个在前向计算中表现良好的分解，其反向传播可能因为需要求逆而引入新的不稳定性。PyTorch 文档中专门讨论了 torch.linalg 的梯度限制，值得仔细阅读。

8.4 随机数值线性代数：大矩阵的新范式

近五年，随机算法在数值线性代数中快速崛起。核心思想是用随机投影把大矩阵降维到小矩阵，在小矩阵上做精确分解。例如随机 SVD：

def randomized_svd(A, k, p=10):
    n = A.shape[1]
    Omega = np.random.randn(n, k + p)
    Y = A @ Omega
    Q, _ = np.linalg.qr(Y)          # 正交基
    B = Q.T @ A                     # 小矩阵
    U_hat, S, Vt = np.linalg.svd(B, full_matrices=False)
    U = Q @ U_hat
    return U[:, :k], S[:k], Vt[:k, :]

Halko、Martinsson 和 Tropp 在 2011 年的 SIAM Review 综述中系统阐述了这一框架，此后被 scikit-learn 的 randomized_svd 等实现广泛采用。其误差有概率界保证，当奇异值衰减快时，k + p 取略大于 k 即可达到接近最优的精度。

本文评述：随机算法的价值在于把 O(n³) 降到 O(n²k)，当 k ≪ n 时优势巨大。但它不是万能的——若奇异值谱衰减缓慢（如白噪声矩阵），随机 SVD 需要 k + p 接近 n 才能保证精度，此时不如用精确算法。

8.5 混合精度：fp16/bf16 下的线性代数

🔒 复制本站文章内容需登录并达到 L3。当前：未登录

理学 MATLAB