Rで行列値関数:行列の平方根

Rで 行列値関数:行列の平方根 を試みます。

1. はじめに：なぜ行列の平方根を考えるのか？

スカラーの世界では、ある数 \(a\) に対して \(x^2 = a\) を満たす数 \(x\) を \(a\) の平方根と呼びます。行列の平方根は、これを行列の世界に拡張したものです。

すなわち、\(n \times n\) の正方行列 \(A\) が与えられたとき、 \[
X^2 = A
\] を満たす行列 \(X\) を見つける問題です。この行列 \(X\) を \(A\) の行列の平方根と呼び、しばしば \(A^{1/2}\) と書かれます。

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

統計学: 多変量正規分布における共分散行列 \(\Sigma\) の分解\(\left(\Sigma = \Sigma^{1/2} (\Sigma^{1/2})^T\right)\)や、変数の標準化（マハラノビス距離の計算）などに使われます。
物理学: 量子力学における演算子の定義や、連続体力学における変形テンソルの計算などで現れます。
制御工学: Riccati方程式など、特定の行列方程式を解く際に利用されます。

2. 行列の平方根の定義

\(n \times n\) の正方行列 \(A\) に対して、

\[
X^2 = XX = A
\]

を満たす \(n \times n\) の正方行列 \(X\) を、\(A\) の平方根行列と呼びます。

スカラーの場合と異なり、行列の平方根には注意すべき点があります。

一意性: 平方根は一つとは限りません。複数存在する場合があります。
存在: 行列によっては、平方根が存在しない場合もあります。
主平方根: 複数ある平方根のうち、その固有値の実部がすべて非負であるものを主平方根 (principal square root) と呼び、通常 \(A^{1/2}\) はこれを指します。

特に、統計学などでよく扱われる半正定値対称行列は、一意な半正定値対称行列の平方根を持つことが知られており、これが応用上非常に重要です。

3. 行列の平方根の主な性質

\(X\) を \(A\) の平方根とします。

\((A^{1/2})^2 = A\) （定義そのものです）
\(A\) と \(A^{1/2}\) は可換です。つまり、\(A A^{1/2} = A^{1/2} A\) が成り立ちます。
もし \(A\) が対称行列ならば、その平方根の中にも対称なものが存在します。
もし \(A\) が半正定値対称行列ならば、一意な半正定値対称行列の平方根が存在します。

注意が必要な性質

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

4. 行列の平方根の計算方法

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

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

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

このとき、平方根行列 \(X\) を \(X = P D^{1/2} P^{-1}\) の形で探します。実際に2乗してみると、 \[\begin{eqnarray}X^2 &=& (P D^{1/2} P^{-1})^2 \\&=& (P D^{1/2} P^{-1})(P D^{1/2} P^{-1}) \\&=& P D^{1/2} (P^{-1}P) D^{1/2} P^{-1} \\&=& P (D^{1/2})^2 P^{-1} \\&=& P D P^{-1}\\ &=& A
\end{eqnarray}\] となり、\(X^2 = A\) が満たされることがわかります。

\(D\) の対角成分が \(\lambda_1, \lambda_2, \dots, \lambda_n\) であるとき、\(D^{1/2}\) は対角成分が \(\pm\sqrt{\lambda_1}, \pm\sqrt{\lambda_2}, \dots, \pm\sqrt{\lambda_n}\) となる対角行列です。

符号の組み合わせの数だけ平方根が存在することになります（\(2^n\) 個、ただし縮退がある場合は減る）。

主平方根は、すべての符号を正（複素数の場合は実部が正）に取ったものに対応します。

【手順】

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

方法2：数値解法（反復法）

ニュートン法などを応用した反復計算で平方根を求める方法もあります。例えば、ドゥンハム-ベイズ法として知られる以下の反復式は、適切な初期値 \(X_0\) から始めると \(A\) の平方根に収束します。 \[
X_{k+1} = \dfrac{1}{2}(X_k + A X_k^{-1})
\]

5. 具体的な計算例

行列 \(A = \begin{pmatrix} 5 & 4 \\ 4 & 5 \end{pmatrix}\) の主平方根を求めてみましょう。

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

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

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

固有ベクトルは正規化しておくと \(P\) が直交行列になり \(P^{-1} = P^T\) となって計算が簡便になりますが、ここでは一般的に解きます。

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

3. D¹/² を計算する

主平方根を求めるため、正の平方根をとります。

\[D^{1/2} = \begin{pmatrix} \sqrt{1} & 0 \\ 0 & \sqrt{9} \end{pmatrix} = \begin{pmatrix} 1 & 0 \\ 0 & 3 \end{pmatrix}\]

4. X = PD¹/²P⁻¹ を計算する \[
\begin{aligned}
X &= P D^{1/2} P^{-1} \\
&= \begin{pmatrix} 1 & 1 \\ -1 & 1 \end{pmatrix} \begin{pmatrix} 1 & 0 \\ 0 & 3 \end{pmatrix} \dfrac{1}{2} \begin{pmatrix} 1 & -1 \\ 1 & 1 \end{pmatrix} \\
&= \dfrac{1}{2} \begin{pmatrix} 1 & 3 \\ -1 & 3 \end{pmatrix} \begin{pmatrix} 1 & -1 \\ 1 & 1 \end{pmatrix} \\
&= \dfrac{1}{2} \begin{pmatrix} 4 & 2 \\ 2 & 4 \end{pmatrix} = \begin{pmatrix} 2 & 1 \\ 1 & 2 \end{pmatrix}
\end{aligned}
\] よって、主平方根は \(X = \begin{pmatrix} 2 & 1 \\ 1 & 2 \end{pmatrix}\) です。

検算

\[X^2 = \begin{pmatrix} 2 & 1 \\ 1 & 2 \end{pmatrix} \begin{pmatrix} 2 & 1 \\ 1 & 2 \end{pmatrix} = \begin{pmatrix} 4+1 & 2+2 \\ 2+2 & 1+4 \end{pmatrix} = \begin{pmatrix} 5 & 4 \\ 4 & 5 \end{pmatrix} = A\]

6. Rによる解法

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

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

`sqrtm()` 関数の利用

sqrtm() 関数に正方行列を渡すことで、その主平方根が計算されます。前のセクションで手計算した行列 \(A = \begin{pmatrix} 5 & 4 \\ 4 & 5 \end{pmatrix}\) について、Rで計算してみましょう。

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

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

# sqrtm() 関数を使って行列の主平方根 X を計算
X <- sqrtm(A)

# 計算結果を表示
cat("\n--- 行列 A の主平方根 X を確認 ---\n")
print(X)

# 検算: X^2 が A と等しいか確認
cat("\n--- 行列 A と求めた行列 A の主平方根 X の2乗は一致するか? ---\n")
all.equal(X %*% X, A)

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

--- 行列 A の主平方根 X を確認 ---
     [,1] [,2]
[1,]    2    1
[2,]    1    2

--- 行列 A と求めた行列 A の主平方根 X の2乗は一致するか? ---
[1] TRUE