矩阵分解为何慢 —— BLAS 线程、缓存、内存布局

含有 SVD 或特征值分解的 op “不知为何很慢”时，原因往往不在算法，而在于 BLAS 的使用方式。本文不是 op 的说明，而是讲述怀疑这些地方的顺序，以及如何测量确认的教材。

全部是 2026-09-06 在我自己的机器上测得的值，出处仅限我自己的实测（24 个逻辑 CPU / OpenBLAS 0.3.31 / numpy 2.4.6 / Windows，停掉其他负载的状态，中位数）。换一台机器数字会变，但趋势的方向应当能够再现。

1. 最有效、却最容易被忽视 —— 线程的过度分配

OpenBLAS 和 MKL 默认都根据逻辑 CPU 数决定线程数。对矩阵乘法来说这是对的。对分解来说则不对。

np.linalg.svd(n x n, full_matrices=False)      单位 ms
n        1t      2t      4t      8t     24t   最快
48    0.165   0.195   0.266   0.344   0.651   1t
96    0.650   0.918   1.081   1.305   2.217   1t
384  19.331  23.572  21.686  28.090  36.243   1t
512   45.9    44.0    45.5    47.3   105.3    2t
1024 351.7   320.0   259.4   262.4   868.8    4t
2048 2675    1950    1600    1344    3263     8t

参考：矩阵乘法 a @ a 的趋势相反（768 时 1t 14.0ms / 8t 2.4ms = 快 5.8 倍）

在 svd / eigh / qr / pinv 这 4 种 × 7 个尺寸中，24 线程获胜的格子连 1 个都没有。因为分解内部的 GEMM 太小，同步的开销超过了计算量。

判断

• 觉得分解慢时，在怀疑算法之前先怀疑线程数。确认只需 1 行：如果 OPENBLAS_NUM_THREADS=1 python your_script.py 变快了，就是这个原因。

• 不要限制整个进程。 矩阵乘法反而会变慢。只在分解周围限制。

• 在 fullseye 中，fsthreads 会自动完成这件事。想关闭时用 FULLSEYE_BLAS_THREADS=off，想固定时用 FULLSEYE_BLAS_THREADS=1。

陷阱 —— 这不只是速度的问题

线程数一变，归约的顺序就会变，因此分解的结果会在低位比特上发生变化。

svd(64x64)   1t vs 24t   bitwise 一致 False   最大差 3.6e-16
eigh(256)    1t vs 24t   bitwise 一致 False   最大差 1.7e-12

线程数由逻辑 CPU 数决定，所以什么都不做的状态在每台机器上给出不同的值。4 核的 CI 与 24 核的开发机上，分解的低位比特不一致。这可能成为“本地能通过、在 CI 上却超出容差”的原因。如果需要可再现性，就显式固定线程数。

2. 看什么来决定上限 —— 短边。不是长边，也不是元素数

在 (m, n) 的 92 个格子（svd 与 lstsq）上，针对每个预测量测量了“只用该规则决定线程上限时损失的时间”。

规则                          损失时间   超过 1.2 倍的格子   最差
什么都不做（保持 24t）          3975ms       79/92 格        19.91x
短边 <512:1 <1024:4 else:8       303ms        5/92 格         1.41x
元素数 <512:1 <1024:4 else:8     235ms       14/92 格         2.50x
长边 <512:1 <1024:4 else:8       241ms       ---              3.89x

判断

按合计，元素数和长边获胜。即便如此仍采用短边。 因为比起合计的差 68ms，最差值的差（1.41 倍还是 2.50 / 3.89 倍）更重要。与其平均快 1 成，不如最差也只多 4 成，这样更便于运维。

道理：LAPACK 的工作量大致为 長辺 × 短辺²（长边 × 短边的平方），线程是否起作用由短边的分块宽度决定。如果按长边把 16384x3 的最小二乘算作“大矩阵”，就会允许 8 个线程，但实际上它是宽度为 3 的 QR。

换算的参考

| 短边 | 允许的线程数 | 依据 |

| --- | --- | --- |

| 小于 32 | 不处理 | 限制机制本身就要 2.4 µs。套在比它更快的分解上反而吃亏 |

| 32〜511 | 1 | 在此范围内没有多线程获胜的格子 |

| 512〜1023 | 4 | 与 1t 不相上下，在上端 4t 开始获胜 |

| 1024 以上 | 8 | 在 2048 时 8t 比 1t 快 2.0 倍。24t 在全范围内最差 |

3. 边长为 2 的幂时会因缓存冲突而变慢（相当于 #pragma pack 的话题）

leading dimension 为 2 的幂时，缓存组会发生冲突而变慢。稠密矩阵库故意错开 lda 来分配内存，原因就在于此。

svd 的单线程有效 GFLOPS（* 为 2 的幂）
n     255    256*    257    264      -> 只有幂次慢 13%
    27.9    26.1    30.4   30.9
n     511    512*    513    520      -> 慢 18%
    41.4    34.6    41.2   44.1
n    1023   1024*   1025   1032      -> 慢 19%
    43.9    35.9    44.4   44.5
n = 64 / 128 时没有差异（能放进缓存，所以不发生冲突）

判断

• 13〜19% 是真实存在的。 只用 2 的幂的尺寸做基准测试，就会把只在那里才慢的值记录为“正常性能”。

• 不过最快的线程数不变（255/256/257/264 全是 1t，1023/1024/1025/1032 全是 4t）。不影响线程方面的判断。

• 测量性能时，要故意偏离 2 的幂再加 1 个点。不偏离的话，“格子的分布”会有偏差 —— 像正方 1023 这样的点连 1 个都没有。

4. 行优先与列优先 —— 转换后再用，大多得不偿失

numpy 是行优先，LAPACK 是列优先，所以每次分解都会发生重排。人们很容易想“那就按列优先传入好了”，但把转换的开销算进去就输了。

lstsq                原样       转换后传入     一开始就列优先   仅转换
(65536,   36)         27.4ms    26.0 (0.95x)   17.6 (0.64x)     6.4ms
(16384,  128)         43.9      50.1 (1.14x)   39.0 (0.89x)     8.5
( 4096,  512)        414.5     676.5 (1.63x)  432.8 (1.04x)     8.8
( 1024, 1024)        172.3     178.6 (1.04x)  177.1 (1.03x)     3.6

先 np.asfortranarray 再传入，5 种形状中有 4 种变慢。因为拷贝要花 6〜9 ms，而节省的时间比这还少。

转置视图不是免费的吗？ —— 免费的只有视图

取 a.T                       0.10 µs     只是交换步长 = 真正免费
np.ascontiguousarray(a.T)    6.5-8.7 ms  转回来就是真正的拷贝
svd(A.T) 后交换 U,V：
  (16384,128) 慢 1.90 倍 / (4096,512) 1.11 倍 / (1024,1024) 1.04 倍
  （奇异值一致。数学上正确，却在运行时间上输了）

由于 A = (Aᵀ)ᵀ，交换 svd(A.T) 的 U 和 V 即可得到 svd(A)，但 LAPACK 会把它当作 (n, m) 的宽矩阵走另一条路径，因此变慢。

判断

• “列优先就快”只在已经是列优先的情况下才成立。

• 只有在“一开始就按列优先组装，而且是竖长”的情况下才有效。在用 np.stack(..., axis=1) 构造设计矩阵的地方，只需用 order="F" 分配内存，竖长时就能获得 3 成左右的收益（(65536,36) 时为 0.64 倍）。

• 除此之外不要动手。转置再转回的设计，开销大于收益。

5. 有意利用连续轴 —— 排序与顺序统计（与分解不同维度的话题）

前面讲的是 LAPACK，但也有积极利用“快的方向”的设计。numpy 是行优先，所以最后一个轴是连续的。同样的运算如果可以选择轴，就靠向连续轴。

np.sort（单位 ms）
形状             dtype     axis=-1(连续)  axis=0(跳跃)   比
(65536,   32)   float32        3.68           18.39        5.00
( 4096,  256)   float32        1.95            9.12        4.67
( 4096,  256)   float64        5.42           13.94        2.57
( 1024, 1024)   float64        5.91           12.27        2.08

陷阱 —— “用小类型就快”并不成立

即使把同一个数组量化，比较排序也不会变快，反而更慢。

(1024,1024) 的 np.sort(axis=-1)
  float64  5.91ms      float32  3.04ms
  uint16  21.97ms      uint8   21.03ms      <- 比 float64 慢 3.6 倍

numpy 针对 float 的排序走的是经过 SIMD 优化的路径，而小整数类型不在这条路径上（在此版本中）。把 dtype 变小本身并不能换来速度。

另 1 个陷阱：上表中只有 (256, 4096) 的 uint8 看起来 axis=0 更快（19.45 对 22.57）。这不是布局的效果，而是1 次排序的长度的效果： axis=0 是对长度 256 的序列排序 4096 次，axis=-1 是对长度 4096 的序列排序 256 次。比较轴时，要先确认长度有没有变。

量化真正换来的东西 —— 可以省去排序

收益不是来自类型，而是来自值域有限。uint8 只有 256 种取值，所以不排序、只计数就够了。

uint8 的整体排序
  n = 16,777,216   np.sort 355.8ms  ->  计数排序 66.7ms   5.33 倍
  n =  2,097,152   np.sort  43.8ms  ->  计数排序  8.5ms   5.17 倍

uint8 的中位数（顺序统计本来就不需要排序）
  n = 8,388,608    np.median 35.6ms -> 直方图 27.6ms  1.3 倍
  n = 1,048,576    np.median  6.4ms -> 直方图  3.7ms  1.7 倍

计数排序只需用 np.bincount(a, minlength=256) 计数，再用 np.repeat(np.arange(256, dtype=np.uint8), counts) 展开即可。中位数只要对累积和做 searchsorted，连展开都不需要。

判断

• 可以选择轴的运算（排序、累积、归约）靠向最后一个轴。2〜5 倍。

• 如果为了速度引入量化，就走计数的路，而不是比较排序。 对 uint8/uint16 直接调用 np.sort，会比 float 还慢。

• 顺序统计（中位数、百分位数、rank filter）在已量化的情况下，可以用直方图的累积和替代。排序就消失了。

没有效果的做法（记录在案）

在对 N 个 D 维描述子做“距离 + argmin”时，比较了把描述子放在行上还是列上，没有出现显著差异（4096×64 时 1.00 倍，8192×128 时 1.06 倍，2048×256 时 0.97 倍）。因为 einsum 吸收了步长。布局的效果仅限于“可以选择轴的运算”，在伴随缩并的运算中不会出现。

6. 怀疑的顺序

含有分解的处理很慢时，按这个顺序确认。越靠上越容易见效，确认成本也越低。

1. 线程数 —— OPENBLAS_NUM_THREADS=1 是否变快（1 行，几秒）

2. 矩阵的大小与形状 —— 短边是多少。竖长还是正方

3. 调用次数 —— 是 1 次调用就慢，还是把小矩阵调用了几万次（如果是后者，限制机制的开销反而更高。在把 3×3 的 eigh 运行 2 万次的路径上，每 1 次都调用 threadpool_limits 会慢 1.9 倍：58.9ms → 112.9ms）

4. 尺寸是否为 2 的幂 —— 基准测试的尺寸是否有偏

5. 内存布局 —— 走到这一步才考虑。而且转换的价值也不大

测量时的注意事项

• threadpoolctl.threadpool_limits(...) 每次都会重新清点已加载的共享库，所以 1 次要 312 µs。持有 1 个 ThreadpoolController 并反复使用则是 2.4 µs。在循环中使用时用后者，而且要在循环外放 1 次。

• 取中位数而不是最小值。最小值会捡到“碰巧快了的 1 次”，不能代表常态的速度。

• 先确认没有其他进程占着 CPU 再测量。本文初版的数字，是在我自己放着不管的 3 个进程占用 3 个核的状态下测得的，比率被高估了 2〜3 倍（48×48 时 83 倍 → 实际是 27 倍）。

7. 从症状到原因

| 症状 | 首先怀疑 | 确认方法 |

| --- | --- | --- |

| 核数多的机器反而更慢 | BLAS 线程过度分配 | 加上 OPENBLAS_NUM_THREADS=1 重新运行。变快了就是这个原因 |

| 本地快而 CI 慢（或者相反） | 逻辑 CPU 数不同 | 在两边输出 threadpoolctl.threadpool_info() 进行比较 |

| 本地能通过、在 CI 上却超出容差 | 线程数改变了归约顺序，低位比特发生变动 | 在两边固定为相同的线程数后重新运行。差异消失就是这个原因 |

| 大量运行小矩阵的处理，加了限制反而变慢 | 限制机制本身的开销（每 1 次 2.4 µs，直接调用 threadpool_limits 则为 312 µs） | 只在循环外放 1 次。短边小于 32 时什么都不做 |

| 竖长矩阵的最小二乘比预想的慢 | 按长边判断为“大矩阵”而允许了多线程 | 按短边（列数）判断。16384x3 是宽度为 3 的 QR |

| 只有基准测试异常慢 | 尺寸是 2 的幂，发生缓存组冲突 | 把 n 与 n±1 并排测量。出现 13〜19% 的差就是这个原因 |

| 改成列优先也没有变快 | np.asfortranarray 的拷贝开销超过了节省 | 单独测量仅转换的时间，与节省量比较 |

| 量化了 np.sort 却没有变快 | 小整数类型不在 SIMD 路径上 | 替换为用 np.bincount 的计数排序进行比较 |

| 堆叠中位数很慢 | np.median(cube, axis=0) 是非连续轴 + 完整的中位数计算 | 张数为奇数时用 np.partition，张数多时把轴移到连续一侧 |

相关

• fsthreads —— fullseye 中实现这一上限的层。分级表与实测记录在其 docstring 中。

• math_metrology / measurement_uncertainty —— 测得值的处理（本文讲的是 “测量之前怀疑速度的顺序”，对象不同）。

外部的一手资料

本文的数值全部是我自己的实测，从外部引用的值连 1 个都没有。以下是规格与实现本身的出处（2026-09-06 确认可访问）。

• threadpoolctl（线程上限的实现。BSD-3）— <https://github.com/joblib/threadpoolctl>

• LAPACK Users' Guide（分解的计算量与分块）— <https://www.netlib.org/lapack/lug/>

• numpy 的线性代数 API — <https://numpy.org/doc/stable/reference/routines.linalg.html>