pythonを使って、libsvmにおける多項式カーネル使用時の写像後の素性の重みベクトルの計算してみたいと思います。
多項式カーネル
多項式カーネルは複数の素性の組み合わせを見ることができると言われています。
libsvmにおける多項式カーネルは以下のようになっています。$g, r, d$は、それぞれ学習時のパラメータです。
K(\mathbf{x}, \mathbf{y}) = (g\space\mathbf{x\cdot y} + r) ^ d
例えば、元々の特徴ベクトルを素性が2つ、すなわち、
\mathbf{x} = [x_{1}, x_{2}]
とすると、$d=2$の時のカーネル関数は、
\begin{aligned}
K(\mathbf{x}, \mathbf{y}) &= (g\space\mathbf{x\cdot y} + r) ^ 2 \\
&= (gx_{1}y_{1}+gx_{2}y_{2} + r)^2 \\
&= (gx_{1}y_{1})^2 + (gx_{2}y_{2})^2 + 2g^2x_{1}x_{2}y_{1}y_{2} + 2grx_{1}y_{1} + 2grx_{2}y_{2} + r^2\\
&=(gx_1^2)(gy_1^2) + (gx_2^2)(gy_2^2) + (\sqrt{2}gx_1x_2)(\sqrt{2}gy_1y_2) \\
&+ (\sqrt{2gr}x_1)(\sqrt{2gr}y_1) +(\sqrt{2gr}x_2)(\sqrt{2gr}y_2) + (r)(r) \\
&= \phi(\mathbf{x}) \cdot \phi(\mathbf{y})
\end{aligned}
すなわち、特徴ベクトル$\phi(\mathbf{x})$は、
\phi(\mathbf{x}) = [gx_1^2, gx_2^2, \sqrt{2}gx_{1}x_{2}, \sqrt{2gr}x_1, \sqrt{2gr}x_2, r]
というように特徴空間を写像します。特に3番目の素性を見てみると、$x_1x_2$の組み合わせ素性を作っていることが分かります。
これを一般化すると多項定理により、
\phi(\mathbf{x}) = \Biggl[ \frac{d!}{\prod_j p_j!}\prod_j b_j^{p_j}, ... \Biggr] (\forall [p_1, p_2, ..., p_J] | \sum_j p_j = d, 0 \leq p_j \leq d)
\phi(\mathbf{x}) = [\phi_1(\mathbf{x}), \phi_2(\mathbf{x}), ..., \phi_I(\mathbf{x})]
\phi_i(\mathbf{x}) = \frac{d!}{\prod_j p_{i,j}!}\prod_j b_j^{p_{i,j}}
\mathbf{p} = \bigl\{ \mathbf{p_i} | \sum_j p_{i,j} = d, 0 \leq p_j \leq d \bigr\}
ただし
b = [\sqrt{g}x_1, \space \sqrt{g}x_2, \space ..., \space \sqrt{g}x_n, \space \sqrt{r}]
$J$は入力ベクトルの素性数+1となっています。
また、$J=3, d=2$の時、
\mathbf{p} = \bigl[ [0, 0, 2], [0, 1, 1], [0, 2, 0], [1, 0, 1], [1, 1, 0], [2, 0, 0] \bigr]
となります。
素性の重み
svmでは識別関数は、分類したいデータを$\mathbf{x}_{d}$とすると、
\begin{aligned}
f(\mathbf{x}_{d}) &= \sum_{i}\alpha_{i}y_{i}K(\mathbf{x}_i, \mathbf{x}_d) - b \\
&=\sum_{i}\alpha_{i}y_{i}\phi(\mathbf{x}_i)\cdot\phi(\mathbf{x}_d) - b
\end{aligned}
と表現できます。ただし、$\mathbf{x}_{i}, \alpha_i, y_i$はそれぞれ、サポートベクトル、ラグランジュ乗数、サポートベクトルのラベルを表しています。
また、
f(\mathbf{x}_d) = \mathbf{w}\cdot\phi(\mathbf{x}_{d})-b
ただし
\mathbf{w} = \sum_{i}\alpha_{i}y_{i}\phi(\mathbf{x}_i)
と表せ、この$\mathbf{w}$が、写像後の特徴ベクトル$\phi(\mathbf{x})$に対する重みベクトルとなります。
wは何に対する重みベクトルなのか
線形カーネル
線形カーネルの場合は、
\phi(\mathbf{x}) = \mathbf{x}
入力ベクトルをそのまま用いるため、入力ベクトルに対応する重みベクトルが計算できます。
例えば、入力ベクトルが
\mathbf{x} = [x_1, x_2]
の場合、$\mathbf{x}$に対する重みベクトルである
\mathbf{w} = [w_1, w_2]
が計算できます。
多項式カーネル
一方、多項式カーネルの場合は、入力ベクトルは異なる特徴空間に写像されるため、写像後の特徴ベクトルに対する重みベクトルが計算できます。
例えば、入力ベクトルが
\mathbf{x} = [x_1, x_2]
の場合、特徴ベクトルは、
\phi(\mathbf{x}) = [gx_1^2, gx_2^2, \sqrt{2}gx_{1}x_{2}, \sqrt{2gr}x_1, \sqrt{2gr}x_2, r]
と写像され、$\phi(\mathbf{x})$に対する重みベクトルである、
\mathbf{w} = [w_1, w_2, w_3, w_4, w_5, w_6]
が計算できます。
libsvmにおける重み計算
libsvmにおいて、訓練後に生成できるモデルファイル内には、サポートベクトル$\mathbf x_i$と、その係数$\alpha_iy_i$ が記述されています。
重みベクトルを求めるには、これらを用いて、
\mathbf{w} = \sum_{i}\alpha_{i}y_{i}\phi(\mathbf{x}_i)
を計算することになります。
線形カーネル
線形カーネルの場合は以下のように計算できます。
from collections import defaultdict
from itertools import izip
from svm import svm_parameter
from svmutil import svm_problem, svm_train
def linear_kernel_mapper(support_vector):
return support_vector
def calc_weight(model, mapping_function):
weight_dict = defaultdict(float)
for support_vector, coef in izip(model.get_SV(), model.get_sv_coef()):
support_vector.pop(-1)
mapped_vector = mapping_function(support_vector)
for feature, value in mapped_vector.items():
weight_dict[feature] += coef[0] * value
return weight_dict
if __name__ == '__main__':
data = [[0, 1], [1, 1], [1, 0]]
label = [-1, 1, -1]
study_data = svm_problem(label, data)
params = svm_parameter('-t 0 -c 1')
model = svm_train(study_data, params)
print calc_weight(model, linear_kernel_mapper)
model.get_SV()はサポートベクトル$\mathbf x_i$のリストを、
model.get_sv_coef()はget_SV()に対応した$\alpha_iy_i$(が[0]に入ったタプル)のリストを返します。
サポートベクトル$\mathbf x_i$は、{素性番号:値}のdictになっています。通常libsvmでは、素性番号は1から始まるため参照する際は注意が必要かもしれません。
linear_kernel_mapper()は線形カーネルの写像関数$\phi(\mathbf x) = \mathbf x$を表しています。
多項式カーネル
多項式カーネルの場合は以下のように計算できます。
from collections import defaultdict, Counter
from itertools import izip, combinations_with_replacement
from math import sqrt, pow, factorial
from operator import mul
from svm import svm_parameter
from svmutil import svm_problem, svm_train
class PolynomialKernelMapper(object):
def __init__(self, g, r, d):
self._g = g
self._r = r
self._d = d
def __call__(self, support_vector):
mapped_vector = defaultdict(float)
features = support_vector.copy()
features['r'] = sqrt(self._r)
for feature_index in features.keys():
features[feature_index] *= sqrt(self._g)
#[A, B, C] = [(AAA), (AAB), (AAC), (ABA), ..., (CCC)]
#のように展開
for combo in combinations_with_replacement(features, self._d):
counter = Counter(combo)
combo_coef = sqrt(factorial(self._d) / reduce(mul, [factorial(count) for count in counter.values()])) #多項定理により多項係数を計算
f_value_multiplied = reduce(mul, [features[f_index] for f_index in combo])
mapped_vector[combo] += combo_coef * f_value_multiplied
return mapped_vector
#再掲
def calc_weight(model, mapping_function):
weight_dict = defaultdict(float)
for support_vector, coef in izip(model.get_SV(), model.get_sv_coef()):
support_vector.pop(-1)
mapped_vector = mapping_function(support_vector)
for feature, value in mapped_vector.items():
weight_dict[feature] += coef[0] * value
return weight_dict
if __name__ == '__main__':
data = [[0, 1], [1, 1], [1, 0]]
label = [-1, 1, -1]
study_data = svm_problem(label, data)
params = svm_parameter('-t 1 -g 1 -r 1 -d 2')
model = svm_train(study_data, params)
print calc_weight(model, PolynomialKernelMapper(g=1, r=1, d=2))
PolynomialKernelMapper()は、多項式カーネルの写像関数$\phi(\mathbf x)$を表しています。
結果は、以下のようになります。
defaultdict(<type 'float'>, {(1, 2): 0.8079544634564579, (2, 'r'): 0.4037817308453665, ('r', 'r'): 0.0, (1, 'r'): 0.4041727326110915, (2, 2): 0.28551679999999996, (1, 1): 0.28579328})
(1, 1)は$x_{1}^2$の項の素性、(1, 2)は$x_{1}x_{2}$の項の素性、(1, r)は$x_{1}$の項の素性といった具合に表しています。
この例では、(1, 2)、すなわち$\sqrt{2g}x_{1}x_{2}$を表す素性の重みが大きくなっています。
おまけ：他のカーネルの場合
他のカーネルの場合、どうなるか考えてみます。
RBF
RBFのカーネル関数は、
K(\mathbf{x}, \mathbf{y}) = \exp(-g||\mathbf{x} - \mathbf{y}||^2)
と定義されます。$g(g \gt 0)$はパラメータです。これを展開すると、
\begin{aligned}
K(\mathbf{x}, \mathbf{y}) &= \exp(-g\mathbf{x}\cdot\mathbf{x} + 2g \mathbf{x}\cdot\mathbf{y} - g\mathbf{y}\cdot\mathbf{y}) \\
&= \exp(-g\mathbf{x}\cdot\mathbf{x})\exp(-g\mathbf{y}\cdot\mathbf{y})\exp(2g\mathbf{x}\cdot\mathbf{y})
\end{aligned}
べき級数展開により、
\begin{aligned}
\exp(2g\mathbf{x}\cdot\mathbf{y}) &= \sum_{d=0} ^ {\infty} \frac{(2g)^d}{d!}( \mathbf{x}\cdot\mathbf{y})^d
\end{aligned}
また多項定理により、
\begin{aligned}
(\mathbf{x}\cdot\mathbf{y})^n &= \Biggl( \frac{n!}{\prod_j p_j!}\prod_j b_j^{p_j} + ... \Biggr) (\forall [p_1, p_2, ..., p_J] | \sum_j p_j = n, 0 \leq p_j \leq n)
\end{aligned}
ただし、
b = [x_1y_1, \space x_2y_2, \space ..., \space x_Jy_J]
$J$は入力ベクトルの素性数です。
2次の入力ベクトル$\mathbf {x} = [x_1, x_2]$に対し、2項までべき級数展開を行ってみると、
\begin{aligned}
&(\mathbf{x}\cdot\mathbf{y})^0 = 1 \\
&(\mathbf{x}\cdot\mathbf{y})^1 = x_1y_1 + x_2y_2 = (x_1, x_2)\cdot(y_1, y_2) \\
&(\mathbf{x}\cdot\mathbf{y})^2 = x_1^2y_1^2 + 2x_1x_2y_1y_2 + x_2^2y_1^2 = (x_1^2, \sqrt{2}x_{1}x_{2}, x_2^2)\cdot(y_2^2, \sqrt{2}y_1y_2, y_2^2)
\end{aligned}
これよりカーネル関数$K(\mathbf x, \mathbf y)$と写像関数$\phi(\mathbf x)$は、
\begin{aligned}
K(\mathbf{x}, \mathbf{y}) &= \exp(-g\mathbf{x}\cdot\mathbf{x}) \exp(-g\mathbf{y}\cdot\mathbf{y}) \sum_{d=0} ^ {2} \frac{(2g)^d}{d!}( \mathbf{x}\cdot\mathbf{y})^d \\
&= \exp(-g\mathbf{x}\cdot\mathbf{x}) \exp(-g\mathbf{y}\cdot\mathbf{y}) \Biggl(1 + 2g(x_1y_1 + x_2y_2) + 2g^2 (x_1^2y_1^2 + 2x_1x_2y_1y_2 + x_2^2y_1^2) \Biggr) \\
&= \exp(-g\mathbf{x}\cdot\mathbf{x}) \exp(-g\mathbf{y}\cdot\mathbf{y}) \\
& \Biggl((1)\cdot(1) + (\sqrt{2g}x_1)(\sqrt{2g}y_1) + (\sqrt{2g}x_2)(\sqrt{2g}y_2) + (\sqrt{2}gx_1^2)(\sqrt{2}gy_1^2) + (2gx_1x_2)(2gy_1y_2) + (\sqrt{2}gx_2^2)(\sqrt{2}gy_1^2) \Biggr) \\
\phi(\mathbf{x}) &= \exp(-g\mathbf{x}\cdot\mathbf{x})\Biggl[1, \sqrt{2g}x_1, \sqrt{2g}x_2, \sqrt{2}gx_1^2, 2gx_1x_2, \sqrt{2}gx_2^2 \Biggr]
\end{aligned}
このように2項までの展開を行った場合は、2次元から6次元に写像されます。多項式カーネルと同様複数の素性の組み合わせを見ていることが分かります。
実際のRBFカーネルでは無限次元に写像され（すなわち全ての素性の組み合わせを見る）、
\phi(x) = \exp(-g\mathbf{x}\cdot\mathbf{x}) \Biggl[ \sqrt{\frac{(2g)^d}{d!}} \sqrt{\frac{d!}{\prod_j p_j!}} {\displaystyle \prod_j b_j^{p_j}}, ... \Biggr] (d=0, 1, 2, ... \infty, \forall [p_1, p_2, ..., p_J] | \sum_j p_j = d, 0 \leq p_j \leq d)
上のコード例と同様にこの写像関数を実装すれば、RBFの重みベクトルも有限項まで出せそうです。