連立一次方程式 $\bm{Ax} = \bm{b}$ を解くことは、科学技術計算の最も基本的な操作です。大学の線形代数では「逆行列 $\bm{A}^{-1}$ を求めて $\bm{x} = \bm{A}^{-1}\bm{b}$ とすればよい」と習いますが、実際の数値計算で逆行列を直接計算することはほとんどありません。なぜなら、逆行列の計算は遅く、数値的にも不安定だからです。
では、大規模な連立方程式を効率的に解くにはどうすればよいでしょうか。その答えがLU分解（LU decomposition）です。行列を下三角行列 $\bm{L}$ と上三角行列 $\bm{U}$ の積に分解すれば、三角行列の連立方程式は前進代入・後退代入で $O(n^2)$ で解けます。
LU分解は、手計算で馴染みのあるガウスの消去法を行列の言葉で再解釈したものです。この対応を理解すると、「なぜガウスの消去法がうまくいくのか」「どこで数値的な問題が起きるのか」が明確になります。
LU分解を理解すると、以下のような場面で活用できます。
- 連立一次方程式の高速な求解: 同じ $\bm{A}$ に対して複数の右辺 $\bm{b}_1, \bm{b}_2, \ldots$ を解く場合、LU分解を1回行えば各 $\bm{b}_i$ は $O(n^2)$ で解ける
- 行列式の計算: $\det(\bm{A}) = \det(\bm{L}) \cdot \det(\bm{U}) = \prod u_{ii}$ で高速に計算できる
- 逆行列の計算: 必要な場合でも、LU分解を経由する方が直接計算より安定
- 数値シミュレーション: 有限要素法や差分法で現れる大規模疎行列の求解
本記事の内容
- ガウスの消去法の行列表現
- LU分解の定義と導出
- 部分ピボット選択（PA = LU）
- 前進代入と後退代入
- 計算量の解析
- Pythonでのスクラッチ実装と数値実験
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
- 行列式の定義・計算方法・幾何学的意味 — 行列式と正則性
- 固有値分解の理論と応用 — 行列分解の基本的な考え方
ガウスの消去法を行列で表現する
ガウスの消去法の復習
ガウスの消去法は、連立方程式の拡大係数行列に行基本操作を繰り返し適用して上三角行列（階段行列）に変換する手法です。具体的には、各列について「ピボット行の定数倍を下の行から引く」操作を繰り返します。
3次の例で見てみましょう。
$$ \bm{A} = \begin{pmatrix} 2 & 1 & 1 \\ 4 & 3 & 3 \\ 8 & 7 & 9 \end{pmatrix} $$
ステップ1: 第1列の消去。第2行から第1行の $4/2 = 2$ 倍を引き、第3行から第1行の $8/2 = 4$ 倍を引きます。
$$ \begin{pmatrix} 2 & 1 & 1 \\ 0 & 1 & 1 \\ 0 & 3 & 5 \end{pmatrix} $$
ステップ2: 第2列の消去。第3行から第2行の $3/1 = 3$ 倍を引きます。
$$ \begin{pmatrix} 2 & 1 & 1 \\ 0 & 1 & 1 \\ 0 & 0 & 2 \end{pmatrix} = \bm{U} $$
これで上三角行列が得られました。ここまでは高校数学でも見かける操作ですが、この「行から定数倍を引く」操作を行列の言葉で書き直すことが、LU分解への鍵です。
消去操作の行列表現
「第 $i$ 行から第 $j$ 行の $m_{ij}$ 倍を引く」操作は、消去行列（elimination matrix）$\bm{E}$ を左から掛けることに対応します。
ステップ1の操作は2つの消去行列で表されます。
$$ \bm{E}_1 = \begin{pmatrix} 1 & 0 & 0 \\ -2 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix}, \quad \bm{E}_2 = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ -4 & 0 & 1 \end{pmatrix} $$
$\bm{E}_1$ は「第2行から第1行の2倍を引く」、$\bm{E}_2$ は「第3行から第1行の4倍を引く」操作です。ステップ2では
$$ \bm{E}_3 = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & -3 & 1 \end{pmatrix} $$
が「第3行から第2行の3倍を引く」操作です。
ガウスの消去法全体は次のように書けます。
$$ \bm{E}_3 \bm{E}_2 \bm{E}_1 \bm{A} = \bm{U} $$
したがって
$$ \bm{A} = \bm{E}_1^{-1} \bm{E}_2^{-1} \bm{E}_3^{-1} \bm{U} $$
ここで重要な事実があります。消去行列の逆行列は符号を反転するだけです。
$$ \bm{E}_1^{-1} = \begin{pmatrix} 1 & 0 & 0 \\ 2 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix}, \quad \bm{E}_2^{-1} = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 4 & 0 & 1 \end{pmatrix}, \quad \bm{E}_3^{-1} = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 3 & 1 \end{pmatrix} $$
さらに、これらの逆消去行列の積は、各乗数を対応する位置にそのまま配置した下三角行列になります。
$$ \bm{L} = \bm{E}_1^{-1} \bm{E}_2^{-1} \bm{E}_3^{-1} = \begin{pmatrix} 1 & 0 & 0 \\ 2 & 1 & 0 \\ 4 & 3 & 1 \end{pmatrix} $$
対角成分より下の要素が、まさにガウスの消去法で使った乗数 $m_{21} = 2$, $m_{31} = 4$, $m_{32} = 3$ になっています。これが「ガウスの消去法は LU分解と同じ操作である」という主張の正体です。
ガウスの消去法を行列の言葉で書き直せたところで、LU分解を一般的に定式化しましょう。
LU分解の定義と存在条件
定義
$n \times n$ 正方行列 $\bm{A}$ が
$$ \begin{equation} \bm{A} = \bm{L}\bm{U} \end{equation} $$
と分解できるとき、これをLU分解と呼びます。ここで $\bm{L}$ は対角成分が全て1の下三角行列（単位下三角行列）、$\bm{U}$ は上三角行列です。
$$ \bm{L} = \begin{pmatrix} 1 & & & \\ l_{21} & 1 & & \\ \vdots & \ddots & \ddots & \\ l_{n1} & \cdots & l_{n,n-1} & 1 \end{pmatrix}, \quad \bm{U} = \begin{pmatrix} u_{11} & u_{12} & \cdots & u_{1n} \\ & u_{22} & \cdots & u_{2n} \\ & & \ddots & \vdots \\ & & & u_{nn} \end{pmatrix} $$
$\bm{L}$ の対角成分を1に固定する規約をDoolittle分解と呼び、最も一般的です。$\bm{U}$ の対角成分を1にする場合はCrout分解と呼びます。
存在条件
全ての正則行列がLU分解を持つわけではありません。LU分解が存在するための必要十分条件は、$\bm{A}$ の全ての首座小行列式（leading principal minor）がゼロでないことです。
つまり、$k = 1, 2, \ldots, n-1$ に対して
$$ \det(\bm{A}_{k}) \neq 0, \quad \bm{A}_{k} = \begin{pmatrix} a_{11} & \cdots & a_{1k} \\ \vdots & \ddots & \vdots \\ a_{k1} & \cdots & a_{kk} \end{pmatrix} $$
が成り立つことが条件です。
直感的には、ガウスの消去法の途中でピボット（対角要素）がゼロにならないことを意味します。ピボットがゼロになると消去操作が止まってしまうため、行の入れ替え（ピボット選択）が必要になります。
一意性
LU分解が存在するとき（つまり全ての首座小行列式がゼロでないとき）、$\bm{L}$ の対角成分を1に固定すれば分解は一意です。
この一意性は、$\bm{L}$ と $\bm{U}$ の要素を直接計算するアルゴリズム（後述）の正当性を保証します。
LU分解の存在条件を見ると、ピボットがゼロになる場合の対処が必要です。ここでピボット選択の概念を導入しましょう。
部分ピボット選択（PLU分解）
なぜピボット選択が必要か
ガウスの消去法でピボット $u_{kk}$ がゼロになると消去が不可能になります。ゼロでなくても、ピボットが非常に小さいと深刻な問題が起きます。
例えば、$a_{11} = 10^{-20}$ のとき、乗数 $m_{21} = a_{21}/a_{11}$ が巨大になり、丸め誤差が増幅されます。これは数値的不安定性と呼ばれ、理論的には正しい結果を、有限精度の浮動小数点演算では得られなくなります。
部分ピボット選択の方法
部分ピボット選択（partial pivoting）は、第 $k$ ステップで第 $k$ 列の $k$ 行目以降で最大の絶対値を持つ要素を探し、その行と第 $k$ 行を入れ替える手法です。
行の入れ替えは置換行列 $\bm{P}$ で表されます。$\bm{P}$ は単位行列の行を入れ替えた行列で、以下の性質を持ちます。
- $\bm{P}$ は直交行列: $\bm{P}^T\bm{P} = \bm{I}$
- $\bm{P}^{-1} = \bm{P}^T$
- $\det(\bm{P}) = \pm 1$
部分ピボット選択付きLU分解は
$$ \begin{equation} \bm{PA} = \bm{LU} \end{equation} $$
と表されます。$\bm{P}$ は全てのステップの行入れ替えをまとめた置換行列です。
重要な事実: 任意の正則行列 $\bm{A}$ に対して、適切な置換行列 $\bm{P}$ を選べば $\bm{PA} = \bm{LU}$ は常に存在します。つまり、部分ピボット選択により存在条件の問題は解消されます。
数値安定性の保証
部分ピボット選択では、乗数 $|l_{ij}| \leq 1$ が保証されます。なぜなら、常にその列で最大の要素をピボットとして選ぶため、他の要素をピボットで割った値の絶対値は1以下になるからです。
これにより丸め誤差の蓄積が抑制され、実用的にほぼ全てのケースで数値的に安定な結果が得られます。NumPyやLAPACKなどの標準的な線形代数ライブラリは、LU分解に常に部分ピボット選択を使用しています。
LU分解の理論がわかったところで、分解結果を使って連立方程式を解く手順を見ていきましょう。
前進代入と後退代入
連立方程式への適用
$\bm{Ax} = \bm{b}$ を解くために、LU分解 $\bm{A} = \bm{LU}$（ピボット付きなら $\bm{PA} = \bm{LU}$）を利用します。
$$ \bm{LUx} = \bm{b} \quad (\text{or} \quad \bm{LUx} = \bm{Pb}) $$
$\bm{y} = \bm{Ux}$ とおくと、2段階で解けます。
ステップ1: 前進代入（forward substitution）で $\bm{Ly} = \bm{b}$ を解く
$$ y_i = b_i – \sum_{j=1}^{i-1} l_{ij} y_j, \quad i = 1, 2, \ldots, n $$
$\bm{L}$ は下三角行列なので、$y_1$ から順に確定していきます。
ステップ2: 後退代入（backward substitution）で $\bm{Ux} = \bm{y}$ を解く
$$ x_i = \frac{1}{u_{ii}}\left(y_i – \sum_{j=i+1}^{n} u_{ij} x_j\right), \quad i = n, n-1, \ldots, 1 $$
$\bm{U}$ は上三角行列なので、$x_n$ から逆順に確定していきます。
なぜ2段階に分けるのか
一見すると回りくどいように思えますが、三角行列の方程式は $O(n^2)$ で解ける点が重要です。LU分解自体は $O(n^3)$ ですが、一度分解してしまえば
- 同じ $\bm{A}$ に対する別の右辺 $\bm{b}’$ は $O(n^2)$ で解ける
- 行列式は $\det(\bm{A}) = \pm \prod_{i=1}^n u_{ii}$ で即座に計算できる
右辺が1つだけなら計算量は同じですが、複数の右辺を扱う場合にLU分解は大きな利点を持ちます。
計算の手順が明確になったところで、計算量を定量的に評価しましょう。
計算量の解析
LU分解の計算量
ガウスの消去法（LU分解）で第 $k$ ステップの操作量を数えます。
第 $k$ ステップでは $(n-k)$ 行に対して、乗数の計算（1回の除算）と $(n-k)$ 個の要素の更新（乗算と減算）が必要です。したがって
$$ \text{flops} = \sum_{k=1}^{n-1} (n-k)(n-k) = \sum_{m=1}^{n-1} m^2 \approx \frac{n^3}{3} $$
乗算と加算をまとめて数えると、LU分解の計算量は $\frac{2}{3}n^3$ flops です。
前進代入・後退代入の計算量
前進代入は
$$ \sum_{i=1}^{n} (i-1) = \frac{n(n-1)}{2} \approx \frac{n^2}{2} $$
後退代入も同様に $\frac{n^2}{2}$ flops です。
比較
| 操作 | 計算量 |
|---|---|
| LU分解 | $\frac{2}{3}n^3$ |
| 前進代入 + 後退代入 | $n^2$ |
| 逆行列計算 | $2n^3$ |
| クラメルの公式 | $O(n \cdot n!)$ |
LU分解は逆行列計算の約1/3の計算量であり、複数の右辺を扱う場合はさらに有利です。クラメルの公式は理論的には美しいですが、$n$ が大きくなると階乗的に爆発するため実用的ではありません。
LDU分解
LU分解の変形として、$\bm{U}$ の対角成分を分離したLDU分解があります。
$$ \bm{A} = \bm{L}\bm{D}\bm{U}’ $$
ここで $\bm{D} = \text{diag}(u_{11}, u_{22}, \ldots, u_{nn})$ は対角行列、$\bm{U}’$ は対角成分が1の上三角行列です。LDU分解は対称行列のコレスキー分解との関係で重要になります。
$\bm{A}$ が対称行列（$\bm{A} = \bm{A}^T$）のとき、LDU分解は $\bm{A} = \bm{L}\bm{D}\bm{L}^T$ という美しい形になります。これはコレスキー分解の基礎です。
理論と計算量の解析が済んだところで、Pythonでスクラッチ実装してみましょう。
Pythonでの実装
LU分解のスクラッチ実装
まず、部分ピボット選択付きLU分解を実装します。
import numpy as np
import matplotlib.pyplot as plt
def lu_decomposition_no_pivot(A):
"""ピボット選択なしのLU分解"""
n = A.shape[0]
L = np.eye(n)
U = A.astype(float).copy()
for k in range(n - 1):
if abs(U[k, k]) < 1e-15:
raise ValueError(f"ゼロピボット: u[{k},{k}] = {U[k,k]}")
for i in range(k + 1, n):
# 乗数の計算
L[i, k] = U[i, k] / U[k, k]
# 行の更新
U[i, k:] -= L[i, k] * U[k, k:]
return L, U
def lu_decomposition_partial_pivot(A):
"""部分ピボット選択付きLU分解 (PA = LU)"""
n = A.shape[0]
L = np.zeros((n, n))
U = A.astype(float).copy()
P = np.eye(n)
for k in range(n - 1):
# ピボット選択: k列のk行以降で最大絶対値を探す
pivot_row = k + np.argmax(np.abs(U[k:, k]))
# 行の入れ替え
if pivot_row != k:
U[[k, pivot_row]] = U[[pivot_row, k]]
P[[k, pivot_row]] = P[[pivot_row, k]]
# Lの既に計算済みの部分も入れ替え
L[[k, pivot_row], :k] = L[[pivot_row, k], :k]
if abs(U[k, k]) < 1e-15:
raise ValueError("行列が特異です")
for i in range(k + 1, n):
L[i, k] = U[i, k] / U[k, k]
U[i, k:] -= L[i, k] * U[k, k:]
L += np.eye(n) # 対角に1を追加
return P, L, U
このコードでは2つの関数を実装しています。lu_decomposition_no_pivot はピボットなしの基本版、lu_decomposition_partial_pivot は部分ピボット選択付きの実用版です。ピボット選択では各ステップで最大絶対値の要素をピボットとして選び、行を入れ替えています。
前進代入と後退代入の実装
def forward_substitution(L, b):
"""前進代入: Ly = b を解く"""
n = len(b)
y = np.zeros(n)
for i in range(n):
y[i] = b[i] - L[i, :i] @ y[:i]
return y
def backward_substitution(U, y):
"""後退代入: Ux = y を解く"""
n = len(y)
x = np.zeros(n)
for i in range(n - 1, -1, -1):
x[i] = (y[i] - U[i, i+1:] @ x[i+1:]) / U[i, i]
return x
def solve_lu(A, b):
"""LU分解を使って Ax = b を解く"""
P, L, U = lu_decomposition_partial_pivot(A)
Pb = P @ b
y = forward_substitution(L, Pb)
x = backward_substitution(U, y)
return x
前進代入では $y_1$ から順に、後退代入では $x_n$ から逆順に値を確定していきます。solve_lu 関数はこれらを組み合わせて連立方程式を解きます。
具体例での検証
np.random.seed(42)
# 3x3の具体例
A = np.array([[2, 1, 1],
[4, 3, 3],
[8, 7, 9]], dtype=float)
b = np.array([1, 1, 1], dtype=float)
# LU分解
P, L, U = lu_decomposition_partial_pivot(A)
print("=== LU分解の結果 ===")
print(f"P =\n{P}\n")
print(f"L =\n{L}\n")
print(f"U =\n{U}\n")
print(f"PA - LU の最大誤差: {np.max(np.abs(P @ A - L @ U)):.2e}")
# 連立方程式の求解
x = solve_lu(A, b)
x_numpy = np.linalg.solve(A, b)
print(f"\nスクラッチ実装の解: {x}")
print(f"NumPyの解: {x_numpy}")
print(f"解の差の最大値: {np.max(np.abs(x - x_numpy)):.2e}")
# 行列式
det_lu = np.prod(np.diag(U)) * np.linalg.det(P)
det_numpy = np.linalg.det(A)
print(f"\nLU分解による行列式: {det_lu:.6f}")
print(f"NumPyの行列式: {det_numpy:.6f}")
この検証コードでは、$\bm{PA} – \bm{LU}$ の誤差がゼロ（機械精度のオーダー）であることを確認し、スクラッチ実装の解とNumPyの解が一致することを検証しています。行列式もLU分解の $\bm{U}$ の対角成分の積（符号は $\bm{P}$ の行列式で調整）から計算できます。
数値安定性の可視化
ピボット選択の重要性を実験で確認しましょう。
fig, axes = plt.subplots(1, 3, figsize=(16, 5))
# (a) ピボット選択の効果
ax = axes[0]
sizes = [10, 20, 50, 100, 200]
errors_no_pivot = []
errors_with_pivot = []
for n in sizes:
# ランダム行列で実験（100回平均）
errs_no, errs_with = [], []
for _ in range(100):
A_rand = np.random.randn(n, n)
b_rand = np.random.randn(n)
x_true = np.linalg.solve(A_rand, b_rand)
try:
L_np, U_np = lu_decomposition_no_pivot(A_rand)
y_np = forward_substitution(L_np, b_rand)
x_np = backward_substitution(U_np, y_np)
errs_no.append(np.max(np.abs(x_np - x_true)))
except ValueError:
errs_no.append(np.nan)
x_pp = solve_lu(A_rand, b_rand)
errs_with.append(np.max(np.abs(x_pp - x_true)))
errors_no_pivot.append(np.nanmedian(errs_no))
errors_with_pivot.append(np.median(errs_with))
ax.semilogy(sizes, errors_no_pivot, 'ro-', markersize=8, linewidth=2,
label='Without pivoting')
ax.semilogy(sizes, errors_with_pivot, 'bs-', markersize=8, linewidth=2,
label='With partial pivoting')
ax.set_xlabel('Matrix size $n$', fontsize=12)
ax.set_ylabel('Median error', fontsize=12)
ax.set_title('Numerical Stability: Pivoting Effect', fontsize=13)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3, which='both')
# (b) 計算量の実測
ax = axes[1]
import time
sizes_timing = [50, 100, 200, 500, 1000]
times_lu = []
times_solve = []
for n in sizes_timing:
A_t = np.random.randn(n, n)
b_t = np.random.randn(n)
start = time.time()
for _ in range(3):
P_t, L_t, U_t = lu_decomposition_partial_pivot(A_t)
times_lu.append((time.time() - start) / 3)
start = time.time()
for _ in range(3):
y_t = forward_substitution(L_t, P_t @ b_t)
x_t = backward_substitution(U_t, y_t)
times_solve.append((time.time() - start) / 3)
ax.loglog(sizes_timing, times_lu, 'ro-', markersize=8, linewidth=2,
label='LU decomposition')
ax.loglog(sizes_timing, times_solve, 'bs-', markersize=8, linewidth=2,
label='Forward + Backward subst.')
# 理論曲線 O(n^3) と O(n^2)
n_ref = np.array(sizes_timing, dtype=float)
ax.loglog(n_ref, times_lu[0] * (n_ref/n_ref[0])**3, 'r--', alpha=0.5,
label='$O(n^3)$ reference')
ax.loglog(n_ref, times_solve[0] * (n_ref/n_ref[0])**2, 'b--', alpha=0.5,
label='$O(n^2)$ reference')
ax.set_xlabel('Matrix size $n$', fontsize=12)
ax.set_ylabel('Time (seconds)', fontsize=12)
ax.set_title('Computational Cost', fontsize=13)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3, which='both')
# (c) L, Uの構造可視化
ax = axes[2]
n_vis = 8
A_vis = np.random.randn(n_vis, n_vis)
P_vis, L_vis, U_vis = lu_decomposition_partial_pivot(A_vis)
# LとUを並べて表示
combined = np.zeros((n_vis, 2 * n_vis + 1))
combined[:, :n_vis] = L_vis
combined[:, n_vis+1:] = U_vis
im = ax.imshow(np.abs(combined), cmap='Blues', aspect='equal')
ax.axvline(x=n_vis - 0.5, color='red', linewidth=2)
ax.set_xticks([n_vis//2 - 0.5, n_vis + n_vis//2 + 0.5])
ax.set_xticklabels(['$L$', '$U$'], fontsize=14)
ax.set_title('Structure of $L$ and $U$', fontsize=13)
plt.colorbar(im, ax=ax, shrink=0.8)
plt.tight_layout()
plt.savefig('lu_decomposition.png', dpi=150, bbox_inches='tight')
plt.show()
このグラフから、LU分解の実用的な特性が読み取れます。
-
左図（ピボット選択の効果）: 部分ピボット選択ありの場合（青）は行列サイズが大きくなっても誤差が機械精度付近に留まりますが、ピボット選択なし（赤）は誤差が大きくなる傾向があります。特に大きな行列では、ピボット選択なしでは数値的に不安定な結果が得られる場合があります
-
中央図（計算量）: LU分解の実行時間は $O(n^3)$ の理論曲線に沿い、前進・後退代入は $O(n^2)$ に沿っています。行列サイズが大きくなるほど、分解と求解の計算量の差が顕著になります
-
右図（L, Uの構造）: $\bm{L}$ は下三角構造（左下に非零要素が集中）、$\bm{U}$ は上三角構造（右上に非零要素が集中）を明確に示しています。この疎な構造がメモリ効率と計算効率をもたらします
実用的な注意点
インプレース計算
実装では $\bm{L}$ と $\bm{U}$ を別々の行列として格納しましたが、実際のライブラリ実装では元の行列 $\bm{A}$ を上書きして $\bm{L}$ と $\bm{U}$ を同じ配列に格納します。$\bm{L}$ の対角成分は常に1なので格納不要であり、対角より下に $\bm{L}$ の成分、対角以上に $\bm{U}$ の成分を入れられます。
バンド行列の場合
$\bm{A}$ がバンド幅 $p$（下バンド幅）と $q$（上バンド幅）のバンド行列の場合、LU分解の $\bm{L}$ はバンド幅 $p$、$\bm{U}$ はバンド幅 $q$ を持ちます。この場合、計算量は $O(npq)$ に削減され、有限差分法で現れる三重対角行列（$p = q = 1$）では $O(n)$ で解けます。
SciPy/NumPyでの使い方
from scipy.linalg import lu, lu_factor, lu_solve
# LU分解
P_sp, L_sp, U_sp = lu(A)
print(f"SciPy LU分解: PA - LU error = {np.max(np.abs(P_sp @ A - L_sp @ U_sp)):.2e}")
# 効率的な求解（factorを使い回す）
lu_piv = lu_factor(A)
x1 = lu_solve(lu_piv, b)
# 複数の右辺を効率的に解ける
b2 = np.array([2, 3, 4], dtype=float)
x2 = lu_solve(lu_piv, b2)
print(f"x1 = {x1}")
print(f"x2 = {x2}")
scipy.linalg.lu_factor は分解結果を圧縮形式で保持し、lu_solve で複数の右辺を効率的に解けます。これがLU分解の真の利点 — 「一度分解すれば何度でも $O(n^2)$ で解ける」を活用するパターンです。
コレスキー分解との関係
LU分解と密接に関連する分解として、正定値対称行列に特化したコレスキー分解があります。ここでは、LU分解からコレスキー分解がどのように導かれるかを見ておきましょう。
対称行列 $\bm{A} = \bm{A}^T$ のLDU分解 $\bm{A} = \bm{L}\bm{D}\bm{U}’$ において、対称性から $\bm{U}’ = \bm{L}^T$ が成り立ちます。したがって $\bm{A} = \bm{L}\bm{D}\bm{L}^T$ と書けます。$\bm{A}$ が正定値なら $\bm{D}$ の対角成分は全て正であるため、$\bm{D}^{1/2} = \text{diag}(\sqrt{d_1}, \ldots, \sqrt{d_n})$ が定義でき、$\tilde{\bm{L}} = \bm{L}\bm{D}^{1/2}$ とおけば
$$ \bm{A} = \tilde{\bm{L}} \tilde{\bm{L}}^T $$
これがコレスキー分解 $\bm{A} = \bm{L}\bm{L}^T$（ここでの $\bm{L}$ は対角成分が正の下三角行列）です。計算量は $\frac{1}{3}n^3$ と、LU分解の約半分です。共分散行列や正規方程式の係数行列 $\bm{A}^T\bm{A}$ のように正定値対称行列が現れる場面では、コレスキー分解を使うことで計算コストと数値安定性の両面で有利になります。
まとめ
本記事では、LU分解の理論と実装について解説しました。
- LU分解は行列 $\bm{A}$ を下三角行列 $\bm{L}$ と上三角行列 $\bm{U}$ の積に分解する手法で、ガウスの消去法の行列表現である
- LU分解が存在するには全ての首座小行列式が非零である必要があるが、部分ピボット選択（$\bm{PA} = \bm{LU}$）により任意の正則行列に適用できる
- 連立方程式は前進代入（$\bm{Ly} = \bm{b}$）と後退代入（$\bm{Ux} = \bm{y}$）の2段階で $O(n^2)$ で解ける
- LU分解の計算量は $\frac{2}{3}n^3$ で、同じ $\bm{A}$ に対する複数の右辺を効率的に処理できる
- 部分ピボット選択は数値安定性を保証するために不可欠であり、実用的なライブラリでは常に使われている
次のステップとして、以下の記事も参考にしてください。
- QR分解とグラム・シュミット法 — 直交行列を使った別の分解
- コレスキー分解の理論と実装 — 対称正定値行列の場合のLU分解
- 行列式の定義・計算方法・幾何学的意味 — LU分解による行列式の計算