numpyやscipyの逆行列の実装について - データサイエンス関連+αのメモ

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

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

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

A X = B

(1)

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

A X = I

(2)

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

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

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

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

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

scipy.linalg.inv¶

docs: scipy.linalg.inv — SciPy v1.12.0 Manual
source: source

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¶

docs: numpy.linalg.inv — NumPy v1.26 Manual
source:
- numpy/numpy/linalg/linalg.py
- numpy/numpy/linalg/umath_linalg.cpp
  - gesv関数を呼び出している

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)

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}

(3)

は、行列式の計算で $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)}

(4)

で、 $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]])

参考¶

幸谷智紀. (2016). LAPACK/BLAS 入門: Linear Algebra PACKage Basic Linear Algebra Subprograms.
渡部善隆「連立1次方程式の直接解法とソフトウェア」
行列計算における高速アルゴリズム