手撕机器学习算法 | 线性回归 | LinearRegression

谭贤芳

已于 2024-07-16 18:31:57 修改

阅读量852

点赞数 12

分类专栏：手撕算法文章标签：机器学习算法线性回归 c++ python

于 2024-07-06 14:58:52 首次发布

本文链接：https://blog.csdn.net/m0_70254540/article/details/140227907

版权

手撕算法专栏收录该内容

7 篇文章 0 订阅

订阅专栏

❗写在前面 ❗：

… … … … … … … … . … … …
1️⃣ 提供两种语言实现：本文提供两种语言实现，分别为Python和C++，分别依赖的核心库为Numpy和Eigen.

2️⃣ 提供关键步的公式推导：写出详尽的公式推导太花时间，所以只提供了关键步骤的公式推导。更详尽的推导，可移步其他博文或相关书籍.

3️⃣ 重在逻辑实现：文中的代码主要以实现算法的逻辑为主，在细节优化上面，并没有做更多的考虑，如果发现有不完善的地方，欢迎指出.

🌞 欢迎关注专栏，相应的内容会持续更新，祝大家变得更强！！
… … … … … … … … . … … …

1、原理推导

目标函数：

$\hat{w}^*=arg\ min(y-X\hat{w})^T(y-X\hat{w})$

令 $L=(y-X\hat{w})^T(y-X\hat{w})$ ，然后对参数 $\hat{w}$ 求导可得梯度：

$\frac{\partial L}{\partial \hat{w}^*}=2X^T(X\hat{w}-y)$ $(△)$

令上式为0就可以解得线性回归的解析解：

$\hat{w}^*=(X^TX)^{-1}X^Ty$

但有些时候，矩阵 $X^TX$ 并不是一个满秩矩阵，但可以通过对 $X^TX$ 添加一个正则化项来使得该矩阵可逆。于是有：

$\hat{w}^*=(X^TX+λI)^{-1}X^Ty$

其中 $λ I$ 为添加的正则化项。

以上是解析解的推导过程，本文不采用这种方式来进行求解参数，而是采用梯度下降的优化算法来求解 $\hat{w}^*$ 。

2、算法实现

🌈Python实现

# 基于梯度下降的线性回归算法
import numpy as np

class LinearRegression:

    def __init__(self):
        self.W = None
        self.b = None 
    
    # 计算损失和梯度
    def _linearLossGradient(self, X, y):
        # 总样本数
        m = X.shape[0]  
        # 当前参数下的预测值
        yHat = np.dot(X, self.W) + self.b
        # 当前参数才的均方误差
        loss = ((yHat - y) ** 2).sum() / m
        # 关于W的梯度
        dW = np.dot(X.T, (yHat - y)) / m
        # 关于b的梯度
        db = (yHat - y).sum() / m

        return yHat, loss, dW, db

    # 梯度下降求解最优解
    def fit(self, X, y, alpha=0.01, epochs=1000, tol=0.01):
        """
        alpha: 学习率
        epochs: 迭代次数
        tol: 停止迭代的损失最大值
        """
        # 参数W的维度 
        dim = X.shape[1]
        # 初始化参数W和截距b
        self.W = np.zeros(dim)
        self.b = 0.0

        # 记录损失的列表
        self.losslist = []

        # 迭代进行梯度下降求解W和b
        for i in range(1, epochs+1):
            _, loss, dW, db = self._linearLossGradient(X, y)
            
            self.losslist.append(loss)

            # 损失小于设定值停止迭代
            if loss < tol:
                break
            
            # 参数沿负梯度方向移动
            self.W -= alpha * dW  
            self.b -= alpha * db
            
            if (i % 500 == 0):
                print(f"epochs {i}: loss = {loss}")
    
    # 模型预测
    def predict(self, X):
        return np.dot(X, self.W) + self.b
    
    # 基于R^2的得预测得分
    def r2score(self, yPred, yTrue):
        return 1 - ((yPred - yTrue) ** 2).sum() / ((yPred - yTrue.mean()) ** 2).sum()

⚡C++实现

#include <iostream>
#include <Eigen/Dense>
#include <algorithm>

using namespace Eigen;
using namespace std;


// 声明参数的结构体
struct Theta
{
    VectorXd W;
    double b;
};

class LinearRegression
{
private:
    Theta theta;  // 模型参数

private:
    // 初始化模型参数
    void initParams(int dim)
    {
        theta.W = VectorXd::Zero(dim);
        theta.b = 0.0;
    }

    // 声明参数和损失的结构体
    struct LossGradient
    {
        VectorXd yHat;
        double loss;
        VectorXd dw;
        double db;
    };

    // 计算损失和梯度
    LossGradient linearLoss(const MatrixXd& X, const VectorXd& y)
    {   
        // 总样本数
        int m = X.rows();
        // 当前参数下的预测值
        VectorXd yHat = (X * theta.W).array() + theta.b;
        // 当前参数才的均方误差
        double loss = (yHat - y).squaredNorm();
        // 关于W的梯度
        VectorXd dw = (X.transpose() * (yHat - y)) / m;
        // 关于b的梯度
        double db = (yHat - y).sum() / m;

        return { yHat, loss, dw, db };
    }

public:
    void fit(const MatrixXd& X, const VectorXd& y, double alpha=0.01, int epochs=5000, double tol=0.01)
    {   
        // 初始化参数
        initParams(X.cols()); 

        for (int i = 0; i < epochs; i++)
        {
            // 计算损失和梯度
            auto result = linearLoss(X, y); 
            
            // 损失小于设定值停止迭代
            if (result.loss < tol)
            {
                break;
            }
            
            // 更新权重W
            theta.W -= alpha * result.dw;
            // 更新偏置b
            theta.b -= alpha * result.db;
        }
    }

    // 获取训练后的参数
    Theta getParams() const
    {
        return theta;
    }

    // 模型预测
    VectorXd predict(const MatrixXd& X)
    {
        return (X * theta.W).array() + theta.b;
    }

    // 基于R^2的得预测得分
    double score(const VectorXd& yPred, const VectorXd& yTrue)
    {
        int m = yPred.size();
        double mse = (yPred - yTrue).squaredNorm() / m;
        double var = (yTrue.array() - yTrue.mean()).square().sum() / m;

        return 1 - (mse / var);
    }
};