広告 プログラミング

【解説＆実装】ガウスの消去法

※本ページには、プロモーション（広告）が含まれています。

悩んでいる人
悩んでいる人

数値解析の講義でガウスの消去法について学習したが、実装しようと思うと手が止まってしまう。

考え方や実装例を教えて欲しい。

こんなお悩みを解決します。

今回は、数値解析を実施する上で、一度は目にすることになる「ガウスの消去法」について解説したいと思います。

これから数値計算を学習し始める学生や、再度勉強し始めた社会人の方にとって、参考になればよいと思っています。

考え方

ガウスの消去法では、行列表示された連立一次方程式の係数行列を上三角行列に変形（前進消去）し、今度は逆方向に解を求めていく（後退代入）という手法をとります。

ここでは、以下の連立一次方程式を例に取り上げ、ガウスの消去法について説明します。

$$
\begin{eqnarray*}
\left(
\begin{array}{ccc}
2 & 1 & -2\\
1 & 1 & -1\\
1 & -2 & 3\\
\end{array}
\right)\left(
\begin{array}{c}
x_{1}\\
x_{2}\\
x_{3}\\
\end{array}
\right) = \left(
\begin{array}{c}
1\\
4\\
-1\\
\end{array}
\right)
\end{eqnarray*}
$$

前進消去

前進消去では、係数行列の対角成分より下側を0にする操作となります。

\((2, 1)\)成分の消去

\((2, 1)\)成分の消去は、以下のように計算します。

  • 1列目に関して、\(\frac{(2, 1)\text{成分}}{(1, 1)\text{成分}}\)を計算します。
  • 2行目から\(\frac{(2, 1)\text{成分}}{(1, 1)\text{成分}} \times 1\text{行目}\)を引きます。

ここで、\(\frac{(2, 1)\text{成分}}{(1, 1)\text{成分}} = \frac{1}{2}\)ですので、2行目を式変形すると、以下のようになります。

$$
\begin{eqnarray*}
&& \left(
\begin{array}{ccc}
2 & 1 & -2\\
1 & 1 & -1\\
1 & -2 & 3\\
\end{array}
\right)\left(
\begin{array}{c}
x_{1}\\
x_{2}\\
x_{3}\\
\end{array}
\right) = \left(
\begin{array}{c}
1\\
4\\
-1\\
\end{array}
\right)\\
& \Leftrightarrow & \left(
\begin{array}{ccc}
2 & 1 & -2\\
1 - \frac{1}{2} \cdot 2 & 1 - \frac{1}{2} \cdot 1 & -1 - \frac{1}{2} \cdot (-2)\\
1 & -2 & 3\\
\end{array}
\right)\left(
\begin{array}{c}
x_{1}\\
x_{2}\\
x_{3}\\
\end{array}
\right) = \left(
\begin{array}{c}
1\\
4 - \frac{1}{2} \cdot 1\\
-1\\
\end{array}
\right)\\
& \Leftrightarrow & \left(
\begin{array}{ccc}
2 & 1 & -2\\
0 & \frac{1}{2} & 0\\
1 & -2 & 3\\
\end{array}
\right)\left(
\begin{array}{c}
x_{1}\\
x_{2}\\
x_{3}\\
\end{array}
\right) = \left(
\begin{array}{c}
1\\
\frac{7}{2} \\
-1\\
\end{array}
\right)
\end{eqnarray*}
$$

\((3, 1)\)成分の消去

同様に、\((2, 2)\)成分の消去は、以下のように計算します。

  • 1列目に関して、\(\frac{(3, 1)\text{成分}}{(1, 1)\text{成分}}\)を計算します。
  • 3行目から\(\frac{(3, 1)\text{成分}}{(1, 1)\text{成分}} \times 1\text{行目}\)を引きます。

ここで、\(\frac{(3, 1)\text{成分}}{(1, 1)\text{成分}} = \frac{1}{2}\)ですので、3行目を式変形すると、以下のようになります。

$$
\begin{eqnarray*}
&& \left(
\begin{array}{ccc}
2 & 1 & -2\\
0 & \frac{1}{2} & 0\\
1 & -2 & 3\\
\end{array}
\right)\left(
\begin{array}{c}
x_{1}\\
x_{2}\\
x_{3}\\
\end{array}
\right) = \left(
\begin{array}{c}
1\\
\frac{7}{2}\\
-1\\
\end{array}
\right)\\
& \Leftrightarrow & \left(
\begin{array}{ccc}
2 & 1 & -2\\
0 & \frac{1}{2} & 0\\
1 - \frac{1}{2} \cdot 2 & -2 - \frac{1}{2} \cdot 1 & 3 - \frac{1}{2} \cdot (-2)\\
\end{array}
\right)\left(
\begin{array}{c}
x_{1}\\
x_{2}\\
x_{3}\\
\end{array}
\right) = \left(
\begin{array}{c}
1\\
\frac{7}{2}\\
-1 - \frac{1}{2} \cdot 1\\
\end{array}
\right)\\
& \Leftrightarrow & \left(
\begin{array}{ccc}
2 & 1 & -2\\
0 & \frac{1}{2} & 0\\
0 & -\frac{5}{2} & 4\\
\end{array}
\right)\left(
\begin{array}{c}
x_{1}\\
x_{2}\\
x_{3}\\
\end{array}
\right) = \left(
\begin{array}{c}
1\\
\frac{7}{2}\\
-\frac{3}{2}\\
\end{array}
\right)
\end{eqnarray*}
$$

\((3, 2)\)成分の消去

\((3, 2)\)成分の消去は、以下のように計算します。

  • 2列目に関して、\(\frac{(3, 2)\text{成分}}{(2, 2)\text{成分}}\)を計算します。
  • 3行目から\(\frac{(3, 2)\text{成分}}{(2, 2)\text{成分}} \times 2\text{行目}\)を引きます。

ここで、\(\frac{(3, 2)\text{成分}}{(2, 2)\text{成分}} = -5\)ですので、3行目を式変形すると、以下のようになります。

$$
\begin{eqnarray*}
&& \left(
\begin{array}{ccc}
2 & 1 & -2\\
0 & \frac{1}{2} & 0\\
0 & -\frac{5}{2} & 4\\
\end{array}
\right)\left(
\begin{array}{c}
x_{1}\\
x_{2}\\
x_{3}\\
\end{array}
\right) = \left(
\begin{array}{c}
1\\
\frac{7}{2}\\
-\frac{3}{2}\\
\end{array}
\right)\\
& \Leftrightarrow & \left(
\begin{array}{ccc}
2 & 1 & -2\\
0 & \frac{1}{2} & 0\\
0 & -\frac{5}{2} - (-5) \cdot \frac{1}{2} & 4 - (-5) \cdot 0\\
\end{array}
\right)\left(
\begin{array}{c}
x_{1}\\
x_{2}\\
x_{3}\\
\end{array}
\right) = \left(
\begin{array}{c}
1\\
\frac{7}{2}\\
-\frac{3}{2} - (-5) \cdot \frac{7}{2} \\
\end{array}
\right)\\
& \Leftrightarrow & \left(
\begin{array}{ccc}
2 & 1 & -2\\
0 & \frac{1}{2} & 0\\
0 & 0 & 4\\
\end{array}
\right)\left(
\begin{array}{c}
x_{1}\\
x_{2}\\
x_{3}\\
\end{array}
\right) = \left(
\begin{array}{c}
1\\
\frac{7}{2}\\
16\\
\end{array}
\right)
\end{eqnarray*}
$$

後退代入

後退代入では、行列の下から上に向かって計算することで、未知数を求めていきます。前進消去した結果をもとに、一度連立一次方程式の形に直して考えたいと思います。

連立一次方程式としての表現

連立一次方程式の形で表現した結果を以下に示します。

$$
\left\{
\begin{array}{rcc}
2 x_{1} + x_{2} -2 x_{3} & = & 1\\
\frac{1}{2}x_{2} + 0x_{3} & = & \frac{7}{2}\\
4x_{3} & = & 16\\
\end{array}
\right.
$$

\(x_{3}\)を求める

連立一次方程式の3行目に注目します。未知変数は\(x_{3}\)のみなので、\(x_{3}\)の値を求めることができます。

$$
\begin{eqnarray*}
& & 4x_{3} = 16\\
& \Leftrightarrow & x_{3} = 4
\end{eqnarray*}
$$

\(x_{2}\)を求める

連立一次方程式の2行目に注目します。\(x_{3}\)の値はすでに得られているため、未知変数\(x_{2}\)を求めることができます。

$$
\begin{eqnarray*}
& & \frac{1}{2}x_{2} + 0 \cdot 4 = \frac{7}{2}\\
& \Leftrightarrow & x_{2} = 7
\end{eqnarray*}
$$

\(x_{1}\)を求める

連立一次方程式の1行目に注目します。\(x_{2}, x_{3}\)の値はすでに得られているため、未知変数\(x_{1}\)を求めることができます。

$$
\begin{eqnarray*}
& & 2 x_{1} + 7 -2 \cdot 4 & = & 1\\
& \Leftrightarrow & 2 x_{1} = 1 - 7 + 8\\
& \Leftrightarrow & x_{1} = 1
\end{eqnarray*}
$$

疑似コード

考え方では、\(3 \times 3\)行列で考えましたが、一般化すると以下のようなアルゴリズムになります。

Input A, b
    # A: a[0][0] ~ a[n - 1][n - 1]
    # b: b[0] ~ b[n - 1]

    # 前進消去
    for k = 0 To n - 1
        for i = k + 1 To n - 1
            w ← a[i][k] / a[k][k]

            for j = k To n - 1
                a[i][j] ← a[i][j] - a[k][j] * w
            end for
        end for
    end for

    # 後退代入
    for i = n - 1 To 0
        sum ← b[i]

        for j = i + 1 To n
            sum ←  sum - a[i][j] * b[j]
        end for
        b[i] ← sum / a[i][i]
    end for

Output b

C言語にする実装

実装結果は、GitHubにあります。必要に応じて参照してください。

https://github.com/yuruto-free/gaussian-elimination

gaussian_elimination.h

#ifndef GAUSSIAN_ELIMINATION_H__
#define GAUSSIAN_ELIMINATION_H__

#include 

#define GE_OK (0)
#define GE_NG (1)

/**
 * \brief ガウスの消去法
 * \param[in]    dim    次元数
 * \param[in]    matrix 係数行列
 *                      サイズ：dim * dim
 * \param[inout] vec    右辺ベクトル
 *                      サイズ：dim
 * \retval       GE_OK  正常終了
 * \retval       GE_NG  異常終了
 */
int32_t GE_gauss_solver(int32_t dim, double *matrix, double *vec);

#endif

gaussian_elimination.c

#include "gaussian_elimination.h"
#include 
#include 
#include 

#define RETURN_OK_GE (0)
#define RETURN_NG_GE (1)
#define MEPS_GE (1e-10)

/**
 * \brief ピボット選択
 * \param[in]  dim    次元数
 * \param[in]  col    対象列
 * \param[in]  matrix 係数行列
 * \param[out] pivot  ピボット値
 * \retval RETURN_OK_GE 正常終了
 * \retval RETURN_NG_GE 異常終了
 */
static int32_t pivot_selection(int32_t dim, int32_t col, double *matrix,
                               int32_t *pivot);
/**
 * \brief 入れ替え処理
 * \param[in]    dim    次元数
 * \param[in]    row    入れ替え元の行
 * \param[in]    pivot  入れ替え先の行
 * \param[inout] matrix 係数行列
 * \param[inout] vec    右辺ベクトル
 * \retval RETURN_OK_GE 正常終了
 * \retval RETURN_NG_GE 異常終了
 */
static int32_t swap(int32_t dim, int32_t row, int32_t pivot, double *matrix,
                    double *vec);

int32_t GE_gauss_solver(int32_t dim, double *matrix, double *vec) {
    int32_t ret = (int32_t)GE_NG;
    int32_t func_val;
    int32_t pivot, idx;
    int32_t row, col;
    double diag, scale, sum;

    if ((NULL != matrix) && (NULL != vec)) {
        /* 前進消去 */
        for (idx = 0; idx < dim - 1; idx++) {
            /* ピボット選択 */
            func_val = pivot_selection(dim, idx, matrix, &pivot);
            if ((int32_t)RETURN_OK_GE != func_val) {
                goto EXIT_GE_SOLVER;
            }
            /* 対角成分の絶対値が最大でない場合、行を入れ替える */
            if (idx != pivot) {
                func_val = swap(dim, idx, pivot, matrix, vec);
                if ((int32_t)RETURN_OK_GE != func_val) {
                    goto EXIT_GE_SOLVER;
                }
            }
            /* 対角成分の最大値を取得 */
            diag = matrix[idx * dim + idx];

            if (fabs(diag) < (double)MEPS_GE) {
                goto EXIT_GE_SOLVER;
            }
            for (row = idx + 1; row < dim; row++) {
                scale = matrix[row * dim + idx] / diag;

                for (col = idx; col < dim; col++) {
                    matrix[row * dim + col] -= matrix[idx * dim + col] * scale;
                }
                vec[row] -= vec[idx] * scale;
            }
        }

        /* 後退代入 */
        for (row = dim - 1; row >= 0; row--) {
            sum = vec[row];

            for (col = row + 1; col < dim; col++) {
                sum -= matrix[row * dim + col] * vec[col];
            }
            vec[row] = sum / matrix[row * dim + row];
        }
        ret = (int32_t)GE_OK;
    }
EXIT_GE_SOLVER:

    return ret;
}

static int32_t pivot_selection(int32_t dim, int32_t col, double *matrix,
                               int32_t *pivot) {
    int32_t ret = (int32_t)RETURN_NG_GE;
    int32_t i;
    int32_t p;
    double max_val, val;

    if ((NULL != matrix) && (NULL != pivot)) {
        /* 対角成分の絶対値が最大値となると仮定 */
        max_val = fabs(matrix[col * dim + col]);
        p       = col;

        /* 対角成分より下側の値のうち、絶対が最大となるものを探索 */
        for (i = col + 1; i < dim; i++) {
            val = fabs(matrix[i * dim + col]);

            if (max_val < val) {
                max_val = val;
                p       = i;
            }
        }
        (*pivot) = p;
        ret      = (int32_t)RETURN_OK_GE;
    }

    return ret;
}

static int32_t swap(int32_t dim, int32_t row, int32_t pivot, double *matrix,
                    double *vec) {
    int32_t ret = (int32_t)RETURN_NG_GE;
    int32_t col;
    double tmp;

    if ((NULL != matrix) && (NULL != vec)) {
        /* 行の入れ替え */
        for (col = 0; col < dim; col++) {
            tmp                       = matrix[row * dim + col];
            matrix[row * dim + col]   = matrix[pivot * dim + col];
            matrix[pivot * dim + col] = tmp;
        }
        /* 右辺ベクトルの入れ替え */
        tmp        = vec[row];
        vec[row]   = vec[pivot];
        vec[pivot] = tmp;
        ret        = (int32_t)RETURN_OK_GE;
    }

    return ret;
}

ここで、部分ピボット選択という処理を追加しています。上記で解説していなかったので、この場で紹介します。

部分ピボット選択

ここでは、以下のような連立一次方程式を考えます。

$$
\begin{eqnarray*}
\left(
\begin{array}{ccc}
0 & 0 & 2\\
3 & 3 & 1\\
-2 & 2 & 0\\
\end{array}
\right) \left(
\begin{array}{c}
x_{1}\\
x_{2}\\
x_{3}\\
\end{array}
\right) = \left(
\begin{array}{c}
1\\
2\\
3\\
\end{array}
\right)
\end{eqnarray*}
$$

上記の連立一次方程式の場合、\((1, 1)\)成分が0となっているため、前進消去ができないです。これを解決するには、対象としている列（ここでは、1列目）のうち、対角成分以降で絶対値が最大となるものと入れ替えれば良いです。

ここでは、2行目が最大となるので、1行目と2行目を入れ替えると以下のようになります。

$$
\begin{eqnarray*}
\left(
\begin{array}{ccc}
3 & 3 & 1\\
0 & 0 & 2\\
-2 & 2 & 0\\
\end{array}
\right) \left(
\begin{array}{c}
x_{1}\\
x_{2}\\
x_{3}\\
\end{array}
\right) = \left(
\begin{array}{c}
2\\
1\\
3\\
\end{array}
\right)
\end{eqnarray*}
$$

1列目に対し、前進消去を行うと、以下のようになります。

$$
\begin{eqnarray*}
\left(
\begin{array}{ccc}
3 & 3 & 1\\
0 & 0 & 2\\
0 & 4 & \frac{2}{3}\\
\end{array}
\right) \left(
\begin{array}{c}
x_{1}\\
x_{2}\\
x_{3}\\
\end{array}
\right) = \left(
\begin{array}{c}
2\\
1\\
\frac{13}{3} \\
\end{array}
\right)
\end{eqnarray*}
$$

同様に、\((2, 2)\)成分が0となっているため、行を入れ替える必要があります。ここでは、2行目と3行目を入れ替えれば良いです。入れ替えた結果は以下のようになります。

$$
\begin{eqnarray*}
\left(
\begin{array}{ccc}
3 & 3 & 1\\
0 & 4 & \frac{2}{3}\\
0 & 0 & 2\\
\end{array}
\right) \left(
\begin{array}{c}
x_{1}\\
x_{2}\\
x_{3}\\
\end{array}
\right) = \left(
\begin{array}{c}
2\\
\frac{13}{3} \\
1\\
\end{array}
\right)
\end{eqnarray*}
$$

スポンサードリンク



-プログラミング
-,