Introduction and Advances in Gaussian Processes

持橋大地
NTT コミュニケーション科学基礎研究所
daichi@cslab.kecl.ntt.co.jp
SVM 2009
Modified for NTT 厚木基礎研
2009–10–13 (Tue)Introduction and Advances in Gaussian Processes – p.1/29

流れ• Gaussian Process とは?
◦ SVM との関係
◦ Dirichlet Process との関係
◦ 例と簡単な応用• Gaussian Process Latent Variable Models (Lawrence 2003)• Gaussian Process Dynamical Models (Hertzmann+ 2005)• Gaussian Process Density Sampler (Adams+ ICML 2008)• Discussion• Literature and Readings

Introduction and Advances in Gaussian Processes – p.2/29

Gaussian Process とは210−1−2−3−8−6−4−202468• 入力 x → y を予測する回帰関数 (regressor) の確率モデル◦N(n)(n)データ D = (x , y ) n=1 が与えられた時, 新しいx(n+1) に対する y (n+1) を予測◦ ランダムな関数の確率分布◦ 連続空間で動く, ベイズ的なカーネルマシン (後で)Introduction and Advances in Gaussian Processes – p.3/29

Gaussian Process とは3210−1−2−3−4−8−6−4−202468• 入力 x → y を予測する回帰関数 (regressor) の確率モデル◦N(n)(n)データ D = (x , y ) n=1 が与えられた時, 新しいx(n+1) に対する y (n+1) を予測◦ ランダムな関数の確率分布◦ 連続空間で動く, ベイズ的なカーネルマシン (後で)Introduction and Advances in Gaussian Processes – p.3/29

Gaussian Process とは3210−1−2−3−4−8−6−4−202468• 入力 x → y を予測する回帰関数 (regressor) の確率モデル◦N(n)(n)データ D = (x , y ) n=1 が与えられた時, 新しいx(n+1) に対する y (n+1) を予測◦ ランダムな関数の確率分布◦ 連続空間で動く, ベイズ的なカーネルマシン (後で)Introduction and Advances in Gaussian Processes – p.3/29

線形回帰モデル• 入力 x から出力 y ∈ R を予測する回帰関数 y = f (x) を求めたい◦ x = (x1 , · · · , xd ) ∈ Rd は時間や任意のベクトル• y = f (x) を, x を一般の関数 φ(x) で変換した上で線形モデルで表してみる◦ 例:

y = wT φ(x)(1)

φ(x) = (1, x1 , · · · , xd , x21 , · · · , x2d )w = (w0 , w1 , · · · , wd , wd+1 , · · · , w2d ) とおくと,
y = wT φ(x)= w 0 + w1 x 1 + · · · + wd x d+ wd+1 x21 + · · · + w2d x2d .

◦ 一般に, 基底関数 φ(x) = (φ1 (x), · · · , φH (x))
Introduction and Advances in Gaussian Processes – p.4/29

GP の導入 (1)• y (1) · · · y (N ) について同時に書くと, 下のように y = Φw と行列形式で書ける (Φ : 計画行列)⎛y (1)⎞⎛

φ1 (x(1) )⎜⎜ φ1 (x(2) )······(x(1) )⎞⎛⎞w1φH⎟⎜ ⎟⎜ (2) ⎟φH (x(2) ) ⎟ ⎜ w2 ⎟⎜y ⎟⎟⎜ . ⎟⎜ . ⎟=⎜..
..
⎟⎜ . ⎟⎜ . ⎟ ⎜.
.
⎠⎜ . ⎟⎝ . ⎠ ⎝
⎜ . ⎟(N)(N)(N)⎟yφ1 (x ) · · · φH (x ) ⎜⎝ .. ⎠wHy

Φ(2)w• 重み w がガウス分布 p(w) = N(0, α−1 I) に従っているとすると, y = Φw もガウス分布に従い,• 平均 0, 分散

 yyT  = (Φw) (Φw)T= α−1 ΦΦT= ΦwwT ΦT(3)の正規分布となる
Introduction and Advances in Gaussian Processes – p.5/29

GP の導入 (2)p(y) = N(y|0, α−1 ΦΦT )(4)は, どんな入力 {xn }N
n=1 についても成り立つ → ガウス過程の定義• どんな入力 (x1 , x2 , · · · , xN ) についても, 対応する出力

y = (y1 , y2 , · · · , yN ) がガウス分布に従うとき, p(y) はガウス過程に従う という.
◦ ガウス過程 = 無限次元のガウス分布◦ ガウス分布の周辺化はまたガウス分布なので, 実際にはデータのある所だけの有限次元• K = α−1 ΦΦT の要素であるカーネル関数k(x, x ) = α−1 φ(x)T φ(x )(5)だけでガウス分布が定まる◦ k(x, x ) は x と x の距離 ; 入力 x が近い→出力 y が近い
Introduction and Advances in Gaussian Processes – p.6/29

GP の導入 (3)• 実際には, 観測値にはノイズ  が乗っている

y = wT φ(x) +  ∼ N(0, β −1 I)=⇒p(y|f ) = N(wT φ(x), β −1 I)(6)• 途中の f = wT φ(x) を積分消去

p(y|x) =p(y|f )p(f |x)df= N(0, C)(7)(8)◦ 二つの独立な Gaussian の畳み込みなので, C の要素は共分散の単純な和:C(xi , xj ) = k(xi , xj ) + β −1 δ(i, j).

(9)◦ GP は, カーネル関数 k(x, x ) とハイパーパラメータ α, β だけで表すことができる.
Introduction and Advances in Gaussian Processes – p.7/29

様々なカーネル
2.5

322

1.5
1

10yy

0.5
0

−0.5
−1−1−1.5−2−2−2.5−5−4−3−2−1012

34−3−55−4−3−2−1x0

12

345x

Gaussian: exp(−(x − x )2 /l)Exponential: exp(−|x − x |/l) (OU process)2.5

3221.511yy0.500−0.5−1−1−1.5−2−2−3−5−4−3−2−10x123452 x−x
Periodic: exp(−2 sin ( 2 )/l2 )−2.5−5−4−3−2−1012345x2 x−x
Periodic(L): exp(−2 sin ( 2 )/(10l)2 )Introduction and Advances in Gaussian Processes – p.8/29

「基底関数」の消去• RBF 基底関数 φ(x) = exp((x − h)2 /r2 ) を考えてみる• 1 次元の場合, h を無限個用意すると

k(x, x ) = σ 2Hh=1∞

φh (x)φh (x )(10)2(x − h)exp −dh→exp −22rr−∞(11)

 2
 2
√
(x − x )(x − x )2≡ θ1 exp −= πr exp −22r
θ22
(x − h)2(12)◦ (x, x ) の RBF カーネルは, 無限個の RBF 基底関数を考えたことと等価.
◦ カーネルのパラメータ θ1 , θ2 は最尤推定で最適化できる
Introduction and Advances in Gaussian Processes – p.9/29

GP の予測• 新しい入力 y new とこれまでの y の結合分布がまた Gaussianになるので,

p(y new |xnew , X, y, θ)(13)p((y, y new )|(X, xnew ), θ)=(14)p(y|X, θ)⎛⎞−1 
y1K knewT −1 ⎠⎝−y∝ exp − ([y, y ]K y)
Tnew2ykk(15)∼ N(kT K−1 y, k − kT K−1 k).

(16)

ここで
◦ K = [k(x, x )], k = (k(xnew , x1 ), · · · , k(xnew , xN )),◦ k = k(xnew , xnew ) + β −1 .
Introduction and Advances in Gaussian Processes – p.10/29

GP の適用例• 『ガウシアンプロセスによる名演奏の学習』寺村&前田+,

IPSJ-MUS 2008 & ICMPC 2008.
◦ 楽譜情報 x →実際の演奏 y の対応を学習◦ 実際には,− x: 楽譜からの素性 (拍子, 音高, 音長, 強弱指定, 主旋律か, …)− y : 楽譜の指定と実際の出力との差 (アタック, リリース,強弱の 3 種類の y )◦ {x(n) , y (n) }N
n=1 のセットから対応関係を学習• 学習のポイント: (x(1) , · · · , x(N ) ) → (y (1) , · · · , y (N ) ) への対応を「まとめて」学習していること◦ カーネルはak(x, x ) =(1 + (x − x )T B(x − x ))c◦ ただし, (x(1) , · · · , x(N ) ) が時系列という仮定は無いのに注意
Introduction and Advances in Gaussian Processes – p.11/29

GP による分類と SVM• y = {+1, −1} のとき, p(y|f ) = σ(y · f ) (logit) または Ψ(y · f )(probit) で分類器になるminimize: − log p(y|f )p(f |X)
1 T −1= f K f−log p(yi |fi )2N(17)

i=1• ソフトマージン SVM では Kα = f とすると,

w=

αi xi → |w|2 = αT Kα = f T K −1 f ゆえ,iN1minimize: |w|2 − C(1 − yi fi )+2i=1N
1= f T K −1 f − C(1 − yi fi )+ .
2(18)

i=1

◦ そっくりだが, SVM は hinge loss なのでスパース.
Introduction and Advances in Gaussian Processes – p.12/29

Loss functions
Relationships between GPs and Other Modelslog(1 + exp(−z))−log Φ(z)max(1−z, 0)2g (z)10−2014.

−(a)0 −(b)zFigure 6.3: (a) A comparison of the hinge error, gλ and gΦ . (b) The -insensitive
error function used in SVR.

• お馴染みの議論かも◦ SVM と ME の関係についても, 以前工藤君がそっくりの議論をしていた• CS 研での議論: 分類したいだけなら, GP classifier は回りく

どすぎる (他の理由があるなら有効)
Introduction and Advances in Gaussian Processes – p.13/29

DP との関係• Gaussian process と Dirichlet process の定義はそっくり[偶然ではない]◦ GP: どんな (x1 , x2 , · · · , x∞ ) をとってきても, 対応する(y1 , y2 , · · · , y∞ ) がガウス分布に従う◦ DP: 空間のどんな離散化 (X1 , X2 , · · · , X∞ ) についても, 対応する離散分布がディリクレ分布
Dir(α(X1 ), α(X2 ), · · · , α(X∞ )) に従う• どちらも, 無限次元の smoother になっている
Introduction and Advances in Gaussian Processes – p.14/29

Gaussian Process Latent Variable Models (GPLVM)• Probabilistic PCA (Tipping&Bishop 1999):p(yn |W, β) =p(yn |xn ,W, β)p(xn )dxn (19)
p(Y|W, β) =p(yn |W, β) → W を最適化n• GPLVM (Lawrence, NIPS 2003): W の方に prior を与えて積分消去

p(W) =DN(wi |0, α−1 I)(20)

i=1p(Y|X, β) =p(Y|X, β)p(W)dW(21)

11−1Ttr(K=exp−YY)
(22)
DN/2
D/2
2(2π)
|K|
Introduction and Advances in Gaussian Processes – p.15/29

GPLVM (2): PPCA の DualD
1DN
log(2π) − log |K| − tr(K−1 YYT )log p(Y|X, β) = −222
K = αXX + βT−1(23)(24)

I

X = [x1 , · · · , xN ]T(25)• X に関して微分すると,∂L= αK−1 YYT K−1 X − αDK−1 X = 0∂X1⇐⇒ X = YYT K−1 X ⇒ X UQ LV TD(26)(27)◦ UQ (N × Q) : YY T の Q 個の最大固有値 λ1 · · · λQ に対応する固有ベクトル◦ L = diag(l1 , · · · , lQ ); li = 1/ λi − 1αDαβ
Introduction and Advances in Gaussian Processes – p.16/29

GPLVM (3) : Kernel 化
D1DN
log(2π) − log |K| − tr(K−1 YYT )log p(Y|X, β) = −222(28)K = αXXT + β −1 I ,X = [x1 , · · · , xN ]T(29)• 自然にカーネル化されている =⇒ 任意のカーネル K を導入 γ
k(xn , xm ) = α exp − (xn − xm )2 + δ(n, m)β −12

•

(30)

∂L= K−1 YYT K−1 − DK−1
∂K
∂L ∂K
∂L
◦
=を適用して微分∂xn,j∂K ∂xn,j◦ Scaled Conjugate Gradient で解ける

GPLVM in MATLAB: http://www.cs.man.ac.uk/˜neill/gplvm/
Introduction and Advances in Gaussian Processes – p.17/29

GPLVM (4): 例
0.3
2.4
0.2

0.2

2.2
0.15
20.1
0.1

1.8
0.05

01.6
0−0.1

1.4−0.05
1.2

−0.2
−0.1

1
−0.15

−0.3

0.8
−0.2

−0.4

0.6
−0.25

−0.5
−0.8

−0.6

−0.4

−0.2

00.2

0.4

0.6

−0.2

−0.1

00.1

0.2

0.3

0.4

Figure 1: Visualisation of the Oil data with (a) PCA (a linear GPLVM) and (b) A GPLVM whichuses an RBF kernel. Crosses, circles and plus signs represent stratified, annular and homogeneousflows respectively. The greyscales in plot (b) indicate the precision with which the manifold wasexpressed in data-space for that latent point. The optimised parameters of the kernel were,and f.

    

• 線形の PPCA(左) より, GP-LVM(右) の方が分離性能が高い◦ ベイズなので, Confidence の分布が同時に得られる• 計算量が問題 (O(N 3 )): active set (サポートベクターみたいなもの) を選んで, データをスパース化して最適化
Introduction and Advances in Gaussian Processes – p.18/29

Gaussian Process Dynamical Model (Hertzmann 2005)

http://www.dgp.toronto.edu/˜jmwang/gpdm/• GPLVM では, 潜在変数 xn に分布がなかった⇓• xn が (GP で) 時間発展するモデル.
◦ 人間の動き (角度ベクトル) 等の時系列データ.
◦ 自然言語の時系列データ?
Introduction and Advances in Gaussian Processes – p.19/29

GPDM (2): Formulation (1)xt = f (xt−1 ; A) + x,t,yt = g(xt ; B) + y,tf ∼ GP(0, Kx )(31)

g ∼ GP(0, Ky )(32)

• 結合確率 p(Y, X|α, β) = p(Y|X, β)p(X|α) を考える.
• 第1項

1−12 Ttr(Kp(Y|X, β) =exp
−YWY )
Y
N
D/2
D/2
2(2π)
|KY ||W|N

(33)

は GPLVM と基本的に同じ.
◦ KY は (この論文では) 普通の RBF カーネル
Introduction and Advances in Gaussian Processes – p.20/29

GPDM (3): Formulation (2)• 第 2 項は Markov 時系列N
p(xt |xt−1 , A, α) p(A|α) dA(34)  t=2Gaussian
11−1Ttr(K= p(x1 )exp−XX−−)
Xd(N−1)/2d2(2π)|KX |p(X|α) = p(x1 )(35)◦ X− = [x2 , · · · , xN ]T とおいた
◦ KX は x1 · · · xN −1 の RBF+線形カーネル

 α
2k(x, x ) = α1 exp − ||x − x ||2 +α3 xT x+α4−1 δ(x, x ) .
2(36)Introduction and Advances in Gaussian Processes – p.21/29

GPDM (4): Formulation(3)p(Y, X, α, β) = p(Y|X, β)p(X|α)p(α)p(β)
−1p(α) ∝αi , p(β) ∝βi−1 .
i(37)(38)i• 対数尤度は1
1−1
T
2 T− log p(Y, X, α, β) = tr(KX X− X− ) + tr(K−1YWY )Y
22d
D+ log |KX | + log |KY | (正則化項)
22
− log |W| +log αj +log βj (39)
j



j



(定数)

−→ 最小化.

(40)
Introduction and Advances in Gaussian Processes – p.22/29

Gaussian Process Density Sampler (1)3333222211110000−1−1−1−1−2−2−2−2−3−3−2−1012(a) x = 1, y = 1, α = 1

3−3−3−2−1012(b) x = 1, y = 1, α = 10

3−3−3−2−10123(c) x = 0.2, y = 0.2, α = 5−3−3−2−10123(d) x = 0.1, y = 2, α = 5• GP は任意の関数の prior →確率密度関数のモデルに使えないか?

1Φ(f (x))π(x)p(x) =Z(f )(41)◦ f (x) ∼ GP(x) ; π(x) : 事前分布◦ Φ(x) ∈ [0, 1] : シグモイド関数− ex. Φ(x) = 1/(1 + exp(−x))
Introduction and Advances in Gaussian Processes – p.23/29

Gaussian Process Density Sampler (2)1Φ(f (x))π(x)p(x) =Z(f )(42)• 生成プロセス: Rejection sampling

1. Draw x ∼ π(x).
2. Draw r ∼ Uniform[0, 1].
3. If r < Φ(g(x)) then accept x; else reject x• Accept された N 個の観測データの背後に, reject された M 個のデータとその場所が存在 (隠れ変数)◦ 分配関数 Z(f ) は求まらないが, Φ(g(x)) は求まる
→ MCMC!
◦ Infinite Mixture とは別の確率密度のモデル化

Introduction and Advances in Gaussian Processes – p.24/29

Gaussian Processes for NLP?
• 言語では出力 y = (y1 , · · · , yN ) は Gaussian ではないので, 直接は使えない• 隠れ変数の (連続的な) 分布に使える可能性は高い◦ Ex: Latent CCA in (Haghighi+, ACL 2008)• 分類もできるが, 特に GP でないといけない理由は希少• 言語+α のモデル化に特に重要だと思われる

Introduction and Advances in Gaussian Processes – p.25/29

まとめ• Gaussian Process…連続的な関数のベイズモデル◦ カーネル関数で定義される, 無限次元のガウス分布◦ 基底関数の空間での線形モデルで, 重みを積分消去したもの• カーネル設計が重要◦ カーネルの違いにより, さまざまな振舞い◦ カーネルのパラメータは, 確率モデルなので最適化できる• 回帰問題だけでなく, 隠れ変数や確率分布のモデルにも使える,確率モデルの Building Block• 計算量が課題→研究は最近様々に進んでいる

Introduction and Advances in Gaussian Processes – p.26/29

Literature• “Gaussian Process Dynamical Models”. J. Wang, D. Fleet,and A. Hertzmann. NIPS 2005.
http://www.dgp.toronto.edu/ jmwang/gpdm/• “Gaussian Process Latent Variable Models for Visualizationof High Dimensional Data”. Neil D. Lawrence, NIPS 2003.
• “The Gaussian Process Density Sampler”. Ryan Prescott

Adams, Iain Murray and David MacKay. NIPS 2008.
• “Archipelago: Nonparametric Bayesian Semi-SupervisedLearning”. Ryan Prescott Adams and Zoubin Ghahramani.
ICML 2009.

Introduction and Advances in Gaussian Processes – p.27/29

Readings• “Gaussian Processes for Machine Learning”. Rasmussenand Williams, MIT Press, 2006.
http://www.gaussianprocess.org/gpml/• “Gaussian Processes — A Replacement for Supervised

Neural Networks?”. David MacKay, Lecture notes at NIPS
1997. http://www.inference.phy.cam.ac.uk/mackay/GP/◦ Videolectures.net: “Gaussian Process Basics.”
http://videolectures.net/gpip06_mackay_gpb/• Pattern Recognition and Machine Learning, Chapter 6.

Christopher Bishop, Springer, 2006.
http://ibisforest.org/index.php?PRML• ガウス過程に関するメモ (1). 正田備也, 2007.

http://www.iris.dti.ne.jp/˜tmasada/2007071101.pdf

Introduction and Advances in Gaussian Processes – p.28/29

Codes• GPML MATLAB codes.

http://www.gaussianprocess.org/gpml/code/• GPLVM MATLAB/C++ codes.

http://www.cs.manchester.ac.uk/˜neill/gplvm/• GPDM code. http://www.dgp.toronto.edu/˜jmwang/gpdm/

Introduction and Advances in Gaussian Processes – p.29/29