ガウス過程と機械学習という本を読み進めたので、そのメモ。
丸写しする訳ではない上、全部載せる気力もないので重要そうなところだけをかいつまんで書いていく。
詳しい導出部分や解説はぜひ、上の書籍を参考にしてみてほしい。
下のリンクから飛ぶことができる。
ガウス過程と機械学習 – サポートページ
定式化
そもそもガウス過程って
ガウス過程って名前が正直ピンと来ないというか中身が想像しにくい内容であったが、それも書籍内で解説してくれていた。
どうやらガウス過程というのは 無限次元のガウス分布 のことを指しており,
確率過程は入力に対してそれに対応する確率変数 y に同時分布 p(y) を与える確率モデルのことをさすらしい。
なので、確率モデルの一種ということであり、やりたいことは出力となる y の期待値(平均)と分散を求めるということになりそうである。
似たようなものとしてベイズ線形回帰があるが、こちらは重みとなる w の期待値と分散を求めることが目的なので若干モデルの構造が異なる。
今回も推定していくモデルは基底関数 ϕi(x) の線型結合で表現され、
y^=wTϕ(x)=w0ϕ0(x)+w1ϕ1(x)+⋯+wHϕH(x)
この基底関数と重みを直接求めることなく、 y^ の期待値と分散をモデルかしてしまうのがガウス過程回帰だ。
ベイズ線形回帰については以下のページでまとめたのできになる人は見てみてほしい。
ベイズ線形回帰をまとめてみる(Python, bayesian regression)
カーネルトリック
先ほども言及したように、ガウス過程回帰は基底関数と重みを直接求めることなく y^ をモデル化する確率モデルとなっている。
まずは y^=(y^1,y^2,⋯,y^N)T の共分散行列 K は重み w を以下のようにモデル化することでこのように表せる。
w∽N(0,λ2I)(1)
K=λ2ΦΦT(2)
(Φ=(ϕ(x1),ϕ(x2),⋯,ϕ(xN))T)(3)
このように、重み w を使わずに表現することが可能。
この時、 K の (n,n’) における要素 Knn’ の値は
Knn’=ϕ(xn)Tϕ(xn’)(4)
と表せ、基底関数ベクトル ϕ を求める必要はなく、 Knn’ を直接求めるようなカーネル関数
k(xn,xn’)
で表現すれば良いということになる。
Knn’ は ϕ(xn) と ϕ(xn’) との近似度を表している(値の傾向が近ければ近いほど共分散の値は大きくなる)ので、カーネル関数も xn と xn’ との類似度(距離)を表すような関数を設定すれば良い。
このように、カーネル関数を用いて共分散行列を表現することを カーネルトリック と呼ぶ。
カーネル関数の種類
最もよく使われるのは RBFカーネル(動径基底関数; Radial Basis Function) と呼ばれるもので、以下のように表される。
k(xn,xn’)=θ1exp(−θ2∣xn–xn’∣2)(5)
他にも 線形カーネル(linear kernel)
k(xn,xn’)=xnTxn’(6)
や 指数カーネル(exponential kernel)
k(xn,xn’)=exp(−θ∣xn–xn’∣)(7)
などを使う。
ガウス過程回帰では、これらの線型結合でカーネル関数を表現することもできる。
k(xn,xn’)=θ0exp(−θ1∣xn–xn’∣2)+θ2xnTxn’+θ3δ(n,n’)(8)
δ(n,n’) は n=n’ の時のみ1を返す関数で、これは観測値 y に期待値 0 、分散 σ2 をもつノイズが乗っていると過程した際の追加分の分散。
そのほかにも様々なタイプのカーネル関数が存在するらしい。
ガウス過程回帰モデル
今までのを一旦まとめると、N個のデータの組 D=(x1,y1),(x2,y2),⋯,(xN,yN) があるとすると
y∽N(0,K)(9)
と y を確率モデル化することができる。
未知の値への予測
では、得られたこの分布から、新たな値 x∗ の新たな予測値 y∗ の期待値と分散を求めていく。
この時、
k∗=(k(x∗,x1),k(x∗,x2),⋯,k(x∗,xN))Tk∗∗=k(x∗,x∗)
とすると、
p(y∗∣x∗,D)=N(k∗TK−1y,k∗∗–k∗TK−1k∗)(10)
これを複数の未知の入力
X∗
が与えられた時の期待値と分散は以下のようになる。
この時、
k∗,k∗∗
を以下のようにすると、


未知の観測値の分布は以下のように表される
p(y∗∣X∗,D)=N(k∗TK−1y,k∗∗–k∗TK−1k∗)(11)
ハイパーパラメータ推定
今までみてきたように、カーネル関数を定めてしまえばあとは手順に沿って計算するだけで未知の値の期待値と分散をそれぞれ求めることができてしまう。
そのため、学習するのはカーネル関数に含まれるパラメータとなる θ となる。
よって式 (8) のようなカーネル関数を考える。
これらの θ を持ったカーネル行列 Kθ を使うと、
p(y∣X,θ)=N(y∣0,Kθ)=(定数)×∣Kθ∣211exp(−21yTKθ−1y)(12)
そのため、これらの対数をとると
log p(y∣X,θ)=–21log∣Kθ∣–21yTKθ−1y+const.=L+const.(13)
このように L をおくと、 L を最大化するような θ を求めれば良いということになる。
したがって L をハイパーパラメータ θ で微分すると
∂θ∂L=–tr(Kθ−1∂θ∂Kθ)+(Kθ−1y)T∂θ∂Kθ(Kθ−1y)(14)
∂θKθ は Kθ の各要素に対してそれぞれ θ で偏微分した値を要素にもつ行列となる。
あとは、確率的勾配降下法(SGD; Stochastic Gradient Descent)を用いてパラメータを何回か更新する。
ここのパラメータ最適化の方法はいくつかあるらしいが自分の馴染みのあるものを使用することにする。
パラメータ更新率 η を使って、
θnew=θ+η∂θ∂L(15)
と各々の θ を更新していく。
実装
カーネルクラス
まずはカーネルを作成するクラスを設定する。
カーネルは以下のものを使用する。
k(xn,xn’)=θ0exp(−θ1∣xn–xn’∣2)+θ2xnTxn’+θ3δ(n,n’)
実装
カーネルクラス
まずはカーネルを作成するクラスを設定する。
カーネルは以下のものを使用する。
k(xn,xn’)=θ0exp(−θ1∣xn–xn’∣2)+θ2xnTxn’+θ3δ(n,n’)
今後を考えて、カーネルマトリックスの微分系( ∂θ∂Kθ や k∗∗ )を求めるものも用意した。
class Kernel(object):
def __init__(self, theta=[], seed=42):
np.random.seed(seed)
if len(theta) != 4:
self.theta = np.random.randn(4)
else:
self.theta = np.array(theta)
def make_kernel(self, X):
"""
inputs :
X : 2d-array. shape=[N, d].
returns :
K : 2d-array. shape=[N,N].
The matrix that represents covariance.
"""
N = X.shape[0]
K = np.zeros(shape=[N,N],dtype=float)
for i in range(N):
for j in range(N):
K[i,j] = self._kernel(X[i], X[j], i, j)
return K
def _kernel(self, x, x_, n=0, n_=1):
"""
inputs :
x : 1d-array. shape=[d]
x_ : 1d-array. shape=[d]
n, n_ : the index of data.
returns:
k : scaler.
"""
k_1 = self.theta[0] * np.exp(- np.sum((x - x_)**2) / self.theta[1])
k_2 = self.theta[2] * np.sum((x - x_)**2)
k_3 = self.theta[3] * (n==n_)
return k_1 + k_2 + k_3
def part_k(self, x, x_, n, n_):
"""
inputs:
x, x_ : 1d-array. shape=[d]
n, n_ : the index of data
returns:
dk : 1d-array. shape=[4].
the result of dk/dθ
"""
dk = np.zeros(4)
dk[0] = np.exp(- np.sum((x - x_)**2) / self.theta[1])
dk[1] = self.theta[0] * np.sum((x-x_)**2) / self.theta[1]**2 * np.exp(- np.sum((x - x_)**2) / self.theta[1])
dk[2] = np.sum((x - x_)**2)
dk[3] = 1.0 * (n==n_)
return dk
def part_K(self, X):
"""
inputs:
X : 2d-array. shape=[N, d].
returns:
dK :3d-array. shape=[4, N, N,].
the result of dK/dθ
"""
N = X.shape[0]
dK = np.zeros(shape=[N, N, 4])
for i in range(N):
for j in range(N):
dK[i, j] =self. part_k(X[i], X[j], i, j)
dK = dK.transpose(2, 0, 1)
return dK
def make_k_star(self, X, X_star):
"""
inputs:
X : 2d-array. The known values
X_star: 2d-array. shape=[M,d]. The unkonwn values.
returns:
k_star: 2d-array. shape=[N, M]
"""
N = X.shape[0]
M = X_star.shape[0]
k_star = np.zeros((N, M))
for i in range(N):
for j in range(M):
k_star[i,j] = self._kernel(X[i], X_star[j])
return k_star
ガウス過程回帰クラス
今度は本体部分のクラス。
part_l関数が ∂θ∂L を求める関数。
回帰してみる
今回は以下の関数
f(x)=exp(x−2.0)+1.2
を予測させることにした。

この結果がこれ。

分散に当たる vが場合によっては負の値になってしまっているので数式的にまずいミスを犯している気がしてならないが、一旦保留。
毎度の事ながら、実装したコードは以下のレポジトリにアップロードしてある。
https://github.com/wildgeece96/gaussian_process
「ガウス過程回帰についてまとめてみる(Python, gaussian process regression)」への1件のフィードバック
コメントは受け付けていません。