Posted at

SVMより良いという噂のGaussian Processによる回帰を実装してみる

More than 3 years have passed since last update.


概要


“Anything you can do easily with an SVM you can do with a Gaussian Process better”


Gaussian Processについて、Zoubin Ghahramani氏

ご本人もpossibly controversialと言っているし少し言い過ぎな気もしますが(何をもって"better"なのか?)、Gaussian Processがどんなものなのか理解するためにとりあえず実装してみました。


Multivariate Gaussian Distribution

Gaussian ProcessではMultivariate Gaussian Distributionを使います。

Multivariate Gaussian Distributionは正規分布の多変量版のようなもので、変数同士の共分散が定義されます。

$n$次元の確率変数がMultivariate Gaussian Distributionは、平均を表す$n$次元ベクトル$\mu$と$n$行$n$列の共分散行列$\Sigma$で定義されます。

この分布からのサンプリングをするには、Cholesky分解$\mathrm{chol}(\cdot)$と標準正規分布から独立同分布でサンプリングしたベクトル$z$を用い、

x = \mu + \mathrm{chol}(\Sigma)z

を計算します。

sample_from_multivariate_gaussian <- function(mean, covariance) {

dim <- length(mean)
covariance.cholesky <- t(chol(covariance))
mean + covariance.cholesky %*% rnorm(dim)
}

試しに2次元でサンプリングした点をプロットしてみます。

\Sigma = \begin{bmatrix}

1 & 0 \\
0 & 4
\end{bmatrix}

\Sigma = \begin{bmatrix}

1 & 0.5 \\
0.5 & 1
\end{bmatrix}

共分散の性質を再現してデータが生成されています。


Gaussian Process


定義

Gaussian Processは$n$次元ベクトル$x$の関数$f(x)$を生成する分布で、

平均関数$m(x)$と共分散関数$k(x_1, x_2)$に対して以下の性質を与えるものです。

E\left[f(x)\right] = m(x)

E\left[\left(f(x_1) - m(x_1)\right)\left(f(x_2) - m(x_2)\right)\right] = k(x_1, x_2)

$m$と$k$で定義される平均値と共分散を用いたMultivariate Gaussian Distributionによって関数値が生成されると解釈することができます。

以降は議論を簡単にするために$m(x) = 0$とします。

またこの記事では以下の共分散関数を使います。

k(x_1, x_2) = \sigma_f \exp \left( - \frac{1}{2 l^2} (x_1 - x_2)^2 \right)


回帰問題への応用

回帰問題では、入力$(x_1, \cdots, x_n)$の一部の$x_i$に対する$f(x_i)$の値が既知で、その他の点に対する関数値が未知で予測の対象と考えます。

既知の成分を$X$、未知の成分を$X^*$と表し、共分散行列を与える関数$K(X_1, X_2)$を以下のように定義します。

\left(K(X_1, X_2)\right)_{ij} = k\left( \left(X_1\right)_{i}, \left(X_2\right)_{j} \right)

参考文献[1]で導出される以下の式で、未知の成分に対するMultivariate Gaussian Distributionを得ます。

\mu = K(X, X^*)K(X, X)^{-1}f(X)

\Sigma = K(X^*, X^*) - K(X, X^*)K(X, X)^{-1}K(X, X^*)

今回は$X = (-1, -0.8, -0.6, \cdots, 1)$とし、うち3点の関数値を既知としました。

# define input space

x <- seq(-1, 1, 0.2)
num_input <- length(x)

# create k matrix
k <- matrix(nrow = num_input, ncol = num_input)
l <- 1
sigma_f <- 0.1
for (row in 1:nrow(k)) {
for (col in 1:ncol(k)) {
x_row <- x[row]
x_col <- x[col]
k[row, col] = sigma_f * exp(-0.5 / l^2 * (x_row - x_col)^2)
}
}

# fix y
known_index <- c(4, 5, 8)
known_y <- c(1, 0.5, 0.25)

# decompose K matrix
all_index <- 1:num_input
unknown_index <- setdiff(all_index, known_index)
k.known.inv <- solve(k[known_index, known_index])
k.unknown <- k[unknown_index, unknown_index]
k.cross <- k[unknown_index, known_index]
k.cross.T <- t(k.cross)
mean <- k.cross %*% k.known.inv %*% known_y
covariance <- k.unknown - k.cross %*% k.known.inv %*% k.cross.T


実行結果

前項で計算したMultivariate Gaussian Distributionとそのサンプリングのアルゴリズムを使って、未知の点をサンプリングしてみます。

y_samples <- matrix(nrow = num_samples, ncol = num_input)

for (sample in 1:num_samples) {
y_samples[sample, unknown_index] <- sample_from_multivariate_gaussian(mean, covariance)
y_samples[sample, known_index] <- known_y
}

Gaussian Processで何が起きているのかを図示してみます。

first_samples <- y_samples %>% head %>% melt

names(first_samples) <- c("sample_index", "x", "y")
ggplot(first_samples, aes(x = x, y = y, color = factor(sample_index))) + geom_path() + geom_point()

各色の線が一度のサンプルで得られた関数値で、既知の点もプロットしています。

次にサンプルの平均値と各点の95%信頼区間をプロットしてみます。

# post process

col_var <- apply(y_samples, MARGIN = 2, FUN = var)
plot_df <- data.frame(
y = colMeans(y_samples),
sd = sqrt(col_var),
index = factor(1:ncol(y_samples)),
x = x,
known = all_index %in% known_index
)

# plot error bar
dodge <- position_dodge(width = 1)
p <- ggplot(plot_df, aes(x = x, y = y, ymax = y + 3*sd, ymin = y - 3*sd))
p + geom_line() + geom_point(aes(color = known), size = 4) + geom_errorbar(position = dodge, width = 0.1)

既知の点の近くの点ほど信頼区間が狭いという直感的に解釈しやすい結果を得られました。


共分散関数のパラメータ

共分散関数の定義に使用したパラメータ$l$は「各点がどれだけ近くの点の影響を受けるか」を表していて、より感覚的には関数の「柔らかさ」のようなものに相当します。例えば$l$を$0.1$倍すると以下の結果が得られます。


まとめ

簡単な回帰問題をGaussian Processで解くことができました。

この手法は分類問題に応用可能で、例えば二値分類の場合は$f(x)$をシグモイドに与えて分類確率を得るといったことが考えられます[1]。

冒頭のSVMとの対比に対しては以下の点を思いました。


  • 計算コストはSVMのほうが有利。$K(X, X)^{-1}$の計算がネック。

  • overfittingしない点でSVMより優秀。これはベイジアンモデル全般に言える。

  • 各点の信頼度が出せるのは多くの最尤法やSVMよりも優れている

  • パラメータの調整が必要な点ではSVMに対して優劣はつけられない


参考資料


  1. Gaussian Process For Machine Learning

  2. Zoubin Ghahramani氏のBayesian Method全般の講義

  3. David MacKay氏によるGaussian Processの講義


(参考)全コード

library(ggplot2)

library(reshape2)
library(magrittr)

sample_from_multivariate_gaussian <- function(mean, covariance) {
dim <- length(mean)
covariance.cholesky <- t(chol(covariance))
mean + covariance.cholesky %*% rnorm(dim)
}

# define input space
x <- seq(-1, 1, 0.2)
num_input <- length(x)

# create consistent K matrix
k <- matrix(nrow = num_input, ncol = num_input)
l <- 1
sigma_f <- 0.1
for (row in 1:nrow(k)) {
for (col in 1:ncol(k)) {
x_row <- x[row]
x_col <- x[col]
k[row, col] = sigma_f * exp(-0.5 / l^2 * (x_row - x_col)^2)
}
}

# fix y
known_index <- c(4, 5, 8)
known_y <- c(1, 0.5, 0.25)

# decompose K matrix
all_index <- 1:num_input
unknown_index <- setdiff(all_index, known_index)
k.known.inv <- solve(k[known_index, known_index])
k.unknown <- k[unknown_index, unknown_index]
k.cross <- k[unknown_index, known_index]
k.cross.T <- t(k.cross)
mean <- k.cross %*% k.known.inv %*% known_y
covariance <- k.unknown - k.cross %*% k.known.inv %*% k.cross.T

# sample
num_samples <- 1024
y_samples <- matrix(nrow = num_samples, ncol = num_input)
for (sample in 1:num_samples) {
y_samples[sample, unknown_index] <- sample_from_multivariate_gaussian(mean, covariance)
y_samples[sample, known_index] <- known_y
}

# post process
col_var <- apply(y_samples, MARGIN = 2, FUN = var)
plot_df <- data.frame(
y = colMeans(y_samples),
sd = sqrt(col_var),
index = factor(1:ncol(y_samples)),
x = x,
known = all_index %in% known_index
)

# plot error bar
dodge <- position_dodge(width = 1)
p <- ggplot(plot_df, aes(x = x, y = y, ymax = y + 3*sd, ymin = y - 3*sd))
p + geom_line() + geom_point(aes(color = known), size = 4) + geom_errorbar(position = dodge, width = 0.1)

# plot samples
first_samples <- y_samples %>% head %>% melt
names(first_samples) <- c("sample_index", "x", "y")
ggplot(first_samples, aes(x = x, y = y, color = factor(sample_index))) + geom_path() + geom_point()