連立方程式Ax=bを解くときは、逆行列を明示的に作らずsolveを使い、残差と条件数を併せて確認する。残差が小さくても、入力の小さなずれで解が大きく動く問題はある。数値手続きが式を満たしたことと、得られた解が入力誤差に対して安定であることを区別しよう。
最小例で確かめる
以下は説明用に作った小さなデータである。実測データや実行速度の測定結果ではない。コード全体をexample.pyとして保存すれば、入力ファイルを別途用意せずに実行できる。assertは、この例で成り立つべき形や値を確認するために入れてある。
import numpy as np
from scipy.linalg import solve
a = np.array([[1., 1.], [1., 1. + 1e-8]])
b = np.array([2., 2. + 1e-8])
x = solve(a, b)
b_changed = b + np.array([0., 1e-8])
x_changed = solve(a, b_changed)
condition = np.linalg.cond(a)
def scaled_residual(a, x, b):
scale = np.linalg.norm(a, 2) * np.linalg.norm(x, 2) + np.linalg.norm(b, 2)
return np.linalg.norm(a @ x - b, 2) / scale if scale else 0.
r0 = scaled_residual(a, x, b)
r1 = scaled_residual(a, x_changed, b_changed)
print("solution:", x.round(6).tolist())
print("perturbed solution:", x_changed.round(6).tolist())
print("condition number:", f"{condition:.2e}")
print("both scaled residuals < 1e-12:", r0 < 1e-12 and r1 < 1e-12)
print("relative input change:", f"{np.linalg.norm(b_changed-b)/np.linalg.norm(b):.2e}")
print("relative solution change:", f"{np.linalg.norm(x_changed-x)/np.linalg.norm(x):.2e}")
try:
solve([[1.,1.],[1.,1.]], [2.,2.]) # 意図的な特異系
except np.linalg.LinAlgError:
print("expected: LinAlgError")
else:
raise AssertionError("singular system was accepted")
np.testing.assert_allclose(x, [1.,1.], atol=1e-6)
np.testing.assert_allclose(x_changed, [0.,2.], atol=1e-6)
assert condition > 1e8 and r0 < 1e-12 and r1 < 1e-12
実行結果
solution: [1.0, 1.0]
perturbed solution: [0.0, 2.0]
condition number: 4.00e+08
both scaled residuals < 1e-12: True
relative input change: 3.54e-09
relative solution change: 1.00e+00
expected: LinAlgError
式を満たすかは残差で確認する
残差はa @ x – bであり、得られたxを元の式へ代入した差である。ただし係数や右辺が大きければ残差の絶対値も大きくなりやすいため、この例では行列とベクトルの2ノルムを使った尺度で割っている。どの残差を報告したかを明示し、単位や値の大きさが違う問題を無条件に同じ閾値で比べない。
solveは係数行列と右辺から数値解を求める。inv(a) @ bという形で逆行列を作る必要はなく、通常は直接求解の方が目的を明確に表せる。入力のshape、dtype、有限性を確認し、複数の右辺がある場合も、それぞれの列が何に対応するかを決めておこう。
条件数は入力に対する敏感さの目安
この行列は二つの行が非常に近く、二つの方程式から別々の情報を取り出しにくい。np.linalg.condの標準では2ノルムの条件数を計算し、例では約4×10**8となる。大きい条件数は、入力にある方向の小さな変化が、解の大きな変化につながり得ることを示す。
条件数は特定の今回の誤差が必ずその倍率で増えるという予測ではなく、悪い方向の感度を表す指標である。また行列の尺度や変数の単位にも影響される。単位を適切にそろえるなどの前処理が必要な場合はあるが、桁を変えるだけで元の問題が十分な情報を持つようになるとは限らない。
二つの小さな残差でも解は動く
右辺の2番目へ1e-8を加えただけで、解はおよそ[1, 1]から[0, 2]へ動く。どちらもそれぞれの右辺に対する残差は小さい。したがって、小さな残差だけを根拠に「解の各成分が多くの桁まで信頼できる」と結論するのは危険である。
この例で表示する入力変化と解の変化は、作った配列に対する決定的な計算であり、測定器の精度や経験的な不確かさを示したものではない。実測問題なら、係数と右辺の誤差の大きさを見積もり、その範囲で解がどれだけ変わるかを調べる必要がある。
特異、悪条件、モデルの不足を分ける
二つの行が完全に同じ特異系では、このsolve呼び出しはLinAlgErrorになる。一方、悪条件でも必ず例外になるわけではない。ライブラリの警告の有無だけで受け入れず、条件数や残差を検査する。しきい値は求める解の精度と入力の精度に合わせて判断する。
未知数より観測が多い問題や、式を厳密には同時に満たせない問題では、最小二乗など問題の定式化を見直すこともある。正則化を使えば安定させられる場合があるが、追加した仮定やバイアスを明示する必要がある。単に例外を避けるために係数へ小さな数を足す方法を万能策として使わない。
掲載例では解の参照値、残差、条件数、特異系の失敗をまとめて検査した。複雑な係数行列でも、まず小さな既知の系で入力と出力の契約を確認し、実問題では推定した値の有効桁を入力誤差の範囲内で解釈するのがよい。
動作確認環境と参考資料
Linux・CPython 3.12.14、NumPy 2.3.5、pandas 2.2.3、SciPy 1.17.0、Matplotlib 3.10.8の環境で掲載コードを実行した。使用するライブラリはコード冒頭のimportを参照してほしい。公式資料の最新版と、この実行確認版は区別している。数値の末尾や表の表示幅は環境によって変わることがある。