数値解析　~Python Tips~

6209 ワード
scipy Python numpy 数値解析 Python テキストリンク

0. はじめに

数理計画の内点法を実装する際に線形方程式$Ax = b$を解くことがあり、Pythonでの実装例を調べた。本記事では線形方程式$Ax = b$を解くモジュールの使用例を記載する。

1. 線形方程式Ax=b

1.0 逆行列

素直に逆行列から演算する。

import numpy as np

A_inv = np.linalg.inv(A)
x = np.dot(A_inv,b)

1.1 LU分解

LU分解は、線形方程式を高速に解く手法である。

import scipy.linalg as linalg

LU = linalg.lu_factor(A)
x = linalg.lu_solve(LU, b)

1.2 QR分解

QR分解は、線形方程式を高速に解く手法である。

import scipy.linalg as linalg
import numpy as np

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

1.3 コレスキー分解

コレスキー分解は、線形方程式を高速に解く手法である。
行列Aが正定値対称行列であるときに適用できる。

import numpy as np

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

2　数値実験

線形方程式$Ax=b$を解いた時の各手法の実行時間を比較する。
ただし、行列$A$は正定値実対称行列である。

A = 
\left(
    \begin{array}{ccc}
      6 & 4 & 1 \\
      4 & 8 & 2 \\
      1 & 2 & 1
    \end{array}
  \right)
, \ \
b = 
\left(
    \begin{array}{c}
      7  \\
      6  \\
      8 
    \end{array}
  \right)

実行時間を下記に示す。コレスキー分解が高速に演算していることがわかる。

手法 実行時間(ms)
inv 1.24
LU分解 1.45
QR分解 0.52
コレスキー分解 0.29

補足　固有値と固有ベクトル

行列$A$の固有値$w$と固有ベクトル$v$を求めるには下記のようにすればよい。

import numpy as np

w, v = np.linalg.eig(A)

#固有値
print(w)
#[11.57822322  2.953963    0.46781378]

#固有ベクトル
print(v)
#[[ 0.59437613  0.80399525  0.01756856]
# [ 0.7780667  -0.5694105  -0.26529962]
# [ 0.20329591 -0.17135727  0.96400594]]
yaf-example