Hastie, T., (2020) Ridge Regularization: An essential concept in data science. Technometrics.
Hastie先生、リッジ正則化についてこれでもかこれでもかこれでもかこれでもかと語り倒すの巻。都合により勉強したい箇所があって目を通したんだけど、正直、疲れた…
1. リッジ回帰と線形回帰
OLS回帰では $$ \hat{\beta} = (X^\top X)^{-1} X^\top y $$だが、\(n \times p\)行列\(X\)の列ランクが\(p\)を下回るとき(ないし条件数が大きいとき)、0(ないし0に近い)固有値を持つ\(X^\top X\)の逆行列をとることになり、よろしくない。そこで対角成分に「リッジ」を足して、$$ \hat{\beta}_\lambda = (X^\top X + \lambda I)^{-1} X^\top y \ \ \ \ldots {\rm (2)}$$ただし\(\lambda > 0\)といたしましょう。というのが、Hoerl & Kennard(1970)が提案したリッジ回帰であった。
最適化問題として捉えると $$ minimize_\beta || y – X\beta ||^2 + \lambda||\beta||^2 \ \ \ \ldots {\rm (3)}$$ である(\(||\cdot||\)はL2ノルム(ユークリディアンノルム))。あるいは、こういうラグランジュ形式じゃなくて上界を使って $$ minimize_\beta ||y-X\beta|| \ {\rm subject \ to} \ ||\beta|| \leq c$$ と書いても同じことである。
\(\lambda\)の値をどうするか？ あくまで数値的な解法としてみるのであれば、0.001とか、\(X^\top X\)の最大の固有値の逆数とかでよい。
[ここでlasso, elastic net, group lassoとの比較。メモ省略]
ベイズ流の見方もできる。$$ (y_i|\beta, X = x_i) \sim x_t^\top \beta + \epsilon_i, \ \ \epsilon_i \ {\rm iid} \ N(0, \sigma^2_\beta I)$$と仮定しよう[原文にはないけど左辺をかっこで包んだ。だってみにくいんだもん]。さらに\(\beta\)も確率変数で \(\beta \sim N(0, \sigma^2_\beta I) \)と仮定しよう。すると負の対数事後分布は \( \lambda = \sigma^2_\epsilon/\sigma^2_\beta \)としたときの式(3)に比例し、事後平均は式(2)となる。事前分散\(\sigma^2_\beta\)が小さいほど、事後平均は( \(\beta\) の事前平均である)0にむかって縮小する。
[速い！速いよ先生!! ゆっくりやりなおすとこういうことであろう。
\(\beta\)の事後分布は $$ p(\beta|y, X) \propto p(\beta) p(y | X, \beta)$$ 正規分布と正規分布の積だから正規分布である。ここで $$ p(\beta) \propto \exp \left(-\frac{||\beta||}{2\sigma^2_\beta} \right)$$ $$ p(y | X, \beta) \propto \exp \left( -\frac{||y-X\beta||}{2\sigma^2_\epsilon} \right)$$ だから $$ p(\beta|y, X) \propto \exp \left( -\frac{||y-X\beta||}{2\sigma^2_\epsilon} – \frac{||\beta||}{2\sigma^2_\beta} \right)$$ この分布の最頻値(正規分布だから平均値でもある)は、この密度関数を最大化する\(\beta\)である。ということは、対数を\(-2 \sigma^2_\epsilon\)倍した \( ||y-X\beta|| + \frac{\sigma^2_\epsilon}{\sigma^2_\beta} ||\beta|| \) を最小化する\(\beta\)である。\( \lambda = \frac{\sigma^2_\epsilon}{\sigma^2_\beta} \)とすれば式(3)と一致する。というわけで、\(\beta\)のリッジ推定量は、ベイズ的にいえば、\(\sigma^2_\epsilon, \sigma^2_\beta\)が既知のときの\(\beta\)の事後最頻値(事後平均)である]
[ここでlassoと比べるコメント。メモ略]
2. リッジ回帰とSVD
たいていの応用場面では、\(\lambda\)をチューニングパラメータと捉えて手元の問題に照らしてよい値を得なければならない。validation用のデータをとっておくとか、LOO-CVとか、\(C_p\)とか、アプローチはいろいろあるけれど、式(3)で二乗誤差損失を考えているからには、とにかくSVDを使って考えるのが効率的である。
\(X\)のSVDを\(X = UDV^\top\)とする。これはfull formであって、\(X\)が\(n \times p\)なら\(U\)は\(n \times n\), \(V\)は\(p \times p\), \(D\)は\(n \times p\)である[つまり\(X\)がランク落ちしていても\(D\)の対角にゼロを埋めてサイズを維持するってことね]。\(D\)の対角要素を大きいほうから\(d_1, \ldots, d_m \geq 0\)とする。
式(2)に代入すると, $$ \hat{\beta}_\lambda = \sum_{d_j > 0} v_j \frac{d_j}{d^2_j + \lambda} \langle u_j, y \rangle$$ である。
[Hastie先生、俺に分からせようという気持ちが全くないね… 仕方がないので自分でゆっくりやりなおすと、$$ \hat{\beta}_\lambda = (X^\top X + \lambda I)^{-1} X^\top y $$ $$ = (V D^\top U^\top U D V^\top + \lambda I)^{-1} V D^\top U^\top y $$ $$ = (V D^2 V^\top + \lambda V V^\top)^{-1} V D^\top U^\top y $$ $$ = (V (D^2 + \lambda I) V^\top)^{-1} V D^\top U^\top y $$ \( (ABC)^{-1} = C^{-1} (AB)^{-1} = C^{-1} B^{-1} A^{-1}\)だから$$ = (V^\top)^{-1} (D^2 + \lambda I)^{-1} V^{-1} V D^\top U^\top y$$ \(V\)は直交行列だから\(V^\top = V^{-1}\)なので $$ = V (D^2 + \lambda I)^{-1} D^\top U^\top y $$ \( D^2 + \lambda I \)は要素\(d^2_j+\lambda\)を持つ対角行列。その逆行列は、要素\(\frac{1}{d^2_j+\lambda}\)を持つ対角行列。これを\(D\)に掛けると、要素\(\frac{d}{d^2_j+\lambda}\)を持つ対角行列。だから$$ \hat{\beta}_\lambda = \sum v_j \frac{d_j}{d_j^2 + \lambda} \langle u_j, y \rangle $$ となる。やれやれ、追いついた…]
この\(\hat{\beta}_\lambda\)をつかってあてはめた値は $$ \hat{y}_\lambda = \sum u_j \frac{d^2_j}{d^2_j + \lambda} \langle u_j, y \rangle$$ となる[ \( \hat{y}_\lambda = X \hat{\beta}_\lambda = U D (D^2 + \lambda I)^{-1} D^\top U^\top y \) だから]。いま\(X\)が列中心化されているとしよう(ふつうそうする。切片に罰則を与えたくないから)。このとき\(u_j\)は\(X\)の主成分である。\(\langle u_j, y\rangle\)は\(y\)の第\(j\)主成分上の座標である。\(\lambda\)を無視すれば、上の式は主成分を基底にとった\(X\)の列空間に\(y\)を写像していることになる。\(\lambda\)は\(\hat{y}_\lambda\)の座標を縮小させている。その程度は\(j\)が進むと共に大きくなる。
[ここで \(p > n\)のときにどうなるかという話。メモ省略]
3. リッジとバイアス-バリアンス・トレードオフ
リッジ回帰の係数は原点に向かって縮小する。\(n > p\)で\( X \)がフルランクだったら、$$ Bias(\hat{\beta}_\lambda) = \sum_j^p v_i \frac{\lambda}{d^2_j + \lambda} \langle v_j, \beta \rangle$$である。いっぽう分散は、\(\sigma^2 = Var(\epsilon)\)として$$ Var(\hat{\beta}_\lambda) = \sigma^2 \sum_j^p \frac{d^2_j}{(d^2_j+\lambda)^2} v_j v_j^\top$$ である。新しい \(x_0\)についての予測のMSEは $$ MSE(x_0, \lambda) = x_0^\top Var(\hat{\beta}_\lambda) x_0 + (x_0^\top Bias(\hat{\beta}_\lambda))^2$$ となる。モデルが正しくても縮小推定量のほうが性能がよいということがありうる。
[シミュレーションの紹介。略]
[James-Stein推定量への一般化。略]
4. リッジとdata augmentation
\(X, y\)が中心化されているとして、\(X\)の下に\(p\)行の\(\sqrt{\lambda}I_p\)を付け加え、\(y\)の下に0を付け加える[たぶんp個の0をくっつけるのであろう]。するとOLS係数は\(\hat{\beta}_\lambda\)となる。[このくだり、読んでしばらく呆然とした。まじか…]
平均0, 共分散\(\lambda I\)の任意の多変量分布から\(n_a\)個の値をランダムドローして\(X\)の下にくっつけ、\(n_a\)個の0を\(y\)にくっつける。オリジナルの点は重み1, 追加した点は重み\(1/n_a\)としてWLS推定すると、係数は\(\hat{\beta}_\lambda\)の近似となる。\( \frac{1}{n_a} X^\top_a X_a \approx \lambda I \)だから。[ええええ…]
\(X\)のそれぞれのデータ点を \(x^\prime_{ij} = x_i + \epsilon_{ij}, \ \ \epsilon_{ij} \sim N(0, \frac{\lambda}{n}I) \)とpurturbしても同じことになる。\(X\)の列数を\(m\)として\( \sum_i \sum_j x^\prime_{ij} x^{\prime\top}_{ij} \approx m(X^\top X + \lambda I) \) だから。[もうなにを言われても驚かんわ…]
もっと現実的な話をすると、DNNでイメージを分類するとき、たとえば訓練データのなかにプードルの画像があるとして、それをちょっと回転したり、サイズを変えたり色を変えたりしてもやっぱりプードルに見える。それらを訓練データに足して過学習を避けることがある。これもリッジ正則化みたいなものである。[話がここまで来るともはや完全にアナロジーだけど、でも勉強になりますです]
5. リッジとドロップアウト正則化
[ニューラルネットのドロップアウト正則化とリッジ回帰の関係。知識不足でよくわからんのでメモ省略]
6. リッジとカーネルトリック
[見出しだけでちょっと笑ってしまった。碩学が興に乗って語り倒すという感じだ。すいません、心の準備ができてません… 読まずに飛ばした]
7. リッジとLOO-CV
SVDをつかうと計算の効率がよいという話をしたけれど、同様に、たとえばk-fold CVをやるときもSVDを使えばk回計算しなくて済む。
n-fold CV すなわちLOO-CVについて考えよう。\( (x_i, y_i) \)を除いて求めたリッジ推定値を\( \hat{\beta}^{(-i)}_\lambda \)と書く。なにも除かずに求めた\(n \times n\)のオペレータ行列[その行列を\( y \)の左から掛けると\( \hat{\beta}_\lambda \)になるような行列のことであろう]を $$ R^\lambda = X(X^\top X + \lambda I)^{-1} X^\top $$ と書く。すると $$ LOO_\lambda = \sum_i (y_i – x^\top_i \hat{\beta}^{(-i)}_\lambda )^2 = \sum_i \frac{ (y_i – x_i^\top \hat{\beta}_\lambda)^2 }{ (1-R^\lambda_{ii})^2 } $$である。つまり、リッジ回帰のすべてのLOO残差[の二乗和]は、オリジナルの残差を \(1/(1-R^\lambda_{ii}) \)だけスケールアップすれば得られるわけだ。SVDを使えば $$ R^\lambda = U S(\lambda) U^\top $$ ただし\(S(\lambda)\)は要素\(d^2_j / (d^2_j + \lambda)\)を持つ対角行列である。[なるほど… まあとにかく、OLS回帰と同様に、LOO残差二乗和の計算にn回の反復はいらないってことね]
[ここからなんだか難しくなるのでメモ省略。Sherman-Morrison-Woodbury identityってなに…]
8. リッジと最小ノルムと二重勾配
ディープ・ラーニング・モデルはイメージ分類のようなSN比の高い領域で支配的である。勾配降下法で訓練誤差が0になるまで訓練した多層のCNNが、テストデータで良い成績を上げることが知られている。そういうネットワークはふつう訓練事例数より多くのパラメータを持っている。つまり通常であれば深刻なoverfittingに分類されるべきものである。しかし、パラメータ数を増やすことで成績が良くなってしまう場合がある。
\(p\)が\(n\)よりはるかに大きいとき、損失を二乗誤差として、\(\beta = 0\)からはじめて勾配降下で訓練誤差が0になるまで降りていくと、最小L2ノルム解になることが知られている。なぜそうなるかというと、勾配 \(\Delta_\beta RSS(\beta) = -2 X^\top (y-X\beta) \)が\(X\)の行空間にあり、よってすべての勾配がそこに残っているからだ[…このくだり、学力不足でよく理解できない…涙涙]。
簡単なデモンストレーション。
モデル\(y = f(x) + \epsilon, \ x \in \mathbb{R}^9\)からデータを生成する。\(x \sim N(0, I_9), \epsilon \sim N(0, \sigma^2) \)とし、\(\sigma\)はSN比が\(Var(f(x))/\sigma^2\)となるように決める。\(f\)は非線形・非加法の関数[としか書いてない…どういうこと?]。標本サイズは\(n = 100\)とし、out-of-sample誤差を測る際にはすごく大きなテスト集合を使う。
で、推定の際にはこういう加法モデルを使う。$$ g(x; d) = \theta_0 + \sum_j^9 \theta^\top_j h_j(x_j; d)$$ ここで\(h_j(\cdot; d)\)は自然スプライン基底関数の長さ\(d\)のベクトル。knotは\(x_j\)の訓練値のuniform quantilesを選ぶ[なんだかよくわかんないけど、まあいいや…]。モデルのパラメータ数は\(9d + 1\)。\(d\)は1から30まで動かす。\(d=11\)のとき、パラメータ数が\(n=100\)に一致する。
横軸にパラメータ数、縦軸に{訓練時, 最小ノルムLS, 最適\(\lambda\)}の予測誤差をとったチャートを見よ。訓練時予測誤差はパラメータ数と共に下がり、100を超えると0になる。最小ノルムLS(100より手前ではOLS)の予測誤差は100に近づくと劇的に増大するが、100を超えるとまた下がる[←なるほど]。リッジ回帰なら(\(\lambda\)が最適なら)予測誤差はほぼフラットで、パラメータ数20くらいのときに最小となる。
最小ノルムLSの曲線が、通常のバイアス-バリアンス・トレードオフを示していない点に注意せよ。なぜパラメータ数100近辺で劇的に増大するかというと、訓練誤差0に到達するために、モデルが訓練データを補完しなければならないからだ。その結果、予測曲線は訓練点のあいだですごくwiggleせざるをえなくなり(wiggleってガクガクと波打つっていう感じかな)、L2ノルムも\(\hat{\theta}\)も大きくなる。これはout-of-sample予測に対しては悪い前兆である。しかしパラメータ数が100を超えると、訓練誤差ゼロにもっとスムーズに到達できるようになり、L2ノルムも下がる。よってout-of-sample予測はかえって改善する。
このように、モデルの複雑さというのはパラメータ数だけでは決まらない。それぞれのパラメータにおける解のL2ノルムが小さくなるポテンシャルも役割を果たす。それぞれのパラメータにおける「objective」は、損失と複雑さを結合したものであり、そのせいで、通常のバイアス-バリアンス・トレードオフが覆われてしまう。
[このくだり、細かいところはよく理解できていないのだが、勉強になりました…ような気がします…]
9. リッジ＋リッジ＝lasso
[力尽きたので読み飛ばした…]
10. 考察
いやあ、リッジ無しのデータサイエンスティスト生活なんて想像もできませんね。云々。
… 哀しいかな、後半はちんぷんかんぷんな箇所が多かったんだけど、まあ一番勉強したかったのは2節、リッジ推定量と計画行列のSVDとの関係だったので、これでよしとしよう。よしとするぞ。
4節のデータ拡張のくだりでびっくりした。正則化っていうのは単にデータを追加しただけだという考え方ができるのか… ここのくだりが一番の収穫であった。