Rで行列値関数:行列の指数関数

1. はじめに：なぜ行列の指数関数を考えるのか？

まず、私たちがよく知っているスカラー（ただの数）の指数関数 \(e^x\) を思い出してみましょう。

この関数は、微分方程式 \(\dfrac{dy}{dx} = ay\) の解が \(y = C e^{ax}\) となるように、微分方程式を解く上で重要な役割を果たします。

行列の指数関数は、これを行列とベクトルの世界に拡張したものです。 特に、以下のような連立線形微分方程式を解くために利用されます。

\[
\dfrac{d\mathbf{x}}{dt} = A\mathbf{x}
\]

ここで、

• \(\mathbf{x}\) は状態変数を並べたベクトル（例：位置と速度）
• \(A\) は係数を並べた正方行列
• \(\dfrac{d\mathbf{x}}{dt}\) は \(\mathbf{x}\) の時間微分

この方程式の解は、初期値を \(\mathbf{x}(0)\) とすると、スカラーの場合と同じ形で、

\[
\mathbf{x}(t) = e^{At} \mathbf{x}(0)
\]

と書くことができます。

この \(e^{At}\) が行列の指数関数です。物理学、制御工学、経済学など、さまざまな分野でシステムの動的な振る舞いを記述するために使われます。

2. 行列の指数関数の定義

\(n \times n\) の正方行列 \(A\) に対して、その指数関数 \(e^A\) は、スカラーの指数関数 \(e^x\) のテイラー展開を形式的に行列に置き換えることで定義されます。

スカラーの場合のテイラー展開は、 \[
e^x = 1 + x + \dfrac{x^2}{2!} + \dfrac{x^3}{3!} + \dots = \sum_{k=0}^{\infty} \dfrac{x^k}{k!}
\] です。これを行列 \(A\) で置き換えると、

\[
e^A = I + A + \dfrac{A^2}{2!} + \dfrac{A^3}{3!} + \dots = \sum_{k=0}^{\infty} \dfrac{A^k}{k!}
\]

となります。ここで、

• \(I\) は単位行列です (\(A^0 = I\) と考えます)。
• \(A^k\) は行列 \(A\) を \(k\) 回掛け合わせたものです。
• この無限級数は、任意の正方行列 \(A\) に対して必ず収束することが知られています。

3. 行列の指数関数の主な性質

スカラーの指数関数と似ている点、異なる点があります。

似ている性質

1. \(e^{\mathbf{0}} = I\) （\(\mathbf{0}\) は零行列）
2. \(e^A\) は常に正則（逆行列を持つ）。そして \((e^A)^{-1} = e^{-A}\)
3. \(\dfrac{d}{dt} e^{At} = A e^{At} = e^{At} A\) (これが微分方程式の解となる理由)
4. \(\det(e^A) = e^{\text{tr}(A)}\) （行列式の値は、Aのトレース(対角成分の和)の指数になる）

異なる性質

スカラーの場合は \(e^{x+y} = e^x e^y\) が常に成り立ちますが、行列では、\(e^{A+B} = e^A e^B\) が成り立つのは、\(A\) と \(B\) が可換（\(AB = BA\)）である場合に限られます。

4. 行列の指数関数の計算方法

5. 具体的な計算例

行列 \(A = \begin{pmatrix} 4 & 1 \\ 3 & 2 \end{pmatrix}\) の \(e^A\) を求めてみましょう。（対角化可能な例）

1. 固有値・固有ベクトルを求める

• 固有多項式: \(\det(A - \lambda I) = (4-\lambda)(2-\lambda) - 3 = \lambda^2 - 6\lambda + 5 = (\lambda-1)(\lambda-5) = 0\)

• 固有値: \(\lambda_1 = 1, \lambda_2 = 5\)

• \(\lambda_1 = 1\) のとき、固有ベクトルは \(\mathbf{v}_1 = \begin{pmatrix} 1 \\ -3 \end{pmatrix}\)

• \(\lambda_2 = 5\) のとき、固有ベクトルは \(\mathbf{v}_2 = \begin{pmatrix} 1 \\ 1 \end{pmatrix}\)

2. P, D, P⁻¹ を作る

• \(P = \begin{pmatrix} 1 & 1 \\ -3 & 1 \end{pmatrix}\)
• \(D = \begin{pmatrix} 1 & 0 \\ 0 & 5 \end{pmatrix}\)
• \(P^{-1} = \dfrac{1}{4} \begin{pmatrix} 1 & -1 \\ 3 & 1 \end{pmatrix}\)

3. e^D を計算する

• \(e^D = \begin{pmatrix} e^1 & 0 \\ 0 & e^5 \end{pmatrix} = \begin{pmatrix} e & 0 \\ 0 & e^5 \end{pmatrix}\)

4. e^A = Pe^DP⁻¹ を計算する

\[
\begin{aligned}
e^A &= P e^D P^{-1} \\
&= \begin{pmatrix} 1 & 1 \\ -3 & 1 \end{pmatrix} \begin{pmatrix} e & 0 \\ 0 & e^5 \end{pmatrix} \dfrac{1}{4} \begin{pmatrix} 1 & -1 \\ 3 & 1 \end{pmatrix} \\
&= \dfrac{1}{4} \begin{pmatrix} e & e^5 \\ -3e & e^5 \end{pmatrix} \begin{pmatrix} 1 & -1 \\ 3 & 1 \end{pmatrix} \\
&= \dfrac{1}{4} \begin{pmatrix} e + 3e^5 & -e + e^5 \\ -3e + 3e^5 & 3e + e^5 \end{pmatrix}
\end{aligned}
\] これが \(e^A\) の計算結果となります。

6. Rによる解法

ここでは、 expm パッケージを利用する方法を紹介します。

library(expm)

A <- matrix(c(4, 3, 1, 2), nrow = 2)

exp_A <- expm(A)
cat("--- expm {expm}により求めた e^A ---\n")
print(exp_A)

--- expm {expm}により求めた e^A ---
         [,1]     [,2]
[1,] 111.9894 36.42372
[2,] 109.2712 39.14200

e_val <- exp(1)
e5_val <- exp(5)
hand_calc_result <- (1 / 4) * matrix(
  c(
    e_val + 3 * e5_val, -e_val + e5_val,
    -3 * e_val + 3 * e5_val, 3 * e_val + e5_val
  ),
  nrow = 2, byrow = TRUE
)
cat("--- 手計算で求めた e^A ---\n")
hand_calc_result
cat("\n--- 両者は一致しているか? ---\n")
all.equal(exp_A, hand_calc_result)

--- 手計算で求めた e^A ---
         [,1]     [,2]
[1,] 111.9894 36.42372
[2,] 109.2712 39.14200

--- 両者は一致しているか? ---
[1] TRUE

まとめ

• 行列の指数関数 \(e^A\) は、スカラーの指数関数のテイラー展開を行列に拡張したものです。
• 最大の応用先は、連立線形微分方程式 \(\dfrac{d\mathbf{x}}{dt} = A\mathbf{x}\) を解くことであり、その解は \(\mathbf{x}(t) = e^{At} \mathbf{x}(0)\) となります。
• \(e^{A+B} = e^A e^B\) は、\(A\) と \(B\) が可換な場合にのみ成立します。

以上です。