Python手写Lasso回归算法

理论

Lasso回归在最小二乘法的基础上加上了一个$l_1$惩罚项

损失函数：$$J(\theta )=\frac{1}{2m}\sum _{i=1}^m(h_{\theta }(x^{(i)})-y^{(i)}{)}^2+\lambda \sum _{j=1}^n|{\theta }_j|$$

相比岭回归可以直接通过矩阵运算得到回归系数相比，LASSO的计算变得相对复杂。由于惩罚项中含有绝对值，此函数的导数是连续不光滑的，所以无法进行求导并使用梯度下降优化。

两种求解Lasso回归的方法：

• 坐标下降法
• 最小角回归

坐标下降法

坐标下降法，是沿着坐标轴的方向去下降。

坐标下降法的数学依据是：

         一个可微的凸函数$J(\theta )$，其中$\theta$是$n\ast 1$的向量，即有$n$个维度。如果在某一点$\overline{\theta }$，使得$J(\theta )$在每一个坐标轴${\overline{\theta }}_i(i=1,2,...,n)$上都是最小值，那么$J({\overline{\theta }}_i)$就是一个全局最小值

于是，我们的优化目标是：在$\theta$的$n$个坐标轴上，对损失函数做迭代的下降，当所有的坐标轴上的${\theta }_i(i=1,2,...,n)$都收敛，此时损失函数最小，此时的$\theta$即为我们要求的结果。

具体算法流程：

1、首先，初始化$\theta$向量，随机取值即可，即为${\theta }^{(0)}$，上面的括号里的数字表示当前迭代的轮数。

2、对于第$k$轮的迭代，我们从${\theta }_1^{(k)}$开始，到${\theta }_n^{(k)}$为止，依次求${\theta }_i^{(k)}$。${\theta }_i^{(k)}$的表达式如下：

${\theta }_i^{(k)}\in \underset{{\theta }_i}{\underbrace{argmin}}J({\theta }_1^{(k)},{\theta }_2^{(k)},...{\theta }_{i-1}^{(k)},{\theta }_i,{\theta }_{i+1}^{(k-1)},...,{\theta }_n^{(k-1)})$

$\theta$向量的$n$个维度的迭代式如下：

${\theta }_i^{(k)}\in \underset{{\theta }_i}{\underbrace{argmin}}J({\theta }_1^{(k)},{\theta }_2^{(k)},...{\theta }_{i-1}^{(k)},{\theta }_i,{\theta }_{i+1}^{(k-1)},...,{\theta }_n^{(k-1)})$

也就是说${\theta }_i^{(k)}$是使得$J({\theta }_1^{(k)},{\theta }_2^{(k)},...{\theta }_{i-1}^{(k)},{\theta }_i,{\theta }_{i+1}^{(k-1)},...,{\theta }_n^{(k-1)})$最小化的${\theta }_i$的值。此时$J(\theta )$只有${\theta }_i^{(k)}$是变量，其余都是常量，所以这就是一个$J(\theta )$关于${\theta }_i^{(k)}$的一元函数，很容易通过求导求得最小值。

如果上面这个式子不好理解，我们具体一点，在第$k$轮，$\theta$向量的$n$个维度的迭代式如下：

${\theta }_1^{(k)}\in \underset{{\theta }_1}{\underbrace{argmin}}J({\theta }_1,{\theta }_2^{(k-1)},...,{\theta }_n^{(k-1)})$

${\theta }_2^{(k)}\in \underset{{\theta }_2}{\underbrace{argmin}}J({\theta }_1^{(k)},{\theta }_2,{\theta }_3^{(k-1)}...,{\theta }_n^{(k-1)})$

$...$

${\theta }_n^{(k)}\in \underset{{\theta }_n}{\underbrace{argmin}}J({\theta }_1^{(k)},{\theta }_2^{(k)},...,{\theta }_{n-1}^{(k)},{\theta }_n)$

3、检查${\theta }^{(k)}$向量和${\theta }^{(k-1)}$向量在各个维度上的变化情况，如果在所有维度上变化都足够小，则${\theta }^{(k)}$为最终结果，否则转入2，继续迭代

最小角回归

在介绍最小角回归之前，需要先看看两个预备算法：

• 前向选择算法
• 前向梯度算法

前向选择算法

前向选择算法的原理是一种典型的贪心算法。要解决的问题是：

        对于$Y=X\theta$这样的线性关系，如何求解系数$\theta$。其中$Y$是$m\ast 1$的向量，$X$是$m\ast n$的矩阵，$\theta$为$n\ast 1$的向量。$m$为样本数量，$n$为特征维度。

把矩阵$X$看成$n$个$m\ast 1$的向量$X_i(i=1,2,...,n)$。在这$n$个向量中选择一个与目标$Y$的余弦距离最大的一个$X_k$，用$X_k$来逼近$Y$，得到下式：
$\overline{Y}=X_k{\theta }_k$，其中${\theta }_k=\frac{<X_k,Y>}{||X_k|{|}_2}$

即$\overline{Y}$是$Y$在$X_k$上的投影。那么，可以定义残差：$Y_{yes}=Y-\overline{Y}$。由于是投影，可知$Y_{yes}$和$X_k$是正交的。再以$Y_{yes}$作为新的因变量，去掉$X_k$后的剩下的自变量的集合$X_i(i=1,2,...,k-1,k+1,...,n)$作为新的自变量集合，重复刚才投影和残差的操作，直到残差为0，或者所有的自变量都用完了，才停止算法。

当$X$只有2维时，如上图所示，和$Y$最接近的是$X_1$，首先在$X_1$上投影，残差如上图长虚线。此时$X_1{\theta }_1$模拟了$Y$，${\theta }_1$模拟了$\theta$（仅仅模拟了一个维度）。接着发现最接近的是$X_2$，此时用残差接着在$X_2$投影，残差为图中短虚线。由于没有其他自变量了，此时$x_1{\theta }_1+x_2{\theta }_2$模拟了$Y$，对应的模拟了两个维度的$\theta$即为最终结果。

此算法对每个变量只需执行一次操作，效率高，运算快。但，当自变量不是正交的时候，每次都在做投影，所以算法只能给出一个局部近似解。这个简单的算法太粗糙，不能直接用于Lasso回归。

前向梯度算法

前向梯度算法和前向选择算法有类似的地方，也是在$n$个$X_i$中选择和目标$Y$最为接近(余弦距离最大)的一个变量$X_k$，用$X_k$来逼近$Y$。但前向梯度算法不是粗暴的用投影，而是每次在最为接近的自变量$X_t$的方向移动一小步，然后再看残差$Y_{yes}$和哪个$X_i$最为接近。此时我们也不会把$X_t$去除，因为我们只前进了一小步，有可能下面最接近的自变量还是$X_t$。如此进行下去，直到残差$Y_{y}es$减小到足够小，算法停止。

当$X$只有2维时，例子如上图，和$Y$最接近的是$X_1$，首先在$X_1$上面走一小段距离，此处$ϵ$为一个较小的常量，发现此时的残差还是和$X_1$最接近。那么接着沿$X_1$走，一直走到发现残差不是和$X_1$最接近，而是和$X_2$最接近，此时残差如上图长虚线。接着沿着$X_2$走一小步，发现残差此时又和$X_1$最接近，那么开始沿着$X_1$走，走完一步后发现残差为0，那么算法停止。此时$Y$由刚才所有的所有步相加而模拟，对应的算出的系数$\theta$即为最终结果。

最小角回归算法

最小角回归对前向梯度和前向选择做了这种，保留了前向梯度算法一定程度上的精确性，同时简化了前向梯度算法一步步迭代的过程：

首先，找到与因变量$Y$最接近或相关度最高的自变量$X_k$，使用类似于前向梯度算法中的残差计算方法，得到新的目标$Y_{y}es$，此时不用和前向梯度算法一样小步小步的走，而是直接向前走直到出现一个$X_t$，使得$X_t$和$Y_{(yes)}$的相关度和$X_k$与$Y_{(yes)}$的相关度是一样的，此时残差$Y_{yes}$就在$X_t$和$X_k$的角平分线上，此时我们开始沿着这个残差角平分线走，直到出现第三个特征$X_p$和$Y_{y}es$的相关度等于${\theta }_t,{\theta }_k$与$Y_{yes}$的一样。将其也加入到$Y$的逼近特征集合中，并用$Y$的逼近特征集合的共同角分线，作为新的逼近方向，循环直到$Y_{yes}$足够小或者所有变量都取完位置。

当$\theta$只有2维时。例子如上图，和$Y$最接近的是$X_1$,首先在$X_1$上走一段距离，直到残差在$X_1$和$X_2$的角平分线上，此时沿着角平分线走，直到残差足够小才停止。此时对应的系数$\beta$即为最终的结果。

最小角回归法是一个适用于高维数据的回归算法，其主要的优点有：

1）特别适合于特征维度n 远高于样本数m的情况。

2）算法的最坏计算复杂度和最小二乘法类似，但是其计算速度几乎和前向选择算法一样

3）可以产生分段线性结果的完整路径，这在模型的交叉验证中极为有用

主要的缺点是：

由于LARS的迭代方向是根据目标的残差而定，所以该算法对样本的噪声极为敏感

Python实现