線形方程式を解く ¶

本稿では SciPy を用いて線形方程式、つまり連立一次方程式を数値計算で解く方法について記す。 SciPy というよりは NumPy だけで事足りるテーマかもしれないが、勢いでここに書く。

とにかく解く ¶

線形代数計算の専門書によると、Cramer の公式がどうだの、Gauss の消去法がどうだの、 LU 分解がどうだのという議論に終始するものだ。しかしここではもう解が求まればよいことにする。

コード的な手順は次のとおりとなる：

係数行列を np.array の形式で定義する。ここでは行列を A とする。
方程式右辺の列ベクトルを np.array の形式で定義する。ここでは b とする。
関数 scipy.linalg.solve を呼び出す。

これで解が存在する限り一次方程式 \(Ax = b\) を解くことができる。

次の方程式を解くためのコード例を示す。

\begin{align*} \left\{ \begin{array}{ccccc} 2x& -y& +2z& =& 8\\ x& -y& -2z& =& -1\\ -2x& +y& -z& =& -6 \end{array} \right. \end{align*}

#!/usr/bin/env python
"""linear_equation.py: Demonstrate solving linear equations by using SciPy.
"""
import numpy as np
from scipy.linalg import solve

# pylint: disable=invalid-name

A = np.array([[2, -1, 2],
              [1, -1, -2],
              [-2, 1, -1]])

b = np.array([[8],
              [-1],
              [-6]])

print("Solve Ax = b: \n", solve(A, b))

私の環境での実行結果は次のとおり（見栄えのために一部空白文字を手動で挿入）。

Solve Ax = b:
 [[ 1.]
  [-2.]
  [ 2.]]

LU 分解 ¶

線形方程式を解く計算法に有用な行列の LU 分解は、関数 scipy.linalg.lu で得られる。ただしこの関数は行列 A を \(A = PLU\) の形、つまり置換行列・下三角行列・上三角行列の積の形に分解する。

次の行列 A を LU 分解するコード例を示す。

\begin{pmatrix} 1 &2 &2\\ 2 &5 &6\\ 3 &8 &12\\ \end{pmatrix}

#!/usr/bin/env python
"""lu.py: Demonstrate LU decomposition.
"""
import numpy as np
from scipy.linalg import lu

# pylint: disable=invalid-name

A = np.array([[1, 2, 2],
              [2, 5, 6],
              [3, 8, 12]])

print("LU decomposition (A = PLU): ")
P, L, U = lu(A)
print("P: \n", P)
print("L: \n", L)
print("U: \n", U)
print("PLU: \n", P @ L @ U)

私の環境での実行結果は次のとおり（見栄えのために一部空白文字を手動で挿入）：

LU decomposition (A = PLU):
P:
 [[ 0.  1.  0.]
  [ 0.  0.  1.]
  [ 1.  0.  0.]]
L:
 [[ 1.          0.          0.        ]
  [ 0.33333333  1.          0.        ]
  [ 0.66666667  0.5         1.        ]]
U:
 [[  3.           8.          12.        ]
  [  0.          -0.66666667  -2.        ]
  [  0.           0.          -1.        ]]
PLU:
 [[  1.   2.   2.]
  [  2.   5.   6.]
  [  3.   8.  12.]]

LU 分解は一意に定まらないものなので、三角行列が期待しているようなきれいな？形にならないかもしれない。

補足 ¶

モジュール scipy.linalg の構成要素について補足する。

行列 A に対して、もし存在すれば逆行列は linalg.inv(A) または A.I で得られる。
行列 A の行列式は linalg.det で求める。
ノルムには関数 linalg.norm を用いる。ノルムの種類を引数で指示する。
最小二乗法には linalg.lstsq を用いる。これは別ページで述べる。
行列のテイラー展開による各種関数もサポート。

SymPy を用いた解法について補足する。

In [1]: solve([2*x - y + 2*z - 8, x - y - 2*z + 1, -2*x + y - z + 6], [x, y, z])
Out[1]: {x: 1, y: -2, z: 2}

In [1]: A = Matrix([[1, 2, 2], [2, 5, 6], [3, 8, 12]])

In [2]: L, U, swp = A.LUdecomposition()

In [3]: L
Out[3]:
Matrix([
[1, 0, 0],
[2, 1, 0],
[3, 2, 1]])

In [4]: U
Out[4]:
Matrix([
[1, 2, 2],
[0, 1, 2],
[0, 0, 2]])

線形方程式を解く¶

とにかく解く¶

LU 分解¶

補足¶

線形方程式を解く ¶

とにかく解く ¶

LU 分解 ¶

補足 ¶