行列指数関数 exp(A) の理論と応用

電子回路のRC直列回路では、電圧の時間変化は微分方程式 $\dot{v} = -\frac{1}{RC}v$ に従い、解は $v(t) = v_0 e^{-t/RC}$ です。では、複数の変数が連立する場合 — 例えば多自由度の振動系や状態空間モデル

$$ \dot{\bm{x}}(t) = \bm{A}\bm{x}(t) $$

の解はどうなるでしょうか。スカラーの場合との類推で $\bm{x}(t) = e^{\bm{A}t}\bm{x}(0)$ と書きたくなりますが、そもそも「行列の指数関数」とは何でしょうか。

行列指数関数 $e^{\bm{A}}$ は、スカラーの指数関数のテイラー展開を行列に拡張した概念です。$e^x = 1 + x + \frac{x^2}{2!} + \cdots$ の $x$ を行列 $\bm{A}$ に置き換えて

$$ e^{\bm{A}} = \bm{I} + \bm{A} + \frac{\bm{A}^2}{2!} + \frac{\bm{A}^3}{3!} + \cdots $$

この無限級数は任意の正方行列に対して収束し、連立微分方程式の解を与えます。

行列指数関数を理解すると、以下のような場面で活用できます。

  • 制御工学: 線形時不変系の状態遷移行列 $e^{\bm{A}t}$
  • 量子力学: 時間発展演算子 $e^{-i\bm{H}t/\hbar}$
  • マルコフ過程: 連続時間マルコフ連鎖の遷移確率 $e^{\bm{Q}t}$
  • ロボティクス: 回転行列のリー代数表現(ロドリゲスの公式)

本記事の内容

  • 行列指数関数の定義と収束性
  • 基本性質
  • 対角化可能な行列での計算
  • ジョルダン標準形での計算
  • 微分方程式との関係
  • 数値計算法(スケーリング&スクエアリング、Padé近似)
  • Pythonでの実装と可視化

前提知識

この記事を読む前に、以下の記事を読んでおくと理解が深まります。

行列指数関数の定義

テイラー展開による定義

$n \times n$ 正方行列 $\bm{A}$ に対して、行列指数関数(matrix exponential)は次のべき級数で定義されます。

$$ \begin{equation} e^{\bm{A}} = \exp(\bm{A}) = \sum_{k=0}^{\infty} \frac{\bm{A}^k}{k!} = \bm{I} + \bm{A} + \frac{\bm{A}^2}{2!} + \frac{\bm{A}^3}{3!} + \cdots \end{equation} $$

ここで $\bm{A}^0 = \bm{I}$ です。

収束性

この級数は任意の正方行列 $\bm{A}$ に対して絶対収束します。行列ノルム(例えばスペクトルノルム $\|\cdot\|_2$ やフロベニウスノルム $\|\cdot\|_F$)で評価すると

$$ \left\|\frac{\bm{A}^k}{k!}\right\| \leq \frac{\|\bm{A}\|^k}{k!} $$

右辺の和 $\sum_{k=0}^{\infty} \frac{\|\bm{A}\|^k}{k!} = e^{\|\bm{A}\|}$ は収束するので、行列級数も収束します。

収束の速さはスカラーの場合と同じです。$\|\bm{A}\|$ が小さければ少ない項数で良い近似が得られますが、$\|\bm{A}\|$ が大きいと多くの項が必要になります。この問題は後で議論する数値計算法で対処します。

ゼロ行列と単位行列

$$ e^{\bm{O}} = \bm{I} $$

スカラーの $e^0 = 1$ に対応します。

$$ e^{c\bm{I}} = e^c \bm{I} $$

対角成分が全て $c$ の行列なので、スカラーの指数関数が各成分に作用します。

行列指数関数の定義を押さえたところで、その性質を調べましょう。

基本性質

行列式との関係

行列指数関数の行列式は、トレースの指数関数で表されます。

$$ \begin{equation} \det(e^{\bm{A}}) = e^{\text{tr}(\bm{A})} \end{equation} $$

証明: 対角化可能な場合、$\bm{A}$ の固有値を $\lambda_1, \ldots, \lambda_n$ とすると $e^{\bm{A}}$ の固有値は $e^{\lambda_1}, \ldots, e^{\lambda_n}$ です。したがって

$$ \det(e^{\bm{A}}) = \prod_{i=1}^n e^{\lambda_i} = e^{\sum_i \lambda_i} = e^{\text{tr}(\bm{A})} $$

この公式はヤコビの公式(Jacobi’s formula)と呼ばれ、$e^{\bm{A}}$ は常に正則であることを示しています($e^{\text{tr}(\bm{A})} > 0$)。

逆行列

$e^{\bm{A}}$ の逆行列は

$$ (e^{\bm{A}})^{-1} = e^{-\bm{A}} $$

これはスカラーの $e^{-x} = (e^x)^{-1}$ に対応します。

転置と共役転置

$$ (e^{\bm{A}})^T = e^{\bm{A}^T}, \quad (e^{\bm{A}})^* = e^{\bm{A}^*} $$

交換可能な行列の場合

$\bm{AB} = \bm{BA}$ のとき($\bm{A}$ と $\bm{B}$ が可換のとき)

$$ \begin{equation} e^{\bm{A}+\bm{B}} = e^{\bm{A}} e^{\bm{B}} \end{equation} $$

注意: 一般に $\bm{AB} \neq \bm{BA}$ のとき $e^{\bm{A}+\bm{B}} \neq e^{\bm{A}}e^{\bm{B}}$ です。これはスカラーの場合(常に可換)との大きな違いであり、行列指数関数の取り扱いを複雑にする原因です。

非可換の場合の正しい公式はベイカー・キャンベル・ハウスドルフの公式(BCH formula)で与えられます。

$$ e^{\bm{A}}e^{\bm{B}} = e^{\bm{A}+\bm{B}+\frac{1}{2}[\bm{A},\bm{B}]+\frac{1}{12}([\bm{A},[\bm{A},\bm{B}]]+[\bm{B},[\bm{B},\bm{A}]])+\cdots} $$

ここで $[\bm{A},\bm{B}] = \bm{AB} – \bm{BA}$ は交換子です。

微分

$$ \frac{d}{dt} e^{\bm{A}t} = \bm{A} e^{\bm{A}t} = e^{\bm{A}t} \bm{A} $$

$\bm{A}$ と $e^{\bm{A}t}$ は常に可換です($e^{\bm{A}t}$ は $\bm{A}$ のべき級数なので)。この性質が連立微分方程式の解に直結します。

基本性質を整理したところで、実際にどう計算するかを見ていきましょう。

対角化による計算

対角化可能な場合

$\bm{A} = \bm{PDP}^{-1}$ と対角化できるとき

$$ \bm{A}^k = \bm{PD}^k\bm{P}^{-1} $$

なので

$$ e^{\bm{A}} = \sum_{k=0}^{\infty} \frac{\bm{A}^k}{k!} = \bm{P} \left(\sum_{k=0}^{\infty} \frac{\bm{D}^k}{k!}\right) \bm{P}^{-1} = \bm{P} e^{\bm{D}} \bm{P}^{-1} $$

対角行列の指数関数は各対角成分に独立に作用します。

$$ \begin{equation} e^{\bm{D}} = \begin{pmatrix} e^{\lambda_1} & & \\ & \ddots & \\ & & e^{\lambda_n} \end{pmatrix} \end{equation} $$

したがって

$$ e^{\bm{A}} = \bm{P} \begin{pmatrix} e^{\lambda_1} & & \\ & \ddots & \\ & & e^{\lambda_n} \end{pmatrix} \bm{P}^{-1} $$

対称行列の場合

$\bm{A} = \bm{A}^T$ なら直交行列 $\bm{Q}$ で $\bm{A} = \bm{QDQ}^T$ と対角化でき

$$ e^{\bm{A}} = \bm{Q} e^{\bm{D}} \bm{Q}^T $$

$e^{\bm{A}}$ も対称行列になります。

計算例

$\bm{A} = \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix}$ の行列指数関数を計算します。

固有値: $\lambda = \pm i$(虚数固有値)

固有ベクトルを使って対角化し計算すると

$$ e^{\bm{A}t} = \begin{pmatrix} \cos t & \sin t \\ -\sin t & \cos t \end{pmatrix} $$

これは回転行列です。反対称行列 $\bm{A}$ の指数関数は直交行列(回転行列)を生成するという重要な性質が表れています。ロボティクスでは、この性質を利用して回転をリー代数(反対称行列)で表現します。

ジョルダン標準形による計算

対角化不可能な場合

$\bm{A}$ が対角化不可能な場合は、ジョルダン標準形 $\bm{A} = \bm{PJP}^{-1}$ を使います。

$$ e^{\bm{A}} = \bm{P} e^{\bm{J}} \bm{P}^{-1} $$

$\bm{J}$ がブロック対角なので $e^{\bm{J}}$ は各ジョルダンブロックの指数関数のブロック対角行列です。

ジョルダンブロックの指数関数

サイズ $k$ のジョルダンブロック $\bm{J}_k(\lambda) = \lambda\bm{I} + \bm{N}$ と分解できます。ここで $\bm{N}$ はべき零行列($\bm{N}^k = \bm{O}$)です。

$\lambda\bm{I}$ と $\bm{N}$ は可換なので

$$ e^{\bm{J}_k(\lambda)} = e^{\lambda\bm{I} + \bm{N}} = e^{\lambda\bm{I}} e^{\bm{N}} = e^{\lambda} e^{\bm{N}} $$

$\bm{N}$ はべき零なのでテイラー展開が有限項で打ち切れます。

$$ e^{\bm{N}} = \bm{I} + \bm{N} + \frac{\bm{N}^2}{2!} + \cdots + \frac{\bm{N}^{k-1}}{(k-1)!} $$

具体的に $k = 3$ の場合

$$ e^{\bm{J}_3(\lambda)t} = e^{\lambda t} \begin{pmatrix} 1 & t & \frac{t^2}{2} \\ 0 & 1 & t \\ 0 & 0 & 1 \end{pmatrix} $$

対角化可能な場合は指数関数のみ $e^{\lambda t}$ が現れますが、ジョルダンブロックの場合は $t$ の多項式($1, t, t^2/2, \ldots$)が追加で現れます。これは微分方程式の解で「共鳴」を起こす項に対応します。

理論的な計算法を見てきましたが、数値計算ではこれらの方法は必ずしも最良ではありません。より実用的な数値計算法を見ましょう。

数値計算法

テイラー展開の直接計算の問題点

定義通りに $e^{\bm{A}} \approx \sum_{k=0}^{N} \frac{\bm{A}^k}{k!}$ を計算する方法は、$\|\bm{A}\|$ が大きいと $N$ を非常に大きく取る必要があり、丸め誤差の蓄積も問題になります。

スケーリング&スクエアリング法

最も広く使われる数値計算法はスケーリング&スクエアリング(scaling and squaring)です。

アイデア: $e^{\bm{A}} = (e^{\bm{A}/2^s})^{2^s}$ を利用します。

  1. スケーリング: $s$ を選んで $\|\bm{A}/2^s\|$ が十分小さくなるようにする
  2. 近似: $e^{\bm{A}/2^s}$ をPadé近似で計算する
  3. スクエアリング: 結果を $s$ 回二乗する

$\|\bm{A}/2^s\|$ が小さければPadé近似が高精度で使え、二乗操作は $O(sn^3)$ で済みます。

Padé近似

Padé近似は有理関数近似であり、テイラー展開より広い範囲で良い近似を与えます。次数 $(p, q)$ のPadé近似は

$$ e^{\bm{X}} \approx [\bm{D}_{pq}(\bm{X})]^{-1} \bm{N}_{pq}(\bm{X}) $$

ここで $\bm{N}_{pq}$ と $\bm{D}_{pq}$ は $\bm{X}$ の多項式です。実用的には $p = q$(対角Padé近似)が使われ、SciPyでは $(13, 13)$ 次のPadé近似を採用しています。

なお、SciPyの expm 関数は内部でこのスケーリング&スクエアリング法を用いており、ユーザーが手動でノルムの大きさを心配する必要はありません。しかし、行列のサイズが非常に大きい場合(数千次元以上)は計算コストが $O(n^3)$ になるため、クリロフ部分空間法などの反復的手法が検討されることもあります。行列が疎(スパース)構造を持つ場合は特に反復法が有利です。

数値計算法の理論を踏まえて、Pythonで実装・可視化しましょう。

Pythonでの実装と可視化

行列指数関数の基本計算

import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import expm

# テイラー展開による行列指数関数
def matrix_exp_taylor(A, N=50):
    """テイラー展開でexp(A)を近似"""
    n = A.shape[0]
    result = np.eye(n)
    term = np.eye(n)
    for k in range(1, N + 1):
        term = term @ A / k
        result += term
    return result

# テスト
A1 = np.array([[1, 2], [0, 3]], dtype=float)
A2 = np.array([[0, 1], [-1, 0]], dtype=float)  # 反対称行列

exp_A1_taylor = matrix_exp_taylor(A1)
exp_A1_scipy = expm(A1)

print("=== 行列指数関数の計算 ===\n")
print(f"exp(A1) [Taylor N=50]:\n{exp_A1_taylor}\n")
print(f"exp(A1) [SciPy]:\n{exp_A1_scipy}\n")
print(f"差: {np.max(np.abs(exp_A1_taylor - exp_A1_scipy)):.2e}\n")

# 性質の検証
print("=== 性質の検証 ===\n")
print(f"det(exp(A1)) = {np.linalg.det(exp_A1_scipy):.6f}")
print(f"exp(tr(A1))  = {np.exp(np.trace(A1)):.6f}")
print(f"差: {abs(np.linalg.det(exp_A1_scipy) - np.exp(np.trace(A1))):.2e}\n")

# 反対称行列 → 直交行列
exp_A2 = expm(A2)
print(f"A2 (antisymmetric):\n{A2}\n")
print(f"exp(A2):\n{exp_A2}\n")
print(f"exp(A2)^T @ exp(A2) =\n{exp_A2.T @ exp_A2}")
print(f"直交性誤差: {np.max(np.abs(exp_A2.T @ exp_A2 - np.eye(2))):.2e}")

テイラー展開とSciPyの計算が一致し、$\det(e^{\bm{A}}) = e^{\text{tr}(\bm{A})}$ が数値的にも確認されます。反対称行列の指数関数が直交行列になることも確認できます。

連立微分方程式の可視化

fig, axes = plt.subplots(2, 3, figsize=(16, 10))

# システム1: 安定ノード(実固有値、負)
A_node = np.array([[-2, 0], [0, -1]])

# システム2: 安定渦(複素固有値、負の実部)
A_spiral = np.array([[-0.5, 2], [-2, -0.5]])

# システム3: 鞍点(実固有値、正と負)
A_saddle = np.array([[1, 0], [0, -2]])

systems = [
    (A_node, 'Stable Node\n$\\lambda = -2, -1$'),
    (A_spiral, 'Stable Spiral\n$\\lambda = -0.5 \\pm 2i$'),
    (A_saddle, 'Saddle Point\n$\\lambda = 1, -2$')
]

t = np.linspace(0, 5, 300)

for col, (A_sys, title) in enumerate(systems):
    eigenvalues_sys = np.linalg.eigvals(A_sys)
    print(f"{title}: eigenvalues = {eigenvalues_sys}")

    # 相平面
    ax = axes[0, col]
    x0_list = [
        [1, 0], [-1, 0], [0, 1], [0, -1],
        [1, 1], [-1, -1], [1, -1], [-1, 1]
    ]

    for x0 in x0_list:
        x0 = np.array(x0, dtype=float)
        trajectory = np.array([expm(A_sys * ti) @ x0 for ti in t])
        ax.plot(trajectory[:, 0], trajectory[:, 1], linewidth=1.2)
        ax.plot(x0[0], x0[1], 'ko', markersize=3)

    # ベクトル場
    xx, yy = np.meshgrid(np.linspace(-2, 2, 10), np.linspace(-2, 2, 10))
    dx = A_sys[0, 0] * xx + A_sys[0, 1] * yy
    dy = A_sys[1, 0] * xx + A_sys[1, 1] * yy
    ax.quiver(xx, yy, dx, dy, alpha=0.2, color='gray')

    ax.set_xlim(-2, 2)
    ax.set_ylim(-2, 2)
    ax.set_aspect('equal')
    ax.set_title(title, fontsize=11)
    ax.set_xlabel('$x_1$', fontsize=11)
    ax.set_ylabel('$x_2$', fontsize=11)
    ax.grid(True, alpha=0.3)

    # 時間発展
    ax = axes[1, col]
    x0 = np.array([1.0, 0.5])
    trajectory = np.array([expm(A_sys * ti) @ x0 for ti in t])
    ax.plot(t, trajectory[:, 0], 'b-', linewidth=2, label='$x_1(t)$')
    ax.plot(t, trajectory[:, 1], 'r-', linewidth=2, label='$x_2(t)$')

    # 固有値による減衰/成長の参考線
    for lam in eigenvalues_sys:
        if np.isreal(lam):
            lam_real = lam.real
            ax.plot(t, np.exp(lam_real * t), 'k--', alpha=0.3, linewidth=1)

    ax.set_xlabel('$t$', fontsize=11)
    ax.set_ylabel('$x(t)$', fontsize=11)
    ax.legend(fontsize=9)
    ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('matrix_exponential_ode.png', dpi=150, bbox_inches='tight')
plt.show()

このグラフから、行列指数関数と微分方程式の解の関係が明確にわかります。

  1. 左列(安定ノード): 固有値が共に負の実数($-2, -1$)なので、全ての軌道が原点に直線的に収束します。上図の相平面で軌道が直線的に原点に向かい、下図で両成分が単調に減衰しています。速い方向($\lambda = -2$)で先に収束し、遅い方向($\lambda = -1$)が支配的になります

  2. 中央列(安定渦): 固有値が複素数($-0.5 \pm 2i$)で実部が負なので、螺旋状に収束します。虚部が振動周期を、実部が減衰率を決めています。下図で振動しながら減衰する様子が確認でき、$e^{-0.5t}$ のエンベロープ内に収まっています

  3. 右列(鞍点): 固有値が正($+1$)と負($-2$)の混合なので、安定方向と不安定方向が共存します。安定多様体に沿って原点に近づいた後、不安定多様体に沿って発散するパターンが上図で見えます。下図では $x_1$ が成長し $x_2$ が減衰しています

テイラー展開の収束の可視化

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# (a) テイラー展開の項数と精度
ax = axes[0]
A_test = np.array([[5, 1], [0, -3]], dtype=float)
exp_A_exact = expm(A_test)

N_values = range(1, 40)
errors = []
for N in N_values:
    exp_approx = matrix_exp_taylor(A_test, N)
    errors.append(np.max(np.abs(exp_approx - exp_A_exact)))

ax.semilogy(list(N_values), errors, 'bo-', markersize=4, linewidth=1.5)
ax.axhline(y=1e-15, color='red', linewidth=1, linestyle='--', label='Machine epsilon')
ax.set_xlabel('Number of Taylor terms $N$', fontsize=12)
ax.set_ylabel('Max error', fontsize=12)
ax.set_title(f'Taylor Series Convergence ($\\|A\\|_2 = {np.linalg.norm(A_test, 2):.1f}$)', fontsize=13)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3, which='both')

# (b) 非可換性: exp(A+B) vs exp(A)exp(B)
ax = axes[1]
A_nc = np.array([[1, 1], [0, 0]], dtype=float)
B_nc = np.array([[0, 0], [1, 1]], dtype=float)

# 可換かチェック
commutator = A_nc @ B_nc - B_nc @ A_nc
print(f"[A, B] = AB - BA =\n{commutator}")
print(f"可換? {np.allclose(commutator, 0)}")

t_vals = np.linspace(0, 2, 50)
diff_norm = []
for t_val in t_vals:
    exp_sum = expm((A_nc + B_nc) * t_val)
    exp_prod = expm(A_nc * t_val) @ expm(B_nc * t_val)
    diff_norm.append(np.linalg.norm(exp_sum - exp_prod, 'fro'))

ax.plot(t_vals, diff_norm, 'r-', linewidth=2)
ax.set_xlabel('$t$', fontsize=12)
ax.set_ylabel('$\\|e^{(A+B)t} - e^{At}e^{Bt}\\|_F$', fontsize=12)
ax.set_title('Non-commutativity: $e^{A+B} \\neq e^A e^B$', fontsize=13)
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('matrix_exp_convergence.png', dpi=150, bbox_inches='tight')
plt.show()

このグラフから、行列指数関数の計算特性が確認できます。

  1. 左図(テイラー展開の収束): 項数 $N$ を増やすにつれて誤差が急速に減少し、$N \approx 25$ で機械精度に達しています。$\|\bm{A}\|_2 \approx 5$ のこの行列では25項で十分ですが、ノルムが大きい行列ではさらに多くの項が必要になり、スケーリング&スクエアリング法の必要性が理解できます

  2. 右図(非可換性): $e^{(\bm{A}+\bm{B})t}$ と $e^{\bm{A}t}e^{\bm{B}t}$ の差のフロベニウスノルムが $t$ の増加とともに大きくなっています。$t = 0$ では両者は $\bm{I}$ で一致しますが、$t$ が大きくなると非可換性の影響が顕著になります。$[\bm{A}, \bm{B}] \neq \bm{O}$ であることがこの乖離の原因です

行列指数関数と回転

反対称行列 $\bm{\Omega}$ の指数関数は直交行列(回転行列)を生成するという重要な性質があります。

3次元の場合、角速度ベクトル $\bm{\omega} = (\omega_1, \omega_2, \omega_3)^T$ に対応する反対称行列

$$ [\bm{\omega}]_\times = \begin{pmatrix} 0 & -\omega_3 & \omega_2 \\ \omega_3 & 0 & -\omega_1 \\ -\omega_2 & \omega_1 & 0 \end{pmatrix} $$

の指数関数 $\bm{R} = e^{[\bm{\omega}]_\times \theta}$ は、$\bm{\omega}$ 軸周りに角度 $\theta = \|\bm{\omega}\|$ だけ回転する回転行列です。この関係はロドリゲスの回転公式として知られています。

$$ e^{[\bm{\omega}]_\times \theta} = \bm{I} + \sin\theta \, [\hat{\bm{\omega}}]_\times + (1 – \cos\theta) \, [\hat{\bm{\omega}}]_\times^2 $$

ここで $\hat{\bm{\omega}} = \bm{\omega}/\|\bm{\omega}\|$ は単位ベクトルです。

この公式はテイラー展開を直接計算するより効率的であり、ロボティクスやコンピュータグラフィックスで広く使われています。

行列指数関数と行列対数関数

行列指数関数の逆操作として行列対数関数 $\log(\bm{A})$ が定義されます。$e^{\bm{X}} = \bm{A}$ を満たす $\bm{X}$ が行列対数です。ただし、行列対数はスカラーの場合と同様に多価関数であり、主値を選ぶ必要があります。

$\bm{A}$ が正則で負の実数の固有値を持たないとき、主対数(principal logarithm)$\text{Log}(\bm{A})$ が一意に定まります。対角化可能な場合は

$$ \text{Log}(\bm{A}) = \bm{P} \begin{pmatrix} \ln\lambda_1 & & \\ & \ddots & \\ & & \ln\lambda_n \end{pmatrix} \bm{P}^{-1} $$

で計算できます。

行列対数関数は、回転行列 $\bm{R}$ から対応する反対称行列(角速度行列)を逆算する場面で重要です。たとえば、ロボティクスでは2つの姿勢 $\bm{R}_1$ と $\bm{R}_2$ の間の「距離」を $\|\text{Log}(\bm{R}_1^T\bm{R}_2)\|_F$ で定義します。これは測地線距離と呼ばれ、回転の大きさを適切に測る指標となります。

また、行列指数関数を繰り返し計算する場面では、マグナス展開(Magnus expansion)が有用です。時変係数の微分方程式 $\dot{\bm{x}}(t) = \bm{A}(t)\bm{x}(t)$ の解は一般には $e^{\int_0^t \bm{A}(\tau)d\tau}$ とはならず($\bm{A}$ が異なる時刻で非可換なため)、マグナス展開による補正項が必要です。

まとめ

本記事では、行列指数関数の理論と計算法について解説しました。

  • 行列指数関数 $e^{\bm{A}} = \sum_{k=0}^{\infty} \frac{\bm{A}^k}{k!}$ は任意の正方行列に対して収束し、連立微分方程式 $\dot{\bm{x}} = \bm{Ax}$ の解 $\bm{x}(t) = e^{\bm{A}t}\bm{x}_0$ を与える
  • $\det(e^{\bm{A}}) = e^{\text{tr}(\bm{A})}$ であり、$e^{\bm{A}}$ は常に正則
  • 対角化可能な場合は $e^{\bm{A}} = \bm{P} e^{\bm{D}} \bm{P}^{-1}$、ジョルダン標準形の場合はべき零行列の有限テイラー展開で計算
  • 非可換性に注意: $e^{\bm{A}+\bm{B}} = e^{\bm{A}}e^{\bm{B}}$ は $\bm{AB} = \bm{BA}$ の場合のみ成立
  • 数値計算ではスケーリング&スクエアリング法Padé近似の組み合わせが標準
  • 反対称行列の指数関数は回転行列を生成し、ロボティクスで重要

次のステップとして、以下の記事も参考にしてください。