カーネル法を勉強し直すその3　〜カーネルの導入とカーネルトリック〜

11月の中旬〜下旬あたり、体調を崩してしまっていたので更新できませんでした。*1 寒さも厳しくなってきましたので、皆様もお気をつけください。

back:カーネル法を勉強し直すその2 〜カーネル法を導くための準備〜 - root3のメモ帳
next:まだ

カーネルの導入
- 前回の復習
- 特徴量の内積の性質
  - 1.対称性
  - 2.正値性
カーネルトリック
- カーネルトリックとは？
- カーネルトリックで嬉しいこと
実装
まとめと次回予告

カーネルの導入

前回の復習

さて、前回の記事では、ロジスティック回帰を別の視点で行うことを考えました。具体的に言うと、以下の3つです。

特徴写像を導入し、
確率の推定式を内積のみで表記し、
特徴量の係数に適切な仮定をおくことで、問題を置き換える

さらにこの事から、特徴量の内積さえ求めておけばロジスティック回帰問題を解くことが可能であることを示しました。

この事実がカーネルの導入に繋がります。

特徴量の内積の性質

前回の記事での例を思い出してください。同心円状に分布した点を分類するために、以下の特徴写像を導入しました。

$\displaystyle \Phi(x_1, x_2) = (1, x_1, x_2, {x_1}^2+{x_2}^2)$

この特徴写像によって変換された特徴量同士の内積は、次のように計算できます。（ただし、 ${\bf x} = (x_1, x_2)$ と表記しています）

$\displaystyle \langle \Phi({\bf x}), \Phi({\bf y}) \rangle \\ = \langle (1, x_1, x_2, {x_1}^2+{x_2}^2), (1, y_1, y_2, {y_1}^2+{y_2}^2) \rangle \\ = 1 + x_1 y_1 + x_2 y_2 + {x_1}^2 {y_1}^2 + {x_1}^2 {y_2}^2 + {x_2}^2 {y_1}^2 + {x_2}^2 {y_2}^2$

さて、上記の関数は、最後の式だけ見ればただのxとyの2変数関数です。このことを明記するand記述を簡単にするために、

$\displaystyle K({\bf x}, {\bf y}) = \langle \Phi({\bf x}), \Phi({\bf y}) \rangle$

と定義しておきます。

さて、 $K({\bf x}, {\bf y})$ は一見ただの2変数関数ですが、特徴量の内積から導出したために、一般的な2変数関数にはない性質を2つ持っています。
これは今回の例のみで成立する訳ではなく、 $K({\bf x}, {\bf y}) = \langle \Phi({\bf x}), \Phi({\bf y}) \rangle$ 全てに対して主張できるものとなっております。

1.対称性

$K({\bf x}, {\bf y})$ は、xとyを入れ替えても値が変わりません。*2
式で書くと以下のようになります。

$\displaystyle \forall {\bf x}, {\bf y}\,\, K({\bf x}, {\bf y}) = K({\bf y}, {\bf x})$

2.正値性

任意の自然数nと、任意のn個の点 ${\bf x_1}, \cdots, {\bf x_n}$ と、任意のn個の実数 $c_1, \cdots, c_n$ に対して、以下の式が成立します。

$\displaystyle \sum_{i, j=1} ^n c_i c_j K({\bf x_i}, {\bf x_j}) \geq 0$

上記の2つの性質をもつ2変数関数を、正定値カーネルと呼ぶことにします。つまり、特徴量の内積から構成した2変数関数 $K({\bf x}, {\bf y}) = \langle \Phi({\bf x}), \Phi({\bf y}) \rangle$ は正定値カーネルであるという事です。

カーネルトリック

ここまででようやくカーネル法の凄いところを説明する準備が整いました。そういう意味ではここまでは準備期間です。

カーネルトリックとは？

ここまでやってきた事で重要なのは以下の2つです。

特徴量の内積がわかれば、データそのものを直接使わなくてもロジスティック回帰問題を解くことができる
特徴量の内積は「正定値カーネル」と呼ばれるクラスに属している

この2つの事実から、偉い人は次のような発想にたどり着きました。

正定値カーネルさえ与えれば、特徴写像を与える事なく回帰問題が解けるのではないだろうか？

もちろんこの主張はこのままではただの仮説にすぎません。なぜなら、「特徴量の内積は正定値カーネルである」は正しくても、その逆である「正定値カーネルは、特徴量の内積になっている」は示されていないからです。

さらに言うと、「正定値カーネルは、特徴量の内積になっている」の証明は非常に難しいことが推測されます。「存在自体を示すこと」が、「存在するものの性質を示すこと」に比べてはるかに難易度が高いのはなんとなく感じて頂けるかと思います。

結論を言うと、この主張は正しいことが示されています。すなわち、正定値カーネルを1つ決めることで、そのカーネルに対応する特徴写像が存在することが証明されています。*3
この事実から、正定値カーネルを1つ決めて、それを回帰に利用することで、特徴写像を決めることなく、何かしらの特徴写像で変換を行った回帰分析を行っていることになります。この一見気持ち悪い事実を、気持ち悪さの意味も込めてカーネルトリックと呼んでいる訳ですね。

繰り返しですが、この気持ち悪い（でも素晴らしい）結果は、「正定値カーネルは、特徴量の内積になっている」という非自明な主張が証明されていることから導かれています。なので、直感で受け入れ難いのは当たり前で、だからこそ「トリック」なんて名前が付いているんですね。多分。

カーネルトリックで嬉しいこと

では、カーネルトリックが使えることで何が嬉しいのでしょうか？それとも、単に非自明な結果を示すことが出来た数学者（？）の自己満足なのでしょうか？

ここで、「その２」で述べた非線形回帰の課題を思い出してみてください。

有効な非線形項がわからない事がある
考えられる非線形項が多すぎる

という事が課題に挙げられていました。

しかし、ここでカーネルトリックの主張を思い出してみてください。

正定値カーネルを1つ決めて、それを回帰に利用することで、特徴写像を決めることなく、何かしらの特徴写像で変換を行った回帰分析を行っていることになります。

…ということは。

いい感じの正定値カーネルさえあれば、いちいち特徴写像の構成に頭を悩ませる必要はないのではないでしょうか。

また、

正定値カーネルを使った計算の計算量はサンプルサイズにしか依存しないので、考慮している非線形項が莫大になったとしても、比較的高速に計算できるのではないでしょうか。

ということで、正定値カーネルを使うと、非線形回帰の問題点が色々と解決しそうな雰囲気があるのです。実際にカーネルトリックを応用して非線形回帰を行う手段は現在確立されており、これをカーネル法と呼んでいる訳ですね。

もちろん「いい感じの正定値カーネルとは何か」みたいな課題は残されています。これは次回のテーマにしたいと思います。

実装

最後に、正定値カーネルを使ったロジスティック回帰を簡単に実装してみましょう。今回は正定値カーネルとして、多項式カーネルを使うことにします。

$\displaystyle K({\bf x}, {\bf y}) = {( {\bf x}^T {\bf y}+ 1 )}^d$

なお、データ加工などは省略しますのでご了承ください。（その２までと同様です）

カーネルの定義

## 多項式カーネルの定義
## x, y:n次元ベクトル(nは任意だが、xとyで同じ次元である必要がある)
## d:累乗の次元(default:d=2)
k_P <- function(x, y, d = 2) {
    (t(x) %*% y + 1)^d
}

カーネル自体は非常に単純なので、これで終わりです。

グラム行列の計算

次に、訓練データに対するグラム行列の計算をします。グラム行列というのは $K({\bf x_i}, {\bf x_j})$ を(i, j)成分に配置した行列です。その２では内積からグラム行列を計算していたのでベクトルのかけ算をするだけでしたが、今回はそうもいかないので素直に計算します。

「Rではforループ使うな」と言われているので*4頑張って工夫しましたが、同じ計算であれば入れ子ループした方がわかりやすいと思います。

## グラム行列の計算
K <- foreach(i = 1:NROW(X_train), .combine = "rbind") %do% {
    apply(X_train, 1,  k_P, x = X_train[i, ])
}
row.names(K) <- NULL # foreachによるインデックスを消去

損失関数の計算

その２と同じです。コードだけ再掲しておきます。

## 損失関数
## c:係数ベクトル、訓練データのサイズと同じ次元
lossfunc <- function(c) {
    require(foreach)
    
    loss <- foreach(i = 1:nrow(K), .combine = "sum") %do% {
        p <- 1 / (1 + exp(-c %*% K[i, ])) # y = 1になる確率
        return(-y_train[i] * log(p) - (1 - y_train[i]) * log(1 - p))
    }
    
    return(loss)
}

最適化

ここもその２と同様です。ちょっと気合い入れて（？）やってみましたが、やっぱり完全分離可能であるからか収束しません。適当なところで打ち切っておきます。

## 損失関数を最小化するcの算出
## 完全分離可能だからか収束しないので、20回目で打ち止め(十分近いところまでは行ってる？)
opt <- optim(par = rep(0, nrow(K)), fn = lossfunc, method = "BFGS", control = list(trace = 1, REPORT = 1, maxit = 20)) 
print(opt)
c <- opt$par

# > print(opt)
# $par
#   [1] -0.228796482 -0.214240013 -0.218561861 -0.252296538 -0.255087166 -0.266654067 -0.166373633 # -0.265606251 -0.235332926 -0.175364375 -0.234709029
#  (中略)
# [144]  0.138411815  0.090848513  0.071024964 -0.004708963  0.031054236  0.086940162  0.144573074
# 
# $value
# [1] 0.02973776
# 
# $counts
# function gradient 
#      59       20 
# 
# $convergence
# [1] 1
# 
# $message
# NULL

精度検証

その２と同じように、決定領域を描いてみます。

f:id:skroot3:20171219235652p:plain

似たような絵ばかりで飽きてきたかと思いますが、今回もうまく分離できていることが確認できました。

まとめと次回予告

今回は、カーネルトリックについて説明しました。個人的には、カーネルトリックの非自明さ、気持ち悪さ、そしてそれゆえのカーネル法の素晴らしさをわかっていただければ満足です。

最後の方は前回とほぼ同様の実装だったので半分手抜きになりましたが、まあ次回また同じようなことをするので大丈夫でしょう（
今回はどちらかというと説明の方がメインです！

次回は、カーネルの選び方について説明します。
ここまでくれば実用的に使える程度の知識を身につけた、と言えるのではないでしょうか。たぶん。

*1:言い訳

*2:内積の対称性から容易に示すことができます

*3:非常に難易度の高い証明です。興味ある方は参考文献(その２参照)を読んでください。

*4:実は今ではそんなに遅くないらしいですね