Python数値計算を最大235倍高速化：1701×1701実データ共分散で反復χ²評価を実測

記事

IT・テクノロジー

Fujikix_Nova

2026/09/02 13:43

1701×1701の実データ共分散行列を使った反復χ²評価で、K=2000では repeated linear solve loop に対して batched Cholesky route が235.34倍高速でした。

科学技術計算やデータ解析では、単純な二乗誤差ではなく、観測点どうしの相関を含む「共分散行列」を使って誤差評価を行う場面があります。

代表的な計算は、次のような covariance-weighted chi-square 評価です。

χ² = rᵀ C⁻¹ r

ここで C は共分散行列、r は観測値とモデル値の差を表す残差ベクトルです。

この計算は、1回だけなら大きな問題にならないこともあります。しかし、パラメータ探索、最適化、MCMC、シミュレーション、グリッドサーチのように、同じ共分散行列 C を固定したまま、残差ベクトル r だけを何千回、何万回と変えて評価する場合、実行時間に大きな差が出ます。

今回、Nova Statistical Core の数値計算検証として、Pantheon+ の実データ由来 1701×1701 dense covariance matrix を用い、反復 χ² 評価の高速化を実測しました。

これは宇宙論パラメータ推定そのものではありません。あくまで、大規模な密共分散行列を使う線形代数計算核のベンチマークです。

検証条件は以下です。

・実データ共分散：Pantheon+ 1701×1701 dense covariance

・行列サイズ：1701×1701

・rank：1701

・condition number：約3143.7

・route 判定：ACCEPT_CHOLESKY

・残差ベクトル：synthetic residual vectors

・比較対象：repeated linear solve loop / cached Cholesky / batched Cholesky matrix-RHS

・宇宙論fit：なし

・Hubble tension に関する主張：なし

まず、同じ残差ベクトル群に対して、毎回 `np.linalg.solve` を行う方式と、Cholesky factor を使ってまとめて評価する方式を比較しました。

K=1000 の直接実測では、結果は以下でした。

repeated linear solve loop：

54.3209 秒

batched Cholesky matrix-RHS：

0.4559 秒

高速化倍率：

119.14倍

次に、K=2000 でも同じ比較を行いました。

repeated linear solve loop：

99.8038 秒

batched Cholesky matrix-RHS：

0.4241 秒

高速化倍率：

235.34倍

この結果から、同じ共分散行列を使って多数の残差ベクトルを評価する場合、毎回 solve を行う実装と比べて、batched Cholesky route が非常に大きな時間短縮を生むことを確認しました。

重要なのは、単に速いだけではなく、数値結果が一致している点です。

K=1000 では、repeated solve と batched Cholesky の最大相対差は、

1.05 × 10⁻¹⁴

K=2000 では、

7.05 × 10⁻¹⁵

でした。

これは倍精度浮動小数点の丸め誤差レベルの一致です。

さらに、batch size optimization も行いました。

total K=10000、つまり1万回の χ² 評価を行う条件で、複数の batch size を比較しました。

比較した batch size は、

50, 100, 200, 300, 500, 1000, 2000, 5000, 10000

です。

この環境では、batch size=2000 が最良でした。

K=10000 の処理時間：

0.6035 秒

1評価あたりの時間：

0.00006035 秒

スループット：

約16,570 evaluations/sec

つまり、1701×1701 の実データ共分散を使った χ² 評価を、1万回まとめて処理しても、約0.6秒で完了しました。

また、内部実装ルートの比較も行いました。

scipy.linalg.cho_solve batched：

0.1359 秒

scipy.linalg.solve_triangular two-pass batched：

0.1036 秒

K=1000 の条件では、two-pass の solve_triangular route が cho_solve よりも速く、数値差は 0.0 でした。

メモリ使用量の見積もりも確認しました。

K=100000、batch size=2000 の場合でも、推定ピーク作業メモリは約96 MiBでした。

K=100000、batch size=10000 の場合でも、推定ピーク作業メモリは約304 MiBでした。

つまり、すべての残差を巨大行列として一括保持するのではなく、batch block 単位で処理することで、速度とメモリの両方を制御できます。

さらに、不正な共分散行列に対する fail-fast guard も検証しました。

nonfinite、non-symmetric、rank-deficient、indefinite、ill-conditioned といったケースでは、反復ループに入る前に遮断または警告できることを確認しました。

これは「速いだけでなく、無駄な計算に入る前に止められる」ことを意味します。

最後に、実運用に近い形として、synthetic model loop も実測しました。

これは宇宙論モデルではなく、軽量な代数モデルを使った模擬パイプラインです。

K=300 の条件で、

naive model + repeated solve：

13.8692 秒

batched model + batched Cholesky：

0.0924 秒

高速化倍率：

150.08倍

数値差：

1.98 × 10⁻¹⁵

という結果でした。

今回確認したことをまとめると、以下です。

・1701×1701 の実データ dense covariance を使用

・source SHA 一致を確認

・rank、固有値、条件数を診断

・正定値行列として Cholesky route を選択

・K=1000 で 119.14倍の直接実測高速化

・K=2000 で 235.34倍の直接実測高速化

・K=10000 を約0.6035秒で処理

・最大相対誤差は 10⁻¹⁴〜10⁻¹⁵ 台

・batch size optimization を実施

・solve_triangular route の高速性も確認

・memory footprint を確認

・fail-fast guard による安全遮断も確認

・synthetic end-to-end model loop で150.08倍の高速化を確認

このような計算は、宇宙論だけでなく、時系列解析、金融リスク解析、地球観測データ、分光データ、材料・電池データなどでも発生します。

特に、

「同じ共分散行列を使い、残差ベクトルだけを大量に変えながら評価する」

という構造では、アルゴリズムと実装ルートの選択が、計算時間に大きな差を生みます。

注意点として、この結果は同一マシン上でのローカルベンチマークであり、速度はCPU、メモリ帯域、BLAS/LAPACK backend、他プロセスの負荷に依存します。また、今回の結果は線形代数計算核および synthetic model loop の検証であり、宇宙論fit、Hubble tension、特定の物理結論を主張するものではありません。

Nova Statistical Core では、単に速いコードを書くのではなく、

・source lock

・行列診断

・rank check

・condition number check

・route selection

・数値一致検証

・timing benchmark

・memory footprint audit

・fail-fast guard

・claim boundary 管理

まで含めて、第三者に説明できる形で計算結果を整理します。