❗写在前面 ❗:
… … … … … … … … . … … …
1️⃣ 提供两种语言实现:本文提供两种语言实现,分别为Python和C++,分别依赖的核心库为Numpy和Eigen.
2️⃣ 提供关键步的公式推导:写出详尽的公式推导太花时间,所以只提供了关键步骤的公式推导。更详尽的推导,可移步其他博文或相关书籍.
3️⃣ 重在逻辑实现:文中的代码主要以实现算法的逻辑为主,在细节优化上面,并没有做更多的考虑,如果发现有不完善的地方,欢迎指出.
🌞 欢迎关注专栏,相应的内容会持续更新,祝大家变得更强!!
… … … … … … … … . … … …
1、原理推导
目标函数:
w
^
∗
=
a
r
g
m
i
n
(
y
−
X
w
^
)
T
(
y
−
X
w
^
)
\hat{w}^*=arg\ min(y-X\hat{w})^T(y-X\hat{w})
w^∗=arg min(y−Xw^)T(y−Xw^)
令 L = ( y − X w ^ ) T ( y − X w ^ ) L=(y-X\hat{w})^T(y-X\hat{w}) L=(y−Xw^)T(y−Xw^),然后对参数 w ^ \hat{w} w^求导可得梯度:
∂ L ∂ w ^ ∗ = 2 X T ( X w ^ − y ) \frac{\partial L}{\partial \hat{w}^*}=2X^T(X\hat{w}-y) ∂w^∗∂L=2XT(Xw^−y) ( △ ) (△) (△)
令上式为0就可以解得线性回归的解析解:
w ^ ∗ = ( X T X ) − 1 X T y \hat{w}^*=(X^TX)^{-1}X^Ty w^∗=(XTX)−1XTy
但有些时候,矩阵 X T X X^TX XTX并不是一个满秩矩阵,但可以通过对 X T X X^TX XTX添加一个正则化项来使得该矩阵可逆。于是有:
w ^ ∗ = ( X T X + λ I ) − 1 X T y \hat{w}^*=(X^TX+λI)^{-1}X^Ty w^∗=(XTX+λI)−1XTy
其中 λ I λI λI为添加的正则化项。
以上是解析解的推导过程,本文不采用这种方式来进行求解参数,而是采用梯度下降的优化算法来求解
w
^
∗
\hat{w}^*
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);
}
};
3、总结
至此,我们分别通过Python和C++实现了线性回归算法的核心逻辑。线性回归包含了最朴素的由自变量到因变量的机器学习建模思想。基于均方误差最小化的最小二乘法是线性回归模型求解的基本方法,通过均方误差和决定系数可以评估线性回归的拟合效果。此外,线性回归模型也是其他各种线性模型的基础。