「アルゴリズムノート」ガウス消元

リード
ガウス消元(gauss elimination)(g a u s s e l i m i n a t i o n)は線形代数のアルゴリズムである.線形方程式群を解くために用いることができ,行列のランクを求めることができ,可逆方程式の逆行列を求めることができる.ガウス消元法の原理は,初等行変換で広がり行列を行階段行列に変換し,その後方程式の解を求めることである.
例題
n n元線形方程式のグループを解く.
ガイド人
二元一次方程式のグループをどのように解くか考えてみましょう.
{ax+by=c (1)dx+ey=f (2) { d x + e y = f   ( 2 ) a x + b y = c   ( 1 )
a,b,c,d,e,f>0 a,b,c,d,e,f>0およびae≠bd a e≠b dを保証する
       (1)−ad×(2)               ( 1 ) − a d × ( 2 )
=>(b−aed) y=c−afd => ( b − a e d )   y = c − a f d
=>y=cd−afbd−ae => y = c d − a f b d − a e
(1)(1)を代入してx=ce−bfbd−ae x=c e−bfbd−aeを得る
消元して持ち込みますよね?
プログラミングシミュレーションを使用します.
アルゴリズムフロー
方程式グループは、a 1,1 x 1+a 1,2 x 2+...+a1,nxn=b1 a 1 , 1 x 1 + a 1 , 2 x 2 + . . . + a 1 , n x n = b 1 a2,1x1+a2,2x2+...+a2,nxn=b2 a 2 , 1 x 1 + a 2 , 2 x 2 + . . . + a 2 , n x n = b 2 ... . . . an,1x1+an,2x2+...+an,nxn=bn a n , 1 x 1 + a n , 2 x 2 + . . . + a n , n x n = b n

i回目の消元を行う場合は、まず任意の行を取り出し、row r o wと記す.

その後、

は、取り出されていない全ての行jについて消元する.具体的なプロセスは次のとおりです.

はd=aj,iarow,i d=a j,i a r o w,iを命令する.

bj←bj−d×bi b j ← b j − d × b i .

aj,k←aj,k−d×arow,k a j , k ← a j , k − d × a r o w , k .

消元が完了するまでxn,n−1,...,1 x n , n − 1 , . . . , 1を順に代入すると,方程式の解が求められる.アルゴリズム時間複雑度O(n 3)O(n 3)

最適化
|arow,i||a r o w,i|を最大にするrow r o wを探すたびに、精度誤差を効果的に低減することができる.
コード#コード#

//        n       
#include 
#include 
#include 
const double eps = 1e-18;
const int maxn = 105;
using namespace std;
typedef double db;
typedef long double ldb;
db t;
int n;
ldb a[maxn][maxn], b[maxn], x[maxn];
int main() {
    scanf("%d", &n);
    for (int i = 1; i <= n; i++) {
        for (int j = 1; j <= n; j++) {
            scanf("%lf", &t);
            a[i][j] = t;
        }
        scanf("%lf", &t);
        b[i] = t;
    }
    for (int i = 1; i <= n; i++) {
        int ind = i;
        for (int j = i + 1; j <= n; j++) {
            if (abs(a[j][i]) > abs(a[ind][i])) {
                ind = i;
            }
        }
        if (abs(a[ind][i]) < eps) {
            puts("-1");
            return 0;
        }
        if (ind != i) {
            swap(b[i], b[ind]);
            for (int j = 1; j <= n; j++) {
                swap(a[i][j], a[ind][j]);
            }
        }
        for (int j = i + 1; j <= n; j++) {
            double d = a[j][i] / a[i][i];
            b[j] -= b[i] * d;
            for (int k = i; k <= n; k++) {
                a[j][k] -= a[i][k] * d;
            }
        }
        /*
        for (int j = 1; j <= n; j++) {
            for (int k = 1; k <= n; k++) {
                printf("%lf ", a[j][k]);
            }
            printf("%lf
", b[j]);
        }
        */
    }
    for (int i = n; i > 0; i--) {
        for (int j = i + 1; j <= n; j++) {
            b[i] -= x[j] * a[i][j];
        }
        x[i] = b[i] / a[i][i];
    }
    for (int i = 1; i <= n; i++) {
        t = x[i];
        printf("x[%d] = %.6lf;
", i, t);
    }
    return 0;
}

練習する
nを求める×n n × nの行列の行列式.
ヒント
行列を上三角行列に変換すればよい.
コード#コード#

//            
#include 
#include 
#include 
const int maxn = 105;
const double eps = 1e-6;
typedef long double ldb;
using namespace std;
double t;
int n;
ldb a[maxn][maxn];
double read() {
    scanf("%lf", &t);
    return t;
}
void write(ldb x) {
    t = x;
    printf("%lf;
", t);
}
int main() {
    scanf("%d", &n);
    for (int i = 1; i <= n; i++) {
        for (int j = 1; j <= n; j++) {
            a[i][j] = read();
        }
    }
    ldb x = 1;
    for (int i = 1; i <= n; i++) {
        //          ，      
        int ind = i;
        for (int j = i + 1; j <= n; j++) {
            if (abs(a[ind][i]) < abs(a[j][i])) {
                ind = j;
            }
        }
        if (ind != i) {
            x = -x;
            for (int j = 1; j <= n; j++) {
                swap(a[i][j], a[ind][j]);
            }
        }
        if (abs(a[i][i]) < eps) {
            puts("0");
            return 0;
        }
        for (int j = i + 1; j <= n; j++) {
            ldb d = a[j][i] / a[i][i];
            for (int k = i; k <= n; k++) {
                a[j][k] -= d * a[i][k];
            }
        }
        printf("Step #%d of Gaussian Elimination
", i);
        for (int j = 1; j <= n; j++) {
            for (int k = 1; k <= n; k++) {
                t = a[j][k];
                printf("%lf%c", t, k == n ? '
' : ' ');
            }
        }
    }
    for (int i = 1; i <= n; i++) {
        x *= a[i][i];
    }
    printf("Det(A) = ");
    write(x);
    return 0;
}

ガウス消元とダイナミックプランニング
例題
【LightOJ1511】Snakes and Ladders
問題解
ヘビがいなければ確率動的計画には後効性がなく,状態遷移方程式を押すとよい.ヘビがいる場合,dp d p値を未知数と見なし,n n元の線形方程式群をリストし,Gauss消元法で解くとよい.
コード#コード#

#include 
#include 
#include 
#include 
typedef long double ldb;
const ldb eps = 1e-6;
const int maxn = 105;
const int n = 100;
using namespace std;
int tasks, m, to[maxn];
ldb a[maxn][maxn], b[maxn], x[maxn];
int main() {
    scanf("%d", &tasks);
    for (int task = 1; task <= tasks; task++) {
        memset(to, 0, sizeof(to));
        for (int i = 1; i <= n; i++) {
            for (int j = 1; j <= n; j++) {
                a[i][j] = 0;
            }
            b[i] = 0, x[i] = 0;
        }
        scanf("%d", &m);
        for (int u, v, i = 1; i <= m; i++) {
            scanf("%d %d", &u, &v);
            to[u] = v;
        }
        for (int i = 1; i <= n; i++) {
            if (i == n) {
                a[i][i] = 1;
            } else if (to[i]) {
                a[i][i] = 1, a[i][to[i]] = -1;
            } else if (i + 6 <= 100) {
                for (int j = i + 1; j <= i + 6; j++) {
                    a[i][j] = -1;
                }
                a[i][i] = 6, b[i] = 6;
            } else {
                a[i][i] = 100 - i;
                for (int j = i + 1; j <= n; j++) {
                    a[i][j] = -1;
                }
                b[i] = 6;
            }
        }
        for (int i = 1; i <= n; i++) {
            int ind = i;
            for (int j = i + 1; j <= n; j++) {
                if (abs(a[j][i]) > abs(a[ind][i])) {
                    ind = j;
                }
            }
            if (i != ind) {
                swap(b[i], b[ind]);
                for (int j = 1; j <= n; j++) {
                    swap(a[i][j], a[ind][j]);
                }
            }
            for (int j = i + 1; j <= n; j++) {
                ldb d = a[j][i] / a[i][i];
                b[j] -= b[i] * d;
                for (int k = 1; k <= n; k++) {
                    a[j][k] -= a[i][k] * d;
                }
            }
        }
        double tmp;
        for (int i = n; i >= 1; i--) {
            for (int j = i + 1; j <= n; j++) {
                b[i] -= a[i][j] * x[j];
            }
            x[i] = b[i] / a[i][i];
            tmp = x[i];
        }
        tmp = x[1];
        printf("Case %d: %.10lf
", task, tmp);
    }
    return 0;
}
/*
100
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
*/

P.S.
実は、このタイプのテーマは反復的な方法で作ることもでき、最後にはテーマの要求の精度を達成することも多い.
まとめ
いくつかの後方性のある確率動的計画(n≦500)(n≦500)に遭遇すると,まずGauss消元を考えるべきである.これは情報学コンテストでよく使われる手段である.