行列指数関数の近似計算精度

定義

行列指数関数は次式で定義される．行列指数関数の性質はwikipediaを参照されたい．

e^A = \sum_{k=0}^\infty \frac{1}{k!}A^k

定義は上記の通りであり簡単に求まりそうであるが，これを厳密に計算しようとするとかなり煩雑である（計算過程で虚数を扱ったりする）．
$e^A$の要素を初等関数で表したかったり，あるいは組込み制御系で計算機が貧弱だったりする場合，精度はほどほどで良いから簡単に扱いたいという場合がある．

途中で打ち切った場合

行列指数関数の最も簡単な近似方法は，定義の式の通りに計算し，計算を途中で打ち切る方法だと考えられる．
これは，マクローリン展開，及び分母多項式を$1$としたときののPade近似と一致する．

この近似方法で問題になるのは近似精度である．
一般論を考えるのは難しいので，実際に値を突っ込んで計算してみる．

一例として，

A = \begin{bmatrix}
0 & 1 \\
-1 & 0 \\
\end{bmatrix}

の場合を考える．$e^A$の厳密な値は次の通り．

e^A = \begin{bmatrix}
  0.540302 & 0.841471\\
 -0.841471 & 0.540302 \\
\end{bmatrix}

途中で打ち切ったものを近似式とし，$f(A,n) := \sum_{k=0}^n \frac{1}{k!}A^k$と表すことにする．
$k=1$~$4,10$のときの近似値は次のようになる．

f(A, 1) = \begin{bmatrix}
  1.0 & 1.0 \\
 -1.0 & 1.0 \\
\end{bmatrix} \\
f(A, 2) = \begin{bmatrix}
  0.5 & 1.0 \\
 -1.0 & 0.5 \\
\end{bmatrix} \\
f(A, 3) = \begin{bmatrix}
  0.5 & 0.833333 \\
 -0.833333 & 0.5 \\
\end{bmatrix} \\
f(A, 4) = \begin{bmatrix}
  0.541667 & 0.833333 \\
 -0.833333 & 0.541667 \\
\end{bmatrix} \\
f(A, 10) = \begin{bmatrix}
  0.540302 & 0.841471 \\
 -0.841471 & 0.540302 \\
\end{bmatrix}

真値と比較すると，上記の$A$では1,2次ぐらいでは十分近似できているとは認めがたい．

近似誤差をフロベニウスノルム$\sqrt{\rm{trace}(AA^{\rm{T}})}=\sqrt{\sum_{i, j} |A_{i,j}|^2}$によって評価する．
横軸に$k$を取り，縦軸を誤差$e^A-f(A,k)$のフロベニウスノルムとしたグラフが次である．2つ目のグラフは縦軸を対数表示している．

今回の$A$では，少なくとも4次くらいまでは計算しておきたい．
これは，行列$A$に大きく依存するので，どこまでで近似するかは必要に応じて適宜求める必要がある．

Aの累乗から求められている点から，Aの各要素の値が十分小さければ，低次で打ち切ってもよいと考えられる．

# Julia 1.6.2

using Plots
using LinearAlgebra

# exp(A)のk次近似関数 (AはSymPy.Symの行列でもOK）
exp_Maclaurin(A::Matrix{T}, k::Int) where T = sum((A^i/factorial(i) for i=0:k))

A = [0. 1.;
     -1. 0.]

p = plot(0:15, [fnorm(exp(A)-exp_approx(A,i)) for i = 0:15], marker=:circle, ms=5, label="")
#p = plot(0:15, [fnorm(exp(A)-exp_approx(A,i)) for i = 0:15], marker=:circle, ms=5, label="", yscale=:log10)
xlabel!("k")
ylabel!("norm(exp(a) - f(a,k))")

Pade近似

スカラーの行列指数関数と同様に，Pade近似を使うこともできる．
行列指数関数の定義式の通り，行列指数関数$e^A$はスカラーの指数関数$e^a$のマクローリン展開と見た目が一致している．
よって，Pade近似の係数の導出も同様の議論ができ，スカラーの場合の係数をそのまま流用できる．

近似式は次のようになる．

\begin {align}
e^X \approx & g(X,m,n) \\
:= & [b_0+b_1X+b_2X^2+\cdots+b_nX^n]^{-1} \\
& [a_0+a_1X+a_2X^2+\cdots+a_mX^m]
\end {align}

$a_i$, $b_i$は次を満たす指数関数のPade近似の係数である．

e^x \approx \frac{a_0+a_1x+a_2x^2+\cdots+a_mx^m}
{b_0+b_1x+b_2x^2+\cdots+b_nx^n}

Pade近似の係数は，手計算するか，ググるか，JuliaのPolynomialsやPythonのmpmathにあるpade関数を使えばよい．

using Polynomials

m = 2   # 分子の次数
n = 3   # 分母の次数

# Pade近似の多項式(分子, 分母)
function exp_pade(A::Matrix{T}, m::Int, n::Int) where T
    # e^Aのマクローリン展開多項式
    buf = Polynomial([1/factorial(k) for k = 0:m+n])

    # 近似関数
    f = Polynomials.PolyCompat.Pade(buf, m, n)

    return f(A)
end

exp_pade(A, m, n)

一般にPade近似はマクローリン展開より良い近似になるようで，その精度はテイラー展開に対して$O(m+n+1)$である．原点から遠い点ではマクロリーン展開より精度は高いだろうが，計算に逆行列が必要になる．

古典制御論におけるむだ時間要素の伝達関数$e^{-Ls}$に対しては，近似した伝達関数をプロパーにしたい＝分子の次数を分母の次数以下にしたいのでPade近似が有用である．とにかく$e^A$の値や式が欲しいという場合には，逆行列があるので使い勝手が悪い．

e^A = \begin{bmatrix}
  0.540302 & 0.841471\\
 -0.841471 & 0.540302 \\
\end{bmatrix} \\

f(A, 5) = \begin{bmatrix}
  0.541667 & 0.841667 \\
 -0.841667 & 0.541667 \\
\end{bmatrix} \\
f(A, 6) = \begin{bmatrix}
  0.540278 & 0.841667 \\
 -0.841667 & 0.540278 \\
\end{bmatrix} \\

g(A, 2, 3) \begin{bmatrix}
  0.540251 & 0.841349 \\
  -0.841349 & 0.540251 \\
\end{bmatrix} \\
g(A, 3, 3) \begin{bmatrix}
  0.54031 & 0.841466 \\
  -0.841466 & 0.54031 \\
\end{bmatrix} \\

\qquad \\

\|f(A,5)-e^A\| = 1.94924 \times 10^{-3} \\
\|f(A,6)-e^A\| = 2.78901 \times 10^{-4} \\
\|g(A,2,3)-e^A\| = 1.87631 \times 10^{-4} \\
\|g(A,3,3)-e^A\| = 1.34915 \times 10^{-5} \\

SymPy(Julia)で行列指数関数を求める

SymPyを用いれば，簡単に求まる場合もある．
上で定義した$A$に非負の定数$t$をかけたものの行列指数関数$e^{At}$は次のように調べられる．
ちなみに，$e^{At}$は制御工学などで頻出し，LTIシステムの状態遷移行列となる．

# Julia 1.6.2

using SymPy

A = [0. 1.;
     -1. 0.]

# t >= 0
t = symbols("t", positive=true)

exp(A*t)

e^{At} = \begin{bmatrix}
  1.0⋅\cos(t) & 1.0⋅\sin(t) \\
 -1.0⋅\sin(t) & 1.0⋅\cos(t) \\
\end{bmatrix}

$t=1$を代入してみると，$e^A$と一致する．
（たまにおかしな解が出る？のであまり信用できていない）