C#でLU分解

LU分解とは、ある行列AをL(下三角行列)とU(上三角行列)の積に分解するものです。

LU分解が役に立つのは、同じ行列に対してベクトルが違う方程式をいくつも解く場合です。LU分解自体はガウスの消去法と同じオーダーの計算量となりますが、LU分解を用いる場合、分解は方程式がいくつあっても1度しか行わないため、その分方程式を解くのにかかる時間は少なくなります。

解くべき方程式を%%$%\boldsymbol{A}\boldsymbol{x}=\boldsymbol{b}%$%%とするとき、%%$%\boldsymbol{A}%$%%をLU分解するとします。このとき、

##$#\begin{align*}
\boldsymbol{A} & =\boldsymbol{L}\boldsymbol{U}\\
\left(\begin{array}{cccc}
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) & =\left(\begin{array}{cccc}
l_{11} & 0 & \cdots & 0\\
l_{21} & l_{22} & \cdots & 0\\
\vdots & \vdots & \ddots & \vdots\\
l_{n1} & l_{n2} & \cdots & l_{nn}
\end{array}\right)\left(\begin{array}{cccc}
u_{11} & u_{12} & \cdots & u_{1n}\\
0 & u_{22} & \cdots & u_{2n}\\
\vdots & \vdots & \ddots & \vdots\\
0 & 0 & \cdots & u_{nn}
\end{array}\right)\\
\left(\begin{array}{cccc}
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) & =\left(\begin{array}{cccc}
l_{11}u_{11} & l_{11}u_{12} & \cdots & l_{11}u_{1n}\\
l_{21}u_{11} & l_{21}u_{12}+l_{22}u_{22} & \cdots & l_{21}u_{1n}+l_{22}u_{2n}\\
\vdots & \vdots & \ddots & \vdots\\
l_{n1}u_{11} & l_{n2}u_{12}+l_{n2}u_{22} & \cdots & \sum_{i=1}^{n}l_{ni}u_{in}
\end{array}\right)
\end{align*}#$##

であることより、%%$%\boldsymbol{L}%$%%と%%$%\boldsymbol{U}%$%%の要素は下のように計算されます。

##$#\begin{align*}
u_{1j} & =\frac{a_{1j}}{l_{11}}\\
l_{i1} & =\frac{a_{i1}}{u_{11}}\\
u_{ij} & =\frac{a_{ij}-\sum_{k=1}^{i}l_{ik}u_{kj}}{l_{ii}}\\
l_{ij} & =\frac{a_{ij}-\sum_{k=1}^{i-2}l_{jk}u_{ki}}{u_{ii}}
\end{align*}#$##

プログラムで実装する際、%%$%\boldsymbol{L}%$%%の対角要素を1と決めておくことにより、%%$%\boldsymbol{L}%$%%と%%$%\boldsymbol{U}%$%%の要素の算出を同時に行い、2つまとめて%%$%\boldsymbol{A}%$%%として保持するという処理が行われます。処理・領域の効率が良いので、プログラムでは大抵この技法を取っています。

これらの行列から%%$%\boldsymbol{x}%$%%を求めるためには、計算の間に新たなベクトル%%$%\boldsymbol{y}%$%%を置くようにすることで、

##$#\begin{align*}
\boldsymbol{A}\boldsymbol{x} & =\boldsymbol{b}\\
\boldsymbol{L}\boldsymbol{U}\boldsymbol{x} & =\boldsymbol{b}\\
\boldsymbol{U}\boldsymbol{x} & =\boldsymbol{y}\\
\boldsymbol{L}\boldsymbol{y} & =\boldsymbol{b}
\end{align*}#$##

という式を使って解いていきます。この処理は以前の処理と似ており、

##$#\begin{align*}
\boldsymbol{L}\boldsymbol{y} & =\boldsymbol{b}\\
\left(\begin{array}{cccc}
l_{11} & 0 & \cdots & 0\\
l_{21} & l_{22} & \cdots & 0\\
\vdots & \vdots & \ddots & \vdots\\
l_{n1} & l_{n2} & \cdots & l_{nn}
\end{array}\right)\left(\begin{array}{c}
y_{1}\\
y_{2}\\
\vdots\\
y_{n}
\end{array}\right) & =\left(\begin{array}{c}
b_{1}\\
b_{2}\\
\vdots\\
b_{n}
\end{array}\right)\\
\left(\begin{array}{c}
y_{1}\\
y_{2}\\
\vdots\\
y_{n}
\end{array}\right) & =\left(\begin{array}{c}
\frac{b_{1}}{l_{11}}\\
\frac{b_{2}-l_{21}y_{1}}{l_{22}}\\
\vdots\\
\frac{b_{n}-\sum_{i=1}^{n-1}l_{ni}y_{i}}{l_{nn}}
\end{array}\right)\\
\left(\begin{array}{c}
y_{1}\\
y_{2}\\
\vdots\\
y_{n}
\end{array}\right) & =\left(\begin{array}{c}
b_{1}\\
b_{2}-l_{21}y_{1}\\
\vdots\\
b_{n}-\sum_{i=1}^{n-1}l_{ni}y_{i}
\end{array}\right)\;\left(\because l_{ii}=1\right)
\end{align*}#$##

によってまず%%$%\boldsymbol{y}%$%%を求め、それから

##$#\begin{align*}
\boldsymbol{U}\boldsymbol{x} & =\boldsymbol{y}\\
\left(\begin{array}{cccc}
u_{11} & u_{12} & \cdots & u_{1n}\\
0 & u_{22} & \cdots & u_{2n}\\
\vdots & \vdots & \ddots & \vdots\\
0 & 0 & \cdots & u_{nn}
\end{array}\right)\left(\begin{array}{c}
x_{1}\\
x_{2}\\
\vdots\\
x_{n}
\end{array}\right) & =\left(\begin{array}{c}
y_{1}\\
y_{2}\\
\vdots\\
y_{n}
\end{array}\right)\\
\left(\begin{array}{c}
x_{1}\\
x_{2}\\
\vdots\\
x_{n}
\end{array}\right) & =\left(\begin{array}{c}
\frac{y_{1}-\sum_{i=2}^{n}u_{1i}y_{i}}{u_{11}}\\
\frac{y_{2}-\sum_{i=3}^{n}u_{2i}y_{i}}{u_{22}}\\
\vdots\\
\frac{y_{n}}{u_{nn}}
\end{array}\right)
\end{align*}#$##

によって%%$%\boldsymbol{x}%$%%を求めます。

式は長々していますが、プログラム自体は大きく複雑にはなっていません。ソースコード等はこちらから

繰返しになりますが、LU分解のプログラムでは2つの行列は同時に、同じように求められます。なお、簡単な見た目ながらも計算時間はかかります。その部分が下の通りです。

			// LU分解
			for (int i = 0; i < Dimension - 1; i++)
			{
				for (int j = i + 1; j < Dimension; j++)
				{
					var s = matrixA[j, i] / matrixA[i, i];
					matrixA[j, i] = s;
					for (int k = i + 1; k < Dimension; k++)
					{
						matrixA[j, k] -= matrixA[i, k] * s;
					}
				}
			}

			// サンプル用にコピー 不要な処理
			for (int i = 0; i < Dimension; i++)
			{
				matrixL[i, i] = 1.0;
				for (int j = 0; j < i; j++)
				{
					matrixL[i, j] = matrixA[i, j];
				}
			}
			for (int i = 0; i < Dimension; i++)
			{
				for (int j = i; j < Dimension; j++)
				{
					matrixU[i, j] = matrixA[i, j];
				}
			}

プログラム中にある通り、matrixLmatrixUをわざわざ用意する必要はなくmatrixAの要素が変化するだけなので、LU分解をしても新たなスペースが必要になることはありません。なんと素晴らしいアルゴリズムなのでしょう!

あとはガウスの消去法のような後退代入のような手法を2回行うだけです。下に、%%$%\boldsymbol{L}\boldsymbol{y}=\boldsymbol{b}%$%%と%%$%\boldsymbol{U}\boldsymbol{x}=\boldsymbol{y}%$%%の処理を順に示します。

			// Ly = b
			// "matrixL"は"matrixA"と変えて同じ
			for (int i = 0; i < Dimension; i++)
			{
				var s = vectorB[i];
				for (int j = 0; j < i; j++)
				{
					s -= matrixL[i, j] * vectorY[j];
				}
				vectorY[i] = s; // vectorY[i] = s / matrixL[i,i]
			}
			// Ux = y
			// "matrixU"は"matrixA"と変えて同じ
			for (int i = Dimension - 1; i >= 0; i--)
			{
				var s = vectorY[i];
				for (int j = i + 1; j < Dimension; j++)
				{
					s -= matrixU[i, j] * solutionX[j];
				}
				solutionX[i] = s / matrixU[i, i];
			}

全体としては下のプログラムになります。(一部省略)

		static void Solve()
		{
			const int Dimension = 4;

			double[,] matrixA = new double[Dimension, Dimension]
			{
				{ 2, 3, 1, 4 },
				{ 4, 1, -3, -2 },
				{ -1, 2, 2, 1 },
				{ 3, -4, 4, 3 },
			};

			double[] vectorB = new double[Dimension]
			{
				10,
				0,
				4,
				6
			};

			double[,] matrixL = new double[Dimension, Dimension];
			double[,] matrixU = new double[Dimension, Dimension];

			double[] vectorY = new double[Dimension];
			double[] solutionX = new double[Dimension];

			Console.WriteLine("初期状態");
			Console.WriteLine(nameof(matrixA));
			WriteMatrix(matrixA);
			Console.WriteLine(nameof(vectorB));
			WriteVector(vectorB);
			Console.WriteLine();

			/* LU分解 */

			Console.WriteLine("LU分解後");
			Console.WriteLine(nameof(matrixA));
			WriteMatrix(matrixA);
			// おまけ
			Console.WriteLine(nameof(matrixL));
			WriteMatrix(matrixL);
			Console.WriteLine(nameof(matrixU));
			WriteMatrix(matrixU);
			Console.WriteLine();

			/* Ly = b */

			Console.WriteLine("Ly = b の後");
			Console.WriteLine(nameof(vectorY));
			WriteVector(vectorY);
			Console.WriteLine();

			/* Ux = y */

			Console.WriteLine("Ux = y の後");
			Console.WriteLine("解");
			WriteVector(solutionX);
		}

計算結果は次の通りです。

初期状態
matrixA
  2.0000   3.0000   1.0000   4.0000 
  4.0000   1.0000  -3.0000  -2.0000 
 -1.0000   2.0000   2.0000   1.0000 
  3.0000  -4.0000   4.0000   3.0000 
vectorB
 10.0000
  0.0000
  4.0000
  6.0000

LU分解後
matrixA
  2.0000   3.0000   1.0000   4.0000 
  2.0000  -5.0000  -5.0000 -10.0000 
 -0.5000  -0.7000  -1.0000  -4.0000 
  1.5000   1.7000 -11.0000 -30.0000 
matrixL
  1.0000   0.0000   0.0000   0.0000 
  2.0000   1.0000   0.0000   0.0000 
 -0.5000  -0.7000   1.0000   0.0000 
  1.5000   1.7000 -11.0000   1.0000 
matrixU
  2.0000   3.0000   1.0000   4.0000 
  0.0000  -5.0000  -5.0000 -10.0000 
  0.0000   0.0000  -1.0000  -4.0000 
  0.0000   0.0000   0.0000 -30.0000 

Ly = b の後
vectorY
 10.0000
-20.0000
 -5.0000
-30.0000

Ux = y の後
解
  1.0000
  1.0000
  1.0000
  1.0000

matrixAだけでmatrixLmatrixUの情報を保持できること、行列×ベクトルの計算が多いことがわかります。

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です