はじめに

LU分解に関する記事はいくつもありますが、省メモリで計算を実行できるというLU分解の特徴に言及している記事がなかったので、改めて省メモリに注目してかつ少し丁寧にLU分解を解説しようと思います。

Reactの勉強がてら簡単LU分解を実行できるサイトを作成してみました。ただ、レンダリングの都合上、以下で述べる省メモリな実装は行っていません。LU分解体験

三角行列

LU分解について説明する前に、三角行列について説明します。LU分解のゴールは、連立方程式の係数を三角行列のみで表すことです。その理由を説明します。まず、簡単な連立方程式を考えてみます。

\left\{
\begin{array}{rrrrrrr}
2x_1 & &      & &      &=&  2 \\
 x_1 &+& 4x_2 & &      &=&  9 \\
4x_1 &-& 3x_2 &+& 3x_3 &=& -5
\end{array}
\right.

この連立方程式の解は、逐次的に解いていくことで簡単に求めることができます。

\begin{align}
x_1 &= \frac{2}{2} = 1 \\
x_2 &= \frac{9 - x_1}{4} = \frac{9-1}{4} = 2 \\
x_3 &= \frac{-5 - 4x_1 + 3x_2}{3} = \frac{-5 - 4 + 6}{3} = -1
\end{align}

ここで、この連立方程式を行列で表すと、

\left(
\begin{matrix}
2 & 0 & 0 \\
1 & 4 & 0 \\
4 & -3 & 3
\end{matrix}
\right)
\left(
\begin{matrix}
x_1 \\
x_2 \\
x_3
\end{matrix}
\right)
=
\left(
\begin{matrix}
2 \\
9 \\
-5
\end{matrix}
\right) \tag{1}

のようになります。連立方程式の係数が、対角成分より上の要素が$0$である下三角行列で表されていることが分かります。この場合、$x_i (i \in 2, 3)$を求める時にはすでに$x_{i-1} (i \in 1,2)$が分かっているので、順に解くことができます。この方法を前進代入と呼びます。次に、下のような連立方程式を考えてみます。

\left\{
\begin{array}{rrrrrrr}
3x_1 &-& 3x_2 &+& 4x_3 &=& -5 \\
     & & 4x_2 &+&  x_3 &=&  9 \\
     & &      & & 2x_3 &=&  2
\end{array}
\right. \\

\left(
\begin{matrix}
3 & -3 & 4\\
0 & 4 & 1 \\
0 & 0 & 2
\end{matrix}
\right)
\left(
\begin{matrix}
x_1 \\
x_2 \\
x_3
\end{matrix}
\right)
=
\left(
\begin{matrix}
-5 \\
9 \\
2
\end{matrix}
\right)

この場合、連立方程式の係数は、対角成分より下の成分が$0$である上三角行列で表されます。この場合は、$x_i (i \in 2, 1)$を求める際には、$x_{i + 1} (i \in 3, 2)$が分かっているので、下三角行列の場合の逆順に解くことができます。この方法を後退代入と呼びます。このように、連立方程式の係数が三角行列で表される場合、簡単なアルゴリズムで解くことができます。ここで、言葉の定義をしておきます。式$(1)$を次ように表現します。

A \boldsymbol{x} = \boldsymbol{b} \\
A = 
\left(
\begin{matrix}
3 & -3 & 4\\
0 & 4 & 1 \\
0 & 0 & 2
\end{matrix}
\right), \ 
\boldsymbol{x} = 
\left(
\begin{matrix}
x_1 \\
x_2 \\
x_3
\end{matrix}
\right), \ 
\boldsymbol{b} =
\left(
\begin{matrix}
2 \\
9 \\
-5
\end{matrix}
\right)

ここで、$A$を系数行列、$\boldsymbol{x}$を解ベクトル、$\boldsymbol{b}$を右辺ベクトルと呼びます。以下では、この表記を一般化し、$A$を$n \times n$行列、$\boldsymbol{x}, \ \boldsymbol{b}$を$n$次元ベクトルとします。

LU分解

上で述べた通り、連立方程式の係数を三角行列で表すことができれば、その方程式ほぼ解けたようなものです。ここで、係数行列$A$を下三角行列$L$と上三角行列$U$の積で表現することを考えます。

LU \boldsymbol{x} = \boldsymbol{b} \\

$U \boldsymbol{x}=\boldsymbol{y}$とすると、$\boldsymbol{y}$については前進代入で求まり、その結果を用いることで、$\boldsymbol{x}$については後退代入によって求められることが分かります。

\begin{align}
L \boldsymbol{y} &= \boldsymbol{b} \\
U \boldsymbol{x} &= \boldsymbol{y}
\end{align}

では、どのようにして$L$と$U$を求めるのでしょうか。本題に入る前に、自由度に注目してみます。$A$の要素は$n^2$個、$L, \ U$の要素はそれぞれ$n(n+1)/2$個あり、合計で$n^2 + n$個の要素があります。したがって、分解後は自由度が$n$大きいですが、これは、$L$の対角成分が全て$1$であるとすることで、自由度が$n^2$となり、$A$に対して$L, \ U$が一意に決めることができます。（ちなみに、$L$の対角成分を$1$とするか、$U$の対角成分を$1$とするかはどちらでも構いません。人によります。）

求め方

では、具体的な求め方を説明しますが、その前に、小行列について説明します。

小行列

次のような行列$M$を考えます。

M = \left(
\begin{matrix}
m_{11} & m_{12} &|& m_{13} & m_{14} & m_{15} & m_{16} \\
m_{21} & m_{22} &|& m_{23} & m_{24} & m_{25} & m_{26} \\
- & - & - & - & - & - & - \\
m_{31} & m_{32} &|& m_{33} & m_{34} & m_{35} & m_{36} \\
m_{41} & m_{42} &|& m_{43} & m_{44} & m_{45} & m_{46} \\
m_{51} & m_{52} &|& m_{53} & m_{54} & m_{55} & m_{56} \\
m_{61} & m_{62} &|& m_{63} & m_{64} & m_{65} & m_{66}
\end{matrix}
\right)

この行列を点線で区切り、それぞれのブロック（小行列）に名前を付けます。

M_{11} = 
\left(
\begin{matrix}
m_{11} & m_{12} \\
m_{21} & m_{22} 
\end{matrix}
\right)
, \ 
M_{12} = 
\left(
\begin{matrix}
m_{13} & m_{14} & m_{15} & m_{16} \\
m_{23} & m_{24} & m_{25} & m_{26}
\end{matrix}
\right) \\

M_{21} =
\left(
\begin{matrix}
m_{31} & m_{32} \\
m_{41} & m_{42} \\
m_{51} & m_{52} \\
m_{61} & m_{62}
\end{matrix}
\right)
, \ 
M_{22} = 
\left(
\begin{matrix}
m_{33} & m_{34} & m_{35} & m_{36} \\
m_{43} & m_{44} & m_{45} & m_{46} \\
m_{53} & m_{54} & m_{55} & m_{56} \\
m_{63} & m_{64} & m_{65} & m_{66}
\end{matrix}
\right)

すると行列$M$は次のように表されます。

M = 
\left(
\begin{matrix}
M_{11} & M_{12} \\
M_{21} & M_{22}
\end{matrix}
\right)

ここで、$6 \times 6$の行列$N$を考え、 $M$と同様にして区切って表してみます。

N = 
\left(
\begin{matrix}
N_{11} & N_{12} \\
N_{21} & N_{22}
\end{matrix}
\right)

そして、$M$と$N$の積は実は次のように表すことができます。

\begin{align}
MN &= 
\left(
\begin{matrix}
M_{11} & M_{12} \\
M_{21} & M_{22}
\end{matrix}
\right)
\left(
\begin{matrix}
N_{11} & N_{12} \\
N_{21} & N_{22}
\end{matrix}
\right) \\
&=
\left(
\begin{matrix}
N_{11}N_{11} + N_{12}N_{21}  & N_{11}N_{12} + N_{12}N_{22} \\
N_{21}N_{11} + N_{22}N_{21} & N_{21}N_{12} + N_{22}N_{22}
\end{matrix}
\right) \\
\end{align}

つまり、小行列を行列の要素のように扱うことができるということです。（実際に計算してみると分かります。）ただし、それぞれの小行列は矛盾なく計算できるように$M$と$N$を分割する必要があります。

求め方本題

前置きが長くなりましたが、ここでやっと行列$L$と$U$を求める手順を示します。まず、$L$と$U$を下のように小行列に分割します。ここで$O$は全ての要素が$0$の行列です。

L = 
\left(
\begin{matrix}
1 & O \\
L_{21} & L_{22}
\end{matrix}
\right), \ 

U = 
\left(
\begin{matrix}
u_{11} & U_{12} \\
O & U_{22}
\end{matrix}
\right)

$L_{22},　U_{22}$はそれぞれ、$(n-1) \times (n-1)$の下三角行列と上三角行列です。また、係数行列$A$についても同様に分割します。

A = 
\left(
\begin{matrix}
a_{11} & A_{12} \\
A_{21} & A_{22}
\end{matrix}
\right)

$A=LU$ですので、

\begin{align}
\left(
\begin{matrix}
a_{11} & A_{12} \\
A_{21} & A_{22}
\end{matrix}
\right)
&=
\left(
\begin{matrix}
1 & O \\
L_{21} & L_{22}
\end{matrix}
\right)
\left(
\begin{matrix}
u_{11} & U_{12} \\
O & U_{22}
\end{matrix}
\right) \\
&=
\left(
\begin{matrix}
u_{11} & U_{12} \\
L_{21}u_{11} & L_{21}U_{12} + L_{22}U_{22}
\end{matrix}
\right)
\end{align}

と表されます。この関係から、

u_{11} = a_{11}, \ U_{12} = A_{12}

を得ます。ここで、$a_{11} \neq 0$と仮定すると、

L_{21} = \frac{A_{21}}{a_{11}}

と$L_{21}$も決まります。これにより、$L$の第１列と$U$の第1行が求まります。さて、ここまではおそらく容易に理解できると思います。問題は、$A_{22} = L_{21}U_{12} + L_{22}U_{22}$をどう処理するかです。ここで、先ほど述べたのように、$L_{21}$と$U_{12}$はすでに求まっているので、右辺第1項はすでに求まっています。次に、右辺第2項ですが、求め方本題の冒頭で述べた通り、$L_{22}$は下三角行列、$U_{22}$は上三角行列です。したがって、第1項を移項して、

\hat{A}_{22} = A_{22} - L_{21}U_{12}

とすれば、

A^\prime_{22} = L_{22}U_{22}

と表現することができます。これは、すなわち、$(n-1) \times (n-1)$の行列$A^\prime_{22}$を下三角行列と上三角行列に分解していることにほかなりません。したがって、この$A^\prime_{22}$に対して、上と同様の手順を施すことで、$(n-2) \times (n-2)$の行列のLU分解の問題に帰着し、さらに、$(n-3) \times (n-3)$の行列の行列のLU分解の問題に帰着し...というように逐次的に解くことができます。ここまで理解できれば、あとはコードに落とし込むだけです。

実装

省メモリな実装

では実際に行列$L$と$U$を求めるコードを書いていくのですが、何も考えなければ、係数行列$A$用の二次元配列、下三角行列$L$用の二次元配列、上三角行列$U$用の二次元配列の三つを変数として持とうと思うのではないでしょうか。実際、LU分解について紹介している記事（例えば、LU分解法をプログラムで書いてみた）では、それぞれの行列用の変数を定義しています。しかし、破壊的な方法ではありますが、$A$を上書きしていくという方法をとれば、必要な変数は1つで十分です。具体的には、$L$の対角成分が$1$であることを利用し、次のように$L$と$U$を1つの行列$V$で管理することを考えます。

L =
\left(
\begin{matrix}
1 & 0 & 0 \\
l_{11} & 1 & 0  \\
l_{21} & l_{22} & 1  \\
\end{matrix}
\right), \

U =
\left(
\begin{matrix}
u_{11} & u_{12} & u_{13} \\
0 & u_{22} & u_{22}  \\
0 & 0 & u_{33}  \\
\end{matrix}
\right) \\

V =
\left(
\begin{matrix}
u_{11} & u_{12} & u_{13} \\
l_{11}  & u_{22} & u_{22}  \\
l_{21} & l_{22} & u_{33}  \\
\end{matrix}
\right)

このような$V$を考え、さらに$A$に$V$を上書きしていくことで、配列の保存に必要なメモリは$A$の分のみとなり、メモリを節約することができます。

実装本題

前置きがかなりのボリュームになってしまいましたが、Pythonでの実装例を紹介しておきます。長々と説明しましたが、実装自体は驚くほど簡単です。

import numpy as np

def main():
    n = int(input())
    A = np.array([np.array(list(map(float, input().split()))) for i in range(n)])

    for i in range(n):
        for j in range(i + 1, n):
            A[j][i] /= A[i][i]
            for k in range(i + 1, n):
                A[j][k] -= A[j][i] * A[i][k]
    return A

if __name__ == "__main__":
    print(main())

部分ピポット選択

以上で基本的なLU分解の解説は終了ですが、$a_{ii} \neq 0$を仮定しているため、係数行列の対角成分に$0$が含まれる場合、ZeroDivisionErrorで計算がストップしてしまいます。また、$0$ではなくても、小さい値で割り算をするのは、精神衛生上もよろしくないですしね。というわけで、なんとかしないといけないのですが、ピボット選択という手法でなんとかします。この手法では、$a_{ii}$に注目する際に、$a_{ii}$から$a_{ni}$の絶対値を比較し、一番大きい値を含む行と$i$行目を入れ替えてから、計算を実行します。こうすることで、0徐算や小さい値での徐算を防ぐことができます。行だけだはなく、列も含んで最大値を探す方法を完全ピボット選択と呼ぶのに対し、上で説明したように行のみに注目する方法を部分ピボット選択と呼びます。部分ピボット選択で実用上は十分です。

概要

部分ピボット選択を実行するに当たり、以下のような行列$P_{k,j}$を考えます。

P_{k,j} =
 
\left(
\begin{matrix}
1 &   &   &   &  &         \\
  & 1 &   &   &  &         \\
  &   & 1 &   &  &         \\
  &   &   & 0 &  &  1      \\
  &   &   &   & \ddots     \\
  &   &   & 1 &  &  0      \\
  &   &   &   &  &  & 1    \\
  &   &   &   &  &  &  & 1 \\
\end{matrix}
\right)

この行列は、第$k$列と第$j$列を入れ替えた単位行列です。この行列を左側からかけると第$k$行と第$j$行を入れ替えることができます。このような行列を置換行列と呼びます。この行列を毎回のLU分解計算時に左側から順次係数行列$A$にかけていくことで、ピボット選択を表現することができます。第$k$段階目の置換行列を$P^{(k)}$と表すと、ピボット選択をしながらLU分解を行うことは、

P^{(n-1)} \dots P^{(2)} P^{(1)} A = LU

と表すことができます。ここで、

P = P^{(n-1)} \dots P^{(2)} P^{(1)}

とおけば、元々解きたい方程式$A \boldsymbol{x} = \boldsymbol{b}$は

PA \boldsymbol{x} = P \boldsymbol{b} \\
\therefore LU \boldsymbol{x} = P \boldsymbol{b}

と表すことができます。したがって、LU分解を実行する際に、置換行列の積を同時に計算し、以下のように上三角方程式・下三角方程式を得ることができます。

L \boldsymbol{y} = P \boldsymbol{b} \\
U \boldsymbol{x} = \boldsymbol{y}

実装

以下、ピボット選択を導入したときのLU分解のコードです。

import numpy as np

def main():
    n = int(input())
    A = np.array([np.array(list(map(float, input().split()))) for i in range(n)])
    P = np.identity(n)
    for i in range(n):
        maxidx = np.argmax(abs(A[i:n, i])) + i
        A[i, :], A[maxidx, :] = A[maxidx, :], A[i, :].copy()
        P[i, :], P[maxidx, :] = P[maxidx, :], P[i, :].copy()
        for j in range(i + 1, n):
            A[j][i] /= A[i][i]
            for k in range(i + 1, n):
                A[j][k] -= A[j][i] * A[i][k]
    return A, P

if __name__ == "__main__":
    A, P = main()

参考

工学のための数値計算

最後に

何か質問があればなんでも聞いてください！あまり数学に詳しいわけではないので、するどいまさかりには対応できないかもしれません。