Scikit-Learn线性回归(六)

1、Lasso回归（套索回归）

本文接上篇：Scikit-Learn线性回归(五)

在上篇中，我们详解了误差与模型复杂度的关系、范数与正则化、Scikit-Learn岭回归及使用案例

我们已经知道，岭回归模型可以解决线性回归模型中矩阵$X^{\prime }X$不可逆的问题，岭回归通过给损失函数添加L2正则化项（惩罚项）使得偏回归系数缩小（趋于0）。但是不管如何减小，岭回归都会始终保留建模时的所有特征（自变量），因此岭回归无法降低模型的复杂度。从广义上来讲，岭回归只是让回归曲线变得更平滑（收敛）了，这一点上，岭回归也算降低了模型的复杂度。岭回归最大的问题在于无法稀疏化特征、不能做特征选择。而Lasso回归模型克服了一点

岭回归为损失函数添加了L2正则化项，而Lasso回归为损失函数添加了L1正则化项。Lasso回归的损失函数定义为：
$$L(\omega )=\sum _{i=1}^m(y_i-y{)}^2+\lambda \sum _{i=1}^m|{\omega }_i|=||Y-X\omega |{|}^2+\lambda ||\omega |{|}_1$$

其中，$\lambda ||\omega |{|}_1=\lambda {\sum }_{i=1}^m|{\omega }_i|$，即L1正则化项，它是目标函数的惩罚项。需要注意的是，正则化项中的回归系数${\omega }_i$为每个自变量对应的回归系数，不包含回归常数项。复杂度参数（惩罚项系数）λ用于控制收缩量

Lasso回归通过使用L1正则化项，它能够在回归系数中引入稀疏性，也就是允许某些系数在优化过程中缩小为零，从而实现特征的选择

与岭回归模型类似，Lasso回归通过添加L1正则化项也可以解决多重共线性（矩阵$X^{\prime }X$不可逆）问题，Lasso回归同样属于缩减性估计，它也使用了缩减（Shrinkage）回归系数的方法。但不同的是，Lasso回归针对不同的自变量，会使其收敛的速度不一样。有的变量很快就趋于0了，有的却会很慢。因此Lasso回归非常适合于做特征选择

Lasso（Least Absolute Shrinkage and Selection Operator）全称为最小绝对收缩和选择算子，Lasso回归（套索回归）是一种线性回归模型的扩展，主要用于解决高维数据中的特征选择问题

2、Lasso回归的解

模型训练的过程就是寻找让损失函数$L(\omega )$最小时的参数解ω。Lasso回归为损失函数添加了一个带惩罚系数λ的ω向量的L1范数作为惩罚项，所以这种正则化方式称为L1正则化。由于惩罚项是关于回归系数ω的绝对值之和，因此惩罚项在零点处是不可导的，导致其代价函数不是处处可导的，所以就没办法通过直接求导的方式来直接得到解析解ω，类似于应用在岭回归上的最小二乘法在此失效。不仅如此，梯度下降法（详见：传送门）也无法搜索出Lasso回归的最优解ω

函数可导的充要条件：函数$f(x)$在$x$点连续且左导数、右导数都存在并相等

为了能够得到Lasso回归的系数ω，可以使用坐标轴下降法（Coordinate Descent）和最小角回归法（Least Angle Regression，LARS）

2.1、坐标轴下降法

坐标轴下降法与梯度下降法类似，都属于迭代算法，不同的是，坐标轴下降法是是沿着坐标轴（维度）下降，而梯度下降法是沿着梯度的负方向下降。坐标下降法的数学思想是：对于m维参数的可微凸函数$J(\omega )$，如果存在一点$\hat{\omega }$，使得函数$J(\omega )$在每个坐标轴上均达到最小值，则$J(\hat{\omega })$就是点$\hat{\omega }$上的全局最小值

坐标轴下降法的原理及Lasso回归解的推导过程参考：https://segmentfault.com/a/1190000041240490

Lasso回归坐标轴下降法算法实现（参考上述链接）如下：

import numpy as np

# Lasso回归，使用坐标下降法求解
# 论文：https://www.cs.cmu.edu/~ggordon/10725-F12/slides/25-coord-desc.pdf
# X：训练数据集 y：目标标签值 alpha：惩罚项系数 n_iter：最大迭代次数 tol：变化量容忍值
# 返回权重系数ω
def coordinate_descent(X, y, alpha=0.1, n_iter=1e4, tol=1e-8):
    # 初始化w为零向量
    w = np.zeros(X.shape[1])
    for i_iter in range(int(n_iter)):
        done = True
        # 遍历所有自变量
        for i in range(0, len(w)):
            # 记录上一轮系数
            weight = w[i]
            # 求出当前条件下的最佳系数
            w[i] = down(X, y, w, i, alpha)
            # 当其中一个系数变化量未到达其容忍值，继续循环
            if np.abs(weight - w[i]) > tol:
                done = False
        # 所有系数都变化不大时，结束循环
        if done: break
    return w

# 坐标下降迭代
def down(X, y, w, index, alpha=0.1):
    """
    cost(w) = (x1 * w1 + x2 * w2 + ... - y)^2 + ... + λ(|w1| + |w2| + ...)
    为简化计算，假设只有w1一个变量，这时其他的值均为常数，带入上式得
    cost(w1) = (a * w1 + b)^2 + ... + λ|w1| + c (a,b,c,λ均为常数)
    展开后
    cost(w1) = aa * w1^2 + 2ab * w1 + λ|w1| + c (aa,ab,c,λ均为常数)
    """
    # 展开后的二次项的系数之和
    aa = 0
    # 展开后的一次项的系数之和
    ab = 0
    for i in range(X.shape[0]):
        # 括号内一次项的系数
        a = X[i][index]
        # 括号内常数项的系数
        b = X[i][:].dot(w) - a * w[index] - y[i]
        # 可以很容易的得到展开后的二次项的系数为括号内一次项的系数平方的和
        aa = aa + a * a
        # 可以很容易的得到展开后的一次项的系数为括号内一次项的系数乘以括号内常数项的和
        ab = ab + a * b
    # 由于是一元二次函数，当导数为零时，函数值最小值，只需要关注二次项系数、一次项系数和λ
    return det(aa, ab, alpha)

# 代价函数求导
def det(aa, ab, alpha=0.1):
    """
    通过代价函数的导数求w，当w=0时，不可导
    det(w) = NaN (w = 0)
    => w = 0
    det(w) = 2aa * w + 2ab + λ = 0 (w > 0)
    => w = - (2 * ab + λ) / (2 * aa)
    det(w) = 2aa * w + 2ab - λ = 0 (w < 0)
    => w = - (2 * ab - λ) / (2 * aa)
    """
    w = - (2 * ab + alpha) / (2 * aa)
    if w < 0:
        w = - (2 * ab - alpha) / (2 * aa)
        if w > 0:
            w = 0