ピボット選択

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

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

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

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

/*!
 * ピボット選択(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;
}

Menu