最近書いたSVARモデルの話や、以前書いたVARモデルの直交化インパルス応答の話と少し関係あるのですが、コレスキー分解と言う行列分解があります。これをNumpyで行う方法を紹介しておこうと言う記事です。
Wikipediaでは正定値エルミート行列について定義が書いてありますが一旦、実行列をメインで考えるのでその直後に書いてある実対称行列の話で進めます。
正定値実対称行列$\mathbf{A}$を下三角行列$\mathbf{L}$とその転置行列$\mathbf{L}^T$の積に分解するのがコレスキー分解です。
NumpyやScipyには専用のメソッドが実装されおり、容易に行うことができます。NumpyでできるならNumpyでやってしまうのがオススメです。(scipyの方はなぜか上三角行列がデフォルトですし。)
参考: numpy.linalg.cholesky — NumPy v1.24 Manual
scipy.linalg.cholesky — SciPy v1.10.1 Manual
では適当な行列用意してやってみましょう。
import numpy as np
# 適当に正定値対称行列を用意する
A = np.array([
[2, 1, 3],
[1, 4, 2],
[3, 2, 6],
])
# 一応正定値なのを見ておく
print(np.linalg.eigvals(A))
# array([8.67447336, 0.39312789, 2.93239875])
# コレスキー分解とその結果を表示
L = np.linalg.cholesky(A)
print(L)
"""
[[1.41421356 0. 0. ]
[0.70710678 1.87082869 0. ]
[2.12132034 0.26726124 1.19522861]]
"""
# 転置行列と掛け合わせることで元の行列になることを確認
print(L@L.T)
"""
[[2. 1. 3.]
[1. 4. 2.]
[3. 2. 6.]]
"""簡単でしたね。
さて、ここまでだとメソッド一つ紹介しただけなのですが、実はここからが今回の記事の本題です。Wikipediaには修正コレスキー分解ってのが紹介されています。これは元の行列を対角成分が全て1の下三角行列$\mathbf{P}$と、対角行列$\mathbf{D}$と、1個目の三角行列の転置行列$\mathbf{P^T}$の3つの積に分解する方法です。(Wikipedia中の式ではこの三角行列もLになってますが、先出のLと被ると以降の説明が書きにくいのでここではPとします。)
$$\mathbf{A}=\mathbf{PDP^T}$$
と分解します。
直交化インパルス応答関数の記事中で使ってるのはこちらなのですが、どうやらこれをやってくれる関数が用意されてないのでやってみました。
いやいや、このメソッドがあるじゃないか、と言うのをご存知の人がいらっしゃいましたら教えてください。
さて、僕がやった方法ですがこれは単純で、コレスキー分解の結果をさらに分解し、
$$\mathbf{L}=\mathbf{PD^{1/2}}$$
とすることを目指します。
ちょっと計算するとすぐわかるのですが、下三角行列$\mathbf{L}$を対角成分が1の下三角行列と対角行列の積に分解しようとしたら、対角行列の方は元の三角行列の対角成分を取り出したものにするしかありません。あとはそれ使って$\mathbf{P}$も計算したら完成です。
対角成分を取り出した行列はnumpy.diagを使います。これ2次元配列に使うと対角成分を1次元配列で取り出し、1次元配列を渡すとそこから対角行列を生成すると言うなかなかトリッキーな動きをするので、対角成分を取り出して対角行列を作りたかったらこれを2回作用させます。平方根も外すために2乗することも忘れないようにすると、$\mathbf{D}$は次のようにして求まり、さらにその逆行列を使って$\mathbf{P}$も出せます。
# 対角行列Dの計算
D = np.diag(np.diag(L))**2
print(D)
"""
[[2. 0. 0. ]
[0. 3.5 0. ]
[0. 0. 1.42857143]]
"""
# Dを用いて、対角成分1の下三角行列Pを求める
P = L@np.linalg.inv(D**0.5)
print(P)
"""
[[1. 0. 0. ]
[0.5 1. 0. ]
[1.5 0.14285714 1. ]]
"""
# PDP^tで元の行列Aが復元できることを見る
print(P@D@P.T)
"""
[[2. 1. 3.]
[1. 4. 2.]
[3. 2. 6.]]
"""はい、これでできました。ブログの都合上、逐一printして出力も表示してるのでめんどくさそうな感じになってますが、実質以下の3行で実装できています。
L = np.linalg.cholesky(A)
D = np.diag(np.diag(L))**2
P = L@np.linalg.inv(D**0.5)