ガウスの消去法

1. 概要

ガウスの消去法（Gaussian Elimination）は、係数行列 $A$ が正方で正則なとき、連立一次方程式 $A\mathbf{x} = \mathbf{b}$ の一意の解を求める最も基本的な直接法である（特異な場合は、消去の結果から解なし・無数の解を判定できる）。拡大係数行列 $[A|\mathbf{b}]$ に行基本変形を施して係数部分を上三角行列に変換し（前進消去）、その後、後退代入により解を求める。

C. F. Gauss (1777–1855) にちなんで名付けられたが、同様の手法は紀元前の中国の『九章算術』にも見られる。

2. 前進消去

$n \times n$ の係数行列 $A$ と右辺ベクトル $\mathbf{b}$ から拡大係数行列を構成し、第 $k$ 列（$k = 1, 2, \ldots, n-1$）について、$k$ 行目より下の全ての行から $k$ 行目の定数倍を引いて $a_{ik}$ を消去する。

3. 後退代入

前進消去により、対角成分 $u_{kk}$ がすべて 0 でない上三角形 $U\mathbf{x} = \mathbf{c}$ が得られたら、最後の行から順に解を求める：

4. 計算過程のアニメーション例

図1・図2（5 元）とは別に、ここではより小さな3変数の連立方程式を例に、前進消去と後退代入の全過程をステップごとに確認する。

拡大係数行列による計算例

部分ピボット選択付き C++ 実装例

#include <vector>
#include <cmath>
#include <stdexcept>
#include <limits>
#include <algorithm>

std::vector<double> gaussianElimination(
    std::vector<std::vector<double>> A,
    std::vector<double> b)
{
    const int n = static_cast<int>(A.size());
    // 特異判定の閾値: 行列の最大要素に比例させる（スケールに左右されない）
    double scale = 0.0;
    for (const auto& row : A) for (double v : row) scale = std::max(scale, std::abs(v));
    const double tol = 100 * std::numeric_limits<double>::epsilon() * scale;
    for (int k = 0; k < n - 1; ++k) {
        // 部分ピボット選択: k 列の絶対値最大行を探す
        int pivot = k;
        for (int i = k + 1; i < n; ++i)
            if (std::abs(A[i][k]) > std::abs(A[pivot][k])) pivot = i;
        if (pivot != k) { std::swap(A[k], A[pivot]); std::swap(b[k], b[pivot]); }
        if (std::abs(A[k][k]) <= tol)
            throw std::runtime_error("singular matrix");
        // 消去
        for (int i = k + 1; i < n; ++i) {
            double m = A[i][k] / A[k][k];
            for (int j = k; j < n; ++j) A[i][j] -= m * A[k][j];
            b[i] -= m * b[k];
        }
    }
    if (std::abs(A[n - 1][n - 1]) <= tol)   // 最後のピボットも検査する
        throw std::runtime_error("singular matrix");
    // 後退代入
    std::vector<double> x(n);
    for (int i = n - 1; i >= 0; --i) {
        double s = b[i];
        for (int j = i + 1; j < n; ++j) s -= A[i][j] * x[j];
        x[i] = s / A[i][i];
    }
    return x;
}

5. ピボット選択

前進消去の第 $k$ ステップにおいて、ピボット $a_{kk}^{(k)}$ がゼロまたはゼロに近い場合、ゼロ除算や大きな丸め誤差が生じる。これを避けるためにピボット選択を行う。

5.1 部分ピボット選択（Partial Pivoting）

第 $k$ 列の $k$ 行目以下で絶対値が最大の要素を選び、行の交換を行う。行の交換は連立方程式の式の順番を入れ替えるだけなので、変数の対応関係は変わらない。実用上はほとんどの場合これで十分であり、LAPACKなどの標準的な数値計算ライブラリでも採用されている。

5.2 完全ピボット選択（Complete Pivoting）

完全ピボット選択では、第 $k$ ステップにおいて 係数行列の $k$ 行目以下かつ $k$ 列目から右の全要素（右辺の列は含めない）から絶対値最大のものを探し、行と列の両方を交換する。

列の交換は変数の順序の入れ替えに相当する。例えば第2列と第3列を交換すると、$x_2$ と $x_3$ の役割が入れ替わる。したがって、全ての消去が終わった後に列の交換履歴を逆順にたどって、変数の順序を元に戻す必要がある。

5.3 有理数演算でのピボット選択

上記のピボット選択戦略は浮動小数点演算を前提としている。有理数（多倍長整数による分数）で厳密にガウスの消去法を実行する場合、丸め誤差は発生しないが、代わりに中間値の膨張（coefficient swell）が問題となる。消去の過程で分子・分母の桁数が爆発的に増大し、計算が極端に遅くなることがある。

この問題を軽減するために、浮動小数点とは違って値の大きさ（絶対値）ではなく、分子・分母の大きさに注目してピボットを選ぶ方法がある。有理数 $p/q$（既約分数）の高さ $H(p/q) = \max(|p|, |q|)$ が小さい（ゼロでない）要素をピボットに選ぶと、そのステップの乗数の分子・分母が小さくなりやすい。ただしこれは局所的な目安（ヒューリスティック）であり、全体の中間値の膨張を最小にする保証はない。例えば $1/1000$ は絶対値こそ小さいが高さは $1000$、$2/3$ は絶対値が大きいが高さは $3$ である。高さ最小ピボットでは（絶対値が大きくても高さの小さい）$2/3$ を選ぶ点が、浮動小数点での絶対値最大ピボットとの違いである。

6. 計算量

$n$ 元連立一次方程式に対するガウスの消去法の計算量は次の通りである。

• 前進消去: $\dfrac{2}{3}n^3 + O(n^2)$ 回の浮動小数点演算
• 後退代入: $n^2 + O(n)$ 回の浮動小数点演算
• 全体: $O(n^3)$

7. LU分解との関係

行交換をせずに消去できる場合（$A$ の先頭主小行列式がすべて 0 でない場合）、ガウスの消去法の前進消去は、行列 $A$ を下三角行列 $L$ と上三角行列 $U$ の積に分解する操作と等価である：

消去の乗数 $m_{ik}$ が $L$ の対角より下の要素となり、消去後の上三角行列が $U$ である。部分ピボット選択で行を交換する場合は、置換行列 $P$ を用いて $PA = LU$、完全ピボット選択で列も交換する場合は $PAQ = LU$ と表す（行を交換するときは、それまでの乗数も同じ行へ並べ替える）。一度 $LU$ 分解を行えば、異なる右辺 $\mathbf{b}$ に対して $O(n^2)$ で解が求まるため、同じ係数行列で複数の右辺を解く場合に効率的である。

8. 注意点

• 要素成長率と数値安定性: 部分ピボット選択付きガウスの消去法の数値安定性は、消去過程での要素成長率 $\rho_n = \max_{i,j,k} |a_{ij}^{(k)}| / \max_{i,j} |a_{ij}|$ に左右される。理論最悪値は $\rho_n \le 2^{n-1}$ で Wilkinson 行列により達成されるが、多くの実用例やランダム行列では経験的に小さく（$n$ 程度以下に）収まることが多い（保証ではない）。完全ピボットでは $\rho_n \le n^{1/2} (2 \cdot 3^{1/2} \cdots n^{1/(n-1)})^{1/2}$ と緩やかに増大する。
• スケーリング感受性: 行のスケールが極端に異なる行列では、絶対値最大という基準だけでは不適切なピボットを選ぶことがある。各行を $\|A_{i,:}\|_\infty$ で正規化してから判定する「スケール付き部分ピボット選択」を採用する実装もある。
• 特異・準特異の検出: $|a_{kk}^{(k)}| \le \tau \cdot \max_{i,j} |a_{ij}|$（$\tau$ は上の実装上のアドバイスと同じく計算機イプシロンの数十〜100 倍程度。行列の次元・スケーリング・必要な精度に応じて選ぶ目安）になったら、行列は数値的に特異とみなし、最小二乗法や SVD など別アルゴリズムに切替えるのが安全である。
• 条件数による精度限界: 計算された解 $\hat{\mathbf{x}}$ の相対誤差は概ね $\kappa(A) \cdot u$（$u$ は丸めの単位）のオーダーであり、$\kappa(A) \approx 1/u$ になると有効桁が完全に失われる。事前に条件数推定（LAPACK の ?gecon 等）を行うのが望ましい。
• 反復改良 (iterative refinement): 残差 $\mathbf{r} = \mathbf{b} - A\hat{\mathbf{x}}$ を高精度で計算し、$A\mathbf{d} = \mathbf{r}$ を $LU$ 分解で解いて $\hat{\mathbf{x}} \leftarrow \hat{\mathbf{x}} + \mathbf{d}$ を 1〜2 回反復すれば、$\kappa(A) \cdot u \ll 1$ の下で精度を取り戻せる。
• 密行列前提: ガウスの消去法は密行列向けで、疎行列に直接適用するとフィルイン（消去で 0 が非 0 に変わる現象）により疎性が失われる。

9. よくある質問

10. 参考資料

• Wikipedia「ガウスの消去法」（日本語版）
• Wikipedia「Gaussian elimination」（英語版）
• G. H. Golub & C. F. Van Loan, Matrix Computations, 4th ed., Johns Hopkins, 2013（§3 で $LU$ 分解、§3.4 で要素成長率と Wilkinson 行列）。
• L. N. Trefethen & D. Bau III, Numerical Linear Algebra, SIAM, 1997（Lecture 20-23 で $LU$ と部分ピボット）。
• N. J. Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed., SIAM, 2002（§9 でガウスの消去法の後退誤差解析、§12 で反復改良）。