机器学习算法原理实现——使用梯度下降求解Lasso回归和岭回归-526互联

本文本质上是在线性回归的基础上进行扩展，加入了正则化而已！

正则化在机器学习中是一种防止过拟合的技术，它通过在损失函数中添加一个惩罚项来限制模型的复杂度。
举一个实际的例子，假设你正在训练一个机器学习模型来预测房价。你有很多特征，如房间数量、地理位置、建筑年份等。如果你的模型过于复杂，例如它尝试拟合每一个训练样本的细微差异，那么它可能在训练数据上表现得很好，但在新的、未见过的数据上表现得很差。这就是过拟合。

为了防止过拟合，你可以使用正则化。在这个例子中，正则化可能会通过对模型的权重施加某种惩罚（例如，使权重的平方和最小）来限制模型的复杂度。这样，模型就不能过于依赖任何一个特征，而是需要考虑所有的特征。这可以帮助模型在新的数据上表现得更好，因为它不会过于依赖训练数据中的特定模式，这些模式可能在新的数据中并不存在。

机器学习中，正则化通常通过在损失函数中添加一个惩罚项来实现。这个惩罚项通常与模型的权重有关。常见的正则化方法有L1正则化和L2正则化。

1. L1正则化（Lasso回归）：在损失函数中添加权重的绝对值的和。公式如下：

L = ∑(y - f(x))^2 + λ∑|w|

其中，y 是真实值，f(x) 是预测值，w 是模型的权重，λ 是正则化参数。

2. L2正则化（岭回归）：在损失函数中添加权重的平方和。公式如下：

L = ∑(y - f(x))^2 + λ∑w^2

其中，y 是真实值，f(x) 是预测值，w 是模型的权重，λ 是正则化参数。

这两种正则化方法都可以有效地防止模型过拟合，但它们的效果和适用场景可能会有所不同。L1正则化可以产生稀疏的权重，即许多权重为零，这可以用于特征选择。L2正则化则会使权重接近零，但不会完全为零，这可以防止权重过大。

如何使用梯度下降法求解呢？

$\text{[math]}$

对分量 $\text{[math]}$ 求偏导可得

$\text{[math]}$

Lasso回归是一种使用L1正则化的线性回归方法，它可以通过梯度下降法进行求解。下面是Lasso回归梯度下降法的求解步骤：

1. 初始化模型参数：初始化权重系数w和偏置b为0或者随机值。

2. 计算预测值：根据当前的权重系数w和偏置b，计算预测值y_hat。

3. 计算损失函数的梯度：计算损失函数关于预测值y_hat的梯度，即对于线性回归问题，计算损失函数关于y_hat的偏导数。

4. 计算L1正则化项的梯度：计算L1正则化项关于权重系数w的梯度，即L1范数的导数。L1范数的导数可以表示为w的符号函数。（高中数学，就是上图里的sign函数）

5. 计算总的梯度：将损失函数的梯度和L1正则化项的梯度进行组合，得到总的梯度。

6. 更新模型参数：根据梯度下降法的更新规则，更新权重系数w和偏置b。例如，w = w - learning_rate gradient_w，b = b - learning_rate gradient_b，其中learning_rate是学习率。

7. 重复步骤2-6，直到达到指定的迭代次数或收敛条件。

总结起来，Lasso回归的梯度下降法求解步骤是：初始化模型参数，计算预测值，计算损失函数的梯度，计算L1正则化项的梯度，计算总的梯度，更新模型参数。通过迭代这些步骤，可以逐步优化模型参数，实现Lasso回归的求解。

好了，可以开始代码编程了！

import numpy as np
 
 
 ### 定义符号函数
def sign(x):
    '''
    输入：
    x：浮点数值
    输出：
    整型数值
    '''
    if x > 0:
        return 1
    elif x < 0:
        return -1
    else:
        return 0


# 利用NumPy对符号函数进行向量化
vec_sign = np.vectorize(sign)
"""
向量化的符号函数 vec_sign 可以对NumPy数组进行元素级别的符号函数计算，
即对数组中的每个元素应用符号函数。这样可以提高计算效率，避免使用循环逐个计算。

import numpy as np

# 定义一个示例数组
arr = np.array([-2.5, 0, 3.7, -1.2, 0.5])

# 使用向量化的符号函数对数组进行计算
result = vec_sign(arr)
print(result)
输出：[-1  0  1 -1  1]
"""
 

### 包括线性回归模型公式、均方损失函数和参数求偏导三部分
def lasso_loss(X, y, w, b, alpha):
    '''
    输入：
    X：输入变量矩阵
    y：输出标签向量
    w：变量参数权重矩阵
    b：偏置
    alpha: L1正则
    输出：
    y_hat：线性回归模型预测值
    loss：均方损失
    dw：权重系数一阶偏导
    db：偏置一阶偏导
    '''
    # 训练样本量
    num_train = X.shape[0]
    # 训练特征数
    num_feature = X.shape[1]
    # print("num_train:", num_train, "num_feature:", num_feature)
    # 线性回归预测值
    y_hat = np.dot(X, w) + b
    # 计算预测值与实际标签之间的均方损失
    loss = np.sum((y_hat-y)**2) / num_train + alpha * np.sum(np.abs(w))
    # 基于均方损失对权重系数的一阶梯度
    dw = np.dot(X.T, (y_hat-y)) / num_train + alpha * vec_sign(w)
    # 基于均方损失对偏置的一阶梯度
    db = np.sum((y_hat-y)) / num_train
    return y_hat, loss, dw, db
 
 
### 初始化模型参数
def initialize_params(dims):
    '''
    输入：
    dims：训练数据的变量维度
    输出：
    w：初始化权重系数
    b：初始化偏置参数
    '''
    # 初始化权重系数为零向量
    w = np.zeros((dims, 1))
    # 初始化偏置参数为零
    b = 0
    return w, b
 
 
### 定义线性回归模型的训练过程
def lasso_train(X, y, learning_rate=0.01, epochs=1000, alpha=0.1):
    '''
    输入：
    X：输入变量矩阵
    y：输出标签向量
    learning_rate：学习率
    epochs：训练迭代次数
    输出：
    loss_his：每次迭代的均方损失
    params：优化后的参数字典
    grads：优化后的参数梯度字典
    '''
    # 记录训练损失的空列表
    loss_his = []
    # 初始化模型参数
    w, b = initialize_params(X.shape[1])
    # 迭代训练
    for i in range(1, epochs):
        # 计算当前迭代的预测值、均方损失和梯度
        y_hat, loss, dw, db = lasso_loss(X, y, w, b, alpha)
        # 基于梯度下降法的参数更新
        w += -learning_rate * dw
        b += -learning_rate * db
        # 记录当前迭代的损失
        loss_his.append(loss)
        # 每500次迭代打印当前损失信息
        if i % 500 == 0:
            print('epoch %d loss %f' % (i, loss))
        # 将当前迭代步优化后的参数保存到字典中
        params = {
            'w': w,
            'b': b
        }
        # 将当前迭代步的梯度保存到字典中
        grads = {
            'dw': dw,
            'db': db
        }
    return loss_his, params, grads
 
 
### 定义线性回归模型的预测函数
def predict(X, params):
    '''
    输入：
    X：测试集
    params：模型训练参数
    输出：
    y_pred：模型预测结果
    '''
    # 获取模型参数
    w = params['w']
    b = params['b']
    # 预测
    y_pred = np.dot(X, w) + b
    return y_pred



from sklearn.datasets import make_regression 
from sklearn.model_selection import train_test_split
from sklearn.linear_model import Lasso
from sklearn.metrics import mean_squared_error

# 生成样本数据
X, y = make_regression(n_samples=100, n_features=10, random_state=0, noise=0.5)
# 划分数据集
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=0)
X_train = np.array(X_train)
X_test = np.array(X_test)
y_train = np.array(y_train).reshape((-1, 1))
y_test = np.array(y_test).reshape((-1, 1))


# 创建Lasso回归模型实例
lasso = Lasso(alpha=0.1)
# 在训练集上拟合模型
lasso.fit(X_train, y_train)
# 在测试集上进行预测
y_pred = lasso.predict(X_test)
# 计算均方误差
mse = mean_squared_error(y_test, y_pred)
print("均方误差：", mse)
# 打印模型的截距和系数
print("截距：", lasso.intercept_)
print("系数：", lasso.coef_)


######################################
# 执行训练示例
loss, params, grads = lasso_train(X_train, y_train, 0.01, 5000, 0.1)
# 获取训练参数
print(params)
y_pred2 = predict(X_test, params)
mse2 = mean_squared_error(y_test, y_pred2)
print("均方误差：", mse2)

输出：
均方误差： 0.2847177370581323 ==》sklearn的输出
截距： [-0.04487816]
系数： [77.33057691 34.0732081 70.06993493 1.36023431 82.10576216 96.62025529
88.44587752 61.26161977 99.22830202 3.55673871]
epoch 500 loss 70.173990
epoch 1000 loss 61.768456
epoch 1500 loss 61.693305
epoch 2000 loss 61.690459
epoch 2500 loss 61.690221
epoch 3000 loss 61.690197
epoch 3500 loss 61.690195
epoch 4000 loss 61.690195
epoch 4500 loss 61.690195
{'w': array([[77.33094424],
       [34.07043903],
       [70.07046341],
       [ 1.35861401],
       [82.1045371 ],
       [96.61784068],
       [88.44665202],
       [61.26265291],
       [99.22846277],
       [ 3.55626159]]), 'b': -0.04567962607950337}
均方误差： 0.2861090798694742 ==》自己手写代码的输出

可以看到，我们自己写的和sklearn的输出已经很接近了！

岭回归的代码：

import numpy as np
 

### 定义LASSO回归损失函数
def ridge_loss(X, y, w, b, alpha):
    '''
    输入：
    X：输入变量矩阵
    y：输出标签向量
    w：变量参数权重矩阵
    b：偏置
    alpha：正则化系数
    输出：
    y_hat：线性模型预测输出
    loss：均方损失值
    dw：权重系数一阶偏导
    db：偏置一阶偏导
    '''
    # 训练样本量
    num_train = X.shape[0]
    # 训练特征数
    num_feature = X.shape[1]
    # 回归模型预测输出
    y_hat = np.dot(X, w) + b
    # L2损失函数
    loss = np.sum((y_hat-y)**2)/(2*num_train) + alpha*(np.sum(np.square(w)))/2 
    # 参数梯度计算
    dw = np.dot(X.T, (y_hat-y)) /num_train + alpha*w/num_train #######注意：这里有/num_train#######
    db = np.sum((y_hat-y)) /num_train
    return y_hat, loss, dw, db 

 
### 初始化模型参数
def initialize_params(dims):
    '''
    输入：
    dims：训练数据的变量维度
    输出：
    w：初始化权重系数
    b：初始化偏置参数
    '''
    # 初始化权重系数为零向量
    w = np.zeros((dims, 1))
    # 初始化偏置参数为零
    b = 0
    return w, b
 
 
### 定义线性回归模型的训练过程
def ridge_train(X, y, learning_rate=0.01, epochs=1000, alpha=0.1):
    '''
    输入：
    X：输入变量矩阵
    y：输出标签向量
    learning_rate：学习率
    epochs：训练迭代次数
    输出：
    loss_his：每次迭代的均方损失
    params：优化后的参数字典
    grads：优化后的参数梯度字典
    '''
    # 记录训练损失的空列表
    loss_his = []
    # 初始化模型参数
    w, b = initialize_params(X.shape[1])
    # 迭代训练
    for i in range(1, epochs):
        # 计算当前迭代的预测值、均方损失和梯度
        y_hat, loss, dw, db = ridge_loss(X, y, w, b, alpha)
        # 基于梯度下降法的参数更新
        w += -learning_rate * dw
        b += -learning_rate * db
        # 记录当前迭代的损失
        loss_his.append(loss)
        # 每500次迭代打印当前损失信息
        if i % 5000 == 0:
            print('epoch %d loss %f' % (i, loss))
        # 将当前迭代步优化后的参数保存到字典中
        params = {
            'w': w,
            'b': b
        }
        # 将当前迭代步的梯度保存到字典中
        grads = {
            'dw': dw,
            'db': db
        }
    return loss_his, params, grads
 
 
### 定义线性回归模型的预测函数
def predict(X, params):
    '''
    输入：
    X：测试集
    params：模型训练参数
    输出：
    y_pred：模型预测结果
    '''
    # 获取模型参数
    w = params['w']
    b = params['b']
    # 预测
    y_pred = np.dot(X, w) + b
    return y_pred



from sklearn.datasets import make_regression 
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error
from sklearn.linear_model import Ridge
from sklearn.preprocessing import StandardScaler

# 生成样本数据
X, y = make_regression(n_samples=100, n_features=10, random_state=0, noise=0.5)
# 归一化数据
scaler = StandardScaler()
X = scaler.fit_transform(X)

# 划分数据集
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=0)
X_train = np.array(X_train)
X_test = np.array(X_test)
y_train = np.array(y_train).reshape((-1, 1))
y_test = np.array(y_test).reshape((-1, 1))

# 创建Ridge回归模型实例
ridge = Ridge(alpha=1.0)
ridge.fit(X_train, y_train)
# 在测试集上进行预测
y_pred = ridge.predict(X_test)
# 计算均方误差
mse = mean_squared_error(y_test, y_pred)
print("均方误差：", mse)
# 打印模型的截距和系数
print("截距：", ridge.intercept_)
print("系数：", ridge.coef_)


######################################
# 执行训练示例
loss, params, grads = ridge_train(X_train, y_train, 0.01, 50000, 1)
# 获取训练参数
print(params)
y_pred2 = predict(X_test, params)
mse2 = mean_squared_error(y_test, y_pred2)
print("均方误差：", mse2)

结果：

均方误差： 7.295579168041554

截距： [-34.3629433]

系数： [[ 75.16219935 30.20047069 66.95939498 1.37559128 81.28261278

98.67605116 79.6055435 52.40872579 100.44228584 3.6482418 ]]

epoch 5000 loss 23291.978297

epoch 10000 loss 23291.978297

epoch 15000 loss 23291.978297

epoch 20000 loss 23291.978297

epoch 25000 loss 23291.978297

epoch 30000 loss 23291.978297

epoch 35000 loss 23291.978297

epoch 40000 loss 23291.978297

epoch 45000 loss 23291.978297

{'w': array([[ 75.16219935],

[ 30.20047069],

[ 66.95939498],

[ 1.37559128],

[ 81.28261278],

[ 98.67605116],

[ 79.6055435 ],

[ 52.40872579],

[100.44228584],

[ 3.6482418 ]]), 'b': -34.362943295912025}

均方误差： 7.295579168052714 ==》和sklearn结果基本一致！！！

注意：

# 归一化数据是可选的，去掉了结果也没多大影响

# scaler = StandardScaler()

# X = scaler.fit_transform(X) # 可选

影响最大当属参数梯度计算：

    # L2损失函数
    loss = np.sum((y_hat-y)**2)/(2*num_train) + alpha*(np.sum(np.square(w)))/2 
    # 参数梯度计算
    dw = np.dot(X.T, (y_hat-y)) /num_train + alpha*w/num_train #######注意：这里有/num_train#######
如果是没有/num_train的话，则结果会有很大差异！为啥？？？网上没有找到相关解释。

GPT4给出解释如下：