複素多変量正規分布からの効率的なサンプリングに関して解説します。Cholesky 分解に基づくサンプリングですので、多変量正規分布に従うサンプリング手法としては標準的なものです。
複素多変量正規分布
複素 N 次元 複素多変量正規分布とは、標本空間 $ \mathbb{C}^N $ 上の確率分布であり、平均（複素 N 次元）μ と分散共分散行列 （複素 N × N 正定値 Hermite 行列）Σ により指定されます。
複素多変量正規分布の確率密度関数
確率密度関数は、
p(z \mid \mu, \Sigma) = \frac{1}{\pi^N |\Sigma|} \, e^{ - (z - \mu)^* \Sigma^{-1} (z - \mu)}
で与えられます。
ここで、行列の右上肩の * は複素共役転置を表し、| | は行列の行列式を表しています。
複素多変量正規分布の回転対称性
複素多変量正規分布には回転対称性があります。例えば、平均 $0 \in \mathbb{C}^N$ の複素多変量正規分布に従う確率変数 Z を、N 次元複素平面の原点を中心に回転させて定義された確率変数
e^{\sqrt{-1} \phi} Z \quad (\phi \in \mathbb{R})
は、確率変数 Z と同じ確率分布に従います。
複素多変量正規分布の応用例
複素多変量正規分布には回転対称性があるので、「風力と風向」・「海面上の波の動き」・「地中の地震波」に代表される２次元の観測量に関して応用があります。また、電気信号は複素数で表記すると便利であることから、電気信号処理にも応用があります。
複素多変量正規分布と実多変量正規分布との関係
複素 N 次元 複素多変量正規分布は、実 2N 次元 実多変量正規分布の特別な場合です。実際、平均 μ・分散共分散行列 Σ の複素多変量正規分布に従う複素 N 次元 確率変数 Z に対して、実 2N 次元 確率変数 (Re(Z), Im(Z)) は、平均
\begin{bmatrix}
\Re{\mu} \\ \Im{\mu}
\end{bmatrix}
\in \mathbb{R}^{2N}
分散共分散行列（$2N \times 2N$ 実対称行列）
\begin{bmatrix}
\frac{1}{2} \Re{\Sigma} & -\frac{1}{2} \Im{\Sigma} \\
\frac{1}{2} \Im{\Sigma} & \frac{1}{2} \Re{\Sigma}
\end{bmatrix}
の実多変量正規分布に従います。
一方で、実 2N 次元 実多変量正規分布は、複素 N 次元 複素多変量正規分布とは限りません。つまり、複素 N 次元量 Z に対して、仮に確率変数 (Re(Z), Im(Z)) が実 2N 次元 実多変量正規分布に従っていたとしても、確率変数 Z が複素多変量正規分布に従う保証はありません。
複素多変量正規分布と疑似分散共分散行列
複素 N 次元確率変数 Z を考えます。簡単のため、Z の平均 μ は 0 とします。
もし、Z が複素多変量正規分布に従うのであれば、その分散共分散行列は
E [Z Z^*] = \Sigma
です。
一方で、Z' を Z の転置行列とするとき、E[Z Z'] は疑似分散共分散行列 (pseudo-variance-covariance matrix) と呼ばれ、Z が複素多変量正規分布に従うのであれば、
E[Z Z'] = 0
となります。このことは、例えば、実 2N 次元 確率変数 (Re(Z), Im(Z)) が従う分布に注意すれば分かります。このことからも、実 2N 次元 実多変量正規分布は、複素 N 次元 複素多変量正規分布とは限らないことが分かると思います。
つまり、「複素正規分布の回転対称性」と「複素正規分布の疑似分散共分散行列が零行列になる性質」との間には深い関係があります。
Cholesky 分解
Cholesky 分解とは、正定値エルミート行列 Σ を下三角行列 L と L の共役転置 L* との積に分解
\Sigma = L L^{*}
することを言います。
複素多変量正規分布からのサンプリング
複素 N 次元量 Z を、複素 N 次元量
W := L^{-1} (Z - \mu) \in \mathbb{C}^N
へと線形変換することを考えます。
確率変数 Z が平均 μ・分散共分散行列 Σ = L L* の複素多変量正規分布に従うのであれば、確率変数 W の確率密度関数は
p(w \mid \mu, \Sigma) = \frac{1}{\pi^N} \, e^{ - w^* w}
と書けることに注意します。つまり、w は標準複素多変量正規分布（分散共分散行列が単位行列）に従います。
したがって、複素多変量正規分布からのサンプリング z は、まず標準複素多変量正規分布から w をサンプリングして、次に z := Lw + μ とすれば良いことが分かります。実際、W を標準複素多変量正規分布とすると、(簡単のため μ = 0 として、)
\begin{align*}
E[(LW) (LW)^*] &= L ( E[W W^*] ) L^* = L L^* = \Sigma \\
E[(LW) (LW)'] &= L ( E[W W'] ) L' = 0
\end{align*}
となることが分かります。つまり、確率変数 Z := LW は分散共分散行列 Σ の複素多変量正規分布に従っています。
標準複素多変量正規分からのサンプリングは、実部と虚部とを独立に分散 1/2 の 1 次元の実正規分布から N 個 + N 個サンプリングすれば良いです。
注意したいのは、「標準複素多変量正規分からのサンプリング」を「標準実多変量正規分からのサンプリング」に置き換えることはできないということです。確かに、W を標準実多変量正規分に従う確率変数とすると、
E[(LW) (LW)^*] = L ( E[W W'] ) L^* = L L^* = \Sigma
となるので、確率変数 Z := LW の分散共分散行列は Σ ですが、その疑似分散共分散行列は
E[(LW) (LW)'] = L ( E[W W'] ) L' = L L' \neq 0
です。つまり、確率変数 Z := LW は複素多変量正規分布には従いません。
LAPACK
Scipy の線形代数ライブラリ scipy.linalg は、内部的には数値線形代数のための数値解析ソフトウェアライブラリである LAPACK (Linear Algebra PACKage) を利用する Low-level LAPACK functions (scipy.linalg.lapack) に基づき実装されています。cf. scipy.linalg.get_lapack_funcs。
LAPACK の複素数型と NumPy/SciPy の複素数型
LAPACK の公式リファレンス Data Types and Precision によると、複素数型に関しては Fortran data type COMPLEX*16 により実装されています。詳細は分かりませんが、多くの環境で倍精度複素数による計算を可能にしているものと思われます。私の理解では、これは NumPy/SciPy の complex128 型に対応する型だと認識しています。したがって、python 側で complex256 型を利用した場合は、結果的に complex128 型に精度が落ちるものと思われます。詳しくは、NumPy の公式リファレンス Scalars を参照してください。
LAPACK による Cholesky 分解
正定値エルミート行列を Cholesky 分解するには、scipy.linalg.cholesky を利用します。今回、Cholesky 分解 Σ = L L* に際しては、下三角行列 L を利用するので lower=True を指定します（lower=False とすることにより上三角行列 U を利用して Cholesky 分解 Σ = U* U することもできます）。内部的には LAPACK の potrf() を呼び出しています。
正定値エルミート行列を Cholesky 分解するには、scipy.linalg.cho_factor を用いることもできます。scipy.linalg.cho_solve は scipy.linalg.cholesky とは異なり返却する下三角行列 L の下三角部分以外の要素（対角成分より左上の要素）に対して 0 を代入しません。したがって、正確には希望する下三角行列 L を返却する訳ではありませんが、次述の scipy.linalg.cho_factor の引数に渡すことにより、Σ x = y 型の線形方程式を解くことができます。scipy.linalg.cholesky も scipy.linalg.cho_factor も、内部的には scipy.linalg._cholesky を呼び出しています。
LAPACK による線形方程式解
線形方程式 $Lw = (z - \mu)$ を w に関して解くには、L が下三角行列であることを前提に scipy.linalg.solve_triangular を用います。下三角行列 L を利用する場合は lower=True を指定するのを忘れないようにしてください。内部的には LAPACK の trtrs() を呼び出しています。
線形方程式 $Lw = (z - \mu)$ を w に関して解くには、scipy.linalg.cho_solve を用いることもできます。引数には scipy.linalg.cho_factor の返却値を利用することができます。もちろん、scipy.linalg.cholesky の返却値を利用することもできます。やはり、下三角行列 L を利用する場合は lower=True を指定するのを忘れないようにしてください。内部的には LAPACK の potrs() を呼び出しています。
cho_factor + cho_solve の組み合わせの方が、cholesky + solve_triangular の組み合わせよりも、計算量が良くなるかについては、検証していません。
複素多変量正規分布サンプラー
コード
ComplexMultivariateNormal のインスタンス生成時に 1 回だけ scipy.linalg.cholesky を呼び出し、分解結果 (下三角行列 L) をインスタンス変数としてインスタンスに保持します。これ以降は、二度と Cholesky 分解してはいけません。
敢えて、対数確率密度 logpdf に関しては cho_solve を、確率密度 pdf に関しては cho_solve を使っています（もちろん、どちらの関数を使用しても問題ないですし、pdf の結果は logpdf の結果を用いて exp(logpdf) と算出するのが自然です）。
複素数型の精度に関しては、結果的に complex128 型へと精度が低下するものと思いますが、complex256 型を利用することとします。
import numpy as np
import scipy as sp
from numpy import random
from scipy import linalg
class ComplexMultivariateNormal:
def __init__(self, mean, cov):
# cov = L @ L.T.conj(); L: lower triangular
self.N = np.shape(cov)[0]
self.mean = mean
self.L_lower = True
self.L = sp.linalg.cholesky(cov, lower = self.L_lower)
self.L_det = np.abs(np.prod(np.diag(self.L)))
def logpdf(self, z):
z = np.array(z, dtype = np.complex256)
w = sp.linalg.solve_triangular(self.L, z - self.mean, lower = self.L_lower)
return - self.N * np.log(np.pi) - 2 * np.log(self.L_det) - (np.conj(w.T) @ w).real
def pdf(self, z):
z = np.array(z, dtype = np.complex256)
w = sp.linalg.cho_solve((self.L, self.L_lower), z - self.mean)
return 1.0 / (np.power(np.pi, self.N) * np.power(self.L_det, 2)) * np.exp(-(np.conj(w.T) @ w).real)
def sample(self, size):
w = np.reshape((np.sqrt(0.5) * np.random.standard_normal(self.N * size)) + (np.sqrt(0.5)*1.0j * np.random.standard_normal(self.N * size)), (self.N, size))
# w = np.reshape(np.random.standard_normal(self.N * size), (self.N, size)) ## NOT CORRECT
return ((self.L @ w).T + np.tile(self.mean, (size, 1)))
def det(self):
return np.power(self.L_det, 2)
分散共分散行列の行列式の導出において、複素 N × N 正定値 Hermite 行列の固有値は正の実数であること利用して計算していることに注意してください。
テストコード
テストケースに乱数を使うのは、本来のテストの在り方として良くないという意見もあるかとは思いますが、今回は乱数を利用してテストケースを書くことにします。
最低限のクリアすべきテストケースとして、例えば
- 経験分散共分散行列が真の分散共分散行列に近づいているかどうか、
- 経験疑似分散共分散行列が真の疑似分散共分散行列 (複素多変量正規分布の場合は零行列)に近づいているかどうか
などを考えることができます。
ここでは、定常な複素自己回帰モデルに従う、連続する N 個観測量の分散共分散行列 ca.variance(xi, N) を真の分散共分散行列 cov として利用することにします。
import unittest
import numpy as np
import complex_multivariate_normal as cmn
import complex_autoregressive as ca
class Test(unittest.TestCase):
def test_cov(self):
xi = [0.5 + 0.5j, -0.5 + 0.5j]
N = 60
cov = ca.variance(xi, N).astype(np.complex128)
mean = np.zeros(N)
cmn_instance = cmn.ComplexMultivariateNormal(mean, cov)
size = 100000
sample = cmn_instance.sample(size)
sample_cov = (sample.T @ np.conj(sample)) / size
self.assertLess(np.sum(np.abs((sample_cov - cov)**2)), 0.001 * np.sum(np.abs((cov)**2)))
def test_pseudo_cov(self):
xi = [0.3 + 0.5j, -0.5 + 0.6j]
N = 60
cov = ca.variance(xi, N).astype(np.complex128)
mean = np.zeros(N)
cmn_instance = cmn.ComplexMultivariateNormal(mean, cov)
size = 100000
sample = cmn_instance.sample(size)
sample_cov = (sample.T @ np.conj(sample)) / size
sample_pseudo_cov = (sample.T @ sample) / size
self.assertLess(np.sum(np.abs((sample_pseudo_cov)**2)), 0.001 * np.sum(np.abs((cov)**2)))
複素 P 次自己回帰モデルの複素 AR 根 (絶対値が 1.0 未満で定常) を xi = [xi_1, ..., xi_P] として、xi から分散共分散行列 cov を計算して ComplexMultivariateNormal のインスタンスを生成します。
import numpy as np
import complex_multivariate_normal as cmn
class ComplexAutoregressive(cmn.ComplexMultivariateNormal):
def __init__(self, mean, cov = None, xi = None):
self.xi = xi
if cov is None:
N = np.shape(mean)[0]
cov = variance(xi, N)
super().__init__(mean, cov)
def variance(xi, N):
## 省略
ComplexMultivariateNormal.sample で標準複素正規分布ではなく標準実正規分布を利用してサンプリングすると、疑似分散共分散行列に関するテスト test_pseudo_cov で AssertionError が出ることを確認してください。
======================================================================
FAIL: test_pseudo_cov (__main__.Test)
----------------------------------------------------------------------
Traceback (most recent call last):
File "test_complex_multivariate_normal.py", line 45, in test_pseudo_cov
self.assertLess(np.sum(np.abs((sample_pseudo_cov)**2)), 0.001 * np.sum(np.abs((cov)**2)))
AssertionError: 130.8844907571622 not less than 2.2458458139014605