今回は、線形方程式の数値計算法の QR 分解 , コレスキー分解 を紹介します。

Python で線形連立方程式を QR 分解, コレスキー分解で解く

QR 分解

QR 分解は LU 分解と同じで線形系の行列方程式を高速に解くための手法で、係数行列 \(\mathbf{A}\) を実直行行列 \(\mathbf{Q}\) と上三角行列 \(\mathbf{R}\) の積に分解して計算します。 LU 分解よりも計算量が多いですが数値計算の安定性がよりよいのが特徴です。とはいえ一般的な行列では LU 分解で十分安定して解けるので、大規模な行列で特異性が心配になった場合に使ってみたらいいかと思います。

簡単に式の変形を示すと、

\begin{equation} \mathbf{A}\boldsymbol{x} = \boldsymbol{b} \nonumber \end{equation}

に \(\mathbf{A} = \mathbf{Q}\mathbf{R}\) を代入すると

\begin{equation} \mathbf{Q}\mathbf{R}\boldsymbol{x} = \boldsymbol{b} \nonumber \end{equation}

となり、\(\mathbf{Q}\) を右辺に移動させて以下のようになります。

\begin{equation} \mathbf{R}\boldsymbol{x} = \mathbf{Q}^{-1}\boldsymbol{b} = \mathbf{Q}^{T}\boldsymbol{b} \nonumber \end{equation}

ここで、\(\mathbf{Q}\) は実直交行列なので \(\mathbf{Q}^{-1}=\mathbf{Q}^{T}\) です。

\(\mathbf{R}\) は上三角行列なので、ガウスの消去法の後退代入が使えて解が求まります(逆行列を計算する必要がない)。やはり SciPy に QR 分解のメソッドがあるので、プログラムはそれを用いて以下のように書くことができます。

import scipy.linalg as linalg
import numpy as np

A = np.array([[6, 4, 1],
              [1, 8, -2],
              [3, 2, 0]])
b = np.array([7, 6, 8])

Q, R = linalg.qr(A)
t = np.dot(Q.T, b)
x = linalg.solve(R, t)

print(x)

[ 4. -2. -9.]

コレスキー分解

コレスキー分解も線形系の方程式を高速に解くための手法です。

LU 分解よりも高速で、使用するメモリも少ない手法ですが、\(\mathbf{A}\) がエルミート行列(または実対称行列)で正定値である場合に使えます。

係数行列 \(\mathbf{A}\) を下三角行列 \(\mathbf{L}\) とその共役転置 \(\mathbf{L}^*\) の積に分解して計算します。

\begin{equation} \mathbf{A}=\mathbf{L}\mathbf{L}^{*} \nonumber \end{equation}

と分解され、方程式は以下となります。

\begin{equation} \mathbf{L}\mathbf{L}^{*}\boldsymbol{x}=\boldsymbol{b} \nonumber \end{equation}

NumPy にコレスキー分解のメソッドがあるので、プログラムはそれを用いて以下のように書くことができます。

import numpy as np

A = np.array([[7, -1, 0, 1],
              [-1, 9, -2, 2],
              [0, -2, 8, -3],
              [1, 2, -3, 10]])
b = np.array([-5, 15, -10, 20])

L = np.linalg.cholesky(A)

t = np.linalg.solve(L, b)
x = np.linalg.solve(L.T.conj(), t)

print(x)

[-0.80693816  1.11613876 -0.3092006   1.76470588]