本指南僅有日文版(它是人寫的散文,這裡不做機械替換)。

行列分解が遅い理由の知識 — 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 かかる。ThreadpoolController を 1 個持って

使い回すと 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>