Ridge Regression - 岭回归

Why 目的

predictors X之间存在严重的多重共线性（multicollinarity, 即自变量之间线性相关-correlation很高）时，会使得least-square（最小二乘法）计算公式β^=(XTX)−1XTY中的R=(XTX) 不可逆（成为了singulation 奇异值 即|XTX| = 0）从而最小二乘法失效。

What 本质是什么

给R内部加上一个惩罚项，使得R成为singulation的可能性大大降低。
本质还是bias – variance tradeoff，这次是牺牲unbias换取小variance。

How 具体是怎么解决的

先回顾一下。||a⃗ ||=(a1−0)2+(a2−0)2+...+(an−0)2−−−−−−−−−−−−−−−−−−−−−−−−−−−−−√表n维向量a⃗ 在笛卡尔坐标系的长度。因此最小二乘法的最小化的目标函数表示为||Y−Xβ||2，得出的解为β^=(XTX)−1XTY。但当XTX存在很严重的多重共线性时，XTX会接近奇异，||β^||值会变得很大，即使||β||本身并不大的情况下。因此我们可以把目标函数变为||Y−Xβ||2+k||β||2,k>0.时，引入的正惩罚项k||β||2使得得出的β^长度不会太大，从而解决了严重多重共线性引入的问题。 注意要求X之前都standardized 了（X−X¯sd(X)）使得X们处于同一单位。
β(k)^=(R+kI)−1X′Y

背景

为什么R=XTX接近奇异时会导致||β^||变很大呢？下面的内容我也看不懂，但日本数学家说过看不懂可以抄。对R作谱分解，得到R=Σ1≤j≤pλjvjvTj, 其中λ为特征值，vj为对应单位特征向量(vTjvj=1)。

E[||β^||2]=E[||β+(β^−β)||2]=||β2||+E[tr(β^−β)T(β^−β)]=||β2||+tr[Cov(β^)]=||β2||+σ2tr[R−1]=||β2||+σ2Σ1≤j≤p1/λj

。当R秩为q<p时非满秩，作为奇异矩阵其最小的p−q个特征值为0。在这种情况下不存在唯一的最小二乘解，normal equation有无数个解且Σp−q+1≤j≤p1/λj=∞。相应的当R接近奇异时Σp−q+1≤j≤p1/λj虽然不是正无穷但也会变得特别大，因此在||β||2并不大时E[||β^||2]依然会很大。通过马尔可夫不等式我们可以知道||β^||2很大的概率还很高。因此我们需要变更最小化目标函数使得||β^||2尽可能小。

记得我们前面说过岭回归牺牲了一部分无偏差性，使得variance变小。现在我们来看看数学推导。让β(k)=E[β(k)^]=E[(R+kI)−1X′Y]=(R+kI)−1Rβ=β−k(R+kI)−1β
Cov(β^(k))=σ2(R+kI)−1R(R+kI)−1
Bias-variance tradeoff: SSE(k) = ||Y - X\hat{\beta}(k)||^2SSE(k)=||Y−Xβ^(k)||2, H(k) = X(R + kI)^{-1}X^TH(k)=X(R+kI)−1XT, \mu(k) = X\beta(k)μ(k)=Xβ(k)
E(SSE(k)) = \sigma^2 tr((I - H(k) )) + ||\mu - \mu(k)||^2E(SSE(k))=σ2tr((I−H(k)))+||μ−μ(k)||2

bias-variance tradeoff同样反映在了组成model的各部分。我们先来看看\beta(k)β(k)在岭回归下的bias-variance tradeoff.
Risk: E(||\beta - \hat{\beta}(k)||^2) = E(||\hat{\beta(k)} - \beta(k) + \beta(k) + \beta||^2) = E(||\hat{\beta(k)} - \beta(k)||^2) + ||\beta(k) - \beta||^2 + 2E[(\hat{\beta}(k) - \beta(k))'(\beta(k) - \beta)]

E(||β−β^(k)||2)=E(||β(k)^−β(k)+β(k)+β||2)=E(||β(k)^−β(k)||2)+||β(k)−β||2+2E[(β^(k)−β(k))′(β(k)−β)]

因为实际上只有\hat{\beta}(k)β^(k)带有YY是随机变量，其余的都是常量。所以E(||\beta - \hat{\beta}(k)||^2) = E(||\hat{\beta(k)} - \beta(k)||^2) + ||\beta(k) - \beta||^2 = variance + bias^2 E(||β−β^(k)||2)=E(||β(k)^−β(k)||2)+||β(k)−β||2=variance+bias2
我们在前面使用谱分解算过variance 部分，E(||β(k)^−β(k)||2)=E[tr(β^+β)T(β^+β)]=tr[Cov(β^)]=σ2Σ1≤j≤pλj/(k+λj).
接着是bias部分，||β(k)−β||2=||k(R+kI)−1β||2=k2Σ(vTjβ)2/(k+λj)2.
让variance部分对k求导，得到variance随着k单调减，而bias是单调增。
因此当k增大时，variance降低而bias增大。任意k>0时，存在E(||β(k)^−β(k)||2)<E(||β(0)^−β(0)||2)，也就是说比普通最小二乘法的 β^ variance小。
把variance-bias结合起来得到当k≈0 时，有E(||β−β^(k)||2)<E(||β−β^(0)||2)。

同理我们还可以对μ(k)=Xβ(k) 做这样的分析。

E(||Xβ−Xβ^(k)||2)=E(||Xβ−Xβ(k)+Xβ(k)+Xβ^(k)||2)=E(||Xβ(k)+Xβ^(k)||2)+||Xβ−Xβ(k)||2=σ2Σ1≤j≤pλ2j/(k+λj)+k2Σ(vTjβ)2λj/(k+λj)2

.
和上节相比分母没变，只是分子有点小变化。

R

X = as.matrix(crime.s.train[, -1])D = t(X) %*% Xeigen(D)eigen(D)[[1]][1]/sum(eigen(D)[[1]])library('MASS')select(lm.ridge(Crime ~ 0 + ., data = crime.s.train,lambda = seq(0,10, .01)))## Model with k = 2.48:k.ridge = 2.48library('ridge')fit.ridge = linearRidge(Crime ~ 0 + ., data = crime.s.train, lambda = k.ridge, scaling = 'none'); summary(fit.ridge)#fitted valuefitted.ridge = predict(fit.ridge, newx = crime.s.train[, 2:13], s = k.ridge)#residualsresid.ridge = crime.s.train[ ,1] - fitted.ridgepar(mfrow = c(1, 3))plot(crime.s.train$Crime, fitted.ridge , xlab = "Crime rate", ylab = " fitted Crime rate", main = "ridge method: Y vs Fitted Y")plot(resid.ridge, xlab = "", main = "ridge method: plot of resid")qqnorm(resid.ridge); qqline(resid.ridge)par(mfrow = c(1, 1))D.ridge = t(X) %*% X + k.ridge * diag(12)VIF.ridge = diag(solve(D.ridge) %*% t(X) %*% X %*% solve(D.ridge)) VIF.ridge