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

Rで 行列値関数:行列の対数 を試みます。

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

行列の対数は、行列の指数関数の逆演算として定義されます。

スカラーの世界で、\(e^x = a\) を満たす \(x\) が対数 \(x = \log a\) であるのと同様に、行列の世界でも、\(n \times n\) の正方行列 \(A\) が与えられたとき、 \[
e^X = A
\] を満たす行列 \(X\) を見つける問題が考えられます。この行列 \(X\) を \(A\) の行列の対数と呼び、\(\log A\) と書きます。

行列の対数は、以下のような分野で重要な役割を果たします。

リー群とリー代数: リー群（行列のなす連続的な群）とその接空間であるリー代数を対応付ける写像として、行列の指数関数と対数関数は中心的な役割を担います。これはロボットのアーム制御や3Dグラフィックスの回転計算などに応用されます。
制御理論: 連続時間システムと離散時間システムの間の変換に利用されます。
統計学・機械学習: 共分散行列の空間（リーマン多様体）上での補間や平均を計算する際に、対数写像が使われます。
医用画像処理: MRIのテンソル画像解析（DTI）などで、テンソルの平均や統計解析を行うために行列対数が用いられます。

2. 行列の対数の定義

\(n \times n\) の正方行列 \(A\) に対して、 \[
e^X = A
\] を満たす \(n \times n\) の正方行列 \(X\) を、\(A\) の対数行列と呼びます。

スカラーの対数と同様に、行列の対数にも注意すべき点があります。

存在: 対数行列は、正則な（逆行列を持つ）行列 \(A\) に対してのみ存在します。
一意性: スカラーの複素対数が \(2\pi i\) の整数倍の不定性を持つのと同様に、行列の対数も一意には定まりません。複数存在する場合があります。
主対数: 複数ある対数のうち、その固有値の虚部が \((-\pi, \pi]\) の範囲にあるものを主対数 (principal logarithm) と呼び、通常 \(\text{Log } A\) と書かれます。

3. 行列の対数の主な性質

\(X = \log A\) とします。

\(e^{\log A} = A\) （定義そのものです）
\(A\) が実行列であっても、\(\log A\) は複素行列になることがあります。これは \(A\) が負の固有値を持つ場合に起こります（例：\(\log(-1) = i\pi\)）。
\(\det(A) = e^{\text{tr}(\log A)}\) という関係が成り立ちます。これは、\(\det(e^X) = e^{\text{tr}(X)}\) から導かれます。

注意が必要な性質

スカラーの場合と異なり、一般的に以下の式は成立しません。 \[
\log(AB) \neq \log A + \log B
\] これが成り立つのは、\(A\) と \(B\) が可換（\(AB=BA\)）であるなど、非常に限定的な条件下のみです。

4. 行列の対数の計算方法

行列の対数を求めるには、いくつかの方法があります。

方法1：行列が対角化可能な場合

行列 \(A\) が対角化可能であるとは、ある正則行列 \(P\) を用いて、 \[
A = PDP^{-1}
\] と書けることを意味します。ここで \(D\) は \(A\) の固有値 \(\lambda_i\) を対角成分に持つ対角行列です。

このとき、対数行列 \(X\) は \(X = P (\log D) P^{-1}\) として計算できます。ここで \(\log D\) は、\(D\) の各対角成分の対数をとった行列です。 \[
\log D =
\begin{pmatrix}
\log \lambda_1 & & 0 \\
& \ddots & \\
0 & & \log \lambda_n
\end{pmatrix}
\] \(A\) が正則であるという条件から、どの固有値 \(\lambda_i\) も \(0\) ではありません。もし負の固有値 \(\lambda_k < 0\) があれば、\(\log \lambda_k = \log(|\lambda_k|) + i\pi\) のように結果は複素数になります。

【手順】

行列 \(A\) の固有値 \(\lambda_i\) と、対応する固有ベクトルを求める。
固有ベクトルを並べて行列 \(P\) を作る。
固有値を対角成分に持つ対角行列 \(D\) を作る。
\(P\) の逆行列 \(P^{-1}\) を求める。
\(D\) の各対角成分の対数を計算して \(\log D\) を作る。
\(X = P (\log D) P^{-1}\) を計算する。

方法2：級数展開を利用する方法

スカラーのテイラー展開 \(\log(1+x) = \sum_{k=1}^{\infty} (-1)^{k-1} \dfrac{x^k}{k}\) を利用する方法もあります。行列 \(A\) のすべての固有値の絶対値が \(1\) に近い場合、\(I-A\) のノルムが \(1\) より小さいときなどに、以下の級数が収束します。 \[
\log A = -\sum_{k=1}^{\infty} \dfrac{(I-A)^k}{k}
\] この方法は特定の条件下での数値計算に用いられます。

5. 具体的な計算例

行列 \(A = \begin{pmatrix} 2 & 1 \\ 1 & 2 \end{pmatrix}\) の主対数を求めてみましょう。

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

固有多項式: \[\det(A - \lambda I) = (2-\lambda)^2 - 1 = \lambda^2 - 4\lambda + 3 = (\lambda-1)(\lambda-3) = 0\]
固有値: \(\lambda_1 = 1, \lambda_2 = 3\)
\(\lambda_1 = 1\) のとき、固有ベクトルは \(\mathbf{v}_1 = \begin{pmatrix} 1 \\ -1 \end{pmatrix}\)
\(\lambda_2 = 3\) のとき、固有ベクトルは \(\mathbf{v}_2 = \begin{pmatrix} 1 \\ 1 \end{pmatrix}\)

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

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

3. log D を計算する

\(\log D = \begin{pmatrix} \log 1 & 0 \\ 0 & \log 3 \end{pmatrix} = \begin{pmatrix} 0 & 0 \\ 0 & \log 3 \end{pmatrix}\)

4. X = P (log D) P⁻¹ を計算する \[
\begin{aligned}
X &= P (\log D) P^{-1} \\
&= \begin{pmatrix} 1 & 1 \\ -1 & 1 \end{pmatrix} \begin{pmatrix} 0 & 0 \\ 0 & \log 3 \end{pmatrix} \dfrac{1}{2} \begin{pmatrix} 1 & -1 \\ 1 & 1 \end{pmatrix} \\
&= \dfrac{1}{2} \begin{pmatrix} 0 & \log 3 \\ 0 & \log 3 \end{pmatrix} \begin{pmatrix} 1 & -1 \\ 1 & 1 \end{pmatrix} \\
&= \dfrac{1}{2} \begin{pmatrix} \log 3 & \log 3 \\ \log 3 & \log 3 \end{pmatrix} = \begin{pmatrix} \dfrac{\log 3}{2} & \dfrac{\log 3}{2} \\ \dfrac{\log 3}{2} & \dfrac{\log 3}{2} \end{pmatrix}
\end{aligned}
\] これが \(A\) の主対数です。

6. Rによる解法

準備：`expm` パッケージのロード

# パッケージをロード
library(expm)

`logm()` 関数の利用

logm() 関数に正方行列を渡すことで、その主対数が計算されます。

前のセクションで手計算した行列 \(A = \begin{pmatrix} 2 & 1 \\ 1 & 2 \end{pmatrix}\) について、Rで計算してみましょう。

# 行列 A を定義
A <- matrix(c(2, 1, 1, 2), nrow = 2)

cat("--- 行列 A を確認 ---\n")
print(A)

# logm() 関数を使って行列の主対数 X を計算
X <- logm(A)

cat("\n--- 行列 A の主対数 X を確認 ---\n")
print(X)

cat("\n--- log(3)/2 の値を確認 ---\n")
log(3) / 2

cat("\n--- 検算: expm(logm(A)) が A に戻るか確認 ---\n")
all.equal(expm(logm(A)), A)

--- 行列 A を確認 ---
     [,1] [,2]
[1,]    2    1
[2,]    1    2

--- 行列 A の主対数 X を確認 ---
          [,1]      [,2]
[1,] 0.5493061 0.5493061
[2,] 0.5493061 0.5493061

--- log(3)/2 の値を確認 ---
[1] 0.5493061

--- 検算: expm(logm(A)) が A に戻るか確認 ---
[1] TRUE