TL;DR
$Y,X$の間に線形な関係
$Y=Xw+\epsilon$
を仮定したとき、$w$を求めるための目的関数は、この記事での目的関数の最終形態として、
正則化係数も正則化項も(X,Y)にかかるノイズも全て成分ごと、サンプリングごとに異なる場合の重みつきリッジ回帰
$$\mathcal L=tr(M(Xw-Y)G(Xw-Y)^T)+tr(\Lambda(w-w_0)G(w-w_0)^T)$$
で表現できる。最小値をとるときの$w$は、
$$w=(X^TMX+\Lambda)^{-1}(X^TMY+\Lambda w_0)$$
となる。ここで、
- $X$は写像済み計画行列
- $Y$も計画行列
- $M$は各サンプルの分散(=誤差二乗)を対角成分においた対角行列の逆行列(精度行列)
- $G$はノイズ$\epsilon$が従う正規分布に対応する精度行列(分散共分散行列の逆行列)。回帰結果には直接影響しないが確率論的に無意味な値ではないため、目的関数上には残しておく。
- $\Lambda$は正則化項だが、確率論的にはウエイト$w$の事前分布が従う正規分布の、行方向の分散共分散行列に対応する精度行列である
- $w_0$は正則化の際に狙いとするウエイト行列。確率論的には前述の事前分布の平均値行列である。
リッジ回帰の解釈でよくあるのが、wで張るロス平面上での説明図だが、それは式にかかれていることを単に図にしているに過ぎないため、個人的にはリッジ回帰は確率論的には事前分布に
$$p(w)=\mathcal N(w|w_0, \Lambda^{-1};G^{-1})$$
を予め想定した際の回帰である、という解釈の方がしっくりくる。この辺は次の記事らへんで、具体例交えてまとめようと思う。
線形回帰とは？
線形回帰とはそれ以上でもそれ以下でもない以下の式だと思っている。
$$\mathbf y=\mathbf X^T\mathbf w$$
$$\mathbf w=\min_{\mathbf w}\sum_i (\mathbf (\mathbf y_i-\mathbf X_i^T\mathbf w)^2+||\lambda\mathbf w||_2)$$
ただし、$\mathbf y$,$\mathbf x$は列ベクトルとし、$$||w||_2$$は全成分の二乗和を意味する。
カーネル法であろうが重回帰だろうが単回帰だろうが、これ1本で十分。
違いは$X$の取り方だけである。
Xは入力データそのままである必要がない
$X$は入力データそのものである必要はなく、入力データからの写像であってもいいわけだ。
たとえば、入力$(x,y)$に対して$X=(x,1)$
というベクトルを考えよう。冒頭の式に入力したら
$$y=w_0x+w_1$$
であり、回帰係数$w_0,w_1$の線形回帰のできあがりである。
複数の説明変数がある重回帰は、そもそも入力をそのまま使ってもいいし、定数項を入力に加えて
$$\hat X=X\oplus 1$$
$$y=\hat X^Tw$$
ってやればいい(ここからは太字でないものもベクトルや行列として扱う)。
$$X=(1, x, x^2, x^3,\dots,x^n)$$
ならば多項式回帰である。
当然、個々のデータ点x_iを使って
$$k(x,x_i)=\exp(-(x-x_i)^T(x-x_i))$$
のようなカーネル関数を定義して
$$X=(k(x,x_0),k(x,x_1),\dots)$$
というベクトルで回帰してもいい。これはカーネル法と呼び、別記事で説明する。
最小二乗法で最適値を得る
準備: 行列のトレースの行列微分
ある行列$A$に依存するスカラー関数$f$について、$A$の各成分で微分したものを並べて行列にする操作を$nabla_A f(A)$として以下の式で定義しよう。
$$\nabla_A f(A)=\begin{pmatrix}
\frac{\partial f}{\partial a_{11}} &
\frac{\partial f}{\partial a_{12}} & \ldots \
\frac{\partial f}{\partial a_{21}} &
\ddots \
\vdots
\end{pmatrix}$$
$f(A)$が行列のトレースを使って
$$f(A)=tr(CABA^TC^T)=\sum_{i,j,k,l,m}C_{ij}A_{jk}B_{kl}A_{ml}C_{im}$$
のようなとき、微分は以下のように計算できる。
$$\nabla_A tr(CABA^TC^T)=\sum_{i,l,m}C_{ij}B_{kl}A_{ml}C_{im}+\sum_{i,j,k}C_{ij}A_{jk}B_{kl}C_{im}=C^TCA(B^T+B)$$
$$\nabla_A tr(CABD^T)=C_{ij}B_{kl}D_{il}=C^TDB^T$$
$$\nabla_A tr(DBA^TC^T)=D_{ij}B_{jk}C_{il}=C^TDB$$
$$\nabla_A tr((CA-D)B(CA-D)^T)=(C^TCA-C^TD)(B^T+B)$$
ここから先は、この微分計算を駆使して計算を進める。
最適値計算
最小二乗誤差は、入出力が多次元である前提で以下で定義できる。
$$\mathcal L=\sum_i (X_i^Tw-Y_i)(-X_i^Tw-Y_i)^T+\sum_{i,j} \lambda w_{ij}^2$$
ここで、$X_i, Y_i$は$i$番目のサンプルに対応している。
サンプルを行方向に並べ、1サンプルで得られる列ベクトルを並べて行列として扱おう。これを計画行列とよぶ。
$$X={x_{ij}}, Y={y_{ij}}$$
計画行列$X,Y$を利用すると、
$$\mathcal L=\sum_{i,j,k,l} (X_{ij}w_{jk}-Y_{ik})(X_{il}w_{lk}-Y_{ik})+\lambda\sum_{i,j}w_{ij}w_{ij}=tr((Xw-Y)(Xw-Y)^T)+\lambda tr(w^Tw)$$
あとは、この式を$w$について微分し、0になるような$w$を求めればよい。
冒頭の準備であらかじめ導出したトレースの微分を使えば、
$$\nabla_w \mathcal L=2(X^TXw-X^TY)+2\lambda w$$
である($A,B,C,D$に$w,I,X,Y$を割り当てれば得られる。$I$は単位行列)。右辺が0であることが$\mathcal L$が最小値のときなので、
$$(X^TX+\lambda I)w-X^TY=0$$
$$\therefore w=(X^TX+\lambda I)^{-1}X^TY$$
となる。Xを冒頭のとおり適する形状として与えてやれば、単回帰、重回帰、カーネル法は同じ式でウエイトを求められる。
最尤法による最適値探索
さて、今度は確率論的考え方で回帰線を求てみよう。以降、xは写像込みの入力だとする。
得られたデータ(x,y)は線形な関係にノイズ$\epsilon$が乗った結果できたとしよう。そのとき、ウエイトが与えられた際に$(x_i,y_i)$が得られる確率は以下である。
$$p(y_i|w)=x_i^Tw+\epsilon$$
である。もし\epsilonが平均0分散vの正規分布に依存したとすると、
$$p(y_i|w)=\mathcal N(y_i|x_i^Tw, vI)$$
と定義できる。
複数のサンプルがすべて得られたときの同時確率は、各データが独立にサンプリングされていると仮定して、以下の積になる。
$$p(Y|X,w)\prod p(y_i|X,w)=\mathcal N(Y|Xw, vI)$$
ただし$Y,X$は計画行列、分散共分散行列$vI$は計画行列の列ベクトルに対するものとする。
右辺を展開すると、
$$p(Y|X,w)=\exp(-v^{-1}tr((Y-Xw)(Y-Xw)^T))$$
$$\ln p(Y|X,w)=-v^{-1}tr((Y-Xw)(Y-Xw)^T)$$
となる。既に得られている$(X,Y)$に対して、変数wを条件として
$$L(w)=p(Y|X,w)$$
を計算したものは尤度というのだった。
この尤度を最大にする$w$が、データ(X,Y)を最もうまく表現する直線だと言えるのではないか？
これを尤度最大化問題と呼ぶ。解いてみよう。
尤度は対数をとっても最大点に変化はないので、対数尤度を最大化する点を探そう。
$$\max_w L(w)=-\min_w \ln L(w)$$
$$=\min_w \frac{tr((Y-Xw)(Y-Xw)^T)}{v}$$
$$=\min_w \frac{tr((Xw-Y)(Xw-Y)^T)}{v}$$
ということで、最小二乗法と係数を除いて同じ最適化式になる。
ようは、最小二乗法で求めた$w$は確率論的考え方で最適値を考えた場合の$w$と同じ式になる、ということが分かる。
なお定数は最適化結果に影響しないため、データに加わるノイズの大きさ$v$は最適化式に依存しないことになる。
一律でないノイズが加わった場合
ここまでの議論ではサンプルに加わるノイズは常に一定である仮定をしたが、サンプリングに対してノイズは一定でない場合が多いだろう。
おまけに多変量な場合は成分ごとに独立していないノイズもあるだろうということで、サンプル同士の独立性は崩さない範囲でノイズの自由度を最大にしてやろう。
このとき、ノイズはサンプルごとに異なる分散共分散行列$V_i$をもつことになる。便宜上分散共分散テンソル(仮称)をVと定義しよう。
すると、尤度最大化の式は以下のように変更される。
$$-\ln L(w)=tr((Xw-Y)V^{-1}(Xw-Y)^T)$$
なお、$V^{-1}$はサンプリング成分ごとに逆行列をとるものとし、行列積もサンプリング成分ごとに行うものとする。
$$(Xw-Y)V^{-1}(Xw-Y)^T$$
はサンプリング成分が三個独立で残るので三階のテンソルであり、trによって対角成分(つまり三個の成分)をそろえた和となる。
これは成分ごとに分解すると、以下の計算と等価である。
$$-\ln L(w)=\sum_{i,j,k,l,m}(X_{ij}w_{jk}-Y_{ik})S_{ikl}(w_{ml}X_{im}-Y_{il})$$
ただし、$S=V^{-1}$とした。
このままではキレイな式にはならないので、$S$に対してちょっと条件を加えよう。具体的には、分散共分散行列はサンプリングごとに変わらないとする。
$$S_{ijk}=M_{i}G_{jk}$$
ここで、$G$は2つの成分方向に対する分散共分散行列の逆行列であり、精度行列ともいう。$M$はサンプリングの分散の逆数(=精度)ベクトルとした。
これで先ほどの式を修正してみよう。
$$-\ln L(w)=\sum_{i,j,k,l,m}(X_{ij}w_{jk}-Y_{ik})M_{i}G_{kl}(w_{ml}X_{im}-Y_{il})=tr(M(Xw-Y)G(Xw-Y)^T)$$
ただし、$M$は対角成分に各サンプルの精度、非対角成分を0とする行列として置き直した。
微分をとると、
$$\nabla_w -\ln L(w)=X^TM^{T}XwG^{T}+X^TMXwG-X^TM^TYG^{T}-X^TMYG$$
いま、$M$と$G$はそれぞれ対称行列なので
$$\nabla_w -\ln L(w)=2X^TMXwG-2X^TMYG$$
となり、勾配0で尤度最大になるので、
$$X^TMXwG=X^TMYTG$$
を解けばウエイト$w$は求まる。
左辺の$(X^TMX)$は正方行列であるので、正則であれば(=解が定まるならば)、wについて解くと以下になる。
$$w=(X^TMX)^{-1}X^TMY$$
このことから分かることは以下である。
- 成分ごとの相関や分散の関係を分散共分散行列で与えても、少なくとも正則化項のないサンプリングなら回帰結果に影響しない
- サンプリングごとに分散が異なる場合は、誤差の-2乗(=精度)が重みとしてかかる
ところで、ノイズが正規分布に従うと仮定するとサンプリングの分散は誤差(標準偏差)の二乗で定義できる。
ゆえに、上記式はサンプリングごとに測定誤差幅が異なる場合の重みつき最小二乗法に対応する。
事前分布を仮定した回帰方法
この記事の最後に、今度は事前分布を仮定した回帰について議論する。
あるデータが得られた際のwの事後確率を考え、それを最大化するようにwを決めることを考えよう。
これは、ベイズの定理を使えば
$$p(w|X,Y)=\frac{p(Y|X,w)p(w)}{Z}$$
で表現できる。ただし、$Z$は周辺尤度であり、
$$Z=\int p(Y|X,w')p(w')dw'$$
で定義できる$w$に依存しない定数である。
前節と同様に$-\ln(p(w))$に対して最小値を求めよう。
$$-\ln p(w|X,Y)=-\ln p(Y|X,w)-\ln p(w)+\ln Z$$
これを$w$に対する勾配が0になるように$w$が決まれば、最適値となる。
$$-\nabla_w\ln p(w|X,Y)=-\nabla_w\ln p(Y|X,w)-\nabla_w p(w)$$
$p(Y|X,w)$はwに対する尤度であるが、これはノイズを前述の分布に即すると考えれば、
$$\nabla_w(-\ln p(Y|X,w)=2X^TMXwG-2X^TMYG$$
と計算できるのだった。
右辺第二項$p(w)$は事前分布である。言い換えると、ウエイトwはデータがない状態ではp(w)に従ってサンプリングされる、と考える。
たとえばウエイトの事前分布を全成分、平均0分散$\lambda$の正規分布とする。
$$p(w)=\mathcal N(w|0, \lambda I)$$
$$-\nabla_w \ln p(w)=2\lambda w$$
これを組み合わせれば事後分布が最大となる$w$は、
$$\nabla_w(-\ln p(w|X,Y)=2X^TMXwG-2X^TMYG+2\lambda w$$
となり、これが0になるときの$w$を求ればよい。一方、この式は簡単な式で解くことはできない。$G$が$w$の右にかかったりかかってなかったりするからだ。
ちょっと事前分布を修正して、
$$\nabla_w(-\ln p(w|X,Y)=2(X^TMXwG-2X^TMYG+2\lambda wG)$$
という具合にできれば、解くことができる。こうなるときの事前分布ほどんな意味をもつか？
さかのぼるように事前分布の対数を決めると、
$$\ln p(w)=-\lambda tr(wGw^T)$$
これは、$w$の列ベクトルuを
$$p(u)=\mathcal N(u|0,G^{-1})$$
からサンプリングすることに相当する。もともと、
$$y=Xw+\epsilon$$
$$p(Y_i)=\mathcal N(Y_i|X_iw, M;G^{-1})$$
でサンプリングする想定で成分間の分散共分散行列$G^{-1}$を定義していた。
列ベクトル$u$とはy成分ごとにかかるウエイトであるので、$y$成分の間に$G^{-1}$の相関があるならウエイトの事前分布もy成分間に同じ相関がでるのは自然な仮定だといえる。
ゆえにウエイトの事前分布を自然な相関を持つ正規分布(広がり$\lambda$)に従うとする。
$$p(w)=\mathcal N(w|0,\lambda^{-1};G^{-1})$$
このとき、最適化は
$$\nabla_w(-\ln p(w|X,Y))=2X^TMXwG-2X^TYG+2\lambda wG$$
となり、解くと
$$w=(X^TMX+\lambda I)^{-1}X^TY$$
となる。つまり事前分布を考慮した最適化の式において成分間の分散共分散行列は無関係となる。
この式はリッジ回帰の最適化条件の式と等価である。
分散共分散行列は自然な正則化項の選び方をすれば最適化結果に依存しないことがわかったので、一旦考慮をやめて単位行列にしよう。
すると、
$$-\ln p(w|X,Y)=tr(M(Xw-Y)(Xw-T)^T)+\lambda tr(ww^T)$$
となり、これは重みつきリッジ回帰
$$\mathcal L=tr(M(Xw-Y)(Xw-Y)^T)+tr(w^Tw)$$
と一致する。(ここで$tr(w^Tw)=tr(ww^T)$)
ラムダ項の一般化
y軸側の分散共分散行列は事前分布側も揃えた方がいい、という話はしたが、もう拡張の余地はないのか？という疑問が残る。そういえば$\lambda$が依然スカラー量のままだが、ここも行列に拡張できないのか？
結論からいえば、できる。
いままでは事前分布を、なかば辻褄あわせ的に列成分同士の相関性について着目してきたが、行成分同士の相関性をまだ議論していなかったので、そちらも考慮しよう。
サンプリングが行成分同士の相関を持つ場合、その精度行列を$A$とすると、
$$-\ln p(w)=tr(w^TAw)$$
となるべきである。トレースの中は行列積は循環できるので、
$$-\ln p(w)=tr(Aww^T)$$
とできる。ついでに列成分側の分散と相関も考慮すると
$$-ln p(w)=tr(AwGw^T)$$
とできる。かけ算の位置的に、$A$が$\lambda$の行列版である、といえる。
つまり、目的関数
$$\mathcal L=tr(M(Xw-Y)G(Xw-Y)^T)+tr(\Lambda wGw^T)$$
を最小化する拡張版リッジ回帰は、確率論的立場で意味を考えるとウエイト行列が従う事前分布を
$$p(w)=\mathcal N(w|0, \Lambda^{-1};G^{-1})$$
と、行成分と列成分の分散がそれぞれ分散共分散行列$\Lambda^{-1},G^{-1}$であるような正規分布に従うと仮定した場合に、データが得られた際の事後分布を最大にするような$w$を求めるような操作だと解釈できる。
なお、事前分布に正規分布でなくラプラス分布
$$p(w)=\frac{1}{Z}\exp(-\lambda|w|)$$
に従うと仮定した場合は、ラッソ回帰に対応する。
$$\mathcal L=M(Xw-Y)(Xw-Y)^T+\lambda|w|$$
また、事前分布が平均値行列を0としない正規分布に従う場合は、平均を$w_0$とすれば
$$\mathcal L=tr(M(Xw-Y)G(Xw-Y)^T)+tr(\Lambda(w-w_0)G(w-w_0)^T)$$
$$\nabla_w \mathcal L=X^TMXwG-X^TMYG+\Lambda (w-w_0)G$$
$$w=(X^TMX+\Lambda)^{-1}(X^TMY+\Lambda w_0)$$
とできる。
今回の記事はここまで。次回は具体例のテストなどを予定。
所感
以前から頭の中で考えていたことだが、こうやって書き下すと予想以上に長くなった。あと何個か、頭の中で温めている話があるので、どんどん出していこうと思う。