ベイズ最適化：ガウシアンプロセス - グラフ機械学習と強化学習について

創薬を行っていくとき、定量的構造活性相関（QSAR）を用いればハイスループットスクリーニングが可能となりますが、本当に欲しいのは所望の特性を満たす化合物です。目標となる活性値や物性値を得るために用いられる方法は、

数理最適化
ベイズ最適化
強化学習

などが考えられます。QSARモデルが得られてしまえば、逆方向の計算をするだけです。例えば、線形回帰で活性相関のあるモデルが得られれば、あとは何らかの数理最適化アルゴリズムを用いればよいです。勾配法, 遺伝的アルゴリズム, 進化計算, PSO, 焼きなまし法など、何でもよいです。決定論的なアルゴリズムだと最適解（局所解）が求まります。

ベイズ的なアプローチでは分布が求まりますので、様々な候補を得ることができます。また、獲得関数を計算すれば次の実験を行うための指針にもなります。ここでは、ベイズ最適化で良く用いられるガウシアンプロセス、ガウス過程についての理論を記載します。また、MCMCによりベイズモデルを求め、ベイズ最適化をするのも良いです。ブートストラップを用いた機械学習モデルも使用することができます。

詳細なサンプルコードはgithubにあります。PRMLの通りに実装しています。

PRML/06_Kernel_Methods at master · Masatsugar/PRML · GitHub

ガウス過程回帰

カーネル関数

$y$ と $x$ の関係を $y=f(x, w)$ とモデル化し $w$ に対する事前分布を定めるという方法でベイズ的な分析が行えます。

しかし、パラメータ $w$ を導入せず、事前分布を直接以下のように定めることができます。

\begin{align} y \sim \mathcal{N}(0, K) \end{align}

ここで、Kは

\begin{align} K_{ij} = k(x_i, x_j) = \frac{1}{\alpha} \Psi(x_i)^{T} \Psi(x_j) \end{align}

で定まるグラム行列です。このとき

\begin{align} k(x, x') = \mathbb{E}[y(x)y(x')] \end{align}

であるので、カーネル関数 $k(x, x')$ を共分散関数と解釈することができます。

例 Gaussian kernelとExponential kernel

\begin{align} k(x, x') &= \exp(\frac{- |x-x'|^{2}}{2 σ^{2}}) \\ k(x,x') &= \exp(-\theta |x-x'|) \end{align}

ソースコード

import numpy as np
import matplotlib.pyplot as plt

def gaussian_kernel(x1, x2):
    return np.exp(-(x1 - x2)**2/2*SIGMA**2)

def exponential_kernel(x1, x2):
    return np.exp(-THETA*abs(x1-x2))

def Gram(kern, x):
    N = x.shape[0]
    K = np.zeros((N, N))
    for i in range(N):
        for j in range(N):
            K[i, j] = kern(x[i], x[j])
    return K

SIGMA = 1.
THETA = 1.
ALPHA = 1.

N = 100
x = np.linspace(-1, 1, N)
kernel = gaussian_kernel
K = Gram(kernel, x)
for i in range(5):
    y = np.random.multivariate_normal(np.zeros(N), K/ALPHA)
    plt.plot(x, y)
plt.xlim(-1, 1)
plt.show()

結果

f:id:udnp:20190330140406p:plain — gaussian kernel & exponential kernel

なお、カーネル法では観測データ $\mathbf{t}$ とすると、レプリゼンタ定理より

\begin{align} y(\mathbf{x}) &= \mathbf{w}^{T} \phi(\mathbf{x}) \\ &=(\mathbf{K} + \lambda \mathbf{I})^{-1} \mathbf{t} \end{align}

という回帰式が得られます。ガウス過程回帰では、これをベイズ的アプローチを加えます。

ガウス過程回帰

観測データの組 $(\mathbf{x}_i, \mathbf{t}_i)$ が

\begin{align} t_i = y_i + \varepsilon_i \end{align} と表されるとします。 $\varepsilon_i$ はノイズを表す確率変数です。 $\varepsilon_i$ は正規分布に従うとします。

\begin{align} p(t_i|y_i) = \mathcal{N}(t_i|y_i, \beta^{-1}) \end{align}

とおきます。 $\beta$ は精度パラメータです。個々のノイズがi.i.dで生じるならば

\begin{align} p(t|y) = \mathcal{N}(t|y, \beta^{-1}\mathbf{I}) \end{align}

と書けるので、さらに、 $y$ をガウス過程としてモデル化し、

\begin{align} p(\mathbf{y}) = \mathcal{N}(0, \mathbf{K}) \end{align}

とします。このとき、観測値の分布は

\begin{align} p(\mathbf{t}) &= \int p(\mathbf{t}|\mathbf{y})p(\mathbf{y})d\mathbf{y} \\ &=\mathcal{N}(\mathbf{t}|0, \mathbf{C}) \hspace{10pt} (\mathbf{C}=\mathbf{K}+\beta^{-1}\mathbf{I}) \end{align}

です。このCの逆行列を求めればよいです。

正規分布の周辺化の公式 \begin{align} p(\mathbf{x})&=\mathcal{N}(\mathbf{x}|\boldsymbol{\mu}, \boldsymbol{\Lambda}^{-1}) \\ p(\mathbf{y}|\mathbf{x})&=\mathcal{N}(\mathbf{y}|\mathbf{A}\mathbf{x}+\mathbf{b}, \mathbf{L}^{-1}) \end{align}

このとき、 $p(y)$ は下式です。 \begin{align} p(\mathbf{y})=\mathcal{N}(\mathbf{y}|\mathbf{A}\boldsymbol{\mu}+\mathbf{b}, \mathbf{L}^{-1}+\mathbf{A}\boldsymbol{\Lambda}^{-1}\mathbf{A}^T) \end{align}

さて、目的は学習データ ${(\mathbf{x}_1, t_1), ...,(\mathbf{x}_n, t_n) }$ が与えられているときに、新しい $\mathbf{x}_{n+1}$ に対する $t_{n+1}$ を予測することです。つまり、 $\mathbf{t}=(t_1,...,t_n)^{T}$ が与えられたとき $t_{n+1}$ の分布を求めればよいです。

\begin{align} p(t_{n+1}|\mathbf{t}) \end{align}

まず、先ほどと同様にして、

\begin{align} p(\mathbf{t}, t_{n+1})=\mathcal{N}(\mathbf{t}, t_{n+1}|\mathbf{0}, \mathbf{C'}) \end{align}

となります。ここで、

\begin{align} \mathbf{C'}&=\begin{pmatrix} \mathbf{C} & \mathbf{k} \\ \mathbf{k}^T & c \end{pmatrix} \\ \mathbf{k} &= (k(\mathbf{x}_1, \mathbf{x}_{n+1}), ...,\mathbf{x}_n, \mathbf{x}_{n+1}), ) \\ c &= k(\mathbf{x}_{n+1}, \mathbf{x}_{n+1}) + \beta^{-1} \end{align}

なので、もとめる $t_{n+1}$ の分布は

\begin{align} p(t_{n+1} | \mathbf{t} ) = \mathcal{N} (t_{n+1} | \mu, σ^{2}) \end{align}

です。ただし、

\begin{align} \mu &=\mathbf{t}^{T} \mathbf{C}^{-1} \mathbf{k} \\ σ^{2} &= c - \mathbf{k}^{T} \mathbf{C}^{-1} \mathbf{k} \end{align}

となります。従って、 $\mathbf{a}=\mathbf{C}^{-1}\mathbf{t}$ とおけば、

\begin{align} y = \sum_{i=1}^n a_i k(\mathbf{x}_i,\mathbf{x}) \end{align}

となります。これはカーネル法の式と同じです。

Hyper parameterの学習

PRMLの本でよく用いられるカーネルは以下の式と記載されています。 \begin{align} k(x, x') = \theta_0 \exp(-\theta_1 |x-x'|^{2} / 2 ) + \theta_2 + \theta_3 xx' \end{align}

$\boldsymbol{\theta}=(\theta_0,\theta_1,\theta_2,\theta_3)$ は事前分布の設定に関わるハイパーパラメータです。これを最尤法を使って求めましょう。第二種の最尤推定、経験ベイズとも呼ばれます。

\begin{align} L(\boldsymbol{\theta}|\mathbf{t}) = p(\mathbf{t}|\boldsymbol{\theta}) \end{align}

を最大化します。

\begin{align} p(\mathbf{t}|\boldsymbol{\theta}) &= \mathcal{N}(\mathbf{t}|\mathbf{0}, \mathbf{C}) \hspace{10pt} (\mathbf{C}=\mathbf{K}+\beta^{-1}\mathbf{I}) \\ &=\frac{1}{(2 \pi)^{n} \sqrt{|\mathbf{C}|}} \exp \left\{ -\frac{1}{2} \mathbf{t}^{T} \mathbf{C}^{-1} \mathbf{t} \right\} \end{align}

であるので、対数尤度は

\begin{align} \ln L(\boldsymbol{\theta}|\mathbf{t})=-\frac{1}{2}\ln |\mathbf{C}|-\frac{1}{2}\mathbf{t}^T\mathbf{C}^{-1}\mathbf{t}-\frac{n}{2}\ln(2\pi) \end{align}

です。これを最大化したいので、勾配を計算します。

\begin{align} \frac{\partial}{\partial \theta_i}\ln L(\boldsymbol{\theta}|\mathbf{t})=-\frac{1}{2}\text{Tr}\bigg(\mathbf{C}^{-1}\frac{\partial \mathbf{C}}{\partial \theta_i}\bigg)-\frac{1}{2}\mathbf{t}^T\mathbf{C}^{-1}\frac{\partial \mathbf{C}}{\partial \theta_i}\mathbf{C}^{-1}\mathbf{t} \end{align}

【解説】 \begin{align} \mathbf{C}\mathbf{C}^{-1} = \mathbf{I} \end{align}

の両辺を $\theta_i$ で微分すればよい。

\begin{align} \frac{\partial \mathbf{C}}{\partial \theta_i}\mathbf{C}^{-1} + \mathbf{C} \frac{\partial \mathbf{C}^{-1}}{\partial \theta_i}= \mathbf{O} \end{align}

ゆえに、

\begin{align} \frac{\partial \mathbf{C}^{-1}}{\partial \theta_i}=\mathbf{C}^{-1}\frac{\partial \mathbf{C}}{\partial \theta_i}\mathbf{C}^{-1} \end{align}

また、 $\mathbf{C}$ の固有値・固有ベクトルを $(\lambda_i, \mathbf{u}_i) \hspace{6pt}(i=1,...,n)$ とする。 $\{\mathbf{u}_i \}$ は正規直交基底をなすようにとる。 $\mathbf{C}$ は対象行列なので、

\begin{align} \mathbf{C} = \sum_i \lambda_i \mathbf{u}_i \mathbf{u}_{i}^{T} \end{align}

と表します。これを $\theta$ で微分すると

\begin{align} \frac{\partial \mathbf{C}}{\partial \theta} = \sum_i \frac{\partial \lambda_i}{\partial \theta} \mathbf{u}_i\mathbf{u}_i^{T} + \sum_i 2\lambda_i \mathbf{u}_i \frac{\partial \mathbf{u}_i^{T}}{\partial \theta} \end{align}

となる。これと、

\begin{align} \mathbf{C}^{-1} = \sum_i \frac{1}{\lambda_i} \mathbf{u}_i \mathbf{u}_i^{T} \end{align}

を掛けて、 $\mathbf{u}_i \mathbf{u}_i^{T}=\delta_{ij}$ を用いると、

\begin{align} \mathbf{C}^{-1}\frac{\partial \mathbf{C}}{\partial \theta} = \sum_i \frac{1}{\lambda_i} \frac{\partial \lambda_i}{\partial \theta} \mathbf{u}_i\mathbf{u}_i^{T} + \sum_i 2 \mathbf{u}_i\frac{\partial \mathbf{u}_i^{T}}{\partial \theta} \end{align}

となります。 $\text{Tr}(a b^{T})=a^{T} b$ より

\begin{align} \text{Tr}\bigg(\mathbf{C}^{-1}\frac{\partial \mathbf{C}}{\partial \theta}\bigg) = \sum_i \frac{1}{\lambda_i} \frac{\partial \lambda_i}{\partial \theta} \mathbf{u}_i^{T}\mathbf{u}_i + \sum_i 2 \mathbf{u}_i^{T} \frac{\partial \mathbf{u}_i}{\partial \theta} \end{align}

となり、 $||\mathbf{u}_i||^2=1$ とこの両辺を微分して得られる

\begin{align} 2u_i^{T} \frac{\partial u_i}{\partial \theta} = 0 \end{align}

を代入することで、

\begin{align} \text{Tr}\bigg(\mathbf{C}^{-1}\frac{\partial \mathbf{C}}{\partial \theta}\bigg) = \sum_i \frac{1}{\lambda_i} \frac{\partial \lambda_i}{\partial \theta} \end{align}

となります。ここで、 $|C|=\prod_i \lambda_i$ であるから、

\begin{align} \frac{\partial}{\partial \theta} \ln |C| &= \frac{\partial}{\partial \theta} \sum_i \ln \lambda_i \ &= \sum_i \frac{1}{\lambda_i}\frac{\partial \lambda_i}{\partial \theta} \end{align}

となります。以上より、

\begin{align} \frac{\partial}{\partial \theta} \ln |C| =\text{Tr} \bigg( \mathbf{C}^{-1} \frac{\partial \mathbf{C}}{\partial \theta} \bigg) \end{align}

です。本来はこのように直接微分して求めるのではなく、計算機上で近似して求めます。

GPソースコード

上記の一連の式を記載すればガウス過程回帰が行えます。

import numpy as np
import matplotlib.pyplot as plt
from scipy import linalg as LA


class GaussianProcess:
    def __init__(self, kernel=None, alpha=1e-10, beta=2):
        self.kernel = kernel
        self.alpha = alpha
        self.beta = beta
        self.train_x = None
        self.train_t = None
        self.K = None
        self.Cinv = None
        self.a = None
        self.N = None

        if kernel is None:
            SIGMA = 1.0
            self.kernel = lambda x1, x2: np.exp(-(x1 - x2) ** 2 / (2 * SIGMA ** 2))

    def fit(self, train_x, train_t):
        # calculate Gram matrix
        N = train_x.shape[0]
        K = np.zeros((N, N))
        for i in range(N):
            for j in range(N):
                K[i, j] = self.kernel(train_x[i], train_x[j])

        C = K + np.identity(N) / self.beta
        Cinv = LA.inv(C)
        a = Cinv.dot(train_t)

        self.train_x = train_x
        self.train_t = train_t
        self.K = K
        self.Cinv = Cinv
        self.a = a
        self.N = N

    # m (x_{N+1}) = \sum_{n=1}^N a_n k(x_n, x_{N+1})
    def predict(self, new_x):
        y = 0
        for i in range(self.N):
            y += self.a[i] * self.kernel(self.train_x[i], new_x)
        return y

    def sigma(self, new_x):
        vec_k = np.array([self.kernel(train_x[i], new_x) for i in range(self.N)])
        kt_Cinv_k = vec_k.T.dot(self.Cinv).dot(vec_k)
        c = self.kernel(new_x, new_x)
        return np.sqrt(c + 1.0/self.beta - kt_Cinv_k)

    # p(t_{n+1} | t) = N( t | mu(x), sigma(x))
    def cond_p(self, new_x, train_t):
        mu = self.predict(new_x)
        sig = self.sigma(new_x)
        return 1/np.sqrt(2*np.pi*sig**2) * np.exp(-(train_t - mu)**2/(2*sig**2))

    # p(t) = \int p(t|y)p(y)dy = N(0, C)
    # K = phi(x_n) phi(x_m)^t
    # C = K/alpha
    def p_t(self):
        for i in range(10):
            y = np.random.multinomial(np.zeros(self.N), self.K/self.alpha)
        return y

    
def rbf_sine_kernel(x1, x2):
    return sum([np.sin(i*np.pi*x1)*np.sin(i*np.pi*x2) for i in range(n_basis)])

THETA0 = 1;THETA1 = 4;THETA2 = 0;THETA3 = 5
def normal_kernel(x1, x2):
    return THETA0 * np.exp(-THETA1*(x1-x2)**2/2) + THETA2 + THETA3*x1*x2

結果

f:id:udnp:20190330134901p:plain — ガウス過程回帰

カーネルの選択が大事になってきます。ベイズ最適化を行うには、獲得関数を計算することになります。

UCB（Upper Confidence Bound）
EI (Expected Improvement
PI (Probability of Improvement)

などがあります。次回、これらをまとめていきたいと思います。