線形な手法とカーネル法（回帰分析） #機械学習 - Qiita

線形な手法とカーネル法（回帰分析） #機械学習 - Qiita 442 いいねしたユーザー一覧へ移動 358 X（Twitter）でシェアする Facebookでシェアする はてなブックマークに追加する more_horiz 記事を削除する close 一度削除した記事は復旧できません。 この記事の編集中の下書きも削除されます。 削除してよろしいですか？ キャンセル 削除する delete info この記事は最終更新日から3年以上が経過しています。 5日目 東北大学 計算機科学研究会 Advent Calendar 2017 @ wsuzume ( Yoshinobu Ogura ) 線形な手法とカーネル法（回帰分析） 機械学習 線形回帰 機械学習入門 カーネル法 多重線形回帰 442 最終更新日 2022年05月02日 投稿日 2017年12月06日 まえがき 線形な手法とはなんぞや。カーネル法とはなんぞや。機械学習について学んだことがある人なら言葉くらいは聞いたことがあるかもしれません。 線形な手法とは、ざっくり言えば一次の連立方程式を解くだけで答えが求まるように問題を設計することです。つまり逆行列さえ計算できれば問題が解けるということなので、とても簡単です。 ただ、線形な手法で解ける問題は限られています。たとえばグラフ上の点を直線で近似するようなことしかできません。けれど世の中そんなまっすぐじゃないですよね。もっとクネクネしたい。クネクネしているのを「非線形」といいますが、非線形な手法は解くのが非常に難しい場合が多いです。 そこで用いるのがカーネル法です。カーネル法は生データに非線形写像を施してやって、それを新しいデータとみなして線形な手法を適用してやることで問題を解きます。つまり「生データをクネクネした空間にうまいこと並べ直してあげれば、その空間の中ではスパッと一直線になってるんじゃないの？」というのがカーネル法の発想です。そしてクネクネした空間で引いた直線を、私たちの住むまっすぐな空間に持ってきてあげたら、きっとクネクネになっていますよね。やったぜ、ほとんど線形な手法を使ってクネクネできた！ カーネル法を使うと、線形な手法で培ってきた知識をほとんどそのまま生かして非線形なデータを扱うことができます。また、カーネル法はいま流行りの深層学習やニューラルネットワークに比べて数学的に整備されており、なぜ上手くいくのかといった点にブラックボックスがありません。 今回は線形回帰と多重線形回帰についてじっくり解説したあと、カーネル回帰へと話を展開します。 注意 この記事では列ベクトルを${\bf x}$と表記し、行ベクトルは${\bf x}^{\mathrm{T}}$のように転置記号を用いて表記します。 今回用いるサンプルコードは こちら です。Google Colaboratory にアップロードしてあるので、Googleアカウントをお持ちであれば実行できるかと思います。 2020/11/28追記 当記事のコメント欄で議論した、より理論的な側面に踏み込んだ内容を整理した記事を公開しました。 カーネルリッジ回帰 線形回帰 ものすごい今さら感がありますが、まずは線形回帰について解説します。 そもそも回帰ってなんぞや。統計とか機械学習の分野で回帰と言ったら回帰分析のことを指します。回帰分析とは、実際に観測されたデータの分布を見て、それを再現するような関数を探し出してあげよう、という方法です。回帰分析によって関数を見つけてあげれば、実際に観測していないデータについてもどんな値を取るか予測できるようになります。データを直線に回帰する場合、これを線形回帰といいます。 どういうものかは実際に見たほうが早いです。下のように、いかにも直線っぽいデータが観測されたとします。 たとえば縦軸を電圧$V$、横軸を電流$I$として物質の電気抵抗$R$を測ろうとしたらこういうグラフが出てきますよね。理論的には$V=RI$なので綺麗に直線になっていてほしいのですが、電気抵抗に限らず実際のデータはなんかよくわからない理由によってゆらぎます。ここから先は議論を一般化するため、縦軸を$y$、横軸を$x$として直線$\hat{y}=f(x)=w_0 + w_1x$の切片$w_0$と傾き$w_1$を求めよ、ということにします（回帰した結果としての出力にはハット記号をつけて$\hat{y}$のように表記し、実際に観測された値$y$と区別することが多いです）。 私たちは$N$個のデータを観測したとします。すなわち$x$に対応した$y$の値のデータの組を、 (x^{(1)}, y^{(1)}), (x^{(2)},y^{(2)}), \ldots, (x^{(N)},y^{(N)}) のように集めてきました。観測されたデータに含まれるゆらぎを考慮しつつも、それっぽい直線を引きたい。じゃあどんな直線が一番それっぽいの？ということを表すために誤差関数$R({\bf w})$というものを定義します（念のため、抵抗$R$とはまったく関係ないです）。 誤差関数$R$の選び方にはいろいろあるのですが、とりあえず二乗誤差$r_{square}$（残差平方和$RSS$ともいいます）を用いることを考えて見ましょう。二乗誤差とは、実際に観測されたデータと回帰直線によって予測された値との差の2乗を取って、それをすべての点について総和したものであり、$w_0$と$w_1$の関数です（${\bf w}=(w_0, w_1)^{\mathrm{T}}$とおけば、$R$はベクトル${\bf w}$の関数です）。 R({\bf w}) = r_{square}(w_0,w_1)= \sum_{i=1}^{N}\left\{y^{(i)}- f(x^{(i)})\right\}^2 = \sum_{i=1}^{N}\left\{y^{(i)}- (w_0 + w_1 x^{(i)})\right\}^2 つまりこの誤差関数$R({\bf w})$を最小化する${\bf w}$を見つけてやれば、二乗誤差という尺度のもとでもっともそれっぽい直線を引けるというわけです。見つけ方は簡単です。$R({\bf w})$は$w_1,w_2$に関する二次式になっているので、偏微分して$0$になる場所が極値になります（$R({\bf w})$は0以上の範囲にある二次関数なので、極値で最小となります）。ここで、 {\bf y} = \left(\begin{matrix} y^{(1)} \\ y^{(2)} \\ \vdots \\ y^{(N)} \end{matrix}\right) ,\; X= \left(\begin{matrix} 1 & x^{(1)} \\ 1 & x^{(2)} \\ \vdots & \vdots\\ 1 & x^{(N)} \end{matrix}\right) とおくと、実は誤差関数は R({\bf w}) = ({\bf y} - X{\bf w})^{\mathrm{T}} ({\bf y} - X{\bf w}) と表すことができます（本当にそうなるか、一度確かめてみると面白いです）。合成関数の微分を用いれば、 \begin{align} \frac{\partial R({\bf w})}{\partial {\bf w}} &= \left(\begin{matrix} \frac{\partial R({\bf w})}{\partial w_0} \\ \frac{\partial R({\bf w})}{\partial w_1} \end{matrix}\right) = \left(\begin{matrix} - 2\sum_{i=1}^{N}1 \cdot \left\{y^{(i)}-(w_0+w_1x^{(i)})\right\} \\ - 2\sum_{i=1}^{N}x^{(i)}\left\{y^{(i)}-(w_0+w_1x^{(i)})\right\} \end{matrix}\right)\\ &= -2X^{\mathrm{T}}({\bf y} - X{\bf w}) ={\bf 0} \end{align} となりますので、これを解いて \begin{align} -2X^{\mathrm{T}}({\bf y} - X{\bf w}) ={\bf 0} \\ \therefore {\bf w} =(X^{\mathrm{T}}X)^{-1}X^{\mathrm{T}}{\bf y} \end{align} となります（ただし$X^{\mathrm{T}}X$に逆行列が存在するときに限ります）。 では実際にやってみましょう。実は、さっき図に示したいかにも直線っぽいデータはガウス分布にしたがった人工的なバラつきを与えて直線の周囲に生成したデータです。 import numpy as np import pandas as pd from matplotlib import pyplot as plt import seaborn as sns # 線形回帰モデル用データセット # a : slope of a line # b : intercept of a line # size : number of data # xlim : domain of variable x # scale : standard deviation def linear_dataset ( a , b , size , xlim = [ 0 , 1 ], scale = None ): x = np . random . uniform ( xlim [ 0 ], xlim [ 1 ], size ) y = a * x + b if scale is not None : noize = np . random . normal ( 0 , scale , size ) y = y + noize df = pd . DataFrame () df [ ' x ' ] = x df [ ' y ' ] = y return df # データを生成 data = linear_dataset ( 0.5 , 0.2 , 20 , scale = 0.03 ) # データを描画 plt . figure ( figsize = ( 6.472 , 4 )) plt . scatter ( data [ ' x ' ]. values , data [ ' y ' ]. values ) plt . xlim (( 0.1 , 1.05 )) plt . ylim (( 0.2 , 0.8 )) plt . show () 傾きは0.5,切片は0.2に設定してあります。上のコードを実行すると大体同じようなデータセットが生成されると思います。そしてこれを線形回帰にかけます。 def linear_regression ( X , Y ): X = np . stack (( np . ones ( X . shape [ 0 ]), X ), axis = 1 ) W = ( np . linalg . inv ( X . T . dot ( X )). dot ( X . T )). dot ( Y ) return W # 重みWを求める W = linear_regression ( data [ ' x ' ]. values , data [ ' y ' ]. values ) こうして求めたWを print(W) などを用いて表示すると求めた重みを直接覗くことができます。私が実行したときは [ 0.19135655 0.50836672] となりました。元々が切片0.2,傾き0.5の直線から生成したデータであることを考えれば、うまくいってますね。実際はこんな出来レースみたいな評価はよろしくないので、求めた直線を元のグラフに描画して確かめます。 X = np . arange ( 0 , 1.5 , 0.01 ) Y = W [ 0 ] + W [ 1 ] * X plt . figure ( figsize = ( 6.472 , 4 )) plt . scatter ( data [ ' x ' ]. values , data [ ' y ' ]. values ) plt . plot ( X , Y ) plt . xlim (( 0.1 , 1.05 )) plt . ylim (( 0.2 , 0.8 )) plt . show () 多重線形回帰 先ほどは入力$x$と出力$y$がいずれも1次元の場合を扱っていたので、直線$\hat{y}=w_0 + w_1x$に回帰させていました。入力がより高次元の${\bf x}=(x_1,x_2,\ldots,x_n)^{\mathrm{T}}$のn次元の入力に一般化された線形回帰を多重線形回帰といい、次の式で表します。 \hat{y} = f({\bf x}) = w_0 + w_1 x_1 + w_2 x_2 + \ldots + w_n x_n この式は行列を用いて \hat{y} = \left(\begin{matrix} 1 & x_1 & x_2 & \cdots & x_n \end{matrix}\right) \left(\begin{matrix} w_0 \\ w_1 \\ w_2 \\ \vdots \\ w_n \end{matrix}\right) と書き直すことができます。さらに出力が${\bf y} =(y_1,y_2,\ldots,y_m)^{\mathrm{T}}$の$m$次元の場合にも Y = \left(\begin{matrix} y_1^{(1)} & y_2^{(1)} & \cdots & y_m^{(1)} \\ y_1^{(2)} & y_2^{(2)} & & \vdots \\ \vdots & & \ddots & \vdots\\ y_1^{(N)} & \cdots & \cdots & y_m^{(N)} \end{matrix}\right) ,\; X = \left(\begin{matrix} 1 & x_1^{(1)} & x_2^{(1)} & \cdots & x_n^{(1)}\\ 1 & x_1^{(2)} & x_2^{(2)} & & \vdots\\ \vdots & & & \ddots&\vdots\\ 1 & x_1^{(N)} & \cdots & \cdots & x_n^{(N)} \end{matrix}\right) W = \left(\begin{matrix} w_{10} & w_{20} & \cdots & w_{m0} \\ w_{11} & w_{21} & & \vdots \\ \vdots & & \ddots & \vdots\\ w_{1n} & \cdots & \cdots & w_{mn} \end{matrix}\right) のように定義することで、 \hat{Y} = f(X)=XW \\ R(W) = (Y-XW)^{\mathrm{T}} (Y-XW) \\ \therefore W = (X^{\mathrm{T}}X)^{-1}X^{\mathrm{T}} Y として、一次元の場合と同様に解くことができます。 カーネル回帰 ちょっとだけ多重線形回帰を思い出してみましょう。多重線形回帰は$x_0=1$とおけば以下の式で表すことができるのでした（$\Sigma$を用いて書き直してあります）。 f({\bf x}) = \sum_{i=0}^{n} w_i x_i カーネル回帰では以下のような式で関数を近似します。 f({\bf x}) = \sum_{i=1}^{N} \alpha_i k({\bf x}^{(i)}, {\bf x}) 多重線形回帰によく似た式になっていますね。$k({\bf x}, {\bf x}')$はカーネル関数と呼ばれるもので、一定の条件さえ満たせば自由に設計することができます（この辺りを詳しく知りたい人は参考文献[1]をご覧ください）。今回はよく使われるガウスカーネルを用います。ガウスカーネルは次のように定義されます。 k({\bf x}, {\bf x}') = exp(-\beta \|{\bf x} - {\bf x}'\|^2) ただし$\|{\bf x} - {\bf x}'\|=\sqrt{\sum_{j=1}^n (x_j - x'_j)^2}$はベクトル同士の距離（$L^2$ノルム）を表しています。$\beta$は$\beta \gt 0$を満たす実数のハイパーパラメータであり、使う人が勝手に決めます（結果を見ながらいろいろいじくってみるのが普通です）。 もう少しカーネル回帰の式を詳しく考察してみましょう。入力が1次元の場合、ガウスカーネルは$k(x,x')=exp(-\beta(x-x')^2)$と書き直すことができます。$x'=0$とおいて$x$を動かしてガウスカーネルをプロットすると、下の図のようになります。 def kernel ( xi , xj , beta = 1 ): return np . exp ( - beta * np . sum (( xi - xj ) ** 2 )) X = np . arange ( - 1.5 , 1.5 , 0.01 ) Y = np . zeros ( len ( X )) for i in range ( len ( X )): Y [ i ] = kernel ( X [ i ], 0 , 10 ) plt . figure ( figsize = ( 6.472 , 4 )) plt . plot ( X , Y ) plt . show () $x'$に値を入れると、このガウスカーネルは左右に動きます。実際は観測されたデータの数だけ$x'$に値を代入して、その$N$個のガウスカーネルを線形に足し合わせたものがカーネル回帰になります。ガウスカーネルを5個くらい適当にズラしつつ適当に重みづけて足し合わせたのが下の図です。 def kernel ( xi , xj , beta = 1 ): return np . exp ( - beta * np . sum (( xi - xj ) ** 2 )) X = np . arange ( - 1.5 , 1.5 , 0.01 ) Y = np . zeros ( len ( X )) centers = [ - 1 , 0 , 0.5 , 0.6 , 1 ] weights = [ 0.2 , 1 , 0.3 , - 1 , 0.5 ] for i in range ( len ( X )): for weight , center in zip ( weights , centers ): Y [ i ] += weight * kernel ( X [ i ], center , 10 ) plt . figure ( figsize = ( 6.472 , 4 )) plt . plot ( X , Y ) plt . show () これ、うまいこと足し合わせればどんな関数でも近似できそうですよね（実際、無限個のガウス関数の足し合わせで、任意の非線形関数を任意の精度で近似することができます）。このうまい足し合わせ方を見つけてやるのがカーネル回帰です。 カーネル回帰の式はやはりベクトルを使って次のように書くことができます。 \hat{y} = \left(\begin{matrix} k({\bf x}^{(1)}, {\bf x}) & k({\bf x}^{(2)}, {\bf x}) & \cdots & k({\bf x}^{(N)}, {\bf x}) \end{matrix}\right) \left(\begin{matrix} \alpha_0 \\ \alpha_1 \\ \vdots \\ \alpha_N \end{matrix}\right) ここまで来れば、なぜ線形回帰についてしつこく丁寧に解説してきたか、察しがつくかと思います。 {\bf y} = \left(\begin{matrix} y^{(1)} \\ y^{(2)} \\ \vdots \\ y^{(N)} \end{matrix}\right) ,\; K=\left(\begin{matrix} k({\bf x}^{(1)}, {\bf x}^{(1)}) & k({\bf x}^{(2)}, {\bf x}^{(1)}) & \cdots & k({\bf x}^{(N)}, {\bf x}^{(1)}) \\ k({\bf x}^{(1)}, {\bf x}^{(2)}) & k({\bf x}^{(2)}, {\bf x}^{(2)}) & & \vdots \\ \vdots & & \ddots & \vdots \\ k({\bf x}^{(1)}, {\bf x}^{(N)}) & \cdots & \cdots & k({\bf x}^{(N)}, {\bf x}^{(N)}) \end{matrix}\right) ,\; {\bf \alpha} = \left(\begin{matrix} \alpha^{(1)} \\ \alpha^{(2)} \\ \vdots \\ \alpha^{(N)} \end{matrix}\right) とおきます（この$K$はグラム行列と呼ばれています）。グラム行列を用いるとカーネル回帰による二乗誤差$R({\bf \alpha})$は \begin{align} R({\bf \alpha}) &= r_{square}({\bf \alpha})= \sum_{i=1}^{N}\left\{y^{(i)}- f({\bf x}^{(i)})\right\}^2 =\sum_{i=1}^{N}\left\{y^{(i)}- \sum_{i=1}^{N} \alpha_i k({\bf x}^{(i)}, {\bf x})\right\}^2 \\ &=({\bf y}- K {\bf \alpha})^{\mathrm{T}}({\bf y}- K {\bf \alpha}) \end{align} と表すことができ、したがって {\bf \alpha} = (K^{\mathrm{T}}K)^{-1}K^{\mathrm{T}}{\bf y} = K^{-1}{\bf y} と解くことができます。線形回帰はせいぜい直線とかまっすぐな平面でしか近似することができませんでした。つまり次のようなデータセットに対して線形回帰を用いてもうまくいかないことはすぐに分かります。 # カーネル回帰モデル用データセット # size : number of data # xlim : domain of variable x # scale : standard deviation def wave_dataset ( size , xlim = [ 0 , 1 ], scale = None ): x = np . random . uniform ( xlim [ 0 ], xlim [ 1 ], size ) y = np . sin ( x ) if scale is not None : noize = np . random . normal ( 0 , scale , size ) y = y + noize df = pd . Da