机器学习算法实践-岭回归和LASSO

前言

继续线性回归的总结, 本文主要介绍两种线性回归的缩减(shrinkage)方法的基础知识: 岭回归(Ridge Regression)和LASSO(Least Absolute Shrinkage and Selection Operator)并对其进行了Python实现。同时也对一种更为简单的向前逐步回归计算回归系数的方法进行了相应的实现。

正文

通过上一篇《机器学习算法实践-标准与局部加权线性回归》中标准线性回归的公式$w = (X^{T}X)^{-1}X^{T}y$中可以看出在计算回归系数的时候我们需要计算矩阵$X^TX$的逆，但是如果该矩阵是个奇异矩阵，则无法对其进行求解。那么什么情况下该矩阵会有奇异性呢?

X本身存在线性相关关系(多重共线性), 即非满秩矩阵。如果数据的特征中存在两个相关的变量，即使并不是完全线性相关，但是也会造成矩阵求逆的时候造成求解不稳定。
当数据特征比数据量还要多的时候, 即$m < n$, 这时候矩阵$X$是一个矮胖型的矩阵，非满秩。

对于上面的两种情况，我们需要对最初的标准线性回归做一定的变化使原先无法求逆的矩阵变得非奇异，使得问题可以稳定求解。我们可以通过缩减的方式来处理这些问题例如岭回归和LASSO.

中心化和标准化

这里先介绍下数据的中心化和标准化，在回归问题和一些机器学习算法中通常要对原始数据进行中心化和标准化处理，也就是需要将数据的均值调整到0，标准差调整为1, 计算过程很简单就是将所有数据减去平均值后再除以标准差:
$$
x_i^{‘} = \frac{x_i - \mu}{\sigma}
$$

这样调整后的均值:
$$\mu^{‘} = \frac{(\sum_{i=1}^{n}x_i)/n - \mu}{\sigma} = 0$$

调整后的标准差:
$$
\sigma^{‘} = \frac{(x_i - \mu)^2/n}{\sigma^2} = \frac{\sigma^2}{\sigma^2} = 1
$$

之所以需要进行中心化其实就是个平移过程，将所有数据的中心平移到原点。而标准化则是使得所有数据的不同特征都有相同的尺度Scale, 这样在使用梯度下降法以及其他方法优化的时候不同特征参数的影响程度就会一致了。

如下图所示，可以看出得到的标准化数据在每个维度上的尺度是一致的(图片来自网络，侵删)

岭回归(Ridge Regression)

标准最小二乘法优化问题:
$$
f(w) = \sum_{i=1}^{m} (y_i - x_{i}^{T}w)^2
$$
也可以通过矩阵表示:
$$
f(w) = (y - Xw)^{T}(y - Xw)
$$
得到的回归系数为:
$$
\hat{w} = (X^{T}X)^{-1}X^{T}y
$$

这个问题解存在且唯一的条件就是$X$列满秩:$rank(X) = dim(X)$.

即使$X$列满秩，但是当数据特征中存在共线性，即相关性比较大的时候，会使得标准最小二乘求解不稳定, $X^TX$的行列式接近零，计算$X^TX$的时候误差会很大。这个时候我们需要在cost function上添加一个惩罚项$\lambda\sum_{i=1}^{n}w_{i}^2$，称为L2正则化。

这个时候的cost function的形式就为:
$$
f(w) = \sum_{i=1}^{m} (y_i - x_{i}^{T}w)^2 + \lambda\sum_{i=1}^{n}w_{i}^{2}
$$

通过加入此惩罚项进行优化后，限制了回归系数$w_i$的绝对值，数学上可以证明上式的等价形式如下:
$$
f(w) = \sum_{i=1}^{m} (y_i - x_{i}^{T}w)^2 \\
s.t. \sum_{i=1}^{n}w_{j}^2 \le t
$$
其中t为某个阈值。

将岭回归系数用矩阵的形式表示:
$$
\hat{w} = (X^{T}X + \lambda I)^{-1}X^{T}y
$$

可以看到，就是通过将$X^TX$加上一个单位矩阵是的矩阵变成非奇异矩阵并可以进行求逆运算。

岭回归的几何意义

以两个变量为例, 残差平方和可以表示为$w_1, w_2$的一个二次函数，是一个在三维空间中的抛物面，可以用等值线来表示。而限制条件$w_1^2 + w_2^2 < t$，相当于在二维平面的一个圆。这个时候等值线与圆相切的点便是在约束条件下的最优点，如下图所示，

岭回归的一些性质

当岭参数$\lambda = 0$时，得到的解是最小二乘解
当岭参数$\lambda$趋向更大时，岭回归系数$w_i$趋向于0，约束项$t$很小

岭回归的Python实现

通过矩阵的形式计算$\hat{w}$, 可以很简单的实现

def ridge_regression(X, y, lambd=0.2):
    ''' 获取岭回归系数
    '''
    XTX = X.T*X
    m, _ = XTX.shape
    I = np.matrix(np.eye(m))
    w = (XTX + lambd*I).I*X.T*y
    return w

岭迹图

可以知道求得的岭系数$w_i$是岭参数$\lambda$的函数，不同的$\lambda$得到不同的岭参数$w_i$, 因此我们可以增大$\lambda$的值来得到岭回归系数的变化，以及岭参数的变化轨迹图(岭迹图), 不存在奇异性时，岭迹图应稳定的逐渐趋向于0。

通过岭迹图我们可以:

观察较佳的$\lambda$取值
观察变量是否有多重共线性

绘制岭迹图

上面我们通过函数ridge_regression实现了计算岭回归系数的计算，我们使用《机器学习实战》中的鲍鱼年龄的数据来进行计算并绘制不同$\lambda$的岭参数变化的轨迹图。数据以及完整代码详见 https://github.com/PytLab/MLBox/tree/master/linear_regression

选取30组不同的$\lambda$来获取岭系数矩阵包含30个不同的岭系数。

def ridge_traj(X, y, ntest=30):
    ''' 获取岭轨迹矩阵
    '''
    _, n = X.shape
    ws = np.zeros((ntest, n))
    for i in range(ntest):
        w = ridge_regression(X, y, lambd=exp(i-10))
        ws[i, :] = w.T
    return ws

绘制岭轨迹图

if '__main__' == __name__:
    ntest = 30
    # 加载数据
    X, y = load_data('abalone.txt')

    # 中心化 & 标准化
    X, y = standarize(X), standarize(y)

    # 绘制岭轨迹
    ws = ridge_traj(X, y, ntest)
    fig = plt.figure()
    ax = fig.add_subplot(111)

    lambdas = [i-10 for i in range(ntest)]
    ax.plot(lambdas, ws)

    plt.show()

上图绘制了回归系数$w_i$与$log(\lambda)$的关系，在最左边$\lambda$系数最小时，可以得到所有系数的原始值(与标准线性回归相同); 而在右边，系数全部缩减为0, 从不稳定趋于稳定；为了定量的找到最佳参数值，还需要进行交叉验证。要判断哪些变量对结果的预测最具影响力，可以观察他们的系数大小即可。

LASSO

岭回归限定了所有回归系数的平方和不大于$t$, 在使用普通最小二乘法回归的时候当两个变量具有相关性的时候，可能会使得其中一个系数是个很大正数，另一个系数是很大的负数。通过岭回归的$\sum_{i=1}^{n} w_i \le t$的限制，可以避免这个问题。

LASSO(The Least Absolute Shrinkage and Selection Operator)是另一种缩减方法，将回归系数收缩在一定的区域内。LASSO的主要思想是构造一个一阶惩罚函数获得一个精炼的模型, 通过最终确定一些变量的系数为0进行特征筛选。

LASSO的惩罚项为:
$$
\sum_{i=1}^{n} \vert w_i \vert \le t
$$

与岭回归的不同在于，此约束条件使用了绝对值的一阶惩罚函数代替了平方和的二阶函数。虽然只是形式稍有不同，但是得到的结果却又很大差别。在LASSO中，当$\lambda$很小的时候，一些系数会随着变为0而岭回归却很难使得某个系数恰好缩减为0. 我们可以通过几何解释看到LASSO与岭回归之间的不同。

LASSO的几何解释

同样以两个变量为例，标准线性回归的cost function还是可以用二维平面的等值线表示，而约束条件则与岭回归的圆不同，LASSO的约束条件可以用方形表示，如下图:

相比圆，方形的顶点更容易与抛物面相交，顶点就意味着对应的很多系数为0，而岭回归中的圆上的任意一点都很容易与抛物面相交很难得到正好等于0的系数。这也就意味着，lasso起到了很好的筛选变量的作用。

LASSO回归系数的计算

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

坐标下降法是每次选择一个维度的参数进行一维优化，然后不断的迭代对多个维度进行更新直到函数收敛。SVM对偶问题的优化算法SMO也是类似的原理，这部分的详细介绍我在之前的一篇博客中进行了整理，参考《机器学习算法实践-SVM中的SMO算法》。

下面我们分别对LASSO的cost function的两部分求解：

RSS部分

$$
RSS(w) = \sum_{i=1}^{m}(y_i - \sum_{j=1}^{n}x_{ij}w_j)^2
$$

求导:
$$
\frac{\partial RSS(w)}{\partial w_k} = -2\sum_{i=1}^{m}x_{ik}(y_i - \sum_{j=1}^{n}x_{ij}w_j) \\
= -2\sum_{i=1}^{m}(x_{ik}y_i - x_{ik}\sum_{j=1, j \ne k}^{n}x_{ij}w_j - x_{ik}^2w_k) \\
= -2\sum_{i=1}^{m}x_{ik}(y_i - \sum_{j=1, j \ne k}^{n}x_{ij}w_{j}) + 2w_k\sum_{i=1}^{m}x_{ik}^2
$$

令$p_k = \sum_{i=1}^{m}x_{ik}(y_i - \sum_{j=1, j \ne k}^{n}x_{ij}w_{j})$, $z_k = \sum_{i=1}{m}x_{ik}^2$ 得到:

$$
\frac{\partial RSS(w)}{\partial w_j} = -2p_k + 2z_kw_k
$$

正则项

关于惩罚项的求导我们需要使用subgradient，可以参考LASSO（least absolute shrinkage and selection operator）回归中如何用梯度下降法求解？

$$
\lambda \frac{\partial \sum_{i=1}^{n}\vert w_j \vert}{\partial w_k} = \begin{cases}
-\lambda & w_k < 0 \\
[-\lambda, \lambda] & w_k = 0 \\
\lambda & w_k > 0
\end{cases}
$$

这样整体的偏导数:
$$
\frac{\partial f(w)}{\partial w_k} = 2z_kw_k - 2p_k + \begin{cases}
-\lambda & w_k < 0 \\
[-\lambda, \lambda] & w_k = 0 \\
\lambda & w_k > 0
\end{cases} \\
= \begin{cases}
2z_kw_k - 2p_k - \lambda & w_k < 0 \\
[-2p_k - \lambda, -2p_k + \lambda] & w_j = 0 \\
2z_kw_k - 2p_k + \lambda & w_k > 0
\end{cases}
$$

令$\frac{\partial f(w)}{\partial w_k} = 0$ 得到

$$
\hat{w_k} = \begin{cases}
(p_k + \lambda/2)/z_k & p_k < -\lambda/2 \\
0 & -\lambda/2 \le p_k \le \lambda/2 \\
(p_k - \lambda/2)/z_k & p_k > \lambda/2
\end{cases}
$$

通过上面的公式我们便可以每次选取一维进行优化并不断迭代得到最优回归系数。

LASSO的Python实现

根据上面代码我们实现梯度下降法并使用其获取LASSO回归系数。

def lasso_regression(X, y, lambd=0.2, threshold=0.1):
    ''' 通过坐标下降(coordinate descent)法获取LASSO回归系数
    '''
    # 计算残差平方和
    rss = lambda X, y, w: (y - X*w).T*(y - X*w)

    # 初始化回归系数w.
    m, n = X.shape
    w = np.matrix(np.zeros((n, 1)))
    r = rss(X, y, w)

    # 使用坐标下降法优化回归系数w
    niter = itertools.count(1)

    for it in niter:
        for k in range(n):
            # 计算常量值z_k和p_k
            z_k = (X[:, k].T*X[:, k])[0, 0]
            p_k = 0
            for i in range(m):
                p_k += X[i, k]*(y[i, 0] - sum([X[i, j]*w[j, 0] for j in range(n) if j != k]))

            if p_k < -lambd/2:
                w_k = (p_k + lambd/2)/z_k
            elif p_k > lambd/2:
                w_k = (p_k - lambd/2)/z_k
            else:
                w_k = 0

            w[k, 0] = w_k

        r_prime = rss(X, y, w)
        delta = abs(r_prime - r)[0, 0]
        r = r_prime
        print('Iteration: {}, delta = {}'.format(it, delta))

        if delta < threshold:
            break

    return w

我们选取$\lambda = 10$, 收敛阈值为0.1来获取回归系数

if '__main__' == __name__:
    X, y = load_data('abalone.txt')
    X, y = standarize(X), standarize(y)
    w = lasso_regression(X, y, lambd=10)

    y_prime = X*w
    # 计算相关系数
    corrcoef = get_corrcoef(np.array(y.reshape(1, -1)),
                            np.array(y_prime.reshape(1, -1)))
    print('Correlation coefficient: {}'.format(corrcoef))

迭代了150步收敛到0.1，计算相对比较耗时:

Iteration: 146, delta = 0.1081124857935265
Iteration: 147, delta = 0.10565615985365184
Iteration: 148, delta = 0.10326058648411163
Iteration: 149, delta = 0.10092418256476776
Iteration: 150, delta = 0.09864540659987142
Correlation coefficient: 0.7255254877587117

LASSO回归系数轨迹

类似岭轨迹，我们也可以改变$\lambda$的值得到不同的回归系数，通过作图可以看到回归系数的轨迹

ntest = 30

# 绘制轨迹
ws = lasso_traj(X, y, ntest)
fig = plt.figure()
ax = fig.add_subplot(111)

lambdas = [i-10 for i in range(ntest)]
ax.plot(lambdas, ws)

plt.show()

得到的轨迹图如下:

通过与岭轨迹图进行对比发现，随着$\lambda$的增大，系数逐渐趋近于0，但是岭回归没有系数真正为0，而lasso中不断有系数变为0.迭代过程中输出如下图:

逐步向前回归

LASSO计算复杂度相对较高，本部分稍微介绍一下逐步向前回归，他属于一种贪心算法，给定初始系数向量，然后不断迭代遍历每个系数，增加或减小一个很小的数，看看代价函数是否变小，如果变小就保留，如果变大就舍弃，然后不断迭代直到回归系数达到稳定。

下面给出实现

def stagewise_regression(X, y, eps=0.01, niter=100):
    ''' 通过向前逐步回归获取回归系数
    '''
    m, n = X.shape
    w = np.matrix(np.zeros((n, 1)))
    min_error = float('inf')
    all_ws = np.matrix(np.zeros((niter, n)))

    # 计算残差平方和
    rss = lambda X, y, w: (y - X*w).T*(y - X*w)

    for i in range(niter):
        print('{}: w = {}'.format(i, w.T[0, :]))
        for j in range(n):
            for sign in [-1, 1]:
                w_test = w.copy()
                w_test[j, 0] += eps*sign
                test_error = rss(X, y, w_test)
                if test_error < min_error:
                    min_error = test_error
                    w = w_test
        all_ws[i, :] = w.T

    return all_ws

我们去变化量为0.005，迭代步数为1000次，得到回归系数随着迭代次数的变化曲线:

逐步回归算法的主要有点在于他可以帮助人们理解现有的模型并作出改进。当构建了一个模型后，可以运行逐步回归算法找出重要的特征，即使停止那些不重要特征的收集。

总结

本文介绍了两种回归中的缩减方法，岭回归和LASSO。两种回归均是在标准线性回归的基础上加上正则项来减小模型的方差。这里其实便涉及到了权衡偏差(Bias)和方差(Variance)的问题。方差针对的是模型之间的差异，即不同的训练数据得到模型的区别越大说明模型的方差越大。而偏差指的是模型预测值与样本数据之间的差异。所以为了在过拟合和欠拟合之前进行权衡，我们需要确定适当的模型复杂度来使得总误差最小。

参考

《Machine Learning in Action》
机器学习中的Bias(偏差)，Error(误差)，和Variance(方差)有什么区别和联系？
Lasso回归的坐标下降法推导
数据什么时候需要做中心化和标准化处理？