線形回帰モデル - yasuhisa's blog

明日発表の分のゼミの資料。PRMLの3.1.2から3.1.5までです。先週のはこの辺に書いている。

今日の日記 - yasuhisa's blog

最小二乗法の幾何学

ここではN=3と固定して考えてみる
- ということなので、3次元空間で考える
- 各軸が $t_1$ 、 $t_2$ 、 $t_3$ で与えられる3次元空間
図についてはここを見る
M < NのケースとしてM=2の場合を考える
- すると次元が一つ落ちるので、ベクトル $\varphi_1$ 、 $\varphi_2$ は(図でいうところの)M=2次元の線形部分空間Sを貼る
- この平面を動き回る
yはベクトル、のパラメータによる線形結合で表わされるので、線形部分空間Sの中のどこにいてもよい
- どこにいてもいいんだけど、どこにいると一番いいんでしょう
  - tとyとの距離が一番近いところがいいよね→tとyとのユークリッド距離
  - 図からも適切なtというのはtの部分空間Sの上への正射影に対応していることが分かる
逆行列を取る関係で、が非正則、非正則に近いとき、数値計算が不安定になる
- この問題は特異値分解(SVD)によって解決することができる

SVDとは何か?

今調べている感じなので、全然知らない
rankがrのの実行列Aはという形に分解できる
- ただし、Uはn次直行行列、Vはm次直行行列、 $\Sigma$ は $\begin{align} ( \array{ D & O_{r,m-r} \\ O_{n-r,r} & O_{n-r,m-r} } ) \end{align}$ 、 $D = \mbox{diag}(\sigma_1,\cdots,\sigma_r),(\sigma_1 \geq \sigma_2 \geq \sigma_r >0)$ であるような行列
- 直行行列は $M^T M = M M^T = I$ を満たす行列
- 対角要素 $\sigma_i(i=1,\cdots,r)$ はAの特異値と呼ばれる
RのSVDによる実行例
- dが $\Sigma$ のdiagなようだ

> example(svd)

svd> hilbert <- function(n) { i <- 1:n; 1 / outer(i - 1, i, "+") }

svd> X <- hilbert(9)[,1:6]

svd> (s <- svd(X))
$d
[1] 1.668433e+00 2.773727e-01 2.223722e-02 1.084693e-03 3.243788e-05
[6] 5.234864e-07

$u
            [,1]       [,2]        [,3]        [,4]        [,5]         [,6]
 [1,] -0.7244999  0.6265620  0.27350003 -0.08526902  0.02074121 -0.004024550
 [2,] -0.4281556 -0.1298781 -0.64293597  0.55047428 -0.27253421  0.092815916
 [3,] -0.3121985 -0.2803679 -0.33633240 -0.31418014  0.61632113 -0.440903754
 [4,] -0.2478932 -0.3141885 -0.06931246 -0.44667149  0.02945426  0.530119859
 [5,] -0.2063780 -0.3140734  0.10786005 -0.30241655 -0.35566839  0.237038375
 [6,] -0.1771408 -0.3026808  0.22105904 -0.09041508 -0.38878613 -0.260449267
 [7,] -0.1553452 -0.2877310  0.29280775  0.11551327 -0.19285565 -0.420944825
 [8,] -0.1384280 -0.2721599  0.33783778  0.29312535  0.11633231 -0.160790254
 [9,] -0.1248940 -0.2571250  0.36542543  0.43884649  0.46496714  0.434599540

$v
           [,1]       [,2]       [,3]        [,4]        [,5]         [,6]
[1,] -0.7364928  0.6225002  0.2550021 -0.06976287  0.01328234 -0.001588146
[2,] -0.4432826 -0.1818705 -0.6866860  0.50860089 -0.19626669  0.041116974
[3,] -0.3274789 -0.3508553 -0.2611139 -0.50473697  0.61605641 -0.259215626
[4,] -0.2626469 -0.3921783  0.1043599 -0.43747940 -0.40833605  0.638901622
[5,] -0.2204199 -0.3945644  0.3509658  0.01612426 -0.46427916 -0.675826789
[6,] -0.1904420 -0.3831871  0.5110654  0.53856351  0.44663632  0.257248908


svd> D <- diag(s$d)

svd> s$u %*% D %*% t(s$v) #  X = U D V'
           [,1]      [,2]      [,3]       [,4]       [,5]       [,6]
 [1,] 1.0000000 0.5000000 0.3333333 0.25000000 0.20000000 0.16666667
 [2,] 0.5000000 0.3333333 0.2500000 0.20000000 0.16666667 0.14285714
 [3,] 0.3333333 0.2500000 0.2000000 0.16666667 0.14285714 0.12500000
 [4,] 0.2500000 0.2000000 0.1666667 0.14285714 0.12500000 0.11111111
 [5,] 0.2000000 0.1666667 0.1428571 0.12500000 0.11111111 0.10000000
 [6,] 0.1666667 0.1428571 0.1250000 0.11111111 0.10000000 0.09090909
 [7,] 0.1428571 0.1250000 0.1111111 0.10000000 0.09090909 0.08333333
 [8,] 0.1250000 0.1111111 0.1000000 0.09090909 0.08333333 0.07692308
 [9,] 0.1111111 0.1000000 0.0909091 0.08333333 0.07692308 0.07142857

svd> t(s$u) %*% X %*% s$v #  D = U' X V
              [,1]          [,2]          [,3]          [,4]          [,5]
[1,]  1.668433e+00  0.000000e+00  4.440892e-16 -2.220446e-16  2.775558e-17
[2,]  4.787837e-16  2.773727e-01  2.775558e-17  8.326673e-17  3.469447e-17
[3,] -4.336809e-19  2.688821e-17  2.223722e-02  8.673617e-18 -3.122502e-17
[4,] -7.261444e-17 -2.734900e-17  1.355253e-18  1.084693e-03 -1.409463e-17
[5,]  4.813392e-17  1.809855e-17  1.922511e-17  1.641550e-18  3.243788e-05
[6,]  4.102411e-17  1.037212e-17 -5.631168e-18 -7.879048e-18 -2.771809e-18
              [,6]
[1,] -1.387779e-16
[2,]  2.081668e-17
[3,] -1.734723e-18
[4,] -1.626303e-18
[5,] -5.141490e-18
[6,]  5.234864e-07

なんでSVDなのか?

Ax=bみたいな問題を考えているとき、 $\mathbf{x} = A^{-1} \mathbf{b}$ みたいにしたいけど、Aは一般に正方行列で、正則であるとは限らない
こういうときに逆行列みたいなの考えたいなあ
- ムーア-ベンローズの一般化逆行列(p140)
これを用いると方程式の解は $\mathbf{x}_0 = V \Sigma^{-1} U^T \mathbf{b} = \sum^r_{i=1} \frac{1}{\sigma_i} (\mathbf{u}^T_i \mathbf{b}) \mathbf{v}_i$ と書ける

逆行列を持たない(丸め誤差の関係で、も含む)時→一次従属な(に近い)関係があるとき
どっかの行がどっかの行とどっかの行の線形和で書き表わせる
- 行列の情報量を削減できるんじゃね?
- 上位の固有値を使うことで低次元で近似しよう
- はの非ゼロな固有値の平方根
  - なんで「非ゼロ」になるの?自明?
  - これか
  - http://books.google.com/books?id=4nmPRu1C4m4C&pg=PA182&lpg=PA182&dq=A%5ET+A+固有値&source=bl&ots=GVaaCSOShv&sig=8f6yJ4BarKcUZSdMHvNEeiutitg&hl=ja&ei=Ayz2See0Norq6gPcqKSGAg&sa=X&oi=book_result&ct=result&resnum=1#PPA182,M1
  - $0 \leq |A \mathbf{v}|^2 = \mathbf{v}^T A^T A \mathbf{v} = \lambda \mathbf{v}^T \mathbf{v} = \lambda$

逐次学習

いわゆるオンライン学習の話
めちゃくちゃデータがでかいときに一度に処理する(バッチ手法)
ではなく、パラメータを逐次更新していく方法
オンライン学習にも色々方法がある
- 前に書いたし、資料も上がっている
- ここでは確率的勾配降下法(stochastic gradient descent)について書いてある
によってパラメータを更新していく
- 二乗和誤差の場合は、この式は $\mathbf{w}^{(\tau + 1)} = \mathbf{w}^{\tau} - \eta ({\mathbf{w}^{(\tau)}}^T \phi_n)\phi_n$ で計算できる

正則化最小二乗法

過学習を防ぐための正則化項
- 過学習については1.1節が詳しい
- 多項式フィッティングでの例→P6とP10
$E_D(\mathbf{w}) + \lambda E_w(\mathbf{w})$
$E_W(\mathbf{w}) = \frac{1}{2} \mathbf{w}^T \mathbf{w}$
$E_D(\mathbf{w}) = \frac{1}{2} \sum^N_{n=1}\{t_n - \mathbf{w}^T \phi(\mathbf{x}_n)\}^2$
$E_D(\mathbf{w}) + \lambda E_w(\mathbf{w}) = \frac{1}{2} \sum^N_{n=1}\{t_n - \mathbf{w}^T \phi(\mathbf{x}_n)\}^2 + \frac{\lambda}{2} \mathbf{w}^T \mathbf{w}$
荷重減衰、パラメータ縮小推定
- 多くのパラメータが0に近づいていく
- P10の表

$\mathbf{w}$ に関する勾配を求めると
$\Phi^T (\mathbf{t} - \Phi \mathbf{w}) + \lambda \mathbf{w} = \Phi^T \mathbf{t} - \Phi^T \Phi \mathbf{w} + \lambda \mathbf{w}$
となり、これを0とすると
$\Phi^T \mathbf{t} = (\Phi^T \Phi + \lambda \mathbf{I}) \mathbf{w}$ となり、 $\mathbf{w} = (\Phi^T \Phi + \lambda \mathbf{I})^{-1} \Phi^T \mathbf{t}$ が得られる。

今は二乗誤差関数について考えたが、一般化して考えることもできる
- $\frac{1}{2} \sum^N_{n=1}\{t_n - \mathbf{w}^T \phi(\mathbf{x}_n)\}^2 + \frac{\lambda}{2} \sum^N_{n=1} |w_j|^q$
実用上はq=1とq=2のケースが多く取り扱われているようで、それぞれL1正則化、L2正則化と呼ばれている

L1正則化とL2正則化のメリットデメリット

L1正則化(Lasso正則化とも呼ばれる)の特徴として

多くの重みが0になる
- 結果が解釈しやすい
特徴数が少なくてL2と精度が同じ
- 推定が高速
微分不可能
- 勾配を使うような最適化手法が直接は使えない

などがある。

ToDo:なんで多くの重みが0になるか

L2正則化の特徴は

全ての重みが0にならない
- Representor Theorem
ちょっとでも性能を上げたい、という時に使われる

凸関数について

を満たすようなもの
- ただし、 $\theta$ は $0 \leq \theta \leq 1$ を満たすような実数
凸関数の和は凸関数
http://www-tkm.ics.nitech.ac.jp/~wadayama/2009-ml/lec02.pdf

出力変数が多次元の場合

tが多次元になった場合を考えよう
- tの要素それぞれに対して別々の基底関数を使ってK個の回帰問題、としてもいいんだけどここではそうでない方法

$\begin{align} \mathbf{y}(\mathbf{x},\mathbf{w}) &= \mathbf{W}^T \phi(\mathbf{x}) \\ &= $ \array{ w_{11} & \cdots & w_{1K} \\ \vdots & \ddots & \vdots \\ w_{M1} & \cdots & w_{MK} } $^T $ \array{\phi_0(\mathbf{x}) \\ \vdots \\ \phi_{M-1}(\mathbf{x})} $ \\ &= $ \array{ w_{11} & \cdots & w_{M1} \\ \vdots & \ddots & \vdots \\ w_{1K} & \cdots & w_{MK} } $ $ \array{\phi_0(\mathbf{x}) \\ \vdots \\ \phi_{M-1}(\mathbf{x})} $ \end{align}$
ここで、1変数でやったときと同じようなことを考える。つまり目標変数が条件付きガウス分布であると仮定する(ただし、ここでは多変数なので多変数のガウス分布であるが)。P145では簡単化のために「等方性」のガウス分布を仮定してある。つまり、 $p(\mathbf{t}|\mathbf{x},\mathbf{W},\beta) = N(\mathbf{t}|\mathbf{W}^T \phi(\mathbf{x}),\beta^{-1} \mathbf{I})$ である。

このように設定した上で、1変数のときと同じように、対数尤度関数について考える。対数尤度関数は次のようになる。

$\begin{align} \log p(\mathbf{t}|\mathbf{x},\mathbf{W},\beta) &= \sum^N_{n=1} \log N(\mathbf{t}|\mathbf{W}^T \phi(\mathbf{x}),\beta^{-1} \mathbf{I}) \\ &= \sum^N_{n=1} \log \frac{1}{(2 \pi)^{\frac{K}{2}}} \frac{1}{|\beta^{-1} \mathbf{I}|^{\frac{1}{2}}} \exp{\{- \frac{1}{2} (\mathbf{t} - \mathbf{W}^T \phi(\mathbf{x}))^T (\beta^{-1} \mathbf{I})^{-1} (\mathbf{t} - \mathbf{W}^T \phi(\mathbf{x}))\}} \\ &= \sum^N_{n=1} \{ \log \frac{1}{(2 \pi)^{\frac{K}{2}}} + \log \frac{1}{|\beta^{-1} \mathbf{I}|^{\frac{1}{2}}} - \frac{1}{2} (\mathbf{t} - \mathbf{W}^T \phi(\mathbf{x}))^T (\beta^{-1} \mathbf{I})^{-1} (\mathbf{t} - \mathbf{W}^T \phi(\mathbf{x})) \} \\ &= \sum^N_{n=1} \log \frac{1}{(2 \pi)^{\frac{K}{2}}} + \sum^N_{n=1} \log \frac{1}{|\beta^{-1} \mathbf{I}|^{\frac{1}{2}}} - \frac{1}{2} \sum^N_{n=1} (\mathbf{t} - \mathbf{W}^T \phi(\mathbf{x}))^T (\beta^{-1} \mathbf{I})^{-1} (\mathbf{t} - \mathbf{W}^T \phi(\mathbf{x})) \\ &= \frac{NK}{2} \log (\frac{1}{2 \pi} ) + \frac{N}{2} \log (\beta) - \frac{\beta}{2} \sum^N_{n=1} (\mathbf{t} - \mathbf{W}^T \phi(\mathbf{x}))^T (\mathbf{t} - \mathbf{W}^T \phi(\mathbf{x})) \\ &= \frac{NK}{2} \log (\frac{1}{2 \pi} ) + \frac{N}{2} \log (\beta) - \frac{\beta}{2} \sum^N_{n=1} || t_n - \mathbf{W}^T \phi(\mathbf{x_n}) ||^2 \end{align}$
…と思ったら本の記述と合わない＞＜PRMLのP145(3.33)式には
$\frac{NK}{2} \log (\frac{\beta}{2 \pi}) - \frac{\beta}{2} \sum^N_{n=1} || t_n - \mathbf{W}^T \phi(\mathbf{x_n}) ||^2$
と書いてありました。多変量のガウス分布の最尤推定が書いてあるP91の(2.118)式と見くらべても自分のやつがあってると思うんだよなー。二項目からKが出てくる気がしない。というわけで本のほうが間違っているんじゃないかと思うんだけど、どうでしょう＞＜。あと「しかしながら、より興味深く、かつ実際によく使われるアプローチは」のところの理由がよく分からなかった。