线性最小二乘拟合数据点效果不佳，求原因及修正方案

问题描述 #

使用以下代码进行2D平面的线性最小二乘拟合时，尽管两种方法得到的系数一致，但拟合平面效果严重失真：

import numpy as np
from scipy import linalg as la
from scipy.linalg import solve

# data
f1 = np.array([1., 1.5, 3.5, 4.])
f2 =  np.array([3., 4., 7., 7.25])
# z = np.array([6., 6.5, 8., 9.])
A= X= np.array([f1, f2]).T

b= y= np.array([0.5, 1., 1.5, 2.]).T

##################### la.lstsq

res= la.lstsq(A,b)[0]
print(res)    
##################### custom lu

#custom OLS 
def ord_ls(X, y):
    A = X.T @ X
    b = X.T @ y
    beta = solve(A, b, overwrite_a=True, overwrite_b=True,
                 check_finite=True)
    return beta

res = ord_ls(X, y)
print(res)

##################### plot

# use the optimized parameters to plot the fitted curve in 3D space.
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# Create 3D plot of the data points and the fitted curve
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(f1, f2, y, color='blue')
x_range = np.linspace(0, 7, 100)
y_range = np.linspace(0, 7,100)

X, Y = np.meshgrid(x_range, y_range)
Z = res[0]*X + res[1]

ax.plot_surface(X, Y, Z, color='red', alpha=0.5)
ax.set_xlabel('feat.1')
ax.set_ylabel('feat.2')
ax.set_zlabel('target')
plt.show()

# [0.2961165  0.09475728]
# [0.2961165  0.09475728]

疑问：拟合平面失真的原因是什么？如何修正？是否需要正则化？是否是特征共线性导致的？

附：需将scipy-0.18.0文档第184页的LU分解代码p, l, u = la.lu(A)翻译成中文。

问题原因 #

1. 核心问题：缺失截距项
   标准线性回归模型应为 z = β₀ + β₁*f1 + β₂*f2，其中β₀是截距项（代表当所有特征为0时的预测值）。你的代码中仅对f1和f2拟合了系数β₁、β₂，完全忽略了截距β₀，导致拟合平面强制过原点，这是效果失真的主要原因。

2. 次要问题：特征高度共线性
   计算f1和f2的相关系数：np.corrcoef(f1, f2)[0,1]结果约为0.999，说明两个特征高度线性相关。共线性会导致系数估计的方差变大，稳定性下降，但当前平面失真的核心原因是缺失截距，而非共线性。

修正方法 #

1. 添加截距项到特征矩阵 #

在特征矩阵X中添加一列全1的常数项，用于拟合截距β₀：

# 修改数据部分：添加截距项
X = np.hstack([np.ones((len(f1),1)), np.array([f1, f2]).T])

2. 修正拟合与绘图代码 #

修正后的完整代码：

import numpy as np
from scipy import linalg as la
from scipy.linalg import solve
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# data
f1 = np.array([1., 1.5, 3.5, 4.])
f2 = np.array([3., 4., 7., 7.25])
y = np.array([0.5, 1., 1.5, 2.]).T

# 添加截距项：特征矩阵第一列为全1
X = np.hstack([np.ones((len(f1), 1)), np.array([f1, f2]).T])

##################### la.lstsq
res_lstsq = la.lstsq(X, y)[0]
print("lstsq结果（截距β0, β1, β2）:", res_lstsq)

##################### custom lu
def ord_ls(X, y):
    A = X.T @ X
    b = X.T @ y
    beta = solve(A, b, overwrite_a=True, overwrite_b=True, check_finite=True)
    return beta

res_custom = ord_ls(X, y)
print("自定义OLS结果（截距β0, β1, β2）:", res_custom)

##################### plot
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(f1, f2, y, color='blue')

x_range = np.linspace(0, 7, 100)
y_range = np.linspace(0, 7, 100)
X_mesh, Y_mesh = np.meshgrid(x_range, y_range)

# 计算拟合平面：Z = β0 + β1*X + β2*Y
Z_mesh = res_custom[0] + res_custom[1]*X_mesh + res_custom[2]*Y_mesh

ax.plot_surface(X_mesh, Y_mesh, Z_mesh, color='red', alpha=0.5)
ax.set_xlabel('feat.1')
ax.set_ylabel('feat.2')
ax.set_zlabel('target')
plt.show()

运行后会得到包含截距的系数（例如[0.05343511 0.31007752 0.08830601]），拟合平面会贴合数据点，解决失真问题。

关于正则化与共线性 #

• 当前问题无需正则化，修正截距后即可得到合理的拟合结果。
• 若后续因共线性导致系数波动过大（比如新增数据后系数变化剧烈），可考虑使用**岭回归（Ridge Regression）**等正则化方法，通过添加L2惩罚项降低系数方差。

LU分解代码翻译 #

p, l, u = la.lu(A)
对矩阵A进行LU分解，返回三个矩阵：置换矩阵p、下三角矩阵l、上三角矩阵u，满足数学关系：p @ A = l @ u。

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