読了:Hastie (2020) リッジ正則化についてこれでもかこれでもかと語り倒す

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であって、Xn \times pならUn \times n, Vp \times p, Dn \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_jXの主成分である。\langle u_j, y\rangleyの第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に分類されるべきものである。しかし、パラメータ数を増やすことで成績が良くなってしまう場合がある。
 pnよりはるかに大きいとき、損失を二乗誤差として、\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 + 1dは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節のデータ拡張のくだりでびっくりした。正則化っていうのは単にデータを追加しただけだという考え方ができるのか… ここのくだりが一番の収穫であった。