1. Home
  2. プログラマーのための線形代数
  3. 行列の分解
  4. LU分解

LU分解

LU分解は行列の分解の中で最も基本的なものです。当ページでは、このLU分解について知っておきたいことをまとめて確認することができます。具体的には以下のことがわかります。

当ページでわかること

  • LU分解とは?
  • LU分解のやり方
  • LU分解可能な行列の条件
  • ほとんどどの行列でも可能なPLU分解とは?
  • PythonでLP分解

早速見ていきましょう。

LU分解とは?

LU 分解は、行列の分解方法の一つで、正方行列 A を下三角行列 L と上三角行列 U の積に分解するというものです。

\[
A=L \cdot U
\]

たとえば以下の行列 A は、ご覧の通り、2 つの三角行列 L と U の積で表すことが可能です。

\[\begin{eqnarray}
\overset{A}{
\begin{pmatrix}
5 & 6 & 7 \\
10 & 20 & 23 \\
15 & 50 & 67
\end{pmatrix}
}
=
\overset{L}{
\begin{pmatrix}
1 & 0 & 0 \\
2 & 1 & 0 \\
3 & 4 & 1
\end{pmatrix}}
\overset{U}{
\begin{pmatrix}
5 & 6 & 7 \\
0 & 8 & 9 \\
0 & 0 & 10
\end{pmatrix}}
\end{eqnarray}\]

このようにある行列を下三角行列と上三角行列に分解するのが LU 分解です。これには次のような利点があります。

LU分解の利点

  • 連立一次方程式を解く際の計算量がガウスの消去法よりも少なくて済む。
  • 逆行列を簡単に求められるようになる。
  • 行列式の値を簡単に求められるようになる。
  • 統計学の線形回帰式の係数を求められるようになる。

要するに LU 分解は、行列のその他の重要な演算や操作をするときに役立つものであるということです。プログラミングでも目にする機会がたまにあります。

LU分解のやり方

続いて LU 分解のやり方を見ていきましょう。ポイントは「一方の三角行列の対角線の要素は 1 とする」という点です。このように決めておくと、行列 L と行列 U の積は次のように表すことができます。

\[\begin{eqnarray}
\overset{A}{
\begin{pmatrix}
5 & 6 & 7 \\
10 & 20 & 23 \\
15 & 50 & 67
\end{pmatrix}
}
&=&
\overset{L}{
\begin{pmatrix}
1 & 0 & 0 \\
l_{21} & 1 & 0 \\
l_{31} & l_{32} & 1
\end{pmatrix}
}
\overset{U}{
\begin{pmatrix}
u_{11} & u_{12} & u_{13} \\
0 & u_{22} & u_{23} \\
0 & 0 & u_{33}
\end{pmatrix}}\\
&=&
\begin{pmatrix}
u_{11} &u_{12} & u_{13} \\
l_{21}u_{11} & l_{21}u_{12}+u_{22} & l_{21}u_{13}+u_{23} \\
l_{31}u_{11} & l_{31}u_{12}+l_{32}u_{22} & l_{31}u_{13}+l_{32}u_{23}+u_{33}
\end{pmatrix}
\end{eqnarray}\]

このことから、数値と記号の対応関係は以下の通りになることがわかります。

  • \(u_{11}=5\)
  • \(u_{12}=6\)
  • \(u_{13}=7\)
  • \(l_{21}u_{11}=10\)
  • \(l_{21}u_{12}+u_{22}=20\)
  • \(l_{21}u_{13}+u_{23}=23\)
  • \(l_{31}u_{11}=15\)
  • \(l_{31}u_{12}+l_{32}u_{22}=50\)
  • \(l_{31}u_{13}+l_{32}u_{23}+u_{33}=67\)

あとはこれを上の行から順番に解いていくだけです。

1 行目は、「\(u_{11}=5\)」、「\(u_{12}=6\)」、「\(u_{13}=7\)」 で既に答えが出ています。

続いて、2 行目はまず「\(l_{21}u_{11}=10\)」から「\(l_{21} \cdot 5 =10\)」となります。そこから「\(l_{21}=2\)」であることが求められます。そして「\(l_{21}u_{12}+u_{22}=20\)」から「\(2 \cdot 6+u_{22}=20\)」となります。そこから「\(u_{22}=8\)」であることが求められます。同じようにして「\(l_{21}u_{13}+u_{23}=23\)」から「\(u_{23}=9\)」であることが求められます。

3 行目も同じです。まず「\(l_{31}u_{11}=15\)」から「\(l_{31}=3\)」であることがわかります。続いて\(l_{31}u_{12}+l_{32}u_{22}=50\)」 から「\(l_{32}=4\)」であることがわかります。最後に\(l_{31}u_{13}+l_{32}u_{23}+u_{33}=67\)」から「\(u_{33}=10\)」であることがわかります。

こうやって上から順番に解いていくことで、全ての記号の値が求められます。以上のことから、行列 \(A\) を分解したときの下方行列 \(L\) と上方行列 \(U\) が求められます。

\[\begin{eqnarray}
\begin{pmatrix}
5 & 6 & 7 \\
10 & 20 & 23 \\
15 & 50 & 67
\end{pmatrix}
=
\begin{pmatrix}
1 & 0 & 0 \\
2 & 1 & 0 \\
3 & 4 & 1
\end{pmatrix}
\begin{pmatrix}
5 & 6 & 7 \\
0 & 8 & 9 \\
0 & 0 & 10
\end{pmatrix}
\end{eqnarray}\]

以上がLU分解のやり方です。

LU分解可能な行列の条件

なおすべての正方行列で LU 分解が可能わけではありません。

ある行列が LU 分解可能なのは、すべての「主座小行列」の行列式の値が 0 ではない場合です。主座小行列とは、行列の左上に位置するすべての部分行列のことです。

たとえば、以下の正方行列があるとします。

\[\begin{eqnarray}
A
=
\begin{pmatrix}
5 & 6 & 7 \\
10 & 20 & 23 \\
15 & 50 & 67
\end{pmatrix}
\end{eqnarray}\]

この場合、次の 3 つが主座小行列になります。

\[\begin{eqnarray}
\begin{pmatrix}
5
\end{pmatrix}, \ \ \
\begin{pmatrix}
5 & 6 \\
10 & 20
\end{pmatrix}, \ \ \
\begin{pmatrix}
5 & 6 & 7 \\
10 & 20 & 23 \\
15 & 50 & 67
\end{pmatrix}
\end{eqnarray}\]

これらのすべての主座小行列の行列式の値が 0 でないときのみ LU 分解が可能です。

ほぼどんな行列でも可能なPLU分解とは?

上述の通りLU 分解では分解ができない場合があります。そこで、ほとんどどのような行列でもLU分解することを可能にしたPLU分解というものがあります。

P は置換行列(ある行列の行や列を入れ替える行列)です。

\[A=PLU\]

たとえば、以下の行列を普通にLU分解してみましょう。

\[\begin{eqnarray}
\overset{A}{
\begin{pmatrix}
0 & 6 & 7 \\
10 & 20 & 23 \\
15 & 50 & 67
\end{pmatrix}
}
&=&
\overset{L}{
\begin{pmatrix}
1 & 0 & 0 \\
l_{21} & 1 & 0 \\
l_{31} & l_{32} & 1
\end{pmatrix}
}
\overset{U}{
\begin{pmatrix}
u_{11} & u_{12} & u_{13} \\
0 & u_{22} & u_{23} \\
0 & 0 & u_{33}
\end{pmatrix}}\\
&=&
\begin{pmatrix}
u_{11} &u_{12} & u_{13} \\
l_{21}u_{11} & l_{21}u_{12}+u_{22} & l_{21}u_{13}+u_{23} \\
l_{31}u_{11} & l_{31}u_{12}+l_{32}u_{22} & l_{31}u_{13}+l_{32}u_{23}+u_{33}
\end{pmatrix}
\end{eqnarray}\]

この行列をLU分解しようとしても、主座小行列の行列式の値が 0 なのでLU分解できないことに気づきます。実際に確かめてみましょう。まず、この行列の積の数値と記号の対応関係は以下の通りになります。

  • \(u_{11}=0\)
  • \(u_{12}=6\)
  • \(u_{13}=7\)
  • \(l_{21}u_{11}=10\)
  • \(l_{21}u_{12}+u_{22}=20\)
  • \(l_{21}u_{13}+u_{23}=23\)
  • \(l_{31}u_{11}=15\)
  • \(l_{31}u_{12}+l_{32}u_{22}=50\)
  • \(l_{31}u_{13}+l_{32}u_{23}+u_{33}=67\)

これを順番に解いていくだけなのですが、今回は \(u_{11}\) が 0 なので、いきなり \(l_{21}u_{11}=10\) が解けないことがわかります。

このようなときに行を入れ替えることで、LU分解を可能にしてくれるのが置換行列です。たとえば、行列 \(A\) に以下の置換行列 \(P\) を左から掛けると、1 行目と 3 行目が入れ替えられます。

\[\begin{eqnarray}
\overset{P}{
\begin{pmatrix}
0 & 0 & 1 \\
0 & 1 & 0 \\
1 & 0 & 0
\end{pmatrix}
}
\overset{A}{
\begin{pmatrix}
0 & 6 & 7 \\
10 & 20 & 23 \\
15 & 50 & 67
\end{pmatrix}
}
=
\overset{PA}{
\begin{pmatrix}
15 & 50 & 67 \\
10 & 20 & 23 \\
0 & 6 & 7
\end{pmatrix}
}
\end{eqnarray}\]

このように置換しておくと、主座小行列の行列式の値が 0 ではなくなるので、LU分解が可能になります。実際に解くと以下のようになります。

\[\begin{eqnarray}
\overset{PA}{
\begin{pmatrix}
15 & 50 & 67 \\
10 & 20 & 23 \\
0 & 6 & 7
\end{pmatrix}
}
=
\overset{L}{
\begin{pmatrix}
1 & 0 & 0 \\
0.667 & 1 & 0 \\
0.333 & 0.8 & 1
\end{pmatrix}
}
\overset{U}{
\begin{pmatrix}
15 & 50 & 60 \\
0 & -13.33 & -21.67 \\
0 & 0 & 2
\end{pmatrix}
}
\end{eqnarray}\]

以上のことから、元々の行列 \(A\) は次のようにPLU分解できるということになります。

\[\begin{eqnarray}
\overset{A}{
\begin{pmatrix}
0 & 6 & 7 \\
10 & 20 & 23 \\
15 & 50 & 67
\end{pmatrix}
}
=
\overset{P}{
\begin{pmatrix}
0 & 0 & 1 \\
0 & 1 & 0 \\
1 & 0 & 0
\end{pmatrix}
}
\overset{L}{
\begin{pmatrix}
1 & 0 & 0 \\
0.667 & 1 & 0 \\
0.333 & 0.8 & 1
\end{pmatrix}
}
\overset{U}{
\begin{pmatrix}
15 & 50 & 60 \\
0 & -13.33 & -21.67 \\
0 & 0 & 2
\end{pmatrix}
}
\end{eqnarray}\]

以上がPLU分解です。

こちらの方法だと、ほとんどの行列を安定してLU分解することが可能であり、システムの安定性が増すため、プログラミング言語の Python では、PLU分解がデフォルトになっています。

PythonでLU分解

LU分解は、Python では SciyPy ライブラリの lu() 関数で実装されています。より具体的には、この関数は PLU分解を実行します。以下の通りです。

In [1]:
# NumPy と SciPy のインポート
import numpy as np
from scipy.linalg import lu

# 正方行列を定義
A = np.array([
    [1,2,3],
    [4,5,6],
    [7,8,9]])
print(A)
[[1 2 3]
 [4 5 6]
 [7 8 9]]
In [2]:
# PLU分解
P,L,U = lu(A)
print(P)
[[0. 1. 0.]
 [0. 0. 1.]
 [1. 0. 0.]]
In [3]:
print(L)
[[1.         0.         0.        ]
 [0.14285714 1.         0.        ]
 [0.57142857 0.5        1.        ]]
In [4]:
print(U)
[[ 7.00000000e+00  8.00000000e+00  9.00000000e+00]
 [ 0.00000000e+00  8.57142857e-01  1.71428571e+00]
 [ 0.00000000e+00  0.00000000e+00 -1.58603289e-16]]
In [5]:
# 元々の行列の再作成
B= P @ L @ U 
print(B)
[[1. 2. 3.]
 [4. 5. 6.]
 [7. 8. 9.]]