numpyやscipyの逆行列の実装について

numpyやscipyの逆行列の実装について#

  • numpyやscipyの内部ではLAPACK（C/C++の線形代数ライブラリ）を使っている

  • LAPACKはだいたいLU分解をしている

逆行列は連立一次方程式の問題に帰着できる#

逆行列\(A^{-1}\)を求めたい場合

\[ A X = B \]

という連立一次方程式を、\(B=I\)として

\[ A X = I \]

とすることで\(X=A^{-1}\)を得ることができる

LAPACKで連立一次方程式を解く#

正則な係数行列をもつ連立一次方程式をLAPACKで解くには、xGESV関数（連立一次方程式を解く関数）もしくはxGETRF関数（LU分解の関数）とxGETRS関数（前進代入・後退代入を行う関数）を使えば解ける（幸谷, 2016）

xGESV関数は\(Ax=b\)を\(x\)について解き、\(b:=x\)として\(b\)に代入する関数

逆行列計算の実装を覗く#

scipy.linalg.inv#

scipyのinvはGETRF、GETRI

# scipyのinv
import numpy as np
from scipy import linalg
A = np.array([
    [1, 2],
    [3, 4]
])
linalg.inv(A)

array([[-2. ,  1. ],
       [ 1.5, -0.5]])

numpy.linalg.inv#

import numpy as np
np.linalg.inv(A)

array([[-2. ,  1. ],
       [ 1.5, -0.5]])

scipyのLAPACK functions#

Low-level LAPACK functions (scipy.linalg.lapack) — SciPy v1.12.0 Manual

xGESV#

DGESV computes the solution to a real system of linear equations

A * X = B,

where A is an N-by-N matrix and X and B are N-by-NRHS matrices.

The LU decomposition with partial pivoting and row interchanges is used to factor A as

A = P * L * U,

where P is a permutation matrix, L is unit lower triangular, and U is upper triangular. The factored form of A is then used to solve the system of equations A * X = B.

from scipy.linalg.lapack import dgesv

# AX = I としてX=A^{-1}を推定する
I = np.eye(A.shape[0])

lu, piv, x, info = dgesv(A, I)

x

array([[-2. ,  1. ],
       [ 1.5, -0.5]])

xGETRF#

LU分解をする関数。

scipy.linalg.lapack.dgetrf — SciPy v1.12.0 Manual

docs: LAPACK: dgetrf

DGETRF computes an LU factorization of a general M-by-N matrix A using partial pivoting with row interchanges. The factorization has the form

A = P * L * U

where P is a permutation matrix, L is lower triangular with unit diagonal elements (lower trapezoidal if m > n), and U is upper triangular (upper trapezoidal if m < n).

from scipy.linalg.lapack import dgetrf

lu, piv, info = dgetrf(A)

lu

array([[3.        , 4.        ],
       [0.33333333, 0.66666667]])

piv

array([1, 1], dtype=int32)

xGETRI#

LU分解の結果から逆行列を取得する関数

scipy.linalg.lapack.dgetri — SciPy v1.12.0 Manual

cpp docs: LAPACK: dgetri

DGETRI computes the inverse of a matrix using the LU factorization computed by DGETRF.

This method inverts U and then computes inv(A) by solving the system inv(A)*L = inv(U) for inv(A).

from scipy.linalg.lapack import dgetri

inv_a, info = dgetri(lu, piv)
inv_a

array([[-2. ,  1. ],
       [ 1.5, -0.5]])

LAPACKはどうやって計算量を下げているのか#

掃き出し法とLU分解+逆行列の計算は計算量的にほぼ一緒で\(O(n^3)\)程度らしい （『数値計算の常識』、安易に逆行列を数値計算するのはやめよう ）

しかし、一般の実数行列につかえるアルゴリズムの中では最善の様子。

LUは疎行列に強いというのもある

逆行列の計算量#

行列式を使う場合#

\(A\)を\(n\)次の正方行列とする。逆行列

\[ A^{-1} = \frac{1}{|A|} \tilde{A} \]

は、行列式の計算で\(O(n^3)\)、

行列式の計算量

行列式は

\[ |A| = \sum_{\sigma \in S_n} \operatorname{sgn}(\sigma) a_{1 \sigma(1)} a_{2 \sigma(2)} \cdots a_{n \sigma(n)} \]

で、\(S_n\)が\(n!\)あり、各項で\(n\)回\(a_{n \sigma(n)}\)を掛けるので\(O(n\times n!)\)

「行列Aのi行から、j行の定数倍を引いても、行列式の値は変わらない。」という性質を使って上三角行列に変形してから行列式を計算する。三角行列への変換は\(O(n^3)\)、上三角行列の行列式は対角成分の積なので\(O(n)\)なので、全体で\(O(n^3)\)になる

（参考：Tech Tips: 行列式の計算）

余因子行列の計算量

\((i, j)\)余因子は「『\(i\)行目と\(j\)列目を除いた行列』の行列式に\((-1)^{i+j}\)をかけたもの」なので、\(n-1\)次行列の行列式を求める計算になり、\(O([n-1]^3)\)

余因子行列は「\((i, j)\)余因子を\(i,j\)成分に持つ行列を転置したもの」なので、\((i, j)\)余因子を\(n\times n\)個計算する必要がある→\(O(n^2 \times [n-1]^3)\)か？

import numpy as np

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

np.linalg.det(A)

-2.0000000000000004

A[1, :] += (-3) * A[0, :] 
A

array([[ 1,  2],
       [ 0, -2]])

参考#