ガウスの消去法

ガウスの消去法(Gaussian elimination)と呼ばれる方法では, 前進消去と後退代入という2段階の計算により解を求める.

一般的な(中学校で習う)解き方の一つに加減法というものがある. たとえば,以下のような連立方程式で考える.

ls_gauss_elimination.eq1.gif

加減法では片方の式の両辺に実数をかけて,係数を合わせて,加減算する. ここでは下の式の両辺を2倍して,上の式から引くことにする.

ls_gauss_elimination.eq2.gif

あとは,ls_gauss_elimination.eq3.gifの結果を1番目の式に代入することで,ls_gauss_elimination.eq4.gifを得る.

これを行列で表すと,

ls_gauss_elimination.eq5.gif

に対して,2行目にls_gauss_elimination.eq6.gifをかけて引くことで,

ls_gauss_elimination.eq7.gif

という形にする.そして,一番下の行から未知数を解いて代入していく.

ls_gauss_elimination.eq8.gif

前半の処理では係数行列ls_gauss_elimination.eq9.gifの対角要素を除いた左下部分がすべて0となる 以下に示すような上三角行列(upper triangle matrix)を求めている.

ls_gauss_elimination.eq10.gif

この処理のことを前進消去(forward elimination)と呼ぶ. 前進消去により,未知数ls_gauss_elimination.eq11.gifls_gauss_elimination.eq12.gifで求められ, このls_gauss_elimination.eq11.gifls_gauss_elimination.eq13.gif行目に代入することでls_gauss_elimination.eq14.gifが得られ, ls_gauss_elimination.eq15.gifls_gauss_elimination.eq16.gif行目に代入することでls_gauss_elimination.eq17.gifが得られる. これをls_gauss_elimination.eq18.gifまで繰り返すことですべての解を得る. この後半部分の処理のことを後退代入(back substitution)と呼ぶ. 前進消去と後退代入により連立1次方程式の解を求める方法をガウスの消去法という.

ガウスの消去法の実装

ガウスの消去法をls_gauss_elimination.eq19.gif元連立1次方程式に適用する. まず,ls_gauss_elimination.eq9.gifに定数ベクトルls_gauss_elimination.eq20.gifを列の最後に追加したls_gauss_elimination.eq21.gifの行列を考える.

ls_gauss_elimination.eq22.gif

この行列は拡大行列と呼ばれる.

拡大行列の要素を

ls_gauss_elimination.eq23.gif

とする.ここで,C,C++で実装することを考えて要素番号を0始まりにしている. 今後,実装の説明では0で始まるインデックスを用いるので注意.

前進消去は次の式で表される.

ls_gauss_elimination.eq24.gif

ここで,

ls_gauss_elimination.eq25.gif

上式によりls_gauss_elimination.eq26.gifls_gauss_elimination.eq27.gifで更新することで上三角行列を得る. この処理のC++のコードを以下に示す.

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
    // 前進消去(forward elimination)
    //  - 対角要素をのぞいた左下要素をすべて0にする(上三角行列にする)
    for(int k = 0; k < n-1; ++k){
        double akk = A[k][k];
        for(int i = k+1; i < n; ++i){
            double aik = A[i][k];
            for(int j = k; j < n+1; ++j){ // 確認のため左下要素が0になるようにj=kとしたが,実際にはj=k+1でよい
                A[i][j] = A[i][j]-aik*(A[k][j]/akk);
            }
        }
    }

2次元配列A[n][n+1]には係数行列と定数ベクトルが格納されている.

次に後退代入の式を以下に示す.

ls_gauss_elimination.eq28.gif

後退代入のC++のコード例を以下に示す.

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
    // 後退代入(back substitution)
    //  - x_nの解はb_n/a_nn,x_nをさらにn-1行の式に代入することでx_(n-1)を求める.
    //  - この作業を1行目まで続けることですべての解を得る.
    A[n-1][n] = A[n-1][n]/A[n-1][n-1];
    for(int i = n-2; i >= 0; --i){
        double ax = 0.0;
        for(int j = i+1; j < n; ++j){
            ax += A[i][j]*A[j][n];
        }
        A[i][n] = (A[i][n]-ax)/A[i][i];
    }

計算された結果は別の配列でなくA[i][n]に格納しているので注意.

ガウスの消去法の全体のコードは以下.

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
/*!
 * ガウスの消去法(ピボット交換なし)
 * @param[inout] A n×nの係数項とn×1の定数項(b)を併せたn×(n+1)の行列.n+1列目に解が入る.
 * @param[in] n n元連立一次方程式
 * @return 1:成功
 */
int GaussElimination(vector< vector<double> > &A, int n)
{
    // 前進消去(forward elimination)
    //  - 対角要素をのぞいた左下要素をすべて0にする(上三角行列にする)
    for(int k = 0; k < n-1; ++k){
        double akk = A[k][k];
        for(int i = k+1; i < n; ++i){
            double aik = A[i][k];
            for(int j = k; j < n+1; ++j){ // 確認のため左下要素が0になるようにj=kとしたが,実際にはj=k+1でよい
                A[i][j] = A[i][j]-aik*(A[k][j]/akk);
            }
        }
    }
 
    // 後退代入(back substitution)
    //  - x_nの解はb_n/a_nn,x_nをさらにn-1行の式に代入することでx_(n-1)を求める.
    //  - この作業を1行目まで続けることですべての解を得る.
    A[n-1][n] = A[n-1][n]/A[n-1][n-1];
    for(int i = n-2; i >= 0; --i){
        double ax = 0.0;
        for(int j = i+1; j < n; ++j){
            ax += A[i][j]*A[j][n];
        }
        A[i][n] = (A[i][n]-ax)/A[i][i];
    }
 
    return 1;
}

標準の2次元配列の代わりにSTLのvectorを使っている.

ピボット選択

ガウスの消去法における前進消去の式をもう一度みてみる.

ls_pivoting.eq1.gif

右辺の第2項にls_pivoting.eq2.gifで割るという計算が含まれている. 元々対角要素に0を含んでいる,もしくは,前進消去の過程で対角要素に0が出てくるとその時点で計算が止まってしまう. また,0でなくても絶対値が非常に小さい数だと桁落ちが起こり,誤差が増大する原因ともなる.

この問題に対処するために,ピボット(pivot:枢軸)選択(pivotal elimination method),もしくは,ピボッティング(pivoting) と呼ばれる方法を説明する. まず,連立方程式において式の順番は交換可能であることに注目する. つまり,ls_pivoting.eq3.gif行目の処理を行う際,その対角要素ls_pivoting.eq2.gifと同じ列でまだ処理していない要素ls_pivoting.eq4.gif の値を調べ,その絶対値が最大の行とls_pivoting.eq3.gif行を入れ替えた後で前進消去の式を適用する. これがピボット選択である. ちなみに,このように行の入れ替えだけ行うことを部分的ピボッティング(partial pivoting)と呼び, さらに列方向についても調べて入れ替えを適用することを完全ピボッティング(complete pivoting)と呼ぶ. ただし,部分的ピボッティングでは未知数の順番は変わらないが, 完全ピボッティングの場合は未知数の順番も変化するので注意.

以下にピボット選択を含めたガウス消去法のコード例を示す.

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
/*!
 * ピボット選択(Pivoting)
 *  - 行入れ替えだけの部分的ピボッティング
 * @param[inout] A n×nの係数項とn×1の定数項(b)を併せたn×(n+1)の行列
 * @param[in] n n元連立一次方程式
 * @param[in] k 対象行
 * @return 1:成功
 */
inline void Pivoting(vector< vector<double> > &A, int n, int k)
{
    // k行目以降でk列目の絶対値が最も大きい要素を持つ行を検索
    int p = k; // 絶対値が最大の行
    double am = fabs(A[k][k]);// 最大値
    for(int i = k+1; i < n; ++i){
        if(fabs(A[i][k]) > am){
            p = i;
            am = fabs(A[i][k]);
        }
    }
    // k != pならば行を交換(ピボット選択)
    if(k != p) swap(A[k], A[p]);
}
 
/*!
 * ガウスの消去法(行に関するピボット選択(部分ピボッティング)あり)
 * @param[inout] A n×nの係数項とn×1の定数項(b)を併せたn×(n+1)の行列.n+1列目に解が入る.
 * @param[in] n n元連立一次方程式
 * @return 1:成功
 */
int GaussEliminationWithPivoting(vector< vector<double> > &A, int n)
{
    // 前進消去(forward elimination)
    //  - 対角要素をのぞいた左下要素をすべて0にする(上三角行列にする)
    for(int k = 0; k < n-1; ++k){
        // ピボット選択
        Pivoting(A, n, k);
 
        double akk = A[k][k];
        for(int i = k+1; i < n; ++i){
            double aik = A[i][k];
            for(int j = k; j < n+1; ++j){
                A[i][j] = A[i][j]-aik*(A[k][j]/akk);
            }
        }
    }
 
    // 後退代入(back substitution)
    //  - x_nの解はb_n/a_nn,x_nをさらにn-1行の式に代入することでx_(n-1)を求める.
    //  - この作業を1行目まで続けることですべての解を得る.
    A[n-1][n] = A[n-1][n]/A[n-1][n-1];
    for(int i = n-2; i >= 0; --i){
        double ax = 0.0;
        for(int j = i+1; j < n; ++j){
            ax += A[i][j]*A[j][n];
        }
        A[i][n] = (A[i][n]-ax)/A[i][i];
    }
 
    return 1;
}

ガウス・ジョルダン法

ガウスの消去法の前進消去を拡張すると正方行列の逆行列を求めることができる. これは,ガウス・ジョルダン法(Gauss-Jordan method),もしくは,掃出し法(sweeping-out method)と呼ばれる.

まず,ls_gauss_jordan.eq1.gifの正方行列ls_gauss_jordan.eq2.gifを考える. ls_gauss_jordan.eq3.gifならばls_gauss_jordan.eq2.gifは正則行列(regular matrix)である (逆にls_gauss_jordan.eq4.gifなら特異行列(singular matrix)). ls_gauss_jordan.eq1.gifの単位行列をls_gauss_jordan.eq5.gifとすると,ls_gauss_jordan.eq2.gifが正則ならば,ls_gauss_jordan.eq6.gifを満足する行列ls_gauss_jordan.eq7.gifが1つだけ存在する. このls_gauss_jordan.eq7.gifを逆行列と呼びls_gauss_jordan.eq8.gifで表される.

さて,ガウスの消去法で対象とした連立1次方程式ls_gauss_jordan.eq9.gifにおいて,ls_gauss_jordan.eq10.gifls_gauss_jordan.eq7.gifに,ls_gauss_jordan.eq11.gifls_gauss_jordan.eq5.gifに置き換えることを考える. そうすると拡大行列は以下のようになる.

ls_gauss_jordan.eq12.gif

ls_gauss_jordan.eq6.gifであるので,両辺にls_gauss_jordan.eq8.gifをかけると

ls_gauss_jordan.eq13.gif

となる. つまり,拡大行列の左側半分が単位行列となるように変換することで,右側半分に逆行列ができる.

ls_gauss_jordan.eq14.gif

このとき,逆行列ls_gauss_jordan.eq8.gifは,

ls_gauss_jordan.eq15.gif

となる. 前進消去を拡張し,これらの処理を行うことで逆行列を求めるのが, ガウス・ジョルダン法である.

ガウス・ジョルダン法の実装

ガウスの消去法における前進消去では上三角行列にするために, ls_gauss_jordan.eq16.gifの条件を満たす要素のみ対象として処理していたが, ガウス・ジョルダン法ではすべての要素に対象を広げる. ガウス・ジョルダンでの前身消去の式は以下である.

ls_gauss_jordan.eq17.gif

ガウス・ジョルダン法で逆行列を求めるコード例を以下に示す.

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
/*!
 * ガウス・ジョルダン法(ピボット選択なし)
 * @param[inout] A n×2nの拡張行列
 * @param[in] n n元連立一次方程式
 * @return 1:成功
 */
int GaussJordan(vector< vector<double> > &A, int n)
{
    // 拡張行列の右半分を単位行列にする
    for(int i = 0; i < n; ++i){
        for(int j = 0; j < n; ++j){
            A[i][j+n] = (i == j ? 1.0 : 0.0);
        }
    }
 
    // ガウス・ジョルダン法(Gauss-Jordan method)で逆行列計算
    for(int k = 0; k < n; ++k){
        double akk = A[k][k];
        // 対角要素を1にするために,k行目のすべての要素をa_kkで割る
        for(int j = 0; j < 2*n; ++j){
            A[k][j] /= akk;
        }
 
        // k列目の非対角要素を0にする
        for(int i = 0; i < n; ++i){
            if(i == k) continue;
            double aik = A[i][k];
            for(int j = 0; j < 2*n; ++j){
                A[i][j] -= A[k][j]*aik;
            }
        }
    }
 
    return 1;
}

参考文献

  • 佐藤次男, 中村理一郎, “よくわかる数値計算 アルゴリズムと誤差解析の実際”, 日刊工業新聞社, 2001.
  • 川上一郎, “数値計算の基礎”, http://www7.ocn.ne.jp/~kawa1/

トップ   編集 凍結 差分 履歴 添付 複製 名前変更 リロード   新規 一覧 検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2012-07-12 (木) 10:58:36