ホロノミック勾配法 (HGM, Holonomic Gradient Method)
高山信毅 (神戸大学)
参考文献等 hgm OpenXM search.

(ある点 x = a の近傍でなめらかな) 多変数関数f (x) = f (x1 , . . . , xn ). 全方向で多項式係数の常微分方程式
ri∑aik (x)∂ik • f = 0,(i = 1, . . . , n)

(1)

k=0k

∂
∂ fを満す. aik : x = (x1 , . . . , xn ) の多項式. ∂i = ∂x, ∂ik • f = ∂xk.
ii

このような性質を満す多変数関数 (の解析接続) を, holonomic 解析関数と呼ぶ.
例: f = exp(−x2 x12 ) は 2 変数の holonomic 解析関数. 実際,∂1 • f = −2x1 x2 f , ∂2 • f = −x12 f .
1/18

(超関数を含む一般の) Holonomic 関数 f は, 全方向で常微分方程式 (1) をみたしさらにもう少し沢山の多項式係数偏微分方程式系を満す関数.
例: H(x) = 1 (x ≥ 0), H(x) = 0 (x < 0) を Heaviside 関数とする.
x∂x • H = 0. exp(−x2 x12 )H(x1 ) は holononic 関数.
Holonomic 関数 f を零化する線形偏微分作用素全部の集合はholonomic ideal (system).
holonomic 系を特殊関数の代わりにしよう (佐藤幹夫の哲学, 1970 年代).
精密に定義するには, イデアル, Hilbert 関数, 等の概念.

Holonomic 関数の性質
1. f , g が holonomic 関数なら, f ± g , fg (定義できれば) もholonomic 関数.
2. f (θ, t) を n + m 個の変数 θ = (θ , . . . , θ ), t = (t , . . . , t )1n1mについての holonomic 関数. パラメータ付き積分 (例: f が(非正規化) 分布密度関数なら Z は正規化定数)∫Z (θ) =f (θ, t)dt1 · · · dtm
Rm

は (積分についてのある technical な条件のもと) n 変数のholonomic 関数.
証明は D-加群の理論. 微分方程式系を導出するアルゴリズムもあり (Zeilberger∫ 1990, Oaku 1997, ...).
例: Z (θ) = R exp(−θt 2 )H(t)dt は holononic 関数. 実際,(2θ∂θ +∫ 1) • Z = 0.
練習: R exp(−y (x − 1)2 )H(x)dx では?
[ [(Dy+(x-1)^2) (x*Dx+2*x*(x-1)*y)] [(x)] [[(x) (y)] []] 0]
integration ::[2 , [-2*e_*y-2*y*Dy-2*y-1 , -2*y*Dy^2-2*y*Dy-3*Dy-1 ]]練習 以下の中で holonomic 関数はどれ?
1. exp(f (x)) where f is a rational function,
2.

1[Hint] Use Th: Any solution of the ordinary diﬀerential equation
sin x(am (x)∂ m + · · · + a0 (x)) • f = 0, ai ∈ C[x], is holomorphic out of thesingular locus {x | am (x) = 0}.
3.

Γ(x), [Hint] Γ(x) has poles at x = −n, n ∈ N0 .

4.

2x ,

5.

H(x) (Heaviside function),6.

x a where a is a constant,

7.

|x|,
∫ +∞

8.

−∞

exp(−xt 6 − t)dt, x > 0.

Holonomic 分布: ホロノミック関数で (非正規化) 分布密度関数が書ける分布. f (θ, t)/Z (θ) が t 空間での分布密度関数.

Holonomic 関数の差分版. 差分作用素 Ei : Ei • f (yi ) = f (yi + 1).
Mellin 変換: Ei 7→ xi , yi 7→ −xi ∂i − 1 = −∂i xi .
F (xi ) 7→ G (yi ) =∫ ∞0yF (xi )xi i dxif (y1 , . . . , yn ) が差分版 holonomic 関数 (holonomic 数列) とは, fの満す線形差分作用素達をそれぞれ Mellin 変換したとき,holonomic ideal をなすこと.
この時 f は全方向に多項式係数の線形漸化式を満す+α.
u u∑
p1 1 p2 2
例: Z (n; p) = u1 +u2 =n,ui ≥0 u1 !u2 ! とおく.
t −n−1を t で展開したときの t −1 の係数
1 − p1 t − p2 t∫t −n−11f (t, p1 , p2 )dt, f =Z (n; p) =n!2πi C
1 − p1 t − p2 t
En • f (n, p) = tf (n, p), ..., ⇒ [(n + 1)En − (p1 + p2 )] • Z (n; p) = 0(Z は二項展開).
例: 離散 A-超幾何分布 (予稿 7 節) は β について差分 holonomic
関数.
5/18n!Z (n; p) =

Pfaﬃan 方程式とは?
[a2 (x)(F =fdfdx)(d •F =
とおくと, dx

d2dx 2+ a1 (x)0−a0 /a2ddx]+ a0 (x) • f = 01−a1 /a2)
F を満たす. これの多変数化.

δi = ∂i , i = 1, . . . , n. α1 = 0 ∈ Nn0 , αi ∈ Nn0 , i = 2, . . . , r をうまく与えると, c 基底ベクトル1δα • Z α2 δ•Z
·F =··r

δα • Zに対して,∂i • F = Pi F , i = 1, . . . , n

微分 Pfaﬃan 系となる r × r 行列 Pi (x) が存在. この方程式を Pfaﬃan (方程式)系とよぶ.
注: holonomic ideal から Pfaﬃan 系への変換はグレブナー基底で原理的には可能. yang.rr など.

例: Von Mises 分布. exp(θ1 cos(t) + θ2 sin(t))/Z (θ),∫ 2πZ (θ) = 0 exp(θ1 cos(t) + θ2 sin(t))dt. (以下の式は数学セミナーの記事 [20] より , 予稿 3 節 Fisher-Bingham 分布).
()01
∂Z=Z,
θ12θ22 −θ12∂θ1
θ12 +θ22
θ1 (θ12 +θ22 )()θ20∂Z
θ1=Z−2θ2θ1 θ2∂θ2
θ2 +θ 2θ 2 +θ 2
1212

ここで Z = (Z , ∂1 • Z )T , ∂i = ∂/∂θi である.

7/18

例, 2 × 2 分割表, 予稿 7 節 離散 A 超幾何分布
0 0 1 1∑u
A =  1 0 1 0 , Z (β; x) = Au=β,u∈Nn xu!
00 1 0 1
F = (Z , E3 • Z )T , D = x1 x4 − x2 x3 , B = β1 − β2 − β3 ,()−β1 x1 x4 +β3 D−x2 D1(B−1)x(B−1)x33
E1 −x4β3−x2

1E2 −x1(−B−x2

β3 Bx3Dβ2 x1 x4 +β3 x2 x3D)

Z (β = 0, x) = 1.
E1 • F (β1 , β2 , β3 ; x) = F (β1 + 1, β2 , β3 ; x)

. HGM (Holonomic Gradient Method) の 3 ステップ
. 正規化定数 Z に対する holonomic 系を求める. さらに

1holonomic 系を Pfaﬃan に変換する. 理論的方法または D-加群のアルゴリズムを用いる.
. 正規化定数とその偏微分または偏差分を必要な rank まであ2る点で計算する. [c 基底ベクトルの初期値を求める.].3 2 で求めた値を Pfaﬃan で必要な点まで延長する. [初期値問題を解く.]
注: rank r は step 3 の計算量を決める.
注: 最尤推定, e.g., [20] では,∑log f (θ, T ) − N log Z (θ) = log

T ∈ 観測値の N 個の集合

∏ f (θ, T )ZTの最大化を θ について行う.
c 基底ベクトル F の数値がわかれば, Z の Hessian やより高い階数の微分の値もわかる.

∂2 F∂xi ∂xj=∂Pi (β, p)∂xj

F + Pi∂F∂xj(=理論的手法による Pfaﬃan の導出と応用例: 予稿 3 節から 8 節.

∂Pi (β, p)∂xj)+ Pi PjF.

(2)

. Fisher-Bingham 分布 (予稿 3 節)

∑22d 次元球面 S d (r ) = {(u1 , . . . , ud+1 ) | d+1
i=1 ui = r , r > 0} の上でパラメータ xij , yi を持つ分布
d+1∑∑1exp xij ui uj +yi ui  δµ(u; x, y , r )δ :=Z (x, y , r )i=11≤i≤j≤d+1

を Fisher-Bingham 分布と呼ぶ (1980 年代から始まる A.T.A.
Wood 等の研究が先駆). ここで δ = δ(r , u) は∫d 2π (d+1)/2Rd+1 δ(r , u) = r Γ((d+1)/2) をみたす球面に台をもつデルタ関数,Z は次で定義される正規化定数∫d+1∑∑
Z (x, y , r ) =exp 
xij ui uj +yi ui  δ(r , u). (3)Rd+1

1≤i≤j≤d+1

i=1(清, 竹村が発端. その後, 小山民雄ら)
10/18

. 例：神戸市の平均風向 (16 方位) のデータ

Figure: 頻度と密度関数

方位
北
北北東
北東
東北東
東
東南東
南東
南南東
南
南南西
南西
西南西
西
西北西
北西
北北西

ラジアン
1.570760
1.178070
0.785380
0.392690
0.000000
5.890349
5.497659
5.104969
4.712279
4.319589
3.926899
3.534209
3.141519
2.748829
2.356139
1.963450

頻度
0.00
0.00
0.24
0.24
0.03
0.12
0.03
0.00
0.02
0.00
0.19
0.070.03
0.00
0.00
0.00

出典: 小山, 中山, 西山, 高山: 高次元における FB 分布の MLE への hgd 法の適用について (スライド), 2011–2013.

. 天球の星, S での Fisher-Bingham 分布による fitting
2

出典: 中山, 西山, 野呂, 小原, 清, 高山, 竹村, Holonomic Gradient Descent and its Application to Fisher-Bingham

Integral, Advances in Applied Mathematics 47 (2011), 639–658,

Q. ところで, 3.2 節の holonomic ideal (system) の中には常微分方
程式はどこにもないですが.
A. これらに微分作用素を左から掛けたりして, 足したり, 引いたりすると常微分方程式を作れます. d = 1 の時,

G=test1(); G[0];((-8*x_1_1^3+24*x_2_2*x_1_1^2+(-8*x_1_2^2-24*x_2_2^2)*x_1_1+8*x_2_2*x_1_2^2+8*x_2_2^3)*x_1+(-16*x_1_2*x_1_1^2+32*x_2_2*x_1_2*x_1_1-16*x_1_2^3-16*x_2_2^2*x_1_2)*x_2)*dx_2^5+((-8*x_1_2*x_1_1+8*x_2_2*x_1_2)*x_1^2+(4*x_1_1^2-8*x_2_2*x_1_1-20*x_1_2^2+4*x_2_2^2)*x_2*x_1+(16*x_1_2*x_1_1-16*x_2_2*x_1_2)*x_2^2+16*x_1_2*x_1_1^3-48*x_2_2*x_1_2*x_1_1^2+(16*x_1_2^3+48*x_2_2^2*x_1_2)*x_1_1-16*x_2_2*x_1_2^3-16*x_2_2^3*x_1_2)*dx_2^4+((-2*x_1_1+2*x_2_2)*x_1^3-8*x_1_2*x_2*x_1^2+((2*x_1_1-2*x_2_2)*x_2^2+8*x_1_1^3+(-24*x_2_2+24)*x_1_1^2+(24*x_1_2^2+24*x_2_2^2-48*x_2_2)*x_1_1-24*x_2_2*x_1_2^2-8*x_2_2^3+24*x_2_2^2)*x_1-4*x_1_2*x_2^3+(48*x_1_2*x_1_1+16*x_1_2^3-48*x_2_2*x_1_2)*x_2)*dx_2^3+(-x_2*x_1^3+(8*x_1_2*x_1_1+(-8*x_2_2+2)*x_1_2)*x_1^2+(-x_2^3+(-4*x_1_1^2+(8*x_2_2-6)*x_1_1+12*x_1_2^2-4*x_2_2^2+6*x_2_2)*x_2)*x_1+(-12*x_1_2*x_1_1+(12*x_2_2-20)*x_1_2)*x_2^2-16*x_1_2*x_1_1^3+(48*x_2_2-16)*x_1_2*x_1_1^2+(-16*x_1_2^3+(-48*x_2_22+32*x_2_2)*x_1_2)*x_1_1+(16*x_2_2+16)*x_1_2^3+(16*x_2_2^3-16*x_2_2^2)*x_1_2)*dx_2^2+(x_1^3+2*x_1_2*x_2*x_1^2+((-2*x_1_1+2*x_2_2-1)*x_2^2-20*x_1_1^2+(-10*x_1_2^2+40*x_2_2-6)*x_1_1+(10*x_2_2++(-x_1_2*x_1^2+(x_2^3+(6*x_1_1-x_1_2^2-6*x_2_2)*x_2)*x_1+(-4*x_1_2*x_1_1+(4*x_2_2+16)*x_1_2)*x_2^2+8*x_1_2+(-4*x_2_2-8)*x_1_2^3+8*x_2_2^2*x_1_2)*dx_2^0[1957] vars(G[0]);[x_1,x_2,x_1_1,x_1_2,x_2_2,dx_2]

Wishart 分布の第一固有値, 行列引数超幾何系 1 F1 (予稿 4. 節)
累積分布関数 Pr[ℓ1 < x], m = 10, n = 12,
Σ−1 = diag(1, 2, . . . , 10)

1.2
by hg

10.8

0.6

0.40.2

出典: 橋口, 沼田
0 , 高山, 竹村, Holonomic gradient method for the distribution function of the largest root of a051015

Wishart matrix, Journal of Multivariate Analysis, 117, (2013) 296-312.

20

. Orthant probability (予稿 5 節)

Orthant
Probability のHGM による計算(小山, 竹村,2012))

. 離散 A 超幾何分布 (予稿 7 節, arxiv:1510.02269)
差分 HGM では厳密計算が可能. 近似公式の精度を調べるなど理論研究にも使える.
.
Theorem (T-栗木-竹村 (arxiv:1510.02269)).
k → +∞ で,Z (kβ; x) ∼x km(2πk)n−dΓ(km + 1) (det(ĀM −1 ĀT ))1/2

M. = diag(m).
2 × 4 分割表. 行和 (4, 19)k, 列和 (9, 5, 3, 6)k,kβ = k(19, 9, 5, 3, 6). x = (1, 1/3, 12, 1/5001, 1, 1, 1, 1). x4 は他より小さい. xi の期待値 ( xi ∂i • Z /Z ) の近似 (by IPS) m =(2.79518, 0.652785, 0.551505, 0.000540425, 6.20482, 4.34722, 2.4485, 5.99946).

期待値の近似 (by IPS) m =(2.79518, 0.652785, 0.551505, 0.000540425, 6.20482, 4.34722, 2.4485, 5.99946).
HGM による厳密値:(2.83214, 0.627808, 0.539555, 0.000496547, 6.16786, 4.37219, 2.46044, 5.9995).
比:(0.98695, 1.03978, 1.02215, 1.08837, 1.00599, 0.994289, 0.995147, 0.999993).

k → +∞ の時
k9200300log Z−568.0127−26598.4556−42685.5415

Approx by Th 5
−569.8179−26598.9446−42685.9149|error|1.80520.48900.3734注: 近似値で確率を計算すると誤差は大きい. 例.
u = (33, 1, 1, 1, 48, 44, 26, 53), k = 9. u を得る正確な確率は3.26465 × 10−7 . 上記の log Z の近似による確率は1.98529 × 10−6 . 約６倍異なる.
17/18

厳密計算なのでこんな数が出力されます

N=200
A=[[1,0,0,1,0,1,0,1,0,1],[0,1,0,1,0,1,0,1,0,1],[0,0,1,-1,0,0,0,0,0,0],[0,0,0,0,1,-1,0,0,0,0],[0,0,0,0,0,0,1,-1,0,0],[0,0,0,0,0,0,0,0,1,-1]]
Beta=[452,412,-37,-23,-13,31]
at ([x1,x2,x3,x4,x5,x6,x7,x8,x9,x10]=[140/411,40/137,25/822,31/411,14/411,17/274,17/822,5/137,10/137,29/82
oohg_native=0, oohg_curl=1
x3 の期待値を計算.
EV([x3])=[484018240471728953822203320553380653219481012643866487201043272204554116427335942534923953734369
86365699839168924385947529623435213755551773022215922104722152504652845614751116627622765024345097422807743057500921935232293131676851615762862014663994664872134693815356637343841938809747418295142613240962333345
3442753508220352031310549167268194351651787783253898660000276995488979059934881671963927282777353837308851/1944222849842515553043842429125888595116006553330637894368400560720768008344952556960403129403576682658482063685905755102313943954044436017805458085864176093731784381898126374058702803535631819651190493876403501
9417725144895331947497817468402087056746060088760317342886715324762007018565160119564515972685383799358741
32090627201429825951569856280808639609886906110220425511570638764915578591464428000430220868340937739443549573932056327206030262721912023810463723569352286063413912998077871191506911]
Time=84562.4

N
級数の直接計算
Macaulay 型行列の方法
06822s (1.89 時間)
61399s (約 17 時間)
Intel Xeon
100 138640s (1 日と約 14.5 時間) 73126s(約 20.3 時間)
200 未計測84562s (約 23.5 時間)
E5-4650 (2.7GHz) with 256G memory, the computer algebra system Risa/Asir(20140528). パッケージ ot hgm ahg.rr.