Automatic Relevance Determination（ARD）

はじめに

今回は、Automatic Relevance Determination（ARD）について解説する。ARDとは、目的変数に対する個々の説明変数の寄与の大きさを見積もる手法である。今回は以下の2つの文脈でARDを行う。

ベイズ推定を用いた回帰によりARDを行う。

ガウス過程を用いた回帰によりARDを行う。

ベイズ推定を用いた回帰によるARD

定式化

目的変数を $X=(x_0,\cdots,x_{N-1})$ 、説明変数を $Y=(y_0,\cdots,y_{N-1})$ とする。ここで、 $x_n\in\mathbb{R}^M,y_n\in\mathbb{R}$ である。潜在変数 $\theta$ を導入し、同時確率分布 $p(X,Y,\theta)$ を考える。この分布にベイズの定理を適用すると次式を得る。

(1) $\begin{equation*} p(\theta|X,Y)=\frac{p(Y|X,\theta)p(\theta)}{p(Y|X)} \end{equation*}$

尤度 $p(Y|X,\theta)$ として

(2) $\begin{eqnarray*} p(Y|X,\theta)&=&\prod_{n=0}^{N-1}p(y_n|x_n,\theta) \\ p(y_n|x_n,\theta)&=&\mathcal{N}(y_n|w^T x_n,\alpha^{-1}) \end{eqnarray*}$

を仮定する。 $\mathcal{N}(y|\mu,\sigma^2)$ は平均 $\mu$ 、標準偏差 $\sigma$ を持つ正規分布を表す。 $\theta=(w,\alpha)$ とした。 $w\in \mathbb{R}^M,\alpha\in\mathbb{R}$ である。これらのパラメータに対する事前確率 $p(\theta)$ として次式を仮定する。

(3) $\begin{eqnarray*} p(\theta) &=& p(w,\alpha)\\ &=& p(w|\lambda)p(\lambda)p(\alpha) \end{eqnarray*}$

ここで、新たなパラメータ $\lambda=(\lambda_0,\cdots,\lambda_{M-1})$ を導入した。各確率分布を次式で定義する。

(4) $\begin{eqnarray*} p(w|\lambda)&=&\mathcal{N}(w|0_M,A^{-1})\\ &=&\prod_{m=0}^{M-1}\mathcal{N}(w_m|0,\lambda_m^{-1}) \\ A&=&diag(\lambda_0,\cdots,\lambda_{M-1}) \\ p(\alpha)&=&{\rm Gam}(\alpha|a,b) \\ p(\lambda_m)&=&{\rm Gam}(\lambda_m|c,d) \\ \end{eqnarray*}$

${\rm Gam}(\lambda|c,d)$ はハイパーパラメータ $c,d$ を持つガンマ分布、 $0_M$ は $M$ 次元のゼロベクトル、 $A$ は $M\times M$ の対角行列である。 $w,\alpha,\lambda$ は推定すべきパラメータであり、 $a,b,c,d$ はエビデンス $p(Y|X)$ が最大になるように決めるハイパーパラメータである。

さて、式(2)は回帰式を

(5) $\begin{eqnarray*} y_n&=&w^T x_n \\ &=&w_0 x_{n0} + \cdots + w_{M-1} x_{nM-1} \end{eqnarray*}$

と仮定していることを意味する（ノイズ部分は無視した）。つまり、 $w_m$ は説明変数 $x_{nm}$ の寄与の大きさを表す。これがベイズ推定によるARDの仕組みである。

仮定した全ての確率分布を式(1)に代入すると、次の事後確率分布を厳密に（解析的に）求めることができる。

$p(w|X,Y)$

$p(\alpha|X,Y)$

$p(\lambda|X,Y)$

事後確率分布 $p(w|X,Y)$ から、個々の説明変数の寄与の大きさを知ることができる。また、次式により未知データ $x_{*}$ に対する $y_{*}$ を予測することができる。

(6) $\begin{eqnarray*} p(y_*|x_*,X,Y)=\int dw\;d\alpha\;d\lambda\;p(y_*|x_*,w,\alpha,\lambda)p(w|X,Y)p(\alpha|X,Y)p(\lambda|X,Y) \end{eqnarray*}$

データ作成

回帰に用いるデータを次の手順で作成した。

$N=130,M=3$

とした。つまり、3次元ベクトル $x_n=(x_{n0},x_{n1},x_{n2})$

を130個用意する。

[0,1)の範囲で乱数を生成し、これを $x_{n0}$

とした。

[0,100)の範囲で乱数を生成し、これを $x_{n1}$

とした。

[0,0.1)の範囲で乱数を生成し、これを $x_{n2}$

とした。

次式により $y_n$

を生成した。

(7) $\begin{equation*} y_n=\cos\left({\frac{x_{n0}\;x_{n2}}{2}\right)+\frac{x_{n0}\;x_{n2}}{10} \end{equation*}$

説明変数 $x_{n1}$ は目的変数 $y_n$ に全く寄与せず、各説明変数のスケールは大きく異なる状況をシミュレートしたものである。訓練データを100個、テストデータを30個とした。前者で学習を行い、後者で検証を行う。

学習

学習の前に以下の前処理を行う。

(8) $\begin{eqnarray*} x_{nm}^{\prime}&=&\frac{x_{nm}-\mu_{m}}{\sigma_{m}}\\ \mu_{m}&=&\frac{1}{N}\sum_{n=0}^{N-1} x_{nm}\\ \sigma_{m}^2&=&\frac{1}{N}\sum_{n=0}^{N-1} (x_{nm}-\mu_m)^2 \end{eqnarray*}$

各軸（各次元）ごとに、平均値を0に標準偏差を1に正規化する。このあと、ARDを行う。Pythonライブラリsklearnに、上の「定式化」で示したロジックを実装したクラスARDRegressionがある。これを用いると、学習部分のコードは以下のようになる（ソースはここのsample_sklearn.pyである）。

def train(x, y):
    print("training ...")
    clf = ARDRegression(n_iter=1000, compute_score=True)
    clf.fit(x, y)
    print("training done!")
    return clf

def train(x, y):

print("training ...")

clf = ARDRegression(n_iter=1000, compute_score=True)

clf.fit(x, y)

print("training done!")

return clf

3行目でARDRegressionオブジェクトを作成し、4行目で学習する。次に結果を示す。

結果

最初に示すのは、 $(w_0,w_1,w_2)$ の平均値と標準偏差である。これは $p(w|X,Y)$ から求まる量である。

期待通りの結果である。 $x_0,x_2$ からの寄与が大きく、 $x_1$ からの寄与はほとんどない。標準偏差は小さ過ぎるので上図では判別できなくなっている。次に示すのは30個のテストデータを予測した結果である。

$|{\rm Truth - Prediction}|$ の平均値は0.266となった。

ガウス過程を用いた回帰によるARD

定式化

出発点は先と同じである。

(9) $\begin{equation*} p(\theta|X,Y)=\frac{p(Y|X,\theta)\;p(\theta)}{p(Y|X)} \end{equation*}$

パラメータ $\theta$ をガウス過程に使うカーネルに関わる部分 $\theta_K$ とそうでない部分 $\theta_L$ の2つに分離する。さらに、 $X$ と $Y$ の間に次式を仮定する。

(10) $\begin{eqnarray*} y_n&=&f(x_n)+\epsilon_n \\ \epsilon_n&\sim&\mathcal{N}(\epsilon_n|0,\sigma^2) \end{eqnarray*}$

ここで、 $\mathcal{N}(\epsilon_n|0,\sigma^2)$ は平均 $0$ 、分散 $\sigma^2$ の正規分布を表す。各観測値は独立同分布から生成されると仮定すると、尤度は

(11) $\begin{eqnarray*} p(Y|X,\theta_L) &=&\prod_{n=0}^{N-1} p(y_n|x_n,\theta_L) \\ &=&\prod_{n=0}^{N-1} \mathcal{N}(y_n|f(x_n),\sigma^2) \end{eqnarray*}$

と書くことができる。ここで、 $\theta_L=(f(\cdot),\sigma)$ とした。 $\sigma$ の事前分布として $p(\sigma)$ を、関数 $f(\cdot)$ の事前分布として次のガウス過程を仮定する。

(12) $\begin{equation*} p(f|\theta_K)=\mathcal{N}(f|0_N,K_{NN}^{\theta_K}) \end{equation*}$

ここで、 $f=\left(f(x_0),\cdots,f(x_{N-1})\right)$ とした。 $0_N$ は $N$ 次元のゼロベクトル、 $K_{NN}^{\theta_K}$ はパラメータ $\theta_K$ を持つ $N\times N$ の共分散行列、その成分 $(K_{NN}^{\theta_K})_{n,n^{\prime}}=k_{NN}^{\theta_K}(x_n,x_{n^{\prime}})$ がカーネル関数となる。カーネル関数として次のものを用いる。

(13) $\begin{equation*} k_{NN}^{\theta_K}(x_n,x_{n^{\prime}})=\exp{\left(-\frac{1}{2}\sum_{m=0}^{M-1}\left(\frac{x_{nm}-x_{n^{\prime}m}}{l_m}\right)^2\right)} \end{equation*}$

ここで、 $\theta_K=(l_0,\cdots,l_{M-1})$ とした。 $l_m$ は各軸のスケール長を表す。この逆数が説明変数の重みに相当する。さらに、 $\theta_K$ の事前分布 $p(\theta_K)$ も導入する。ここまでの結果をまとめると、式(9)は次のようになる。

(14) $\begin{equation*} p(f,\sigma,\theta_K|X,Y)=\frac{p(Y|X,\sigma,f)\;p(\sigma)\;p(f|\theta_K)\;p(\theta_K)}{p(Y|X)} \end{equation*}$

$f$ の事前分布 $p(f|\theta_K)$ がガウス過程であることに注意する。上式をベイズ推定により評価することもできるが、ここでは右辺分子を取り出しMAP推定を行う。右辺分子の対数を取ると

(15) $\begin{equation*} L=\log{\left(p(Y|X,\sigma,f)\;p(\sigma)\;p(f|\theta_K)\;p(\theta_K)\right)} \end{equation*}$

を得る。これを最大にするような $\theta_K$ と $\sigma$ を求める。

(16) $\begin{equation*} (\theta_K^{MAP},\sigma^{MAP})={\rm arg}\max_{\theta_K,\sigma}L \end{equation*}$

この2つのパラメータについては点推定であることに注意する（確率を求めるのではない）。未知の値 $x_{\star}$ に対する $f_\star=f(x_\star)$ の確率分布は次式から計算される。

(17) $\begin{equation*} p(f_\star|\theta_K^{MAP})=\int df\;p(f_\star|f,\theta_K^{MAP})\;p(f|\theta_K^{MAP}) \end{equation*}$

右辺の積分内の最初の項 $p(f_\star|f,\theta_K^{MAP})$ は $p(f(x_\star)|f(x_0),\cdots,f(x_{N-1}),\theta_K^{MAP})$ を表し、式(12)から解析的に求めることができる条件付き確率分布である。上式から $f_{\star}$ の分布が分かれば、未知の値 $x_{\star}$ に対する $y_\star$ の値は次式から決定することができる。

(18) $\begin{eqnarray*} f_\star&\sim&p(f_\star|\theta_K^{MAP}) \\ \epsilon_\star&\sim&\mathcal{N}(\epsilon_\star|0,\left(\sigma^{MAP}\right)^2) \\ y_\star&=&f_\star+\epsilon_\star \end{eqnarray*}$

以上で定式化できたので、実際に計算を行う。

学習

上の定式化で示したロジックは、ガウス過程のためのPythonライブラリGPyに実装されている。先と同じデータを使用し、先と同じ前処理を施したあと、以下を実行する（ソースはここのsample_gp_multi_dim.ipynbである）。

# カーネルを定義する。
kernel = GPy.kern.RBF(input_dim=n_features, variance=0.01, ARD=True)
# kernel = GPy.kern.Matern52(n_features, ARD=True)

# ガウス過程を用いた回帰を行う。
model = GPy.models.GPRegression(train_xs, train_ys, kernel)

# 最適化（MAP推定）を行う。
model.optimize(messages=True, max_iters=1e5)

# カーネルを定義する。

kernel = GPy.kern.RBF(input_dim=n_features, variance=0.01, ARD=True)

# kernel = GPy.kern.Matern52(n_features, ARD=True)

# ガウス過程を用いた回帰を行う。

model = GPy.models.GPRegression(train_xs, train_ys, kernel)

# 最適化（MAP推定）を行う。

model.optimize(messages=True, max_iters=1e5)

2行目でカーネル関数（式(13)に相当する）を定義し、6行目と9行目でガウス過程による回帰を行う（ベイズ推定ではなくMAP推定が行われる）。

結果

最初に重み $(1/l_0,1/l_1,1/l_2)$ の値を示す。

期待通りの結果である。 $x_0,x_2$ からの寄与が大きく、 $x_1$ からの寄与はほとんどない。 $l_m$ は点推定される量であるから、「ベイズ推定を用いた回帰によるARD」の場合とは異なり標準偏差は存在しないことに注意する。次に、テストデータを予測した結果を示す。

標準偏差は小さ過ぎて上図では判別できなくなっている。 $|{\rm Truth - Prediction}|$ の平均値は0.212となった。「ベイズ推定を用いた回帰によるARD」のときは0.266であったから改善されている。ガウス過程の場合、非線形回帰となるため精度が上がるのだろう。

まとめ

今回は、2つの文脈でAutomatic Relevance Determination（ARD）を見た。

ベイズ推定を用いた回帰によるARD
ガウス過程を用いた回帰によるARD

1の回帰は線形回帰であるが、2の回帰は非線形回帰である。今回のデータに対しては、予測精度は2の方が高い。また、1の計算では $M\times M$ の行列を、2の計算では $N\times N$ の行列を扱う。したがって、今回のような $M<N$ の場合は、1の方が高速に計算できる。

Kumada Seiya

仕事であろうとなかろうと勉強し続ける、その結果”中身”を知ったエンジニアになれる

One Class Support Vector Machine(One Class SVM)入門

点群からの法線の導出

Automatic Relevance Determination（ARD）

はじめに

ベイズ推定を用いた回帰によるARD

定式化

データ作成

学習

結果

ガウス過程を用いた回帰によるARD

定式化

学習

結果

まとめ

Kumada Seiya

最近の記事

レイアウト解析

ChatGPTを用いた画像異常検知

クォータニオンと回転行列の関係

OpenCVで深層学習

OpenCVのcv::Matに関する小ネタ集

テンソルネットワークの入り口（その2）

線形回帰をPythonで数式から逃げずに実装してみた

PyMCによるMarkov Chain Monte Carlo

CatBoostのランク学習（Learning to rank）をためそう

Learning Chainer with Examples 〜深層学習への入門〜

LightGBMによるFIFA18の選手データでの市場価値推定

Doc2Vecでpythonスクリプトを学習

Google Vision APIでOCR

アーカイブ

カテゴリー