手撸机器学习算法 - 岭回归

系列文章目录：

• 感知机
• 线性回归
• 非线性问题
• 多项式回归
• 岭回归

算法介绍

今天我们来一起学习一个除了线性回归、多项式回归外最最最简单的回归算法：岭回归，如果用等式来介绍岭回归，那么就是：\(岭回归 = 多项式回归 + 惩罚项\)，\(多项式回归 = 线性回归 + 多项式特征构建\)，从上述等式可以看到，所谓学习岭回归，只需要学习多项式和惩罚项即可，由于之前我们已经学习过多项式回归了，因此现在的重点是惩罚项或者叫正则项；

从多项式回归到岭回归

岭回归是在多项式回归的基础上增加了惩罚项，准确的说法是：在多项式回归的优化函数上增加了约束条件用于限制算法的假设空间，以应对模型的过拟合问题，下面我们分别看看如何增加约束条件、为什么可以防止过拟合、约束条件对推导的影响；

算法推导

既然岭回归是在多项式回归的基础上实现的，那么我们就以一个二元二次多项式回归为例子：
\(w_0*x_0^2+w_1*x_1^2+w_2*x_0*x_1+w_3*x_0+w_4*x_1+b\)
假设现在通过上述模型拟合数据显示过拟合，一般的做法是将模型从二阶降低到一阶（降阶可以减少特征数），则模型变为：
\(w_3*x_0+w_4*x_1+b\)
这个降阶的方式可以为手动指定w0、w1、w2为0来实现，对于多项式回归来说，它唯一控制模型复杂度的就是阶数，阶数越大，特征越多，模型越复杂，反之则越简单，但是这种控制方法难免显得不够灵活平滑，如果我们期望更平滑的降低复杂度的方法呢，这时就需要通过惩罚项来实现；

如何增加约束条件

增加约束的方式也很简单，从公式上看就是增加了服从条件，如下对条件W增加的约束，使得W的可取范围为半径为r的圆内：

\[y = w_0*x + w_1*b \\ s.t. ||W||^2 < r^2, W=(w_0, w_1) \]

为什么可以防止过拟合

对于上述约束，我们可以这样理解它，在没有加约束之前W=(w0 w1)的所有取值为整个二维平面上的点，而\(||W||^2 < r^2\)将W限制在原点为中心，半径为r的圆内，由于它减少了W的可取范围，因此起到了降低算法的假设空间（或者说是算法复杂度）的效果，也就可以作为一个有效的惩罚项；

约束条件下的公式推导

首先我们回顾下线性回归的公式推导，首先优化目标如下：

\[argmin \frac{1}{N}\sum_{i=1}^{N}(w*x_i+b-y_i)^2 \]

在岭回归中，优化目标增加了约束条件，如下：

\[argmin \frac{1}{N}\sum_{i=1}^{N}(w*x_i+b-y_i)^2 \\ s.t. ||W||^2 < r^2, W=(w_0, w_1) \]

通过拉格朗日将约束条件转为函数的一部分，通过添加拉格朗日乘子，注意下述公式的λ作为超参数，因此看作常量，如下：

\[argmin \frac{1}{N}\sum_{i=1}^{N}(Wx_i-y_i)^2 + λ(||W||^2-r^2) \]

对上述公式针对W求导并令其为零有：

\[\frac{1}{N}N(2x_i^Tx_iW-2x_i^Ty_i+0)+2λW = 0 \\ 2x_i^Tx_iW + 2λW = 2x_i^Ty_i \\ (x_i^Tx_i+λI)W = x_i^Ty_i, I为单位阵 \\ W = (x_i^Tx_i+λI)^{-1}x_i^Ty_i \]

代码实现

岭回归对象初始化

可以看到，当岭回归的拉格朗日乘子λ为0时，岭回归退化为多项式回归：

if self.lambdaVal == 0:
    return super(RidgeRegression,self).train()

参数W计算

xTx = self.X.T @ self.X
I = np.eye(xTx.shape[0])
self.w = np.linalg.inv(xTx + self.lambdaVal*I) @ self.X.T @ self.y
self.w = self.w.reshape(-1)
self.w,self.b = self.w[1:],self.w[0]

运行结果

下面使用5阶多项式回归来观察看使用惩罚项与不使用的区别，可以很容易的看到，由于惩罚项参数λ的存在，使得同样为5阶的模型，惩罚系数越大，模型越趋于简单；

不使用惩罚项

分别使用λ=0.1、1、10作为惩罚系数

全部代码

import numpy as np
import matplotlib.pyplot as plt
from 线性回归最小二乘法矩阵实现 import LinearRegression as LR
from 多项式回归 import PolynomialRegression as PR

'''
惩罚项：亦称为罚项、正则项，用于限制模型复杂度，在公式上可看到是增加了某个约束条件，即：subject to xxxx；

以多项式回归理解惩罚项：
对于二元二次多项式回归，所有假设空间可能为：w0*x0^2+w1*x1^2+w2*x0*x1+w3*x0+w4*x1+b
当二阶多项式过拟合时，通常考虑退回到一阶，即线性回归，假设空间为：w0*x0+w1*x1+b
这种退化可以看到是对二阶多项式增加了约束条件：w0=0,w1=0,w2=0
因此对于多项式回归，任意低阶都可以看作是其高阶+惩罚项的组合结果

惩罚项的意义：通过对公式增加灵活的约束条件，可以更平滑的控制模型复杂度，只要约束条件是有意义的，那么它就降低了原假设空间的大小，例如对于线性回归w0*x0+b，W=(w0 w1)，即W的可取范围为整个二维平面，如果增加约束条件w0^2+w1^2<r^2，则W的取值范围为二维平面上以r为半径的圆内，而W决定了线性回归的假设空间大小，因此通过约束条件得以降低假设空间大小的目的；

岭回归 = 线性回归 + 优化目标(argmin MSE)上增加约束条件(s.t. ||w||^2<=r^2)
'''

class RidgeRegression(PR):
    def __init__(self,X,y,degrees=1,lambdaVal=0):
        super(RidgeRegression,self).__init__(X,y,degrees)
        self.lambdaVal = lambdaVal

    def train(self):
        if self.lambdaVal == 0:
            return super(RidgeRegression,self).train()
        xTx = self.X.T @ self.X
        I = np.eye(xTx.shape[0])
        self.w = np.linalg.inv(xTx + self.lambdaVal*I) @ self.X.T @ self.y
        self.w = self.w.reshape(-1)
        self.w,self.b = self.w[1:],self.w[0]
        return self.w,self.b

def pain(pos=141,xlabel='x',ylabel='y',title='',x=[],y=[],line_x=[],line_y=[]):
    plt.subplot(pos)
    plt.title(title)
    plt.xlabel(xlabel)
    plt.ylabel(ylabel)
    plt.scatter(x,y)
    plt.plot(line_x,line_y)

if __name__ == '__main__':
    rnd = np.random.RandomState(3)
    x_min, x_max = 0, 10
    
    def pain(pos=141,xlabel='x',ylabel='y',title='',x=[],y=[],line_x=[],line_y=[]):
        plt.subplot(pos)
        plt.title(title)
        plt.xlabel(xlabel)
        plt.ylabel(ylabel)
        plt.scatter(x,y)
        plt.plot(line_x,line_y)
    
    # 上帝函数 y=f(x)
    def f(x):
        return x**5-22*x**4+161*x**3-403*x**2+36*x+938
    
    # 上帝分布 P(Y|X)
    def P(X):
        return f(X) + rnd.normal(scale=30, size=X.shape)
    
    # 通过 P(X, Y) 生成数据集 D
    X = rnd.uniform(x_min, x_max, 10)   # 通过均匀分布产生 X
    y = P(X)                            # 通过 P(Y|X) 产生 y

    X,y = X.reshape(-1,1),y.reshape(-1,1)
    x_min,x_max = min(X),max(X)


    for pos,deg in zip([331,332,333],[2,5,10]):
        model = PR(X=X,y=y,degrees=deg)
        w,b = model.train()
        print(f'最小二乘法的矩阵方式结果为：w={w} b={b}')
        line_x = [x_min+(x_max-x_min)*(i/100) for i in range(-1,102,1)]
        line_y = [model.predict(x) for x in line_x]
        pain(pos,'X','y','DEG='+str(deg),X[:,0],y[:,0],line_x,line_y)
    for pos,deg,lambdaVal in zip([334,335,336],[5,5,5],[0.1,1,10]):
        model = RidgeRegression(X=X,y=y,degrees=deg,lambdaVal=lambdaVal)
        w,b = model.train()
        print(f'最小二乘法的矩阵方式结果为：w={w} b={b}')
        line_x = [x_min+(x_max-x_min)*(i/100) for i in range(-1,102,1)]
        line_y = [model.predict(x) for x in line_x]
        pain(pos,'X','y','DEG='+str(deg)+', λ='+str(lambdaVal),X[:,0],y[:,0],line_x,line_y)
    for pos,deg,lambdaVal in zip([337,338,339],[10,10,10],[0.1,1,10]):
        model = RidgeRegression(X=X,y=y,degrees=deg,lambdaVal=lambdaVal)
        w,b = model.train()
        print(f'最小二乘法的矩阵方式结果为：w={w} b={b}')
        line_x = [x_min+(x_max-x_min)*(i/100) for i in range(-1,102,1)]
        line_y = [model.predict(x) for x in line_x]
        pain(pos,'X','y','DEG='+str(deg)+', λ='+str(lambdaVal),X[:,0],y[:,0],line_x,line_y)
    
    plt.show()

最后

相对于多项式回归，由于岭回归有惩罚项的存在，因此它可以更加肆无忌惮的使用更高阶的多项式而不需要太担心过拟合的问题，因此理论上多项式回归能做到的，岭回归可以做的更好，当然了由于参数λ的存在，岭回归需要调的参数也更多了；