基底関数を使った線形回帰とカーネル回帰

https://zenn.dev/inaturam/articles/34d7d3ce828e76

上の記事ではデータ行列の特徴次元に対してパラメータに貢献度を表現させる線形モデルを考えたが、基底関数を使いデータのレコードを1つの値に飛ばすという考えから線形モデルを見ていく。
線形回帰における違いとしては、前者は $[xn∈RD]=X∈RN×D,w∈RD[\bold{x}_n \in \mathbb{R}^{D}] = \bold{X} \in \mathbb{R}^{N×D}, \bold{w} \in \mathbb{R}^{D}$ によりデータのd次元目の貢献度を推定していたのに対し、後者は $[φm(xn)∈R]=Φ∈RN×M,w∈RM[φ_m(\bold{x}_n) \in \mathbb{R}] = \bold{Φ} \in \mathbb{R}^{N×M}, \bold{w} \in \mathbb{R}^{M}$ に注目し、m番目の基底関数によって飛ばした後の空間の貢献度を推定する。

e.g.
データをそのまま扱う場合の線形回帰: $y(x,w)=∑d=1Dwdxdy(\bold{x},\bold{w}) = \sum_{d=1}^D w_d x_d$
↑データの特徴次元ごとの重み(パラメータ)が推定される。パラメータだけでなく入力変数$xx$も線形なので、回帰関数は直線になる。
基底関数を使う場合の線形下記: $y(x,w)=∑m=1Mwmφm(x)y(\bold{x},\bold{w}) = \sum_{m=1}^M w_m φ_m(\bold{x})$
↑データの特徴次元は関係なくなり、$φφ$と$ww$の線形結合によって関数が表される。例えば基底関数に多項式 $φm(x)=xmφ_m(x) = x^m$ などを選べば、多項式曲線を使った曲線が得られる。

これを応用すると、基底関数そのもので回帰を行うカーネル回帰を考えることができる。

線形回帰

まず、扱うデータとして特徴次元 $DD$ を持つ入力変数 $xn\bold{x}_n$ と、目的変数 $tnt_n$ を設定する。

$$X=[xn∈RD]∈RN×D=[x1,1…xN,1⋮⋱⋮xD,1…xN,D]t=[tn]∈RD \begin{align*} \bold{X} & = [ \bold{x}_n \in \mathbb{R}^{D} ] \in \mathbb{R}^{N×D} \\ & = \begin{bmatrix} x_{1,1} & \dots & x_{N,1} \\ \vdots & \ddots & \vdots \\ x_{D,1} & \dots & x_{N,D} \\ \end{bmatrix} \\ \bold{t} & = [ t_n ] \in \mathbb{R}^{D} \end{align*}$$

このデータに対して、次の式で表される、パラメータ $w∈RM\bold{w} \in \mathbb{R}^M$ を持つ 線形モデル をfitさせることを考える。

$$y(x,w)=f(∑m=1Mwmφm(x))=w⊤φ(x) \begin{align*} y(\bold{x},\bold{w}) & = f( \sum_{m=1}^{M} w_m φ_m(\bold{x}) ) \\ & = \bold{w}^\top φ(\bold{x}) \end{align*}$$

簡単に記述するため、今は活性化関数 $f:R↦Rf: \mathbb{R} \mapsto \mathbb{R}$ を恒等写像として無視した。また、バイアスパラメータ $w1w_1$ のために $m=1m=1$ では特別に基底関数 $φm:RD↦Rφ_m: \mathbb{R}^D \mapsto \mathbb{R}$ を $φ1(x)=1φ_1(\bold{x}) = 1$ としてダミー運用し、 $φ:RD↦RDφ: \mathbb{R}^D \mapsto \mathbb{R}^D$ を使った内積として書いている。

このパラメータ $w\bold{w}$ をデータから学習してデータに $y(x,w)y(\bold{x},\bold{w})$ を近似させることが 線形回帰 である。

単純な線形回帰では次の方法でモデルを最適化し、パラメータを求める。

1. データが観測できた時に、モデルがどれだけ正しくデータを生成できるかを示す確率(厳密には確率ではない)、 尤度(likelihood) $p(t∣X,w)p(\bold{t}|\bold{X}, \bold{w})$ を定義する。
2. 尤度を最大化するパラメータを求める(最尤推定)。最大値の条件である $∂∂wp(t∣X,w)=0\frac{∂}{∂\bold{w}} p(\bold{t}|\bold{X}, \bold{w}) = 0$ より、最適なパラメータ $wML\bold{w}_{\text{ML}}$ が求まる。

最初に尤度を定める。今回、目的変数 $tt$ が回帰関数 $yy$ にガウスノイズ $ε∼N(ε∣0,β−1)ε \sim \mathcal{N}(ε|0,β^{-1})$ を付加したものと考える。そうすると、 $tt$ は平均 $yy$ 、分散 $β−1β^{-1}$ の正規分布に従うということになる。

$$t=y(x,w)+ε∼p(t∣x,w,β)=N(t∣y(x,w),β−1) \begin{align*} t & = y(\bold{x},\bold{w}) + ε \\ & \sim p(t| \bold{x}, \bold{w}, β) \\ & = \mathcal{N}(t| y(\bold{x},\bold{w}), β^{-1}) \end{align*}$$

データが独立同分布から生成されたと仮定すると、尤度は以下のように書ける。

$$p(t∣X,w,β)=∏n=1NN(tn∣y(xn,w),β−1)=∏n=1N12πβ−1e(tn−w⊤φ(xn))22β−1ln⁡p(t∣X,w,β)=N2ln⁡β−N2ln⁡2π−β2∑n=1N(tn−w⊤φ(xn))2⏟:=βE(w) \begin{align*} p(\bold{t}| \bold{X}, \bold{w}, β) & = \prod_{n=1}^N \mathcal{N}(t_n| y(\bold{x}_n,\bold{w}), β^{-1}) \\ & = \prod_{n=1}^N \frac{1}{\sqrt{2π β^{-1}}} e^{\frac{(t_n - \bold{w}^\top φ(\bold{x}_n))^2}{2 β^{-1}}} \\ \ln p(\bold{t}| \bold{X}, \bold{w}, β) & = \frac{N}{2}\ln β - \frac{N}{2}\ln 2π - \underbrace{\frac{β}{2} \sum_{n=1}^N ( t_n - \bold{w}^\top φ(\bold{x}_n) )^2}_{:= β E(\bold{w})} \end{align*}$$

$EE$ は二乗和誤差で、最小二乗法でこれを最小化することを考えるのは、尤度最大化から来ている。
尤度の最大値は、尤度(=対数尤度)の勾配が0になる点であるから、以下の条件(正規方程式)を解くことで最尤解である回帰パラメータ $wML\bold{w}_{\text{ML}}$ が求まる。

$$0=∂∂wln⁡p(t∣X,w,β)∈RM⇒0=∂∂w12∑n=1N(tn−w⊤φ(xn))2⇒0=∑n=1N(tn−w⊤φ(xn))φ(xn)⊤⇒0=∑n=1Ntnφ(xn)⊤−∑n=1Nw⊤φ(xn)φ(xn)⊤⇒wML=(Φ⊤Φ)−1Φ⊤t \begin{align*} 0 & = \frac{∂}{∂\bold{w}} \ln p(\bold{t}| \bold{X}, \bold{w}, β) \in \mathbb{R}^M \\ ⇒ 0 & = \frac{∂}{∂\bold{w}} \frac{1}{2} \sum_{n=1}^N ( t_n - \bold{w}^\top φ(\bold{x}_n) )^2 \\ ⇒ 0 & = \sum_{n=1}^N ( t_n - \bold{w}^\top φ(\bold{x}_n) ) φ(\bold{x}_n)^\top \\ ⇒ 0 & = \sum_{n=1}^N t_n φ(\bold{x}_n)^\top - \sum_{n=1}^N \bold{w}^\top φ(\bold{x}_n) φ(\bold{x}_n)^\top \\ ⇒ \bold{w}_{\text{ML}} & = ( \bold{Φ}^\top\bold{Φ} )^{-1} \bold{Φ}^\top \bold{t} \end{align*}$$

計画行列 $Φ\bold{Φ}$ は次のようになる。

$$Φ=[φ1(x1)…φM(x1)⋮⋱⋮φ1(xN)…φM(xN)]∈RN×M \bold{Φ} = \begin{bmatrix} φ_1 (\bold{x}_1) & \dots & φ_M (\bold{x}_1) \\ \vdots & \ddots & \vdots \\ φ_1 (\bold{x}_N) & \dots & φ_M (\bold{x}_N) \\ \end{bmatrix} \in \mathbb{R}^{N×M}$$

以上が単純な線形回帰の最尤推定である。この方法は問題があり、データの特徴次元列に相関がある場合など、擬似逆行列 $(Φ⊤Φ)−1Φ⊤( \bold{Φ}^\top\bold{Φ} )^{-1} \bold{Φ}^\top$ の逆行列部分 $Φ⊤Φ\bold{Φ}^\top\bold{Φ}$ が非正則に近い場合にパラメータが極端に大きな値になることがある。これは最尤推定の 過学習 の根源であり、逆行列を正則化することで過学習を防止する必要がある。

このために、正則化を含めてモデルを最適化し、パラメータを求める。

1. 尤度に加えて、パラメータが従うであろう分布として 事前分布 を設定し、ベイズの定理により 事後分布 を求める。
2. 事後分布最大化(MAP推定)により尤度と事前分布の積の最大値を求める。

先述の通り、問題は $w\bold{w}$ が大きな値を取ることであるため、「 $w\bold{w}$ は平均0付近にあってほしい」という信念をを事前分布に与える。

$$w∼p(w∣α)=N(w∣0,α−1I)=1(2πα−1)M2ew⊤w2α−1 \begin{align*} \bold{w} & \sim p(\bold{w}|α) \\ & = \mathcal{N}(\bold{w}| \bold{0}, α^{-1} \bold{I}) \\ & = \frac{1}{(2πα^{-1})^{\frac{M}{2}}} e^{\frac{\bold{w}^\top\bold{w}}{2α^{-1}}} \end{align*}$$

この時、事後分布はベイズの定理より次のように尤度と事前分布の積に比例する。

$$p(w∣X,t,α,β)=p(t∣X,w,β)p(w∣α)p(X,t,w,α,β)⏟const∝p(t∣X,w,β)p(w∣α) \begin{align*} p(\bold{w}|\bold{X},\bold{t},α,β) & = \frac{p(\bold{t}|\bold{X},\bold{w},β)p(\bold{w}|α)}{\underbrace{p(\bold{X},\bold{t},\bold{w},α,β)}_{\text{const}}} \\ & ∝ p(\bold{t}|\bold{X},\bold{w},β)p(\bold{w}|α) \end{align*}$$

最尤推定と同様に、尤度の最大値の条件から正規方程式を導くと、 正則化項 $λ2∣∣w∣∣2\frac{λ}{2}||\bold{w}||^2$ が現れることがわかる。ただし両辺 $ββ$ で割って $λ=αβλ=\frac{α}{β}$ で置いた。

$$0=∂∂wln⁡p(w∣X,t,α,β)⇒0=∂∂w(β2∑n=1N(tn−w⊤φ(xn))2+α2w⊤w)⇒wMAP=(Φ⊤Φ+λI)−1Φ⊤t \begin{align*} 0 & = \frac{∂}{∂\bold{w}} \ln p(\bold{w}|\bold{X},\bold{t},α,β) \\ ⇒ 0 & = \frac{∂}{∂\bold{w}} (\frac{β}{2} \sum_{n=1}^N ( t_n - \bold{w}^\top φ(\bold{x}_n) )^2 + \frac{α}{2} \bold{w}^\top\bold{w}) \\ ⇒ \bold{w}_{\text{MAP}} & = ( \bold{Φ}^\top\bold{Φ} + λ\bold{I} )^{-1} \bold{Φ}^\top \bold{t} \end{align*}$$

事前分布から得られる正則化項が逆行列部分をフルランクに補強するため、行列計算が安定して極端な値が出にくくなる。
このように、正則化二乗和誤差を最小化する線形回帰を Redge回帰 と言う。

カーネル回帰

双対表現: Redge回帰の問題設定において、回帰パラメータ $w\bold{w}$ を直接扱う代わりに、回帰関数そのものを、データと基底関数の計画行列と目的変数を使って表す方法。回帰関数の全てがカーネルで表されるため、全ての問題をカーネル関数を通じて扱えるようになる。

まず、回帰パラメータについて解き、別の書き方を導入する。

$$0=∂∂wln⁡p(w∣X,t,α,β)⇒0=∂∂w(12∑n=1N(tn−w⊤φ(xn))2+λ2w⊤w)⇒w=∑n=1N1λ(tn−w⊤φ(xn))⏟:=anφ(xn)⇒w=Φ⊤a \begin{align*} 0 & = \frac{∂}{∂\bold{w}} \ln p(\bold{w}|\bold{X},\bold{t},α,β) \\ ⇒ 0 & = \frac{∂}{∂\bold{w}} (\frac{1}{2} \sum_{n=1}^N ( t_n - \bold{w}^\top φ(\bold{x}_n) )^2 + \frac{λ}{2} \bold{w}^\top\bold{w}) \\ ⇒ \bold{w} & = \sum_{n=1}^N \underbrace{\frac{1}{λ} ( t_n - \bold{w}^\top φ(\bold{x}_n) )}_{:= a_n} φ(\bold{x}_n) \\ ⇒ \bold{w} & = \bold{Φ}^\top \bold{a} \end{align*}$$

この結果を正則化二乗和誤差 $E(w)E(\bold{w})$ に代入する。

$$E(w)=12∑n=1N(tn−w⊤φ(xn))2+λ2w⊤w)⇒E(a)=12a⊤KKa−a⊤Kt+12t⊤t+λ2a⊤Ka \begin{align*} E(\bold{w}) & = \frac{1}{2} \sum_{n=1}^N ( t_n - \bold{w}^\top φ(\bold{x}_n) )^2 + \frac{λ}{2} \bold{w}^\top\bold{w}) \\ ⇒ E(\bold{a}) & = \frac{1}{2} \bold{a}^\top\bold{KKa} - \bold{a}^\top\bold{Kt} + \frac{1}{2} \bold{t}^\top\bold{t} + \frac{λ}{2} \bold{a}^\top\bold{Ka} \end{align*}$$

ただし $K=ΦΦ⊤\bold{K} = \bold{ΦΦ}^\top$ で、 $Ki,j=k(xi,xj)K_{i,j} = k(\bold{x}_i,\bold{x}_j)$ となっている。
これらの結果から $a=(K+λI)−1t\bold{a} = (\bold{K} + λ\bold{I})^{-1} \bold{t}$ が導かれ、回帰モデルを書き直すと次のようになる。

$$y(x)=w⊤φ(x)=a⊤Φφ(x)=[k(x1,x)⋮k(xN,x)](K+λI)−1t \begin{align*} y(\bold{x}) & = \bold{w}^\top φ(\bold{x}) \\ & = \bold{a}^\top \bold{Φ} φ(\bold{x}) \\ & = \begin{bmatrix} k(\bold{x}_1, \bold{x}) \\ \vdots \\ k(\bold{x}_N, \bold{x}) \\ \end{bmatrix} (\bold{K} + λ\bold{I})^{-1} \bold{t} \end{align*}$$

この結果から、最小二乗法の解は双対性の概念からカーネル関数 $k(x,x′)k(\bold{x}, \bold{x}')$ のみのよって表現できることが示された。
これを用いた回帰として、カーネル関数 $k(x,x′)k(\bold{x}, \bold{x}')$ が、2点の距離のみに依存している(動径基底)ものを使った RBFネットワーク がある。RBFネットワークでは基底関数 $φ(x)φ(\bold{x})$ に「与えられたデータ点からの距離」を表す関数を選ぶことで、回帰関数 $y(x)y(\bold{x})$ をfitする。

まず、データと目的変数の同時分布 $p(x,t)p(\bold{x},t)$ と、データ点を中心に置かれる同時分布の密度関数 $q(x,t)q(\bold{x},t)$ を用いることで、Parzen(カーネル密度)推定を行う。

$$p(x,t)=1N∑n=1Nf(x−xn,t−tn) p(\bold{x},t) = \frac{1}{N} \sum_{n=1}^N f(\bold{x}-\bold{x}_n, t-t_n)$$

ここで回帰関数 $y(x)=E[t∣x]y(\bold{x}) = \mathbb{E}[t|\bold{x}]$ を調べると、以下のように計算できる。

$$E[t∣x]=∫tp(t∣x)dt=∫tp(x,t)dt∫p(x,t)dt=∑n=1N∫tq(x−xn,t−tn)dt∑n′=1N∫q(x−xn′,t−tn′)dt↓q(x)=∫q(x,t)dt=∑n=1Nq(x−xn)∑n′=1Nq(x−xn′)⏟:=k(x,x′)tn \begin{align*} \mathbb{E}[t|\bold{x}] & = \int t p(t|\bold{x}) dt \\ & = \frac{ \int t p(\bold{x}, t) dt }{ \int p(\bold{x}, t) dt } \\ & = \frac{ \sum_{n=1}^{N} \int t q(\bold{x} - \bold{x}_n, t - t_n) dt }{ \sum_{n'=1}^{N} \int q(\bold{x} - \bold{x}_{n'}, t - t_{n'}) dt } \\ & ↓ q(\bold{x}) = \int q(\bold{x}, t) dt \\ & = \sum_{n=1}^{N} \underbrace{\frac{ q(\bold{x} - \bold{x}_n) }{ \sum_{n'=1}^{N} q(\bold{x} - \bold{x}_{n'}) }}_{:=k(\bold{x},\bold{x}')} t_n \\ \end{align*}$$

この時 $k(x,x′)=q(x−xn)/∑n′=1Nq(x−xn′)k(\bold{x},\bold{x}') = q(\bold{x} - \bold{x}_n)/\sum_{n'=1}^{N} q(\bold{x} - \bold{x}_{n'})$ は距離を意味するカーネルになっており、 $x\bold{x}$ と $xn\bold{x}_n$ が近いほど大きな重みが与えられるため、データ点の周りに大きな密度を持つ予測分布が得られる。

---

カーネル回帰とその一般形であるガウス過程について、より詳しくは以下の記事で実装と解説をした。

https://zenn.dev/inaturam/articles/9d6f6428c4ce0a

以上の知識を整理すれば、NNの数理について学ぶ準備ができる気がする。
PRMLは素晴らしい本である。