【Python】行列指数関数・行列対数関数

行列指数関数・行列対数関数の Python 上の実装の話

# 行列指数関数

$x(t): \mathbb R \to \mathbb R^m$ についての微分方程式

```math
\dot{x}(t) = Ax(t) \quad \left(A\in\mathbb C^{m\times m}\right)
```

の解は $x(t) = \mathop{\mathrm {Exp}}(tA)x(0)$ で与えられる．
ここで，$\mathop{\mathrm {Exp}}:\mathbb C^{m\times m} \to \mathbb C^{m\times m}$ は指数関数のマクローリン展開 $\exp(x) = \sum_{n=0}^\infty \frac1{n!}x^n$ と同様に，

```math
\mathop{\mathrm{Exp}}A = \sum_{n=0}^\infty \frac1{n!}A^n
```

で定義される行列指数関数である．
このように，線形システムの解などにおいて，行列指数関数の計算が必要となる．

## SciPy による実装

[scipy.linalg.expm](https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.expm.html) に行列指数関数の実装がある．
例として，次の行列 $A$ の行列指数関数 $\mathop{\mathrm{Exp}}A$ を考える．

```math
A = \left(\begin{array}{cc}
\lambda & \alpha \\
0 & \mu
\end{array}\right)
```

$A^n$ が

```math
\begin{align*}
A^n &= \left(\begin{array}{cc}
\lambda^n & \alpha \left(\lambda^{n-1} + \lambda^{n-2}\mu + \cdots + \mu^{n-1}\right) \\
0 & \mu^n
\end{array}\right)\\
&= \left(\begin{array}{cc}
\lambda^n & \alpha \frac{\lambda^n - \mu^n}{\lambda - \mu} \\
0 & \mu^n
\end{array}\right)
\end{align*}
```

となるので，行列指数関数 $\mathop{\mathrm{Exp}}A$ は

```math
\begin{align*}
\mathop{\mathrm{Exp}}A &= \left(\begin{array}{cc}
\sum_{n=0}^\infty\frac1{n!}\lambda^n & \frac{\alpha}{\lambda - \mu}\sum_{n=0}^\infty\frac1{n!}\left(\lambda^n - \mu^n\right) \\
0 & \sum_{n=0}^\infty\frac1{n!}\mu^n
\end{array}\right)\\
&= \left(\begin{array}{cc}
\exp(\lambda) & \alpha\frac{\exp(\lambda) - \exp(\mu)}{\lambda - \mu} \\
0 & \exp(\mu)
\end{array}\right)
\end{align*}
```

で計算できる．$\lambda=1, \alpha=2, \mu=3$ のときの SciPy の計算結果を以下に示す．

```python
import numpy as np
from scipy.linalg import expm


a = np.array([[1., 2.], [0., 3.]])
expa = np.array([[np.exp(1), np.exp(3) - np.exp(1)], [0., np.exp(3)]])
expm(a) - expa
# array([[0.00000000e+00, 3.55271368e-14],
#        [0.00000000e+00, 0.00000000e+00]])
```

## 対称行列の行列指数関数
$A\in\mathbb R^{m\times m}$ が対称行列の場合，直交行列 $U\in\mathbb R^{m\times m}\ \left(U^\top U = UU^\top = I\right)$ を用いて

```math
A = U\Lambda U^\top
```

と固有値分解できる．ここで，$\Lambda = \mathop{\mathrm{diag}}(\lambda_1,\lambda_2,\ldots,\lambda_m)$ は対角成分が固有値の対角行列である．このとき，

```math
\begin{align*}
A^n &= U\Lambda U^\top U\Lambda U^\top \cdots U\Lambda U^\top\\
&= U\Lambda \left(U^\top U\right)\Lambda \left(U^\top \cdots U\right)\Lambda U^\top \\
&= U \Lambda^n U^\top
\end{align*}
```

より，行列指数関数 $\mathop{\mathrm{Exp}}A$は

```math
\begin{align*}
\mathop{\mathrm{Exp}}A &= U\left(\sum_{n=0}^\infty \frac1{n!}\Lambda^n\right)U^\top \\
&= U \left(\begin{array}{cccc}
\sum_{n=0}^\infty \frac1{n!}\lambda_1^n & & &\\
& \sum_{n=0}^\infty \frac1{n!}\lambda_2^n & &\\
& & \ddots &\\
& & & \sum_{n=0}^\infty \frac1{n!}\lambda_m^n
\end{array}\right)U^\top\\
&= U \left(\begin{array}{cccc}
\exp(\lambda_1) & & &\\
& \exp(\lambda_2) & &\\
& & \ddots &\\
& & & \exp(\lambda_m)
\end{array}\right)U^\top
\end{align*}
```

で計算できる．$\exp(\lambda_i) > 0$ より，$\mathop{\mathrm{Exp}}A$ は正定値対称行列である．

SciPy による実装は非対称行列にも対応したアルゴリズムなため，対称行列に対しては上記の手法を用いたほうが速く計算できる．このアルゴリズムは多様体上の最適化ライブラリ Pymanopt に実装されている．細かい実装は [Github](https://github.com/pymanopt/pymanopt/blob/1.0.0/pymanopt/tools/multi.py) を参照．

ここでは，Scipy と Pymanopt の計算時間を比較する．

```python
import numpy as np
from scipy.linalg import expm
from pymanopt.tools.multi import multiexp


np.random.seed(0)
a = np.random.randn(10, 10)
a = 0.5 * (a + a.T)

%time expma = expm(a)
# CPU times: user 1.58 ms, sys: 1.85 ms, total: 3.44 ms
# Wall time: 3.28 ms
%time expa = multiexp(a, sym=True)
# CPU times: user 592 µs, sys: 194 µs, total: 786 µs
# Wall time: 855 µs

np.testing.assert_array_almost_equal(expma, expa)
```

# 行列対数関数

$\mathop{\mathrm{Exp}}A = X$ がみたされるとき，$A = \mathop{\mathrm{Log}}X$ を行列対数関数という．$\exp(x) = \exp(x + 2in\pi)$ であるので，$\mathop{\mathrm{Log}}$ は一意ではない．そこで，固有値の虚部が $(-\pi, \pi)$ に入るものを principal logarithm と呼ぶ．

## SciPy による実装

[scipy.linalg.logm](https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.logm.html) に行列対数関数の実装がある．
例として，行列指数関数の例で用いた行列 $A$ を用いる ($\lambda=1, \alpha=2, \mu=3$)．

```math
\begin{align*}
A &= \left(\begin{array}{cc}
\lambda & \alpha \\
0 & \mu
\end{array}\right)\\
\mathop{\mathrm{Exp}}A &= \left(\begin{array}{cc}
\exp(\lambda) & \alpha\frac{\exp(\lambda) - \exp(\mu)}{\lambda - \mu} \\
0 & \exp(\mu)
\end{array}\right)
\end{align*}
```

```python
import numpy as np
from scipy.linalg import logm


expa = np.array([[np.exp(1), np.exp(3) - np.exp(1)], [0., np.exp(3)]])
logm(expa)
# array([[1., 2.],
#        [0., 3.]])
```

## 正定値対称行列の行列対数関数
$X\in\mathbb R^{m\times m}$ が正定値対称行列の場合，直交行列 $V\in\mathbb R^{m\times m}\ \left(V^\top V = VV^\top = I\right)$ を用いて

```math
X = V\Sigma V^\top
```

と固有値分解できる．ここで，$\Sigma = \mathop{\mathrm{diag}}(\sigma_1,\sigma_2,\ldots,\sigma_m)$ は対角成分が<font color="Red">正</font>の固有値の対角行列である．このとき，

```math
\begin{align*}
X &= V \left(\begin{array}{cccc}
\exp(\log(\sigma_1)) & & &\\
& \exp(\log(\sigma_2)) & &\\
& & \ddots &\\
& & & \exp(\log(\sigma_m))
\end{array}\right)V^\top\\
&= \mathop{\mathrm{Exp}}\left(V \left(\begin{array}{cccc}
\log(\sigma_1) & & &\\
& \log(\sigma_2) & &\\
& & \ddots &\\
& & & \log(\sigma_m)
\end{array}\right)V^\top\right)
\end{align*}
```

とかけるので，行列対数関数 $\mathop{\mathrm{Log}}X$は

```math
\begin{align*}
\mathop{\mathrm{Log}}X = V \left(\begin{array}{cccc}
\log(\sigma_1) & & &\\
& \log(\sigma_2) & &\\
& & \ddots &\\
& & & \log(\sigma_m)
\end{array}\right)V^\top
\end{align*}
```

で計算できる．

SciPy による実装は非対称行列にも対応したアルゴリズムなため，正定値対称行列に対しては上記の手法を用いたほうが速く計算できる．このアルゴリズムも行列指数関数と同様に多様体上の最適化ライブラリ Pymanopt に実装されている．細かい実装は [Github](https://github.com/pymanopt/pymanopt/blob/1.0.0/pymanopt/tools/multi.py) を参照．

ここでは，SciPy と Pymanopt の計算時間を比較する．

```python
import numpy as np
from scipy.linalg import logm
from pymanopt.tools.multi import multilog


np.random.seed(0)
x = np.random.randn(10, 10)
x = x @ x.T + np.identity(10) * 1e-8

%time logmx = logm(x)
# CPU times: user 8.49 ms, sys: 3.44 ms, total: 11.9 ms
# Wall time: 8.98 ms
%time logx = multilog(x, pos_def=True)
# CPU times: user 550 µs, sys: 713 µs, total: 1.26 ms
# Wall time: 784 µs

np.testing.assert_array_almost_equal(logmx, logx)
```

# まとめ
対称行列に対しては Pymanopt (`pymanopt.tools.multi`) を，非対称行列には SciPy (`scipy.linalg`) を使いましょう．

* [Pymanopt](https://www.pymanopt.org/)
* [SciPy](https://docs.scipy.org/doc/scipy/reference/)

## 参考文献

* [C. B. Moler and C. F. Van Loan, Nineteen dubious ways to compute the exponential of a matrix, SIAM Rev., 20 (1978), pp. 801–836.](https://epubs.siam.org/doi/pdf/10.1137/1020098?casa_token=bnCWIRyIsUMAAAAA:2rI8Vtrk70dosbw9Xt2Jp4ra_QNmahKSKW6PRRPLSIYDKJjrdMZOaEtkI3WA6yP0pfHPYr8pEhh_rA)
* [Awad H. Al-Mohy and Nicholas J. Higham (2012) “Improved Inverse Scaling and Squaring Algorithms for the Matrix Logarithm.” SIAM Journal on Scientific Computing, 34 (4). C152-C169. ISSN 1095-7197](http://eprints.maths.manchester.ac.uk/1687/1/paper11.pdf)