はじめに

　連立一次方程式を解くための道具（ライブラリ）は、有名どころでは BLAS (Basic Linear Algebra Subprograms) や LAPACK (Linear Algebra PACKage) があるし、Python ならば NumPy や SymPy でサクッと解けるだろう。

　だが、数値計算屋の端くれとしては最低限押さえておきたいアルゴリズムだし、それなのにしばらく離れていると「あれ？　どうだったっけ？」と情けないことになるので、備忘録として残しておこうと思う。

　Python でのサンプルプログラムと、私なりの LU 分解の直感的な考え方を。

　なお、私は LU 分解の手法として、上三角行列 $U$ の対角成分を $1$ にするやり方、いわゆる Crout 型に慣れているので以下のアルゴリズムはそれに従っている。巷では下三角行列 $L$ の対角成分を $1$ にする方法、いわゆる Doolittle 型が溢れているような気もする。なので、参考にする際は、その辺りは注意願いたい。

サンプルプログラム

　連立一次方程式

\begin{align}
Ax&=b \tag{1}
\end{align}

を解くための Python プログラムを以下に示す。
　このサンプルプログラムには以下の関数が記述されている。

関数	概要
`dot`	行列 A とベクトル x の積を返す関数
`lu_decomposition_left_1`	LU分解 (シンプルな Left looking)
`lu_decomposition_left_2`	LU分解 (一般的な Left looking)
`lu_decomposition_right`	LU分解 (Right looking)
`solve_lu`	LU分解を用いて連立一次方程式を解く
`solve_gauss_jordan`	ガウス・ジョルダン法を用いて連立一次方程式を解く
`solve_gauss_elimination`	ガウス消去法を用いて連立一次方程式を解く

　なお、これらのアルゴリズムを記述しておくことが主目的のため、分かりやすさを重視してシンプルにしているので、エラー処理は最低限だし、ピボッティング (pivoting: pivot選択) などは入っていない。

solver.py

# -*- coding: utf-8 -*-

def dot(A, x):
    """
    行列 A とベクトル x の積

    Args:
      A: 行越
      x: ベクトル

    Returns:
      b: ベクトル Ax
    """
    n = len(A)
    b = [0.0] * n
    for i in range(n):
        for j in range(n):
            b[i] += A[i][j] * x[j]
    return b


def lu_decomposition_left_1(A):
    """
    LU分解 (Left looking, ijk)

    Args:
      A: 連立一次方程式の係数行列

    Returns:
      LU: 下三角&上三角行列 (Crout 型: u_ii=1)
    """

    n = len(A)
    LU = [a[:] for a in A]  # A を壊さないようにコピーしておく

    for i in range(n):
        for j in range(n):
            for k in range(min(i, j)):
                LU[i][j] -= LU[i][k] * LU[k][j]
        a = LU[i][i]
        if a == 0:
            raise ValueError("A[i][i] is 0.")
        for j in range(i+1, n):
            LU[i][j] /= a
    return LU


def lu_decomposition_left_2(A):
    """
    LU分解 (Left looking, ikj)

    Args:
      A: 連立一次方程式の係数行列

    Returns:
      LU: 下三角&上三角行列 (Crout 型: u_ii=1)
    """

    n = len(A)
    LU = [a[:] for a in A]  # A を壊さないようにコピーしておく

    for i in range(n):
        for k in range(i):
            for j in range(k+1, n):
                LU[i][j] -= LU[i][k] * LU[k][j]
        a = LU[i][i]
        if a == 0:
            raise ValueError("A[i][i] is 0.")
        for j in range(i+1, n):
            LU[i][j] /= a
    return LU


def lu_decomposition_right(A):
    """
    LU分解 (Right looking)

    Args:
      A: 連立一次方程式の係数行列

    Returns:
      LU: 下三角&上三角行列 (Crout 型: u_ii=1)
    """

    n = len(A)
    LU = [a[:] for a in A]  # A を壊さないようにコピーしておく

    for k in range(n):
        a = LU[k][k]
        if a == 0:
            raise ValueError("A[i][i] is 0.")
        for j in range(k+1, n):
            LU[k][j] /= a
        for i in range(k+1, n):
            for j in range(k+1, n):
                LU[i][j] -= LU[i][k] * LU[k][j]
    return LU


def solve_lu(A, b, lu_decomposition):
    """
    LU分解を用いて連立一次方程式を解く。

    Args:
      A: 連立一次方程式の係数行列
      b: 連立一次方程式の右辺

    Returns:
      x: 連立一次方程式の解
    """

    LU = lu_decomposition(A)
    n = len(LU)
    x = b[:]
    for i in range(n):
        for k in range(i):
            x[i] -= LU[i][k] * x[k]
        x[i] /= LU[i][i]

    for i in reversed(range(n)):
        for k in range(i + 1, n):
            x[i] -= LU[i][k] * x[k]

    return x


def solve_gauss_jordan(A, b):
    """
    ガウス・ジョルダン法を用いて連立一次方程式を解く。

    Args:
      A: 連立一次方程式の係数行列
      b: 連立一次方程式の右辺

    Returns:
      x: 連立一次方程式の解
    """

    n = len(A)
    x = b[:]
    K = [a[:] for a in A]  # A を壊さないようにコピーしておく

    for i in range(n):
        a = K[i][i]
        if a == 0:
            raise ValueError("A[i][i] is 0.")
        for j in range(i, n):
            K[i][j] /= a
        x[i] /= a
        for ii in range(n):
            if i == ii:
                continue
            a = K[ii][i]
            for j in range(i, n):
                K[ii][j] -= a * K[i][j]
            x[ii] -= a * x[i]
    return x


def solve_gauss_elimination(A, b):
    """
    ガウス消去法を用いて連立一次方程式を解く。

    Args:
      A: 連立一次方程式の係数行列
      b: 連立一次方程式の右辺

    Returns:
      x: 連立一次方程式の解
    """

    n = len(A)
    x = b[:]
    K = [a[:] for a in A]  # A を壊さないようにコピーしておく

    for i in range(n):
        a = K[i][i]
        if a == 0:
            raise ValueError("A[i][i] is 0.")
        for j in range(i, n):
            K[i][j] /= a
        x[i] /= a

        for ii in range(i + 1, n):
            a = K[ii][i]
            for j in range(i, n):
                K[ii][j] -= a * K[i][j]
            x[ii] -= a * x[i]

    for i in reversed(range(n)):
        for j in range(i + 1, n):
            x[i] -= K[i][j] * x[j]

    return x


if __name__ == '__main__':
    A = [[2., 1., 1., -1.],
         [1., 2., -1., 2.],
         [0., 1., 2., -2.],
         [-2., 1., 0.,  3.]]
    x = [1., -2., 13., -4.]
    b = dot(A, x)
    print('A=', A)
    print('x=', x)
    print('b=', b)

    x = solve_lu(A, b, lu_decomposition_left_1)
    print('x(LU L1)=', x)

    x = solve_lu(A, b, lu_decomposition_left_2)
    print('x(LU L2)=', x)

    x = solve_lu(A, b, lu_decomposition_right)
    print('x(LU R )=', x)

    x = solve_gauss_jordan(A, b)
    print('x(Gss J)=', x)

    x = solve_gauss_elimination(A, b)
    print('x(Gss E)=', x)

　上記のプログラムには、以下の行列 $A$ およびベクトル $b$ における答え $x$ を求めるテスト例題が入っている。

\begin{align}
A&=
\begin{pmatrix}
2 & 1 & 1 & -1\\
1 & 2 & -1 & 2\\
0 & 1 & 2 & -2\\
-2 & 1 & 0 & 3
\end{pmatrix} \tag{2}\\
b&=\begin{pmatrix}
17\\ -24\\ 32\\ -16
\end{pmatrix} \tag{3}
\end{align}

　答え $x$ は以下になる。

x=\begin{pmatrix}
1 \\ -2 \\ 13 \\ -4
\end{pmatrix} \tag{4}

$ python solver.py
A= [[2.0, 1.0, 1.0, -1.0], [1.0, 2.0, -1.0, 2.0], [0.0, 1.0, 2.0, -2.0], [-2.0, 1.0, 0.0, 3.0]]
x= [1.0, -2.0, 13.0, -4.0]
b= [17.0, -24.0, 32.0, -16.0]
x(LU L1)= [0.9999999999999996, -1.9999999999999991, 13.0, -4.000000000000001]
x(LU L2)= [0.9999999999999996, -1.9999999999999991, 13.0, -4.000000000000001]
x(LU R )= [0.9999999999999996, -1.9999999999999991, 13.0, -4.000000000000001]
x(Gss J)= [1.0000000000000018, -2.0000000000000004, 13.0, -4.000000000000001]
x(Gss E)= [0.9999999999999996, -1.9999999999999991, 13.0, -4.000000000000001]

　五通り（Left lookingのLU分解が二つ、Right lookingのLU分解、ガウスジョルダン法、ガウス消去法）での答えが出力されるが、どれも正しく計算されていることが確認できる。
　（多少の数値誤差はご愛嬌、かな）

LU分解について

　いろいろ検索してみると、LU分解について小難しい理論で解説されていることが多い気がするので、ここではざっくりとした、直感的な解説をしてみる。（私なりの）

　LU分解とは、$N×N$ の行列 $A$ を下三角行列 $L$ と上三角行列 $U$ に分解することである。それはまさに言葉通りであるが、これを式で表すと以下のように書ける。

A = LU \tag{5}

　これを行列表記で表すならば以下のようになる。

\begin{matrix}
\begin{pmatrix}
a_{11} & \cdots & a_{1n} \\
\vdots & \ddots & \vdots \\
a_{n1} & \cdots & a_{nn} \\
\end{pmatrix}\\A
\end{matrix}
=
\begin{matrix}
\begin{pmatrix}
l_{11} & & \text{0} \\
\vdots & \ddots & \\
l_{n1} & \cdots & l_{nn}
\end{pmatrix}\\L
\end{matrix}
\begin{matrix}
\begin{pmatrix}
u_{11} & \cdots & u_{1n} \\
& \ddots & \vdots \\
\text{0} & & u_{nn} \\
\end{pmatrix}\\U
\end{matrix} \tag{6}

　このとき行列 $A$ の成分 $a_{ij}$ は以下のように書ける。

\begin{eqnarray}
a_{ij} &=& \sum_{k}^{\min(i,j)} l_{ik} u_{kj} \tag{7}
\end{eqnarray}

　具体的な例（例えば $a_{53}$）を挙げるならば以下のようになる。

\begin{eqnarray}
a_{53} &=& l_{51} u_{13} + l_{52} u_{23} + l_{53} u_{33} \tag{8}
\end{eqnarray}

　ここで $u_{33}=1$ 、すなわち上三角行列 $U$ の対角成分を $1$ と置くと、

\begin{eqnarray}
l_{53} &=& a_{53} - (l_{51} u_{13} + l_{52} u_{23}) \tag{9}
\end{eqnarray}

　と書ける。すなわち下三角行列 $L$ の成分 $l_{ij}$ は、

l_{ij} = a_{ij} - \sum_{k}^{j-1} l_{ik} u_{kj} \label{lij}\tag{10}

　となる。同様にもう一つ例（$a_{35}$）を挙げると以下になる。

a_{35} = l_{31} u_{15} + l_{32} u_{25} + l_{33} u_{35} \tag{11}

　ここから $u_{35}$ を求めると以下になる。

u_{35} = \frac{1}{l_{33}} \left(a_{35} - (l_{31} u_{15} + l_{32} u_{25}) \right) \tag{12}

　すなわち、上三角行列 $U$ の成分 $u_{ij}$ は以下となる。

u_{ij} = \frac{1}{l_{ii}} \left( a_{ij} - \sum_{k}^{i-1} l_{ik} u_{kj} \right) \label{uij}\tag{13}

　あとは、これらに従って下三角行列 $L$ と上三角行列 $U$ を求めていけばよい。が、ここで改めて式 $(\ref{lij})$ と式 $(\ref{uij})$ を確認してみる。

\left\{\begin{array}{ll}
l_{ij} = a_{ij} - \sum_{k}^{j-1} l_{ik} u_{kj} & (i \ge j) \\
u_{ij} = \frac{1}{l_{ii}} \left( a_{ij} - \sum_{k}^{i-1} l_{ik} u_{kj} \right) & (i \lt j)
\end{array}\right.\tag{14}

　これをもう少しまとめてみる。

\begin{align}
lu_{ij} &= t\left( a_{ij} - \sum_{k}^{\min(i,j)-1} l_{ik} u_{kj} \right)\\
t &=\left\{\begin{array}{ll}
1 & (i \ge j) \\
\frac{1}{l_{ii}} & (i \lt j)
\end{array}\right.\tag{15}
\end{align}

　上記は当然ながら上三角行列 $U$ の対角成分を $1$ とした場合であり、Crout 型と呼ばれる。逆に下三角行列 $L$ の対角成分を $1$ とした場合は、

\begin{align}
t &=\left\{\begin{array}{ll}
\frac{1}{u_{ii}} & (i \gt j) \\
1 & (i \le j)
\end{array}\right.\tag{15}
\end{align}

　となり、こちらは Doolittle 型と呼ばれる。

　上記の式から、添字 $i,j,k$ でループすることはすぐに分かると思う。このループを回す順番で計算効率や並列効率が変わってくる。そしてそれは配列のメモリの持ち方にも依存する。
　とはいえ、それが顕著に影響してくるのは行列 $A$ の対角数 $N$ が大きい場合の話だ。最近のコンピューターの性能を考えると、例えば 100 以上とか。
　$N$ が小さいうち、例えば 20 以下とか、ならば気にするほどの違いは現れないかと思う。
　逆に $N$ が巨大にならば、もっと別なことを考える必要が出てくるかもしれない。密行列のままでいいのか、疎行列で考えたほうがいいのか、直接法でいいのか、それとも反復法を考えるべきなのか、などなど。

終わりに

　連立一次方程式を解く基本的なアルゴリズムを示した。
　が、最初にも書いたが、最近は便利なライブラリが溢れているので、普通に計算するならばそっちを使うほうが断然良い。計算効率も段違いに良いし。

　そういう意味では、本記事は自分向けの備忘録ではあるが、それでも、もし誰かの何かの足しになるなら幸いである。