Pythonで共分散行列の逆行列を直接作らない方がよい理由

Pythonで科学計算やデータ解析をしていると、χ²計算、尤度計算、MCMC、回帰分析などで共分散行列を扱う場面があります。

このとき、よくある実装として

C^{-1}

を明示的に作り、

r.T @ inv(C) @ r

のように計算するコードがあります。

数式としては正しく見えますが、数値計算の実装としてはあまり推奨されません。理由は、逆行列を明示的に作ることで、計算コストや数値誤差が増えやすくなるためです。

より安定した方法は、共分散行列 C を Cholesky 分解し、

C = L L.T

として、残差ベクトル r に対して

L w = r

を解く方法です。

その上で、

χ² = w.T @ w

として計算します。

この方法は、明示的に C^{-1} を作らずに、同じ共分散重み付きχ²を評価できます。

Pythonでは、たとえば以下のように書けます。

mport numpy as np

from scipy.linalg import solve_triangular

def eval_whitened_chi2(L, y_obs, y_pred):

    r = np.asarray(y_obs) - np.asarray(y_pred)

    w = solve_triangular(L, r, lower=True, check_finite=False)

    return float(np.dot(w, w))

ここで L は、あらかじめ

L = np.linalg.cholesky(C)

として計算しておきます。

この方法の利点は以下です。

・明示的な逆行列を作らない  

・数値的に安定しやすい  

・同じ共分散行列を何度も使うMCMCや尤度計算に向いている  

・処理の意味が明確になる  

・第三者が確認しやすいコードになる  

特に、MCMCやベイズ推定では、同じ共分散行列に対して何千回、何万回もモデル評価を行うことがあります。このような場合、最初にCholesky分解を一度だけ行い、その後は三角行列solveを繰り返す構成にすると、実装が整理されます。

私のサービスでは、Python・NumPy・SciPyを使った科学計算コードについて、以下のような確認を行っています。

・χ²計算や尤度関数の確認  

・共分散行列処理の見直し  

・明示的な逆行列計算の回避  

・NumPy / SciPy の高速化  

・MCMCや研究用コードの再現性整理  

・計算結果のログ、README、監査レポート作成  

単に「動くコード」ではなく、「数値的に安全か」「第三者が再現できるか」「結果としてどこまで主張できるか」を整理することを重視しています。

研究用コード、PoC、GitHub公開前の確認、論文用の計算コード整理などでお困りの場合は、お気軽にご相談ください。

Python

NumPy

SciPy

数値解析

科学計算

MCMC

データ解析