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

行列指数関数

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

\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$ と同様に，

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

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

SciPy による実装

scipy.linalg.expm に行列指数関数の実装がある．
例として，次の行列 $A$ の行列指数関数 $\mathop{\mathrm{Exp}}A$ を考える．

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

$A^n$ が

\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$ は

\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 の計算結果を以下に示す．

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)$ を用いて

A = U\Lambda U^\top

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

\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$は

\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 を参照．

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

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 に行列対数関数の実装がある．
例として，行列指数関数の例で用いた行列 $A$ を用いる ($\lambda=1, \alpha=2, \mu=3$)．

\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*}

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)$ を用いて

X = V\Sigma V^\top

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

\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$は

\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 を参照．

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

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) を使いましょう．

参考文献