機械学習プロフェッショナルシリーズ – ガウス過程と機械学習読書メモ

機械学習技術人工知能技術デジタルトランスフォーメーション技術確率的生成モデル本ブログのナビ深層学習技術ノンパラメトリックベイズとガウス過程

サマリー

ガウス過程 (Gaussian Process; GP) は、確率論に基づくノンパラメトリックな回帰やクラス分類の手法であり、連続的なデータのモデリングに利用される確率過程の一種となる。ガウス過程のアプローチでは、データの生成プロセスを確率分布でモデル化し、データ間の関係性を表現するためのカーネル (カーネル関数、カーネル関数行列) を用いて、データ点間の相関を表現し、データ点の間の関係性や不確実性を推定する。ガウス過程のカーネルは、様々な形状を持つことができ、データの特性に合わせてカスタマイズすることができるため、データに適合するような柔軟なモデルを構築することができる。ここではこのガウス過程について機械学習プロフェッショナルシリーズ「ガウス過程と機械学習」をベースに述べる。

今回は読書メモを記載する

機械学習プロフェッショナルシリーズ – ガウス過程と機械学習

「圧倒的に柔軟なベイズ的回帰モデルであるガウス過程の日本初の入門書。基礎の線形回帰から始め、ガウス過程の原理をゼロからていねいに解説。教師なし学習、実応用など最近の話題まで紹介した。さあ、はじめよう！」

第０章たった５分でガウス過程法が分かってしまう

概要
ガウス過程法は、ガウス過程の理論を用いてパターン認識を行う手法

0.1 第１ステップ：機械学習って何？
機械学習とは、人工的に作られた関数関係

0.2 第２ステップ：回帰と最小二乗法
パラメータ空間とデータ空間

0.3 第３ステップ：確率モデリングとベイズ推定
ベイズ推定でのパラメータ空間の事後分布
事前分布として何もわからない状態だと、雲が広がりデータもばらつく
事後分布として実データの影響を反映するとデータは収束する
パラメータが収束

0.4 第４ステップ：ガウス分布と共分散
一次元ガウス分布と二次元ガウス分布
2つ以上の変数がガウス分布に従うとき、これらをまとめたベクトルの分布を多変量ガウス分布
多変量ガウス分布は中心位置と共分散行列で定義される
2つの確率変数の互いに独立でないと共分散がゼロで無くなる

0.5 第５ステップ：ガウス過程とガウス過程回帰
機械学習とは関数f(x)を学習で求めること
データに合うようにf(x)を求めることが回帰問題の目的
回帰ではf(x)の範囲、つまり関数モデルを適切に決めることが重要
ガウス過程回帰で求めたいものは関数f(x)の事後分布の雲
ベイズ推定で求める
関数f(x)の任意個数の入力点x1,…,xnに対して、その出力f(x1),…,f(xN)の分布がN次元の多変量ガウス分布
概念図
ガウス過程回帰の例
コラム：関数の雲とガウス過程
関数の雲が得られると
わからなさの程度がわかる
自信のある領域とない領域の違いがわかる
モデル選択や特徴選択ができる

第１章線形回帰モデル

概要
ガウス過程は線形回帰モデルの無次元への拡張
ある入力xが与えられた時、対応する出力yの値を予測することを回帰(regression)という
連続値が回帰、離散値の場合は識別あるいは分類(classification)

1.1 単回帰

実数xから実数yへの回帰
N個の(x,y)の観測値のペア
D={(x1,y1), (x2,y2), …, (xN,yN)}
最も単純な回帰モデル(単回帰)
y(xから予測されたy値)=a +bx (a,b ∈R)
式1
式2
最小になるのはa,bについての微分が0になる時
1
2
公式単回帰の正規方程式
(例)D={(3,2),(2,4),(-1,1)}の時
1
2

1.2 重回帰とベクトル表現

入力xが一次元でなく、D次元のベクトルxT = (x1, x2, … ,xD)
転置T
縦ベクトルとしての表現
y(xTから予測されたy)=𝛚0 + 𝛚1×1 + 𝛚2×2 +…+ 𝛚DxD
重回帰(multiple regression)
2次元の例
y=1 + 1.2×1 + 0.6×2
(yn – yn(XTから予測されたy)) = (yn (𝛚0 + 𝛚1×1 + 𝛚2×2 +…+ 𝛚DxD))2
y(xTから予測されたy)=𝛚0 + 𝛚1×1 + 𝛚2×2 +…+ 𝛚DxD
N個のデータに対しては
式
Xを計画行列(design matrix)と呼ぶ
ベクトルの一次式の微分の公式
ベクトルの二次形式の微分
(例)

1.3 線形回帰モデル

xの関数𝛟1(x),𝛟2(x),…,𝛟H(x)をH個準備
y(xから予測されたyの値)= 𝛚0 + 𝛚1𝛟1(x) + 𝛚2𝛟2(x) +….+ 𝛚H𝛟H(x)
H=2で𝛟1(x)=x, 𝛟2(x)=x2
y(xから予測されたyの値) = 𝛟𝙬
例
基底ベクトルと特徴ベクトル
𝜱(X) = (1, x, x2, six)
y(xから予測されたy)=WT𝜱(X)
パラメータwに関しては一次式(線形):線形回帰モデル(linaer regression model)
= 𝜔0𝜙0(x) + 𝜔1𝜙1(x) +…+ 𝜔H𝜙H(x)
基底関数(base function)
𝜙0(x), 𝜙1(x), …, 𝜙H(x)
特徴ベクトル(feature vector)
𝜙(X) = (𝜙0(x), 𝜙1(x), …, 𝜙H(x))T
線形回帰モデルの解
例
𝜙(x)=(1, x, x2, six, cosx)T
𝜙(x)が”tokyo”の場合、基底関数𝜙(x)として「xの文字列長」「xが大文字で始まるか(始まれば1、羽島せなければ0)」「xが動詞か」のような関数を考えると𝜙(x)=(5, 1, 0)Tのようになる

1.4＊リッジ回帰

解をうまく計算できないケースがある
逆行列(𝜱T𝜱)-1が存在しない
計算が不安定になる
𝜱(x)のあるi番目の要素が、データ点全てで他のj番目の要素と同じ、またはその定数倍
𝜱の各列は線型独立ではなくなる
例
y = WTX = 𝜔0 + 𝜔1×1 + 𝜔2×2
X1T = (2, 4), X2T = (3, 6), X3T = (4, 8)
行列表現
X
XTX
逆行列は存在しない
完全に定数倍でなくてもほとんど同じ場合
逆行列
𝜔の絶対値が大きくなって不安定
係数Wiの絶対値が極端に大きくて好ましくないとき
正則化(regularization)
係数ベクトルの最適化を行うときに、係数が極端な値を取らないための制約を与える
|W|2 = WTW
E = (y-Xw)2 + αWTW
リッジ回帰の解
コラム：相関係数と回帰モデル
HSIC(Hilbert-Schmidt Independence Criterrion)
非線形な相関を高精度に捉える
カーネル法を使えば実数値やベクトルだけでなく、文字列やグラフ、木といった一般的なデータ構造についてもカーネルが定義できるため、それらの間の相関を捉えることができる
カーネルはガウス過程でも使える

第２章ガウス分布

2.1 ガウス分布とは

概要
回帰モデルの「当てはまりの良さ」は二乗誤差を基準として天下り的に定義
線形回帰モデルを確率モデルとして考える
確率モデルとして考えるということは、データが生成される確率を与えること
xが与えられた時にyを予測するのが回帰モデル
Yがとる確率は?
実際のデータは直線や平面上に全て乗るのではなく、それから少しずれる
ずれεが、誤差の分布として最も標準的な平均0,分散σ2のガウス分布 (Gaussian distribution) p(ε)=N(0,σ2)に従う場合
ガウス分布の確率密度
平均μ、分散σ2
形状
先頭の定数は、xについての積分を1にして確率密度にするための正規化定数
確率密度関数の本体は exp(・・)
比例を表す記号∝を用いて、定数を除いて
平均0、分散1の正規分布を特に、標準正規分布(standard normal distribution)
誤差の分布のかたちが不明な場合、誤差分布のモデルとして正規分布が広く用いられる
必ずしも正規分布に従わない互いに独立な誤差でも
その挿話は中心極限定理により正規分布に近づく
誤差の分布形状が正規分布に従わないことがわかっていても、多くの場合はそう悪い近似にはならない
平均と分散を定めた時、正規分布は情報量最小化(=エントロピー最大化)である
標本サイズが小さい場合でも安定して良い近似が得られる
特別な理由がない限り「正規分布」は標準として扱うべきモデル
誤差分布にコーシー分布を仮定して、外れ値に対して頑健な回帰モデルを作る工夫もある

2.1.1 ガウス分布からのサンプリング

標準正規分布N(0,1)からのサンプルxは (0,1]からの一様乱数を発生させる関数rand()があれば生成することができる(Box-Muller法)
公式
サンプル
平均𝜇、分散𝜎2のガウス分布N(𝜇, 𝜎2)では、 xのスケールが𝜎(標準偏差)倍で、値が𝜇倍だけシフト
X’ = 𝜇 + 𝜎x

2.1.2 線形回帰の確率モデル

xが与えられた時のyの平均は線形回帰
WTX
実際のyの値はこれにノイズεがのったもの
y=WTX + ε
p(ε) = N(0,σ2)
線形回帰モデルへ適用すると
p(y) = N(WTx, 𝜎2)
N個のデータ点のペアDがあった時
Xとyについて
DのN個のペア(xn, yn)が互いに独立であるとすると、Xからyを予測する確率は
log xは単調増加関数なので、 p(y|X)の最大化とlog p(y|X)の最大化は等価
上式の対数を取れば
二乗誤差による線形回帰と同じ

2.2 重みの事前分布とリッジ回帰

回帰係数wの成分として極端な値を避けるリッジ回帰も確率モデルから導くことができる
Wの各要素wiの絶対値が極端に大きくなるのを防ぐ
Wiがそれぞれ独立に平均0,分散λ2の正規分布に従う確率変数
リッジ回帰の式
αの値は学習データのフィッティングについて、観察ノイズと係数のどちらを重視するかのバランスを決める意味を持つ

2.3 多変量ガウス分布

概要
多変量のガウス分布について
多変量ガウス分布の確率密度関数
D次元のベクトルX=(X1,..,Xc)が平均μ、共分散行列Σのガウス分布N(x|μ,Σ)にしたがっている時
μ=(μ1,…,μD)はXの期待値を表す平均ベクトル
μ=𝔼[x] ΣはDxDの共分散行列 (covariance matrix)
その(i,j)要素がxiとxjの共分散を表す
Σij = 𝔼[(xi-𝔼[xi])(xj-𝔼(xj)] = 𝔼[xi,xj]- 𝔼[xi]𝔼[xj] Σの要素である全てのxiと様の組み合わせ
共分散行列はまとめて Σ=𝔼[xxT] – 𝔼[x]𝔼[x]T と行列計算で表せる
ガウス分布では共分散行列Σは常にΣ-1の形で現れるため
𝛬 = Σ-1を精度行列 (precision matrix)と呼ぶ
(μ,Σ)ではなく(μ, 𝞚)とすることもある
Σは(x-μ)の共分散を表す
行列式|Σ-1|はΣ-1の春空間の大きさを表す
先頭の正規化定数
μは分布を平行移動するだけなので、ガウス分布の形状はΣで決まる
関数形例
対応式例
多変量ガウス分布の共分散行列
Σが対角行列の場合
任意の2つの次元が互いに独立で、多変量ガウス分布の確率密度は一次元ガウス分布の積となる
Σが対角行列でない場合
xの分布はその成分の間で相関を持つ
例1 正の相関
N(0,Σ)から100個の点を発生させる
x1とx2に強い相関がある
例2 負の相関
N(0,Σ)から100個の点を発生させる
例3 相関がない場合

2.3.1 多変量ガウス分布と線形変換

平均0のガウス分布に従うベクトル x 〜 N(0,Σ) ∝exp(-1/2xTΣ-1x)
を行列Aで線形変換してy=Axとしても、 yはやはりガウス分布に従う

2.3.2 多変量ガウス分布からのサンプリング

平均0、共分散Σの多変量ガウス分布からのサンプル y 〜 N(0,Σ)は
MATLABやOctaveでは”mvnrnd”で表せる
RではMASSパッケージのmvrnormで生成
自分で生成する場合にはまずΣを Σ=LLT のようにコレスキー分解する
コレスキー分解
対称行列Σを下三角行列、すなわち対角成分だけ値が埋まっている行列Lを用いて分解する
NATLAB,OctaveやRでは”cool”
N(0,Σ)に従うベクトルyを生成するには
まずx=(x1,x2,…xn)をランダムに生成
y=Lxと変換すれば良い

2.3.3 多変量ガウス分布の周辺化

ベクトルxが多変量ガウス分布に従っている時、 xの一部の次元を周辺化しても残りの次元は多変量ガウス分布に従う
例 2次元の多変量ガウス分布を1次元に周辺化
D次元のベクトルxが多変量ガウス分布x〜N(μ,Σ)にしたがっている時
xの最初のL次元をx1,残りのD-L次元をx2とおくと
μとΣも同様にL次元で分割すると
p(x)=p(x1,x2)をx2に関して周辺化したx1の分布は
公式3.2 多変量ガウス分布の周辺化
p(x) = ∫ p(x1, x2) dx2 = N(μ1, Σ11)
x2に対して周辺化することは、平均μ2や共分散行列の一部Σ12,Σ21, Σ22を見なかったことにすること

2.3.4 多変量ガウス分布の条件付き分布

ベクトルxが多変量ガウス分布に従う時、 xの一部の次元を固定した時の分布、すなわちある次元で「切った」切り口の分布も、やはりガウス分布になる
例 2次元のガウス分布p(x1,x2)において x1=c(定数)で切った時のx2の分布p(x2|x1)
Cによって異なる平均、分散のガウス分布となる
D次元のガウス分布に従うベクトルx 〜 N(μ, Σ)を x1,とx2に分けて考える
x1を固定した時のx2の分布 (同時分布p(x1,x2)をx1で切った時のx2の上付き分布)
公式2.4 多変量ガウス分布の条件付き確率
p(x2|x1) =N( μ2 + Σ21Σ11-1(x1 – μ1), Σ22 – Σ21Σ11-1Σ12)
公式の意味
条件付き分布の平均 : μ2 + Σ21Σ11-1(x1 – μ1)
観測値x1が与えられた時、それと相関を持つx2の新しい平均
条件付き分布の平均は元の平均μ2に、 x1のμ1からのずれ(x1-μ1)を行列Σ21Σ11-1で変換して加えたもの
x1の精度行列Σ11-1は
x1の分散が小さい、すなわち観測が正確な場合に大きくなる
x1の分散が大きい、すなわち観測があやふやな場合に小さくなる
Σ11-1が小さい(観測制度が高いほど)、「ずれ」x1-μ1は大きく拡大される
さらにx1とx2の共分散行列Σ21によってx2の空間に投射される
分散が大きいほど(観測の精度が低いほど)、x1はわかっても情報は少なく、x2が更新される量は減る
条件付きの分散 : Σ22 – Σ21Σ11-1Σ12
観測値x1が得られて、新しい分散は元の分散Σ22から Σ21Σ11-1Σ12だけ減少する
x1の精度行列Σ11-1が高い時は大きな値が、精度が低い時には小さな値が引かれる

第３章ガウス過程

概要
ガウス過程は線形回帰モデルの重みを積分消去したもの
無限次元のガウス分布
ガウス過程はランダムな関数を生成する確率分布

3.1 線形回帰モデルと次元の呪い

基底関数とその重み付けによる関数の表現
問題
入力xの次元が小さい場合でしか使えない
基底関数を配置すると次元が大きくなると指数的に必要な基底関数が増加する
次元ののろい(curse of dimensionlity)

3.2 ガウス過程

概要
xが高次元の場合でも柔軟な回帰モデルを実現するにはどうすれば良いか?
線形回帰モデルのパラメータWについて期待値をとって、モデルから積分消去する
入力xについて、xの特徴ベクトルをΦ(x)=(Φ0(x),Φ1(x),Φ2(x),…,ΦH(x))T
各特徴の重みW=(W0,W1,…,WH)による線形回帰モデル
ŷ = WTΦ(x) = w0φ0(x) + w1φ1(x) +….+wHφH(x)
行列式
ŷ = ΦW
Φ:計画行列
X1,…,XNが与えられれば定数行列
W:重みベクトル
ガウス分布に従う
Y 〜 N(0, λ2ΦΦT)
Wは期待値がとられて消えている
xやΦ(x)の次元がどんなに高くても、高次元の重みを求める必要はなく、 yの分布はデータ数Nに依存する共分散行列λ2ΦΦTで決まる
Ŷ = (ŷ1,…,ŷN)は観測データy = (y1,…,yN)に近づけるあてはめ値
Wを平均が0で分散がλ2Iのガウス分布から生成されたものとする
W 〜 N(0, λ2I)
定義3.1 ガウス過程の定義
どんなN個の入力の集合(x1,x2,…xN)についても
対応する出力y=(y1,y2,…,yN)の同時分布p(y)が多変量ガウス分布に従う時
xとyの関係はガウス過程(Gaussian process)に従う
過程
ガウス過程が確率過程(stochastic process)の一種であることから来ている
確率過程とは
入力の集合(x1,x2,…,xN)に対応する確率変数の集合(y1,y2,…,yN)に同時分布p(y1,y2,…,yN)を与える確率モデル
Nをいくらでも大きくとれる
無限次元のガウス分布

3.2.1 ガウス過程の意味

Y 〜 N(0, λ2ΦΦT)の意味は?
共分散行列をK=λ2ΦΦTとおく
Φ(X)=(Φ0(x),…,ΦH(x))Tとして
Xnとサミゃの特徴ベクトルΦ(xn)とΦ(xn’)の内積の定数倍が、共分散行列Kの(n,n’)要素Knn’になっている
多変量ガウス分布において2つの変数の間の共分散が大きい
似た値を取りやすい
Knn’が大きい
特徴ベクトルの空間においてXnとXn’が似ている
対応するynとyn’もにた値を持つ
性質3.2 ガウス過程の直観的な性質
入力xが似ていれば、出力yもにている
入力xが似ている
特徴ベクトルΦ(x)を通じた内積が大きいこと
出力yが似ている
ガウス分布での共分散が大きいために、近い値を取りやすい
公式3.3 平均0のガウス過程
y 〜 N(0,K)
観測データはあらかじめ平均を引いておけば平均を0にできるので、まずはこの式で検討

3.2.2 カーネルトリック

yの分布は共分散行列Kの要素 Knn’ = Φ(xn)TΦ(xn’) だけで定まる
特徴ベクトルを明示的に求める必要はない
カーネル関数(kernel function)
k(xn, xn’) = Φ(xn)TΦ(xn’)
カーネル行列(kernel matrix) または Φのグラム行列(Gram matrix)
カーネル関数k(xn,xn’)から計算される共分散
カーネル関数をk(X, x’) = (XTX’+1)2と定義すると
X=(X1,X2)T、X’=(X1,X’2)Tの時
k(x,x’)=(x1x’1 + x2x’2 +1)2
特徴ベクトルを上式として定義する
k(x,x’)=Φ(X)TΦ(x’)
結果を出すためにΦ(x)を書き下す必要はない
カーネルトリック(kernel trick)
特徴ベクトルΦ(x)を直接表現せずに、カーネル関数だけで内積を計算する
Φ(x)が無限次元でも、カーネルトリックにより簡単に計算できる
Kはガウス分布の共分散行列
ガウス過程として成立するためには
Kは対称行列かつ逆行列K-1が存在し、固有値が全て正となる正定値行列

3.2.3 ガウス過程の定義

ガウス過程とは、関数f(*)を出力する箱
無次元のガウス分布
関数fを定めるには、あらゆる入力x1,x2,x3,…における出力値f(x1),f(x2),f(x3),..を定めれば良い
概念図
入力空間X上の無限個の点x1,x2,x3,…に対して、その間の類似度を表すカーネル行列Kが計算できる
これを共分散行列とするガウス分布から無次元のベクトルf(f(x1),f(x2),f(x3),…)が得られる
無限個の入力X=(x1,x2,x3,…)に対応してガウス過程に従う無次元のベクトルをf=(f(x1),f(x2),f(x3),…)と書き、 f〜N(μ,K)で平均がμのガウス過程を表す
fのうち、観測された有限個のf(x1),f(x2),f(x3),..に限って古海苔を周辺化し、有限個のガウス分布としたもの
周辺分布 (marginal distribution)
同時確率分布から、一部の確率変数を消去した確率分布
条件付き確率分布は、特定の確率変数を特定の値に制限したもの
同時確率分布から周辺分布を作ることを周辺化(marginalization)
例
同時確率(joint probability)
二つの事象がどちらも起こる確率
P(X∩Y) = P(X,Y)
周辺確率(marginal probability)
他の事象に関わりなく一つの事象だけの確率
P(X)とかP(Y)
同時確率を全て足し合わせて、P(X)を求めることをYについての周辺化という
加法定理の適用
例:3点X=(x1,x2,x3)上でのガウス過程
x1,x2,x3の間の近さ(=共分散)に応じた楕円体型のガウス分布N(0,K)が定まり、そこからのサンプルとしてf=(f(x1),f(x2),f(x3))が得られる
定義3.4ガウス過程の正確な定義
どんな自然数Nについても、入力x1,x2,x3,…,xN∈Xに対応する出力のベクトルf(f(x1),f(x2),…,f(xN))が平均μ=(μ(x1),μ(x2),…,μ(xN)), Knn’= k(xn,xn’)を要素とする行列Kを共分散行列とするガウス分布N(μ, K)に従う時 Fきガウス過程に従うといいこれをf 〜 GP(μ(x), k(x,x’))と書く
カーネル法とガウス過程
ガウス過程はベイズ統計の立場から見たカーネル法
カーネル関数と平均関数があれば、fの事前分布であるガウス過程が定まり
データが求められればfの事後分布を求めることができる
ガウス過程にもとづく予測分布は
期待値がカーネルリッジ回帰の結果と一致する
これに分散を加えてガウス分布としたものになる
分類のためのSVM(サポートベクトルマシン) のガウス過程による等価物
RVM (Relevance VectorMachine: 関連ベクトルマシン)
RVMはSVMより少ないサポートベクトル数で高い性能を持つ、優れた分類機
RVMは計算量が原則的にO(Nβ)
SVMと異なり学習問題が凸でない
あまり使われていない

3.2.4 ガウス過程からのサンプル

ガウスカーネル(Gaussian kernel) または動径基底関数(radial basis function, RBF)カーネル
ガウス過程で入力xの間の類似度を図るために最もよく使われるカーネル関数
Θ1,θ2∈Rはカーネル関数の性質を求めるパラメータ
概略
RBFカーネルはxとx’の距離に応じて、値がガウス分布のように指数的に減衰する関数
このRBFを用いてランダムにサンプルfを生成した結果
サンプルの間隔を変えることで共分散行列の値も変化
ガウス過程からのサンプルfには無限個の点があり(=無次元)、並べると曲線になる
そのうち一部を抜き出すと右の図のようになる
少数のパラメータを用いた線形回帰モデルのような、従来のパラメトリックモデルでは得られない
RBFカーネルを使って生成した2次元のガウス過程

3.3 ガウス過程とカーネル

3.3.1 RBFカーネルと基底関数

はじめに
RBFカーネルを用いたガウス過程による曲線、あるいは一般に超曲面は、基底関数を無限個重み付けた線形モデルと見ることができる
ガウス過程ではWを積分消去して、データのある場所だけの有限次元で結果のガウス過程を表現できる
データのある場所での関数の線形和によって、その事後分布がガウス分布として表現される
カーネルと関数形
ガウス過程の事後期待値は、その関数が使ったカーネル関数の重みつき和になる
RBFカーネルに限らずガウス過程の一般的な性質
期待値の式
予測したい点x*を変数とするカーネル関数k(x*,・)の、学習データ点での重みつき和になっている

3.3.2 さまざまなカーネル

RBFカーネル以外にも、 xとx’の近さを表すカーネル関数には様々なものがある
線形カーネル(linear kernel)
k(x,x’)=xTx’
ガウス過程からのサンプル
特徴ベクトルが恒等写像(自分自身をそのまま返す変換) Φ(x)=x であることを意味している
xの中にバイアス項に対応する定数x0=1が含まれているとすれば、重回帰と等価
指数カーネル(exponential kernel)
k(x,x’)=exp(-|x-x’|/θ)
ガウス過程からのサンプル
ブラウン運動と関連
Ornstein-Uhlenbeck過程とも呼ばれる
周期カーネル(periodic kernel)
k(x,x’)=exp(θ1cos(|x-x’|/θ2))
ガウス過程からのサンプル
極座標で単位炎上を動く点の偏角と見た場合
周期カーネルでは、 xとx’の距離|x-x’|が大きくても、この距離が周期θ2の2π倍であれば、カーネルの値は大きくなるため、「にている」と見做され、 yの値が周期的に近い値となる
ガウスカーネル
k(x,x’)=exp(-|x-x’|2/θ)
ガウス過程からのサンプル
性質3.5 ガウス過程回帰と線形回帰
ガウス過程回帰は線形回帰モデルを含有している
カーネルの組み合わせ
上記のカーネルは組み合わせて使うこともできる
2つのカーネルの和および積は半正定値を持つ、正しいカーネル関数となる
例
線形カーネルと周期カーネルの組み合わせ
「全体の一次関係+周期性」を持った回帰を行うことができる
Maternカーネル
その他のカーネル
カーネルは実数ベクトルだけではなく、文字列やグラフ、木構造などについても定義できる
複雑なデータ構造の上に、類似度であるカーネル構造を定義するだけで、その空間で考えることができることは、カーネル法の大きな利点
文字カーネル
Fisherカーネル
HMMの周辺化カーネル

3.3.3 観測ノイズ

3.4 ガウス過程回帰モデル

はじめに
ガウス過程に基づいて具体的にどうやって回帰問題を解くのか?
N個の観測値(入力x∈𝝌と出力y∈ℝのN個のペア) D={(x1,y1),(x2,y2),…,(xN,yN))
xとのyの間に y=f(x) の関係がある
この関数fが平均0のガウス過程 f 〜 GP(0,k(x,x’))から生成されている
観測値のノイズもカーネル関数の中に含まれる
y=(y1,y2,…yN)Tとおけばこのyはガウス分布に従い、入力の全てのペア(xn,xn’)についてカーネル関数k(x,x’)を用いて

3.4.1 ガウス過程回帰の予測分布

x*でのy*の値はどうやれば計算できるか?
線形モデルでは重みwを使ってy*=WTΦ(X*)と予測
ガウス過程の場合はWは積分消去されて存在しない
yに新しいy*含めたものを新しくy’=(y1,…,yN,y*)
x1,…,xN及びx*から計算される(N+1)x(N+1)次元のカーネル関数行列をK’
y’〜N(0,K’)
ガウス過程モデルの新しいデータが入った時
新しい入力x*と学習データの入力の類似度 (カーネル関数の値)を並べたベクトル
k*=(k(x*,x1), k(x*,x2),…,k(x*,xN))T
x*の自分自身との類似度
k**=k(x*, x*)
ガウス分布に従うベクトルが上式の時
ベクトルの一部y1が与えられた時の残りのy2は、ガウス分布p(y2|y1) =N(μ2+ Σ21Σ11-1(y1-μ1), Σ22-Σ21Σ11-1Σ12)
μ1=μ2=0の時
p(y2|y1) = N(Σ21Σ11-1y1, Σ22-Σ21Σ11-1Σ12)
公式3.6 ガウス過程の予測分布
p(y*|x*,D) = N(k*TK-1y, k** – k*TK-1k*)
y*の期待値はガウス分布の平均
公式3.7 ガウス過程の予測分布の期待値
𝔼[y*|x*,D] = k*TK-1y
予測値が複数ある場合
入力したい点がX*=(x1*,…,xM*)とM個ある場合
予測したい出力をy*=(y1*,…,yM*)
yとy*の同時分布
公式3.8 ガウス過程の予測分布(予測値が複数個の場合)
p(y*|X*,D) = N(k*TK-1y, k** – k*TK-1k*)
ガウス過程の予測結果
予測分布の計算
実際の予測分布の期待値はどのようにして計算するのか?
公式では上式のような行列とベクトルの積の計算をする
NxNのカーネル行列Kの逆行列K-1を求めることが必要

3.4.2 ガウス過程回帰の計算

基本アルゴリズム
ガウス過程回帰の基本定義
入力xと出力yのペアからなる学習データ(x1,y1),(x2,y2),…(xN,yN)
入力xと入力x’の間の類似度、すなわちガウス分布の共分散を与えるカーネル関数k(x,x’)
学習データX=(x1,…,XN)及びy=y1,…,yN)が与えられた時
テストデータの新しい入力x*に対する出力y*を p(y*|x*,X,Y) = N(kTK-1Y, k** – k*TK-1k*) として求める
カーネル分布が与えられれば一意に定まる
ガウス過程においては「学習」はカーネルのハイパーパラメータのみ
シンプルなガウスカーネル例
ガウスカーネルによる共分散行列の効率的な計算アルゴリズム
XはD次元の縦ベクトルがN個並んだDxN次元の行列
X.^2は行列Xの各要素ごとの2乗
repeat(X,a,b)は行列Xをa行b列に複製する関数
ガウス誤差を含んだ場合のガウスカーネル

3.4.3 ガウス過程回帰の要素表現

ガウス過程の予測分布の期待値
精度行列を𝜦= K-1とすると
その要素𝜦nn’を用いて
ガウス過程のyの予測の期待値は、各学習データの各ynを上式で重みつき和
予測量が学習データの重みつき和で表される
SVMなどのカーネル法と似ている
重み自体の最適化をする
ガウス過程では重み自体の最適化はしない(カーネルのパラメータ以外は学習しない)
リプリゼンター定理
ガウス過程かいきの予測値の(ニューラル)ネットワーク的な表現
y*は単独のx*だけではなく、 x1,…,xN全体に影響される

3.5 ガウス過程回帰のハイパーパラメータ推定

ガウス過程のハイパーパラメーターを推定するためにはどのようにすれば良いか?
ハイパーパラメータをまとめてθ=(θ1,θ2,θ3)とおくと、カーネルはθに依存する
カーネルの式
Kから計算されるカーネル行列Kもθに依存してKθとなる
この時の学習データの確率は
対数をとると
確率を最大にするθ→右の式を最大にするθを求める
求める方法
尤度関数を最大にするように、θを少しづつ変えて探索するMCMC法
勾配法(gradient method)を用いる
∂Kθ/∂θは、カーネル行列Kθの各要素を注目しているパラメータθ∈θで微分した行列
SCG法やL-BFGS法等の最適化ルーチンを用いれば、ハイパーパラメータを最適化できる
SCG(scaled Conjugate Gradient)法は、ガウス過程の最適化において薬使われてきた最適化法
BFGS法やL-BFGS法より安定して計算できる
MATLABやpythonのライブラリにある
L-BFGS法での計算例
データを増やすと最適なハイパーパラメータも変化して、よりデータに合った柔軟な回帰関数が得られる
方法3.9 ガウス過程のハイパーパラメータの最適化
ハイパーパラメータ推定と局所解
勾配法ではθは一般的に局所解で、真の大域的最適解とは限らない
そのような場合は、MCMCを用いて探索することで、局所解を避けて大域最適解に近づくことができる
カーネルの選択
カーネルは組み合わせて使うことも可能
例
各カーネルの重みθ1,θ2,θ3及び誤差の分散θ5を同時に学習することが可能
特定の一つのカーネルを事前に選択する必要はない
例
陸上の100mの世界記録のデータの複数のカーネルモデルでの計算
タイムが短縮される全体的な傾向→線形カーネル
細かイーン堂をガウスカーネルで表現し、ベイズ的にデータの量に応じた確信度(分散)の回帰もて゜るが得られる
学習データの確率
エビデンス(evidence)

3.6＊ガウス過程回帰の一般化

はじめに
さらなる一般化
観測モデルのバリエーション
ガウス過程による関数値fから観測値yが生成される観測モデル
条件つき確率p(y|f)により表される
観測モデルにガウス分布を想定
p(y|f)=N(f, σ2I)
観測モデルにはガウス分布でない場合も考えられる
解析的に解を求めることができない
MCMCやラプラス近似、期待値伝搬法(EP)法、変分ベイズ等の近似的な推論が必要

3.6.1 ロバストなガウス過程回帰

コーシー分布 (Cauchy distribution)
ガウス分布と比較して分布の裾がずっと厚くなる
外れ値(outlier)への対応
Cauchy(x|x0,γ)=1/π(γ +(x-x0)2/γ)
外れ値がある場合のガウス回帰モデル
ガウス誤差を使うとガウス過程回帰、回帰関数が外れ値に引きずられる
コーシー誤差では外れ値に引きずられない
コーシー分布は外れ値に対して頑健(ロバス)である
事後分布
単純なガウス分布ではないので誤差の項をカーネルの中に含むことができない
Fの推定には変分ベイズやMCMC法のような近似推論法が必要
楕円スライスサンプリングを用いてMCMCでサンプリングしたfの例
ガウス分布の標準偏差σがさらにガンマ分布に従うとして場合
σについての期待値をとって拡張した分布
コーシー分布はt分布で自由度が1の特別な場合
ガウス過程を同様に拡張したt過程(Student-t process)も提案されている
誤差ではなく、f自体が外れ値をとることを許す確率モデル

3.6.2 ガウス過程識別モデル

観測値yが何かのイベントが起きた(1)、起こらない(0)のような二値データである場合
Fが与えられたときにy=1となる確率p(y=1|f)を p(y=1|f) = σ(f)と表す
ガウス過程を用いたこうした識別モデルを、 “ガウス過程識別モデル”と呼ぶ
σ(x)はS字の関数
ガウス過程識別モデルのためのロジスティック関数
確率モデルとして最も素直なのは “プロビット関数(probit function)” とも呼ばれる正規分布の累積密度関数
プロビット関数は正規分布の積分を含むため、代わりにロジスティック関数(logistic function)が広く用いられる
σ(1.7x)はΨ(x)に極めて近い関数となる
誤差 σ(1.7x) – Ψ(x)は全域にわたって0.0.1以下となる
Yが与えられた時のfの分布は
単純なガウス分布とはならない
近似的推論が必要
ツールキットGpyでEP法を使って計算したガウス過程識別モデルの推定結果
ガウス過程識別モデルを、少数の関連ベクトルを選ぶことでスパースに行う
RVM(関連ベクトルマシン)
ガウス過程識別モデルは、単に分類だけをしたいときは、モデルも推定も複雑なのでメリットは少ない
他の確率モデルとの組み合わせによりfが得られ、それに基づいて観測値yが得られる場合に必要なモデル

3.6.3 ポアソン回帰モデル

機械の故障回数や素粒子の崩壊など、観測値yが2,4,0,1,…のように非負の自然数をとる場合
ポアソン分布 (Poisson distribution)
ポアソン分布は離散値y=0,1,2,3,…に対する確率分布
期待値は実数値のパラメータλと一致し、 𝔼[Po(y|λ)] = λ
ポアソン分布の例
ポアソン分布は、素粒子の崩壊のように、一定の小さい確率で起きる事象を多数解同時に考えた場合に、全体で起きる事象の回数が従う分布
λ=efと書けるとすると
p(y|λ) =Po(y|λ) =λye-λ/y! =p(y|f) = (ef)ye-ef/y! =exp(fy -ef)/y!
事後分布はガウス分布ではない
例: 一次元で横に並んだ50区間における、特定の植物の個体数を表すデータをガウス過程によるポアソン分布及びラプラス近似を用いて推定した結果
ここでのMCMCはガウス過程モデルの(低次元の)ハイパーパラメータを探索するもの
コラム：ニューラルネットワークとガウス過程
1層のニューラルネットは隠れ素子数→∞の極限でのガウス過程と等価
ニューラルネットワークの代わりに、ガウス過程を考えることで、ニューラルネットワークにおける多重の重みの最適化が不要になる
何が学習されるか予測できないニューラルネットワークと違って、カーネル関数を通じて問題な関する事前知識を表現したり、時系列やグラフのように自明にベクトルかできない対象を見通しよく扱うことができる
ニューラルネットワークのパラメータを事前分布からランダムにサンプリングした時の、入力(x1,x2)=(-0.2,0.4)に対する主査力の分布
各点一つのニューラルネットワークを表す
隠れ素子の数Hが増えるにつれ、出力の分布がガウス分布に近づく
深層学習のようにネットワークが多層になっている場合
共分散行列を再帰的に計算することで、等価なガウス過程(NNGP, Neural Network GP)を求められる
NNGPは、ニューラルネットワークのような SGDなどで重みを明示的に学習することなく、解析的に表されるガウス過程で同等の性能を持つことが締められた
ニューラルネットワークがガウス過程に漸近するためには、独立で有限の分散を持つ同時分布に従ってさえいれば良い

第４章確率的生成モデルとガウス過程

はじめに
確率生成モデル
「現実世界における観測Yは、何らかの確率分布p(Y)からのサンプリング Y〜p(Y)によって得られたもの」という仮説
最も基本的なモデル
Y = f(x) + ε
入力xを定数として固定
関数f(*)、観測値y、観測ノイズεを確率変数遠く
確率生成モデリングの基礎になる考え方
尤度関数の導入は「目の前の観測値yは、確率変数だ」という発想の転換であり、最尤推定法の基礎になる
事前確率・事後確率の導入は「知りたいもの(隠れ変数fもしくは未知パラメータしー)は確率変数だ」という発想の転換であり、ベイズ推定の基礎になる

4.1 確率変数と確率的生成モデル

4.1.1 確率変数Xと確率分布p(X)

確率変数(random variable)
その値が試行ごとに確率的な決まる変数
例:確率変数Xが理想的なサイコロの出す目とする
Xの撮る値の集合は{1,2,3,4,5,6}
それぞれの出現確率はp(X=1)=1/6,…,p(X=6)=1/6
全ての場合の確率の和は1になる
確率分布(probability distribution)
一般的に
確率変数Xの性質は、 Xがとり得る値の集合と、 Xがこれからの値をとる確率を定める “確率分布(probability distribution)p(x)” から定まる
確率分布p(X)は、Xの実現値x∈Xの出現確率p(X=x)を全ての実現値について定める
例
確率変数Xを、的の中心を狙って投げたダーツの矢の刺さった位置を的の中心から測った距離
実現値がとりうる値は実数 x∈ℝ
xの出現確率は確率密度関数p(x)によって表される
xが特定の範囲内に含まれる確率は積分で表される
確率変数がとりうる全ての値に関する確率の積分は1になる
平均μ、分散σ2の一次元ガウス分布の確率密度関数
確率過程(stochastic process)
関数f(x)を確率的に生成する生成器
定義4.1 確率過程
任意のN個の入力値x1,…,xN∈Xに対して
N個の出力値fN=(f1,…,fN)=(f(x1),…,f(xN) , f(xn)∈yの同時確率p(fN)=p(f1,…,fN)を与えることができる時
この関係f(*)を確率過程と呼ぶ
ここでNは任意の自然数、XとYはそれぞれ入力値と出力値がとりうる値の集合
定義4.2 ガウス過程
確率f(*)において、同時確率p(f1,…,fN)がN次元のガウス分布として得られる時
f(*)をガウス過程と呼ぶ
メモ1
変数xの確率分布がパラメータμに依存する条件先分布であることを強調したい場合
p(x|μ)=N(y|μ、σ2)、もしくはx|μ〜N(μ、σ2)とかく

4.1.2 同時確率p(X,Y)と周辺化

複数の確率変数の同時分布 (joint distribution)を考える
例:確率変数XとYの同時分布をp(X,Y)と書く
2つの確率変数X,Yの組み合わせ(X,Y)を新たに1つの確率変数と考えた時の確率分布
3つ以上でも同様
Xがとりうる値の集合を𝓧、Yがとりうる値の集合𝓨
XとYの組み合わせ(X,Y)がとりうる値の集合を𝓧x𝓨
𝓧と𝓨の直積集合(Cartesian product)
(X,Y) ∈ 𝓧x𝓨
2つ以上も同様
定義4.3 周辺化と周辺分布
複数の確率変数X,Yの同時分布の確率密度関数p(X,Y)から
以下のような積分計算によってXの確率密度関数p(X)を得る操作を周辺化(marginalization)
p(X) = ∫ p(X,Y)dY
こうしてできた確率分布p(X)を周辺分布(marginal distribution)と呼ぶ
「Yに関して周辺化積分する」とか、「周辺化によりYを消す」という
Yがとりうる値の集合が𝓨が離散的な時は、積分ではなく全ての和をとる
定義4.4 条件先確率
条件付き分布(conditional distribution) p(Y|X)は
Xの値が既知もしくは所与とした時の、Yの確率分布
条件Xは必ずしも確率変数である必要はない
条件Xに関しては積分にも特に制約はない
p(Y|X)がYの確率密度関数である時 ∫ p(Y|X)dY =1が必ず成り立つ
条件付き分布、同時分布、周辺分布の間に成り立つ関係
P(Y|X)p(X) = p(X, Y) = p(X|Y)p(Y)
定義4.5 ベイズの定理
p(Y|X) = p(X|Y)p(Y) / p(X)
例: サイコロ・ダーツの連鎖的モデル
サイコロを振り、その出目dによって的の位置座標μを既知の関数関係μ=μ(d)によって決めそこを狙ってダーツを投げ、刺さった一の座標の水平座標をx=X1∈ ℝとする
xの確率分布を確率密度関数p(x)の形でどうかくか?
xの生成過程はサイコロの出目の確率分布p(d)とダーツなげ結果の位置座標の確率分布p(x|d)の仏で決まる
サイコロの出目の確率分布
p(d=1)=1/6, …,p(d=6)=1/6の一様分布
位置μを狙って投げたダーツの刺さる位置
Μを中心として適当な分散σ2を持つ正規分布に従う
p(x|μ) = N(x|μ,σ2)
ダーツで狙う位置座標μは、サイコロの出目に依存
p(x|μ) = p(x|μ(d)) = N(x|μ(d),σ2)
XとDの同時分布p(x,d)
2つの確率分布の積
p(x,d)=p(x|d)p(d)
同時確率p(x,d)をdに対して周辺化することで求めたい密度関数p(x)を得る
具体的な密度関数
複数のガウス分布の重みつき平均の形で書くことのできる確率分布
混合ガウス分布 (mixture Gaussian distribution)
手順の抽象化
1 未知の値を確率変数で表す (例 xとδ)
2 個々の値が生成される確率的な過程をそれぞれ確率分布で表す (例 p(d)とp(x|d))
3 全ての確率変数の同時分布を表す (例 p(x,d) =p(x|d)p(d))
4 必要な周辺分布を求める (例 p(x))
確率的生成モデル (probabilistic generative modeling)
観測される値の生成過程を確率分布で表す作業
例4 ガウス過程回帰モデル

4.1.3 独立性：p(X,Y)=p(X)p(Y)と条件付き独立性：p(X,Y|Z)=p(X|Z)p(Y|Z)

複数の確率変数の関係を調べる時、最も基本的かつ重要な接質は独立性と条件付き独立性
定義4.6 確率変数の独立性
2zの確率変数XとYに関する条件付き分布、同時分布、周辺分布の間に、石かの性質が成り立つ時
確率変数XとYは独立(independent)であるという
p(X,Y)=p(X)p(Y)
この関係はp(X)=p(X|Y)と等価
さらにp(Y)=p(Y|X)と等価
定義4.7 確率変数の条件先独立性
2つの確率変数X,Yと条件Zが以下の関係を満たす時
確率変数X,Yは条件Zのもとで条件付き独立(conditionally independence)であるという
p(X,Y|Z) = p(X|z)p(Y|Z)
例5 3つ以上の確率変数の独立性や条件付き独立性
p(A,B,C|D,E)=p(A|D,E)p(B,C|D,E)が成り立つ時
条件(B,C)の組み合わせをF、条件(D,E)の組み合わせをGと呼び直すことで
p(A,F|G)=p(A|G)p(F|G)
条件(D,E)のもとでAと(B,C)が条件付き独立
例6 多変量ガウス分布における独立性
例7 ガウス過程回帰における条件付き独立性
p(y|f) =p(y1|f)x…xp(yN|f)
p(y|f)=N(f,σ2IN)
fを条件とした条件付き独立
p(y)≠p(y1)x…p(yN)
p(y)に関しては成分間の独立は成り立たない
例8 独立同分布と条件付き独立性
同じ確率分布を持つ確率変数から、繰り返し標本を生成することによって複数の実現値を得ること
独立同分布(independent and identically distributed, iid)

4.1.4 ガウス過程回帰モデルのグラフィカルモデル

グラフィカルモデル(graphical model)は、確率変数の間の独立性及び条件付き独立性の有無を一覧性の良い図で可視化するための記法
例5のケース
グラフィカルモデルでの確率変数の抽象化あるいはパネル表示
例9 線形回帰の確率生成モデル
グラフィカルモデル
例10 ガウス過程回帰の確率生成モデル
グラフィカルモデル
fが全て矢印の無い無向きリンクで繋がっている

4.2 最尤推定とベイズ推定

4.2.1 確率的生成モデルと最尤推定

観測Yの確率生成モデル (probabilistic generative model)
「現実世界における観測Yは、何らかの確率分布p(Y)からのサンプリングY〜p(Y)によって得られたものだ」
確率生成モデルは仮説
確率生成モデルは同じ対象に対して複数のモデルを持っても良い
p(Y)=P1(Y),p(Y)=P2(Y),…
仮説間の違いをパラメータθで表して、確率生成モデルを条件付き確立p(Y|θ)で表すこともある
p(Y|θ)をパラメトリックモデル (parametric model)と呼ぶ
パラメータの推定(estimation)
観測Yに基づいてパラメータθを決定すること
最尤推定 (maximum likelihood estimation)
パラメータを含む確率的生成モデルp(Y|θ)のパラメータθを
観測データの尤度関数(likelihood function) L(θ)=p(Y|θ)が最大になるように定める方法
尤度関数L(θ)
確率関数p(Y|θ)において観測データYを定数と扱うことで、パラメータθの関数と見立てるもの
尤度関数を最大にするパラメータをパラメータの最尤推定値と呼ぶ
Arg maxθ L(θ)は尤度関数L(θ)が最大値をとるようなパラメータθの値
例11 共分散行列が既知であるガウス分布の平均の最尤推定
観測されたN個のD次元縦ベクトルyn ∈ ℝD, n=1,…,Nが
それぞれ平均μと共分散行列Σを持つガウス分布モデルから生成されたものと仮定
確率的生成モデル
yn 〜 N(μ, Σ)
尤度関数は未知の平均ベクトルμの関数として書ける
両辺の対数をとる
対数尤度をμに関して微分したものがゼロになる点(留意点)
最尤推定量
例12 サイコロを使ったダーツ投げモデルの最尤推定
確率的生成モデル
狙ったダーツの位置座標θ=(μ(1),…,μ(6)}が未知パラメータの時の最尤推定?
尤度関数
対数尤度関数

4.2.2 確率的生成モデルとベイズ推定

ベイズ推定では「知りたいものは(未知パラメータθ)は、確率変数だ」と考える
ベイズ推定において未知パラメータθの「推定」とは、観測Xを得ることによって確率変数θの分布を更新すること
確率変数の分布とは
確率密度関数
更新前の分布p(θ)を事前確率(prior もしくは apriori probability)分布
観測Xに基づいた更新後の分布p(θ|X)を事後確率(posterior 若しくは a posteriori probability)分布
ベイズ推定の手順
尤度関数をモデリングする。
未知変数θのもとで観測Yが確率的に生成される過程を、条件付き確率p(Y|θ)で書き下す
これを未知変数θの関数L(θ)=p(Y|θ)とする
事前確率をモデリングする
未知変数θのとり得る値の範囲と、その主観的な可能性の度合いを、確率分布p(θ)で表す
ベイズの定理を適用して未知変数θの事後確率p(θ|Y)を求める
p(θ|Y) = p(Y|θ)p(θ) / p(Y)
p(Y)は、p(Y) = ∫ p(Y|θ)p(θ)dθ で定義される
観測値Yの実現値が所与である時、θによらない数値にもなるが
p(Y) は正規化定数(normalizationconstant) 若しくは、周辺尤度(marginalized likelihood)と呼ばれる
ベイズ推定と最尤推定は尤度関数に基づいて未知のパラメータを知ろうとする点では同じ
「ベイズ推定において未知のパラメータや予測分布が確率分布の形で得られる」ことが大きな特徴
ガウス過程法では未知ベクトルの値について平均と共分散を持つガウス分布が得られる
例13 共分散行列が既知であるときのガウス分布のベイズ推定
観測されたD次元縦ベクトルyn∈ℝD,n=1,…,N
上記の確率的生成モデルが平均μと共分散Σを持つガウス分布と仮定
パラメータのうちΣは既知
未知パラメータμの事前確率はガウス分布で与えられる
yn 〜 N(μ, Σ), n=1,…,N
μ 〜 N(μ0, Σμ0)
μ0は既知のD次元縦ベクトルであり、事前分布の中心
DxDの大きさを持つ対称行列Σμ0も既知であり、事前分布の共分散行列を表す
ベイズ推定の目的は、未知パラメータμの事後分布p(μ|Y)を求めること
ベイズの定理(の対数表現より)
log p(μ|Y) = log p(Y|μ) + log p(μ) – log p(Y)
右辺の第一項は対数尤度関数そのもの
右辺の第二項は事前確率
右辺の第三項はμによらない定数
上式の変形
Constはμに依らない定数をまとめたもの
事後確率p(μ|Y)は正規分布N(μ,Σμ)の形で求まる
分散の逆数
精度(precision)
共分散行列の逆行列
精度行列(precision matrix)
ガウス分布を共分散行列Σではなく、精度行列S=Σ-1を用いる
事後確率のパラメータ
事後確率の精度行列Sμは、事前確率に由来する項Sμ0ょ除けば、観測の確率生成モデルの精度行列のN倍になる
観測ynが加わるごとに、事後確率の精度行列Sμの値がSだけ増えていく

4.3 確率分布の表現

はじめに
確率分布を求めるとはどういうことか?
具体的にどのような作業をすれば良いのか

4.3.1 ノンパラメトリックモデルとは

確率変数の確率分布を計算機上の数値によって表現する方法は大きく分けて2つある
パラメトリック(parametric)な手法
推定対象が従う確率分布が一定次元のパラメータで指定される確率分布
このパラメータを求めることによって確率分布を決定する
代表
ガウス分布
ベータ分布
確率分布がパラメトリックな確率分布に従っていると仮定
パラメータを数値的に求めることがすなわち確率分布の推定である
事前確率、尤度確率、事後確率などをパラメトリックな確率分布で表すこと
パラメトリックな確率モデリング
ノンパラメトリック(nonparametric)な手法
推定対象について既知のパラメトリックな確率分布を仮定せずに確率分布を求める方法
典型例
ヒストグラムを求める方法
カーネル密度推定法
K近傍法
ガウス過程法はノンパラメトリックな方法の一種
推定対象となる未知の関数f(x)の確率分布としていかなるパラメトリックな確率分布も仮定しない

4.3.2 確率分布を標本で表現する

確率分布p(x)を計算機上で表現して可視化する直接的な方法
xがスカラーであれば密度関数p(x)や分布関数F(x)=∫p(t)dtの関数形状を示す
間接的な方法
確率分布から生成した適当な個数の標本X1,…,XNを可視化する
確率分布の様々な表現(1次元)
確率分布の様々な表現(2次元)
有限個の標本を用いて確率分布の可視化を行う有用性
確率密度関数が解析的な形で書けない(若しくは書き表すことが難しい)
標本を計算機で生成するのは容易である
例14 解析的な形で書き表せない確率分布の可視化
範囲(-1,1)から等確率で値dが生成される一様分布と、 dを中心とした分散σ2の正規分布を連鎖系で繋ぐ
d 〜 Unif(-1, 1)
x 〜 N(d, σ2)
周辺分布p(x)はどのような形になるか?
標本による可視化
(A) 確率変数xの散布図
(B) 1000点の標本のヒストグラム
(C)青線はガウス分布カーネルを用いたカーネル密度推定の結果赤の破線は、真の確率密度関数
標本は統計量を求めるためにも使える
関数f(x)=x2が与えられていて
期待値𝔼[f(x)] = ∫ f(x)p(x)dx
分散𝕍[f(x)] = ∫(f(x) – 𝔼[f(x)])2p(x)dx
重み付き標本を用いる工夫
確率分布p(x)をN個の重み付き標本(wnxn), n=1,..,Nで表す
期待値や分散は重み付き平均を用いて計算できる
重み付きサンプリングのアルゴリズム
例15 解析的な形で書き表せない確率分布に関するベイズ推定
例14で作ったモデルのパラメータσ2が未知
事前確率p(σ2)=Unif(0.1,2)と確率変数xについてN個のi.i.dな標本XN=(x1,…,xN)が与えられている
事後確率はどのようにして求めるのか?
パラメータσ2を解析的な形で書き表すのは難しい
事後確率に従うM個の重み付き標本wmσ2, m=1,…,Mを生成することは可能
ベイズの定理から
p(σ2|XN) = p(σ2)p(XN|σ2) / p(XN)
尤度
エビデンス
XNとσ2の値を入力したときに、尤度p(XN|σ2)を数値的に求めることができるなら
事後確率分布p(σ2|XN)の重み付き標本を得ることが可能
事後確率分布の重み付き標本を得るアルゴリズム
隠れ変数が含まれる場合に事後確率分布の重み付き標本を得るアルゴリズム
ベイズ推定の結果を標本や重み付き標本の形で得る方法
モンテカルロ法(Monte Carlo method)
標本系列の作り方に工夫の加わったマルコフ連鎖モンテカルロ(Markov chain Monte Carlo method MCMC法)
例16 カーネル密度推定
カーネル密度推定は、確率分布推定のための最も簡単な方法の一つ
ノンパラメトリックな方法の代表
確率変数x∈ℝdが確率密度関数p(x)に従うとき
B個の標本x1,…,xBを用いて密度関数の推定値Ṗ(x)を上記のように得る
hk(x)≥0, k=1,…,Kはk番目の基底関数
∫ hk(x)dx = 0 の正規化条件を満たす
一つ一つの基底関数をパラメトリックな形で表現し、パラメータを決め、個数Kを適当に決めることで様々な確率密度p(x)をṖ(x)で近似できる
カーネル密度推定は、上記の特別な場合
B個の標本それぞれに対応する基底関数を以下のように定義した場合
hb(x) = h((xb – x)/ σ)
h(x)は正値をとる任意の関数:カーネル関数
σ:バンド幅(基底の半径を制御)
カーネル関数はb=1,…,Bによらず、原点にピークを持つ確率密度関数を適当に与える
具体例
ガウス分布カーネル
h(x) = N(x|0,Id)
エパネクニコフカーネル (Epanechnikov kernel)
h(x) = 3/4(1-∥x∥2)⫿(∥x∥<1)
⫿(*)は()の中がなりたてば1,そうでなければ0ょとる指数関数
確率密度関数として滑らかで見栄えの良いものを得るためにガウス分布カーネルが使われる
見栄えを度外視して密度関数推定の精度について論理的な根拠が必要な場合にエパネクニコフカーネルが使われる
例17 ニューラルネットワークを用いた分布推定法
適当な次元Dのガウス分布から生成されたベクトルパターン x 〜 N(0, ID)
上記をニューラルネットワーク f(x;w)に通して変換する確率生成モデル y= f(x;w)
上記を通して画像や音声などを表す高次元ベクトルパターンyの確率分布を表現
画像分類のための畳み込みニューラルネットワーク
画像生成のための虐畳み込みニューラルネットワーク
コラム：ブラウン運動とガウス過程
ブラウン運動する粒子のx座標をX(t) ∈ ℝ
時刻t+∆t,(∆t>0)に位置との間の差分はガウス分布を用いて X(t+∆t) – x(t) 〜 N(0, ∆tσ2)
W(t):ウィナー過程
W(t+∆t) – W(t) 〜 N(0, ∆t)
一次元のガウス過程
μ(t)=0
共分散関数

第５章ガウス過程の計算法

概要
データ点数Nが大きい場合のガウス過程の計算
カーネル行列やその逆行列を求めるための計算コストがボトルネック
直接観察されない隠れ変数を間に挟む「補助変数法」

5.1 ガウス過程回帰の計算コスト

ガウス過程回帰の主目的
予測分布に基づく出力yの期待値を求めること
計算式
ベクトルk*、行列Kを計算してメモリに格納する
逆行列K-1を計算してメモリに格納する
k*TK-1yを計算してメモリに格納する
計算コスト
メモリ消費量
演算量
ガウス過程回帰の計算コスト
影響を与える因子は、教師データとして与えられる入力点数N,角入力の次元数D
メモリ消費量
数値データを格納するのに必要なメモリ量
O(ND)
ベクトルk*と行列Kを格納するのに、それぞれO(N),O(N2)
演算量
行列KのNxN個の要素を全部計算するのに必要な演算数のオーダー
O(N2D)
逆行列K-1の計算に必要なオーダー
O(N3)
ベクトルk*と行列K-1の赤の計算
O記法(Big-O notation)
NやDなどのスケール因子と計算量の関係の強さを示す表示法
ガウス過程を普通に計算すると、標本サイズNが大きくなるつれ計算コストが嵩む
O(N3)の計算
N=1000で1秒
N=2000では8秒
N=10000では1000秒(17分)
ガウス過程回帰の計算量を考える上でのボトルネック
逆行列計算に必要なメモリ消費量O(N2)と演算量O(N3)
ガウス過程回帰計算の工夫の効果
補助変数法
M次元(M<N)の補助変数を導入することで、ガウス過程ほうで必要な逆行列計算の対象をMxM行列に圧縮する
メモリ消費量のオーダーをO(NM + M2)に、演算量のオーダーをO(NM2 + M3)に減らす
変分法
変分法を導入することで、ガウス過程のハイパーパラメータ最適化が効率化できる
ハイパーパラメータ最適化は、ハイパーパラメータを逐次的に更新するアルゴリズムの 1ステップごとにカーネル行列の更新を必要とする
変分法を導入することで、データをミニバッチ(mini batch)と呼ばれる100個や500個の小口サイズにわけて、逐次的に学習する計算方法(“ミニバッチ学習の概要とアルゴリズム及び実装例“で述べているミニバッチ学習)が可能となる
補助入力てんが格子状に配置されている徳べな場合に格子状構造を用いて大幅に計算効率をアップする
気象などの空間データや、画像・音声・動画などの時系列データを扱う問題で、空間を密に埋めるような入力や出力を扱うときに威力を発揮

5.2 補助変数法

はじめに
補助変数法 (inducing variable method)
ガウス過程回帰法の計算量削減のために有効な近似法

5.2.1 部分データ法

部分データ法 (subset of data approximation , SOD)
Nが大きすぎてNxNの行列計算を避けたいときに使える簡単なアイデア
普通のガウス過程回帰
回帰モデルy=f(x)+εに従って N個の入力点X=(x1,…,xn), xn∈ℝDにおける観測値Y=y=(y1,…,yN)T, yn∈ℝの裏に
ガウス過程に従う関数f(*)とその出力値f(f1,…,fN)T=(f(x1),…,f(xN))T ∈ ℝNを想定
ガウス過程を平均μ(x)とカーネル関数k(x,x’)で定めたとき
ガウス過程回帰モデルは未知数fおよびyに関して以下のガウス過程を想定
f 〜 N(μN, KNN)
μNは第n成分がμ(xN)であるようなN次元縦ベクトル
KNNは第(n,n)成分がk(xn, xn’)であるようなNxN行列
yn 〜 N(fn, σ2)
新しい入力点x*における関数出力値f*=f(x*)と観測値y*を考える
fとf*の同時確率は以下のガウス分布で表される
μ*=μ(x*)およびk**=k(x*,x*)はそれぞれスカラー
kN*は第n成分がk(xn,x*)であるようなN次縦ベクトル
観測ノイズの分散がσ2であることから
予測分布p(f*|y) = N(ḟ*,σ*2)の平均と分散は以下のようになる
f* = k*N(KNN + σ2IN)-1
σ*2 = k** – k*NT(KNN + σ2IN)-1kN*
部分データ法
N個の入力点から M個の入力点(M<N)を全データの分布をよく代表してくれるように選び出す
MxMの行列の逆行列計算によりガウス過程の予測分布が得られる
部分データ法は、選び出されたM個の入力点からなる部分データが全データをよく代表してくれるなら
O(N3)よりもはるかに小さい計算コストO(M3)で同程度の精度が得られる

5.2.2 補助入力点と補助変数法の計算

部分デー法に洗練を加えたもの
補助変数法では全データをよく代表してくれるようにM点からなる部分データを作り出して、これを使って試験入力点における環数値を効率的に推定する
部分データ法ではN-M点文のデータを捨てていた
補助変数法では全データを有効に使うことで精度を保つ
補助変数法では、関数は8*)の定義域内に補助入力点(inducing point)と呼ばれるM個の仮想的な入力点Z=(z1,..,zM)を適切に配置する
適切な補助点が既に得られたものと仮定
補助入力点zmにおける出力値f(zm)を “補助変数(inducing variable)”
上記をm=1,…,Mに関して縦ベクトルにまとめたものを “補助変数ベクトル(inducing vector) u=(u1,…,uM)T
ガウス過程に従うことがわかっている未知の関数f(x)の予測分布を求める手順
1.　事前確率の定義
入力点Xにおける出力値fと、補助入力点Zにおける出力値(補助変数ベクトル)uと、試験入力x*における出力値f*の事前確率p(f,u,f*)を定義する
2. 事後確率を導出
入力点Xにおける観測値yに基づいて、補助変数の事後確率p(u|y)を求める
3. 予測分布を導出
出力値f*の予測分布p(f*|y)を、補助変数の事後分布p(u|y)を用いて近似する
p(f*|y) = ∫ p(f*|f)p(f|y)df ≈ ∫ p(f*|u)p(u|y)du
補助変数法には様々な流儀がある
FITCと呼ばれる方法を採用
補助変数法における事前確率の定義
f(*)が平均関数u(x)と共分散関数k(x,x’)で定義されたガウス過程
出力fと補助変数uの同時事前確率は、上記のガウス分布となる
μM=(μ(z1),…,μ(zM))TはM次元縦ベクトル
補助変数uに対する平均ベクトル
共分散行列Kは(N+M)x(N+M)の実数値正方行列
ベクトルfの分散共分散行列KNN、 uの分散共分散行列KMM、 fとuの共分散行列KNMがKに含まれる
補助変数法では、新しい点x*における予測分布を以下のように近似したい
p(f*|y) = ∫ p(f*|f)p(f|y)df ≈ ∫ p(f*|u)p(u|y)du)
補助点Zが適切に配置されていて、新しい点x*が入力点x1,…,xNのいずれかの周辺にあるときに上式が成り立つならば
新しい点の代わりに各入力点xnを入れても以下の式が成り立つ
p(fn|y) ≈ ∫ p(fn|u)p(u|y)du
補助変数法では以上の考察に基づいて、関数fの確率生成モデルとして以下の3段階の生成過程を経る
公式5.1 補助変数法が想定する3段階生成過程
(1) u 〜 N(0M, KMM)
(2) fn|u 〜 N(kMnTKMM-1u, kn – kMnTKMM-1kMn), n=1,…,N
(3) yn|fn 〜 N(fn, σ2), n=1,…,N
ここでkn = k(xn,xn), kMn = (k(z1,xn),…,k(zM,xn))T
3段階生成過程のうち最初の2段階分を一つの同時分布にまとめる
p(f|u)p(u) = p(f,u)
N+M次元のガウス分布
共分散行列
補助変数法ではその非対角成分全てを無視してゼロに置き換え
対角成分kn=k(xn,xn), n=1,…,Nだけを残す
行列にゼロ要素が多いとき、その行列を”疎行列(sparse matrix)” もしくはスパース行列と呼ぶ
補助変数法の別名
疎ガウス過程sparse Gaussian process)
グラフィカルモデルでの解釈
オリジナルの確率的生成モデルのグラフィカルモデルでは、 fNの全ての成分間に無向リンクが存在
補助変数法のグラフィカルモデルでは、 fNの成分間で互いに直接のリンクを持ち無くなり、その代わりに補助変数uを通じた繋がりになる
新規入力点x*に対応する出力f*も補助変数uを通して観測値y1,..,yNとつながる
補助変数の事後確率の導出

5.2.3 補助変数法の計算コスト

観測点数N=10000であればGPで800メガバイト、補助変数法を使うと補助入力点M=100で9メガバイト

5.2.4 補助入力点の配置

補助変数法の確率的生成モデルはガウス過程回帰の確率的生成モデルの「大胆な近似」
適切な近似になっているかどうか
補助入力点Z=(z1,…,zM)における、関数値f(zm)が全データをよく代表しているか否かに依存する
補助入力点の置き方や密度を変えた場合の近似の比較
関数が滑らかであれば、入力点数の1/10程度の補助入力点で、平均と分散をぴったり近似できる
補助変数法は、補助点の間で補助点の値を滑らかに補完することで関数を近似する
滑らかでない複雑な形の場合は、補助点を多く取らないと良い近似はできない
補助入力点数の数の増加を抑えつつ近似精度を高く保つために
補助入力点を適切に配置する工夫が必要になる
推定する関数f(x)が滑らかな領域は補助入力点の密度を低くする
滑らかでない場合は高めに設定流する
実際の補助入力点の作り方
入力データ点の中から、無作為に一定割合を選び出した「部分データ」を補助入力点とする
入力データ点に関するクラスタリング処理(K平均法など)で、代表点を求めて補助入力点とする
上記方法で選んだものを初期値として、補助入力点位置が最適になるように自動調整する
補助入力点を格子状に並べる
格子状配置の特殊性を用いることで計算量を大幅に削減する方法がある

5.3 変分ベイズ法と確率的勾配法

はじめに
変分ベイズ法(Variational Bayesian method,VB法)に基づいて、確率的勾配法(stochastic gradient method)アルゴリズムを導出する
未知隠れ変数やパラメータが複雑な階層構造を持つようなモデルのベイズ推定問題を数値最適化問題に帰着させる方法
確率勾配法は、数値最適化問題をパラメータの逐次的な更新により解くアルゴリズム
巨大なデータに基づいてニューラルネットワークなどの複雑なパラメトリックモデルのパラメータ決定を効率化するための必須の方法
データ数Nが大きい場合のガウス過程回帰モデルの補助変数法においてハイパーパラメータθの最適化を行いたい場合に、計算効率を改善できる

5.3.1 変分ベイズ法と独立分解仮定

変分ベイズ法はベイズ推定の近似計算法
ガウス過程回帰に限らない一般の機械学習を広く含む状況設定
観測Yの確率的生成モデルp(Y|w,θ)が未知変数wとハイパーパラメータθを条件とした条件つき確率の形で与えられ
未知変数の事前確率p(w)が与えられている場合
ベイズ推定の目的
未知変数の事後確率p(w|Y,θ)とハイパーパラメータθを求めること
事後確率はハイパーパラメータのもとでベイズの定理(左式)を満たすように求める
p(w|θ) = p(Y|w,θ)p(w|θ) / p(Y|θ)
θ* = arg max p(Y|θ)
p(Y|θ) = ∫ p(Y|w,θ)p(w|θ)dw
エビデンス(周辺尤度)
変分ベイズ法でのアプローチ
事後確率分布p(w|Y,θ)を変分事後分布(variational posterior) q(w)によって近似
解析者が変分事後分布q(w)を既知のパラメトリックな事後分布の形で天下り的に与える
例
q(w) = N(w|m,S)のように、パラメータθw = (m,S)に依存するガウス分布とする
次に、変分事後分布を真の事後分布に近づけるべくパラメータの調整を行う
変分事後分布q(w)の意味について
変分ベイズ法アルゴリズムの立場から
変分事後分布は仮の事後分布である
仮の事後分布はアルゴリズム開始時に適当に初期化され、ステップごとに真の事後分布に近づいていく
最終的に求めたい分布に収束する
確率的生成モデリングの立場から
変分事後分布は事後分布のモデルを意味する
事前分布p(w)と尤度関数p(Y|w,θ)がわかりやすい形をしていても、事後分布がわかりやすい形になるとは限らない
解析者は扱いやすいパラメトリックな確率密度関数の形を天下り的に設定し、これを変分事後分布とする
解析者は事後分布に対して任意に制約を加えて、パラメトリックな変分事後分布にする
変分事後分布は仮説の一部
変分ベイズ法の計算目的
wに関する2つの確率分布q(w)とp(w|Y,θ)の違いを最小化すること
二つの確率分布の違いをKL情報量 (Kullback-Leibler divergence) KL[q(w)||p(w|Y,θ)] を用いて測る
定義5.5 KL情報量
確率密度関数q(w)とp(w)で定義された 2つの確率分布間のKL情報量は以下で表される
KL[q(w)||p(w)] = ∫ q(w)log (q(w)/p(w)) dw
任意の確率密度関数p,qに対して KL[q(w)||p(w)] ≥ 0が成り立ち
全てのwにおいて q(w) = p(w)であるときに、 KL[q(w)||p(w)] = 0
一般にKL[q(w)||p(w)] = KL[p(w)||q(w)]は成り立たない
KL情報量をqとpの距離とは呼べない
距離の代わりに”擬距離”と呼ぶ
変分ベイズ法ではq(w)を求めるために、事後分布q(w)の汎関数(functional)の形でエビデンス上階(evidence lower bound, EBLO)もしくは、エビデンスの変分下界(variational lower bound of evidence)を定義してこれを変分事後分布に関して最大化する
定義5.6 エビデンスの変分下界
Fθ{q(w)} = ∫ q(w) log (p(Y|w,θ)p(w|θ) / q(w) dw
なぜエビデンスの変分下界と呼ばれるのか?
エビデンスp(Y|θ)に関する不等式 Fθ{q(w)} ≤ p(Y|θ) が任意の確率密度q(w)について成立し、
上式の等号はq(w)≡p(w|Y,θ)であるときに成立するから
Fθがエビデンスの変分下界であることや KL情報量がゼロ以上であることの証明はイェンセンの不等式(Jensen’s inequality)を用いる
公式5.7 イェンセンの不等式
確率変数x∈ℝと、xの上に凸な関数f(x)に関して、 f(𝔼[x]) ≥ 𝔼[f(x)] が成立する
例
log xは上に凸な関数
log ∫ p(x)f(x)dx ≥ ∫ p(x) log f(x) dx
変分事後分布q(w)の形状がパラメータθwの関数で表される場合に、エビデンスの変分下界は変分事後分布のパラメータθwと、確率的生成モデルのハイパーパラメータθの関数として表現できる F(θw,θ) = Fθ{q(w)}
左辺はθw,θの関数、右辺は確率密度関数q(w)の汎関数(かつθの関数)
エビデンス下界関数
変分ベイズ法では、ベイズ推定を構成する互いに絡み合う2つの問題
事後分布p(w|Y,θ)の推定
ハイパーパラメータθの最適化
エビデンス下界関数を目的関数とした、たった一つの数値最適化問題に帰着させる
[θw*,θ*] = arg max F(θw, θ)
エビデンス下界関数を最大にするような θw*,θが求まったとき
これに基づく変分事後確率q*(w)は事後確率p(w|Y,θ*)の最も良い近似になる
特に、パラメトリックな確率分布q(w)=q(w|θw)がパラメータ次第で真の事後確率を表現可能な場合は
q*(w)=q(w|θw)は厳密な事後確率と一致する
同時にθ*はエビデンスを最大にするハイパーパラメータになる
変分ベイズ法における変分事後確率q(w)の与え方は?
解析者が天下りに与える仮説の一部
独立分解仮説 (factorization model)
変分ベイズにおける変分事後確率q(w)のモデルを作る方法として、頻繁に使われる便利な考え方
例
事後確率p(w|X,θ)において、未知パラメータw=(w1,…,w5)が5次元ベクトルであるとき
変分事後確率において、 q(w)=q(w1,w2)q(w3,w4,w5)を仮定することを、独立分解過程という
独立分解過程に含まれる特別な場合
平均場近似 (mean field approximation)がある
q(w)=q(w1)q(w2)q(w3)q(w4)q(w5)のようにスカラー変数ごとに独立な変分事後確率を考える
独立分解過程では、特定の変数ペア間の相関を無視する

5.3.2 変分ベイズ法を補助変数法に適用する

変分ベイズ法を、ガウス過程回帰の保持変数法に適用
補助変数法に出現する変数を以下のカテゴリに分解
観測される確率変数
観測値y=(y1,…,yN)T
既知定数
入力点X = (x1,…,XN)
未知定数
補助入力点Z = (z1,…,zM)
共分散関数k(x,x’:θ)を調整するパラメータ
推定対象となる確率変数
入力点における関数出力値d=(f1,…,fN)T = (f(x1),f(x2),…f(xN))
補助入力点における関数出力値u=(u1,…,uM)T = (f(z1),…,f(zM))
変分ベイズ法に当てはめる
「推定対象」となる確率変数を全部まとめたものがパラメータw=(f,u)
未知定数を全部まとめたものがハイパーパラメータθ
観測される確率変数がそのままyに対応
既知定数は定数なので無視
仮定5.8 変分ベイズ法による補助変数法のための独立分解仮定

5.3.3 ミニバッチと確率的勾配法

変分ベイズ法と確率的勾配法が威力を発揮するのは、データ点数Nが巨大な時 (具体的には数万から数億個のデータ点が存在している場合)
確率的勾配法 (stochastic gradient descent)
N個の入出力ペア(xn,yn),n=1,…Nを N7=100,500個程度の小口サイズに分ける
小単位の一つ一つをミニバッチと呼ぶ
ミニバッチを用いてKMNを禁じ計算して購買方向を求める
全データを用いて計算したものより少し誤差のある勾配が決まる
ミニバッチごとに誤差の生じる方向は異なるが、全てのミニバッチで平均すると誤差の平均(バイアス)は無視できる
ステップごとに異なるミニバッチによって勾配を求めて、これによりパラメータ更新を進める
通常の勾配法と比較して、収束までに必要な時間が少なくて済む
ミニバッチごとに勾配法こうが確率的に揺らぐので、局所解にトラップされにくい

5.4 格子状補助入力点配置にもとづくガウス過程法計算

はじめに
補助入力点を規則正しく格子状に配置することで、ガウス過程法の計算コストを大幅に低減する
多次元の場合は補助入力点を多数配置しないと制度が得られない
格子状に補助入力点を配置することで計算を効率化する

5.4.1 クロネッカー法

クロネッカー法 (Kronecker method)
多次元格子が一次元格子の直積で表されることと
行列のクロネッカー積の性質を利用して
カーネル行列に関わる演算量を低減する
例
入力点xが、2次元空間上の格子点x=(x(1),x(2))∈𝓧1x𝓧2から与えられている
𝓧1, 𝓧2はそれぞれ縦・横孔子の座標に対する50要素、70要素の集合
xは50×70の格子点上に配置されている
カーネル関数が格子次元ごとのカーネル関数の積の形になっていることが必要
2次元ならk(xn,xn’)=k(1)(xn(1),xn'(1))k(2)(xn(2),xn'(2))のように書くことができる場合
ガウスカーネルならば自動的にこの性質を満たす
全補助入力点xm, m=1,…,Mに関する共分散行列KMM(この例では3500×3500行列)は
格子次元ごとに構成した共分散行列K1(この例では50×50行列)と K2(この例では70×70行列)のクロネッカー積KMM= K1⊗K2で表すことができる
クロネッカー積 (Kronecker product)
任意形状行列(正方行列に限らない)同士で定義される行列積
例
2×3行列A=(aij)と2×2行列B=(bij)のクロネッカー積は
メモリ量が素では3500×3500=1225万
これをK1,K2に分けて保存すると 50×50 + 70×70 = 7400程度で十分
クロネッカー法の肝
クロネッカー積に展開されたカーネル行列の固有値分解
NxN行列の固有値分解K=P𝚪PTには、O(N2)のメモリ消費、O(N3)の演算量
公式5.9 クロネッカー積の固有値

5.4.2 テプリッツ法

テプリッツ法 (Toeplitz method)
一次元格子において補助入力点が等間隔に並んでおり、なおかつカーネル関数が2つの入力値xiとxjの差分で書ける形 k(xi,xj) = h(xi – xj)になっている場合に使える手法
テプリッツ法とクロネッカー法との組み合わせで、多次元の等間隔格子を取り扱うことができる
例
x=0, 0.1,0.2,…,10のような等間隔格子とガウスカーネルを使うような場合
テプリッツ行列 (Toeplitz matrix)
正方行列A(aij)であって aij = ai-1,j-1が成り立つような行列
対角一定行列 (diagonal-constant matrix)
第一行目の行ベクトルと第一列目の列ベクトルの成分が決まれば、その他全ての成分がわかる
巡回行列 (circulant matrix)
テプリッツ行列の特別な場合
MxM巡回行列とは、正方行列B=(bij)であって bij=bi-1,j-1およびbMj=b1,j+1が成り立つもの
テプリッツ行列と巡回行列の例
巡回行列には、一行目の対する離散フーリエ変換によって、その行列の固有値が得られる特徴がある

5.4.3 局所的カーネル補間

局所カーネル補間 (local kernel interporation)
任意の入力点xにおけるカーネル関数の値を得るために、 x近傍の格子点を最小限だけピックアップして格子点上の値に重みつづけして保管する
一次元格子であれば近傍格子点を2点二次元なら高々4点三次元ならば高々8点の補完で任意点の値をきめる
入力点Xとテスト点x*がとりうる値の範囲を広くカバーするようにして格子上補完点を配置する必要がある

5.4.4 KISSGP法とその演算量

クロネッカー法、テプリッツ法、局所カーネル補間の3つのテクニックはそれぞれ互いに補間的な関係にある
KISS-GP法
これらを同時に使用する実装
実装ツール
ガウス過程法の様々なアルゴリズムのMATLAB実装をまとめたGPMLというツールボックス
PythonのGPyTorchというツールボックス

第６章ガウス過程の適用

概要
ガウス過程の特徴を有効に使った実用例として、空間統計学、ベイズ最適化の課題をみる

6.1 クリギングと空間統計学

限られた回数のボーリング検査の結果から採石場全体に含まれる鉱物の総量を求める方法
ボーリング検査
ある場所のどのくらいの深さでどんな品質の鉱石が得られるかがわかる
クリーグの手法: クリンギング(Kringing)
座標x1,…,xNにおいて実際に計測された値f1,…,fNの線形総和を用いて
任意の空間座標xにおいて計測されるであろう値f(x)を常識のように求める
ガウス過程回帰法と、原理的には同じ
バリオグラムとセミバリオグラム
共分散関数のモデリングを行うために、半分散やセミバリオグラムといった言葉を使う
N点x1,…,xNにおける計測値y1,…,yNが与えられる時
任意の2点(xn,xn’)の組み合わせに対応する半分散(semivariation)を、位置xn,xn’における観測値yn,yn’の際の大きさで定義する
Snn’ = 1/2(yn – yn’)2
半分散snn’と観測転換の距離rnn’=∥xn – xn’∥の関係を散布図の形で可視化した図
バリオグラム(variogram:半分散雲)
セミバリオグラム(semi-variogram)
半分散snn’の期待値を距離rnn’で表した関数γ(r)
経験セミバリオグラム (empirical semi-variogram)
関数形状の解釈
ナゲット(nugget)
ひと固まりのものを意味する
同一箇所の繰り返し計測時の半分散の期待値
γ(0)の値をナゲット
シル(sill)
距離が十分に離れたときの半分散の値
γ(r),r >> 0の値をシル
レンジ(range)
距離がこれ以上離れても相関の大きさが変わらないと言える最小範囲
ナゲットが小さくてレンジとシルが大きい時
目的関数は狭い領域のなかでは均一性が高いが、広い範囲で見ると変化が大きい
広い範囲でボーリング作業を行なって高品質鉱床の大きな塊を探すことに意味がある
ナゲットが大きい
狭い領域の中でばらつきが大きい
狭い範囲で少しづつ位置をずらしながら何度もボーリングすることで高品質の小さな塊を少しづつ見つけて探すべき
半分散が必ずしも距離rだけの関数で表されない、より一般的な場合のモデリング
バリオグラムγ(xn,xn’)を2点の観測いち座標xn,xn’の関数の形でモデリングする

6.2 ベイズ最適化

6.2.1 ベイズ最適化とは

ベイズ最適化 (Bayesian optimization)
少数標本と最低限の仮定に基づいて確率的な予測を行うことができるガウス過程回帰法の特徴をフル活用した応用技術
目的
未知の目的関数f(x)が最大値をとるような入力x* = arg max f(x)を求めること
ユーザーは未知目的関数の形状の手かがりを得るために、入力値X=x1,x2,…を任意に決めて実験をする
実験によりコストがかかるので、実験コストを節約しつつ、最適なx*に近づくために、毎回の実験での適切なxiを決める
例
3種類の肥料(A,B,C)をそれぞれ適切な分量x=(xA,xB,xC)だけまいて麦の収量f(x)を最大化したい
新型自動車の流線型の鼻先形状パラメータxを適切に決めて、空気抵抗f(x)を最小化したい模型を3Dプリンタで出力して、風洞実験を行なって空気抵抗を求めるのに一回あたり50万円と3日間がかかる
深層ニューラルネットワークの構造パラメータxを適切に決めて、学習の精度f(x)を最適化したいパラメータを定めて行う実験1回あたり 3時間と計算機使用料60円がかかる
ベイズ最適化以前では”実験計画法(experimental design)と呼ばれる手法があった
検討対象となる要因x1,x2…がそれぞれ離散的な値x1∈{A,B,C},x2∈{大、中、小},.. をとる場合を想定
ベイズ最適化では、 N回の実験の結果として N点分の入出力データ(x1,y1),…,(xN,yN)が得られた時ガウス過程法に基づいて関数f(x)の推定を事後確率の形で求める N+1回目の実験では、これを用いてXN+1をどこに定めれば良いか決める
関数f(x)の推定結果は、事後確率の期待値μ(x)と標準偏差σ(x)の形で求められる
真の最適値を求めるためには
期待値μ(x)が最も大きな値をとる暫定最適解の周り
実験が手薄な場所(標準偏差σ(x)が大きな場所)
μとσを適当に組み合わせて定義した獲得関数(acquisition　functon)a(x)が最大になるように XN+1を決める
公式6.1 ベイズ最適化における獲得関数の定義方法
獲得関数a(x)はガウス過程で表された関数f(x)の事後確率の期待値μ=μ(x)と標準偏差σ=σ(x)を組み合わせた形で定義される
信頼性上限関数 (upper confidence bound)は以上で定義される
推定の期待値μに対して、標準偏差σに比例した大きさのマージンを加えた獲得関数を定義することで、データ密度の低い箇所周辺を好んで探索する
改善期待度 (expected improvement)は以上で定義される
ΤはこれまでN回の観測yN=f(xN)の中で最大のyを表す
t=(μ – τ)/σは、各xにおいてμとτの差を標準偏差で正規化したもの
I( f > τ )はf>τが成り立つときに1、成り立たない時に0の値をとる指数関数
Φ(),φ()はそれぞれ標準正規分布の累積分布関数と密度関数である
UCB基準によるベイズ最適化の計算ステップが進んでいく様子
獲得関数(緑色)が最大になるようなxにおける間数値f(x)を逐次的に実験で求めて推定値を更新しながら f(x)を最大にするxに近づいていく

6.2.2＊関連度自動決定(ARD)

ベイズ最適化では、対象となる関数の性質が未知であることが前提
特に、関数の入力が多次元ベクトルである場合に、 ARD(Automatic RelevanceDetermination、関連自動決定) という仕組みを用いることが有効
ARDによれば多次元ベクトルを入力とする回帰において出力に影響を持つベクトル成分を自動的に決定できる
観測点Xi=(Xi1,…,XiD)のようなD次元ベクトルで表される場合
上式のようなARDカーネルを考える
このカーネルの性質を表すハイパーパラメータθ=(𝛈1,…,𝛈D)はD次元ベクトル
このベクトルの第d成分𝛈dは正の実数値をとり、出力yの変動と観測点ベクトル夢xの第d成分の変動との関連性(relevance)の度合いを決める
𝛈dの値が0であるほか、他と比べて極端に小さい時、ガウス過程からサンプリングされた関数f(xi)の変動は入力xiの第d成分の変動の影響を受けなくなる
ハイパーパラメータ𝛈d, d=1,…,Dを適切に決めることによって入力ベクトルの各成分が関数の変動に対して与える影響の大小を設定することができる
データに基づいてこれらのハイパーパラメータを推定することで、入力ベクトルの各成分が与える影響の代償を推定することができる
ガウス過程回帰による高次元データ解析において、上記のカーネル関数をさらに拡張した公式が一般よく使われている
公式6.2 ARDカーネル関数

6.2.3＊行列微分の公式とARDアルゴリズムの導出

ARDアルゴリズムは、標準的なARDカーネル関数kを用いるほかは、ガウス過程回帰モデルのハイパーパラメータ決定法と大きく変わらない
ARD(自動的に次元選択を行なってくれる工夫)と、 Maternカーネル(共分散関数の多様性を用意してくれる工夫)の組み合わせ
ベイズ最適化の実応用のため、ソフトウェアライブラリが各種公開されている
Pyhton
GPyOpt
GPFlow
GPyTorch
R
rBayesuanOptimization

第７章ガウス過程による教師なし学習

概要
ガウス過程を使うことにより、観測値Yが潜在変数Xから生成されているとする潜在変数モデルにおいて
XからYへの写像を非線形にすることができる
同時に問題を数学的に見通しよく定義することができる
ニューラルネットワークによる潜在変数モデルを、より数学的に考えることと等価
ガウス過程に従う潜在変数のサンプリング方法についても説明

7.1 ガウス過程潜在変数モデル (GPLVM)

概要
入力xと出力yで、yだけがありxは未知のケース
例
xをwebページのクリック履歴とすると
似たユーザーは似たページを見ることが多い
ユーザーの特徴を表す入力xが与えられているわけでは無い
n番目のデータyがyn={yn(1),yn(2),…,yn(D)}T ∈ ℝD
これがN個ある。
Y={y1,y2, …, yN}T
観測値の全体YはNxD次元行列
Ynd = yn(d)はn番目のデータのd次元目の値
Yの各次元dはそれぞれ別の意味を持つ
dごとに、近いy同士は近い値になる
dごとのN次元の縦ベクトルが、それぞれ共通の未知のN個の入力X=x1,x2,…,xNからのガウス過程回帰で生成されると仮定
y(d)の平均が0であると仮定
単純なガウス分布による誤差を考える
各次元d=1,…Dについての式
f(d) 〜 N(0, Kx)
KxはXの各要素館のカーネル行列
(N,n’)要素はk(n,n’)
yn(d) = f(d)(xn) + ε, ε 〜 N(0,σ2) (n=1,…,N)
まとめると
y(d) 〜 N(0, Kx + σ2I)
観測点yは×、◯、△で表した次元ごとにそれぞれのガウス過程に従っている
その裏に▲で表した潜在的な入力×が存在している
×が近ければ、対応するベクトル(×、◯、△)同士も近くなる
データ全体のYの確率はy(1),y(2),…,y(D)の確率の積
各xがK次元の標準ガウス分布に従うと仮定
Kはxの存在する潜在空間の次元数で、Dとは無関係に決めることができる
最適化することさえできれば、 xがベクトルである必要はない
カーネル関数で類似度を測れるものであればなんでも良い
文字間のカーネルを考えて「潜在文字列」を(例えばMCMC法などで)最適化すれば、原理的には「文字から観測値が生成されるGPLVM」なども考えることができる
XとYの同時確率は上式となる
これを最大化する潜在的なXを見つけることが、学習問題の目標
ガウス仮定潜在変数モデル (Gaussian Process Latent Variable Model, GPLVM)
例:人間のポーズをGPLVMによって非線形に2次元空間に埋め込んだ(K=2)もの
Xの事前分布を考えずに直接Xを求めることも可能だが
データによってxが極端な値になることを防ぐことができないため
事前分布をおくほうが望ましい

7.1.1 GPLVMの生成モデル

GPLVMではデータYは一般的に次のように生成されたと考える
(1) For n= 1,…,N 潜在的な入力xn 〜 N(0,Ik)を生成
(2) X=x1,…,xNから、その間のカーネル行列Kxを計算
(3) For d= 1,…,D 関数値f(d) 〜 N(0, Kx)を生成観察モデルp(y|f)に従って、観測値y(d) 〜 p(y|f(d))を生成
ガウス分布による誤差を使う場合は 2つのステップをまとめて
(3)’ For d = 1,…,D 観測値y(d) 〜 N(0, Kx + σ2I)を生成
全ての観測値y(d)=(y1(d),…,yN(d))を一度に生成

7.1.2 GPLVMの目的関数

上式からどうやってxを求めるのか?
p(Y|X)の項に注目する
見やすさのためにσ2Iの項をカーネルに含めて、 Kx+σ2Iを改めてKxとする
式の変形
計算上Kx-1=𝚲と略記するとyについて
上式の破線部は
行列A,Bについて
tr(AB)は上式であることを用いた
行列AとBについて、 tr(ATB)は行列AとBの「内積」を意味する
2つの行列が一致するときに大きな値をとる
以上より、XからYを生成する確率は
Xの情報はKxの中に埋め込まれている
左式の目的関数を最大化することは
観測データの相関行列YYTをグラム行列Kx-1が出来る限りよく近似できるように、 Kxの元になる入力X=(x1,x2,…,xN)を最適化していると考えられる
tr(AB)=tr(BA)なので
上式の変換は以上になる
多次元ガウス分布の式と同様の式で書ける

7.1.3 GPLVMの学習

XからYを生成する確率
右の対数をとると、観測データの対数尤度は
Kxに含まれる、左の式を最大化するXを求めるには、基本的には勾配法を用いる
行列の微分公式から、Kxが対称行列であることを用いると上式になる
よって、右式LのKxでの微分は
微分の連鎖則から
標準的なガウスカーネルを考えると
微分はk(xn,xn’)またはk(xn’,xn)の項だけに現れる
k(xn,xn’)(n’=1,…,N)の場合だけを考えて 2倍すれば良い
最終的な計算部分
具体例
MATLAB実装のGPLVM実装の、送油データの最初の200個の潜在座標の計算結果
主成分分析の2次元写像と比べて、 GPLVMの非線形な圧縮では、本来存在するクラスがよりよく分離されている
GPLVMの目的関数は凸ではなく多数の局所解があるため、勾配法による計算は初期値に依存する
GPLVMの計算はPyhtonのGPyなどのパッケージでも行える

7.2 ガウス過程潜在変数モデルの性質

カーネルPCAとの関係
GPLVMは、主成分分析をガウス過程を用いてカーネル化したものと捉えることもできる
カーネル法によるカーネル主成分分析(カーネルPCA)と何が違うのか?
カーネルPCAでは「潜在変数」を考えずに、代わりに観測データyを投射φ(y)によって高次元の特徴空間に飛ばし、そこで通常の主成分分析を行う
GPLVMの利点と欠点
GPLVMは確立モデルであるため、データの確率が高くなる、すなわちデータをよりよく説明するパラメータ (カーネルのパラメータなど)を観測データの確率を最大化することで選ぶことができる
欠点
最適化問題が凸ではなく、主成分分析のように大域解が容易には求まらない
xとyの間に線形モデルのような明示的な関係がなく、「xが似ていればyも似ている」というガウス過程のみによって結ばれているため、データYを表現できるXには無限の可能性がある
最適化が困難
ベイズガウス潜在変数モデル
ベイズガウス過程潜在モデル (Bayesian GPLVM)
変分ベイズ法を用いたXの事後分布と、ハイパーパラメータを学習するモデル
例
GPyに含まれるGPy.modles.BayesianGPLVM() を使った送油データの学習結果
通常の次元圧縮と違って、潜在空間の「確信度」を付与することができる
明るさが確信度
データ空間Yにおける分散

7.3＊ガウス過程潜在変数モデルの拡張

7.3.1 無限ワープ混合モデル(iWMM)
7.3.2 ガウス過程力学モデル(GPDM)

7.4＊潜在的ガウス過程のサンプリング

7.4.1 ポアソン点過程とCox過程
7.4.2 楕円スライスサンプリング

付録A 付録

A.1 行列の分割と逆行列
A.2 巡回行列の固有値が離散フーリエ変換で得られること
A.3 行列微分の計算
文献案内
ガウス過程のためのソフトウェア
あとがき