在数值计算这条线上待久了会发现一个有意思的现象很多人在纸上写矩阵求逆写的是 (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%。操作浮点运算量主导项相对误差量级内存增量典型用途LU 分解 回代求解(\frac{2}{3}n^3 2n^2)(\kappa(A)\cdot\varepsilon)(n^2)一般方阵求解Cholesky 两次三角求解(\frac{1}{3}n^3 2n^2)(\kappa(A)\cdot\varepsilon)更优常数(n^2/2)对称正定显式求逆经 LU(\approx 2n^3)数倍于求解(2n^2)确需逆矩阵Householder QR 求解(\frac{4}{3}n^3)对病态更鲁棒(n^2)病态、超定SVD / 伪逆(\approx 4n^3) 至 (20n^3)可截断可控(3n^2) 以上秩亏、降噪Newton-Schulz 迭代每步 (4n^3)两次矩阵乘依赖初值二次收敛(2n^2)高并行、低精度起步表里的 (\kappa(A)) 是条件数(\varepsilon) 是机器精度双精度约 (2.2\times10^{-16})单精度约 (1.2\times10^{-7})半精度约 (9.8\times10^{-4})。这个误差估计不是吓唬人的条件数 (10^8) 的双精度矩阵解出来大概只剩 8 位有效数字同样的矩阵如果条件数到了 (10^{16})双精度下结果基本全是噪声。所以看到一个矩阵第一件事是估条件数而不是急着选算法。运算量这块也要有个直观感受。n1000 时(\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 的一半在 n5000 这种规模下能省掉将近半分钟。第二是内存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_finiteFalse) x lu_solve((lu, piv), b, check_finiteFalse) # 手动做行平衡再分解对量级差异大的矩阵很关键 scale np.abs(A).max(axis1) 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)那么 (Ax 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-SchmidtMGS运算量更小约 (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|_11}|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_{k1} 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, iters6, tol1e-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, dtypeA.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 精度分级策略与流水线设计在一个完整的求解流水线里不同阶段的精度需求是不一样的一刀切用双精度会浪费大量带宽。我一般按下面的方式分级阶段推荐精度理由数据输入与存储低精度或单精度省带宽和显存误差在后续被平均掉矩阵乘累加单精度或双精度累加误差会累积不能省分解LU/Cholesky单精度起视条件数升级主元选择对精度敏感三角回代与分解同精度混用会引入额外误差残差计算双精度精化的效果完全依赖残差的准确性精化修正双精度修正量小但必须准这套配置的核心思想是把精度花在误差会被放大的环节省在误差会被稀释的环节。残差计算必须高精度因为它本身就是用小差值算大信息任何舍入都会直接毁掉精化的效果。流水线设计上还有一个容易忽略的点同步点的位置。每次精度切换都意味着一次数据格式转换和可能的同步这在加速卡上代价不低。合理的做法是把同一精度的操作尽量聚在一起减少格式来回转换。我见过一个实现每做一次矩阵乘就转一次精度结果转换开销占了总时间的四成把算法本身省下来的全赔进去了。6. 排查实录结果不对和速度不对的两种现场6.1 结果不对按这个顺序查数值问题的排查最忌讳乱试必须有固定顺序。我自己的顺序是这样的第一步验证输入。检查矩阵有没有 NaN 或 Inf有没有全零行或全零列量级差异是否过大。全零列意味着矩阵奇异任何求逆方法都会失败这时候该做的是降维或者加正则项而不是换个算法继续试。量级差异大的话先做行平衡很多看起来病态的矩阵平衡之后条件数就降下来了。第二步估条件数。这个数字能直接告诉你结果的可信度。条件数 (10^{16}) 的双精度矩阵不管用什么算法都拿不到有意义的结果这时候要做的是正则化或者降维不是调参数。第三步算残差。用 (|Ax - b|/(|A||x| |b|)) 这个相对残差衡量。如果残差小但解和预期差距大说明是矩阵本身病态导致的多解问题不是算法问题。如果残差大才是算法实现有问题。第四步做小规模对照。把矩阵缩到 3×3 或 4×4用手算或者高精度库算一遍和你的实现对比。这一步能揪出绝大多数索引错误、转置错误和缩放方向错误。我的经验是八成以上的算法不对其实是实现里的索引或方向错了。第五步换算法交叉验证。用两个独立的方法比如 LU 和 SVD 伪逆算同一个问题结果一致才可信。如果两者差异很大说明至少有一个的稳定性不够通常是 LU 那边出了问题。6.2 速度不对先看 roofline 再看实现性能问题同样有固定顺序先定位瓶颈类型再优化代码否则很容易优化了不重要的部分。先算算术强度判断是算力受限还是带宽受限。如果是带宽受限优化方向是减少搬运合并访存、提高缓存命中率、减少临时数组、把中间结果留在片上。如果是算力受限优化方向是提高并行度、用更快的指令比如矩阵乘的专用指令、减少分支。然后看并行效率。用加速卡时检查批次大小是否足够大、每个线程的工作量是否平衡、有没有线程在等待同步。小矩阵批量求解时批次太小会导致调度开销占主导批次够大但每个线程做太多事又会让寄存器溢出。这个平衡点要实测。最后看有没有隐藏的同步和内存分配。这两项是最常见的时间黑洞而且在分析器里不一定显眼。我的习惯是在关键路径上禁用动态内存分配全部改成预分配复用同步点则尽量合并能异步的绝不阻塞。6.3 速查表症状最可能的原因处理方式解里出现 Inf 或 NaN遇见零主元或矩阵奇异开部分主元检查零行零列残差小但解明显不对矩阵病态多解加正则项或截断小奇异值残差随规模增大而暴涨条件数过高精度不够做行平衡升级精度加精化逆矩阵和转置搞混索引顺序或存储布局错用小规模对照验证结果与库函数差几个数量级缩放还原方向错检查缩放作用在矩阵还是未知量上GPU 版本比 CPU 还慢批次太小或访存未合并增大批次改数据布局为行主序连续迭代法不收敛预条件子不匹配或矩阵非正定换预条件改用适用更广的 Krylov 方法Newton-Schulz 发散初值缩放不够用范数乘积的倒数缩放检查首次残差更新后精度逐渐劣化增量更新误差累积周期性全量重算或改用稳定的更新形式半精度结果完全不可用条件数超出半精度能力先估条件数超阈值直接走高精度这张表里的每一条我都至少遇到过一次其中缩放还原方向错和半精度条件数超限这两条是在生产环境里出过问题的代价都不小。前者是因为测试用例恰好对称掩盖了方向错误后者是因为开发时的手工数据条件数很低上线后真实数据条件数高了三个数量级。我个人在实际操作中的体会是矩阵求逆这件事算法选择的重要性其实排不到第一排第一的是搞清楚问题的条件数和真实需求。很多时候花两天优化算法不如花两个小时估一下条件数、确认下游到底需要逆矩阵还是只需要解收益大得多。另外一个实用的习惯是把每次求解的条件数估计、分解类型、相对残差这三个数字记进日志出问题时这三个数字基本能直接指出方向比翻代码快得多。
排查与解决方案)
