カーネル法入門: ノート2

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

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

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

以前カーネル法を勉強したときの記事: パターン認識と機械学習 下(第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) の標本分散共分散行列を固有値分解しなければならないところを、fH の部分空間 H_0 から探せば十分ということに着目して f を表現し、最適化問題から  \Phi(X_i) を排除してしまうとは、なかなかトリッキーですね…。

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

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

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

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

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

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

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

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

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

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

 \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ページ中ほどの式になるのはわかるんじゃないかな。

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

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

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

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

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

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

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

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

  • データ数 N を上回るほど、データの次元 m が巨大であるとき。
  • データセットが線型従属なとき(例えば、同じ点が重複して含まれているとか)。
こんな状況だと、誤差2乗和を最小にする \hat{a} は無数に出てきちゃうんだよね。そうだな…以下の PRML の 146 ページに Figure 3.4 があるよね。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 は正則になりません。が、ピンポイントでそうなることはあまりないのではないかと思っています。
f:id:cookie-box:20190101155733p:plain:w60

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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)