C# で LU 分解

LU 分解とは、ある行列 $\boldsymbol{A}$ を $\boldsymbol{L}$（下三角行列）と $\boldsymbol{U}$（上三角行列）の積に分解するものです。 LU 分解が役に立つのは、同じ行列に対してベクトルが違う方程式をいくつも解く場合です。 LU 分解自体はガウスの消去法と同じオーダーの計算量となりますが、 LU 分解を用いる場合、分解は方程式がいくつあっても 1 度しか行わないため、その分方程式を解くのにかかる時間は少なくなります 解くべき方程式を $\boldsymbol{A}\boldsymbol{x}=\boldsymbol{b}$ とするとき、 $\boldsymbol{A}$ を LU 分解するとします。このとき、 $$\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}$$

であることより、 $\boldsymbol{L}$ と $\boldsymbol{U}$ の要素は下のように計算されます。 $$\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}$$

プログラムで実装する際、 $\boldsymbol{L}$ の対角要素を $1$ と決めておくことにより、 $\boldsymbol{L}$ と $\boldsymbol{U}$ の要素の算出を同時に行い、2 つまとめて $\boldsymbol{A}$ として保持するという処理が行われます。処理・領域の効率が良いので、プログラムでは大抵この技法を取っています。 これらの行列から $\boldsymbol{x}$ を求めるためには、計算の間に新たなベクトル $\boldsymbol{y}$ を置くようにすることで、 $$\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}$$

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

$$\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}$$

によってまず $\boldsymbol{y}$ を求め、それから $$\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}$$

によって