関連ベクトルマシンをPythonコードを交えて紹介

はじめに

PRML(Pattern Recognition and Machine Learning)で関連ベクトルマシンについて学んだ内容をまとめて、実際のデータを使って学習しました。主に7章2節の内容です。学習アルゴリズムの全体像を示すことを意識したので、導出は端折っている部分が多いです。詳しい解の導出は、PRMLの本文や演習問題にあります。そのほかに参考になるものは最後にまとめてあります。

関連ベクトルマシンのアルゴリズム

変数一覧

• 訓練データの数 $N$
• 基底関数の数 $M$
• 訓練データの説明変数 $\mathbf{X} = \left\{\mathbf{x}_1, \mathbf{x}_2, \cdots \mathbf{x}_N\right\}, \mathbf{x}_n = (x_{n, 1}, x_{n, 2}, \cdots, x_{n_D})^T$
• 訓練データの目的変数 $\mathbf{t} = \left\{t_1, t_2, \cdots t_N\right\}$
• 基底関数 $\mathbf{\phi}(\mathbf{x}_n) = (\phi_1(\mathbf{x}_n), \cdots \phi_{M}(\mathbf{x}_n))^T$
• 重み $\mathbf{w} = (w_1, \cdots, w_{M})^T$
• 新たな説明変数 $\mathbf{x}$

また以下の計画行列$\mathbf{\Phi}$を使います。

$$\mathbf{\Phi} = \begin{pmatrix} \phi_0(\mathbf{x}_1) & \phi_1(\mathbf{x}_1) & \cdots & \phi_{M-1}(\mathbf{x}_1) \\ \phi_0(\mathbf{x}_2) & \phi_1(\mathbf{x}_2) & \cdots & \phi_{M-1}(\mathbf{x}_2) \\ \vdots & \vdots & \ddots & \vdots \\ \phi_0(\mathbf{x}_N) & \phi_1(\mathbf{x}_N) & \cdots & \phi_{M-1}(\mathbf{x}_N) \\ \end{pmatrix}$$

目的とモデル

最終的な目標は以下の予測分布を求めることです。訓練データから、新たな説明変数に対する目的変数の分布を求めます。

$$\begin{aligned} p(t | \mathbf{x}, \mathbf{t}, \mathbf{X}) \end{aligned}$$

今回は以下のような線形回帰モデルを使用します。

$$\begin{aligned} p(t | \mathbf{x}, \mathbf{w}, \beta) = N(t | \mathbf{w}^T\mathbf{\phi(\mathbf{x}}), \beta^{-1}) \end{aligned}$$

これ以降、見やすさのために$\mathbf{x}$と$\mathbf{X}$は省略します。

事前分布と事後分布

事前分布を一つのハイパーパラメータ$\alpha$を用いて、

$$\begin{aligned} p(\mathbf{w}|\alpha) = \mathcal{N}(\mathbf{w}|\mathbf{0}, \alpha^{-1}\mathbf{I}) \end{aligned}$$

とした場合、以下の記事のようになります。
ベイズ線形回帰とエビデンス近似
今回は、それぞれのパラメータ$w_i$に対して異なる$\alpha_i$を考えています。導出は似ている部分が多いです。

パラメータ$\mathbf{w}$についての事前分布と尤度関数は以下の通りです。

$$\begin{aligned} p(\mathbf{t} | \mathbf{w}, \beta) &= N(\mathbf{t} | \mathbf{\Phi}\mathbf{w}, \beta^{-1}\mathbf{I}) \\ p(\mathbf{w} | \mathbf{\alpha}) &= \prod^M_{i=1}\mathcal{N}(w_i|0, \alpha_i^{-1}) \end{aligned}$$

ここで、$\mathbf{\alpha} = (\alpha_1, ... , \alpha_M)^T$です。よって、事後分布は次のようになります。

$$\begin{aligned} p(\mathbf{w} | \mathbf{t}, \mathbf{\alpha}, \beta) &= N(\mathbf{w} | \mathbf{m}, \mathbf{\Sigma}) \\ where \quad \mathbf{m} &= \beta\mathbf{\Sigma}\mathbf{\Phi^T}\mathbf{t} \\ \mathbf{\Sigma} &= (\mathbf{A} + \beta\mathbf{\Phi^T}\mathbf{\Phi})^{-1} \end{aligned}$$

ここで、$\mathbf{A} = diag(\alpha_i)$です。

予測分布とエビデンス近似

予測分布は以下のようにかけます。

$$\begin{aligned} p(t | \mathbf{t}) = \int\int\int p(t | \mathbf{w}, \beta) p(\mathbf{w} | \mathbf{t}, \mathbf{\alpha}, \beta) p(\mathbf{\alpha}, \beta | \mathbf{t}) d\mathbf{w}d\mathbf{\alpha} d\beta \end{aligned}$$

しかし、この積分をそのまま計算するのはこんなんです。そこで、まず$\mathbf{w}$について周辺化して得られた周辺尤度における$\mathbf{\alpha}, \beta$の最尤推定量$\mathbf{\alpha}^*, \beta^*$を求めて、その値で固定して計算することで、近似的に予測分布を求めます。この手法をエビデンス近似と言います。

$$\begin{aligned} p(t|\mathbf{t}) &\simeq p(t|\mathbf{t}, \mathbf{\alpha}^*, \beta^*) \\ &= \int p(t|\mathbf{w}, \beta^*)p(\mathbf{w}|\mathbf{t}, \mathbf{\alpha}^*, \beta^*)d\mathbf{w} \\ &= \int \mathcal{N}(t|\mathbf{w}^T\mathbf{\phi}(\mathbf{x}), (\beta^*)^{-1})\mathcal{N}(\mathbf{w}|\mathbf{m}, \mathbf{\Sigma})d\mathbf{w} \\ &= \mathcal{N}(t|\mathbf{m}^T\mathbf{\phi}(\mathbf{x}), \sigma^2(\mathbf{x})) \\ where \quad \sigma^2(\mathbf{x}) &= (\beta^*)^{-1} + \mathbf{\phi}(\mathbf{x})^T\mathbf{\Sigma}\mathbf{\phi}(\mathbf{x}) \end{aligned}$$

ここで、$\mathbf{m}, \mathbf{\Sigma}$では、$\mathbf{\alpha}^*, \beta^*$が使われています。

$\alpha_i^* \rightarrow \infty$となるとき、それに対応して$w_i \rightarrow 0$となり、$\phi_i$はモデルから除外されます。多くの場合、大部分の特徴量ベクトルがモデルから除外され、予測に影響を与えなくなります。この性質をスパース性といいます。

$\alpha, \beta$の最尤推定

$\mathbf{\alpha}, \beta$に関する事後分布を求めます。ここでは事前分布は一定と考えます。

$$\begin{aligned} p(\mathbf{\alpha}, \beta | \mathbf{t}) &\propto p(\mathbf{t} | \mathbf{\alpha}, \beta) p(\mathbf{\alpha}, \beta) \\ &\propto p(\mathbf{t} | \mathbf{\alpha}, \beta) \\ &= \int p(\mathbf{t} | \mathbf{w}, \beta)p(\mathbf{w} | \mathbf{\alpha}) d\mathbf{w} \\ &= \int N(\mathbf{t} | \mathbf{\Phi}\mathbf{w}, \beta^{-1}\mathbf{I}) N(\mathbf{w} | \mathbf{0}, \mathbf{A}^{-1}) d\mathbf{w} \\ &= \mathcal{N}(\mathbf{t}|\mathbf{0}, \mathbf{C}) \\ where \quad \mathbf{C} &= \beta^{-1}\mathbf{I} + \mathbf{\Phi}\mathbf{A}^{-1}\mathbf{\Phi}^T \end{aligned}$$

この分布の対数$lnp(\mathbf{t} | \mathbf{\alpha}, \beta)$を考えて、$\mathbf{\alpha}, \beta$について最大化するのが目標です。それぞれについて微分して0とおくと、次の更新式が得られます。

$$\begin{aligned} \alpha_i^{new} &= \frac{\gamma_i}{m_i^2} \\ (\beta^{new})^{-1} &= \frac{||\mathbf{t} - \mathbf{\Phi}\mathbf{m}||^2}{N - \sum_i \gamma_i} \\ where \quad \gamma_i &= 1 - \alpha_i\Sigma_{ii} \end{aligned}$$

この式に従って、$\mathbf{\alpha}, \beta$の値が収束するまで更新を繰り返せばいいです。

もう一つの定式化

$\mathbf{\alpha}$について別の更新式を考えます。$\mathbf{\varphi}_i$を$\mathbf{\Phi}$のi列目とします。$\mathbf{C}$のうち、$\mathbf{\varphi}_i$に依存しない項をまとめて$\mathbf{C}_{-i}$とおきます。また、$lnp(\mathbf{t} | \mathbf{\alpha}, \beta)$のうち、$\mathbf{\varphi}_i$に依存しない項をまとめて$L(\mathbf{\alpha}_{-i})$とおきます。このとき対数周辺尤度は以下のようになります。

$$\begin{aligned} lnp(\mathbf{t} | \mathbf{\alpha}, \beta) &= L(\mathbf{\alpha}_{-i}) + \lambda(\alpha_i) \\ where \quad \lambda(\alpha_i) &= \frac{1}{2}\left[ln\alpha_i - ln(\alpha_i + s_i) + \frac{q_i^2}{\alpha_i + s_i}\right] \\ s_i &= \mathbf{\varphi}_i^T\mathbf{C}_{-i}^{-1}\mathbf{\varphi}_i \\ q_i &= \mathbf{\varphi}_i^T\mathbf{C}_{-i}^{-1}\mathbf{t} \end{aligned}$$

ここで、スパース性係数$s_i$は、基底ベクトル$\mathbf{\varphi}_i$がモデルにすでに存在する基底ベクトルと重なっている程度を表す指標と見なすことができます。クオリティ係数$q_i$は$\mathbf{\varphi}_i$を除いたモデルの誤差と$\mathbf{\varphi}_i$の整合性の尺度です。$q_i$に対して$s_i$が大きいと、$\mathbf{\varphi}_i$がモデルから除外される可能性が高くなります。$q_i^2 \lt s_i$のときは、$\alpha_i \rightarrow \infty$とすれば周辺尤度が最大となります。$q_i^2 \gt s_i$のときの増減表は以下のようになりますが、このときは$\alpha_i = \frac{s_i^2}{q_i^2 - s_i}$とすれば周辺尤度が最大となります。

$\alpha_i$	$0$	$\cdots$	$\frac{s_i^2}{q_i^2-s_i}$	$\cdots$	$\frac{3s_i^2 + s_i\sqrt{s_i(s_i + 8q_i^2)}}{4(q_i^2 - s_i^2)}$	$\cdots$
$\frac{d\lambda(\alpha_i)}{d\alpha_i}$	$\times$	$+$	$0$	$-$	$-$	$-$
$\frac{d^2\lambda(\alpha_i)}{d\alpha_i^2}$	$\times$	$-$	$-$	$-$	$0$	$+$

実装

環境

• Google Colaboratory
• Python 3.7.13
• numpy 1.21.6
• pandas 1.3.5
• sklearn 1.0.2

import文

import numpy as np
import pandas as pd
from sklearn.metrics import r2_score
from sklearn.model_selection import train_test_split

np.random.seed(28)
np.set_printoptions(precision=3)

使用するデータセット

scikit-learnの「糖尿病データセット」というトイデータセットを使います。このデータセットでは、年齢、性別、BMI、平均血圧、6つの血清検査測定値という10の説明変数から、1年後の糖尿病の進行状況を予測します。データ数は442です。このうちの8割に当たる353個を訓練用データとし、残りをテスト用データとします。また、データは標準化しておきます。

input_variables = ['age', 'sex', 'bmi', 'map', 'tc', 'ldl', 'hdl', 'tch', 'ltg', 'glu']

data_url = 'https://www4.stat.ncsu.edu/~boos/var.select/diabetes.tab.txt'
raw_df = pd.read_csv(data_url, sep="\s+")
diabetes_df = raw_df.set_axis(input_variables + ['y'], axis=1)

t_origin_df = diabetes_df.loc[:, 'y']
t_df = (t_origin_df - t_origin_df.mean()) / t_origin_df.std()
x_origin_df = diabetes_df.loc[:, input_variables]
x_df = (x_origin_df - x_origin_df.mean()) / x_origin_df.std()

train_x_df, test_x_df, train_t_df, test_t_df = train_test_split(x_df, t_df, test_size=0.2)

train_x = train_x_df.to_numpy(copy=True)
train_t = train_t_df.to_numpy(copy=True)
test_x = test_x_df.to_numpy(copy=True)
test_t = test_t_df.to_numpy(copy=True)

diabetes_df.head()

index	age	sex	bmi	map	tc	ldl	hdl	tch	ltg	glu	y
0	59	2	32.1	101.0	157	93.2	38.0	4.0	4.8598	87	151
1	48	1	21.6	87.0	183	103.2	70.0	3.0	3.8918	69	75
2	72	2	30.5	93.0	156	93.6	41.0	4.0	4.6728	85	141
3	24	1	25.3	84.0	198	131.4	40.0	5.0	4.8903	89	206
4	50	1	23.0	101.0	192	125.4	52.0	4.0	4.2905	80	135

基底関数

基底関数には各データ点を中心とするガウスカーネルを用いています。よって、基底関数の数は訓練用データの数と同じく353個です。

$$\begin{aligned} \phi_i(\mathbf{x}) &= k(\mathbf{x}, \mathbf{x}_i) \\ where \quad k(\mathbf{x}, \mathbf{x}_i) &= exp(-\gamma||\mathbf{x} - \mathbf{x}_i||^2) \end{aligned}$$

def gauss_kernel(gamma=0.1):
    gamma = gamma

    def closure(x, x_star):
        return np.exp(-gamma * np.dot(x - x_star, x - x_star))

    return closure

関連ベクトルマシン

ここでは、もう一つの定式化で示した更新アルゴリズムを用いています。

class RVMReg:
    def __init__(self, kernel):
        self.kernel = kernel
        self.xs = None
        self.ts = None
        self.alphas = None
        self.beta = None
        self.m = None

    def fit(self, xs, ts, iter_max=1000):
        INF = 1e18
        self.xs = xs
        self.ts = ts
        N = len(xs)
        Phi = np.array([np.array([self.kernel(x1, x2) for x1 in xs]) \
                        for x2 in xs])
        self.alphas = np.full(N, INF)
        self.beta = 0.1 * np.var(self.ts)
        
        for iter in range(iter_max):
            params = np.hstack([self.alphas, self.beta])
            update_flg = True

            # アルファの更新
            for i in range(N):
                phi_i = Phi[:, i]
                if update_flg:
                    inv_A = np.diag(1 / self.alphas)
                    C = np.identity(N) / self.beta + np.dot(np.dot(Phi, inv_A), Phi.T)
                    inv_C = np.linalg.inv(C)
                Qi = np.dot(np.dot(phi_i.T, inv_C), ts)
                Si = np.dot(np.dot(phi_i.T, inv_C), phi_i)
                qi = self.alphas[i] * Qi / (self.alphas[i] - Si)
                si = self.alphas[i] * Si / (self.alphas[i] - Si)
                if qi**2 > si:
                    update_flg = True
                    self.alphas[i] = si**2 / (qi**2 - si)
                else:
                    if self.alphas[i] == INF:
                        update_flg = False
                    else:
                        update_flg = True
                    self.alphas[i] = INF

            # ベータの更新
            A = np.diag(self.alphas)
            Sigma = np.linalg.inv(A + self.beta * np.dot(Phi.T, Phi))
            self.m = self.beta * np.dot(np.dot(Sigma, Phi.T), ts)
            denominator = ts - np.dot(Phi, self.m)
            denominator = np.dot(denominator.T, denominator)
            gammas = np.array([1 - self.alphas[j] * Sigma[j, j] for j in range(N)])
            numerator = N - np.sum(gammas)
            self.beta = numerator / denominator

            # 収束したら終了
            if np.allclose(params, np.hstack([self.alphas, self.beta])):
                print(f"Converged iteration: {iter+1}, The number of non-infinite alphas: {np.sum(self.alphas != INF)}")
                break

            if iter == iter_max - 1:
                print("No convergence was achieved")

    def predict(self, x_star):
        phi_x = np.array([self.kernel(x, x_star) for x in self.xs])

        return np.dot(self.m, phi_x)

結果

ガウスカーネルのパラメータは、$\gamma=0.018$としています。

model = RVMReg(gauss_kernel(0.018))
model.fit(train_x, train_t)
pred_train = np.array([model.predict(x) for x in train_x])
pred_test = np.array([model.predict(x) for x in test_x])
print(f"train_r2: {r2_score(train_t, pred_train)}, test_r2: {r2_score(test_t, pred_test)}")

Converged iteration: 23, The number of non-infinite alphas: 6
train_r2: 0.5370609799956878, test_r2: 0.5183172061767674

結果の部分にThe number of non-infinite alphas: 6とありますが、これは$\alpha_i \ne \infty$であるパラメータの数を示しています。特徴量の数は全部で353個ありますが、そのうち6個の特徴量ベクトルのみが予測に使われ、残りの347個の特徴量ベクトルはモデルから除外されていることになります。

参考

PRMLでReferenceとしてある論文です。もう一つの定式化で示したアルゴリズムが詳しくのっています。この記事では、この論文と全く同じように実装をしているわけではありません。
Tipping, M. E. and A. Faul (2003). Fast marginal likelihood maximization for sparse Bayesian models.

この記事では、もう一つの定式化で示した更新式を用いましたが、$\alpha, \beta$の最尤推定で示した更新式を用いて実装しているサイトです。
PRML第7章 回帰問題に対する関連ベクトルマシン Python実装