第68回 プログラミングについて『連立１次方程式を解こう その2』
前回はガウス・ジョルダン法で行列を解いてみました。今回はもう少し行列の計算を考えてみましょう。まず次の連立方程式をみて下さい。
3x+2y+ z=1 ー（１）
y+2z=2 ー（２）
2z=3 ー（３）
この連立方程式のときには多分迷わずに、式（３）を解いて z を求め、この z の値を式（２）に代入して y を求め、最後に y と z を式（１）に代入して x を求めると思います。この方程式を行列で表現してみましょう。
┌ ┐┌ ┐ ┌ ┐
│3,2,1││x│ │1│
│0,1,2││y│＝│2│
│0,0,2││z│ │3│
└ ┘└ ┘ └ ┘
このような行列を三角行列と言い、上半分が０でないので上三角行列と呼びます。前回の４元の解を求めたプログラムを変形して、この解法で計算してみましょう。
#define DIM 4
double a[DIM][DIM+1]={
2.0, 5.0, 8.0, 3.0, 10.0,
4.0, 2.0, 3.0, 7.0, 25.0,
8.0, 6.0, 9.0, 4.0, 30.0,
9.0, 4.0, 3.0, 8.0, 45.0
};
main()
{
int i,j,k;
double t;
/* 上三角行列を作る */
for(i=0;i<DIM-1;i++){
for(j=i+1;j<DIM;j++){
t=a[j][i]/a[i][i];
for(k=0;k<=DIM;k++)a[j][k]-=a[i][k]*t;
}
disp(a);
getchar();
}
/* 対角以外の要素を０にする */
for(i=DIM-1;i>0;i--){
for(j=i-1;j>=0;j--){
t=a[j][i]/a[i][i];
for(k=0;k<=DIM;k++)a[j][k]-=a[i][k]*t;
}
disp(a);
getchar();
}
/* 対角の要素を１にする */
for(i=0;i<DIM;i++){
t=a[i][i];
for(j=0;j<=DIM;j++)a[i][j]/=t;
}
disp(a);
}
/* 行列の表示 */
disp(a)
double a[DIM][DIM+1];
{
int i,j;
for(i=0;i<DIM;i++){
for(j=0;j<=DIM;j++)printf("%10.4f ",a[i][j]);
printf("\n");
}
}
このプログラムを実行してみると、
2.0000 5.0000 8.0000 3.0000 10.0000
0.0000 -8.0000 -13.0000 1.0000 5.0000
0.0000 -14.0000 -23.0000 -8.0000 -10.0000
0.0000 -18.5000 -33.0000 -5.5000 0.0000
2.0000 5.0000 8.0000 3.0000 10.0000
0.0000 -8.0000 -13.0000 1.0000 5.0000
0.0000 0.0000 -0.2500 -9.7500 -18.7500
0.0000 0.0000 -2.9375 -7.8125 -11.5625
2.0000 5.0000 8.0000 3.0000 10.0000
0.0000 -8.0000 -13.0000 1.0000 5.0000
0.0000 0.0000 -0.2500 -9.7500 -18.7500
0.0000 0.0000 0.0000 106.7500 208.7500
ここで三角行列ができました。
2.0000 5.0000 8.0000 0.0000 4.1335
0.0000 -8.0000 -13.0000 0.0000 3.0445
0.0000 0.0000 -0.2500 0.0000 0.3162
0.0000 0.0000 0.0000 106.7500 208.7500
2.0000 5.0000 0.0000 0.0000 14.2506
0.0000 -8.0000 0.0000 -0.0000 -13.3958
0.0000 0.0000 -0.2500 0.0000 0.3162
0.0000 0.0000 0.0000 106.7500 208.7500
2.0000 0.0000 0.0000 -0.0000 5.8782
0.0000 -8.0000 0.0000 -0.0000 -13.3958
0.0000 0.0000 -0.2500 0.0000 0.3162
0.0000 0.0000 0.0000 106.7500 208.7500
ここで対角以外の要素が０になりました。
1.0000 0.0000 0.0000 -0.0000 2.9391
0.0000 1.0000 0.0000 0.0000 1.6745
0.0000 0.0000 1.0000 -0.0000 -1.2646
0.0000 0.0000 0.0000 1.0000 1.9555
対角の要素を１にして解が求まりました。
計算回数は、ｎが４のときには６０回、ｎが５のときには １２０回、ｎが１００のときには９９万９千９百回の計算回数になり、前回のアルゴリズムに比べ、だいぶ改善されています。この解法はガウス法というものです。
三角行列になったらこんなに計算回数が減るのであれば、行列を上三角行列と下三角行列の積に表すことができると、何かいいことがありそうです。
┌ ┐┌ ┐ ┌ ┐
│Ａ00、Ａ01、Ａ02││Ｘ0│ │Ｖ0│
│Ａ10、Ａ11、Ａ12││Ｘ1│＝│Ｖ1│
│Ａ20、Ａ21、Ａ22││Ｘ2│ │Ｖ2│
└ ┘└ ┘ └ ┘
例えば、上の行列を次のようにできたとします。
┌ ┐┌ ┐┌ ┐ ┌ ┐
│Ｌ00、 ０、 ０││Ｕ00、Ｕ01、Ｕ02││Ｘ0│ │Ｖ0│
│Ｌ10、Ｌ11、 ０││ ０、Ｕ11、Ｕ12││Ｘ1│＝│Ｖ1│
│Ｌ20、Ｌ21、Ｌ22││ ０、 ０、Ｕ22││Ｘ2│ │Ｖ2│
└ ┘└ ┘└ ┘ └ ┘
Ｌ Ｕ Ｘ Ｖ
この行列は、
Ｌ・Ｕ・Ｘ＝Ｖ
と論理的にあらわせますので、この式のＵ・ＸをＹと置き換えると、
Ｌ・Ｙ＝Ｖ ー（１）
Ｕ・Ｘ＝Ｙ ー（２）
となり、式１でＹを求め、式２でＸを求めることができます。ということは、ＬもＵも三角行列ですので、簡単に解を求めることができるという訳です。
問題はどうやって下三角行列と上三角行列に分解するかです。
┌ ┐ ┌ ┐┌ ┐
│Ａ00、Ａ01、Ａ02│ │Ｌ00、 ０、 ０││Ｕ00、Ｕ01、Ｕ02│
│Ａ10、Ａ11、Ａ12│＝│Ｌ10、Ｌ11、 ０││ ０、Ｕ11、Ｕ12│
│Ａ20、Ａ21、Ａ22│ │Ｌ20、Ｌ21、Ｌ22││ ０、 ０、Ｕ22│
└ ┘ └ ┘└ ┘
Ｌ Ｕ
なのですが、このままだとＬとＵは一意にきめることができません。そこで、
┌ ┐ ┌ ┐┌ ┐
│Ａ00、Ａ01、Ａ02│ │Ｌ00、 ０、 ０││ １、Ｕ01、Ｕ02│
│Ａ10、Ａ11、Ａ12│＝│Ｌ10、Ｌ11、 ０││ ０、 １、Ｕ12│
│Ａ20、Ａ21、Ａ22│ │Ｌ20、Ｌ21、Ｌ22││ ０、 ０、 １│
└ ┘ └ ┘└ ┘
Ｌ Ｕ
とＵの対角を１になるような下三角行列と上三角行列に分解することにします。すると、
Ａ00＝Ｌ00＊１
Ａ10＝Ｌ10＊１
Ａ20＝Ｌ20＊１
Ａ01＝Ｌ00＊Ｕ01
Ａ11＝Ｌ10＊Ｕ01＋Ｌ11
Ａ21＝Ｌ20＊Ｕ01＋Ｌ21
Ａ02＝Ｌ00＊Ｕ02
Ａ12＝Ｌ10＊Ｕ02＋Ｌ11＊Ｕ12
Ａ22＝Ｌ20＊Ｕ02＋Ｌ21＊Ｕ12+Ｌ22
ということになりますので、
Ｌ00＝Ａ00
Ｌ10＝Ａ10
Ｌ20＝Ａ20
Ｕ01＝Ａ01／Ｌ00
Ｌ11＝Ａ11ーＬ10＊Ｕ01
Ｌ21＝Ａ21ーＬ20＊Ｕ01
Ｕ02＝Ａ02／Ｌ00
Ｕ12＝（Ａ12ーＬ10＊Ｕ02）／Ｌ11
Ｌ22＝Ａ22ーＬ20＊Ｕ02ーＬ21＊Ｕ12
となります。ここで、Ｕ02を求める順番をＵ01の後にすると、
Ｌ00＝Ａ00
Ｌ10＝Ａ10
Ｌ20＝Ａ20
Ｕ01＝Ａ01／Ｌ00
Ｕ02＝Ａ02／Ｌ00
Ｌ11＝Ａ11ーＬ10＊Ｕ01
Ｌ21＝Ａ21ーＬ20＊Ｕ01
Ｕ12＝（Ａ12ーＬ10＊Ｕ02）／Ｌ11
Ｌ22＝Ａ22ーＬ20＊Ｕ02ーＬ21＊Ｕ12
となります。ということは、
┌ ┐ ┌ ┐
│ │ 、 ０ 、 ０ │ │１、（２）─→ │
│（１）、（３）、 ０ │ │０、 １、（４）│
│ ↓ 、 ↓ 、（５）│ │０、 ０、 １ │
└ ┘ └ ┘
Ｌ Ｕ
の順にＬとＵを交互に縦横縦横と求めていけばいいことになります。このように行列を下三角行列と上三角行列に分解することを、ＬＵ分解（多分ＬｏｗｅｒとＵｐｐｅｒの意味だと思います）といいます。またこの分解方法はクラウト法というものです。
ここでＬＵ分解を行なうプログラムを紹介したいのですが、もう１つプログラムを作る前に考えておくことがあります。行列を次のように分解するのは理論的にはそれでいいのですが、これをプログラムにすると、余分なメモリを用意しなければならず、特策とは言えないのです。
┌ ┐ ┌ ┐┌ ┐
│Ａ00、Ａ01、Ａ02│ │Ｌ00、 ０、 ０││ １、Ｕ01、Ｕ02│
│Ａ10、Ａ11、Ａ12│＝│Ｌ10、Ｌ11、 ０││ ０、 １、Ｕ12│
│Ａ20、Ａ21、Ａ22│ │Ｌ20、Ｌ21、Ｌ22││ ０、 ０、 １│
└ ┘ └ ┘└ ┘
Ｌ Ｕ
そこで、上の行列を見て気付くことは、ＬとＵの０以外で重なっているところは、Ｕの対角のところ、つまり１のところです。ということは、上の行列を、
┌ ┐ ┌ ┐
│Ａ00、Ａ01、Ａ02│ │Ｌ00、Ｕ01、Ｕ02│
│Ａ10、Ａ11、Ａ12│＝│Ｌ10、Ｌ11、Ｕ12│
│Ａ20、Ａ21、Ａ22│ │Ｌ20、Ｌ21、Ｌ22│
└ ┘ └ ┘
ＬＵ
とし、対角のＵの部分は暗黙的に１があるものとできることになります。このようにすることで元の行列以外に余分なメモリーを用意する必要がなくなります。
それでは、ＬＵ分解する部分を作ってみましょう。
int i,j,k;
for(i=0;i<DIM;i++){
for(j=i;j<DIM;j++){ /* Ｌを作る */
for(k=0;k<i;k++)a[j][i]-=(a[j][k]*a[k][i]);
}
for(j=i+1;j<DIM;j++){ /* Ｕを作る */
for(k=0;k<i;k++)a[i][j]-=(a[i][k]*a[k][j]);
a[i][j]/=a[i][i];
}
}
これだけのものになってしまうのです（実は理論をこのプログラムにするのには苦労しました）。定数 DIMは行列の次数を示すものです。これで元の行列がＬとＵを含んだ行列に分解されます。
それでは、ＬＵ分解によるプログラムの全文を作ってみましょう。
#define DIM 4
double a[DIM][DIM+1]={
2.0, 5.0, 8.0, 3.0, 10.0,
4.0, 2.0, 3.0, 7.0, 25.0,
8.0, 6.0, 9.0, 4.0, 30.0,
9.0, 4.0, 3.0, 8.0, 45.0
};
main()
{
int i,j,k;
double t;
disp(a);
getchar();
/* ＬＵ分解 */
for(i=0;i<DIM;i++){
for(j=i;j<DIM;j++){
for(k=0;k<i;k++)a[j][i]-=(a[j][k]*a[k][i]);
}
for(j=i+1;j<DIM;j++){
for(k=0;k<i;k++)a[i][j]-=(a[i][k]*a[k][j]);
a[i][j]/=a[i][i];
}
}
disp(a);
getchar();
/* 行列のＬ成分（下三角行列）から中間値を求める */
for(i=0;i<DIM;i++){
a[i][DIM]/=a[i][i];
a[i][i ] =1.0;
for(j=i+1;j<DIM;j++){
a[j][DIM]-=(a[i][DIM]*a[j][i]/a[i][i]);
a[j][i ]-=(a[i][i ]*a[j][i]/a[i][i]);
}
}
disp(a);
getchar();
/* 行列のＵ成分（上三角行列）から最終値を求める */
for(i=DIM-1;i>0;i--){
for(j=i-1;j>=0;j--){
a[j][DIM]-=(a[i][DIM]*a[j][i]/a[i][i]);
a[j][i ]-=(a[i][i ]*a[j][i]/a[i][i]);
}
}
disp(a);
getchar();
}
/* 行列の表示 */
disp(a)
double a[DIM][DIM+1];
{
int i,j;
for(i=0;i<DIM;i++){
for(j=0;j<=DIM;j++)printf("%10.4f ",a[i][j]);
printf("\n");
}
}
このプログラムを実行してみましょう。
2.0000 5.0000 8.0000 3.0000 10.0000
4.0000 2.0000 3.0000 7.0000 25.0000
8.0000 6.0000 9.0000 4.0000 30.0000
9.0000 4.0000 3.0000 8.0000 45.0000
最初の行列の内容。
2.0000 2.5000 4.0000 1.5000 10.0000
4.0000 -8.0000 1.6250 -0.1250 25.0000
8.0000 -14.0000 -0.2500 39.0000 30.0000
9.0000 -18.5000 -2.9375 106.7500 45.0000
ＬＵ分解した後の行列の状態。
1.0000 2.5000 4.0000 1.5000 5.0000
0.0000 1.0000 1.6250 -0.1250 -0.6250
0.0000 0.0000 1.0000 39.0000 75.0000
0.0000 0.0000 0.0000 1.0000 1.9555
下三角行列から中間値を求めた後の状態。 ↑
これが中間値
1.0000 0.0000 0.0000 0.0000 2.9391
0.0000 1.0000 0.0000 0.0000 1.6745
0.0000 0.0000 1.0000 0.0000 -1.2646
0.0000 0.0000 0.0000 1.0000 1.9555
上三角行列から最終値を求めた後の状態。
計算回数は、ｎが４のときには２６回、ｎが５のときには ５０回、ｎが１００のときには３３万８千３百回の計算回数になり、これまでのアルゴリズムに比べかなり改善されています。ＬＵ分解を行なう方法は計算速度が一般的に速いので、この方法が多く使われています。
連立方程式の求め方にはこれまでに説明してきたもの以外にも様々な方法がありますので、興味のある方は書籍などを参考にして下さい。
かなり長くなってしまいましたが、もう３つばかりお話することがあります。１つは、今回の一番最初の連立方程式が次のような順序だったときです。
2z=3
3x+2y+ z=1
y+2z=2
このときには、
┌ ┐┌ ┐ ┌ ┐
│0,0,2││x│ │3│
│3,2,1││y│＝│1│
│0,1,2││z│ │2│
└ ┘└ ┘ └ ┘
となってしまい、行列の対角に０が入ってしまいます。これでは計算が出来ませんので行を入れ替えて対角が０にならないようにする必要があります。また、これまでの説明では計算（除算をするとき）しようとする値が０のときの対処もしなければなりません。
次に、計算の精度の問題です。以前にも説明しましたが、実数計算を行なうときには大きな値と小さな値との演算を行なうと桁落ちをしてしまい、精度が低下してしまうのです。例えば、ＬＵ分解を行なうときにはＵ成分を生成するときに対角の値で除算を行ないますが、この対角の値が小さいと誤差を増幅してしまうことになりますので、ＬＵ分解を行なう前に出来る限り対角の値が大きくなるように行列の入れ替えを行なう必要があります。このような操作は、ピボッティング操作と呼ばれています。
最後に、これまでの説明で使用したプログラムの例は、配列の次元を固定としていました。あくまでもサンプルですのでこれで十分だったのですが、実際のプログラムでは次元を固定できないことがあります。ＦＯＲＴＲＡＮのように可変配列などを使用できる言語では結構簡単なのですが、Ｃ言語にはそのような機能はありません。ということは、行列の計算を行なうにはなんらかの工夫を行なう必要があるのです。工夫と言ってもたいしたことではないのですが、行列（配列）の要素のアドレスをコンパイラに任せるのではなく、プログラム側で計算することです。
例えば、N * M の要素をもる一次元の配列 A があったとき、この配列を N 行 M 列の配列として任意の要素のアドレスは次のようになります。
仮想要素 アドレス
[0][0] A
[0][1] A+sizeof(type A)
:
:
[0][j] A+sizeof(type A)*j
:
:
[1][j] A+sizeof(type A)*M + sizeof(type A)*j
[2][j] A+sizeof(type A)*M*2 + sizeof(type A)*j
:
:
[i][j] A+sizeof(type A)*M*i + sizeof(type A)*j
となりますので一般形としては i 行 j 列の要素のアドレスは、
A+sizeof(type A)*(M*i + j)
となります。ここで type A は配列 A の変数の型とします。
ということで、今回はこれぐらいにしましょう。行列の計算はときには必要になることがありますので、ある程度の知識を持っておくことをお勧めします。
ではまた次回。
