ガウス過程に基づく分類モデルの学習ノート

はじめに

• 日本語の書籍が入手できなかったので、以下の資料を参考にしました。
• Gaussian Processes for Machine Learning (http://www.gaussianprocess.org/gpml/chapters/)
• 「3 Classification」の再現コードとなります。
• ラプラス近似のみ。Expectation propagationは除く。

Linear Models for Classification

入力データ$x$に対して$y=1$となる確率を$S(x^tw)$とする。$(y=\pm 1)$
ここで$S()$は、潜在変数$f=x^tw$を確率に変換するシグモイド関数である。
$S(-f)=1-S(f)$を満たすものとする。$y=-1$となる確率は$S(-x^tw)$となる。
$w$の事前分布を$N(0,\Sigma_p)$とする。
回帰モデルの場合と同様、$w$の事後分布を求めて、予測モデルから積分消去する。
テストデータ$x_∗$に対して$y_∗=1$となる確率は以下のようになる。

\begin{align}
p(y|x,w)&=S(yf)\quad (f=x^tw)\\
\log p(w|X,y)&\sim -\frac{1}{2}w^t\Sigma_p^{-1}w+\sum\log S(y_if_i)\\
p(y_*=1|x_*,X,y)&=\int p(y_*=1|x_*,w)p(w|X,y)dw\\
&=\int S(x_*^tw)p(w|X,y)dw
\end{align}

以下のサンプルでは、$w$は二次元で定数項はなしとする。$\Sigma_p=I_2$とする。
$S(f)=\frac{1}{1+\exp(-f)}$とする(ロジスティック関数)。
予測確率を求める積分は数値積分する(二次元)。
事前分布・尤度・事後分布・予測確率の等高線図を作成。

xp = c(-5,-1,-0.5,1,1,5)
yp = c(1,-5,-0.5,0,5,4)
X = rbind(xp,yp)
y = c(-1,-1,1,-1,1,1)

my_sigmoid = function(x){
  1 / (1 + exp(-x))
}

par(mfrow=c(2,2),mar=c(4,4,2,2))

my_draw_par = function(FUN,main,w1,w2){
  levels = FUN(w1,w2) - c(1,4,9,16)/4/2
  xgrid = seq(from=-2,to=2,by=0.01)
  ygrid = seq(from=-2,to=2,by=0.01)
  value = outer(xgrid,ygrid,FUN)
  contour(xgrid,ygrid,value,main=main,levels=levels,drawlabels=FALSE)
  abline(v=0,h=0,lty=2)
  points(w1,w2,pch=4)
}
my_prior = function(w1,w2){
  -(w1^2+w2^2) / 2
}
my_draw_par(my_prior,"prior",0,0)

my_likelihood = function(w1,w2){
  colSums(log(my_sigmoid(t(X) %*% rbind(w1,w2) * y)))
}
my_likelihood_optim = function(par){
  my_likelihood(par[1],par[2])
}
ret = optim(c(0,0),my_likelihood_optim,control=list(fnscale=-1))
ret
my_draw_par(my_likelihood,"likelihood",ret$par[1],ret$par[2])

my_posterior = function(w1,w2){
  my_likelihood(w1,w2) + my_prior(w1,w2)
}
my_posterior_optim = function(par){
  my_posterior(par[1],par[2])
}
ret = optim(c(0,0),my_posterior_optim,control=list(fnscale=-1))
ret
my_draw_par(my_posterior,"posterior",ret$par[1],ret$par[2])

# Numerical integration with posterior
w1_sample = ret$par[1]+seq(from=-2,to=2,by=0.2)
w2_sample = ret$par[2]+seq(from=-2,to=2,by=0.2)
v_sample = outer(w1_sample,w2_sample,my_posterior)
v_sample = exp(v_sample-max(v_sample))
v_sample = v_sample / sum(v_sample)

my_prob = function(x,y){
  apply(rbind(x,y),2,function(xy){
    fun = function(w1,w2){my_sigmoid(xy[1]*w1+xy[2]*w2)}
    sum(outer(w1_sample,w2_sample,fun) * v_sample)
  })
}
my_contour = function(){
  levels = c(0.1,0.25,0.5,0.75,0.9)
  xgrid = seq(from=-6,to=6,by=0.05)
  ygrid = seq(from=-6,to=6,by=0.05)
  value = outer(xgrid,ygrid,my_prob)
  contour(xgrid,ygrid,value,main="probability",levels=levels)
  abline(v=0,h=0,lty=2)
  points(X[1,],X[2,],type="p",pch=ifelse(y==1,1,4))
}
my_contour()

Gaussian Process Classification

潜在変数$f$がGPに従うとすると予測モデルは以下のようになる。
ここで$K(a,b)$は、$f$の共分散行列を定めるカーネル関数である。
$K=K(X,X),K_*=K(X,x_*)$と略記する。

\begin{align}
p(f|X)&=N(0,K)\\
p(y|f)&=\prod S(y_if_i)\\
p(f|X,y)&\sim p(f|X)p(y|f)\\
p(f_*|X,x_*,f)&=N(K_*^tK^{-1}f,K(x_*,x_*)-K_*^tK^{-1}K_*)\\
p(f_*|X,y,x_*)&=\int p(f_*|X,x_*,f)p(f|X,y)df\\
p(y_*=1|X,y,x_*)&=\int S(f_*)p(f_*|X,y,x_*)df_*
\end{align}

このままでは計算が困難なため、事後分布$p(f|X,y)$を正規分布で近似する。(ラプラス近似)

The Laplace Approximation for the Binary GP Classifier

ラプラス近似では、事後分布の対数値の最大値を求めて、
これを中心とする二次のテーラー展開で事後分布を近似する。
事後分布の対数値$\Psi(f)$を最大化する$\hat{f}$をニュートン法で求める。

\begin{align}
\Psi(f)&\triangleq\log p(y|f)+\log p(f|X)\\
\nabla\Psi(f)&=\nabla\log p(y|f)-K^{-1}f=0\\
\nabla^2\Psi(f)&=\nabla^2\log p(y|f)-K^{-1}=-W-K^{-1}\\
f^{\mathrm{new}}&=f-(-W-K^{-1})^{-1}(\nabla\log p(y|f)-K^{-1}f)\\
&=(K^{-1}+W)^{-1}(Wf+\nabla\log p(y|f))\\
\Psi(f)&\simeq\Psi(\hat{f})-\frac{1}{2}(f-\hat{f})^t(K^{-1}+W)(f-\hat{f})\\
p(f|X,y)&\simeq N(\hat{f},(K^{-1}+W)^{-1})
\end{align}

$W=-\nabla^2\log p(y|\hat{f})$であり、データが互いに独立であることから、対角行列となる。
ここで$p(f_*|X,x_*,f)=N(K_*^tK^{-1}f,K(x_*,x_*)-K_*^tK^{-1}K_*)$なので、
$p(f_*|X,y,x_*)=\int p(f_*|X,x_*,f)p(f|X,y)df$もまた正規分布となる。
$f_*$の平均値と分散は以下のようになる。

\begin{align}
E[f_*|X,y,x_*]&=\iint f_*p(f_*|X,x_*,f)p(f|X,y)df_*df\\
&=\int K_*^tK^{-1}fp(f|X,y)df\\
&=K_*^tK^{-1}\hat{f}=K_*^t\nabla\log p(y|\hat{f})\\
Var(f_*|X,y,x_*)&=K(x_*,x_*)-K_*^t(K+W^{-1})^{-1}K_*
\end{align}

テストデータ$x_∗$に対して$y_∗=1$となる確率は以下のようになる。

\begin{align}
p(y_*=1|X,y,x_*)&=\int S(f_*)p(f_*|X,y,x_*)df_*\\
&=\int S(f_*)N(f_*|\mu_*,\sigma_*^2)df_*
\end{align}

周辺尤度の近似値は以下のようになる。

\begin{align}
\log p(y|X)&=-\frac{1}{2}\hat{f}^tK^{-1}\hat{f}+\log p(y|\hat{f})-\frac{1}{2}\log|I+WK|
\end{align}

計算を安定化させるために行列を変形。

\begin{align}
(K^{-1}+W)^{-1}&=K-KW^{\frac{1}{2}}B^{-1}W^{\frac{1}{2}}K\\
B&=I+W^{\frac{1}{2}}KW^{\frac{1}{2}}\\
(K+W^{-1})^{-1}&=W^{\frac{1}{2}}B^{-1}W^{\frac{1}{2}}\\
|I+WK|&=|B|
\end{align}

A Toy Problem & Laplace Approximation

カーネル関数$K(a,b)=\sigma_k^2\exp\left(-\frac{|a-b|^2}{2l^2}\right)$
シグモイド関数$S(f)=\frac{1}{1+\exp(-f)}$
予測確率を求める積分は数値積分する(一次元)。

サンプルデータ

xp = c(0.18,0.41,0.47,0.57,0.64,0.65,0.67,0.78,0.86,0.89,0.11,0.13,0.14,0.19,0.23,0.28,0.36,0.41,0.46,0.79)
yp = c(0.26,0.63,0.15,0.78,0.67,0.53,0.38,0.8 ,0.6 ,0.79,0.88,0.12,0.42,0.62,0.76,0.5 ,0.28,0.45,0.88,0.71)
X = rbind(xp,yp)
y = c(-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,1,1,1,1,1,1,1,1,1,1)

カーネル関数とシグモイド関数

my_kernel = function(X1,X2,k_sigma,k_scale){
  apply(X2,2,function(a){
    k_sigma^2*exp(colSums(-(X1-a)^2)/2/k_scale^2)
  })
}
my_kernel_var = function(X,k_sigma,k_scale){
  rep(k_sigma^2,ncol(X))
}
my_sigmoid = function(y,f){
  1 / (1 + exp(-y * f))
}
my_sigmoid_log_1st = function(y,f,prob){
  y * (1-prob)
}
my_sigmoid_log_2nd = function(y,f,prob){
  prob * (1-prob)
}

ラプラス近似

my_laplace_calc = function(k_scale,k_sigma){
  K11 = my_kernel(X,X,k_sigma,k_scale)
  n = length(y)
  f = seq(from=0,to=1,length.out=n)
  while(TRUE){
    prob = my_sigmoid(y,f)
    V1 = my_sigmoid_log_1st(y,f,prob)
    V2 = my_sigmoid_log_2nd(y,f,prob)
    V2q = sqrt(pmax(V2,0))
    B = diag(n) + t(K11 * V2q) * V2q
    L = chol(B) # t(L) %*% L
    b = V2 * f + V1
    a = b - backsolve(L,forwardsolve(t(L),(K11 %*% b) * V2q)) * V2q
    f0 = f
    f = drop(K11 %*% a)
    if(max(abs(f - f0))<1e-5)break
  }
  list(f=f,prob=prob,V1=V1,V2=V2,V2q=V2q,L=L,a=a)
}

予測確率の等高線図

カーネル関数のパラメータ$\sigma_k,l$を変化させて予測確率の等高線図を作成。

par(mfrow=c(2,2),mar=c(4,4,2,2),oma=c(0,0,0,0))

plot(X[1,],X[2,],type="p",xlab="",ylab="",pch=ifelse(y==1,1,20),xlim=c(0,1),ylim=c(0,1))

my_contour = function(k_sigma,k_scale){
  calc = my_laplace_calc(k_scale,k_sigma)
  evidence = -sum(calc$a * calc$f) / 2 + sum(log(calc$prob)) - sum(log(diag(calc$L)))
  print(evidence)
  my_prob = function(x,y){
    xy = rbind(x,y)
    K21 = my_kernel(xy,X,k_sigma,k_scale)
    K22 = my_kernel_var(xy,k_sigma,k_scale)
    v = forwardsolve(t(calc$L),t(K21) * calc$V2q)
    f_mean = drop(K21 %*% calc$V1)
    f_var = K22 - colSums(v^2)
    # Numerical integration
    u_sample = seq(from=-4,to=4,by=0.1)
    f_sample = f_mean + sqrt(f_var) %*% t(u_sample)
    colSums(t(my_sigmoid(1,f_sample)) * dnorm(u_sample)) * 0.1
  }
  levels = c(0.25,0.5,0.75)
  xgrid = seq(from=0,to=1,by=0.01)
  ygrid = seq(from=0,to=1,by=0.01)
  value = outer(xgrid,ygrid,my_prob)
  main = paste(sprintf("%.3g",c(k_sigma,k_scale)),collapse=",")
  contour(xgrid,ygrid,value,levels=levels,main=main)
  points(X[1,],X[2,],type="p",pch=ifelse(y==1,1,20))
  legend("bottomright",legend=sprintf("evidence=%.3g",evidence),bty="n")
}
my_contour(3,0.1)
my_contour(3,0.2)
my_contour(3,0.3)

Binary Handwritten Digit Classification Example

カーネル関数$K(a,b)=\sigma_k^2\exp\left(-\frac{|a-b|^2}{2l^2}\right)$
シグモイド関数$S(f)=\Phi(f)$
シグモイド関数としてプロビット関数を使用する場合、予測確率を解析的に計算できる。
テストデータ$x_∗$に対して$y_∗=1$となる確率は以下のようになる。

\begin{align}
p(y_*=1|X,y,x_*)&=\int \Phi(f_*)N(f_*|\mu_*,\sigma_*^2)df_*
=\Phi\left(\frac{\mu_*}{\sqrt{1+\sigma_*^2}}\right)
\end{align}

サンプルデータの入手

USPS handwritten digit data
http://www.gaussianprocess.org/gpml/data/

library("R.matlab")
mat_data <- readMat('usps_resampled.mat')

index3 = which(mat_data$train.labels[3+1,] == 1)
index5 = which(mat_data$train.labels[5+1,] == 1)
pattern3 = mat_data$train.patterns[,index3]
pattern5 = mat_data$train.patterns[,index5]
y_train = c(rep(1,length(index3)),rep(-1,length(index5)))
x_train = cbind(pattern3,pattern5)

index3 = which(mat_data$test.labels[3+1,] == 1)
index5 = which(mat_data$test.labels[5+1,] == 1)
pattern3 = mat_data$test.patterns[,index3]
pattern5 = mat_data$test.patterns[,index5]
y_test = c(rep(1,length(index3)),rep(-1,length(index5)))
x_test = cbind(pattern3,pattern5)

カーネル関数とシグモイド関数

my_kernel = function(X1,X2,k_sigma,k_scale){
  apply(X2,2,function(a){
    k_sigma^2*exp(colSums(-(X1-a)^2)/2/k_scale^2)
  })
}
my_kernel_var = function(X,k_sigma,k_scale){
  rep(k_sigma^2,ncol(X))
}
my_sigmoid = function(y,f){
  pnorm(y * f)
}
my_sigmoid_log_1st = function(y,f,prob){
  y * dnorm(f) / prob
}
my_sigmoid_log_2nd = function(y,f,prob){
  a = dnorm(f) / prob
  a * (a + y * f)
}

ラプラス近似と周辺尤度

my_laplace_calc = function(k_scale,k_sigma){
  K11 = my_kernel(x_train,x_train,k_sigma,k_scale)
  n = length(y_train)
  f = seq(from=0,to=1,length.out=n)
  while(TRUE){
    prob = my_sigmoid(y_train,f)
    V1 = my_sigmoid_log_1st(y_train,f,prob)
    V2 = my_sigmoid_log_2nd(y_train,f,prob)
    V2q = sqrt(pmax(V2,0))
    B = diag(n) + t(K11 * V2q) * V2q
    L = chol(B) # t(L) %*% L
    b = V2 * f + V1
    a = b - backsolve(L,forwardsolve(t(L),(K11 %*% b) * V2q)) * V2q
    f0 = f
    f = drop(K11 %*% a)
    if(max(abs(f - f0))<1e-5)break
  }
  list(f=f,prob=prob,V1=V1,V2=V2,V2q=V2q,L=L,a=a)
}
my_evidence = function(k_scale,k_sigma){
  calc = my_laplace_calc(k_scale,k_sigma)
  lml = -sum(calc$a * calc$f) / 2 + sum(log(calc$prob)) - sum(log(diag(calc$L)))
  K21 = my_kernel(x_test,x_train,k_sigma,k_scale)
  K22 = my_kernel_var(x_test,k_sigma,k_scale)
  v = forwardsolve(t(calc$L),t(K21) * calc$V2q)
  f_mean = drop(K21 %*% calc$V1)
  f_var = K22 - colSums(v^2)
  info = -mean(pnorm(y_test*f_mean/sqrt(1+f_var),log=T))
  info0 = -sum(log(table(y_train)/length(y_train)) * (table(y_test)/length(y_test)))
  info = (info0 - info)/log(2)
  error = sum(ifelse((y_test*f_mean)>0,0,1))
  c(lml,info,error)
}

Marginal likelihood Maximization

周辺尤度を最大化するパラメータを算出。資料に記載されている値とは少し差がある。

my_evidence_optim = function(par){
  k_scale = exp(par[1])
  k_sigma = exp(par[2])
  K11 = my_kernel(x_train,x_train,k_sigma,k_scale)
  n = length(y_train)
  f = seq(from=0,to=1,length.out=n)
  while(TRUE){
    prob = my_sigmoid(y_train,f)
    V1 = my_sigmoid_log_1st(y_train,f,prob)
    V2 = my_sigmoid_log_2nd(y_train,f,prob)
    V2q = sqrt(pmax(V2,0))
    B = diag(n) + t(K11 * V2q) * V2q
    L = chol(B) # t(L) %*% L
    b = V2 * f + V1
    a = b - backsolve(L,forwardsolve(t(L),(K11 %*% b) * V2q)) * V2q
    f0 = f
    f = drop(K11 %*% a)
    if(max(abs(f - f0))<1e-5)break
  }
  -sum(a * f) / 2 + sum(log(prob)) - sum(log(diag(L)))
}
ret_optim = optim(c(2.85,2.35),my_evidence_optim,control=list(fnscale=-1))

> ret_optim
$par
[1] 2.755177 2.277937

> my_evidence_optim(ret_optim$par)
[1] -98.86856

> my_evidence_optim(c(2.85,2.35))
[1] -98.96438

Marginal likelihood and Information

並列処理パッケージを使用。計算には7分くらいかかる。

library(parallel)

cores = detectCores()
cl = makeCluster(cores,type="PSOCK")

clusterExport(cl,c("x_train","y_train","x_test","y_test"))
clusterExport(cl,c("my_sigmoid","my_sigmoid_log_1st","my_sigmoid_log_2nd"))
clusterExport(cl,c("my_kernel","my_kernel_var","my_laplace_calc","my_evidence"))

par(mfrow=c(1,2),mar=c(4,4,2,2))
my_contour = function(){
  xgrid = seq(from=1.125,to=5.125,by=0.25)
  ygrid = seq(from=-0.375,to=5.125,by=0.25)
  xygrid = as.matrix(expand.grid(xgrid,ygrid))
  ret = parApply(cl,xygrid,1,function(xy){my_evidence(exp(xy[1]),exp(xy[2]))})
  lml = matrix(ret[1,],nrow=length(xgrid))
  info = matrix(ret[2,],nrow=length(xgrid))
  error = matrix(ret[3,],nrow=length(xgrid))

  par(mfrow=c(2,2),mar=c(4,4,2,2))
  levels = c(-100,-105,-115,-130,-150,-200)
  contour(xgrid,ygrid,lml,xlab="log(l)",ylab="log(sf)",levels=levels,main="Log marginal likelihood")
  abline(v=c(1,2,3,4,5),h=c(0,1,2,3,4,5),lty=2)
  points(x=2.85,y=2.35,pch=3)
  levels = c(0.84,0.8,0.7,0.5,0.25)
  contour(xgrid,ygrid,info,xlab="log(l)",ylab="log(sf)",levels=levels,main="Information about test targets in bits")
  abline(v=c(1,2,3,4,5),h=c(0,1,2,3,4,5),lty=2)
  plot(x=0,xlim=c(min(xgrid),max(xgrid)),ylim=c(min(ygrid),max(ygrid)),type="n",
    xlab="log(l)",ylab="log(sf)",main="Number of test misclassifications")
  for(i in seq(length(xgrid))){
    for(j in seq(length(ygrid))){
      if(error[i,j]>=100)next
      text(x=xgrid[i],y=ygrid[j],labels=error[i,j],cex=0.8)
    }
  }
}
my_contour()

stopCluster(cl)

Histogram of the latent means

par(mfcol=c(2,1),mar=c(4,4,2,2))
k_scale = exp(2.85)
k_sigma = exp(2.35)
ret_calc = my_laplace_calc(k_scale,k_sigma)
hist(ret_calc$f,breaks=seq(from=-9.3,to=7.5,by=0.2),xlab="latent means,f",main="Training set latent means")

K21 = my_kernel(x_test,x_train,k_sigma,k_scale)
f_mean = drop(K21 %*% ret_calc$V1)
hist(f_mean,breaks=seq(from=-9.3,to=7.6,by=0.2),xlab="latent means,f",main="Test set latent means")