カーネル法を勉強し直すその2　〜カーネル法を導くための準備〜

カーネル法勉強し直し企画の第二弾です。
今回から3回程度でカーネル法の仕組みについて解説します。

注意点として、数学的な詳しい議論は避ける予定です。*1もし数学的な詳しい議論を読みたい方は、以下の本を読んでみると良いかと思います。

カーネル法入門―正定値カーネルによるデータ解析 (シリーズ多変量データの統計科学)

作者: 福水健次
出版社/メーカー: 朝倉書店
発売日: 2010/11/01
メディア: 単行本
クリック: 19回
この商品を含むブログ (10件) を見る

今更ですが、本連載もこの本の1章〜3章を参考にさせていただいています。

back:カーネル法を勉強し直すその1 〜導入〜 - root3のメモ帳
next:カーネル法を勉強し直すその3 〜カーネルの導入とカーネルトリック〜 - root3のメモ帳

ロジスティック回帰で非線形分類をする
非線形項入りロジスティック回帰の問題点
- 有効な非線形項がわからない事がある
- 考えられる非線形項が多すぎる
非線形項の追加を別の視点で考察する
別の視点でロジスティック回帰を解く
- 実装
- 精度検証
結局何が違ったのか
まとめと次回予告

ロジスティック回帰で非線形分類をする

いきなりこんな見出しで「前回と言ってること違うじゃねーか！」とツッコまれそうですが、実はロジスティック回帰でも非線形分類は可能です。発想は非常に明快かつ単純です。非線形項を自分で与えてあげれば良いのです。*2

まずは前回も提示したサンプルデータをもう一度出します。

f:id:skroot3:20171018220853p:plain

このプロットをみて、「同心円状に同じクラスが分布しているな」という事はパッとみてわかるかと思います。そこからさらに言い換えをすると、

同心円状に同じクラスが分布している
↓
クラスが中心からの距離に依存している
↓
中心からの距離 $\sqrt{{x_1}^2+{x_2}^2}$ を変数として組み込めば良さそう

という事になります。なお、今回の問題ではあくまで2クラスを分類できれば良いので、距離 $\sqrt{{x_1}^2+{x_2}^2}$ そのものではなく、距離の二乗 ${x_1}^2+{x_2}^2$ を利用する事にします。

実際にやってみる前に、今回も学習データと検証データに分割しておきます。まずは必要なパッケージとデータの読み込みから。
※前回は導入だったのでコード紹介はところどころ省略していましたが、今回からは基本的にほぼ全て掲載します。

library(dplyr)
library(magrittr)
library(ggplot2)

library(caret)

# loading
data = read.csv("data/sample_data.csv") %>% 
    mutate(class = as.factor(class))

ちなみに、dataは以下のような構造で入っています。

head(data)
#           x1         x2 class
# 1  0.5748705  0.6174553     0
# 2 -0.7772104 -0.6852452     0
# 3 -0.8842989 -0.4700317     0
# 4 -0.8157662 -0.7874295     0
# 5  0.6256141 -0.7546584     0
# 6 -0.5796135 -0.6355427     0

次に、データの分割をします。caretパッケージのcreateDataPartitionという関数を使います。また、分割前は必ずseedを固定しておきましょう。固定しておかないと、再現性がなくなってしまいます。seedの値はお好みで構いません。*3

# train_test split
set.seed(1732)

train_rate <- 0.75 # 学習データの割合
train_ind <- createDataPartition(data$class, p = train_rate, list = FALSE, times = 1)

train_indは学習に使うデータのインデックスを表しています。例えば、train_ind = c(1, 3)であれば、1番目と3番目のデータを学習データとして使う、といった具合です。

非線形項の追加

では実際に非線形項である距離の二乗 ${x_1}^2+{x_2}^2$ を追加した上で、ロジスティック回帰を行います。前回の最初にやったロジスティック回帰と異なるのは、説明変数に距離の二乗が追加されている事のみです。

# 特徴量の構成
# x1^2 + x2^2という変数を付け加える
data %<>% mutate(z = x1^2 + x2^2) 

head(data)
#           x1         x2 class         z
# 1  0.5748705  0.6174553     0 0.7117271
# 2 -0.7772104 -0.6852452     0 1.0736170
# 3 -0.8842989 -0.4700317     0 1.0029144
# 4 -0.8157662 -0.7874295     0 1.2855197
# 5  0.6256141 -0.7546584     0 0.9609022
# 6 -0.5796135 -0.6355427     0 0.7398663

zという名前で、距離の二乗を追加しました。

距離の二乗を追加したところで、先ほどのtrain_indを用いで学習データと検証データに分けておきましょう。

data_train <- data[train_ind, ]
data_test <- data[-train_ind, ]

モデリング

次に、このデータを用いてロジスティック回帰を行います。式の書き方は前回と同様ですが、zを右辺に追加する必要がある事にだけ注意してください。*4

mod_glm <- glm(formula = class ~ x1 + x2 + z, data = data_train, family = "binomial")
print(mod_glm)

# Call:  glm(formula = class ~ x1 + x2 + z, family = "binomial", data = data_train)
# 
# Coefficients:
# (Intercept)           x1           x2            z  
#     -68.029        1.617        2.006       29.207  
# 
# Degrees of Freedom: 149 Total (i.e. Null);  146 Residual
# Null Deviance:       207.9 
# Residual Deviance: 2.538e-09     AIC: 8

※上記データでglmを実行するとおそらく警告メッセージが出てきますが、完全分離可能な故のメッセージです。今回においては特に害はありません。

標準化していないので厳密な事は言えませんが、zの係数が非常に大きく、有効に働いている事がわかります。また、AICが8と、非常に小さな値になっています。前回の単純なロジスティックモデルはAIC=209.6でした。

精度検証

このモデルの精度を求めるために、確率0.5を閾値とした時の混合行列を出してみます。caretパッケージの、confusionMatrixという関数を使います。

pred_glm <- if_else(predict(mod_glm, newdata = data_test) >= 0, 1, 0)
confmat_glm <- confusionMatrix(data = pred_glm, reference = data_test$class)
print(confmat_glm)

# Confusion Matrix and Statistics
# 
#           Reference
# Prediction  0  1
#          0 25  0
#          1  0 25
#                                     
#                Accuracy : 1      
# （以下略）

精度100%。テストデータに対しても、完全分離できています。zの有無だけでこれほどの違いが生じるという事がおわかりいただけたかと思います。

最後に、決定領域もプロットしてみましょう。前回と同じようにmeshを作成し、プロットします。これに関しては前回と同様ですので、結果のみ表示します。

f:id:skroot3:20171025203705p:plain

綺麗な円形の境界ができています。

非線形項入りロジスティック回帰の問題点

さて、ロジスティック回帰でも非線形要素を取り入れた分類ができました。めでたし、めでたし。

……ではもちろん終わりません。これでめでたくなれるほどデータ分析は甘くはないでしょう。次のステップへ進む前に、このやり方の問題点を整理しておきます。

有効な非線形項がわからない事がある

先ほどの例では、プロットすれば一目でわかるような特徴がありました。しかし、一般的なデータ分析では特徴がすぐにわかるとは限りませんし、むしろわからない事の方が多いでしょう。そもそもプロットすらできないかもしれません。

考えられる非線形項が多すぎる

「わからないなら考えられるもの全て入れてしまえ」という発想もあるかもしれません。では、実際に入れてみたらどうなるでしょうか。
仮に変数の数は2つ、非線形項も多項式のみを考えたとします。この場合、

2次まで： $x_1, x_2, {x_1}^2, {x_2}^2, x_1x_2$ の5個
3次まで： $x_1, x_2, {x_1}^2, {x_2}^2, x_1x_2, {x_1}^3, {x_1}^2x_2, x_1{x_2}^2, {x_2}^3$ の9個
4次まで：14個

といったように、変数の数はどんどん増えていきます。さらに変数の種類が増えればさらに爆発的に増えます。それらを全て含めて計算するのは非常にコストが高いです。

自分で明示的に非線形項を与える事は現実的ではない、という事がわかりました。
そして、前回書いたように、その問題を解決してくれるのがカーネル法です。ここからようやく、カーネル法の仕組みに触れていきたいと思います。*5

非線形項の追加を別の視点で考察する

カーネル法そのものに入る前に、「非線形項を追加する」という行為を、別の視点で考察してみようと思います。しばらく再び天下り式になりますが、お付き合いください。

特徴写像の定義

次の式は、ロジスティック回帰の確率を推定するための式です。

$\displaystyle \ln\frac{p}{1-p} = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \beta_3 ({x_1}^2+{x_2}^2)$

今回は距離の二乗もあったので、変数の数を3つにしています。

これまでやってきたように、距離の二乗を追加することで、変数の数を増やしています。これを、特徴量の抽出というニュアンスで捉えてみます。すなわち、次のように考えます。

元のデータ $(x_1, x_2)$ から、回帰に使う特徴量群 $(1, x_1, x_2, {x_1}^2+{x_2}^2)$ に変換している

元データから特徴量への変換法則を写像だと考えることで、特徴写像が定義されます。また、写影した後の空間を特徴空間と呼ぶ事にします。具体的に書くのであれば、今回の場合は次のような2次元ベクトル(in元の空間)から4次元ベクトル(in特徴空間)への写像です。

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

「なぜこのような言い換えをしておくと良いのか」を、この段階で理解する事は難しいかもしれません。今の所は「数学おきまりの定式化」だと思ってください。

内積表示

これ以降は、距離の二乗を明示的に使わずに特徴写像の形で表記したいと思います。

もう一度確率を推定する式を思い出してみます。今回は、特徴写像で写影した後の成分で考えてみましょう。すると、以下のような表記になります。

$\displaystyle \ln\frac{p}{1-p} = \theta_1 \Phi_1({\bf x})+ \theta_2 \Phi_2({\bf x}) + \cdots + \theta_n \Phi_n({\bf x})$

ここで、表記を簡潔にするために説明変数のベクトルを ${\bf x}$ と表記しています。また、 $\Phi_i({\bf x})$ はi番目の特徴量です。係数についても、元の $\beta$ と異なることを明示するために別の記号に置き換えました。

この式は、内積を利用することでよりシンプルな表記となります。 $\displaystyle \ln\frac{p}{1-p} = \langle {\bf \theta}, {\bf \Phi(x)}\rangle$

これをpについて解くことで、pの推定式は以下のようになります。 $\displaystyle p = \frac{1}{1 + \exp{(-\langle {\bf \theta}, {\bf \Phi(x)}\rangle})}$

問題の置き換え

最後の仕上げとして、係数パラメータ ${\bf \theta}$ の推定の問題を別の問題に置き換えます。

仮に、 ${\bf \theta}$ が各サンプルデータの特徴量の線形和で表せるとしましょう。すなわち、

$\displaystyle {\bf \theta} = \sum_{j=1}^N c_j \Phi({\bf x_j})$

と表せると仮定します。 ${\bf \theta}$ も $\Phi({\bf x_j})$ も特徴空間の次元のベクトルだという事に注意してください。これを先ほど導いた内積表示で現れている内積に代入すると、

$\displaystyle \langle {\bf \theta}, {\bf \Phi(x)}\rangle = \sum_{j=1}^N c_j \langle {\bf \Phi(x_j)}, {\bf \Phi(x)}\rangle$

となります。これは見た目以上にすごい結果を主張をしています。というのも、ロジスティック回帰において最適化する関数(対数尤度)には、学習データに関わる項が、係数と学習データの特徴量の内積 $\langle {\bf \theta}, {\bf \Phi(x_i)}\rangle$ の形でしか現れないため、上記の式でそれらすべてを置き換える事で、特徴量の係数の推定問題を線形結合の係数 ${\bf c}$ の推定問題に置き換えることが出来ることを主張しているからです。

もちろん、最初においた「各サンプルデータの特徴量の線形和で表せる」という仮定が適切である必要があります。しかしこの仮定が適切であることも容易にわかります。なぜなら、 ${\bf \theta}$ の線形和以外の成分は、内積に影響を与えないからです。*6

別の視点でロジスティック回帰を解く

ここまでで、

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

ということをしました。先ほども述べたように、ロジスティック回帰の最適化する関数には、学習データが関わる項は内積の形でしか現れないため、 ${\bf c}$ を推定すれば良い事になります。具体的な形は長くなってしまうので掲載しませんが、元々が $\langle {\bf \theta}, {\bf \Phi(x_i)}\rangle$ という項からなる関数であることから、

$\displaystyle \langle {\bf \theta}, {\bf \Phi(x_i)}\rangle = \sum_{j=1}^N c_j \langle {\bf \Phi(x_j)}, {\bf \Phi(x_i)}\rangle$

という形の式からなる被最適化関数になることはなんとなーくわかるでしょう。つまり、特徴量の内積さえわかってしまえば良いわけです。

実装

さて、数式ばかりになってしまったので、最後に実装してみましょう。ここはカーネル法における重要な部分なので、自分で実装してみます。内積がわかってしまえば良いということで、先に特徴量の内積を計算しておきます。

今回のサンプルデータにおいては、特徴写像が定数項(1)と距離の二乗を付け加えるだけのものです。距離の二乗はすでにデータに含まれているので、定数項のみ付け加えておきます。*7また、後々のために特徴量と目的変数は分離しておきます。

X_train <- data_train %>% 
    mutate(const = 1) %>% 
    select(-class) %>% 
    as.matrix

y_train <- data_train %>% 
    select(class) %>% 
    as.matrix %>% 
    as.numeric

X_test <- data_test %>% 
    mutate(const = 1) %>% 
    select(-class) %>% 
    as.matrix

y_test <- data_test %>% 
    select(class) %>% 
    as.matrix %>% 
    as.numeric

特徴量の内積は、行列の積を利用して容易に求められます。

K <- X_train %*% t(X_train)

次に被最適化関数を定義します。ロジスティック回帰では対数尤度を最大化したいのですが、Rの仕様によりデフォルトは最小化なので、負の対数尤度を定義します。

### 損失関数の定義
### y:観測した目的変数(1 or 0)
### K:内積の行列
lossfunc <- function(c) {
    require(foreach)
    
    loss <- foreach(i = 1:dim(K)[1], .combine = "sum") %do% {
        p <- 1 / (1 + exp(-c %*% K[i, ])) # y=1になる確率
        -y_train[i] * log(p) - (1 - y_train[i]) * log(1 - p)
    }
    
    return(loss)
}

最後に、optimを使って最適化しましょう。なお、最適化手法についてはここでは深く立ち入りません。*8完全分離可能な問題からうまく収束してくれませんので、いい感じのところで止めておきます。

opt <- optim(par = rep(0, 150), fn = lossfunc, method = "BFGS", control = list(trace = 1, REPORT = 1, maxit = 20))
c <- opt$par
print(c)

# [1] -0.448785638 -0.453780897 -0.412366529 -0.493923658 -0.510731854 -0.472685169 -0.395180071 -0.456999698 -0.525874473 -0.439156389 -0.499728880
# [12] -0.407851381 -0.341378967 -0.431868743 -0.427072459 -0.366841003 -0.356438580 -0.415034210 -0.457071461 -0.353744542 -0.462231882 -0.382534156
# （以下略）

さて、print(c)まで実行することで、係数が大量に表示されたと思います。今回学習に使っているデータは150個なので、cも150成分存在します。

精度検証

最後に、検証データに対する予測結果と、決定領域を求めておきます。予測関数を定義すれば、あとは前半と同様です。

# 予測用関数
# x : 特徴ベクトル
predfunc <- function(x) {
    K_newdata <- as.numeric(x %*% t(X_train))
    1 / (1 + exp(-c %*% K_newdata))
}

f:id:skroot3:20171025223854p:plain

最初の例と同様に、100%の精度かつ決定領域が円形になりました。*9

結局何が違ったのか

前半の解き方と後半の解き方で、得られた結果はほぼ同様でした。これは当然の結果で、なぜならあくまでも問題を解く際の視点を変えただけで、「与えられた特徴量から、2クラスを分離する」という問題自体は何も変わっていないからです。それであれば、なぜわざわざややこしい、面倒なやり方で解いたのでしょうか？

この解き方のメリットは、2つあります。

特徴空間の次元をいくら増やしても、学習に必要なデータは増えない
最終的に必要なのが、特徴量そのものではなく特徴量の内積である

1つ目について。前半の方法では、特徴量を増やせば増やすほど入力データが増加していきました。一方、後半の方法は、学習データ同士の内積があれば良いので、いくら次元が高くても必要な計算量はサンプルサイズにしか依存しません。*10

2つ目について。こちらに関しては「だから何？」と思うかもしれません。普通に考えれば、特徴量の内積を計算するためには、具体的な特徴量を経由するほかないはずです。

しかし、実は特徴量を経由せずに、特徴量の内積を計算する方法があります。もしそれが実現するのであれば、内積を計算するために特徴量の具体形を求める必要がなくなります。このことに注目したのが、カーネル法です。

まとめと次回予告

だいぶ長くなってしまいましたが、今回は、

特徴写像
内積表示
被最適化関数の置き換え

ということを主に説明しました。抽象的な議論ばかりでよくわからないところが多かったかと思います。こんな記事よりも丁寧に解説している良書・良サイトはたくさんありますので、そちらを読んでみてください（丸投げ）

次回は、ついにカーネルを導入します。そして、「カーネル法といえばこれ！」と言っても過言ではない、カーネルトリックについて説明します。

*1:僕自身、そんなに詳しくありません。一通り証明は読んでおきたいですが

*2:「それ非線形回帰(分類)じゃねーから」というツッコミも無しです(

*3:1732の由来は察してください

*4:そもそも「class ~ .」と書けば「class以外の全変数で回帰」の意味になるので間違いがないですが、zがある事を明示するためにあえて指定しています

*5:書いてて思いましたが、ここまで前回の内容でよかったですね

*6:直交補空間という概念が関わっています

*7:定数項が最後の列にありますが、気にしないでくて良いです。気になる人は並び替えてください

*8:立ち入れるほど詳しくない、それどころか無知of無知です

*9:実はちょっと閾値をいじっていますが、あくまでも同等の分類になるかどうかの検証なので許してください何でもしm

*10:もちろん特徴量の次元が増えれば内積の計算の負荷は高くなりますが、あくまで線形オーダーです