non-central complex Wishart 行列の最大固有値の CDF の数値計算
Fadil Habibi Danufane, 小原功任 (金沢大学), 高山信毅 (神戸大学), Constantin Siriteanu (大阪大学)

2017 年 11 月 09 日

問題の背景多アンテナ無線通信システム (MIMO=multiple-input andmultiple-output) における通信不能確率 (outage probability) の数値評価法の確立．

NT 個の送信用アンテナと NR 個の受信用アンテナをもつ MIMO
無線通信システムのチャンネル行列 H は，NT × NR 複素数値ランダム行列 1 であって，H の列ベクトルは共分散 Σ の複素数値ガウス分布にしたがう. E [X ] = M とする 2 . 行列 S = Σ−1 HH ∗は，非心 Wishart 分布にしたがうことが知られている. このとき，
s = min{NT , NR }, t = max{NT , NR } とおき，ϕ1 < ϕ2 < · · · < ϕsを行列 S の非負固有値とする．

12cmvnorm
平均は列毎に異なって良い累積確率密度関数
定理 ([KA], 2003)
行列 Σ−1 MM ∗ が正固有値 0 < λ1 < λ2 < · · · < λs をもつならば，S の最大固有値 ϕs について，閾値 x に関する通信不能確率(outage probability) の累積確率密度関数は∑exp(− si=1 λi )∏Pr (ϕs ≤ x) =det Ψ(x)Γ(t − s + 1)s 1≤i<j≤s (λj − λi )となる．ここで，Ψ(x) は行列値関数であり，その (i, j) 成分は，
∫ xt−iHt−s+1 (x, λj ) =y t−i exp(−y ) 0 F1 ( ; t − s + 1; y λj ) dy0[KA] M. Kang, M. S. Alouini, Largest Eigenvalue of Complex
Wishart Matrices and Performance Analysis of MIMO MRC
Systems, IEEE Journal on Selected Areas in Communications 21(2003), 418–426.

合流型超幾何級数 0 F1
定義 合流型超幾何級数 0 F1 を以下で定める．
0 F1 (; c; x) =∞∑n=0xn(c)n n!

ここで，(c)n = Γ(c + n)/Γ(c) はポッホハマーの記号．

また，次の微分方程式をみたす．
{()}
ddxx+c −1 −x •u =0dxdx

累積確率密度関数 (cont.)
問題

累積確率密度関数Pr (ϕs ≤ x) =∑exp(− si=1 λi )∏det Ψ(x)Γ(t − s + 1)s 1≤i<j≤s (λj − λi )を数値的に評価せよ．

行列 Ψ(x) の各成分を数値的に評価

先行研究 ([KA])Ψ(x) の各成分を，クンマー超幾何関数 1 F1 や変形ベッセル関数などで表す．
=⇒ 級数展開

0 F1 との関係

補題∂ k
H (x, y ) = x k e −x 0 F1 ( ; n; xy )∂x n

命題 (e.g. [DM])0 F1 (n; z) ∼√√Γ(2n − 1)
exp(2 z)(4 z)−n+1/2 ,Γ(n − 1/2)

(z → +∞)[DM] A. Dixit, V. H. Moll, The integrals in Gradshteyn andRyzhik. Part 28: The confluent hypergeometric function andWhittaker functions. Scientia, Series A:Mathematical Sciences 26(2015), 49–62.

漸近的評価定理 領域 x ≫ y ≫ 0 において，x → +∞ となるとき，Hnk (x, y ) は漸近的にHnk (x, y ) ∼Γ(2n − 1) √ −n+1/2(4 y )·2·exp(y ) exp(−P(s0 ))Γ(n − 1/2)である．ただし，√P(s) = −(s − y )2 + (1 + 2k − n + 1/2) log s√1 √s0 =( y + y + 2(1 + 2k − n + 1/2))2とする．

(2π′′P (s0 ))1/2

級数表示

定理
関数 Hnk (x, y ) は，以下の級数表示をもつ．
∞Hnk (x, y ) =(−1)q (k + 1)p+qx k+1 ∑x p+q y p
k +1(n)p (1)p (1)q (k + 2)p+qp,q=0

原点に十分近い点では，級数によって近似値を得る

微分方程式系

定理 関数 u = Hnk (x, y ) は，微分方程式系{θy (θy + n − 1) + y (θx − θy − k − 1)} • u = 0,
(θx − θy − k − 1 + x) θx • u = 0.
∂∂をみたす．ただし，θx = x ∂x,θy = y ∂yである．
この方程式系のホロノミックランクは 4 である．

ホロノミック勾配法 (HGM)

ホロノミック関数 u(x1 , . . . , xn ) がランク r の斉次線形偏微分方
程式系をみたすとする．このときグレブナー基底を用いると，うまい微分作用素 t2 , . . . , tr によって，方程式系を∂
F = Pi (x1 , . . . , xn )F∂xi(i = 1, 2, . . . , n)と書き換えることができる (パッフィアン方程式)．ここで，F = t (u, t2 • u, . . . , tr • u) である．

パッフィアン方程式と Runge-Kutta 法を用いた数値計算法をホロノミック勾配法という．

不安定性

Hnk (x, y ) はランク 4 の方程式系をみたすので，グレブナー基底を用いて，∂∂
F = P1 (x, y )F ,F = P2 (x, y )F∂x∂yの形の常微分方程式が導かれる．しかしこの方程式からは，長い区間では，精密な数値的結果を得られない．
目的の関数 Hnk (x, y ) よりも増大度の大きな解が存在不安定性, 例
−1 1 0dF =  0 −1 1  Fdt00 0
解は F 1 = (exp(−t), 0, 0)T , F 2 = (0, exp(−t), 0)T ,F 3 = (1, 1, 1)T . t = 0 での初期値 (1, 0, 0)T で F1 が解.
初期値にすこしエラーを加える: (1, 10−30 , 10−30 )T . すると,tvalue F1 by RK diﬀerence F1 − F1150 1.92827e-22
9.99959e-31
60 8.75556e-27
1.00000e-30
70 1.39737e-30
1.00000e-30
80 1.00002e-30
1.00000e-30

rank 3 方程式の不安定性∂kk −x F ( ; n; xy ) より, H k は,0 1n∂x Hn (x, y ) = x e

L = ((x∂x − k + x)(x∂x − k + x + n − 1) − yx)∂x をみたす. 3
命題 x → +∞ で次の漸近行動を持つ解達が基底.

h1 = 1h2 = exp(−x + 2(xy )1/2 )((xy )−1/2 )−2k+n−1/2 (1 + O((xy )−1/2 ))h3 = exp(−x − 2(xy )1/2 )((xy )−1/2 )−2k+n−1/2 (1 + O((xy )−1/2 ))漸近評価によれば, Hnk は x について O(1). しかし, O(1) ≪ h2 4 .

3formal sol(L,[dx,x],x=infinity);
講演後の訂正: これは誤り. rank 3 の方程式は安定だが, rank 4 の方程式が不安定
4

rank 4 方程式の不安定性
命題 y → +∞ で次の漸近行動を持つ解達が基底.

h1 = (xy )−1/2(1/2+n) exp(−2(xy )1/2 )(1 + O(1/y 1/2 )),h2 = y −k−1 (1 + O(1/y )),h3 = (xy )−1/2(1/2+n) exp(2(xy )1/2 )(1 + O(1/y 1/2 )),h4 = y 1−n+k exp(y )(1 + O(1/y )),計算実験

yRatio1000 7.36595030875893e-452
2000 2.64621603289928e-881
3000 2.67723893601667e-1311
Ratio = (H110 (1/2, y ))/(y 1−n+k exp(y )).
Hnk には h4 成分は無さそう.

よりよい HGM

補題

∂ kH (x, y ) = x k e −x 0 F1 ( ; n; xy )∂x n

補題 f (x, y ) = 0 F1 ( ; n; xy ) は 2 階の斉次線形微分方程式{()}∂∂
xx+ n − 1 − xy f (x, y ) = 0∂x∂xをみたす．

よりよい HGM (cont.)∂定理 ベクトル値関数 F (x, y ) = t (Hnk (x, y ), f (x, y ), x ∂xf (x, y ))は，一階微分方程式

0 x k e −x0∂
1 0
F = P(x, y )F ,P(x, y ) = 0x∂x1−n0yxをみたす．

▶ 初期値は原点近傍の級数表示により求める．

ゲージ変換とは∂ ′∂F = PF, F′ =F∂x∂x
G (x) を行列 (ゲージ変換行列と呼ぶ).
F = GHF′ = G ′ H + G H′ なので,()H′ = −G −1 G ′ + G −1 PG H
解の漸近展開の主部をゲージ変換行列に選ぶ. 解 H は O(1) となり数値解析しやすい.

ゲージ変換u(t, z) = exp(−2tz)f (t 2 , z 2 ) とおく．
∂定理 ベクトル値関数 H(t, z) = t (Hnk (t 2 , z 2 ), u(t, z), t ∂tu(t, z))
は，一階微分方程式∂H = RH∂tをみたす．ここで01
R(t, z) =0t02t 2(k+1) e −t +2tz
0−2tz(2n − 1)2

0
.
1−4tz + 2(n − 1)

数値実験 1

0.6
0.4
0.00.2

CDF of the outage prob.

0.81.0行列 Σ−1 MM ∗ の固有値が {0.1, 0.2, 0.3, 0.4, 0.5} のとき，通信不能確率の累積確率密度関数のグラフは下図のようになる．
((NT , NR ) = (5, 5) が最上部の破線，(9, 5) が最下部の実線)

1012141618x

図: (NT , NR ) = (5, 5), (6, 5), . . . , (9, 5), λ = {0.1, 0.2, 0.3, 0.4, 0.5}

WiFi の解析では...
20NT NR ∼ y = z 2 = λi . つまり y は大きい.
たとえば 5 × 100 で y ∼ 20NT NR = 104 .
ymax∼ 1 or 10 or . . .
ymin
さらに挑戦問題: y ∼ 108 .
しかし, たとえばH10 (108 , 108 ) ≃ 7.749665152657736 × 1043429447
double の最大値 1.797693 × 10308 .
⇒ (1) bigfloat 対応の RK. または (2) さらなる Gauge 変換でdouble で RK.

さらに Gauge 変換 : double 対策

00000  + O(1/t),
R(t, z) ∼  00 −2z(2n − 1) −4z(安定).

g (t) = exp(−t 2 + 2zt + 2(k + 1) log t)()√4(k + 1)
z1+ 1+tm =2z2
Gauge 変換行列diag(g (tm ), 1, 1),t > tmdiag(g (t), 1, 1),t ≤ tma26d-x.c 5 : abs error 10−40 , rel error 10−20 . 初期値は数値積分 6 .
初期値は上の Gauge 変換で shift した値. 内部的に初期値 (∼ 10−200 ) が(1, 1, 1) となるように定数 Gauge 変換.
6NIntegrate[tˆk*Exp[-t]*HypergeometricPFQ[{},{n},t*y],{t,0,x}];5

5 × 100 データ

λ = y = (1000, 9850, 9950, 10050, 18900)
x
18291,18292,18293,18294,18295,
18296,

通信失敗確率 Outage Probability (3 桁程度有効)
9.5131498898356095823246796729448732748e-05
9.7129526345525128081447073098646520549e-05
9.9166868841088132645047214949297442091e-05
0.00010124430281005984620939531913849296543
0.00010336247277432610901154418144577980193
0.00010552219747580335704519387529275255356

18559,
18560,
18561,
18562,
18563,
18564,

0.0098282592564031513975306126792867854095
0.0099656058110287087344059017182109596642
0.010104621126295227580655729963893042756
0.010245322099469326847761622248107088847
0.010387724710229051062645433015770748285
0.010531845711868293978689856959196384411

18565,

0.010677701832260576777072934579875051019

5 × 100 アンテナ77

横軸は ∼ log10 (1/x)

誤差の検出[2053] set_xmin_xmax(5,100);
Lambda=[ 9750 9850 9950 10050 10150 ][2055] auto_ppp(10200);
0.00044597210091649080492296557923849273023[2056] A26_error;[略.
[Sig_digits,8],[Valid_digits,[-1,-2,-2,-2,-2,-2,-2,-3,-3,-2[Assume Sig_digits,14],[then Valid_digits,[5,4,4,4,4,4,4,3,行列の成分に誤差を加えて det の値の変動を計算. これで有効桁を見積もる. 88double では Hnk の精度は 14 桁程度が限界

HnK の HGM, まとめ
1. 微分方程式を見つけよ.
2. 計算したい解は一番増大度の大きい解か? (安定性).
3. Yes なら global に計算可能. No なら local にしか計算できな
い (でも exterpolation には有用).
4. 安定な微分方程式を見つけよ (存在の保証はない).
5. double 範囲の計算のためにうまい Gauge 変換を見つける.
6. (構想) bigfloat または有理数での効率的かつ高信頼性
Runge-Kutta.
7. (構想) Borel-Laplace 変換したものは安定かも. これと一重数値積分を組み合わせる. などなど解析の手法を動員.
8. (構想) 増大度の大きな解を用いて方程式を上三角化, または,解の主部を行列形式で求めて Gauge 変換. 増大度の小さい解を捕まえる?

References[1] F. H. Danufane, K. Ohara, N. Takayama, C.Siriteanu9 :Holonomic Gradient Method for the Distribution Function of
the Largest Root of Complex Non-central Wishart Matrices,preprint. (arXiv:1707.02564 [math.ST])
上記の論文の参考文献表も参照.

9ver 2 より