カーネル法入門：ノート2 - クッキーの日記

以下の本を読みます。キャラクターは架空のものです。解釈の誤りは筆者に帰属します。お気付きの点がありましたらコメント等でご指摘いただけますと幸いです。

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

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

以前カーネル法を勉強したときの記事：パターン認識と機械学習下（第6章：その6〈終〉） - クッキーの日記

前回：ノート1 ／次回：まだ

f:id:cookie-box:20190101155733p:plain:w60

前回のあらすじをまとめると、

データを分析するには（つまり、ある変数とある変数が依存し合っているかを調べたり、それに基づいて回帰したり分類したりするには）、線形な依存性のみを調べるのでは不十分なこともあるので、非線形な依存性も調べたい。
言い換えると、データを $\Phi : \Omega \to H$ なる特徴写像で特徴空間 $H$ に写像することで非線形な特徴も引き出して、この特徴空間で分析をしたい。
しかし、あらゆる高次項を特徴に追加したのでは、組み合わせが爆発してしまい、主成分分析のために標本共分散行列の固有値分解をしたり、回帰分析のためにデータ行列の擬似逆行列を求めたりすることができない。
が、少なくとも特徴空間における主成分分析は、行列の各成分が $k(X_i, X_j) = \langle \Phi(X_i), \Phi(X_j) \rangle$ であるような行数と列数がデータサイズの（有限の大きさの）行列の2乗の固有値分解をすることで達成できる。ので、 $k(X_i, X_j)$ の値さえ評価できればよいことになる。ここでポイントとなったのは、

特徴空間における第1主軸 $f$ は $\displaystyle f = \sum_{i=1}^N a_i \Phi(X_i)$ の形でかける（もしかけない成分があったら、かけない成分をゼロにしてその分かける成分を引き伸ばした方が各データの $f$ との内積の分散が大きくなる）。
そうすると解くべき最大化問題の中にデータは $k(X_i, X_j)$ の形でしか出てこない。

…普通に考えれば $\Phi(X_i)$ の標本分散共分散行列を固有値分解しなければならないところを、 $f$ は $H$ の部分空間 $H_0$ から探せば十分ということに着目して $f$ を表現し、最適化問題から $\Phi(X_i)$ を排除してしまうとは、なかなかトリッキーですね…。

f:id:cookie-box:20190101160814p:plain:w60

カーネルトリックというくらいだからね。

しかし、 $f$ は $H_0$ から探せば十分、というのは通常の主成分分析でも成り立つ話ですよね？通常の主成分分析の言葉でいうと、第1主軸 $u$ を求めるのに、 $u$ をあらゆる向きの $m$ 次元単位ベクトルから探すのではなくて、 $\displaystyle u = \sum_{i=1}^N a_i X_i$ でかけるものを探せば十分、ということになりますよね。

まさしくそうだけど、通常の主成分分析では固有値分解すべき標本分散共分散行列が $m \times m$ の大きさだったのが、カーネルPCAとして解くと固有値分解すべき行列が $N \times N$ の大きさになるから、どちらがありがたいのかという話だね。

あー、元データの次元数 $m$ よりデータサイズ $N$ の方が大きければかえって計算量が大きくなってしまいますね…。8ページからは主成分分析からまた変わって、リッジ回帰という例がありますね。リッジ回帰とは、線形モデルを最小2乗フィッティングするのにパラメータベクトルのL2ノルムの2乗に比例する正則化項を付けて、「誤差2乗和を小さくするパラメータベクトルを求めよ。ただし、パラメータベクトルのL2ノルムの2乗は小さく抑えること」といったもののようですね。その解 $\hat{a}$ は…「容易に求められるように」？何気ない「容易に」が読者を傷つけるのでは？

本だと尺の都合があるからね。最小化したいのは以下だよね。

$\displaystyle \sum_{i=1}^N || Y_i - a^\top X_i ||^2 + \lambda || a ||^2$ これを $a$ で偏微分しようか。ベクトルで微分するときの公式をつかうよ。 $\begin{split} \displaystyle \frac{\partial}{\partial a} \left( \sum_{i=1}^N || Y_i - a^\top X_i ||^2 + \lambda || a ||^2 \right) &= 2 \sum_{i=1}^N ( Y_i - a^\top X_i ) X_i + 2 \lambda a \\ &= 2 \sum_{i=1}^N Y_i X_i - 2 \sum_{i=1}^N (a^\top X_i) X_i + 2 \lambda a \\ &= 2 X^\top Y - 2 X^\top Xa + 2 \lambda a \end{split}$ あとはこれが $0$ になるように $a$ を解けばいい。8ページ中ほどの式になるのはわかるんじゃないかな。

まあそこまで説明してもらえばわかります。その式の下の、リッジ回帰はどんなときに利用されるかというので、「 $X^\top X$ が退化していたり」というのは何ですか？

退化というのは $X^\top X$ が正則でない、逆行列をもたないという意味だね。もし今の回帰問題で正則化項がなかったら、つまり、 $\lambda = 0$ だったら、回帰係数の解は $\hat{a} = (X^\top X)^{-1} X^\top Y$ になるけど、 $X^\top X$ が逆行列をもたなかったらこれは求まらない。

確かに $X^\top X$ が逆行列をもっていなかったら $\hat{a}$ は求まりませんね。いえでも、解けないというのはどういうことなんでしょう？それに、解けないから正則化項 $\lambda || a||^2$ を導入して解決、というのもどういうことなのか…。

$X^\top X$ が正則でないというのは、 $X^\top X$ の階数が $m$ に満たない、さらに言い換えると $X$ の階数が $m$ に満たないって状況だね。そうなっちゃうのは以下のようなときかな。

データ数 $N$ を上回るほど、データの次元 $m$ が巨大であるとき。
データセットが線型従属なとき（例えば、同じ点が重複して含まれているとか）。

こんな状況だと、誤差2乗和を最小にする $\hat{a}$ は無数に出てきちゃうんだよね。そうだな…以下の PRML の 146 ページに Figure 3.4 があるよね。

Christopher Bishop at Microsoft Research

Figure 3.4 の左側の図を例にすると、これはちゃんと $X^\top X$ が退化していないケースで、誤差2乗和（青い等高線）が最小になる点があるんだよね（青い点）。誤差2乗和はこの青い点を底にしたお椀になってる。でも、 $X^\top X$ が退化しているとこれみたいなお椀にならない。スノボのハーフパイプみたいにずーっと谷底が続いてたり、最悪まっさらな平面で「全体的に底」って感じになっちゃう。そうなると最小点を選べなくて困る。その解決策として、「じゃあ別のお椀をもってくればいいじゃん」というのが正則化項だ。Figure 3.4 の左側の図の赤い等高線が原点を底にした「L2ノルムのお椀」だね。このお椀を重ね合わせれば最小点が求まるって寸法だね（※）。もちろん、元々の $X^\top X$ が退化していない場合でも正則化のお椀を重ねていい。Figure 3.4 は誤差2乗和の青いお椀と正則化の赤いお椀を重ねて、重ねてできた新しいお椀の底になるのが $w^\ast$ ってことだね。

※	厳密には、もし $X^\top X$ が固有値 $-\lambda$ をもっていたら、正則化項を導入しても $X^\top X + \lambda I_N$ は正則になりません。が、ピンポイントでそうなることはあまりないのではないかと思っています。

なるほど…誤差2乗和のみでは回帰係数が選びきれないときでも選ぶ基準を新たに与えたんですね。原点を底にしたお椀ということは「選べなかったら原点になるべく近い解にしといて」って感じですよね…？

そうだね。ちなみに、L2ノルムのお椀は Figure 3.4 の左側の赤い等高線のように丸いけど、L1ノルムのお椀はその右側のようにひし形をしてるね。お椀っていうかひっくり返した四角錐だね。だから、L1ノルムで正則化するラッソ回帰では一部の成分がゼロになるような解が求まりやすい。この図でも $w_1 = 0$ の点が解になってるね。

じゃあラッソ回帰は「なるべく平面 $w_1 = 0, \; w_2 = 0, \; \cdots$ に貼り付く解にしといて」って感じですかね。…まあラッソ回帰は置いておいて、今回はリッジ回帰を特徴空間でやりたいわけです。特徴空間でのリッジ回帰の目的関数は式 (1.9) ですね。また特徴空間 $H$ での内積とノルムが出てきていますね…。

データ空間でやっていた「 $X_i$ の各成分を回帰係数 $a$ で足し合わせてそれを予測値としようとする」というのは、特徴空間 $H$ でやろうとすると「 $\Phi(X_i)$ とベクトル $f$ との内積をとってそれを予測値をしようとする」ということだからね。

しかし $\Phi(X_i)$ というベクトルはいま得体の知れないもので、カーネルPCAのときは最終的に $\Phi(X_i)$ は排除されたのですよね。リッジ回帰ではどうなるのでしょう。リッジ回帰でも解は $\displaystyle f = \sum_{j=1}^N c_j \Phi(X_j)$ の形でかけるとありますね。

もし $f$ に直交補空間の成分 $f_\bot$ があったらそれを丸々削った方が正則化項の $|| f||^2$ が抑えられるからね。

全くですね。ということは…今回もまた最小化したい式の中にデータは $k(X_i, X_j)$ の形でしか出てきませんね。 $\Phi(X_i)$ は出てきません。

最小化したい式の中に出てこないだけじゃなく、予測値を計算するときにも $\Phi(X_i)$ が出てこないことが肝心だね。 $\Phi(X_i)$ を直接取り扱いたくないから。でも未知の点 $x$ に対して「 $\Phi(x)$ とベクトル $f$ との内積をとるということをするとやっぱりデータは内積 $k(X_j, x) = \langle \Phi(X_j), \Phi(x) \rangle$ の形でしか出てこない。

えっと、まとめると、特徴空間でそのままリッジ回帰する（予測したい未知データの特徴 $\Phi(x)$ と内積を取るべき特徴空間内のベクトル $f$ そのものをもとめる）のではなく、 $f$ が $\displaystyle f = \sum_{j=1}^N c_j \Phi(X_j)$ の形でかけることに注目して $c_j$ を求める問題にしたのですよね。うーん、カーネルPCAでもカーネルリッジ回帰でも「解はこの形でかけるじゃん」というところがどうも降って湧いたようなんですが…。

最適化問題中の $\Phi(X_j)$ を意地でも $\langle \Phi(X_i), \Phi(X_j) \rangle$ の形にするために、 $f$ から $\Phi(X_i)$ を絞り出したって考えればそんなに降って湧いてもないんじゃないかな。

絶対に $\Phi(X_i)$ を直接取り扱いたくないのだという強い意思といったところですか…。ところで、カーネルPCAとカーネルリッジ回帰ではデータを特徴写像 $\Phi$ で特徴空間 $H$ にとばしましたが、 $\Phi$ や $H$ の詳細には立ち入っていませんでした。 $H$ についてわかっているのは、ベクトル空間であることと、内積 $\langle \cdot, \cdot \rangle$ をもつということだけですね。 $\Phi$ もどんな写像なのかまるでわかりません。直接取り扱っていないので。常にカーネル関数 $k$ を通して取り扱っています。…この $H$ と $\Phi$ は適当に決めていいんでしょうか。

ではないね。まずカーネル関数 $k(X_i, X_j) = \langle \Phi(X_i), \Phi(X_j) \rangle$ の値が評価できないといけないからこれが発散するような $\Phi$ だと困ると思うし…そもそも私たちは直接 $\Phi$ を決めないし…まあ [tex\Phi] とは何かは次章だね。

そもそも非線形な依存性を調べたくて $\Phi$ などという特徴写像を持ち出したのですから、どのような $\Phi$ が可能なのかが肝心な気がしますが…。ところで今回カーネルリッジ回帰を学びましたが、python の scikit-learn でカーネルリッジ回帰ができるようなので実行してみます。iris データの全品種ごちゃ混ぜで、がく片の長さからがく片の幅を予測するという何がしたいのかわからない例ですが…カーネル関数 $k$ には色々なものが指定できるし、自分で実装したものを渡すこともできそうです。組み込みのものを利用するときはおそらく以下の pairwise.kernel_metrics にあるものが指定できると思います。以下のスクリプトではいじっていませんが、パラメータも変更できます。

sklearn.kernel_ridge.KernelRidge ‐ scikit-learn 0.20.3 documentation

sklearn.metrics.pairwise.kernel_metrics ‐ scikit-learn 0.20.3 documentation

f:id:cookie-box:20190328234405p:plain:w480

もともとが線形関数でも非線形関数でもフィッティングできなさそうなのでそもそも例として適当ではありませんが、線形フィッティングは駄目そうですね。 $k$ にpolyカーネル $k(x, y) = (x^\top y)^3$ を選んだときよりもrbfカーネル $k(x, y) = \exp(-\gamma || x-y ||)^2$ を選んだときの方が点の密度が大きいところを通れているといえばいるような。

（ノート3があれば）つづく

スクリプト

import numpy as np
from sklearn.kernel_ridge import KernelRidge
from sklearn.datasets import load_iris

%matplotlib inline
import matplotlib.pyplot as plt
from pylab import rcParams
rcParams['figure.figsize'] = 8, 6
rcParams['font.size'] = 16

iris = load_iris()
#x = np.array(iris['data'][:,0][iris['target']==0]) # setosa の sepal length (cm)
#y = np.array(iris['data'][:,1][iris['target']==0]) # setosa の sepal width (cm)
x = np.array(iris['data'][:,0]) # sepal length (cm)
y = np.array(iris['data'][:,1]) # sepal width (cm)

print(x.shape)
print(y.shape)

x = x.reshape(-1, 1)
y = y.reshape(-1, 1)

print(x.shape)
print(y.shape)

# ----- オリジナル -----
plt.scatter(x, y, label='original')
# ----- linearカーネル -----
kr = KernelRidge(kernel='linear')
kr.fit(x, y)
plt.scatter(x, kr.predict(x), label='predict(linear)')
# ----- polyカーネル -----
kr = KernelRidge(kernel='poly')
kr.fit(x, y)
plt.scatter(x, kr.predict(x), label='predict(poly)')
# ----- rbfカーネル -----
kr = KernelRidge(kernel='rbf')
kr.fit(x, y)
plt.scatter(x, kr.predict(x), label='predict(rbf)')

plt.legend()
plt.show()

(150,)
(150,)
(150, 1)
(150, 1)