最小固有値の分布の精密計算とカイ2乗近似ppt

固有値分布論の動向
– 最小固有値分布の精密計算と χ2 近似 –
橋口博樹 ∗
宮崎英子 ∗∗
∗∗
埼玉大学
∗∗
杉山高一 ∗∗
中央大学
2009/12/5, 6
1 / 35
概要
✓
✏
1. ダイジェスト（２つの成果）
2. 固有値分布論の広がり
3. ゾーナル多項式からジャック多項式へ
4. 最小固有値の精密計算と χ2 近似
✒
✑
2 / 35
✓
(ダイジェスト 1/2) χ2 近似と精密分布の比較
(
)
W ∼ W4 51, diag(52 , 42 , 32 , 12 ) での最小固有値 ℓ4 /51 の分布
✒
✏
✑
2
2
1.5
1.5
1
1
0.5
0.5
1
0.5
1.5
2
0.5
(a) 精密分布
(b) 自由度 48 の
ϵ = max Fℓm/n(x) − Fχ2
0≤x<∞
1
n−m+1
1.5
χ2
2
分布
(x) ≈ 6.93 × 10−16
• k = 92 までの次数のゾーナル多項式が必要であり，1 次から 92 次ま
での分割の総数
∑92
j=1
j
#P4 は約 1.67 × 106
3 / 35
ダイジェスト (2/2) 計算量
m 変量， k 次までのジャック多項式を生成
Koef and Edelman(2006) のジャック多項式生成の計算量
(
)
√
O m exp(2π 2k/3)
Hashiguchi et al. (2000) のゾーナル多項式生成をジャック多
項式の生成へ拡張し，計算量を導出
(
√
m3
O
exp(2π 2k/3)
k
)
m2 < k なら，Hashiguchi et al. (2000) の拡張版が効率的
4 / 35
✓
✏
固有値分布論のイントロダクション
✒
✑
5 / 35
固有値分布論の広がり
今後の動向，展望 (「ランダム行列の広がり」[永尾他 (2007)])
1. 量子カオス系のエネルギー準位統計と数論における L 関数の
ゼロ点分布
2. 組合せ論とそれに関連した非平衡統計力学への応用
3. 複雑ネットワークに関係した社会学や生物学への応用や通信
理論などの工学的応用
工学への応用
Alex Grant (2002), Rayleigh Fading Multi-Antenna Channels,
EURASIP Journal on Applied Signal Processing, 316–329
Sugiyama (1966), (1967), (1972)
Lin et al. (2008), Ergodic Capacity of LTE Downlink Multiuser MIMO
Systems, Proc. IEEE Communications Society, ICC 2008
Sugiyama (1967)
RATNARAJAH et al. (2005), EIGENVALUES AND CONDITION
NUMBERS OF COMPLEX RANDOM MATRICES, SIAM J. MATRIX
ANAL. APPL., 441–456
Sugiyama (1967), (1970), (1972)
6 / 35
✓
はじめに
✏
固有値分布論とゾーナル多項式
✒
✑
ゾーナル多項式は，多変量正規母集団下で共分散行列の分
布，特に固有値分布を記述する重要な多項式．
James (1960): 一般線形群と直交群に関する表現論を利用し
て導出．
共分散行列の固有値関する多くの統計量が，ゾーナル多項式
を使って導出 [James (1964), Constantine (1963)]．
しかし数値計算の困難性から，1960, 1970 年代での高次ゾー
ナル多項式の計算実現は，Sugiyama(1979) の２変数で 200
次の場合のみ
Sugiura(1990): ２変量の場合に最大固有値，最小固有値の分
布の数値計算
１２変量，１２次までは，McLaren (1976) でプログラム実装
３変量以上のゾーナル多項式の数値計算は，Hashiguchi and
Niki (1997)，Hashiguchi, Nakagawa and Niki (2000) で実装．
Hashgiuchi and Niki(2006): 最大，最小固有値の分布の数値
計算
7 / 35
はじめに (2/2)
✓
✏
ゾーナル多項式からジャック多項式へ
✒
✑
Jack(1970) は，James(1960) のゾーナル多項式を拡張しパラ
メータ α 付きの対称式を導入
α = 2 のとき，ジャック多項式はゾーナル多項式
α = 1 のとき，シュアー多項式 (複素ゾーナル多項式)
α = −1 のとき，拡延対称式
Jack (1972) では，拡延対称式からベキ和多項式への相互変換
方法が記述
⇐ 中川，仁木 (1991) の多変量拡延対称式，多変量ベキ和
多項式の相互変換アルゴリズムの一変量版
Stanley(1989): ジャック多項式の組合せ論的な性質を解明
⇒ Koef and Edelman(2006): ジャック多項式の計算方法，固
有値分布の計算
8 / 35
Jack(1970): ジャック多項式の定義 (1/2)
Y : m × m 対称行列
y1 , . . . , y m: Y の固有値
k
Pm
: m 個以下の要素からなる k 次の分割全体
P23 = {(3), (2, 1)}
[注意 (3, 0) = (3)]
..............................................................
k
J κ (Y; α) for κ = (κ1 , κ2 , . . . , κ m) ∈ P m
: ジャック多項式
J κ′ (Y; α) = α m
∂m
I − t 1 Y κ1 − t 2 Y κ2 − · · · − t mY κm
∂t 1 ∂t 2 · · · ∂t m
−1/α
( t 1 = t 2 = · · · = t m = 0)
..............................................................
9 / 35
Stanley(1989): ジャック多項式の定義 (2/2)
任意の分割 κ に対して、以下の 3 つを満たす対称式（ジャック多
項式） J κ = J κ (Y; α) が唯一存在する．
(P1) (直交性) κ
µ ならば，〈J κ , J µ 〉 = 0
∑
(P2) (上三角性) J κ = µ vκµ mµ と書いたとき、 µ ≤ κ でないな
らば，vκµ = 0 である．
(P3) (正規性) | κ |= m、 J κ における y1 y2 · · · y m の係数 vκ(1m) は m!
に等しい
✓
✏
J κ (Y; α) をジャック多項式という．
✒
(注意) J κ (Y; α) は Y の固有値 y1 y2 · · · y m の同次対称式．
✑
10 / 35
ジャック多項式の性質 [Stanley(1989)]
漸化式 [ J κ (y1 , . . . , y m; α) = J κ (Y; α)]
✓
J κ (y1 , . . . , y m; α) =
∑
µ≤κ
✏
|κ/µ|
J µ (y1 , . . . , y m−1 ; α) y m
βκµ ,
✒
✑
m 変数 k 次のジャック多項式は，それより小さい次数 ( k − 1
以下)，少ない変数 ( m − 1) のジャック多項式で計算できる
ラマヌジャンの k 次分割に関する公式
( √
)
1
P kk ∼
exp
π
2k/3
√
4k 3
11 / 35
ジャック多項式生成の計算量: Koef and Edelman(2006)
{
}
k
ジャック多項式 J κ (Y; α) | κ ∈ P m
の計算量
すべての分割のジャック多項式数 ( s) × 一つのジャック多項式
生成 ( sm): ✓
(
(
))
√
s2 m = O m exp 2π 2k/3
✒
✏
✑
ただし， s は k までの分割の総数（の近似）
s =
k
∑
j=1
j
Pm ≤
k
∑
j=1
j
j
P ∼k
(
( √
))
= O exp π 2k/3
( √
)
1
exp
π
2k/3
√
4k 3
12 / 35
基本対称式の定義
✓
(目標)
✒
ジャック多項式 J κ (Y; α) を基本対称式で展開
✏
✑
e1 , . . . , e m: 基本対称式
e1
=
y1 + y2 + · · · + y m
e2
em
=
=
y1 y2 + y2 y3 + · · · + y m−1 y m, . . . ,
y1 y2 · · · y m
k
Eκ (Y): 分割 κ = (κ1 , κ2 , . . . , κ m) ∈ P m
の基本対称式
Eκ (Y) = eκ1 −κ2 eκ2 −κ3 · · · eκm−1 −κm eκmm
1
2
m−1
Eκ (Y) の次数は k:
Y : 2 次 E(2,1) (Y) = (y1 + y2 )2−1 (y1 y2 ) = y21 y2 + y1 y22
13 / 35
✓
q[κ, µ] の漸化式とその項数
J κ (Y; α) =
∑
✏
q[κ, µ] Eµ (Y)
k
µ∈P m
✒
✑
もし κ ≻ µ かつ d(κ; α)
d(µ; α),
m = 2 の場合，2 項間の漸化式
α
(ν1 + 2) (ν1 + 1) q[κ, (µ1 + 1, µ2 − 1)]
q[κ, µ] =
d(µ; α) − d(κ; α)
m = 3 の場合，4 項間の漸化式
α
{
q[κ, µ] =
d(µ; α) − d(κ; α)
(ν1 + 2)(ν1 + 1) q[κ, (µ1 + 1, µ2 − 1, µ3 )]
+ (ν2 + 2)(ν2 + 1) q[κ, (µ1 , µ2 + 1, µ3 − 1)]
+ 3 (ν1 + 1)(ν2 + 1) q[κ, (µ1 + 1, µ2 , µ3 − 1)]
}
where νi = µi − µi+1 (i = 1, 2, 3),
14 / 35
主定理
ジャック多項式を基本対称式で展開したときの係数の漸化式
の項数を T とする，


t(t + 1)(4t − 1)



1+
m = 2t



6
T = 


t(t + 1)(4t + 5)



m = 2t + 1
1 +
6
( 3)
∼ O m
{
}
k
ジャック多項式 J κ (Y; α) | κ ∈ P m
の計算量
✓
k (
∑
j=1
✒
✏
) ( j )
j
k 2
≤ km3 P m
m3 P m × P m
(
( √
)
m3
exp π 2k/3
∼ O
k
)
✑
15 / 35
計算量比較
m 変量， k 次までのジャック多項式を生成
Koef and Edelman(2006) のジャック多項式生成の計算量
(
)
√
O m exp(2π 2k/3)
Hashiguchi et al. (2000) のゾーナル多項式生成をジャック多
項式の生成へ拡張し，計算量を導出
(
√
m3
O
exp(2π 2k/3)
k
)
m2 < k なら，Hashiguchi et al. (2000) の拡張版が効率的
16 / 35
✓
✏
最小固有値分布の精密計算と χ2 近似
✒
✑
最小固有値分布の精密計算
χ2 近似
17 / 35
もう一つのジャック多項式
✓
(tr Y) k =
∑
✏
Cκ(α) (Y)
k
κ∈P m
✒
✑
J κ (Y; α) と Cκ(α) (Y) の関係
Cκ(α) (Y) =
ただし，
jκ =
α|κ| (|κ|)!
∏
jκ
J κ (Y; α)
h∗κ (i, j) hκ∗ (i, j)
(i, j)∈κ
h∗κ (i, j) = k′ − i + α(ki − j + 1), hκ∗ (i, j) = k′ − i + 1 + α(ki − j) は
j
j
それぞれ the upper と lower hook length at (i, j) ∈ κ である．
18 / 35
最小固有値分布の数値計算
α=2
=⇒
zonal polynomial
W m(n, Σ) での最
大，最小固有値
の分布計算
complex zonal
polynomial
CW m(n, Σ) での
最小固有値の分
布計算
(α)
Jack function: Cκ
α=1
=⇒
W m(n, Σ) での最小固有値の分布計算とその限界
CW m(n, Σ) での最小固有値の分布計算とその限界
✓
✒
✏
最小固有値の分布の χ2 近似
✑
19 / 35
Khatri (1972): ウィシャート行列の最小固有値分布
W ∼ W m(n, Σ) ( n ≥ m).
ℓ m: ウィシャート行列 W の最小固有値
j = (n − m − 1)/2: 非負整数
✓
✏
jm
( x
) ∑
x k ∑∗
−1
Pr[ℓ m < x] = 1 − exp − tr Σ
Cκ (Σ −1 )
k k!
2
2
κ
k=0
✒
∑
ただし，和の記号 ∗ は，
✑
κ
{
k
κ = (κ1 , . . . , κ m) ∈ P m
| κ1 ≤ j
}
の条件 (κ1 ≤ j) を満たす分割にわたって和を取る．
20 / 35
Ratnarajah et al. (2005): 複素ウィシャートの最小固有値分布
W ∼ CW m(n, Σ) (n ≥ m)．
ℓ˜m: W の最小固有値
✓
✏
(
)
Pr[ℓ˜m < x] = 1 − exp −x tr Σ −1
m(n−m)
∑
k=0
✒
ただし，Cκ (Y) =
{
Cκ(1) (Y)．
∑∗
κ
x k ∑∗
Cκ (Σ −1 ),
k! κ
✑
は
}
k
| κ1 ≤ n − m .
κ = (κ1 , . . . , κ m) ∈ P m
の条件下で和をとる．
21 / 35
✓
✏
最小固有値 ℓ4 , ℓ˜4 の密度関数
✒
2.5
n = 11
2
2.5
n = 21
1.5
1
1
0.5
0.5
1
n = 10
2
1.5
0.5
✑
1.5
2
Figure: W4 (n, diag(1, 2, 3, 4))
n = 20
0.5
1
1.5
2
Figure: CW 4 (n, diag(1, 2, 3, 4))
22 / 35
数値計算の限界に関する考察 (1/2)
分割の最大次数
✓
実数:
✒
m(n − m − 1)
2
✏
複素: m(n − m)
✑
( x
)
Pr[ℓ m < x] = 1 − exp − tr Σ −1
2
(n−m−1)/2
∑ m
(
)
Pr[ℓ˜4 < x] = 1 − exp −x tr Σ −1
m(n−m)
∑
k=0
k=0
x k ∑∗
2 k k!
Cκ (Σ −1 )
κ
x k ∑∗
Cκ (Σ −1 ),
k! κ
23 / 35
数値計算の限界に関する考察 (2/2)
k
|P m
|
k = 100
k
k
t = |P m
|(|P m
| + 1)/2
(approx.)
m=2
51
1326
(1.3 × 103 )
m=3
884
391170
(3.9 × 105 )
m=4
8039
32300703
(3.2 × 107 )
t ≤ 105
Real
Complex
k
m(n − m − 1)/2 ≤ k
m(n − m) ≤ k
m=2
k ≥ 100
n ≥ 103
n ≥ 52
m=3
k ≤ 70
n ≤ 50
n ≤ 26
m=4
k ≤ 35
n ≤ 23
n ≤ 12
24 / 35
✓
✏
最小固有値分布の精密計算と χ2 近似
✒
✑
最小固有値分布の精密計算
χ2 近似
25 / 35
最小固有値分布の χ2 近似
W ∼ W m(n, Σ), Σ = diag(λ1 , . . . , λm)
1.
1 −1
(λ + λ−1 + · · · + λ−1
m ) とおく．このとき ℓ m の密度関
2
2 1
数は，(2) 式で近似できる．
λ=
1
f (ℓ m) =
2.
Γ
(
λ 2 (n−m+1)
1
(n −
2
m + 1)
1
(n−m+1)−1
2
) ℓm
e−λℓm .
(2)
ℓ m/λm の分布は，λm が他の固有値と離れ n が十分大きい場
合，自由度 n − m + 1 の χ2 分布に従う．
✓
f (ℓ m) =
✒
1
1
(2λm) 2 (n−m+1)
1
(n−m+1)−1
2
ℓm
[
exp −
ℓm
2λm
✏
]
.
(3)
✑
26 / 35
✓
χ2 近似と精密分布の比較
(
)
W ∼ W4 51, diag(52 , 42 , 32 , 12 ) での最小固有値 ℓ4 /51 の分布
✒
✏
✑
2
2
1.5
1.5
1
1
0.5
0.5
1
0.5
1.5
2
0.5
(a) 精密分布
(b) 自由度 48 の
ϵ = max Fℓm/n(x) − Fχ2
0≤x<∞
1
n−m+1
1.5
χ2
2
分布
(x) ≈ 6.93 × 10−16
• k = 92 までの次数のゾーナル多項式が必要であり，1 次から 92 次ま
での分割の総数
∑92
j=1
j
#P4 は約 1.67 × 106
27 / 35
最小固有値分布の χ2 近似
W ∼ W4 (51, diag(52 , 42 , 32 , 1)) の最小固有値 ℓ4 /51
パーセント点
累積確率
α
精密分布 (xα )
χ2
Fχ2 (xα )
48
0.99
0.95
0.05
0.01
1.437
1.271
0.645
0.549
1.444
1.277
0.648
0.552
48
0.989
0.947
0.04784
0.00947
✓
✏
χ2 近似は 99%，95%，5%，1% パーセント点で 3 桁の精度が
あることが分かる．
✒
✑
28 / 35
χ2 近似とモンテカルロシミュレーション
シミュレーション条件
n = 100
m = 4，Σ = diag(42 , 32 , 1.52 , 1)，diag(42 , 32 , 22 , 1)，
diag(52 , 42 , 32 , 1)
x∗α : シミュレーションでの 100α パーセント点
(106 個の並び替え)
Fχ2
✓
n−m+1
(x∗α ): χ2
n−m+1
での累積確率
✏
シミュレーション回数が 106 であるので，その精度は 3 桁程
度である．χ2 とシミュレーション結果が 2, 3 桁合っている．
✒
97
✑
29 / 35
最小固有値分布の χ2 近似シミュレーション結果 (1/4)
Table: 106 回のシミュレーションでのパーセント点 ( x∗α ) と χ2 との比較
97
(1)
W ∼ W4 (100, Σ) での W の最小固有値 ℓ4 のパーセント点
Σ
diag(52 , 42 , 32 , 1)
diag(42 , 32 , 22 , 1)
α
x∗α
Fχ2 (x∗α /n)
x∗α
Fχ2 (x∗α /n)
97
0.99
0.975
0.95
0.90
0.50
0.10
0.05
0.02
0.01
132.182
125.962
120.794
115.019
96.130
79.460
75.129
71.518
67.427
0.9898
0.9743
0.9487
0.8976
0.4941
0.0975
0.0486
0.0243
0.0096
97
131.772
125.569
120.419
114.667
95.846
79.213
74.876
71.283
67.208
0.9891
0.9729
0.9462
0.8936
0.4859
0.0940
0.0465
0.0232
0.0091
x∗ : 106 回のシミュレーションでのパーセント点
30 / 35
最小固有値分布の χ2 近似シミュレーション結果 (2/4)
Table: 106 回のシミュレーションでのパーセント点 ( x∗α ) と χ2 との比較
97
(1)
W ∼ W4 (100, Σ) での W の最小固有値 ℓ4 のパーセント点
Σ
diag(42 , 32 , 22 , 1)
diag(42 , 32 , 1.52 , 1)
α
x∗α
Fχ2 (x∗α /n)
x∗α
Fχ2 (x∗α /n)
97
0.99
0.975
0.95
0.90
0.50
0.10
0.05
0.02
0.01
131.772
125.569
120.419
114.667
95.846
79.213
74.876
71.283
67.208
0.9891
0.9729
0.9462
0.8936
0.4859
0.0940
0.0465
0.0232
0.0091
97
130.922
124.803
119.719
114.032
95.359
78.863
74.563
70.971
66.943
0.9876
0.9698
0.9413
0.8859
0.4719
0.0893
0.0439
0.0217
0.0085
x∗ : 106 回のシミュレーションでのパーセント点
31 / 35
最小固有値分布の χ2 近似シミュレーション結果 (3/4)
Table: 106 回のシミュレーションでのパーセント点 ( x∗α ) と χ2 との比較
91
(2)
W ∼ W10 (100, Σ) での W の最小固有値 ℓ10 のパーセント点
λ1 = 102 , λ2 = 92 , λ3 = 82 , λ4 = 72 , λ5 = 62 , λ6 = 52
diag(λ1 , . . . , λ6 , 4.52 , 42 , 32 , 1) diag(λ1 , . . . , λ6 , 42 , 32 , 22 , 1)
100α
x∗α
Fχ2 (x∗α /n)
x∗α
Fχ2 (x∗α /n)
91
0.99
0.975
0.95
0.90
0.50
0.10
0.05
0.02
0.01
125.007
118.962
113.921
108.318
90.017
73.909
69.731
66.225
62.322
0.9895
0.9738
0.9476
0.8959
0.4905
0.0957
0.0475
0.0235
0.0093
91
124.486
118.505
113.505
107.924
89.699
73.675
69.507
66.014
62.098
0.9886
0.9720
0.9447
0.8912
0.4811
0.0924
0.0456
0.0225
0.0088
x∗ : 106 回のシミュレーションでのパーセント点
32 / 35
最小固有値分布の χ2 近似シミュレーション結果 (4/4)
Table: 106 回のシミュレーションでのパーセント点 ( x∗α ) と χ2 との比較
91
(2)
W ∼ W10 (100, Σ) での W の最小固有値 ℓ10 のパーセント点
λ1 = 102 , λ2 = 92 , λ3 = 82 , λ4 = 72 , λ5 = 62 , λ6 = 52
diag(λ1 , . . . , λ6 , 42 , 32 , 22 , 1) diag(λ1 , . . . , λ6 , 42 , 32 , 1.52 , 1)
100α
x∗α
Fχ2 (x∗α /n)
x∗α
Fχ2 (x∗α /n)
91
0.99
0.975
0.95
0.90
0.50
0.10
0.05
0.02
0.01
124.486
118.505
113.505
107.924
89.699
73.675
69.507
66.014
62.098
0.9886
0.9720
0.9447
0.8912
0.4811
0.0924
0.0456
0.0225
0.0088
91
123.644
117.755
112.825
107.289
89.222
73.318
69.185
65.739
61.901
0.9870
0.9689
0.9397
0.8831
0.4668
0.0874
0.0429
0.0211
0.0083
x∗ : 106 回のシミュレーションでのパーセント点
33 / 35
まとめ
1. Hashiguchi, Nakagawa and Niki (2000) の real zonal
polynomials の計算方法を Jack functions の場合へ拡張し，
complex zonal polynomials を 2 変数から 4 変数で求めた．
2. Jack functions を基本対称式で展開するための計算量を求
めた．
(
)
O
√
m3
exp(2π 2k/3)
k
Koef and Edelman (2006) の方法と比べ， m2 < k なら，
Hashiguchi et al. (2000) の拡張版が効率的
3. 最小固有値分布の精密計算を行った．また分割の総数により
数値計算の限界について考察を行った．
最大固有値の分布計算には一般化超幾何関数の計算が必要と
なるので，数式処理言語での開発を今後行う予定である．
4. 最小固有値分布の χ2 を求め，精密分布，シミュレーション
との比較を行った．
34 / 35
Conclusions
本研究では次のことを示した．
1. Hashiguchi, Nakagawa and Niki (1998) の real zonal
polynomials の計算方法を Jack functions の場合へ拡張し，
complex zonal polynomials を 2 変数から 4 変数で求めた．分
割の次数では 100 以下である．
2. Jack functions を基本対称式で展開するための漸化式の項数
が，O(m3 ) であることを示した．
monomial symmetric functions での展開では，漸化式の項数
は O(k2 ) である．
Jack 多項式生成のための計算量を見積もり，Koef and
Edelman (2006) と比較する必要がある．
3. 固有値分布の数値計算を行った．2 変数で実数のウィシャー
ト行列では 100 次以上の自由度のものを計算できるが，4 変
数では 30 程度が計算可能であることが分かった．複素版で
は実数版の半分で 4 変数で 20 程度である．4 変数で自由度
20 の場合に 68 次までの分割が必要で実際に計算できている
ことも確認した．
最大固有値の分布計算には一般化超幾何関数の計算が必要と
なるので，数式処理言語での開発を今後行う予定である．
35 / 35

Download Report