LU分解を使ってn元連立1次方程式を解く

LU分解を使ってn元連立1次方程式を解くプログラムをC++で作成した。

LU分解とは、行列を下三角行列L(Lower triangle matrix)と上三角行列U(Upper triangle matrix)の積LUに分解することである。

例えば

$A=\begin{pmatrix}1&2\\3&4\end{pmatrix}$

をLU分解すると

$L=\begin{pmatrix}1&0\\3&1\end{pmatrix}$
$U=\begin{pmatrix}1&2\\0&-2\end{pmatrix}$

となる。（これはあくまで一例であり、この場合A=LUとなるL,Uは無限に存在する）

このLU分解を使ってn元連立1次方程式を解く。

n元連立1次方程式を行列を使って表せば

$A\bf{x}=\bf{b}$ ( $A$ はn次の正方行列 $\bf{x}$ , $\bf{b}$ はn次の列ベクトル) … (1)

となる。

このとき $A$ を

$A=LU$ ( $L$ , $U$ はn次の正方行列)

とLU分解する。

すると、(1)の方程式は

$LU\bf{x}=\bf{b}$ … (2)

と書き換えられる。

ここで

$\bf{y}=U\bf{x}$

とおけば、(2)より、

$L\bf{y}=\bf{b}$

となる。

このとき、 $\bf{L}$ は下三角行列なので前進代入により、簡単に $\bf{y}$ が求まる。
さらに、 $\bf{U}$ は上三角行列なので後退代入により、こちらも簡単に $\bf{x}$ が求まる。

このときの計算量のオーダは $O(n^2)$ である。n元連立1次方程式の解法として有名なガウスの消去法のオーダは $O(n^3)$ なので、LU分解を用いた方が高速である（オーダだけで計算量を完全に比較できるわけではないが、ここでは簡易のためにそうしておく）。

しかし、そのオーダはLU分解に必要なコストを除いたものである。一般に、LU分解をするのに必要な計算量のオーダはガウスの消去法と同じく $O(n^3)$ なので、結果として同じオーダになってしまう。これではLU分解を用いる意味が無いように思える。

もちろん、そのような事は無い。

(1)のような方程式が複数与えられ、その $A$ がすべて同じ値だったとする。
例えばこんな問題である。

$\begin{pmatrix}1&2\\3&4\end{pmatrix}\begin{pmatrix}x_1\\x_2\end{pmatrix}=\begin{pmatrix}1\\4\end{pmatrix}$
$\begin{pmatrix}1&2\\3&4\end{pmatrix}\begin{pmatrix}x_1\\x_2\end{pmatrix}=\begin{pmatrix}5\\6\end{pmatrix}$
$\begin{pmatrix}1&2\\3&4\end{pmatrix}\begin{pmatrix}x_1\\x_2\end{pmatrix}=\begin{pmatrix}9\\2\end{pmatrix}$

この場合、ガウスの消去法で解くならば全ての問題に対し $O(n^3)$ のオーダの計算量がかかる。
それに対し、LU分解後の行列 $\bf{LU}$ は $\bf{A}$ に対して決まるため、一度計算してしまえば以降同じ $\bf{LU}$ を使いまわせる。そのため、最初の一問以外は $O(n^2)$ で済む。

これがLU分解を用いる利点である。

LU分解の具体的な方法については省略する（下記のプログラムに書いてはあるが）。

プログラムと問題

プログラム

main.cpp

#include <iostream>

using namespace std;

void LUdecomposition(int dimension, double** A, double** L, double** U){//LU分解
	for(int i=0;i<dimension;i++){
		for(int j=0;j<dimension;j++){
			L[i][j]=0.0;
			U[i][j]=0.0;
			if(i==j)L[i][j]=1.0;
		}
	}
	double sum;
	for(int i=0;i<dimension;i++){
		for(int j=0;j<dimension;j++){
			if(i>j){
				/*L成分を求める*/
				sum=0.0;
				for(int k=0;k<j;k++){
					sum+=L[i][k]*U[k][j];
				}
				L[i][j]=(A[i][j]-sum)/U[j][j];
			}else{
				/*U成分を求める*/
				sum=0.0;
				for(int k=0;k<i;k++){
					sum+=L[i][k]*U[k][j];
				}
				U[i][j]=A[i][j]-sum;
			}
		}
	}
}

void Lforwardsubs(int dimension, double** L, double* b, double* y){//前進代入
	for(int i=0;i<dimension;i++){
		y[i]=b[i];
	}
	for(int i=0;i<dimension;i++){
		y[i]/=L[i][i];
		for(int j=i+1;j<dimension;j++){
			y[j]-=y[i]*L[j][i];
		}
	}
}

void Ubackwardsubs(int dimension, double** U, double* y, double* x){//後退代入
	for(int i=0;i<dimension;i++){
		x[i]=y[i];
	}
	for(int i=dimension-1;i>=0;i--){
		x[i]/=U[i][i];
		for(int j=i-1;j>=0;j--){
			x[j]-=x[i]*U[j][i];
		}
	}
}

int main(){
	int dimension,problem_num;
	double **A,**L,**U,**b,**x,*y;

/*入力を受け取る*/
	cin >> dimension;
	A=new double*[dimension];
	L=new double*[dimension];
	U=new double*[dimension];
	y=new double[dimension];
	for(int i=0;i<dimension;i++){
		A[i]=new double[dimension];
		L[i]=new double[dimension];
		U[i]=new double[dimension];
	}

	for(int i=0;i<dimension;i++){
		for(int j=0;j<dimension;j++){
			cin >> A[i][j];
		}
	}

	cin >> problem_num;
	b=new double*[problem_num];
	x=new double*[problem_num];
	for(int i=0;i<problem_num;i++){
		b[i]=new double[dimension];
		x[i]=new double[dimension];
	}

	for(int i=0;i<problem_num;i++){
		for(int j=0;j<dimension;j++){
			cin >> b[i][j];
		}
	}

/*問題を解く*/
	LUdecomposition(dimension,A,L,U);
	for(int i=0;i<problem_num;i++){
		Lforwardsubs(dimension,L,b[i],y);
		Ubackwardsubs(dimension,U,y,x[i]);
	}

/*解答を出力する*/
	for(int i=0;i<problem_num;i++){
		for(int j=0;j<dimension;j++){
			if(j!=0)cout << " "; 
			cout << x[i][j];
		}
		cout << endl;
	}

/*メモリの解放*/
	for(int i=0;i<dimension;i++){
		delete [] A[i];
		delete [] L[i];
		delete [] U[i];
	}
	for(int i=0;i<problem_num;i++){
		delete [] b[i];
		delete [] x[i];
	}
	delete [] A;
	delete [] L;
	delete [] U;
	delete [] b;
	delete [] x;
	delete [] y;
}

プログラムへの入力は

方程式の次数
係数A(次数の2乗だけ入力)
bの個数
b(上で入力した数だけ入力)

という書式で行う。

下の入力ファイル例を見ればわかるだろう。

入力ファイル(例)

problem.txt

出力ファイル

result.txt

※7/3 少しだけ手直し