Numpy & Scipy - 1.6 The Solution of Matrix Equation (General Matrices)

行列式 (Determinant) の計算

from scipy import linalg

• scipy の linalg をインポートしよう。

• linalg.det(Matrix) を使うと行列式を計算して返してくれる。
• 基本的なアルゴリズムは LU 分解 (LU Decomposition) を使用する。$\quad A = LU$
• また、行列式の計算は以下のように実行される。

\[det(A) = det(L*U) = det(L)*det(U) = det(U)\]

• $L$ の対角成分はすべて 1 であるため $det(L) = 1$ となり、$det(U)$ だけを計算すれば良い。
• したがって、一般的なケース (General Case) の行列式を求める際には、線形代数で学んだクラメールの公式 (Cramer’s Rule) のような方法は使用しない。

A1 = np.array([[1,5,0],[2,4,-1],[0,-2,0]])
det = linalg.det(A1)

print(det)

A2 = np.array([[1,-4,2],[-2,8,-9],[-1,7,0]])
det = linalg.det(A2)

print(det)

-2.0

15.0

逆行列より Ax=b を解く方がよい理由

$A^{-1}$ を求めて解を得る場合

$A^{-1}$ の計算量 : $\sim n^3$ $A^{-1}b$ の計算量 : $\sim n^2$

$A\mathbf{x} = \mathbf{b}$ の方程式を解いて解を得る場合

ガウスの消去法 (Gauss Elimination) : $\sim n^3$

LU 分解 (LU Decomposition) :

• $A = LU$ に分解する計算量 : $\sim n^3$
• $LU\mathbf{x} = \mathbf{b}$ の方程式を解く計算量 : $\sim n^2$

• どちらも計算量は似ているように見えるが、以下の理由から逆行列を使わず $A\mathbf{x} = \mathbf{b}$ の方程式を解く方法が好まれる。

数値精度 (Numerical Accuracy)

コンピュータでは浮動小数点 (floating point)、すなわち限られた精度で数値が表現される。逆行列を使って計算すると、方程式を直接解く方法に比べて近似誤差が多くなり、不正確な結果が生じる。また、後退誤差解析 (backward error analysis) の観点からも逆行列の使用は避けられることが多い。

疎行列 (Sparse Matrix) の逆行列と LU 分解の比較

行列のほとんどの要素が 0 の場合 (疎行列)、その逆行列はどうなるだろうか？

• 疎行列の逆行列もほぼ疎 (sparse) ではないかと考えたが、現実はそう甘くない。

一方、L と U は依然として疎 (sparse) であることが多い。

• まとめると、特に $A^{-1}$ の逆行列が必要な場合でない限り、$A\mathbf{x} = \mathbf{b}$ を解く際に逆行列を求めないようにしよう！

逆行列 (Inverse) の計算

• 逆行列を求める目的が $A\mathbf{x} = \mathbf{b}$ を解くことであれば、もう一度考え直す必要がある。
• linalg.inv(matrix) を使って逆行列を求める。

• ここでパラメータとして入力した行列の逆行列が存在しない場合、すなわち 特異行列 (singular matrix) であればエラーが発生する。

• 基本的なアルゴリズムは LU 分解を使用し $\quad LUA^{-1} = I$ を解くが、後退代入 (backward phase / back substitution) の過程のみを経る。

A1 = np.array([[1,2,1],[2,1,3],[1,3,1]])
inv_a1 = linalg.inv(A1)

prt(inv_a1, fmt="%0.2f")

 8.00, -1.00, -5.00
-1.00,  0.00,  1.00
-5.00,  1.00,  3.00

行列方程式 Ax = b を解く

• $A\mathbf{x} = \mathbf{b}$ を解くこと、すなわち解を求めることは $\mathbf{x}$ を求めることと同じであるため、

• x = linalg.solve(A, b, assume_a="gen") 関数を使えば良い。
• ここで assume_a に入力できるパラメータは 4 種類ある。

gen 汎用ソルバー。行列のタイプ (対称、エルミートなど) が不明な場合に主に使用する。 LU 分解で解く。$\quad A = LU$

sym 対称行列 (symmetric matrix) または複素対称行列 (complex symmetric matrix) に使用できる。複素対称はエルミート (Hermitian) とは異なることに注意しよう。 対称行列の性質は $\; A = A^T \;$ である。 対角ピボット (diagonal pivoting) $\quad A = LDL^T \mbox{(block diagonal)}$ を使用する。

her エルミート行列 (Hermitian matrix) に使用する。$\quad A = A^$ (共役転置) 対角ピボット (diagonal pivoting) $\quad A = LDL^$ で解く。

pos 正定値行列 (positive definite matrix) に使用する。二次形式 (quadratic form) $\mathbf{x}^TA\mathbf{x} > 0$ を意味する。 アルゴリズムはコレスキー分解 (Cholesky Decomposition) $\quad A = R^TR = LL^T$ を使用する。

• ここで注意すべき点は、assume_a を誤って設定してもエラーは表示されず、誤った結果が返ってくることに気をつけること。そのため通常は gen を主に使用する。
• 行列の性質をよく理解していれば assume_a オプションを使用すると良い。(演算速度の面でメリットがある)

b = np.ones((3,), dtype=np.float64)

# singular (not invertible)
A_sing = np.array([[1,3,4],[-4,2,-6],[-3,-2,-7]])

# gen (general)
A_gen = np.array([[0,1,2],[1,0,3],[4,-3,8]])

# sym (symmetric)
A_sym = np.array([[1,2,1],[2,1,3],[1,3,1]])

# sym complex (symmetric complex)
A_sym_c = np.array([[1, 2-1j, 1+2j],[2-1j, 1, 3],[1+2j, 3,1]])

# her (hermitian)
A_her = np.array([[1, 2+1j, 1-2j],[2-1j, 1, 3],[1+2j, 3,1]])

# pos (positive definite)
A_pos = np.array([[2, -1, 0],[-1,2,-1],[0,-1,2]])

• singular
• 特異行列 (singular matrix) であるためエラーが発生する。

x = linalg.solve(A_sing, b)

numpy.linalg.LinAlgError: Matrix is singular.

• gen
• $Ax = b \rightarrow Ax - b = 0$ であることを確認するために最後の行を追加している。

x = linalg.solve(A_gen, b)

prt(x, fmt="%0.5e")
print()
prt(A_gen @ x - b , fmt="%0.5e")

 1.00000e+00,  1.00000e+00,  0.00000e+00

 0.00000e+00,  0.00000e+00,  0.00000e+00

• sym
• sym と gen の違いは解法アルゴリズムにある。
• 求められる x は近似値であることを念頭に置いておこう。gen による LU 分解と sym による対角ピボットで計算した結果には誤差が生じる。

x1 = linalg.solve(A_sym, b) # default 로는 gen
x2 = linalg.solve(A_sym, b, assume_a="sym")

prt(x1, fmt="%0.5e")
print()
prt(A_sym @ x1 - b , fmt="%0.5e")

prt(x2, fmt="%0.5e")
print()
prt(A_sym @ x2 - b , fmt="%0.5e")

 # gen
 2.00000e+00,  0.00000e+00, -1.00000e+00

 0.00000e+00,  0.00000e+00,  0.00000e+00

# sym
 2.00000e+00, -2.77556e-17, -1.00000e+00

 0.00000e+00,  0.00000e+00,  0.00000e+00

• her

x = linalg.solve(A_her, b, assume_a="her")

prt(x, fmt="%0.5e")
print()
prt(A_her @ x - b , fmt="%0.5e")

( 1.11111e-01+1.11111e-01j), ( 3.33333e-01-1.11111e-01j), ( 1.11111e-01+1.11022e-16j)

( 0.00000e+00+2.77556e-17j), (-4.44089e-16+1.94289e-16j), ( 2.22045e-16-1.66533e-16j)

• ここで e-16 は $10^{-16}$ を意味するため、0 に近い値である。

• pos

x1 = linalg.solve(A_pos, b, assume_a="gen")
x2 = linalg.solve(A_pos, b, assume_a="pos")

prt(x1, fmt="%0.5e")
print()
prt(A_pos @ x1 - b , fmt="%0.5e")

prt(x2, fmt="%0.5e")
print()
prt(A_pos @ x2 - b , fmt="%0.5e")

# gen
 1.50000e+00,  2.00000e+00,  1.50000e+00

 0.00000e+00,  2.22045e-16, -4.44089e-16

# pos
 1.50000e+00,  2.00000e+00,  1.50000e+00

-8.88178e-16,  1.55431e-15, -4.44089e-16

三角行列ソルバー (Triangular Matrix Solver)

• 同様に $A\mathbf{x} = \mathbf{b}$ を解くが、行列 $A$ が下三角行列 (lower triangular matrix) または上三角行列 (upper triangular matrix) の場合にのみ使用できる。
• 後退代入 (backward phase / back substitution) の計算のみが必要である。
• x = linalg.solve_triangular(A, b, lower=False)
• lower パラメータは下三角行列の場合は True、上三角行列の場合は False を指定する。

A = np.array([[1,0,0,0],[1,4,0,0],[5,0,1,0],[8,1,-2,2]], dtype=np.float64)
b = np.array([1,2,3,4],dtype=np.float64)

x = linalg.solve_triangular(A, b, lower=True)

prt(x, fmt="%0.5e")

 1.00000e+00,  2.50000e-01, -2.00000e+00, -4.12500e+00

求めた解は正確か？

• $A\mathbf{x} = \mathbf{b}$ を解いて得た $\mathbf{x}$ の解は、数値計算による近似値である。

• $A\mathbf{x}$ と $\mathbf{b}$ が十分に近いか？ == $A\mathbf{x} - \mathbf{b}$ が十分に 0 に近いか？ と同じ問いである。

• したがって、np.allclose 関数を使い、$A\mathbf{x} - \mathbf{b}$ と np.zeros で作成したゼロ行列を比較すれば良い。

A = np.array([[2,-1,0],[-1,2,-1],[0,-1,2]], dtype=np.float64)
b = np.ones((3,), dtype=np.float64)

x = linalg.solve(A, b, assume_a="pos")

prt(x, fmt="%0.5e")
print()
prt(A@x-b, fmt="%0.5e")

zr = np.zeros((3,), dtype=np.float64)

bool_close = np.allclose(A@x-b, zr)
print(bool_close)

 1.50000e+00,  2.00000e+00,  1.50000e+00

-8.88178e-16,  1.55431e-15, -4.44089e-16
True

例題 1.

$A = \begin{bmatrix} 1 & 1/2 & 1/3 & \dots & 1/9 & 1/10 \\ 1/2 & 1/3 & 1/4 & \dots & 1/10 & 1/11 \\ 1/3 & 1/4 & 1/5 & \dots & 1/11 & 1/12 \\ \vdots & \ddots & \ddots & & & \\ 1/9 & 1/10 & & & 1/17 & 1/18 \\ 1/10 & 1/11 & \dots & \dots & 1/18 & 1/19 \end{bmatrix} $

linalg.hilbert(10) を使って 10x10 のヒルベルト行列 (Hilbert matrix) を作成し、A に格納する。 A の逆行列を inv_A に格納する。 x1 = inv_A @ b で求める。 x2 を linalg.solve を使って求める (gen)。 A@x1 - b を小数点以下 15 桁まで出力する (floating)。 A@x2 - b を小数点以下 15 桁まで出力する (floating)。 A@x1 - b と A@x2 - b を allclose を使ってゼロベクトルと比較する。

• 解答
• ヒルベルト行列 (Hilbert matrix) は数値解析的に非常に扱いが難しい行列である。サイズが大きくなるほど解くのが困難になる。
• 理論として学んだが、逆行列の形で解くのは精度の面で劣ることが確認できる。

A = linalg.hilbert(10)
b = np.ones((10,), dtype=np.float64)

# Hilbert 행렬
prt(A, fmt="%0.5f")

A_inv = linalg.inv(A)

x1 = A_inv @ b

print()

x2 = linalg.solve(A, b, assume_a="gen")

prt(A @ x1 - b, fmt="%0.15e")
print()
prt(A @ x2 - b, fmt="%0.15e")
print()

zr = np.zeros((10,), dtype=np.float64)

bool_close1 = np.allclose(A @ x1 - b, zr)
print(bool_close1)

bool_close2 = np.allclose(A @ x2 - b, zr)
print(bool_close2)

• ヒルベルト行列の出力結果

  1.00000,  0.50000,  0.33333,  0.25000,  0.20000,  0.16667,  0.14286,  0.12500,  0.11111,  0.10000
   0.50000,  0.33333,  0.25000,  0.20000,  0.16667,  0.14286,  0.12500,  0.11111,  0.10000,  0.09091
   0.33333,  0.25000,  0.20000,  0.16667,  0.14286,  0.12500,  0.11111,  0.10000,  0.09091,  0.08333
   0.25000,  0.20000,  0.16667,  0.14286,  0.12500,  0.11111,  0.10000,  0.09091,  0.08333,  0.07692
   0.20000,  0.16667,  0.14286,  0.12500,  0.11111,  0.10000,  0.09091,  0.08333,  0.07692,  0.07143
   0.16667,  0.14286,  0.12500,  0.11111,  0.10000,  0.09091,  0.08333,  0.07692,  0.07143,  0.06667
   0.14286,  0.12500,  0.11111,  0.10000,  0.09091,  0.08333,  0.07692,  0.07143,  0.06667,  0.06250
   0.12500,  0.11111,  0.10000,  0.09091,  0.08333,  0.07692,  0.07143,  0.06667,  0.06250,  0.05882
   0.11111,  0.10000,  0.09091,  0.08333,  0.07692,  0.07143,  0.06667,  0.06250,  0.05882,  0.05556
   0.10000,  0.09091,  0.08333,  0.07692,  0.07143,  0.06667,  0.06250,  0.05882,  0.05556,  0.05263
  

• 逆行列による解と np.solve (gen) による解の比較

-5.384686879750245e-05, -4.824035864248177e-05, -4.365437722742005e-05, -3.984553430069759e-05, -3.663714443313815e-05, -3.390060580654719e-05, -3.154054138576612e-05, -2.948524204493541e-05, -2.767985240692550e-05, -2.608172012175114e-05

 1.768403201651836e-10,  2.589735093039280e-10,  1.027888885118955e-11, -5.332401187274627e-13, -2.094682205466825e-10, -1.216354794664198e-10,  8.731149137020111e-11,  4.546962806273314e-11,  5.050826423769195e-11, -6.052347512053302e-11

False
True