玩不转scikit-learn（二）：Ridge Regression

为了能更好的掌握机器学习，自己边学习边自己写一写算法。

为了验证自己算法实现的准确性，我会和scikit-learn的对应算法进行对比。恩，水平有限，尽力分析下为什么scikit-learn的对应算法会比自己实现的朴素算法好。

希望可以对于算法本身和scikit-learn的实现有一定的深度理解。。希望吧

代码用Python3在jupyter notebook(ipython notebook)上写的。代码+图的完整实现我放在GitHub上了：地址如下

https://github.com/lixintong1992/Machine_Learning

求求各位在GitHub上star我。。。

---

Ridge Regression

岭回归是一种改良的最小二乘估计法，实质上是一种改良的最小二乘估计法，通过放弃最小二乘法的无偏性，以损失部分信息、降低精度为代价获得回归系数更为符合实际、更可靠的回归方法，对病态数据的拟合要强于OLS。本质是在自变量信息矩阵的主对角线元素上人为地加入一个非负因子

岭回归形式

由于直接套用线性回归可能产生过拟合，我们需要加入正则化项，如果加入的是L2正则化项，就是Ridge回归

$$E_\hat{w}=(y-X\hat{w})^T(y-X\hat{w}) + {\alpha\hat{w}\hat{w}^T}$$

这样子对上一篇blog中提到的共线性（collinearity）就有更好的鲁棒性（robust）

当线性回归模型中存在多个相关变量时，它们的系数确定性变差并呈现高方差。比如说，在一个变量上的一个很大的正系数可能被在其相关变量上的类似大小的负系数抵消

---

分析

岭回归是对最小二乘回归的一种补充，它损失了无偏性，来换取高的数值稳定性，从而得到较高的计算精度

$$\hat{w}=(X^TX+\alpha I)^{-1}X^Ty$$

当$(X^TX)$的行列式接近于0时，我们将其主对角元素都加上一个数$\alpha$，可以使矩阵为奇异的风险大降低。

特征数远大于样本数量时，不可以用最小二乘法但是可以用岭回归。

从多变量回归的变量选择来说，普通的多元线性回归要做的是变量的剔除和筛选，而岭回归是一种shrinkage的方法，就是收缩。这是什么意思呢， 比如做普通线性回归时候，如果某个变量t检验不显著，我们通常会将它剔除再做回归，如此往复（stepwise)，最终筛选留下得到一个我们满意回归方程，但是在做岭回归的时候，我们并没有做变量的剔除，而是将这个变量的系数beta向0”收缩“，使得这个变量在回归方程中的影响变的很小。 与普通的多元线性回归相比，岭回归的变化更加smooth，或者说continuous。从这点上来说活，岭回归只是shrinkage methods中的一种，大家常说的lasso回归（貌似叫套索回归）其实也属于这种方法。

---

求解

求解复杂度和原始的最小二乘法一样

sckit-learn中可以选择多种求解器：
- ‘auto’自动选择求解器
- ‘svd’ 利用奇异值分解X求解，对于奇异矩阵相较于“cholesky”更稳定
- ‘cholesky’利用的是scipy.linalg.solve可以获得闭式解（closed-form solution）
- ‘sparse_cg’共轭梯度下降，相比于“cholesky”在大数据上更合适（可以设置结果精确度和最大迭代次数）
- ‘lsqr’scipy.sparse.linalg.lsqr，这是最快的，用的是迭代
- ‘sag’随机平均梯度下降，如果样本数和特征数都很大，则表现会很快。只有特征在同一个量级才好用，可以事先处理特征

---

与普通线性回归的对比

针对特征线性相关的情况，岭回归可以比较好的求解
首先生成一个Hilbert matrix作为特征矩阵

X = 1. / (np.arange(1, 11) + np.arange(0, 10)[:, np.newaxis])  # numpy amazing!!
np.linalg.det(X)

X的行列式值2.1644528904064253e-53非常接近于0

y = np.ones(10)
y[:5]=0

随意初始化一个y值

首先我们用线性回归求解

model = linear_model.LinearRegression()
model.fit(X, y)
model.coef_

array([ -1.09983550e+06,   7.73522419e+07,  -1.30690544e+09,
         8.99906111e+09,  -2.98892606e+10,   4.77497011e+10,
        -2.36132637e+10,  -2.87825730e+10,   4.16709711e+10,
        -1.49119035e+10])

求得的系数比较爆炸啊。。
model.score(X,y)
0.99999999999572542

score是scikit-learn中的一个简单预测评价，1最好。不细说了。

此时是对于训练样本的完美拟合

model.intercept_
89644.031877040863
偏置也是十分的过分

---

对于岭回归

clf = linear_model.Ridge(alpha=0.01)
clf.fit(X, y)
clf.coef_

array([ 1.30443212, -1.89068259, -2.11762635, -1.92007014, -1.66880907,
       -1.43915365, -1.24415185, -1.08185894, -0.94716034, -0.8349426 ])

系数正常了不少

clf.score(X,y)
0.72515881280298078
对于训练集的拟合也差了不少

clf.intercept_
1.7828471405439079
偏置正常了许多

岭回归的系数正常了许多，但是代价是score低了不少，也算是防止过拟合了吧。

---

惩罚因子$\alpha$对于系数的影响实验

与最小二乘的直观对比

有了二范数的正则项，果然稳定了不少

惩罚因子$\alpha$对训练的可视化结果

惩罚系数越大，训练的误差越大

---

scikit-learn应用

交叉验证确定参数

一般来说，只要我们觉得数据有线性关系，用LinearRegression类拟合的不是特别好，需要正则化，可以考虑用Ridge类。但是这个类最大的缺点是每次我们要自己指定一个超参数 $\alpha$，然后自己评估$\alpha$的好坏，比较麻烦，一般用RidgeCV类来跑Ridge回归，RidgeCV类对超参数 $\alpha$ 使用了交叉验证，来帮忙我们选择一个合适的 $\alpha$ 。在初始化RidgeCV类时候，我们可以传一组备选的 $\alpha$ 值，10个，100个都可以。RidgeCV类会帮我们选择一个合适的 $\alpha$ 。免去了我们自己去一轮轮筛选 $\alpha$ 的苦恼。

alphas = np.logspace(-6, 6, 200)
reg = linear_model.RidgeCV(alphas=alphas)
reg.fit(X, y)
reg.alpha_

9.9999999999999995e-07
这样子就能从中确定最好的alpha_。这里用的是Generalized Cross-Validation (GCV)，an efficient form of leave-one-out cross-validation。具体不懂。。。

---

自己的实现

scikit-learn中岭回归求解方法多样：
‘auto’自动选择求解器
‘svd’ 利用奇异值分解X求解，对于奇异矩阵相较于“cholesky”更稳定
‘cholesky’利用的是scipy.linalg.solve可以获得闭式解（closed-form solution）
‘sparse_cg’共轭梯度下降，相比于“cholesky”在大数据上更合适（可以设置结果精确度和最大迭代次数）
‘lsqr’scipy.sparse.linalg.lsqr，这是最快的，用的是迭代
‘sag’随机平均梯度下降，如果样本数和特征数都很大，则表现会很快。只有特征在同一个量级才好用，可以事先处理特征

这么多求解方法。。我实现最简单的闭式解（closed-form solution）吧。。
就是根据下面的公式就可以了

$$\hat{w}=(X^TX+\alpha I)^{-1}X^Ty$$

GitHub上有具体的实现，结果总是和scikit-learn差一点点，不明原因差不多吧。。。主要还是solver='cholesky'用的是scipy.linalg.solve
源代码相关部分如下：

下面试试看scipy.linalg.solve

scipy.linalg.solve求解的是$ax=b$，其中$a$是方阵，求解$x$

$$(X^TX+\alpha I)\hat{w}=X^Ty$$

最后放弃了。。不知道为什么总是差一点。。。反正就差一点的说。。

---

结尾

https://github.com/lixintong1992/Machine_Learning

GitHub上求star~写blog不容易呀~