C# で LU 分解

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

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

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

A=LU(a11a12a1na21a22a2nan1an2ann)=(l1100l21l220ln1ln2lnn)(u11u12u1n0u22u2n00unn)(a11a12a1na21a22a2nan1an2ann)=(l11u11l11u12l11u1nl21u11l21u12+l22u22l21u1n+l22u2nln1u11ln2u12+ln2u22i=1nlniuin)\begin{aligned} \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{aligned}

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

u1j=a1jl11li1=ai1u11uij=aijk=1ilikukjliilij=aijk=1i2ljkukiuii\begin{aligned} 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{aligned}

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

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

Ax=bLUx=bUx=yLy=b\begin{aligned} \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{aligned}

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

Ly=b(l1100l21l220ln1ln2lnn)(y1y2yn)=(b1b2bn)(y1y2yn)=(b1l11b2l21y1l22bni=1n1lniyilnn)(y1y2yn)=(b1b2l21y1bni=1n1lniyi)(lii=1)\begin{aligned} \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{aligned}

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

Ux=y(u11u12u1n0u22u2n00unn)(x1x2xn)=(y1y2yn)(x1x2xn)=(y1i=2nu1iyiu11y2i=3nu2iyiu22ynunn)\begin{aligned} \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{aligned}

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

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

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

69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
// 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 回行うだけです。下に、 Ly=b\boldsymbol{L}\boldsymbol{y}=\boldsymbol{b}Ux=y\boldsymbol{U}\boldsymbol{x}=\boldsymbol{y} の処理を順に示します。

110
111
112
113
114
115
116
117
118
119
120
// 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]
}
127
128
129
130
131
132
133
134
135
136
137
// 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];
}

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

Solve()
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
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();
0
/* LU分解 */
100
101
102
103
104
105
106
107
108
Console.WriteLine("LU分解後");
Console.WriteLine(nameof(matrixA));
WriteMatrix(matrixA);
// おまけ
Console.WriteLine(nameof(matrixL));
WriteMatrix(matrixL);
Console.WriteLine(nameof(matrixU));
WriteMatrix(matrixU);
Console.WriteLine();
0
/* Ly = b */
122
123
124
125
Console.WriteLine("Ly = b の後");
Console.WriteLine(nameof(vectorY));
WriteVector(vectorY);
Console.WriteLine();
0
/* Ux = y */
139
140
141
Console.WriteLine("Ux = y の後");
Console.WriteLine("解");
WriteVector(solutionX);

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

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
初期状態
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 の情報を保持できること、行列×ベクトルの計算が多いことがわかります。