連立一次方程式は、科学技術計算のあらゆる場面で現れる最も基本的な問題です。ガウスの消去法（Gaussian elimination）は、連立一次方程式を体系的に解くための標準的なアルゴリズムであり、数値計算の基礎を成す手法です。
本記事では、ガウスの消去法のアルゴリズムを丁寧に解説し、Pythonでスクラッチ実装するところまでを扱います。
本記事の内容
- 連立一次方程式の行列表現
- ガウスの消去法のアルゴリズム（前進消去と後退代入）
- 部分ピボット選択
- Pythonでのスクラッチ実装
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
連立一次方程式の行列表現
連立一次方程式
$$ \begin{cases} a_{11}x_1 + a_{12}x_2 + \cdots + a_{1n}x_n = b_1 \\ a_{21}x_1 + a_{22}x_2 + \cdots + a_{2n}x_n = b_2 \\ \vdots \\ a_{n1}x_1 + a_{n2}x_2 + \cdots + a_{nn}x_n = b_n \end{cases} $$
は、行列とベクトルを用いて
$$ \bm{A}\bm{x} = \bm{b} $$
と書けます。ここで $\bm{A} = (a_{ij})$ は $n \times n$ 係数行列、$\bm{x} = (x_1, \dots, x_n)^T$、$\bm{b} = (b_1, \dots, b_n)^T$ です。
拡大係数行列は、$\bm{A}$ の右側に $\bm{b}$ を付け加えた行列 $[\bm{A} | \bm{b}]$ です。
$$ [\bm{A} | \bm{b}] = \left(\begin{array}{cccc|c} a_{11} & a_{12} & \cdots & a_{1n} & b_1 \\ a_{21} & a_{22} & \cdots & a_{2n} & b_2 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ a_{n1} & a_{n2} & \cdots & a_{nn} & b_n \end{array}\right) $$
ガウスの消去法
ガウスの消去法は、拡大係数行列に行の基本変形を施して上三角行列に変換し（前進消去）、その後に後退代入で解を求めるアルゴリズムです。
行の基本変形
- 2つの行を入れ替える
- ある行を0でない定数倍する
- ある行に他の行の定数倍を加える
これらの操作は連立方程式の解を変えません。
前進消去（Forward Elimination）
拡大係数行列を上三角形に変換する過程です。
第 $k$ ステップ（$k = 1, 2, \dots, n-1$）では、第 $k$ 行を使って第 $k+1, k+2, \dots, n$ 行の第 $k$ 列の成分を0にします。
第 $i$ 行（$i > k$）に対して、
$$ \text{第 } i \text{ 行} \leftarrow \text{第 } i \text{ 行} – \frac{a_{ik}}{a_{kk}} \times \text{第 } k \text{ 行} $$
ここで $a_{kk}$ をピボット（pivot）と呼びます。
後退代入（Back Substitution）
上三角行列になった後、最後の行から順に変数の値を求めます。
$$ x_n = \frac{b_n’}{a_{nn}’} $$
$$ x_k = \frac{1}{a_{kk}’}\left(b_k’ – \sum_{j=k+1}^{n} a_{kj}’ x_j\right) \quad (k = n-1, n-2, \dots, 1) $$
具体例
次の連立方程式を解きます。
$$ \begin{cases} 2x + y – z = 8 \\ -3x – y + 2z = -11 \\ -2x + y + 2z = -3 \end{cases} $$
拡大係数行列:
$$ \left(\begin{array}{ccc|c} 2 & 1 & -1 & 8 \\ -3 & -1 & 2 & -11 \\ -2 & 1 & 2 & -3 \end{array}\right) $$
ステップ1: 第1列の消去
$R_2 \leftarrow R_2 + \frac{3}{2} R_1$, $R_3 \leftarrow R_3 + R_1$
$$ \left(\begin{array}{ccc|c} 2 & 1 & -1 & 8 \\ 0 & 1/2 & 1/2 & 1 \\ 0 & 2 & 1 & 5 \end{array}\right) $$
ステップ2: 第2列の消去
$R_3 \leftarrow R_3 – 4 R_2$
$$ \left(\begin{array}{ccc|c} 2 & 1 & -1 & 8 \\ 0 & 1/2 & 1/2 & 1 \\ 0 & 0 & -1 & 1 \end{array}\right) $$
後退代入:
$$ \begin{align} -z &= 1 \Rightarrow z = -1 \\ \frac{1}{2}y + \frac{1}{2}(-1) &= 1 \Rightarrow y = 3 \\ 2x + 3 – (-1) &= 8 \Rightarrow x = 2 \end{align} $$
解: $(x, y, z) = (2, 3, -1)$
部分ピボット選択
ピボット $a_{kk}$ が0またはほぼ0のとき、前進消去は失敗するか数値的に不安定になります。これを回避するために部分ピボット選択（partial pivoting）を行います。
第 $k$ ステップで、第 $k$ 列の $a_{kk}, a_{k+1,k}, \dots, a_{nk}$ のうち絶対値最大のものを探し、その行と第 $k$ 行を入れ替えます。
Pythonでの実装
import numpy as np
def gaussian_elimination(A, b):
"""部分ピボット選択付きガウスの消去法"""
n = len(b)
# 拡大係数行列を作成
Ab = np.hstack([A.astype(float), b.reshape(-1, 1).astype(float)])
# 前進消去
for k in range(n - 1):
# 部分ピボット選択
max_row = k + np.argmax(np.abs(Ab[k:, k]))
if max_row != k:
Ab[[k, max_row]] = Ab[[max_row, k]]
if abs(Ab[k, k]) < 1e-12:
raise ValueError("行列が特異です。")
for i in range(k + 1, n):
factor = Ab[i, k] / Ab[k, k]
Ab[i, k:] -= factor * Ab[k, k:]
# 後退代入
x = np.zeros(n)
for k in range(n - 1, -1, -1):
x[k] = (Ab[k, -1] - np.dot(Ab[k, k+1:n], x[k+1:n])) / Ab[k, k]
return x
# テスト
A = np.array([[2, 1, -1],
[-3, -1, 2],
[-2, 1, 2]], dtype=float)
b = np.array([8, -11, -3], dtype=float)
x_custom = gaussian_elimination(A, b)
x_numpy = np.linalg.solve(A, b)
print("自作ガウス消去法:", x_custom)
print("NumPy solve: ", x_numpy)
print("一致確認:", np.allclose(x_custom, x_numpy))
消去過程の可視化
import numpy as np
import matplotlib.pyplot as plt
def gaussian_elimination_steps(A, b):
"""消去の各ステップを記録するバージョン"""
n = len(b)
Ab = np.hstack([A.astype(float), b.reshape(-1, 1).astype(float)])
steps = [Ab.copy()]
for k in range(n - 1):
max_row = k + np.argmax(np.abs(Ab[k:, k]))
if max_row != k:
Ab[[k, max_row]] = Ab[[max_row, k]]
for i in range(k + 1, n):
factor = Ab[i, k] / Ab[k, k]
Ab[i, k:] -= factor * Ab[k, k:]
steps.append(Ab.copy())
return steps
A = np.array([[2, 1, -1],
[-3, -1, 2],
[-2, 1, 2]], dtype=float)
b = np.array([8, -11, -3], dtype=float)
steps = gaussian_elimination_steps(A, b)
fig, axes = plt.subplots(1, len(steps), figsize=(4 * len(steps), 3))
for idx, (step, ax) in enumerate(zip(steps, axes)):
im = ax.imshow(step, cmap='RdBu_r', vmin=-12, vmax=12, aspect='auto')
for i in range(step.shape[0]):
for j in range(step.shape[1]):
ax.text(j, i, f'{step[i,j]:.1f}', ha='center', va='center', fontsize=9)
ax.set_title(f'Step {idx}')
ax.set_xticks([])
ax.set_yticks([])
plt.suptitle('Gaussian Elimination Steps', fontsize=14)
plt.tight_layout()
plt.savefig('gaussian_elimination_steps.png', dpi=150, bbox_inches='tight')
plt.show()
計算量の確認
import numpy as np
import time
import matplotlib.pyplot as plt
sizes = [10, 50, 100, 200, 500]
times_custom = []
times_numpy = []
for n in sizes:
A = np.random.randn(n, n)
b = np.random.randn(n)
start = time.time()
for _ in range(3):
gaussian_elimination(A, b)
times_custom.append((time.time() - start) / 3)
start = time.time()
for _ in range(3):
np.linalg.solve(A, b)
times_numpy.append((time.time() - start) / 3)
fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(sizes, times_custom, 'bo-', label='Custom Gaussian Elimination')
ax.plot(sizes, times_numpy, 'r^-', label='NumPy linalg.solve')
ax.set_xlabel('Matrix size $n$')
ax.set_ylabel('Time (s)')
ax.set_title('Gaussian Elimination: Computation Time')
ax.legend()
ax.grid(True, alpha=0.3)
ax.set_yscale('log')
plt.tight_layout()
plt.savefig('gaussian_elimination_time.png', dpi=150, bbox_inches='tight')
plt.show()
ガウスの消去法の計算量は $O(n^3)$ です。NumPyの linalg.solve はLAPACKのLU分解を内部で使用しており、同じ $O(n^3)$ ですが最適化されているため高速です。
まとめ
本記事では、連立一次方程式とガウスの消去法について解説しました。
- ガウスの消去法: 前進消去で上三角行列に変換し、後退代入で解を求める
- 部分ピボット選択: 数値的安定性のために絶対値最大のピボットを選ぶ
- 計算量: $O(n^3)$
- Pythonでのスクラッチ実装により、アルゴリズムの仕組みを確認できる
次のステップとして、以下の記事も参考にしてください。