ガウスの消去法†
ガウスの消去法(Gaussian elimination)と呼ばれる方法では，
前進消去と後退代入という2段階の計算により解を求める．
一般的な(中学校で習う)解き方の一つに加減法というものがある．
たとえば，以下のような連立方程式で考える．
加減法では片方の式の両辺に実数をかけて，係数を合わせて，加減算する．
ここでは下の式の両辺を2倍して，上の式から引くことにする．
あとは，
の結果を1番目の式に代入することで，
を得る．
これを行列で表すと，
に対して，2行目に
をかけて引くことで，
という形にする．そして，一番下の行から未知数を解いて代入していく．
前半の処理では係数行列
の対角要素を除いた左下部分がすべて0となる
以下に示すような上三角行列(upper triangle matrix)を求めている．
この処理のことを前進消去(forward elimination)と呼ぶ．
前進消去により，未知数
は
で求められ，
このを
行目に代入することで
が得られ，
を
行目に代入することで
が得られる．
これを
まで繰り返すことですべての解を得る．
この後半部分の処理のことを後退代入(back substitution)と呼ぶ．
前進消去と後退代入により連立1次方程式の解を求める方法をガウスの消去法という．
ガウスの消去法の実装†
ガウスの消去法を
元連立1次方程式に適用する．
まず，に定数ベクトル
を列の最後に追加した
の行列を考える．
この行列は拡大行列と呼ばれる．
拡大行列の要素を
とする．ここで，C,C++で実装することを考えて要素番号を0始まりにしている．
今後，実装の説明では0で始まるインデックスを用いるので注意．
前進消去は次の式で表される．
ここで，
上式により
を
で更新することで上三角行列を得る．
この処理のC++のコードを以下に示す．
1
2
3
4
5
6
7
8
9
10
11
| | 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){ A[i][j] = A[i][j]-aik*(A[k][j]/akk);
}
}
}
|
2次元配列A[n][n+1]には係数行列と定数ベクトルが格納されている．
次に後退代入の式を以下に示す．
後退代入のC++のコード例を以下に示す．
1
2
3
4
5
6
7
8
9
10
11
| | 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
| |
int GaussElimination(vector< vector<double> > &A, int n)
{
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){ A[i][j] = A[i][j]-aik*(A[k][j]/akk);
}
}
}
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を使っている．
ピボット選択†
ガウスの消去法における前進消去の式をもう一度みてみる．
右辺の第2項に
で割るという計算が含まれている．
元々対角要素に0を含んでいる，もしくは，前進消去の過程で対角要素に0が出てくるとその時点で計算が止まってしまう．
また，0でなくても絶対値が非常に小さい数だと桁落ちが起こり，誤差が増大する原因ともなる．
この問題に対処するために，ピボット(pivot:枢軸)選択(pivotal elimination method)，もしくは，ピボッティング(pivoting)
と呼ばれる方法を説明する．
まず，連立方程式において式の順番は交換可能であることに注目する．
つまり，
行目の処理を行う際，その対角要素と同じ列でまだ処理していない要素
の値を調べ，その絶対値が最大の行と行を入れ替えた後で前進消去の式を適用する．
これがピボット選択である．
ちなみに，このように行の入れ替えだけ行うことを部分的ピボッティング(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
| |
inline void Pivoting(vector< vector<double> > &A, int n, int 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]);
}
}
if(k != p) swap(A[k], A[p]);
}
int GaussEliminationWithPivoting(vector< vector<double> > &A, int n)
{
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);
}
}
}
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)と呼ばれる．
まず，
の正方行列
を考える．
ならばは正則行列(regular matrix)である
(逆に
なら特異行列(singular matrix))．
の単位行列を
とすると，が正則ならば，
を満足する行列
が1つだけ存在する．
このを逆行列と呼び
で表される．
さて，ガウスの消去法で対象とした連立1次方程式
において，
をに，
をに置き換えることを考える．
そうすると拡大行列は以下のようになる．
であるので，両辺にをかけると
となる．
つまり，拡大行列の左側半分が単位行列となるように変換することで，右側半分に逆行列ができる．
このとき，逆行列は，
となる．
前進消去を拡張し，これらの処理を行うことで逆行列を求めるのが，
ガウス・ジョルダン法である．
ガウス・ジョルダン法の実装†
ガウスの消去法における前進消去では上三角行列にするために，
の条件を満たす要素のみ対象として処理していたが，
ガウス・ジョルダン法ではすべての要素に対象を広げる．
ガウス・ジョルダンでの前身消去の式は以下である．
ガウス・ジョルダン法で逆行列を求めるコード例を以下に示す．
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
| |
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);
}
}
for(int k = 0; k < n; ++k){
double akk = A[k][k];
for(int j = 0; j < 2*n; ++j){
A[k][j] /= akk;
}
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;
}
|
参考文献†
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
