工程矩阵求逆全指南：LU/QR、迭代法与硬件加速

在数值计算这条线上待久了，会发现一个有意思的现象：很多人在纸上写矩阵求逆，写的是 (A^{-1} = \frac{A^*}{|A|})；而在工程里真正落地的时候，几乎没人这么干。原因不复杂——伴随矩阵法的复杂度是 (O(n!)) 级别的行列式展开，n 稍大一点就彻底不可用了。工程上说的"矩阵求逆"，实际上指的是一整套围绕分解、迭代、结构复用和硬件映射建立起来的方案体系，核心诉求只有两个：数值上别崩，时间上别炸。这篇内容就是把这两条诉求拆开，讲清楚四类主流方法各自的适用边界、稳定性来源、参数怎么定、坑在哪里。不管你是刚接手一个最小二乘拟合的新手，还是已经在写批量求解 kernel 的老手，应该都能从中找到可以直接抄走的部分。全文按"选型逻辑—分解法—正交化—迭代与结构复用—硬件加速—排查"的顺序展开，中间穿插了我自己在实际项目里踩过的坑和量化过的数据。

1. 把"求逆"当成一次工程选型，而不是一次公式代入

1.1 先问一句：你到底需要逆矩阵，还是只需要解

这是所有讨论的起点，也是我见过最多人栽跟头的地方。大部分工程问题的真实形式是 (Ax = b)，需要的是 (x)，而不是 (A^{-1}) 这个矩阵本身。这两件事在计算上的代价差了一个数量级：对 n 阶方阵，一次 LU 分解是 (\frac{2}{3}n^3) 次浮点运算，回代求解是 (2n^2)；而显式求逆大约是 (2n^3) 加上额外的数据搬运。换句话说，先求逆再相乘，运算量是直接求解的三倍左右，而且精度更差。

精度为什么会更差？因为显式求逆相当于把 (A^{-1}) 的每一个元素都算到机器精度，然后拿去做乘法再引入一次舍入。而直接求解只把误差集中在最终的解向量上。业内的经验法则是：显式求逆的相对误差大致是直接求解的若干倍，条件数越大，差距越明显。所以只要你的下游只需要 (x)，就不要生成逆矩阵。

那什么时候真的需要逆？我能想到的典型场景有三类。第一类是逆矩阵本身有物理或统计含义，比如协方差矩阵的逆就是精度矩阵，要在后续做不确定性量化。第二类是要对同一个 (A) 反复求解成千上万个不同的 (b)，这时候分解一次、复用多次是合理的，但仍然不需要显式求逆，缓存分解结果就够了。第三类是低维矩阵（n 小于几十）嵌在内层循环里，求逆和解方程的差异可以忽略，代码可读性优先。除此之外的绝大多数情况，都应该写成"分解 + 求解"的形式。

还有一种中间态值得提一句：如果你需要的是逆矩阵的对角线（比如在计算杠杆值、方差膨胀因子时），那也没必要全量求逆。逆的对角元可以通过分解后的三角求解逐列获取，代价是 (O(n^3)) 但常数小得多，内存占用从 (n^2) 降到可以逐列回收。这个技巧在统计计算里很常用，但很多新手不知道。

1.2 代价与精度的双重账本

选型不能凭感觉，得先有个量化的账本。我把常见操作的复杂度和精度特性列在下面，这张表基本能覆盖日常决策的 80%。

表里的 (\kappa(A)) 是条件数，(\varepsilon) 是机器精度（双精度约 (2.2\times10^{-16})，单精度约 (1.2\times10^{-7})，半精度约 (9.8\times10^{-4})）。这个误差估计不是吓唬人的：条件数 (10^8) 的双精度矩阵，解出来大概只剩 8 位有效数字；同样的矩阵如果条件数到了 (10^{16})，双精度下结果基本全是噪声。所以看到一个矩阵，第一件事是估条件数，而不是急着选算法。

运算量这块也要有个直观感受。n=1000 时，(\frac{2}{3}n^3) 约等于 (6.7\times10^8) 次浮点运算，在一颗能跑 10 GFLOPS 的核上大约 0.07 秒；n 翻倍到 2000，运算量变成 8 倍，直接到 0.5 秒以上。这也是为什么大规模问题必须往并行和硬件加速上走，算法本身已经没有多少优化空间了。

1.3 四类方法的边界划分

我习惯把工程上的求逆方案分成四类，划分依据不是数学上的分类，而是"它在什么条件下更稳、更快"。这四类的思路完全不同，混用容易出问题。

第一类是分解法，包括 LU、Cholesky、LDL 分解等，本质是把矩阵拆成容易求逆的三角或对角结构，稳定性靠主元策略保证。第二类是正交化路线，包括 QR 和 SVD，用正交变换把问题变成良态的形式，代价高但鲁棒性强，是病态矩阵的兜底方案。第三类是迭代与结构复用，包括各类定常迭代、Krylov 方法、Newton-Schulz 迭代，以及分块、Schur 补、低秩更新这类"只改一点点"的技巧。第四类是硬件加速，它本身不是一种数学方法，而是把前三类方法映射到特定计算单元上的工程手段，重点是精度分级、批次组织和访存优化。

前两类是"一次性把矩阵拆干净"，适合中等规模、需要高精度、可以接受 (O(n^3)) 的场景。第三类是"不拆，或者只拆局部"，适合稀疏、超大规模、或者矩阵会连续微调的场景。第四类则是横切关注点，任何一类方法都可能需要它。下面几节我会逐类拆开讲。

2. 分解法：数值稳定的地基

2.1 LU 与 Cholesky 的取舍逻辑

分解法里最先要做出的判断是：能不能用 Cholesky。如果能，就别用 LU。原因有三点，都是实打实的差异。第一是运算量，Cholesky 只需 (\frac{1}{3}n^3)，是 LU 的一半，在 n=5000 这种规模下能省掉将近半分钟。第二是内存，Cholesky 只要存储下三角，内存减半，这在 GPU 上直接决定了能不能一次装进显存。第三是稳定性，对称正定矩阵做 Cholesky 是不需要主元的，过程天然稳定，不会出现主元过小导致的数值爆炸。

判断矩阵是否对称正定，实际工程里有几种做法。理论上可以做一遍 Cholesky 试探，如果遇到负数开根就失败退出，但这会浪费一次分解的代价。更常见的做法是看问题的来源：协方差矩阵、刚度矩阵、正规方程里的 (A^TA)、图拉普拉斯矩阵，这些天然是对称半正定的，加上一点对角正则就严格正定了。如果是 (A^TA) 这种形式，我通常会额外提醒一句：不要真的去构造 (A^TA) 再分解，那样会把条件数平方，精度损失一半。正确做法是对 (A) 直接做 QR 或 SVD。

LU 的适用面更广，任何非奇异方阵都能用，代价是要处理主元。这里有个实操细节：很多人写代码时用"无主元 LU"，在教科书例子上跑得很好，一到真实数据就崩。原因是真实矩阵经常出现主元位置上是零或者极小值的情况，不做行交换会直接除出无穷大。所以工业级实现里，LU 几乎总是带部分主元（partial pivoting）的。

2.2 主元策略：为什么"部分主元"几乎是默认答案

主元策略的选择本质上是稳定性与并行性的权衡，三种主流策略各有明确的使用场景。

部分主元（每列选绝对值最大的行交换到对角位置）是默认答案，因为它把增长因子控制在 (2^{n-1}) 这个理论上界内，实践中几乎不会触发。代价是需要行交换，这会破坏数据的局部性，在并行环境里还要维护一个置换向量。这个置换向量千万别丢，后续求解、求行列式、求逆都要用到它。

完全主元（行列都选最大）稳定性更好，但搜索成本是 (O(n^3)) 级别的额外比较，而且列交换会打乱变量顺序，恢复起来麻烦。除非你在处理极端病态的矩阵，否则不值得。

对角主元的出发点是保持稀疏性。对稀疏矩阵来说，任何行交换都可能引入大量的填充元（fill-in），让原本 (O(n)) 的非零元变成 (O(n^2))。这时候宁可牺牲一点数值稳定性，也要保住稀疏结构。这也是稀疏直接解法里各种重排序算法（最小度排序、嵌套剖分）存在的原因。

还有一个容易被忽略的点：矩阵缩放。如果矩阵各行各列的量级差异很大（比如一行是 (10^6)，另一行是 (10^{-6})），即使条件数不算太差，浮点表示和主元选择也会出问题。标准的做法是做行、列平衡，把矩阵缩放成行列范数接近的形式，解完再缩放回来。这个步骤在做数值计算库时是默认开启的，但自己手写实现时经常忘掉。

import numpy as np from scipy.linalg import lu_factor, lu_solve # 带部分主元的 LU，返回的 piv 是置换索引 lu, piv = lu_factor(A, check_finite=False) x = lu_solve((lu, piv), b, check_finite=False) # 手动做行平衡再分解，对量级差异大的矩阵很关键 scale = np.abs(A).max(axis=1) scale[scale == 0] = 1.0 A_bal = A / scale[:, None] lu2, piv2 = lu_factor(A_bal) y = lu_solve((lu2, piv2), b) x2 = y / scale # 注意解的还原方向

上面这段代码里最后一步的方向很容易搞错。我见过不止一次有人写成y * scale，在小规模测试上碰巧对，一上真实数据就全错。记法很简单：如果 (A = D \cdot A')，那么 (A'x = D^{-1}b)，解出 (x) 之后不需要再做任何变换。只有当缩放作用在未知量上时才需要还原，这两种情形一定要分清。

注意：对秩亏矩阵或严重病态矩阵，LU 分解不会报错，而是安静地给出一个巨大的解。所以分解完成后必须检查主元的最小绝对值，如果和最大主元差了几个数量级，就该转到 QR 或 SVD 路线了。

2.3 迭代精化：把丢掉的精度捡回来

迭代精化（iterative refinement）是我认为被严重低估的一个技巧，实现成本极低，收益却很直接。基本思路是：分解得到近似解 (x_0)，然后算残差 (r = b - Ax_0)，用同一个分解去解修正量 (Ad = r)，再令 (x_1 = x_0 + d)，重复几次。

关键在于残差的计算精度。如果 (r) 用和分解相同的精度算，精化效果有限；如果残差用更高的精度计算（比如分解用单精度、残差用双精度），效果会非常显著。这就是所谓的"混合精度迭代精化"，在双精度分解的场合对应的是"补偿求和"或者扩展精度累加。

收益有多大？对良态矩阵，一次精化可以把误差降低一到两个数量级；对条件数 (10^{10}) 左右的矩阵，两到三次精化基本能恢复到接近机器精度的水平。而代价只是每次 (O(n^2)) 的矩阵向量乘和一次三角回代，相对于 (O(n^3)) 的分解完全可以忽略。

实现上有个细节：判断收敛不要看修正量 (d) 的范数，要看残差 (r) 的范数相对于 (|A||x| + |b|) 的比值。因为 (d) 本身可能因为分母的病态而偏大，用它判断会误判。还有一个更实际的建议：设定最大迭代次数（3 到 5 次就够），不要写成 while 循环，避免在极端情况下卡死。

3. 正交化路线：病态矩阵的兜底

3.1 QR 求逆的实际代价

QR 分解的思路是把 (A) 写成 (QR)，其中 (Q) 是正交矩阵，(R) 是上三角。因为 (Q^{-1} = Q^T)，所以 (A^{-1} = R^{-1}Q^T)，问题就变成了上三角矩阵的求逆，这个有直接的递推公式。代价上，Householder QR 是 (\frac{4}{3}n^3)，大约是 LU 的两倍，但如果只求解不做显式求逆，回代部分和 LU 一样。

QR 值得用的理由是稳定性。正交变换不放大误差，所以 QR 的误差增长因子通常只和 (\kappa(A)) 成正比，而 LU 还要额外乘上主元增长因子。对条件数在 (10^6) 到 (10^{10}) 之间的矩阵，QR 给出的结果往往比 LU 多一到两位有效数字。这个差距在最小二乘、正交化、特征值计算里是决定性的。

QR 有一个变体叫修正 Gram-Schmidt（MGS），运算量更小（约 (2n^3)），在并行环境下也更好组织。但要注意，经典 Gram-Schmidt 在数值上是不稳定的，必须用修正版本，或者在 MGS 之后再补一次重正交化。我见过一些为了追求速度改用经典 GS 的实现，在正交性要求严格的场景（比如构造基向量）里直接失败。

3.2 SVD 伪逆与截断阈值的定量选择

SVD 是最后的兜底方案。(A = U\Sigma V^T)，伪逆就是 (A^+ = V\Sigma^+U^T)，其中 (\Sigma^+) 是把非零奇异值取倒数。它的独特价值在于能处理秩亏和病态：你可以主动丢掉那些小于阈值的奇异值，得到一个"最合理的近似逆"，而不是一个数值爆炸的结果。

阈值怎么定是最实际的问题。先用一个纯数学的下界：(\sigma_i / \sigma_{\max} < n\varepsilon) 的奇异值在数值上就是零，必须截断，否则倒数会变成 (10^{16}) 量级。但纯数学阈值往往不够，因为工程问题里的"小奇异值"可能代表真实但不可靠的方向。我的经验是分两种情形：

• 如果矩阵来自干净的数据、已知精确秩，用 (\text{tol} = \max(m,n)\cdot\varepsilon\cdot\sigma_{\max}) 就够，这是大多数线性代数库的默认做法。
• 如果矩阵来自带噪测量，阈值应该和噪声水平挂钩。一个可操作的估计是用奇异值谱的"肘部"位置：把 (\sigma_i) 降序排列，找相邻比值的突变点，或者用 (\sigma_i^2) 的累积能量曲线取 99% 或 99.9% 的截断。

截断的副作用必须说清楚：丢掉小奇异值意味着你放弃了那些方向上的信息，得到的解是有偏的。在反问题里这等价于正则化，是有意为之；但如果你只是想求个精确逆，截断就意味着结果不再是 (A^{-1})。所以做这一步之前，一定要确认下游能接受这个近似。

3.3 条件数估计：不用算 SVD 也能知道有多糟

算 SVD 的代价是 (O(n^3)) 且常数很大，为了估个条件数专门跑一次不划算。实际工程里用条件数估计器，基于 1-范数或无穷范数，通过几次矩阵向量乘就能给出一个量级正确的估计，代价只有 (O(n^2))。

原理是利用 (|A^{-1}|1) 的定义：它等于 (\max{|x|_1=1}|A^{-1}x|_1)。这个最大值不好精确求，但可以迭代搜索：取一个使 (|A^{-1}x|) 大的方向，解一次线性系统得到 (y = A^{-1}x)，用 (y) 的符号构造新的 (x)，重复几次，通常三四次就能收敛到真实值的几倍以内。这种估计在判断"要不要换算法"时完全够用。

实操建议：在预处理阶段就对矩阵做一次条件数估计，把结果打进日志。这个数字在后期排查问题时价值巨大，能直接告诉你异常结果是算法选错了还是数据本身有问题。我现在的习惯是在任何求解入口都记录 (\kappa) 的估计值、矩阵维度和实际使用的分解类型，出问题时看这三项基本能定位到方向。

4. 迭代与结构复用：大矩阵、稀疏矩阵和"只改一点点"

4.1 稀疏场景：别求逆，求分解

稀疏矩阵的第一原则是：永远不要显式求逆。原因很直接，稀疏矩阵的逆通常是稠密的。一个三对角矩阵的逆是全稠密的，一个来自网格离散的稀疏矩阵的逆，非零元比例会随规模增长迅速趋近 100%。所以对稀疏问题，要么用稀疏直接分解（LU、Cholesky 的稀疏版本，配合重排序控制填充），要么用迭代解法。

稀疏直接分解的核心矛盾是填充元。分解过程中原本是零的位置会变成非零，填充量取决于消元顺序。这就引出了重排序算法：最小度排序倾向于每次消去度数最小的节点，嵌套剖分则通过递归分割把填充限制在分割面上。选哪个要看矩阵结构：来自二维网格的问题，嵌套剖分通常更优；来自不规则图结构的问题，最小度更稳。

稀疏迭代解法则完全不碰分解，只做矩阵向量乘。它的收敛速度取决于矩阵的谱分布，对条件数差的矩阵需要预条件。Krylov 子空间方法（共轭梯度、GMRES、BiCGSTAB 这一类）是主流选择，其中共轭梯度只适用于对称正定，其他情形要换。这类方法的实现细节很多，但核心的一句忠告是：先看收敛曲线再调参数，如果残差曲线来回震荡，八成是预条件子选错了，而不是迭代次数不够。

4.2 Newton-Schulz：为并行而生的无分解迭代

Newton-Schulz 迭代是我在做高并行求逆时最常用的方案，它的形式极其简单：

[ X_{k+1} = X_k(2I - AX_k) ]

这个式子的来源是牛顿法解方程 (f(X) = X^{-1} - A = 0)，但不用记推导，记住性质就行：如果初值满足 (|I - AX_0| < 1)，迭代二次收敛，误差平方级下降。也就是说，误差从 (10^{-2}) 到 (10^{-4}) 到 (10^{-8}) 到 (10^{-16})，大概四到五次迭代就能到双精度极限。

它最迷人的地方是：整个迭代只有矩阵乘法，没有任何分解、没有三角回代、没有主元选择。矩阵乘法是 GPU 上优化最彻底的算子，能跑满算力，而且天然可以批处理。所以在大批量小矩阵求逆、或者需要在 GPU 上避免分支发散的场景，它比 LU 快得多。

初值怎么取是唯一的难点。最常用的构造是 (X_0 = \alpha A^T)，其中 (\alpha = 1/(|A|1|A|\infty))，这个缩放的目的是让 (AX_0) 的特征值落在 ((0, 2)) 区间内，保证收敛。也可以用更精细的估计，比如用 Frobenius 范数或者矩阵的迹来缩放。如果矩阵是正定的，还可以直接取 (X_0 = I/\text{tr}(A)) 之类的粗估，反正迭代会自己修正。

def newton_schulz_inv(A, iters=6, tol=1e-14): n = A.shape[0] # 用 1-范数和 inf-范数构造缩放因子，保证谱半径条件 alpha = 1.0 / (np.linalg.norm(A, 1) * np.linalg.norm(A, np.inf)) X = alpha * A.T.copy() I = np.eye(n, dtype=A.dtype) for k in range(iters): AX = A @ X # 残差范数就是收敛判据，比看 X 的变化更可靠 res = np.linalg.norm(I - AX, 1) if res < tol: break X = X @ (2.0 * I - AX) return X, k, res

这段代码里，2.0 * I - AX和最后一次乘法是两次 GEMM，每次 (2n^3) 运算量，所以单步是 (4n^3)，比 LU 的单次 (2n^3/3) 贵不少。但因为全是 GEMM，在加速卡上的实际吞吐能高出十几倍，综合下来仍然占优。另一个好处是它不需要转置和置换，内存访问模式非常规整。

注意：Newton-Schulz 对初值敏感，如果 (|I - AX_0| \ge 1)，迭代会发散而不是收敛。所以在生产代码里必须检查第一次迭代后的残差，如果它比初始残差还大，就要退回分解法。这个检查我一般会写死，不做成可选。

4.3 分块、Schur 补与低秩更新

现实中的矩阵经常是"大部分不变、小部分变了"。这种场景下重算整个逆是巨大的浪费，正确做法是利用结构做增量更新。

分块求逆是最基础的工具。把矩阵切成四块：

[ \begin{bmatrix} A & B \ C & D \end{bmatrix}^{-1} = \begin{bmatrix} S^{-1} & -S^{-1}BD^{-1} \ -D^{-1}CS^{-1} & D^{-1} + D^{-1}CS^{-1}BD^{-1} \end{bmatrix}, \quad S = A - BD^{-1}C ]

其中 (S) 就是 Schur 补。这个公式的价值在于：如果 (D) 已经有逆（或者已经分解好），那么更新只需要算 (S) 并求它的逆，规模从 (n) 降到 (n/2) 左右，运算量降到原来的四分之一。这在递归分块算法里是核心，也是并行分块求解的理论基础。

低秩更新（Woodbury 恒等式）解决的是另一类问题：矩阵只在一个低秩方向上被修改。

[ (A + UCV)^{-1} = A^{-1} - A^{-1}U(C^{-1} + VA^{-1}U)^{-1}VA^{-1} ]

如果 (U, V) 是 (n\times k) 且 (k \ll n)，那么中间的求逆只需要对 (k\times k) 的矩阵做，代价从 (O(n^3)) 降到 (O(n^2k))。典型的应用场景是：在线学习中协方差矩阵的增量更新、传感器网络里逐个加入观测、以及各种递推最小二乘。我用这个技巧做过一个实时更新的场景，单次更新从 40 毫秒降到 0.3 毫秒，效果非常直接。

但这里有个精度陷阱必须提醒：Woodbury 公式在数值上是不稳定的，因为 (A^{-1}U) 这些中间量可能放大误差，尤其当 (A) 本身病态时。稳定的做法是改用 QR 更新或者直接把新数据拼进去重新分解，只在 (k) 很小（比如小于 10）且 (A) 良态时才用 Woodbury。我踩过这个坑，一次更新误差累积到第三十几次时结果完全失真，最后是加了周期性全量重算才解决的。

4.4 预条件子：迭代法的真正门槛

迭代法的收敛速度由矩阵的谱决定，而预条件子的作用就是把谱"压紧"，让迭代次数和问题规模脱钩。没有预条件的 Krylov 方法，迭代次数通常随 (\sqrt{\kappa}) 增长，条件数 (10^6) 就要几千次迭代，完全不可接受。

常用的预条件子分几个层次。最简单的对角预条件（Jacobi），只取对角线，代价近乎为零，对对角占优的矩阵效果不错。进阶一点的是不完全分解（ILU、不完全 Cholesky），保留分解中的部分填充元，在稀疏结构上做近似分解。最强的是多重网格类的预条件，对来自偏微分方程离散的矩阵几乎能做到迭代次数与网格规模无关。

选预条件子的实操建议是：从最便宜的试起，逐步加码。先用对角预条件跑一下，看迭代次数是否可接受；不行换 ILU，调填充等级；再不行才考虑多重网格。不要一上来就上最复杂的方案，因为预条件本身也有构造代价，而且越复杂的方案对矩阵结构的假设越多，换个问题就可能失效。

还有一个常被忽略的点：预条件子必须是对称正定的（如果原方法要求），否则共轭梯度会失效。不完全 Cholesky 在丢弃填充元时可能产生负主元，需要做修正（比如对角偏移）。这个细节在实现里经常被跳过，导致原本该收敛的迭代随机发散。

5. 硬件加速：把算法映射到硅片上

5.1 小矩阵海量批次：batched 求解的收益曲线

在图形渲染、物理仿真、深度学习里，最常见的形态不是一个大矩阵，而是几万个 3×3、4×4、8×8 的小矩阵。这种场景下，单个矩阵的求逆用不了多少算力，瓶颈在于组织方式。

关键洞察是：单线程处理一个小矩阵浪费了绝大部分算力，而批处理能把访存和计算重叠起来。具体做法是把同一批次的所有矩阵连续存放，让每个线程或每个线程块负责一个矩阵，全批次一起推进。这样做的另一个好处是消除了分支发散——所有矩阵走完全相同的代码路径。

收益有多明显？在一个 3×3 求逆的场景里，我做过对比：朴素的循环调用标量代码，每秒能处理几十万个；改成批处理内核、每个线程处理一个矩阵、用寄存器存全部元素之后，吞吐提升到每秒上千万个，接近两个数量级。差异的来源不是算法，而是消除了函数调用开销、实现了访存合并、让计算单元一直有活干。

这里有个参数值得琢磨：每个线程处理几个矩阵。处理 1 个，并行度高但每个线程的计算量太少，调度开销占比大；处理 4 到 8 个，可以让指令级并行填满流水线，同时保持访存合并。这个数没有通用最优值，和矩阵尺寸、寄存器压力、缓存大小都有关。我的经验是 3×3 和 4×4 用 1 到 2 个，8×8 到 16×16 用 1 个，再大就要考虑把矩阵放进共享内存、用整个线程块协作。

5.2 混合精度 + 迭代精化：FP16 分解到底靠不靠谱

半精度的动态范围小、尾数只有 10 位，直接拿来做分解，结果基本没法用。但"分解用低精度、修正用高精度"这套组合，在很多场景下确实能跑通，前提是条件数不能太大。

逻辑是这样的：低精度分解给出的近似解误差量级约为 (\kappa(A)\cdot\varepsilon_{\text{half}})。如果 (\kappa(A) = 10^2)，误差约 (10^{-1})，虽然很大但还在"能作为迭代初值"的范围内；然后做两步高精度迭代精化，每步误差降低到原来的 (\kappa\varepsilon) 倍，两步之后就能逼近双精度的解。但如果 (\kappa(A) = 10^5)，低精度解已经是噪声，精化根本救不回来。

所以在决定用不用低精度之前，必须先知道条件数。这也是我在第 3 节强调条件数估计的原因——它不只影响算法选择，还直接影响精度分级策略。

实际实现的推荐路数是：输入用低精度存储（省带宽和显存）、累加用单精度或双精度（保精度）、分解用单精度（折中）、精化用双精度。这套配置在多数中等条件数的问题上能拿到接近双精度的结果，同时吞吐提升两到四倍。如果非要上低精度分解，务必加上"结果校验"环节：算一遍残差范数，超过阈值就退回高精度重算。

5.3 访存受限的现实：为什么求逆加速比比 GEMM 难看

很多人会疑惑：矩阵乘法在加速卡上能提速几十倍，为什么求逆只提速几倍？答案在算术强度上。

算术强度是浮点运算次数与访存字节数的比值。矩阵乘法是 (O(n^3)) 运算对 (O(n^2)) 数据，强度是 (O(n))，n 大的时候完全算力受限，加速卡的高算力能充分发挥。而矩阵求逆里的三角回代是 (O(n^2)) 运算对 (O(n^2)) 数据，算术强度只有 (O(1))，属于典型的访存受限。加速卡的算力再高，数据供不上也是白搭。这就是 Amdahl 定律在数值计算里的具体体现。

换算成数字：假设加速卡的峰值算力是 10 TFLOPS 级别，内存带宽是 1.5 TB/s，那么转折点的算术强度大约是 6.7 flops/byte。低于这个值的操作，性能上限由带宽决定，算力再涨也没用。三角求解的算术强度远低于这个值，所以它在加速卡上的实际利用率通常只有峰值的百分之几到十几。

这个认识带来两个实践上的推论。第一，不要指望靠换硬件把三角求解加速十倍，瓶颈在带宽不在算力，优化方向应该是减少数据搬运：把矩阵分块装进共享内存、复用已经读入的数据、避免不必要的临时数组。第二，优先选择算术强度高的算法形式。这也是 Newton-Schulz 在大规模并行场景下受欢迎的原因——它全程都是 GEMM，算术强度高，能把加速卡吃满。

5.4 精度分级策略与流水线设计

在一个完整的求解流水线里，不同阶段的精度需求是不一样的，一刀切用双精度会浪费大量带宽。我一般按下面的方式分级：

这套配置的核心思想是：把精度花在误差会被放大的环节，省在误差会被稀释的环节。残差计算必须高精度，因为它本身就是"用小差值算大信息"，任何舍入都会直接毁掉精化的效果。

流水线设计上还有一个容易忽略的点：同步点的位置。每次精度切换都意味着一次数据格式转换和可能的同步，这在加速卡上代价不低。合理的做法是把同一精度的操作尽量聚在一起，减少格式来回转换。我见过一个实现每做一次矩阵乘就转一次精度，结果转换开销占了总时间的四成，把算法本身省下来的全赔进去了。

6. 排查实录：结果不对和速度不对的两种现场

6.1 结果不对：按这个顺序查

数值问题的排查最忌讳乱试，必须有固定顺序。我自己的顺序是这样的：

第一步，验证输入。检查矩阵有没有 NaN 或 Inf，有没有全零行或全零列，量级差异是否过大。全零列意味着矩阵奇异，任何求逆方法都会失败，这时候该做的是降维或者加正则项，而不是换个算法继续试。量级差异大的话先做行平衡，很多"看起来病态"的矩阵平衡之后条件数就降下来了。

第二步，估条件数。这个数字能直接告诉你结果的可信度。条件数 (10^{16}) 的双精度矩阵，不管用什么算法都拿不到有意义的结果，这时候要做的是正则化或者降维，不是调参数。

第三步，算残差。用 (|Ax - b|/(|A||x| + |b|)) 这个相对残差衡量。如果残差小但解和预期差距大，说明是矩阵本身病态导致的多解问题，不是算法问题。如果残差大，才是算法实现有问题。

第四步，做小规模对照。把矩阵缩到 3×3 或 4×4，用手算或者高精度库算一遍，和你的实现对比。这一步能揪出绝大多数索引错误、转置错误和缩放方向错误。我的经验是，八成以上的"算法不对"其实是实现里的索引或方向错了。

第五步，换算法交叉验证。用两个独立的方法（比如 LU 和 SVD 伪逆）算同一个问题，结果一致才可信。如果两者差异很大，说明至少有一个的稳定性不够，通常是 LU 那边出了问题。

6.2 速度不对：先看 roofline 再看实现

性能问题同样有固定顺序，先定位瓶颈类型再优化代码，否则很容易优化了不重要的部分。

先算算术强度，判断是算力受限还是带宽受限。如果是带宽受限，优化方向是减少搬运：合并访存、提高缓存命中率、减少临时数组、把中间结果留在片上。如果是算力受限，优化方向是提高并行度、用更快的指令（比如矩阵乘的专用指令）、减少分支。

然后看并行效率。用加速卡时，检查批次大小是否足够大、每个线程的工作量是否平衡、有没有线程在等待同步。小矩阵批量求解时，批次太小会导致调度开销占主导；批次够大但每个线程做太多事，又会让寄存器溢出。这个平衡点要实测。

最后看有没有隐藏的同步和内存分配。这两项是最常见的时间黑洞，而且在分析器里不一定显眼。我的习惯是在关键路径上禁用动态内存分配，全部改成预分配复用；同步点则尽量合并，能异步的绝不阻塞。

6.3 速查表

这张表里的每一条我都至少遇到过一次，其中"缩放还原方向错"和"半精度条件数超限"这两条，是在生产环境里出过问题的，代价都不小。前者是因为测试用例恰好对称，掩盖了方向错误；后者是因为开发时的手工数据条件数很低，上线后真实数据条件数高了三个数量级。

我个人在实际操作中的体会是：矩阵求逆这件事，算法选择的重要性其实排不到第一，排第一的是搞清楚问题的条件数和真实需求。很多时候花两天优化算法，不如花两个小时估一下条件数、确认下游到底需要逆矩阵还是只需要解，收益大得多。另外一个实用的习惯是把每次求解的条件数估计、分解类型、相对残差这三个数字记进日志，出问题时这三个数字基本能直接指出方向，比翻代码快得多。