LU分解による連立方程式の解法
LU分解による連立方程式の解法はコンピュータを用いた数値計算で一般的に利用されている手法である。
条件が整えば計算量が抑えられるなど効率的なアルゴリズムとして知られている。
LU分解による連立方程式の解の導出には以下の三つのステップにより実行される
LU分解で行列を分解した後の処理はガウスの消去法に順ずる処理になっている。
LU分解
LU分解とは与えられた行列\(A\)に対して、\(A\)を下三角行列\(L\)と上三角行列\(U\)の積の形式になる様に分解することである。
\[
A=LU
\]
\[A=
\left[
\begin{array}{}
a_{11} &a_{12} & \cdots &a_{1n} \\
a_{21} &a_{22} & \cdots &a_{2n} \\
\vdots & \vdots & \ddots & \vdots \\
a_{n1} &a_{n2} & \cdots &a_{nn} \\
\end{array}
\right]
\Rightarrow
L=
\left[
\begin{array}{}
a_{11}' & 0 & \cdots &0 \\
a_{21}' &a_{22}' & \cdots &0 \\
\vdots & \vdots & \ddots & \vdots \\
a_{n1}' &a_{n2}' & \cdots &a_{nn}' \\
\end{array}
\right]
,
U=
\left[
\begin{array}{}
a_{11}' &a_{12}' & \cdots &a_{1n}' \\
0 &a_{22}' & \cdots &a_{2n}' \\
\vdots & \vdots & \ddots & \vdots \\
0 & 0 & \cdots &a_{nn}' \\
\end{array}
\right]
\]
これにより連立方程式を含め、再利用可能な形式となり繰り返し計算を行う上で効率が上がったり、簡単に行列式や逆行列を求められるなど多くのメリットがある。
LU分解の仕方については様々なものがあるが、連立方程式への応用については上三角行列\(U\)の対角項を全て１にする様に分解する。
具体的な例をあげると次の行列は以下の様に分解される
\[A=
\left[
\begin{array}{}
1 &2 & 4 \\
2 &7 & 23 \\
4 &13 & 47 \\
\end{array}
\right]
\Rightarrow
L=
\left[
\begin{array}{}
1 &0 & 0 \\
2 &3 & 0 \\
4 &5 & 6 \\
\end{array}
\right]
,
U=
\left[
\begin{array}{}
1 &2 & 4 \\
0 &1 & 5 \\
0 &0 & 1 \\
\end{array}
\right]
\]
以下から具体的なLU分解法について述べる。
次の3元1次連立方程式をLU分解で解く場合、まず係数行列のみに着目する。
\[
\left\{
\begin{array}{1}
x_{1}+2x_{2}+4x_{3}=11\\
2x_{1}+7x_{2}+23x_{3}=43\\
4x_{1}+13x_{2}+47x_{3}=85\\
\end{array}
\right.
\]
連立方程式の係数行列を\(A\)とすると
\[
A=
\left[
\begin{array}{}
1 &2 & 4 \\
2 &7 & 23 \\
4 &13 & 47 \\
\end{array}
\right]
\]
まず\(L\)の１列目に\(A\)の対角項と垂直成分を代入する。
\[
L=
\left[
\begin{array}{}
1 &0 & 0 \\
2 &0 & 0 \\
4 &0 & 0 \\
\end{array}
\right]
\]
次に\(U\)１行目に対角項と水平成分両方の対角項で割った値を代入する（この場合係数は１なのでそのまま）
\[
U=
\left[
\begin{array}{}
1 &2 & 4 \\
0 &0 & 0 \\
0 &0 & 0 \\
\end{array}
\right]
\]
対角項の斜め下の小行列から対角項を除いた\(L\)の垂直成分ベクトルと\(U\)の水平成分ベクトルの積を引いたものをとあらたな\(A’\)行列とする。
\[
A'=
\left[
\begin{array}{}
7 &23 \\
13 &47 \\
\end{array}
\right]
-
\left[
\begin{array}{}
2 \\
4 \\
\end{array}
\right]
\left[
\begin{array}{}
2 &4 \\
\end{array}
\right]
=
\left[
\begin{array}{}
7 &23 \\
13 &47 \\
\end{array}
\right]
-
\left[
\begin{array}{}
4 &8 \\
8 &16 \\
\end{array}
\right]
=
\left[
\begin{array}{}
3 &15 \\
5 &31 \\
\end{array}
\right]
\]
\(L\)に行列\(A’\)の対角項と垂直成分を代入する、\(U\)に行列\(A’\)の対角項と垂直成分を対角項で割った値を代入する。
\[
L=
\left[
\begin{array}{}
1 &0 & 0 \\
2 &3 & 0 \\
4 &5 & 0 \\
\end{array}
\right]
U=
\left[
\begin{array}{}
1 &2 & 4 \\
0 &1 & 5 \\
0 &0 & 0 \\
\end{array}
\right]
\]
\(A'\)の小行列（この場合は値）から水平成分と垂直成分の積（同じく値どうしの積）を引いたものを\(A''\)とする（この行列も値）
\[
A''=
31-(5 \times 5)=6
\]
\(A''\)の値そのものを\(L\)に代入して、さらに値を値で割った値、１を\(U\)に代入するとLU分解が完成する。
\[
L=
\left[
\begin{array}{}
1 &0 & 0 \\
2 &3 & 0 \\
4 &5 & 6 \\
\end{array}
\right]
U=
\left[
\begin{array}{}
1 &2 & 4 \\
0 &1 & 5 \\
0 &0 & 1 \\
\end{array}
\right]
\]
前進代入
ここからはLU分解をした後の解の導出について考える。
与えられた方程式を\(A\)を係数行列、\(x\)を未知数のベクトル、\(b\)を定数項ベクトルとして
\[A\mathbf{x}=\mathbf{b}\]
と置く。AはLU分解によってLとUの積に分解されたので
\[LU\mathbf{x}=\mathbf{b}\]
ここで、未知数ベクトル\(y\)を導入して形式的に
\[U\mathbf{x}=\mathbf{y}\]
と置くと
\[L\mathbf{y}=\mathbf{b}\]
となり\(\mathbf{y}\)の方程式となる。
\(L\)と\(\mathbf{b}\)は既知なので\(\mathbf{y}\)について解け、
\(\mathbf{y}\)が解けると上式(\(U\mathbf{x}=\mathbf{y}\))より\(\mathbf{x}\)が求まる。
ここでyを求める過程で前進代入（前進消去ではないことに注意）、xを求める過程で後退代入を利用する。
まずは
\[L\mathbf{y}=\mathbf{b}\]
から解く。
前進代入の式は以下の通りである
\[
y_i=\frac{b_i-\sum^{i−1}_{k=1}l_{ik}y_k}{l_{ii}}\\
(i=1,2,\dots n)
\]
既知となった\(y_{1}\)から順に\(y_{2},y_{3}\)と代入しながら全ての\(y_i\)を得る。
具体的に求めると
\[
\left[
\begin{array}{ccc|c}
1 & 0 & 0 &11\\
2 & 3 & 0 &43\\
4 & 5 & 6 &85\\
\end{array}
\right]
\]
これを前進代入で解いて
\[
y_{1}=11 \\
y_{2}=\frac{43-2\times11}{3}=7\\
y_{3}=\frac{85-4\times11-5\times7}{6}=1
\]
\[
y_{1}=11, y_{2}=7, y_{3}=1
\]
後退代入
前進代入で求めたyをつかって
\[U\mathbf{x}=\mathbf{y}\]
を解く。後退代入の式は以下の通りである。
\[
y_i=y_i-\sum^{n}_{k=i+1}u_{ik}y_k\\
(i=n,n-1,\dots 1)
\]
既知となった\(x_{n}\)から順に\(x_{n-1},x_{n-2}\)と代入しながら全ての\(y_i\)を得る。
具体的に求めると
\[
\left[
\begin{array}{ccc|c}
1 & 2 & 4 &11\\
0 & 1 & 5 &7\\
0 & 0 & 1 &1\\
\end{array}
\right]
\]
これを解いて
\[
x_{3}=1 \\
x_{2}=7-5\times1=2\\
x_{1}=11-2\times2-4\times1=3
\]
\[\therefore x_{1}=3, x_{2}=2, x_{3}=1 \]
以上がLU分解による連立方程式の解法である。
コード
public double[] LUDecomposition(double[][] a) {
int n=a.length;
double[][] l=new double[n][n];
double[][] u=new double[n][n];
//LU分解
for(int k=0;k<n;k++) {
double alpha=1.0/a[k][k];
for(int i=k;i<n;i++) {
l[i][k]=a[i][k];
u[k][i]=a[k][i]*alpha;
}
for(int i=k+1;i<n;i++) {
for(int j=k+1;j<n;j++) {
a[i][j]=a[i][j]-(l[i][k]*u[k][j]);
}
}
}
//前進代入
double[] y=new double[n];
for(int i=0;i<n;i++) {
y[i]=a[i][n];
for(int j=0;j<i;j++) {
y[i]-=l[i][j]*y[j];
}
y[i]/=l[i][i];
}
//後退代入
double[] x=new double[n];
for(int i=n-1;i>=0;i--) {
x[i]=y[i];
for(int j=n-1;j>=i+1;j--) {
x[i]-=u[i][j]*x[j];
}
}
return x;
}
サンプルコード