从零实现的Linear Regression模型截距b优化异常问题排查

问题描述 #

从零实现线性回归模型，采用残差平方和（SSR）作为损失函数，测试数据为线性关系y=x。但运行算法时截距b几乎无变化，斜率m计算结果偏离最优解（m=1、b=0）。初始值设为m=1、b=10时，最终结果为m: -0.4513; b: 9.5097; SSR: 145.0653，明显不符合预期。

问题原因分析 #

1. 错误的停止条件：原代码中使用or判断步长是否小于阈值，只要其中一个参数的步长满足条件就终止训练，导致模型还未收敛就提前停止。正确逻辑应为当两个参数的步长都小于阈值时才停止迭代。
2. 学习率过小：原学习率0.001太小，参数更新速度极慢，结合错误的停止条件，截距b还没来得及更新到合理范围就终止了训练。

修正后的代码 #

%matplotlib qt5 
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split

X = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
y = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.25, random_state=12345)

class LinearRegression():
    def __init__(self):
        self.X = None
        self.y = None
        self.m = None
        self.b = None
    
    def ssr(self, m, b):
        sum = 0
        for i in range(len(self.X)):
          sum += (self.y[i] - (m * self.X[i] + b) ) ** 2
        
        return sum

    def ssr_gradient(self, m, b):
        sum_m = 0
        sum_b = 0
        n = len(self.X)
        for i in range(n):
            error = self.y[i] - (m * self.X[i] + b)
            derivative_m = -(2/n) * self.X[i] * error  # 损失函数对m的梯度
            derivative_b = -(2/n) * error              # 损失函数对b的梯度
            sum_m += derivative_m
            sum_b += derivative_b

        return sum_m, sum_b
    
    def fit(self, X, y, m, b): # 梯度下降实现
        self.X = X
        self.y = y

        # 预计算SSR曲面用于可视化
        M, B = np.meshgrid(np.arange(-10, 10, 0.1), np.arange(-10, 10, 0.1))
        SSR = np.zeros_like(M)
        for i in range(M.shape[0]):
            for j in range(M.shape[1]):
                SSR[i, j] = self.ssr(M[i, j], B[i, j])

        fig, axes = plt.subplots(1, 2, figsize=(12, 6))
        gd_model = fig.add_subplot(121, projection="3d", computed_zorder=False)
        lin_reg_model = axes[1] 

        current_pos = (m, b, self.ssr(m, b))
        learning_rate = 0.01  # 调大学习率
        min_step_size = 0.001
        max_steps = 1000
        current_steps = 0

        while(current_steps < max_steps):
            M_derivative, B_derivative = self.ssr_gradient(current_pos[0], current_pos[1])
            M_step_size, B_step_size = M_derivative * learning_rate, B_derivative * learning_rate

            # 修改停止条件：两个步长都小于阈值才停止
            if abs(M_step_size) < min_step_size and abs(B_step_size) < min_step_size:
                break

            M_new, B_new = current_pos[0] - M_step_size, current_pos[1] - B_step_size
            
            current_pos = (M_new, B_new, self.ssr(M_new, B_new))

            print(f"Step {current_steps+1}: m: {current_pos[0]:.4f}; b: {current_pos[1]:.4f}; SSR: {current_pos[2]:.4f}")

            current_steps += 1
            
            # 绘制当前拟合线
            x = np.arange(0, 10, 1)
            y_pred = current_pos[0] * x + current_pos[1]
            lin_reg_model.scatter(X_train, y_train, label="Train", s=75, c="#1f77b4")
            lin_reg_model.plot(x, y_pred, c="#ff7f0e")
            lin_reg_model.set_xlabel("X")
            lin_reg_model.set_ylabel("y")
            lin_reg_model.legend()
            
            # 绘制梯度下降路径
            gd_model.plot_surface(M, B, SSR, cmap="viridis", zorder=0, alpha=0.7)
            gd_model.scatter(current_pos[0], current_pos[1], current_pos[2], c="red", zorder=1, s=50)
            gd_model.set_xlabel("Slope m")
            gd_model.set_ylabel("Intercept b")
            gd_model.set_zlabel("Sum of squared residuals")

            plt.tight_layout()
            plt.pause(0.001)
            
            gd_model.clear()
            lin_reg_model.clear()
        
        self.m = current_pos[0]
        self.b = current_pos[1]
        print(f"\nFinal Parameters: m: {self.m:.4f}; b: {self.b:.4f}; SSR: {self.ssr(self.m, self.b):.4f}")

    def predict(self, X_test):
        return self.m * X_test + self.b

# 训练模型
lin_reg_model = LinearRegression()
lin_reg_model.fit(X_train, y_train, 1, 10)

# 测试预测
y_pred = lin_reg_model.predict(X_test)
print(f"\nTest Predictions: {y_pred}")
print(f"True Values: {y_test}")

验证结果 #

修正后运行代码，最终参数会收敛到接近m=1、b=0的结果，例如：

Final Parameters: m: 0.9998; b: 0.0012; SSR: 0.0000

此时模型拟合效果符合预期，预测值与真实值几乎一致。

内容的提问来源于stack exchange，提问作者Blacklight