創薬を行っていくとき、定量的構造活性相関（QSAR）を用いればハイスループットスクリーニングが可能となりますが、本当に欲しいのは所望の特性を満たす化合物です。目標となる活性値や物性値を得るために用いられる方法は、
などが考えられます。QSARモデルが得られてしまえば、逆方向の計算をするだけです。例えば、線形回帰で活性相関のあるモデルが得られれば、あとは何らかの数理最適化アルゴリズムを用いればよいです。勾配法, 遺伝的アルゴリズム, 進化計算, PSO, 焼きなまし法など、何でもよいです。決定論的なアルゴリズムだと最適解（局所解）が求まります。
ベイズ的なアプローチでは分布が求まりますので、様々な候補を得ることができます。また、獲得関数を計算すれば次の実験を行うための指針にもなります。ここでは、ベイズ最適化で良く用いられるガウシアンプロセス、ガウス過程についての理論を記載します。また、MCMCによりベイズモデルを求め、ベイズ最適化をするのも良いです。ブートストラップを用いた機械学習モデルも使用することができます。
詳細なサンプルコードはgithubにあります。PRMLの通りに実装しています。
PRML/06_Kernel_Methods at master · Masatsugar/PRML · GitHub
ガウス過程回帰
カーネル関数
と
の関係を
とモデル化し
に対する事前分布を定めるという方法でベイズ的な分析が行えます。
しかし、パラメータを導入せず、事前分布を直接以下のように定めることができます。
\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}
であるので、カーネル関数を共分散関数と解釈することができます。
- 例 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()
結果
なお、カーネル法では観測データ とすると、レプリゼンタ定理より
\begin{align} y(\mathbf{x}) &= \mathbf{w}^{T} \phi(\mathbf{x}) \\ &=(\mathbf{K} + \lambda \mathbf{I})^{-1} \mathbf{t} \end{align}
という回帰式が得られます。ガウス過程回帰では、これをベイズ的アプローチを加えます。
ガウス過程回帰
観測データの組が
\begin{align}
t_i = y_i + \varepsilon_i
\end{align}
と表されるとします。はノイズを表す確率変数です。
は正規分布に従うとします。
\begin{align} p(t_i|y_i) = \mathcal{N}(t_i|y_i, \beta^{-1}) \end{align}
とおきます。は精度パラメータです。個々のノイズがi.i.dで生じるならば
\begin{align} p(t|y) = \mathcal{N}(t|y, \beta^{-1}\mathbf{I}) \end{align}
と書けるので、さらに、をガウス過程としてモデル化し、
\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}
このとき、は下式です。
\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}
さて、目的は学習データが与えられているときに、新しい
に対する
を予測することです。つまり、
が与えられたとき
の分布を求めればよいです。
\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}
なので、もとめるの分布は
\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}
となります。従って、とおけば、
\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}
は事前分布の設定に関わるハイパーパラメータです。これを最尤法を使って求めましょう。第二種の最尤推定、経験ベイズとも呼ばれます。
\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}
の両辺をで微分すればよい。
\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}
また、の固有値・固有ベクトルを
とする。
は正規直交基底をなすようにとる。
は対象行列なので、
\begin{align} \mathbf{C} = \sum_i \lambda_i \mathbf{u}_i \mathbf{u}_{i}^{T} \end{align}
と表します。これをで微分すると
\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}
を掛けて、を用いると、
\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}
となります。より
\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}
となり、とこの両辺を微分して得られる
\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}
となります。ここで、であるから、
\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
結果
カーネルの選択が大事になってきます。ベイズ最適化を行うには、獲得関数を計算することになります。
- UCB（Upper Confidence Bound）
- EI (Expected Improvement
- PI (Probability of Improvement)
などがあります。次回、これらをまとめていきたいと思います。
