Pythonで学ぶ数値計算のアルゴリズムと実装 第6回：連立一次方程式の直截解法：ガウス消去法とピボット選択

Python
2026/01/01
この記事は約18分で読めます。

こんにちは、JS2IIUです。
数値計算の連載も中盤に差し掛かりました。第5回では、行列の「健康状態」を測るバロメーターである条件数について学び、数値計算がいかに誤差の影響を受けやすいかを理解しました。今回からは、いよいよその行列を用いて「連立一次方程式 \(Ax = b\) を解く」という、数値計算の最重要課題に挑みます。

機械学習の重み学習、物理シミュレーション、最適化問題など、現代のエンジニアリングにおいて $Ax = b$ を解く場面は数え切れないほど存在します。多くのプログラマーは np.linalg.solve や torch.linalg.solve といった関数を一行書くだけで済ませていますが、その内部でどのようなドラマが繰り広げられているのかを知ることは、計算効率の向上や、不安定な計算結果のデバッグにおいて極めて重要です。

今回は、最も古典的でありながら強力な直截解法（Direct Method）である「ガウス消去法（Gaussian Elimination）」のアルゴリズムを、一歩ずつ丁寧に実装していきましょう。今回もよろしくお願いします。

1. 連立一次方程式の解法

連立一次方程式を解くアルゴリズムは、大きく分けて2つのカテゴリーがあります。

ひとつは「直截解法（Direct Method）」です。これは、浮動小数点の丸め誤差を無視すれば、あらかじめ決まった有限の手順で「正確な解」が得られる手法です。今回扱うガウス消去法がその代表例です。

もうひとつは「反復解法（Iterative Method）」です。これは適当な初期値から始めて、徐々に解に近づけていく手法で、非常に大規模な（行列がスカスカな）問題に適しています。これについては連載の後半（第13回など）で詳しく扱います。

小規模から中規模の「密な行列（要素のほとんどが0でない行列）」に対しては、直截解法が最も信頼できる選択肢となります。

2. 前準備：三角行列を係数行列とする方程式

行列 \(A\) がそのままの形では解きにくくても、もし「上三角行列（対角線より下がすべて0の行列）」であれば、解を求めるのは非常に簡単です。

例えば、以下のようなシステムを考えてみましょう。
$$
\begin{cases}
a_{11}x_1 + a_{12}x_2 + a_{13}x_3 = b_1 \\
0x_1 + a_{22}x_2 + a_{23}x_3 = b_2 \\
0x_1 + 0x_2 + a_{33}x_3 = b_3
\end{cases}
$$

一番下の行を見れば、\(x_3 = b_3 / a_{33}\) とすぐに解けます。その結果を二番目の式に代入すれば \(x_2\) が求まり、さらに一番目の式に代入すれば \(x_1\) が求まります。このプロセスを「後退代入（Back Substitution）」と呼びます。

後退代入のPyTorch実装

ガウス消去法の後半戦で必要になるため、まずはこの「解きやすい形」の解法を実装しておきましょう。

Python
import torch

def back_substitution(U, b):
    """
    上三角行列 U とベクトル b に対して Ux = b を解く

    Args:
        U (torch.Tensor): n x n の上三角行列
        b (torch.Tensor): n 次元のベクトル

    Returns:
        torch.Tensor: 解ベクトル x
    """
    n = len(b)
    x = torch.zeros_like(b)

    # 下の行から順番に解いていく
    for i in range(n - 1, -1, -1):
        if U[i, i] == 0:
            raise ValueError("行列が正則ではありません（対角成分に0があります）")

        # すでに求まった x_j (j > i) の影響を b から引く
        sum_val = torch.dot(U[i, i+1:], x[i+1:])
        x[i] = (b[i] - sum_val) / U[i, i]

    return x

# テスト
U_test = torch.tensor([[2.0, 1.0], [0.0, 3.0]])
b_test = torch.tensor([8.0, 9.0])
# 3x = 9 -> x=3, 2x + 1(3) = 8 -> 2x=5 -> x=2.5
print(f"Back Substitution Result: {back_substitution(U_test, b_test)}")
Plaintext
Back Substitution Result: tensor([2.5000, 3.0000])

3. ガウス消去法の基本：前進消去

ガウス消去法の「前進消去（Forward Elimination）」は、連立一次方程式 \(Ax = b\) の係数行列 \(A\) を「上三角行列」（対角線より下がすべて0の行列）に変形するためのプロセスです。これにより、未知数を下の行から順に一つずつ解いていく「後退代入」が可能になります。

この変形は「行の基本変形」（ある行に他の行の定数倍を加える操作）を繰り返すことで実現します。具体的には、各列ごとに、その列の先頭（対角成分）以外の下側の要素をすべて0にしていきます。

なぜ \(A\) だけでなく \(b\) にも同じ操作をするのか？
それは、\(Ax = b\) の両辺に同じ行操作を施さないと、方程式の意味が変わってしまうからです。計算の便宜上、\(A\) と \(b\) を横に並べた「増大行列（Augmented Matrix）」として扱い、同時に操作します。

【手順】

  1. まず1列目について、1行目を基準（ピボット行）とし、2行目以降の各行から「1行目の定数倍」を引いて、1列目の2行目以降の要素をすべて0にします。
    • 具体的には、2行目の1列目の値を0にしたい場合、2行目から「(2行目の1列目 ÷ 1行目の1列目) × 1行目」を引きます。
    • 3行目、4行目…も同様に処理します。
  2. 次に2列目について、2行目を基準（ピボット行）とし、3行目以降の2列目の要素を0にします。
    • 3行目の2列目を0にしたい場合、3行目から「(3行目の2列目 ÷ 2行目の2列目) × 2行目」を引きます。
    • 4行目以降も同様です。
  3. これを \(n-1\) 列目まで繰り返します。
    • つまり、k列目ではk行目を基準に、k+1行目以降のk列目の要素を0にしていきます。

この操作を繰り返すことで、\(A\) は上三角行列の形になり、\(b\) もそれに対応した形に変形されます。

【ポイント】

  • 各ステップで「上から下へ」順番に処理することで、下三角部分がどんどん0になっていきます。
  • 途中で割り算が発生するため、ピボット（基準となる要素）が0や極端に小さい場合は注意が必要です（この対策が「ピボット選択」）。
  • 増大行列を使うことで、\(A\)と\(b\)を同時に正しく変形できます。

このようにして、連立一次方程式を「解きやすい形」に持ち込むのが前進消去の本質です。

4. 数値的安定性の要：ピボット選択（Pivoting）

前進消去の過程では、各ステップで「ピボット」と呼ばれる基準要素 $a_{kk}$ で他の行を割る操作が必ず発生します。
もし \(a_{kk}\) が 0 だった場合、その列の消去ができず、プログラムはクラッシュします。
また、\(a_{kk}\) が0でなくても極端に小さい値（例：\(10^{-15}\)）だと、割り算によって誤差が大きく増幅され、計算結果が信頼できなくなります（第5回で学んだ「条件数の悪化」と同じ現象）。

このような「ゼロ除算」や「誤差の増幅」を防ぐために導入されるのが「ピボット選択（Pivoting）」です。

【ピボット選択の仕組みと意味】

  • 各ステップ（k列目）で、これから消去操作を行う列の中から「絶対値が最大の要素」を持つ行を探します。
  • その行を現在のk行目と入れ替え、必ず大きな値で割り算を行うようにします。
  • これにより、割り算による誤差の増幅を最小限に抑え、計算の安定性が劇的に向上します。

【部分ピボット選択と完全ピボット選択】

  • 本文で説明しているのは「部分ピボット選択（Partial Pivoting）」で、各列ごとに下側の行だけを見て最大値を探します。
  • さらに厳密な「完全ピボット選択（Complete Pivoting）」では、行だけでなく列も入れ替えて、全体で最大の要素をピボットにします（ただし計算コストが増えるため、実務では部分ピボットが主流です）。

【なぜピボット選択が重要なのか】

  • ピボット選択をしないと、たまたま対角成分が0や極小値になるだけで計算が破綻したり、丸め誤差が爆発して全く正しくない解が出ることがあります。
  • 実際の数値計算ライブラリ（NumPy, PyTorch, MATLABなど）でも、ガウス消去法やLU分解の内部では必ずピボット選択が組み込まれています。
  • ピボット選択は「計算の安定性」を守るための最重要テクニックの一つです。

このように、ピボット選択はガウス消去法を「実用的なアルゴリズム」にするために不可欠な工夫です。数値計算の現場では、必ず意識しておきましょう。

5. 実践：PyTorchによるガウス消去法の完全実装

それでは、ピボット選択を含むガウス消去法の全体像を実装しましょう。

Python
def gaussian_elimination(A, b):
    """
    ガウス消去法（部分ピボット選択付き）を用いて Ax = b を解く

    Args:
        A (torch.Tensor): n x n の係数行列
        b (torch.Tensor): n 次元のベクトル

    Returns:
        torch.Tensor: 解ベクトル x
    """
    # Aとbをコピーして破壊的変更を避ける
    A = A.clone().float()
    b = b.clone().float()
    n = len(b)

    # --- 前進消去 (Forward Elimination) ---
    for k in range(n):
        # 1. 部分ピボット選択: k列目の中で絶対値が最大の行を探す
        pivot_idx = torch.argmax(torch.abs(A[k:, k])) + k

        # 行列が正則か（解けるか）の判定
        if torch.abs(A[pivot_idx, k]) < 1e-12:
            raise ValueError("行列は正則でない（または非常に悪条件）ため解けません")

        # 2. 行の入れ替え (A と b 両方)
        if pivot_idx != k:
            # Aのk行目とpivot_idx行目をスワップ
            A[k], A[pivot_idx] = A[pivot_idx].clone(), A[k].clone()
            # bのk要素目とpivot_idx要素目をスワップ
            b[k], b[pivot_idx] = b[pivot_idx].clone(), b[k].clone()

        # 3. 各行から引き算して下を0にする
        for i in range(k + 1, n):
            factor = A[i, k] / A[k, k]
            # i行目の計算
            A[i, k:] -= factor * A[k, k:]
            b[i] -= factor * b[k]

    # --- 後退代入 (Back Substitution) ---
    # 前進消去で A は上三角行列 U になっている
    x = back_substitution(A, b)

    return x

# 実験: 条件数が少し悪い行列でテスト
A_sample = torch.tensor([[1.0, 1.0, 1.0],
                         [2.0, 3.0, 5.0],
                         [4.0, 0.0, 5.0]])
b_sample = torch.tensor([6.0, 15.0, 14.0])

# 自作関数で解く
x_sol = gaussian_elimination(A_sample, b_sample)
print(f"Manual Gaussian Result: {x_sol}")

# PyTorchの組み込み関数と比較
x_torch = torch.linalg.solve(A_sample, b_sample)
print(f"PyTorch solve Result:  {x_torch}")
Plaintext
Manual Gaussian Result: tensor([3.3077, 2.5385, 0.1538])
PyTorch solve Result:  tensor([3.3077, 2.5385, 0.1538])

コードの解説

  • ピボット選択: torch.argmax を使って、現在の列以降で最も大きな要素を見つけています。
  • メモリ格納: 行列の行スワップの際、PyTorchの高度なインデックス指定 A[[k, pivot_idx]] = A[[pivot_idx, k]] を使うことで、簡潔に記述しています。
  • 正則性判定: ピボット要素が極めて小さい場合、その行列は実質的に「逆行列を持たない」とみなしてエラーを投げています。

6. 誤差と計算量の現実

このアルゴリズムの計算量を考えてみましょう。
前進消去のプロセスには、2重のループの中で各行の引き算（さらにループ）があるため、全体として \(O(n^3)\) の計算量が必要になります。

  • \(n=100\) なら約100万回の演算。
  • \(n=10,000\) なら約1兆回の演算。

このように、\(n\) が大きくなると計算時間は爆発的に増えます。また、第5回で触れた通り、ステップ数が増えるほど「丸め誤差」が蓄積していきます。そのため、数万次元を超えるような巨大な連立方程式を解く場合は、この直截解法ではなく、もっとメモリ効率が良く誤差の蓄積が少ない手法（反復法やスパース行列用アルゴリズム）を選択する必要があります。

7. まとめ

第6回では、連立一次方程式を解くための最も基本的な武器である「ガウス消去法」を学びました。

  • 前進消去によって行列を「解きやすい」上三角行列に変形する。
  • 後退代入によって下の変数から順に値を確定させる。
  • ピボット選択は、ゼロ除算を避けるだけでなく、数値的な精度を保つために絶対に欠かせない。
  • 計算量は \(O(n^3)\) であり、行列のサイズに対してコストが急激に増大する。

今回実装したプロセスは、実は次回扱う「行列分解（LU分解）」の基礎にもなっています。行列を一度「分解」しておけば、右辺の \(b\) が変わるたびに \(O(n^3)\) の重い計算をやり直す必要がなくなるのです。

次回は、より洗練された行列の扱い方である「LU分解とコレスキー分解」について解説します。エンジニアが実際のコードで最もよく目にする「行列の形」に迫ります。お楽しみに。

最後まで読んでいただきありがとうございます。

連載記事、各回へのリンク

Pythonで学ぶ数値計算のアルゴリズムと実装 第1回：多項式計算とテイラー展開 | アマチュア無線局JS2IIU
https://js2iiu.com/2025/12/31/python_numerical_01/

Pythonで学ぶ数値計算のアルゴリズムと実装 第2回：多項式補間とチェビシェフ補間の戦略 | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/01/python_numerical_02/

Pythonで学ぶ数値計算のアルゴリズムと実装 第3回：積分の基礎：中点則、台形則からシンプソン則まで | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/01/python_numerical_03/

Pythonで学ぶ数値計算のアルゴリズムと実装 第4回：高度な数値積分：二重指数型（DE）公式と広義積分 | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/01/python_numerical_04/

Pythonで学ぶ数値計算のアルゴリズムと実装 第5回：計算の「正しさ」を評価する：ノルムと条件数 | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/01/python_numeric_05/

Pythonで学ぶ数値計算のアルゴリズムと実装 第6回：連立一次方程式の直截解法：ガウス消去法とピボット選択 | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/01/python_numerical_06/

Pythonで学ぶ数値計算のアルゴリズムと実装 第8回：スパース性を活かす：帯行列の高速計算アルゴリズム | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/02/python_numerical_08/

Pythonで学ぶ数値計算のアルゴリズムと実装 第9回：最小二乗法とQR分解：ハウスホルダー変換による安定化 | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/02/python_numerial_09/

Pythonで学ぶ数値計算のアルゴリズムと実装 第10回：非線形方程式の探索：2分法からニュートン法まで | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/03/python_numerical_10/

Pythonで学ぶ数値計算のアルゴリズムと実装 第11回：行列の本質を捉える：固有値問題の数値解法 | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/03/python_numerical_11/

Pythonで学ぶ数値計算のアルゴリズムと実装 第12回：現象の変化を追う：ルンゲ・クッタ法による常微分方程式の解法 | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/03/python_numerical_12/

Pythonで学ぶ数値計算のアルゴリズムと実装 第13回：空間の熱と振動を解く：偏微分方程式の差分近似と反復解法 | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/03/python_numerical_13/

コメント

タイトルとURLをコピーしました