C# で LU 分解

Tech
Share

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(a11a12⋯a1na21a22⋯a2n⋮⋮⋱⋮an1an2⋯ann)=(l110⋯0l21l22⋯0⋮⋮⋱⋮ln1ln2⋯lnn)(u11u12⋯u1n0u22⋯u2n⋮⋮⋱⋮00⋯unn)(a11a12⋯a1na21a22⋯a2n⋮⋮⋱⋮an1an2⋯ann)=(l11u11l11u12⋯l11u1nl21u11l21u12+l22u22⋯l21u1n+l22u2n⋮⋮⋱⋮ln1u11ln2u12+ln2u22⋯∑i=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=aij−∑k=1ilikukjliilij=aij−∑k=1i−2ljkukiuii\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(l110⋯0l21l22⋯0⋮⋮⋱⋮ln1ln2⋯lnn)(y1y2⋮yn)=(b1b2⋮bn)(y1y2⋮yn)=(b1l11b2−l21y1l22⋮bn−∑i=1n−1lniyilnn)(y1y2⋮yn)=(b1b2−l21y1⋮bn−∑i=1n−1lniyi)  (∵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(u11u12⋯u1n0u22⋯u2n⋮⋮⋱⋮00⋯unn)(x1x2⋮xn)=(y1y2⋮yn)(x1x2⋮xn)=(y1−∑i=2nu1iyiu11y2−∑i=3nu2iyiu22⋮ynunn)\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];
    }
}

プログラム中にある通り、 matrixL や matrixU をわざわざ用意する必要はなく 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 だけで matrixL と matrixU の情報を保持できること、行列×ベクトルの計算が多いことがわかります。