ガウスの消去法 - わかばテクノロジ

線形代数
ガウスの消去法; Gaussian Eliminate

概要

ガウスの消去法は連立方程式の解、逆行列の計算などに用いられるアルゴリズムである。
このアルゴリズムは前進消去、後退代入という2ステップからなり、
前進消去は行の入れ替え(ピボット選択と入れ替え)と、ある行に別の行を何倍にして変数の一部の係数を0にする(消去)からなる。
後退代入はある行の下の行を用い、対角成分だけを1、他の成分を0にする(単位行列)にする処理である。

行列クラス本体はこちら

疑似逆行列の定義

ガウスの消去法では逆行列または疑似逆行列を計算するために用いられる。
\(A\)の疑似逆行列\(A^{+}\)は以下を満たす。ここで\(A^{*}\)を\(A\)の随伴行列(実数を成分とする行列では転置行列と等しい)とする。
・\(A\)と\(A^{+}\)は互いに広義可逆元
\( \quad A A^{+} A = A, \ A^{+} A A^{+} = A^{+} \)
・\(A A^{+}\)と\(A^{+} A\)はエルミート行列
\( \quad (A A^{+})^{*} = A A^{+}, \ (A^{+} A)^{*} = A^{+} A \)

疑似逆行列の性質

疑似逆行列は以下のような性質を持つ
\( \quad A = (A^{+})^{+} \)
\( \quad (A^T)^{+} = (A^{+})^T \)

ソースコード


namespace Algebra {
    /// <summary>行列クラス</summary>
    public partial class Matrix {
        /// <summary>ガウスの消去法</summary>
        public static void GaussianEliminate(Matrix m, ref Matrix v) {
            if(!IsSquare((Matrix)m) || m.Rows != v.Rows) {
                throw new ArgumentException();
            }

            int i, j, k, p, n = m.Rows;
            double pivot, inv_mii, mul, swap;

            for(i = 0; i < n; i++) {
                pivot = Math.Abs(m[i, i]);
                p = i;

                //ピボット選択
                for(j = i + 1; j < n; j++) {
                    if(Math.Abs(m[i, j]) > pivot) {
                        pivot = Math.Abs(m[i, j]);
                        p = j;
                    }
                }

                //ピボットが閾値以下ならばMは正則行列でないので逆行列は存在しない
                if(pivot < 1.0e-12) {
                    v = Invalid(v.Rows, v.Columns);
                    return;
                }

                //行入れ替え
                if(p != i) {
                    for(j = 0; j < n; j++) {
                        swap = m[i, j];
                        m[i, j] = m[p, j];
                        m[p, j] = swap;
                    }

                    for(j = 0; j < v.Columns; j++) {
                        swap = v[i, j];
                        v[i, j] = v[p, j];
                        v[p, j] = swap;
                    }
                }

                // 前進消去
                inv_mii = 1 / m[i, i];
                m[i, i] = 1;
                for(j = i + 1; j < n; j++) {
                    m[i, j] *= inv_mii;
                }
                for(j = 0; j < v.Columns; j++) {
                    v[i, j] *= inv_mii;
                }

                for(j = i + 1; j < n; j++) {
                    mul = m[j, i];
                    m[j, i] = 0;
                    for(k = i + 1; k < n; k++) {
                        m[j, k] -= m[i, k] * mul;
                    }
                    for(k = 0; k < v.Columns; k++) {
                        v[j, k] -= v[i, k] * mul;
                    }
                }
            }

            // 後退代入
            for(i = n - 1; i >= 0; i--) {
                for(j = i - 1; j >= 0; j--) {
                    mul = m[j, i];
                    for(k = i; k < n; k++) {
                        m[j, k] = 0;
                    }
                    for(k = 0; k < v.Columns; k++) {
                        v[j, k] -= v[i, k] * mul;
                    }
                }
            }
        }
    }
}