コンピュータで $\bm{Ax} = \bm{b}$ を解くとき、得られた解はどの程度信頼できるでしょうか。数学的には正確に解けるはずの問題でも、浮動小数点演算の丸め誤差により、結果が大きくずれることがあります。
この「解の信頼性」を左右するのが行列の条件数（condition number）です。条件数が小さい行列は「良条件」（well-conditioned）で、入力の微小な摂動に対して解も微小にしか変化しません。逆に条件数が大きい行列は「悪条件」（ill-conditioned）で、わずかな丸め誤差が解に大きく増幅されます。
例えば、ヒルベルト行列 $H_{ij} = \frac{1}{i+j-1}$ は悪条件行列の典型例です。$10 \times 10$ のヒルベルト行列の条件数は約 $10^{13}$ にも達し、倍精度浮動小数点（約16桁の有効桁数）で解くと有効桁数が3桁しか残りません。
条件数を理解すると、以下のような場面で活用できます。
- 数値計算の信頼性評価: 連立方程式の解の精度を事前に見積もる
- アルゴリズム選択: 正規方程式 vs QR分解 vs SVDの使い分け
- 前処理（プリコンディショニング）: 悪条件行列の改善
- 実験データの解析: 計測誤差が結果に与える影響の定量化
本記事の内容
- 条件数の定義と直感的理解
- 連立方程式の摂動解析
- 条件数と固有値・特異値の関係
- 悪条件行列の例と原因
- 対処法（前処理、正則化、高精度演算）
- Pythonでの実装と数値実験
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
- 特異値分解(SVD)の導出と応用 — 特異値の概念
- LU分解の理論と実装 — 連立方程式の数値解法
条件数とは何か
直感的な理解
高速道路の料金所を想像してください。入口で1円の誤差（丸め誤差）が生じたとき、出口での精算額に何円の誤差が出るでしょうか。単純な足し算なら1円のまま（条件数1）ですが、複雑な変換を経ると誤差が増幅されることがあります。
連立方程式 $\bm{Ax} = \bm{b}$ でも同じです。$\bm{b}$ に微小な摂動 $\delta\bm{b}$ が加わったとき、解の摂動 $\delta\bm{x}$ はどれだけ大きくなるでしょうか。この「誤差の増幅率の最大値」が条件数です。
数学的定義
正則行列 $\bm{A}$ の条件数（condition number）は、行列ノルム $\|\cdot\|$ に対して
$$ \begin{equation} \kappa(\bm{A}) = \|\bm{A}\| \cdot \|\bm{A}^{-1}\| \end{equation} $$
で定義されます。$p$-ノルムに対する条件数は $\kappa_p(\bm{A})$ と書きます。
最も一般的に使われるのは2-ノルム条件数です。
$$ \begin{equation} \kappa_2(\bm{A}) = \|\bm{A}\|_2 \cdot \|\bm{A}^{-1}\|_2 = \frac{\sigma_{\max}}{\sigma_{\min}} \end{equation} $$
ここで $\sigma_{\max}$ と $\sigma_{\min}$ は $\bm{A}$ の最大特異値と最小特異値です。条件数は最大特異値と最小特異値の比という非常に直感的な表現を持ちます。
条件数の基本性質
- $\kappa(\bm{A}) \geq 1$（等号は $\bm{A}$ が直交行列のスカラー倍のとき）
- $\kappa(c\bm{A}) = \kappa(\bm{A})$（スカラー倍で変わらない）
- $\kappa(\bm{A}^T) = \kappa(\bm{A})$
- $\kappa(\bm{I}) = 1$
- $\kappa(\bm{A})$ が大きいほど悪条件
$\kappa(\bm{A}) = 1$ の行列は最も良条件であり、直交行列（回転・鏡映）がこれに該当します。直交行列による変換は誤差を全く増幅しないため、数値計算で重宝されます。QR分解やSVDが数値的に安定なのは、直交行列を積極的に活用しているためです。
条件数の定義を踏まえて、連立方程式の解がどの程度敏感かを定量的に評価しましょう。
摂動解析
右辺の摂動
$\bm{Ax} = \bm{b}$ において、$\bm{b}$ が $\bm{b} + \delta\bm{b}$ に摂動されたとき、解は $\bm{x} + \delta\bm{x}$ に変化します。
$$ \bm{A}(\bm{x} + \delta\bm{x}) = \bm{b} + \delta\bm{b} $$
$\bm{A}\delta\bm{x} = \delta\bm{b}$ より $\delta\bm{x} = \bm{A}^{-1}\delta\bm{b}$。ノルムを取ると
$$ \|\delta\bm{x}\| \leq \|\bm{A}^{-1}\| \|\delta\bm{b}\| $$
一方、$\bm{Ax} = \bm{b}$ から $\|\bm{b}\| \leq \|\bm{A}\| \|\bm{x}\|$ なので
$$ \frac{1}{\|\bm{x}\|} \leq \frac{\|\bm{A}\|}{\|\bm{b}\|} $$
2つの不等式を組み合わせると
$$ \begin{equation} \frac{\|\delta\bm{x}\|}{\|\bm{x}\|} \leq \kappa(\bm{A}) \frac{\|\delta\bm{b}\|}{\|\bm{b}\|} \end{equation} $$
これが摂動の基本不等式です。$\bm{b}$ の相対誤差が条件数倍に増幅されて $\bm{x}$ の相対誤差になる、という意味です。
係数行列の摂動
$\bm{A}$ 自体が摂動された場合（$(\bm{A} + \delta\bm{A})(\bm{x} + \delta\bm{x}) = \bm{b}$）、同様の解析から
$$ \frac{\|\delta\bm{x}\|}{\|\bm{x} + \delta\bm{x}\|} \leq \kappa(\bm{A}) \frac{\|\delta\bm{A}\|}{\|\bm{A}\|} $$
が得られます。
有効桁数の見積もり
浮動小数点演算の丸め誤差は $\bm{b}$ や $\bm{A}$ に約 $\epsilon_{\text{mach}} \approx 10^{-16}$（倍精度）の相対誤差を与えます。したがって、解の相対誤差は
$$ \frac{\|\delta\bm{x}\|}{\|\bm{x}\|} \approx \kappa(\bm{A}) \cdot \epsilon_{\text{mach}} $$
失われる有効桁数は
$$ \text{失われる桁数} \approx \log_{10}(\kappa(\bm{A})) $$
例えば $\kappa(\bm{A}) = 10^8$ なら、倍精度の16桁のうち8桁が失われ、有効桁数は約8桁です。$\kappa(\bm{A}) = 10^{16}$ なら全ての桁が失われ、結果は全く信頼できません。
摂動解析の結果を実際に確認しましょう。次のセクションでは、条件数と特異値の関係を深掘りします。
条件数と特異値
SVDによる理解
$\bm{A} = \bm{U\Sigma V}^T$ とSVD分解すると、$\|\bm{A}\|_2 = \sigma_1$（最大特異値）、$\|\bm{A}^{-1}\|_2 = 1/\sigma_n$（最小特異値の逆数）なので
$$ \kappa_2(\bm{A}) = \frac{\sigma_1}{\sigma_n} $$
特異値が全て等しい（$\sigma_1 = \sigma_2 = \cdots = \sigma_n$）なら $\kappa = 1$ で最良条件。特異値の「スプレッド」が大きいほど悪条件です。
対称行列の場合
対称行列 $\bm{A} = \bm{A}^T$ では特異値 = 固有値の絶対値なので
$$ \kappa_2(\bm{A}) = \frac{|\lambda_{\max}|}{|\lambda_{\min}|} $$
固有値の範囲が広いほど条件が悪くなります。
行列積の条件数
$$ \kappa(\bm{AB}) \leq \kappa(\bm{A}) \cdot \kappa(\bm{B}) $$
行列を掛けるごとに条件が悪化する可能性があります。これが、正規方程式 $\bm{A}^T\bm{A}\hat{\bm{x}} = \bm{A}^T\bm{b}$ の条件数が $\kappa(\bm{A}^T\bm{A}) = \kappa(\bm{A})^2$ になる理由です。
悪条件行列の例
ヒルベルト行列
$$ H_{ij} = \frac{1}{i + j – 1}, \quad H_4 = \begin{pmatrix} 1 & 1/2 & 1/3 & 1/4 \\ 1/2 & 1/3 & 1/4 & 1/5 \\ 1/3 & 1/4 & 1/5 & 1/6 \\ 1/4 & 1/5 & 1/6 & 1/7 \end{pmatrix} $$
ヒルベルト行列は正定値対称行列ですが、条件数がサイズとともに指数的に増加します。$n = 10$ で $\kappa \approx 10^{13}$、$n = 15$ で $\kappa \approx 10^{17}$（倍精度の限界を超える）です。
行が非常に似ているため（$\frac{1}{i+j-1}$ の値が似通う）、行列のランクが「数値的にほぼ不足」している状態です。
ヴァンデルモンド行列
$$ V_{ij} = x_i^{j-1}, \quad V = \begin{pmatrix} 1 & x_1 & x_1^2 & \cdots \\ 1 & x_2 & x_2^2 & \cdots \\ \vdots & & & \ddots \end{pmatrix} $$
多項式補間で現れる行列で、ノード $x_i$ が密に配置されると急速に悪条件になります。
スケーリングの不均一な行列
$$ \bm{A} = \begin{pmatrix} 10^6 & 1 \\ 1 & 10^{-6} \end{pmatrix} $$
異なるスケールの成分が混在すると条件が悪化します。特徴量のスケーリング（標準化）が機械学習で重要な理由の一つです。
悪条件行列の原因を理解したところで、Pythonで数値実験をして理論を確認しましょう。
Pythonでの実装と数値実験
条件数と解の精度
import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import hilbert
np.random.seed(42)
fig, axes = plt.subplots(1, 3, figsize=(16, 5))
# (a) ヒルベルト行列の条件数
ax = axes[0]
sizes = range(2, 16)
cond_numbers = []
for n in sizes:
H = hilbert(n)
cond_numbers.append(np.linalg.cond(H))
ax.semilogy(list(sizes), cond_numbers, 'ro-', markersize=8, linewidth=2)
ax.axhline(y=1/np.finfo(float).eps, color='blue', linewidth=1.5, linestyle='--',
label='$1/\\epsilon_{mach}$')
ax.set_xlabel('Matrix size $n$', fontsize=12)
ax.set_ylabel('$\\kappa_2(H_n)$', fontsize=12)
ax.set_title('Condition Number of Hilbert Matrix', fontsize=13)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3, which='both')
# (b) 条件数 vs 解の誤差
ax = axes[1]
sizes_test = [3, 5, 7, 9, 11, 13]
actual_errors = []
predicted_errors = []
conds = []
for n in sizes_test:
H = hilbert(n)
x_true = np.ones(n)
b = H @ x_true
# 数値的に解く
x_computed = np.linalg.solve(H, b)
rel_error = np.linalg.norm(x_computed - x_true) / np.linalg.norm(x_true)
actual_errors.append(rel_error)
cond = np.linalg.cond(H)
conds.append(cond)
predicted_errors.append(cond * np.finfo(float).eps)
ax.loglog(conds, actual_errors, 'ro', markersize=10, label='Actual error')
ax.loglog(conds, predicted_errors, 'b--', linewidth=2, label='$\\kappa \\cdot \\epsilon_{mach}$')
# 1:1 reference line
cond_range = np.logspace(1, 18, 100)
ax.loglog(cond_range, cond_range * np.finfo(float).eps, 'b--', alpha=0.3)
ax.set_xlabel('$\\kappa_2(A)$', fontsize=12)
ax.set_ylabel('Relative error $\\|\\delta x\\| / \\|x\\|$', fontsize=12)
ax.set_title('Condition Number vs Solution Error', fontsize=13)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3, which='both')
# (c) 特異値分布
ax = axes[2]
for n, color in [(5, 'blue'), (8, 'green'), (12, 'red')]:
H = hilbert(n)
sv = np.linalg.svd(H, compute_uv=False)
ax.semilogy(range(1, n+1), sv, f'{color}o-', markersize=6, linewidth=1.5,
label=f'$n={n}$, $\\kappa={np.linalg.cond(H):.0e}$')
ax.set_xlabel('Singular value index', fontsize=12)
ax.set_ylabel('$\\sigma_i$', fontsize=12)
ax.set_title('Singular Values of Hilbert Matrix', fontsize=13)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3, which='both')
plt.tight_layout()
plt.savefig('condition_number.png', dpi=150, bbox_inches='tight')
plt.show()
このグラフから、条件数と数値精度の関係が明確に確認できます。
-
左図（ヒルベルト行列の条件数）: サイズ $n$ が増えるにつれて条件数が指数的に増大しています。$n \geq 13$ では条件数が $1/\epsilon_{\text{mach}} \approx 10^{16}$（青破線）を超え、倍精度では全く信頼できない領域に入ります
-
中央図（条件数 vs 誤差）: 実際の解の相対誤差（赤丸）が理論的な上界 $\kappa \cdot \epsilon_{\text{mach}}$（青破線）にほぼ沿っています。条件数が $10^{13}$ 程度で有効桁数は約3桁、$10^{16}$ を超えると誤差が1を超え（有効桁数0）、解が全く意味を成さなくなります
-
右図（特異値分布）: ヒルベルト行列の特異値が急激に減衰しています。$n = 12$ では最大特異値と最小特異値の差が15桁以上あり、これが巨大な条件数の原因です。特異値の「崖」が数値的なランク不足を引き起こしています
摂動解析の可視化
fig, axes = plt.subplots(1, 2, figsize=(14, 6))
# (a) 良条件 vs 悪条件の摂動応答
n_pert = 100
# 良条件: ランダム直交行列 + 穏やかなスケーリング
Q, _ = np.linalg.qr(np.random.randn(3, 3))
D_good = np.diag([3, 2, 1])
A_good = Q @ D_good @ Q.T # cond ~ 3
# 悪条件: 大きなスケーリング差
D_bad = np.diag([1000, 1, 0.001])
A_bad = Q @ D_bad @ Q.T # cond ~ 10^6
for idx, (A, title, ax) in enumerate([
(A_good, f'Well-conditioned ($\\kappa = {np.linalg.cond(A_good):.1f}$)', axes[0]),
(A_bad, f'Ill-conditioned ($\\kappa = {np.linalg.cond(A_bad):.1e}$)', axes[1])
]):
x_true = np.array([1, 1, 1], dtype=float)
b_true = A @ x_true
relative_perturbations = []
relative_errors = []
for _ in range(n_pert):
# bをランダムに摂動
delta_b = np.random.randn(3)
delta_b = delta_b / np.linalg.norm(delta_b) # 方向はランダム
eps = 1e-6 * np.linalg.norm(b_true)
delta_b *= eps
b_pert = b_true + delta_b
x_pert = np.linalg.solve(A, b_pert)
rel_pert_b = np.linalg.norm(delta_b) / np.linalg.norm(b_true)
rel_error_x = np.linalg.norm(x_pert - x_true) / np.linalg.norm(x_true)
relative_perturbations.append(rel_pert_b)
relative_errors.append(rel_error_x)
ax.scatter(relative_perturbations, relative_errors, alpha=0.5, s=20,
color='steelblue')
# 理論上界
cond = np.linalg.cond(A)
max_pert = max(relative_perturbations)
ax.plot([0, max_pert], [0, cond * max_pert], 'r--', linewidth=2,
label=f'Upper bound: $\\kappa \\cdot \\delta b/b$')
ax.plot([0, max_pert], [0, max_pert], 'g--', linewidth=1.5,
alpha=0.5, label='1:1 line')
ax.set_xlabel('$\\|\\delta b\\| / \\|b\\|$', fontsize=12)
ax.set_ylabel('$\\|\\delta x\\| / \\|x\\|$', fontsize=12)
ax.set_title(title, fontsize=12)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('perturbation_analysis.png', dpi=150, bbox_inches='tight')
plt.show()
このグラフから、条件数による誤差増幅の違いが直感的にわかります。
-
左図（良条件）: $\bm{b}$ の相対摂動と $\bm{x}$ の相対誤差がほぼ1:1の関係にあります。条件数が3程度なので、入力の誤差が最大でも3倍にしか増幅されず、散布点が緑の1:1ラインの近くに集中しています
-
右図（悪条件）: 同じ大きさの摂動に対して、$\bm{x}$ の相対誤差が最大で $\kappa \approx 10^6$ 倍に増幅されています。赤い上界線と散布点の位置関係から、摂動の方向によって増幅率が大きく異なることもわかります。最悪のケースでは理論上界に近い増幅が起きています
異なるアルゴリズムの安定性比較
fig, ax = plt.subplots(1, 1, figsize=(10, 6))
sizes = [5, 8, 10, 12, 14]
errors_normal = []
errors_qr = []
errors_svd = []
conds = []
for n in sizes:
H = hilbert(n)
x_true = np.ones(n)
b = H @ x_true
conds.append(np.linalg.cond(H))
# 方法1: 正規方程式 (H^T H x = H^T b)
try:
x_normal = np.linalg.solve(H.T @ H, H.T @ b)
errors_normal.append(np.linalg.norm(x_normal - x_true) / np.linalg.norm(x_true))
except np.linalg.LinAlgError:
errors_normal.append(np.nan)
# 方法2: QR分解
Q, R = np.linalg.qr(H)
x_qr = np.linalg.solve(R, Q.T @ b)
errors_qr.append(np.linalg.norm(x_qr - x_true) / np.linalg.norm(x_true))
# 方法3: SVD（最も安定）
x_svd = np.linalg.lstsq(H, b, rcond=None)[0]
errors_svd.append(np.linalg.norm(x_svd - x_true) / np.linalg.norm(x_true))
ax.loglog(conds, errors_normal, 'ro-', markersize=10, linewidth=2,
label='Normal equations ($\\kappa^2$)')
ax.loglog(conds, errors_qr, 'bs-', markersize=10, linewidth=2,
label='QR decomposition ($\\kappa$)')
ax.loglog(conds, errors_svd, 'g^-', markersize=10, linewidth=2,
label='SVD ($\\kappa$)')
# 理論参考線
cond_ref = np.logspace(1, 18, 100)
ax.loglog(cond_ref, cond_ref**2 * np.finfo(float).eps, 'r--', alpha=0.3,
label='$\\kappa^2 \\epsilon_{mach}$')
ax.loglog(cond_ref, cond_ref * np.finfo(float).eps, 'b--', alpha=0.3,
label='$\\kappa \\epsilon_{mach}$')
ax.set_xlabel('$\\kappa_2(A)$', fontsize=12)
ax.set_ylabel('Relative error', fontsize=12)
ax.set_title('Stability Comparison: Normal Eq. vs QR vs SVD', fontsize=13)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3, which='both')
plt.tight_layout()
plt.savefig('stability_comparison.png', dpi=150, bbox_inches='tight')
plt.show()
このグラフから、アルゴリズム選択が数値精度に与える影響が明確に確認できます。
-
正規方程式（赤）: 条件数の二乗に比例して誤差が増大しています。$\kappa(\bm{A}^T\bm{A}) = \kappa(\bm{A})^2$ なので、条件が二重に悪化します。$\kappa \geq 10^8$ 程度で実用的な精度を失います
-
QR分解（青）: 条件数に線形に比例する誤差で、正規方程式より大幅に安定です。$\kappa$ が $10^{12}$ 程度まで実用的な精度を維持しています
-
SVD（緑）: QR分解と同程度の安定性を持ち、さらにランク不足の場合にも対応できる最も頑健な方法です
このように、同じ問題を解くアルゴリズムでも数値安定性が大きく異なります。悪条件な問題では QR分解やSVDを使うべきであり、正規方程式は避けるべきです。
正則化による条件数の改善
ティホノフ正則化の効果をPythonで確認しましょう。
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# (a) 正則化パラメータ vs 条件数
ax = axes[0]
n_reg = 10
H = hilbert(n_reg)
lambdas = np.logspace(-15, 0, 100)
cond_regularized = []
for lam in lambdas:
HtH_reg = H.T @ H + lam * np.eye(n_reg)
cond_regularized.append(np.linalg.cond(HtH_reg))
ax.loglog(lambdas, cond_regularized, 'b-', linewidth=2)
ax.axhline(y=np.linalg.cond(H.T @ H), color='red', linewidth=1.5, linestyle='--',
label=f'$\\kappa(H^TH)$ = {np.linalg.cond(H.T @ H):.1e}')
ax.axhline(y=1/np.finfo(float).eps, color='green', linewidth=1, linestyle=':',
label='$1/\\epsilon_{mach}$')
ax.set_xlabel('$\\lambda$', fontsize=12)
ax.set_ylabel('$\\kappa(H^TH + \\lambda I)$', fontsize=12)
ax.set_title('Regularization Effect on Condition Number', fontsize=13)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3, which='both')
# (b) バイアスと条件数のトレードオフ
ax = axes[1]
x_true = np.ones(n_reg)
b = H @ x_true
lambdas_test = np.logspace(-12, -1, 50)
biases = []
cond_vals = []
for lam in lambdas_test:
x_reg = np.linalg.solve(H.T @ H + lam * np.eye(n_reg), H.T @ b)
biases.append(np.linalg.norm(x_reg - x_true) / np.linalg.norm(x_true))
cond_vals.append(np.linalg.cond(H.T @ H + lam * np.eye(n_reg)))
ax.loglog(lambdas_test, biases, 'r-', linewidth=2, label='Relative error')
ax_twin = ax.twinx()
ax_twin.loglog(lambdas_test, cond_vals, 'b--', linewidth=2, label='$\\kappa$')
ax_twin.set_ylabel('$\\kappa(H^TH + \\lambda I)$', fontsize=12, color='blue')
ax.set_xlabel('$\\lambda$', fontsize=12)
ax.set_ylabel('Relative error', fontsize=12, color='red')
ax.set_title('Bias-Variance Tradeoff in Regularization', fontsize=13)
ax.grid(True, alpha=0.3, which='both')
lines1, labels1 = ax.get_legend_handles_labels()
lines2, labels2 = ax_twin.get_legend_handles_labels()
ax.legend(lines1 + lines2, labels1 + labels2, fontsize=9, loc='center right')
plt.tight_layout()
plt.savefig('regularization_condition.png', dpi=150, bbox_inches='tight')
plt.show()
このグラフから、正則化と条件数の関係が確認できます。
-
左図（正則化の効果）: 正則化パラメータ $\lambda$ を増やすと条件数が劇的に改善されます。$\lambda = 0$ では $\kappa \approx 10^{26}$（赤破線）ですが、$\lambda = 10^{-6}$ で $\kappa \approx 10^{12}$ 程度まで下がり、$\lambda = 1$ ではほぼ1に近づきます
-
右図（バイアスとのトレードオフ）: $\lambda$ が小さいと条件数が大きく数値誤差が支配的になり、$\lambda$ が大きいとバイアス（正則化による真の解からの乖離）が支配的になります。最適な $\lambda$ は両者のバランスが取れる領域にあり、交差検証で選択するのが一般的です
対処法
前処理（プリコンディショニング）
悪条件行列 $\bm{A}$ に対して、条件数を改善する前処理行列 $\bm{M}$ を使って
$$ \bm{M}^{-1}\bm{Ax} = \bm{M}^{-1}\bm{b} $$
を解きます。$\bm{M}$ が $\bm{A}$ の良い近似なら $\kappa(\bm{M}^{-1}\bm{A}) \ll \kappa(\bm{A})$ となり、収束が改善されます。
前処理行列の選び方は問題の構造に依存します。代表的な手法としては、不完全LU分解（ILU）やヤコビ前処理（対角要素の逆数を使う最も単純な方法）があります。理想的には $\bm{M} \approx \bm{A}$ かつ $\bm{M}^{-1}$ の計算コストが低いことが求められます。反復法（共役勾配法やGMRESなど）では、前処理の質が収束速度を支配するため、適切な前処理行列の設計が実用上きわめて重要です。
スケーリング
行や列のスケールが大きく異なる場合、対角スケーリング
$$ \bm{D}_R^{-1}\bm{A}\bm{D}_C^{-1}(\bm{D}_C\bm{x}) = \bm{D}_R^{-1}\bm{b} $$
で条件を改善できます。機械学習での特徴量の標準化はこの一形態です。
正則化
ランク不足に近い行列に対しては、ティホノフ正則化（リッジ回帰）
$$ (\bm{A}^T\bm{A} + \lambda\bm{I})\bm{x} = \bm{A}^T\bm{b} $$
で条件数を $\kappa \approx \frac{\sigma_1^2 + \lambda}{\sigma_n^2 + \lambda}$ に改善できます。$\lambda > 0$ により最小特異値が底上げされます。
高精度演算
倍精度（64ビット）で不十分な場合、四倍精度（128ビット）や任意精度演算を使う方法があります。Pythonではmpmathライブラリが任意精度演算を提供しています。ただし、高精度演算は計算コストが大幅に増加するため、まずスケーリングや前処理で条件数を改善し、それでも不十分な場合の最終手段として用いるのが実践的です。
条件数の応用: 固有値問題と最小二乗法
固有値問題の条件数
連立方程式だけでなく、固有値問題にも条件数の概念があります。行列 $\bm{A}$ の固有値 $\lambda$ に対する摂動 $\delta\lambda$ は、対応する左右の固有ベクトル $\bm{y}, \bm{x}$ を使って
$$ |\delta\lambda| \leq \frac{\|\delta\bm{A}\|}{|\bm{y}^H\bm{x}|} $$
と評価されます。$s(\lambda) = 1/|\bm{y}^H\bm{x}|$ を固有値の条件数と呼びます。対称行列では $\bm{y} = \bm{x}$（直交固有ベクトル）なので $s(\lambda) = 1$ であり、全ての固有値が良条件です。一方、非対称行列では固有ベクトルが「ほぼ平行」になることがあり、固有値が非常に敏感になります。
最小二乗問題の条件数
最小二乗問題 $\min\|\bm{b} – \bm{Ax}\|$ の解の感度は、$\bm{A}$ の条件数だけでなく残差の大きさにも依存します。具体的には
$$ \frac{\|\delta\hat{\bm{x}}\|}{\|\hat{\bm{x}}\|} \leq \kappa(\bm{A}) \left(\frac{\|\delta\bm{A}\|}{\|\bm{A}\|} + \tan\theta \cdot \kappa(\bm{A}) \frac{\|\delta\bm{A}\|}{\|\bm{A}\|}\right) $$
ここで $\theta$ は $\bm{b}$ と列空間 $C(\bm{A})$ のなす角度です。残差が小さい（$\theta \approx 0$、つまりモデルがデータをよく説明している）場合は $\kappa$ に線形に比例しますが、残差が大きい場合は $\kappa^2$ に比例し、正規方程式と同程度に敏感になります。
このことは、回帰モデルの当てはまりが悪い場合（残差が大きい場合）、パラメータ推定が特に不安定になることを意味しています。モデルの改善（説明変数の追加や非線形項の導入）が精度向上に直結する数学的根拠がここにあります。
反復精度改良法
悪条件な連立方程式に対する実用的な手法として反復精度改良法（iterative refinement）があります。LU分解で得た近似解 $\bm{x}_0$ に対して
- 残差を計算: $\bm{r}_0 = \bm{b} – \bm{A}\bm{x}_0$（高精度で計算）
- 補正量を求める: $\bm{A}\bm{d}_0 = \bm{r}_0$（同じLU分解を再利用）
- 解を更新: $\bm{x}_1 = \bm{x}_0 + \bm{d}_0$
このステップを繰り返すと、条件数が大きくても解が改善される場合があります。残差の計算を倍精度より高い精度（拡張精度）で行うことがポイントです。
反復精度改良法は、追加の計算コストが比較的小さいにもかかわらず大きな精度改善が期待できる手法です。LU分解は最初に一度だけ行えばよく、反復ごとに必要なのは行列-ベクトル積と前進・後退代入のみです。そのため、LU分解のコスト $O(n^3)$ に対して各反復は $O(n^2)$ で済みます。実用的には2〜3回の反復で十分な精度が得られることが多く、LAPACK の dgesv 系ルーチンでもオプションとして反復精度改良が実装されています。
以上のように、条件数の理論は連立方程式の直接解法から固有値問題、最小二乗法、さらには反復法の収束解析に至るまで、数値線形代数のあらゆる場面で基盤となる概念です。条件数を把握しておくことで、計算結果を盲目的に信頼するのではなく、その信頼性を定量的に評価し、必要に応じて適切な対処法を選択できるようになります。
まとめ
本記事では、条件数と数値的安定性の理論について解説しました。
- 条件数 $\kappa(\bm{A}) = \|\bm{A}\| \cdot \|\bm{A}^{-1}\| = \sigma_{\max}/\sigma_{\min}$ は入力誤差の最大増幅率
- 摂動解析: $\|\delta\bm{x}\|/\|\bm{x}\| \leq \kappa(\bm{A}) \cdot \|\delta\bm{b}\|/\|\bm{b}\|$ により、失われる有効桁数は $\log_{10}\kappa$ と見積もれる
- ヒルベルト行列は悪条件行列の典型例で、条件数がサイズとともに指数的に増大する
- 正規方程式は条件数が二乗になるため、悪条件な問題ではQR分解やSVDを使うべき
- 前処理、スケーリング、正則化、高精度演算が悪条件への対処法
次のステップとして、以下の記事も参考にしてください。
- LU分解の理論と実装 — ピボット選択と安定性
- QR分解とグラム・シュミット法 — 安定な求解法
- 特異値分解(SVD)の導出と応用 — 最も安定な手法
- 正則化の理論 — 悪条件への対処