Stein's Unbiased risk estomator : Steinの補題とリスク推定 からの続き。

前稿のまとめ

\[ y_i=\mu_i+\sigma_nz_i\quad(i=1,\dots,n)\]

によって生成されている観測値から未知の\(\mu_i\)を推定するとき、そのリスク\(R=\mathbb{E}[(\mu-\widehat\mu)^2]/n\)の不偏推定が

\[ \widehat{R}=-\sigma_n^2+\frac1n\sum_{i=1}^n(y_i-\widehat\mu_i)^2+\frac{2\sigma_n^2}{n}\sum_{i=1}^n\frac{\partial\widehat\mu_i(y)}{\partial y_i}\]

によって計算でき、これをSURE : Stein's unbiased risk estimatorといった。

回帰におけるSURE

\(n\)個の観測点\(y_1,\dots,y_n\)が次のように生成されている場合を考える。

\[ y_i=f^*(x_i)+\sigma_nz_i\]

ここに\(x_1,\dots,x_n\in\mathbb{R}^p\)は説明変数や入力変数を示していて、この関係式は回帰における問題設定でよく用いられる形になっている。

\(n\)個の説明変数が既知であり固定されている場合に\(\mu_i=f^*(x_i)\)とすればこれはSUREにおける問題設定と同じ枠組みによって説明できる。何かしらの回帰モデルによって\(f^*\)の推定\(\widehat f\)を求められたとすると、\(\widehat\mu_i=\widehat{f}(x_i)\)によってSUREにおける\(\widehat\mu\)を指定するように\(\widehat\mu\)の部分も修正する。

これは真の平均\(\mu\)に関する未知性を\(f^*\)という関数に押し付けたものと考えてよい。唯一異なるのは観測\(y\)に付随する\(x_1,\dots,x_n\)という既知の情報が増えたことだけである。\(y=(y_1,\dots,y_n)\)としたのと同様に\(X=(x_1^\top,\dots,x_n^\top)^\top\)という\(n\times p\)行列を考える。これを計画行列（design matrix）という。

SUREにおいて\(y\)という観測値によってのみ決定されていた\(\widehat\mu(y)\)という推定は、この場合には\(X\)にもよることが自然であり、\(\widehat f\)という\(f^*\)の推定も同様に\(y,X\)に依存する。

\[ \widehat\mu_i=\widehat\mu_i(y,X)=\widehat{f}(x_i;y,X)\]

線形回帰におけるリスク推定

ここまではパラメトリックに限定しない一般的な回帰の話。ここから特にパラメトリック回帰モデルの中でも単純な線形回帰モデルを例に挙げる。線形回帰モデルでは、

\[ f^*\in\mathcal{F}:=\{f(x)=x^\top\beta\;|\;\beta\in\mathbb{R}^p\}\]

を仮定して最もデータにフィットするような\(\widehat{f}\)を\(\mathcal{F}\)の中から選択する。

\[ \widehat{f}(x)=x^\top\widehat{\beta},\quad\widehat\beta:=\mathop{\rm argmin}_{\beta\in\mathbb{R}^p}[(y-X\beta)^2]\]

というようにデータに対する二乗誤差を最小化するような\(\widehat{\beta}\)を用いる方法がまず考えられる。このような\(\widehat\beta\)は陽に求まることがわかっていて、

\[ \widehat\beta=(X^\top X)^{-1}X^\top y\]

従って

\[ \widehat{f}(x_i;y,X)=x_i^\top(X^\top X)^{-1}X^\top y\quad(i=1,\dots,n)\]

であり\(i=1,\dots,n\)の部分をまとめて表記すると結局

\[ \widehat\mu(y,X)=X(X^\top X)^{-1}X^\top y\]

がわかる。これを用いて\(\widehat{R}\)の第3項を計算してみよう。今\(H=X(X^\top X)^{-1}X^\top\)とおくとこれは\(n\times n\)行列である。\(\widehat\mu\)の\(i\)行目を\(y_i\)で微分してその総和を取るという操作は、\(H\)の対角和（トレース）を計算することに他ならない。従って

\[ \begin{align}\widehat{R}&=-\sigma_n^2+\frac1n\sum_{i=1}^n(y_i-Hy|_i)^2+\frac{2\sigma_n^2}{n}\mathrm{tr}H\\&=-\sigma_n^2+\frac1n((I_n-H)y)^2+\frac{2\sigma_n^2}{n}\mathrm{tr}H\end{align}\]

とかける。実は\(\mathrm{tr}H\)という量は簡単なトレースの性質\(\mathrm{tr}(AB)=\mathrm{tr}(BA)\)を使えば簡単に計算できる。

\[ \mathrm{tr}H=\mathrm{tr}(X(X^\top X)^{-1}X^\top)=\mathrm{tr}((X^\top X)(X^\top X)^{-1})=\mathrm{tr}(I_p)=p\]

従ってこの場合の第3項はモデルにおける説明変数の個数に等しい。モデル選択に影響する\(H\)と\(y\)に関係するところだけに注目すると以下のようになる。

\[ \widehat{R}=\frac1n\sum_{i=1}^n((I_n-H)y)^2+2\sigma_n^2\frac{p}{n}\]

ちなみに\(H\)には名前があってhat matrixといわれる。第1項は1サンプル当たりの訓練誤差を表し第2項は1サンプル数当たり説明変数の個数を表している。第1項が\(y\)を補空間に射影したときの2乗和の形になっているのも面白い（実際、訓練誤差和はこのようにして表すことが可能）。

この例から、SUREにおけるリスク推定量がAICと本質的には同じような量を示していることが理解できる。何を出発として導出したかによって式形は変わってくることに注意が必要である（ここにある式形だけにとらわれてはいけない！）。

リッジ回帰の場合

Stein's unbiased risk estimatorのより面白い例は実はLASSOなのだがここではパラメータによってSUREの各項のトレードオフを簡単に見ることのできるリッジ回帰を見てみる。リッジ回帰では考えるモデルの種類は同じだが、\(\widehat{f}\)の選択に次のルールを用いる。

\[ \widehat{f}(x)=x^\top\widehat{\beta}_\lambda,\quad\widehat\beta_\lambda:=\mathop{\rm argmin}_{\beta\in\mathbb{R}^p}[(y-X\beta)^2+n\lambda\beta^2]\]

\(\lambda\)をパラメータとする正則化項が含まれていることがさきほどと異なる。これについても\(\widehat\beta_\lambda\)は陽に求めることができて、

\[ \widehat\beta_\lambda=(X^\top X+n\lambda I_p)^{-1}X^\top y\]

がその解になる。まったく同じ議論でもってhat matrixは

\[ H_\lambda=X(X^\top X+n\lambda I_p)^{-1}X^\top\]

となることもわかり、このトレースを計算することでSUREの第3項が求められる。これにはトレースの性質の他に固有値を用いる必要があって、今\(A=X^\top X\)という行列は対称行列、\(A^\top=A\)だから\(P^\top P=I_p\)なる直交行列によって\(P^\top AP=D\)と対角化することが出来る。\(D\)は行列\(A\)の固有値\(\mu_1,\dots,\mu_p\)を対角成分に並べた行列である。さらに\(A=X^\top X\)という形式の行列が正定値であるということから\(\mu_i\ge 0\)もいえる。特異値分解を考えても良いがここでは固有値でいく。

これと、トレース内で行列積の順番を自由に変えられるという性質を用いると

\[ \begin{align}\mathrm{tr}H_\lambda&=\mathrm{tr}(X[P(D+n\lambda I_p)P^\top]^{-1}X^\top)\\&=\mathrm{tr}(XP^{-\top}(D+n\lambda I_p)^{-1}P^{-1}X^\top)\\&=\mathrm{tr}(X^\top X(P^\top P)^{-1}(D+n\lambda I_p)^{-1})\\&=\mathrm{tr}(X^\top XP^\top P(D+n\lambda I_p)^{-1})\\&=\mathrm{tr}(P^\top(X^\top X) P(D+n\lambda I_p)^{-1})\\&=\mathrm{tr}(D(D+n\lambda I_p)^{-1})\\&=\sum_{i=1}^p\frac{\mu_i}{\mu_i+n\lambda}\end{align}\]

ということがわかる。（ほんとは前半こんなに式変形はいらない。最後のところは非ゼロ要素が対角成分しかないので普通の数と同じように計算するだけだ。）

\(\lambda=0\)で\(\mathrm{tr}H=p\)に帰着することも確認できる。さらに、この項が\(\lambda\)が大きくなるにつれて減少していくことも確認できる。この時、もし訓練誤差の項も\(\lambda\)が大きくなるにつれて減少するのであればその範囲内でいくらでも\(\widehat{R}\)を小さくできるということになるがどうだろうか。

\[ \frac1n((I_n-H_\lambda)y)^2\]

今\(\lambda=0\)のときにこの値が最小値をとることがわかっているから、ごく特殊な場合を除いて\(\lambda\)が大きくなるにつれて大きくなるはずだ。だからトレードオフが発生していることが理解できる。

モデルの自由度

線形回帰および、そのリッジ正則化版であるリッジ回帰においてSUREの最後の項でhat matrixのトレースを計算した。実はこの量はモデルの自由度（degrees of freedom）\(\mathrm{df}(\widehat\mu)\)として知られている。正則化が無い場合、単に

\[ \mathrm{df}(\widehat\mu)=p\]

と説明変数の数、つまり言葉のイメージ通りフリーパラメータの数となり自由度は整数以外取りようがないように思える。ところがリッジ回帰で正則化のある場合をみると

\[ \mathrm{df}(\widehat\mu)=\sum_{i=1}^p\frac{\mu_i}{\mu_i+n\lambda}\]

というように自由度という量が整数以外も取る例があることがわかる。この量はモデルの複雑さを示す量であり、ざっくり訓練誤差とのトレードオフや過剰損失などをコントロールしていると見える。ここでは線形回帰に絞ったがパラメトリック、ノンパラメトリック関わらず定義できる普遍的な量だ。

特にLASSOの場合は

\[ \widehat{\mathrm{df}}[\widehat\mu]=\text{# of numzero coefficients}\]

ということがしられておりこの辺の話はたくさんある。

さらに続く。

参考

Trevor Hastie, Robert Tibshirani, Jerome Friedman, The Elements of Statistical Learning: Data Mining, Inference, and Prediction. Second Edition. https://web.stanford.edu/~hastie/ElemStatLearn/. pp.61-68, Ridge regression.