【Python】逆行列との行列積 $X^{-1}Y$ の計算

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

# 基礎事項
Linear Algebra and its Application[^1] という本から計算量を引用．
$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) で測る．

|  | flops |
|:---:|---:|
| 行列積 $XY$ の計算 | $2m^2n$ |
| 逆行列 $X^{-1}$ の計算 | $2m^3$ |
| LU 分解 $X=LU$ の計算 | $\frac23m^3$ |
| 線型方程式 $LW=Y$ の計算 | $m^2 \times n$ |
| 線型方程式 $UZ=W$ の計算 | $m^2 \times n$ |
| Cholesky 分解[^2] $X=LL^\top$ の計算 | $\frac13m^3$ |


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

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

# Python の実装比較
線型方程式を扱うソルバーとして，[`numpy.linalg.solve`](https://numpy.org/doc/stable/reference/generated/numpy.linalg.solve.html) と [`scipy.linalg.solve`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.solve.html) があるのでその辺りを比較する．

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

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

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

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

## 一般の行列 X, Y
```python
def generate(n):
    return np.random.randn(n, n)
```
で行列を生成する．比較する計算式は以下の通り．

```python
z = la.inv(x) @ y
z = la.solve(x, y)
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/394475/9e0b2393-d610-2e36-939b-e032b4a8a4cc.png)

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

## 正定値対称行列 X
```python
def generate(n):
    x = np.random.randn(n, n)
    return x @ x.T
```
で正定値対称行列を生成する．比較する計算式は以下の通り．

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

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/394475/0d89d45a-60df-1662-70cc-a4baee00a819.png)

`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
  
```python
# ひとつめ
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
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/394475/0c5a8853-96c6-8b89-52fe-37147f67a745.png)

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

* $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
  
```python
# ひとつめ
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
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/394475/046a87f8-5a95-4994-3059-bbcf604fa2da.png)

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


[^1]: [David C. Lay, Steven R. Lay and Judi J. McDonald, Linear Algebra with its Applications, 5th edition, Pearson, 2016](https://www.amazon.co.jp/Linear-Algebra-Its-Applications-5th/dp/032198238X)
[^2]: Cholesky 分解は正定値対称行列にのみ有効