今回は、リッジ回帰と特異値分解 (Singular Value Decomposition, SVD) の関係に加え、パラメーター最適化における一般化交差確認 (Generalized Cross-Validation, GCV) についてまとめました 。間違いなどありましたら、ご指摘いただけると幸いです。
最小二乗線形回帰 (OLS) とリッジ回帰 λ \lambda λ リッジ回帰とSVD λ \lambda λ Cross-Validation の効率的な計算 Generalized Cross-Validation
最小二乗線形回帰 (OLS) とリッジ回帰
例えば、説明変数 X ∈ R n × p X \in \R^{n \times p} X ∈ R n × p y ∈ R n \boldsymbol{y} \in \R^{n} y ∈ R n β ∈ R p \boldsymbol{\beta} \in \R^{p} β ∈ R p ε ∈ R n ∼ N ( 0 , σ 2 I ) \boldsymbol{\varepsilon} \in \R^{n} \sim \mathcal{N}(0,\sigma^2 I) ε ∈ R n ∼ N ( 0 , σ 2 I )
y = X β + ε ( 1 ) \boldsymbol{y} = X \boldsymbol{\beta} + \boldsymbol{\varepsilon} \quad (1) y = X β + ε ( 1 ) この回帰係数の最小二乗推定量 β ^ \hat{\boldsymbol{\beta}} β ^
β ^ = ( X ⊤ X ) − 1 X ⊤ y ( 2 ) \hat{\boldsymbol{\beta}} = (X^\top X)^{-1} X^\top \boldsymbol{y} \quad (2) β ^ = ( X ⊤ X ) − 1 X ⊤ y ( 2 ) しかし、X X X X ⊤ X X^\top X X ⊤ X ridge λ \lambda λ
β ^ λ = ( X ⊤ X + λ I ) − 1 X ⊤ y ( 3 ) \hat{\boldsymbol{\beta}}_{\lambda} = (X^\top X + \lambda I)^{-1} X^\top \boldsymbol{y} \quad (3) β ^ λ = ( X ⊤ X + λ I ) − 1 X ⊤ y ( 3 ) (3) 式が、Hoerl and Kennard, 1970 によって提案されたリッジ回帰 の方法です 。これを最適化問題として書くと次のようになります。
argmin β ∥ y − X β ∥ 2 2 + λ ∥ β ∥ 2 2 \underset{\boldsymbol{\beta}}{\text{argmin }} \| \boldsymbol{y}-X\boldsymbol{\beta} \|_2^2 + \lambda \|\boldsymbol{\beta}\|_2^2 β argmin ∥ y − X β ∥ 2 2 + λ ∥ β ∥ 2 2 ここで、∥ ∥ 2 2 \|\|_2^2 ∥ ∥ 2 2 ℓ 2 \ell_2 ℓ 2 ℓ 2 \ell_2 ℓ 2 λ ∥ β ∥ 2 2 \lambda \|\boldsymbol{\beta}\|_2^2 λ ∥ β ∥ 2 2 β \boldsymbol{\beta} β β ^ λ \hat{\boldsymbol{\beta}}_{\lambda} β ^ λ
− 2 X ⊤ ( y − X β ^ ) + 2 λ β ^ = 0 ( X ⊤ X + λ I ) β ^ = X ⊤ y β ^ = ( X ⊤ X + λ I ) − 1 X ⊤ y \begin{aligned}
-2X^\top(\boldsymbol{y}-X\hat{\boldsymbol{\beta}}) + 2\lambda \hat{\boldsymbol{\beta}} &= 0 \\
(X^\top X + \lambda I)\hat{\boldsymbol{\beta}} &= X^\top \boldsymbol{y} \\
\hat{\boldsymbol{\beta}} &= (X^\top X + \lambda I)^{-1} X^\top \boldsymbol{y}
\end{aligned} − 2 X ⊤ ( y − X β ^ ) + 2 λ β ^ ( X ⊤ X + λ I ) β ^ β ^ = 0 = X ⊤ y = ( X ⊤ X + λ I ) − 1 X ⊤ y λ \lambda λ ハイパーパラメーターとなる λ \lambda λ 数値計算的には、λ = 0.001 \lambda=0.001 λ = 0.001 λ \lambda λ
Validation dataset で各 λ \lambda λ
Training dataset 内で Cross-validation をして各 λ \lambda λ
C p C_p C p
ここでは、λ \lambda λ リッジ回帰とSVDの関係 を始めに考えます。
リッジ回帰とSVD
X X X 特異値分解 (SVD) を X = U D V ⊤ X = UDV^\top X = U D V ⊤ U = ( u 1 , … , u n ) ∈ R n × n , U = (\boldsymbol{u}_1, \dots, \boldsymbol{u}_n) \in \R^{n \times n}, U = ( u 1 , … , u n ) ∈ R n × n , V = ( v 1 , … , v n ) ∈ R p × p V = (\boldsymbol{v}_1, \dots, \boldsymbol{v}_n) \in \R^{p \times p} V = ( v 1 , … , v n ) ∈ R p × p D ∈ R n × p D \in \R^{n \times p} D ∈ R n × p d 1 , … , d m ( m = min ( n , p ) ) d_1, \dots, d_{m}\ (m=\min(n,p)) d 1 , … , d m ( m = min ( n , p ))
これを式 (3) に代入すると、リッジ回帰推定量 β ^ λ \hat{\boldsymbol{\beta}}_{\lambda} β ^ λ 。
β ^ λ = ( X ⊤ X + λ I ) − 1 X ⊤ y = ( ( U D V ⊤ ) ⊤ U D V ⊤ + λ I ) − 1 ( U D V ⊤ ) ⊤ y = ( V D ⊤ U ⊤ U D V ⊤ + λ V V ⊤ ) − 1 V D ⊤ U ⊤ y = V ( D ⊤ D + λ I ) − 1 D ⊤ U ⊤ y = V ( d 1 d 1 2 + λ O ⋱ O d m d m 2 + λ ) U ⊤ y = ∑ d j > 0 v j d j d j 2 + λ ⟨ u j , y ⟩ \begin{aligned}
\hat{\boldsymbol{\beta}}_{\lambda} &= (X^\top X + \lambda I)^{-1} X^\top \boldsymbol{y} \\
&= ((UDV^\top)^\top UDV^\top + \lambda I )^{-1} (UDV^\top)^\top \boldsymbol{y} \\
&= (VD^\top U^\top UDV^\top + \lambda V V^\top)^{-1} VD^\top U^\top \boldsymbol{y} \\
&= V (D^\top D + \lambda I )^{-1} D^\top U^\top \boldsymbol{y} \\
&= V \begin{pmatrix}
\frac{d_1}{d_1^2 + \lambda} & & O \\
& \ddots & \\
O & & \frac{d_{m}}{d_{m}^2 + \lambda}
\end{pmatrix} U^\top \boldsymbol{y} \\
&= \sum_{d_j \gt 0} \boldsymbol{v}_j \frac{d_j}{d_j^2 + \lambda} \langle \boldsymbol{u}_j, \boldsymbol{y} \rangle
\end{aligned} β ^ λ = ( X ⊤ X + λ I ) − 1 X ⊤ y = (( U D V ⊤ ) ⊤ U D V ⊤ + λ I ) − 1 ( U D V ⊤ ) ⊤ y = ( V D ⊤ U ⊤ U D V ⊤ + λV V ⊤ ) − 1 V D ⊤ U ⊤ y = V ( D ⊤ D + λ I ) − 1 D ⊤ U ⊤ y = V ⎝ ⎛ d 1 2 + λ d 1 O ⋱ O d m 2 + λ d m ⎠ ⎞ U ⊤ y = d j > 0 ∑ v j d j 2 + λ d j ⟨ u j , y ⟩ ただし、⟨ ⋅ , ⋅ ⟩ \langle\cdot,\cdot\rangle ⟨ ⋅ , ⋅ ⟩ y ^ λ \hat{\boldsymbol{y}}_\lambda y ^ λ
y ^ λ = X β ^ λ = ( U D V ⊤ ) β ^ λ = U D ( D ⊤ D + λ I ) − 1 D ⊤ U ⊤ y = ∑ d j > 0 u j d j 2 d j 2 + λ ⟨ u j , y ⟩ \begin{aligned}
\hat{\boldsymbol{y}}_\lambda &= X \hat{\boldsymbol{\beta}}_{\lambda} = (UDV^\top) \hat{\boldsymbol{\beta}}_{\lambda} \\
&= U D (D^\top D + \lambda I )^{-1} D^\top U^\top \boldsymbol{y} \\
&= \sum_{d_j \gt 0} \boldsymbol{u}_j \frac{d_j^2}{d_j^2 + \lambda} \langle \boldsymbol{u}_j, \boldsymbol{y} \rangle
\end{aligned} y ^ λ = X β ^ λ = ( U D V ⊤ ) β ^ λ = U D ( D ⊤ D + λ I ) − 1 D ⊤ U ⊤ y = d j > 0 ∑ u j d j 2 + λ d j 2 ⟨ u j , y ⟩ λ \lambda λ リッジ回帰の係数は、λ \lambda λ X X X n > p n \gt p n > p β ^ λ \hat{\boldsymbol{\beta}}_{\lambda} β ^ λ
Bias ( β ^ λ ) = E [ β ^ λ ] − β = E [ ( X ⊤ X + λ I ) − 1 X ⊤ y ] − β = ( X ⊤ X + λ I ) − 1 X ⊤ E [ y ] − β = ( X ⊤ X + λ I ) − 1 X ⊤ X β − β = β − λ ( X ⊤ X + λ I ) − 1 β − β = λ V ( D ⊤ D + λ I ) − 1 V ⊤ y = ∑ j = 1 p v j λ d j 2 + λ ⟨ v j , y ⟩ \begin{aligned}
\text{Bias}(\hat{\boldsymbol{\beta}}_{\lambda}) &= \mathbb{E}[\hat{\boldsymbol{\beta}}_{\lambda}]-\boldsymbol{\beta} \\
&= \mathbb{E}[(X^\top X + \lambda I)^{-1} X^\top \boldsymbol{y}]-\boldsymbol{\beta} \\
&= (X^\top X + \lambda I)^{-1} X^\top \mathbb{E}[\boldsymbol{y}]-\boldsymbol{\beta} \\
&= (X^\top X + \lambda I)^{-1} X^\top X \boldsymbol{\beta}-\boldsymbol{\beta} \\
&= \boldsymbol{\beta} - \lambda (X^\top X + \lambda I)^{-1} \boldsymbol{\beta}-\boldsymbol{\beta} \\
&= \lambda V (D^\top D + \lambda I )^{-1} V^\top \boldsymbol{y} \\
&= \sum_{j=1}^p \boldsymbol{v}_j \frac{\lambda}{d_j^2 + \lambda} \langle \boldsymbol{v}_j, \boldsymbol{y} \rangle
\end{aligned} Bias ( β ^ λ ) = E [ β ^ λ ] − β = E [( X ⊤ X + λ I ) − 1 X ⊤ y ] − β = ( X ⊤ X + λ I ) − 1 X ⊤ E [ y ] − β = ( X ⊤ X + λ I ) − 1 X ⊤ X β − β = β − λ ( X ⊤ X + λ I ) − 1 β − β = λV ( D ⊤ D + λ I ) − 1 V ⊤ y = j = 1 ∑ p v j d j 2 + λ λ ⟨ v j , y ⟩ となって、λ > 0 \lambda \gt 0 λ > 0 λ \lambda λ n < p n \lt p n < p
一方、W λ = ( X ⊤ X + λ I ) − 1 X ⊤ X W_\lambda=(X^\top X + \lambda I)^{-1} X^\top X W λ = ( X ⊤ X + λ I ) − 1 X ⊤ X β ^ λ \hat{\boldsymbol{\beta}}_{\lambda} β ^ λ
Var [ β ^ λ ] = Var [ W λ β ^ ] = W λ Var [ β ^ ] W λ ⊤ = W λ σ 2 ( X ⊤ X ) − 1 W λ ⊤ = σ 2 { ( X ⊤ X + λ I ) − 1 X ⊤ X } ( X ⊤ X ) − 1 { ( X ⊤ X + λ I ) − 1 X ⊤ X } ⊤ = σ 2 ( X ⊤ X + λ I ) − 1 X ⊤ X { ( X ⊤ X + λ I ) − 1 } ⊤ = σ 2 ∑ j = 1 p d j 2 ( d j 2 + λ ) 2 v j v j ⊤ \begin{aligned}
\text{Var}[\hat{\boldsymbol{\beta}}_{\lambda}] &= \text{Var}[W_\lambda \hat{\boldsymbol{\beta}}] = W_\lambda \text{Var}[\hat{\boldsymbol{\beta}}] W_\lambda^\top \\
&= W_\lambda\sigma^2 (X^\top X)^{-1} W_\lambda^\top \\
&= \sigma^2 \{(X^\top X + \lambda I)^{-1} X^\top X\} (X^\top X)^{-1}\{(X^\top X + \lambda I)^{-1} X^\top X\}^\top \\
&= \sigma^2 (X^\top X + \lambda I)^{-1} X^\top X \{(X^\top X + \lambda I)^{-1}\}^\top \\
&= \sigma^2 \sum_{j=1}^p \frac{d_j^2}{(d_j^2 + \lambda)^2} \boldsymbol{v}_j \boldsymbol{v}_j^\top
\end{aligned} Var [ β ^ λ ] = Var [ W λ β ^ ] = W λ Var [ β ^ ] W λ ⊤ = W λ σ 2 ( X ⊤ X ) − 1 W λ ⊤ = σ 2 {( X ⊤ X + λ I ) − 1 X ⊤ X } ( X ⊤ X ) − 1 {( X ⊤ X + λ I ) − 1 X ⊤ X } ⊤ = σ 2 ( X ⊤ X + λ I ) − 1 X ⊤ X {( X ⊤ X + λ I ) − 1 } ⊤ = σ 2 j = 1 ∑ p ( d j 2 + λ ) 2 d j 2 v j v j ⊤ と誤差の分散 σ 2 \sigma^2 σ 2 λ \lambda λ これらがバイアス・バリアンスのトレードオフ に繋がります 。
すなわち、λ \lambda λ
補足：最小二乗推定量のバリアンス
Var [ β ^ ] = Var [ { β ^ − E [ β ] } { β ^ − E [ β ] } ⊤ ] = E [ { ( X ⊤ X ) − 1 X ⊤ y − β } { ( X ⊤ X ) − 1 X ⊤ y − β } ⊤ ] = ( X ⊤ X ) − 1 X ⊤ E [ Y Y ⊤ ] X ( X ⊤ X ) − 1 − β β ⊤ = ( X ⊤ X ) − 1 X ⊤ { X β β ⊤ X ⊤ + Var ( ε i ) I } X ( X ⊤ X ) − 1 − β β ⊤ = β β ⊤ + σ 2 ( X ⊤ X ) − 1 − β β ⊤ = σ 2 ( X ⊤ X ) − 1 \begin{aligned}
\text{Var}[\hat{\boldsymbol{\beta}}] &= \text{Var}[\{\hat{\boldsymbol{\beta}}-\mathbb{E}[\boldsymbol{\beta}]\}\{\hat{\boldsymbol{\beta}}-\mathbb{E}[\boldsymbol{\beta}]\}^\top]\\
&= \mathbb{E}[\{(X^\top X)^{-1} X^\top \boldsymbol{y} -\boldsymbol{\beta}\}\{(X^\top X)^{-1} X^\top \boldsymbol{y} -\boldsymbol{\beta}\}^\top]\\
&= (X^\top X)^{-1} X^\top \mathbb{E}[YY^\top] X (X^\top X)^{-1} - \boldsymbol{\beta}\boldsymbol{\beta}^\top \\
&= (X^\top X)^{-1} X^\top \{X\boldsymbol{\beta}\boldsymbol{\beta}^\top X^\top + \text{Var}(\varepsilon_i)I\} X (X^\top X)^{-1} - \boldsymbol{\beta}\boldsymbol{\beta}^\top\\
&= \boldsymbol{\beta}\boldsymbol{\beta}^\top + \sigma^2 (X^\top X)^{-1} - \boldsymbol{\beta}\boldsymbol{\beta}^\top \\
&= \sigma^2 (X^\top X)^{-1}
\end{aligned} Var [ β ^ ] = Var [{ β ^ − E [ β ]} { β ^ − E [ β ] } ⊤ ] = E [{( X ⊤ X ) − 1 X ⊤ y − β } {( X ⊤ X ) − 1 X ⊤ y − β } ⊤ ] = ( X ⊤ X ) − 1 X ⊤ E [ Y Y ⊤ ] X ( X ⊤ X ) − 1 − β β ⊤ = ( X ⊤ X ) − 1 X ⊤ { X β β ⊤ X ⊤ + Var ( ε i ) I } X ( X ⊤ X ) − 1 − β β ⊤ = β β ⊤ + σ 2 ( X ⊤ X ) − 1 − β β ⊤ = σ 2 ( X ⊤ X ) − 1
Cross-Validation の効率的な計算
λ \lambda λ λ \lambda λ Cross-Validation の利用があります。一般には、k-fold Cross-Validation を利用して分割した k k k C p C_p C p
n-fold Cross-Validation (= Leave One Out Cross-Validation, LOOCV ) の誤差 LOO λ \text{LOO}_\lambda LOO λ n n n X = ( x 1 , … , x n ) ⊤ , X=(\boldsymbol{x}_1, \dots, \boldsymbol{x}_n)^{\top}, X = ( x 1 , … , x n ) ⊤ , y = ( y 1 , … , y n ) \ \boldsymbol{y}=(y_1, \dots, y_n) y = ( y 1 , … , y n )
LOO λ = ∑ i = 1 n ( y i − x i ⊤ β ^ λ ( − i ) ) 2 = ∑ i = 1 n ( y i − x i ⊤ β ^ λ ) 2 ( 1 − R i i λ ) 2 ( 4 ) \begin{aligned}
\text{LOO}_\lambda &= \sum_{i=1}^n (y_i - \boldsymbol{x}_i^\top \hat{\boldsymbol{\beta}}_{\lambda}^{(-i)})^2 \\
&= \sum_{i=1}^n \frac{(y_i - \boldsymbol{x}_i^\top \hat{\boldsymbol{\beta}}_{\lambda})^2}{(1-R_{ii}^\lambda)^2} \quad (4)
\end{aligned} LOO λ = i = 1 ∑ n ( y i − x i ⊤ β ^ λ ( − i ) ) 2 = i = 1 ∑ n ( 1 − R ii λ ) 2 ( y i − x i ⊤ β ^ λ ) 2 ( 4 ) ここで、β ^ λ ( − i ) \hat{\boldsymbol{\beta}}_{\lambda}^{(-i)} β ^ λ ( − i ) i i i R λ R^\lambda R λ リッジ演算子行列 です。
R λ = X ( X ⊤ X + λ I ) − 1 X ⊤ R^\lambda = X(X^\top X + \lambda I)^{-1}X^\top R λ = X ( X ⊤ X + λ I ) − 1 X ⊤ したがって、LOO λ \text{LOO}_\lambda LOO λ β ^ λ \hat{\boldsymbol{\beta}}_{\lambda} β ^ λ R λ R^\lambda R λ
R λ = V ( D ⊤ D + λ I ) − 1 D ⊤ U ⊤ = U ( d 1 2 d 1 2 + λ O ⋱ O d r d r 2 + λ ) U ⊤ = U S λ U \begin{aligned}
R^\lambda &= V (D^\top D + \lambda I )^{-1} D^\top U^\top \\
&= U \begin{pmatrix}
\frac{d_1^2}{d_1^2 + \lambda} & & O \\
& \ddots & \\
O & & \frac{d_{r}}{d_{r}^2 + \lambda}
\end{pmatrix} U^\top = U S_\lambda U
\end{aligned} R λ = V ( D ⊤ D + λ I ) − 1 D ⊤ U ⊤ = U ⎝ ⎛ d 1 2 + λ d 1 2 O ⋱ O d r 2 + λ d r ⎠ ⎞ U ⊤ = U S λ U Generalized Cross-Validation
(4) 式で、tr [ R λ ] = ∑ i = 1 n R i i λ \text{tr}[R^\lambda]=\sum_{i=1}^n R^\lambda_{ii} tr [ R λ ] = ∑ i = 1 n R ii λ
R i i λ ≃ 1 n tr [ R λ ] R^\lambda_{ii} \simeq \frac{1}{n}\text{tr}[R^\lambda] R ii λ ≃ n 1 tr [ R λ ] という近似を利用すると、誤差 GCV λ \text{GCV}_\lambda GCV λ この誤差を利用する方法は、Generalized Cross-Validation と呼ばれています 。
GCV λ = ∑ i = 1 n ( y i − x i ⊤ β ^ λ ) 2 ( 1 − 1 n tr [ R λ ] ) 2 ( 5 ) \text{GCV}_\lambda = \sum_{i=1}^n \frac{(y_i - \boldsymbol{x}_i^\top \hat{\boldsymbol{\beta}}_{\lambda})^2}{(1-\frac{1}{n}\text{tr}[R^\lambda])^2} \quad (5) GCV λ = i = 1 ∑ n ( 1 − n 1 tr [ R λ ] ) 2 ( y i − x i ⊤ β ^ λ ) 2 ( 5 ) この方法は、個々の対角要素 R i i λ R^\lambda_{ii} R ii λ tr [ R λ ] \text{tr}[R^\lambda] tr [ R λ ] GCV λ \text{GCV}_\lambda GCV λ λ \lambda λ 漸近的に最適である と示されています
Reference
Trevor, T. (2020). Ridge Regularization: An Essential Concept in Data Science. Technometrics. 62:4. 426-433. url
Hastie, T. (2020). Ridge Regularizaton: an Essential Concept in Data Science. arXiv preprint arXiv:2006.00371. url
van Wieringen, W. N. (2015). Lecture notes on ridge regression. arXiv preprint arXiv:1509.09169. url
Hoerl, A. E., & Kennard, R. W. (1970). Ridge regression: Biased estimation for nonorthogonal problems. Technometrics, 12(1), 55-67. url
Li, K. C. (1986). Asymptotic optimality of C L C_L C L url
荒木孝治 2013 罰則付き回帰とデータ解析環境R 公益社団法人日本オペレーションズ・リサーチ学会 url
読了：Hastie (2020) リッジ正則化についてこれでもかこれでもかと語り倒す url
