逆行列との行列積 $X^{-1}Y$ の計算方法・計算量のメモ
逆行列を直接計算しないで線型方程式を解きましょうというお話の予定

基礎事項

Linear Algebra and its Application¹ という本から計算量を引用．
$X\in\mathbb{R}^{m \times m}$, $Y\in\mathbb{R}^{m \times n}$, $Z = X^{-1}Y \in\mathbb{R}^{m \times n}$ とする．
計算量を floating point operations (flops) で測る．

• 逆行列 $X^{-1}$ を計算 → 行列積 $X^{-1}Y$ を計算すると $2m^3 + 2m^2n$ flops
• LU 分解 $X=LU$ を計算 → 線型方程式 $XZ=Y$ を解くと $\frac23m^3 + 2m^2n$ flops

線型方程式を解いたほうが速そう

Python の実装比較

線型方程式を扱うソルバーとして，numpy.linalg.solve と scipy.linalg.solve があるのでその辺りを比較する．

使用モジュールはこんな感じ．

import numpy as np
import numpy.linalg as la
import scipy.linalg as sla

計算時間は 10 回平均を使用．
行列の次元は $m=n$ として，10 ~ 1,000 を取るよう生成した．

dims = np.round(np.logspace(1, 3, 100)).astype(int)

一般の行列 X, Y

def generate(n):
    return np.random.randn(n, n)

で行列を生成する．比較する計算式は以下の通り．

z = la.inv(x) @ y
z = la.solve(x, y)

np.linalg.solve のほうが速い．

正定値対称行列 X

def generate(n):
    x = np.random.randn(n, n)
    return x @ x.T

で正定値対称行列を生成する．比較する計算式は以下の通り．

z = la.solve(x, y)
z = sla.solve(x, y, assume_a='pos')

scipy.linalg.solve のパラメータ assume_a を 'pos' とすることで正定値対称行列に特化したソルバーを使ってくれるらしい．
numpy で良さそう．

3 つの行列積

$X$, $Y$ がともに正定値対称であるとき，$X^{-1}YX^{-1}$, $YX^{-1}Y$ の計算を考える．
scipy.linalg.solve_triangular で上 (下) 三角行列の線型方程式を効率的に解くことができる．

• $X^{-1}YX^{-1}$ の計算
  
  • 逆行列 $X^{-1}$ を計算 → 行列積 $X^{-1}YX^{-1}$ を計算すると $6m^3$ flops
  • Cholesky 分解 $X=LL^\top$ を計算 → 線型方程式 $XZ=Y$ を解き $Z = X^{-1}Y$ を得る → 線型方程式 $XW = Z^\top$ を解き $W = X^{-1}YX^{-1}$ を得ると $\left(\frac13 + 2 + 2\right)m^3 \approx 4.33m^3$ flops
  • Cholesky 分解 $X=L_XL_X^\top$ を計算 → Cholesky 分解 $Y=L_YL_Y^\top$ を計算 → 線型方程式 $XZ = L_Y$ を解き $Z = X^{-1}L_Y$ を得る → 行列積 $ZZ^\top = X^{-1}YX^{-1}$ を計算すると $\left(\frac13 + \frac13 + 2 + 2\right)m^3 \approx 4.67m^3$ flops

# ひとつめ
xinv = la.inv(x)
z = xinv @ y @ xinv

# ふたつめ (1)
z = sla.solve(x, y, assume_a='pos')
z = sla.solve(x, z.T, assume_a='pos')

# ふたつめ (2)
xc = la.cholesky(x)
w = sla.solve_triangular(xc, y, lower=True)
w = sla.solve_triangular(xc.T, w, lower=False)
w = sla.solve_triangular(xc, w.T, lower=True)
z = sla.solve_triangular(xc.T, w, lower=False)

# みっつめ
yc = la.cholesky(y)
z = sla.solve(x, yc, assume_a='pos')
z = z @ z.T

どれもそこまで差はなさそう．

• $YX^{-1}Y$ の計算
  • 逆行列 $X^{-1}$ を計算 → 行列積 $YX^{-1}Y$ を計算すると $6m^3$ flops
  • 線型方程式 $XZ = Y$ を解き $Z = X^{-1}Y$ を得る → 行列積 $YZ = YX^{-1}Y$ を計算すると $\left(\frac13 + 2 + 2\right)m^3 \approx 4.33m^3$ flops
  • Cholesky 分解 $X=LL^\top$ を計算 → 線型方程式 $LZ=Y$ を解き $Z = L^{-1}Y$ を得る → 行列積 $Z^\top Z = YX^{-1}Y$ を計算すると $\left(\frac13 + 1 + 2\right)m^3 \approx 3.33m^3$ flops

# ひとつめ
xinv = la.inv(x)
z = y @ xinv @ y

# ふたつめ
z = sla.solve(x, y, assume_a='pos')
z = y @ z

# みっつめ
xc = la.cholesky(x)
z = sla.solve_triangular(xc, y, lower=True)
z = z.T @ z

Cholesky 分解を用いたものが一番速い．