如何用最小二乘法拟合参数三次多项式曲线及技术问询

疑问解答 #

1. 分别拟合x/y vs 整体拟合r(u) #

两种方案都可行，适用场景不同：

• 分别拟合x=f(u)和y=g(u)：操作简单，属于线性最小二乘问题，求解效率高、结果稳定。但该方法完全依赖原参数u的分布，若原u的参数化不合理（比如不均匀、和曲线几何特性无关），拟合后的曲线可能出现参数化扭曲。
• 整体拟合r(u)：若指基于原参数u直接拟合三维向量的三次多项式，本质和分别拟合x/y等价；若指重新参数化拟合（即不固定原u，重新寻找最优参数t来匹配曲线几何形态），则属于非线性最小二乘问题，能更好地保持曲线的几何连续性，但求解复杂度更高。

2. 直接拟合r(u)的可行性 #

直接拟合r(u)分为两种情况：

• 若沿用原参数u，直接对向量r(u)拟合三次多项式，本质是对x、y分量分别拟合，和第一种方案无区别。
• 若要重新定义参数t（脱离原u的约束），直接拟合r(t)的三次多项式，则需要同时求解多项式系数和每个数据点对应的t值，属于非线性拟合场景，适合对曲线几何形态要求更高的场景。

3. 线性vs非线性最小二乘的选择 #

• 线性最小二乘：当固定原参数u，对x、y分量分别拟合三次多项式时，目标函数是关于多项式系数的线性函数，可直接通过矩阵求解（无需迭代），属于线性问题。
• 非线性最小二乘：当进行重新参数化拟合时，待估参数包括多项式系数和每个数据点的t值，目标函数是非线性的，需要迭代求解（如Gauss-Newton方法），属于非线性问题。你之前的误解是混淆了“参数化曲线的非线性”和“最小二乘问题的线性性”——只有当待估参数和目标函数存在非线性耦合时，才是非线性最小二乘。

4. scipy.least_squares的使用 #

• 线性拟合场景下，没必要用least_squares，numpy.linalg.lstsq更高效直接。
• 非线性重新参数化拟合场景下，least_squares完全可用（支持Levenberg-Marquardt方法，属于Gauss-Newton的变种），关键是正确定义残差函数并初始化参数。你之前尝试失败大概率是残差函数定义错误或参数初始化不合理。

方法论指导 #

1. 优先尝试线性方法：如果对参数化的一致性要求不高，先采用分别拟合x/y的线性最小二乘，快速得到结果，验证拟合效果。
2. 非线性方法作为补充：当线性拟合结果的参数化不符合需求（比如曲线出现不必要的扭曲），再考虑重新参数化的非线性拟合，注意参数初始化要尽可能接近最优解（比如先用线性拟合结果作为系数初始值），避免陷入局部最优。

Python代码实现 #

方案1：线性最小二乘（分别拟合x/y） #

import numpy as np
import matplotlib.pyplot as plt

# 生成带噪声的模拟数据
u = np.linspace(0, 1, 50)
x_true = 2*u + 3*u**2 - u**3
y_true = 1 - u + 4*u**3
x_noisy = x_true + np.random.normal(0, 0.05, size=u.shape)
y_noisy = y_true + np.random.normal(0, 0.05, size=u.shape)

# 构造三次多项式的设计矩阵
def build_design_matrix(u_vals):
    return np.vstack([np.ones_like(u_vals), u_vals, u_vals**2, u_vals**3]).T

# 拟合x分量
X = build_design_matrix(u)
x_coeffs, _, _, _ = np.linalg.lstsq(X, x_noisy, rcond=None)
x_fit = X @ x_coeffs

# 拟合y分量
y_coeffs, _, _, _ = np.linalg.lstsq(X, y_noisy, rcond=None)
y_fit = X @ y_coeffs

# 可视化结果
plt.figure(figsize=(8, 6))
plt.plot(x_true, y_true, label='真实曲线', color='#1f77b4')
plt.scatter(x_noisy, y_noisy, label='带噪声数据', color='#ff4b5c', s=15)
plt.plot(x_fit, y_fit, label='拟合曲线', color='#2ca02c', linestyle='--')
plt.legend()
plt.xlabel('x')
plt.ylabel('y')
plt.title('线性最小二乘拟合（分别拟合x/y）')
plt.show()

print("x分量拟合系数：", x_coeffs.round(4))
print("y分量拟合系数：", y_coeffs.round(4))

方案2：非线性最小二乘（重新参数化拟合） #

import numpy as np
from scipy.optimize import least_squares
import matplotlib.pyplot as plt

# 生成带噪声的模拟数据（原参数u随机分布）
u = np.random.uniform(0, 1, 50)
u.sort()
x_true = 2*u + 3*u**2 - u**3
y_true = 1 - u + 4*u**3
x_noisy = x_true + np.random.normal(0, 0.05, size=u.shape)
y_noisy = y_true + np.random.normal(0, 0.05, size=u.shape)
data_points = np.vstack([x_noisy, y_noisy]).T

# 定义残差函数：参数包含多项式系数和每个点的t值
def calc_residual(params, data):
    num_points = data.shape[0]
    # 前8个参数是x、y的三次多项式系数：A0,A1,A2,A3,B0,B1,B2,B3
    poly_coeffs = params[:8]
    # 后num_points个参数是每个数据点对应的t值
    t_vals = params[8:]
    
    # 计算拟合的x、y坐标
    x_fit = poly_coeffs[0] + poly_coeffs[1]*t_vals + poly_coeffs[2]*t_vals**2 + poly_coeffs[3]*t_vals**3
    y_fit = poly_coeffs[4] + poly_coeffs[5]*t_vals + poly_coeffs[6]*t_vals**2 + poly_coeffs[7]*t_vals**3
    
    # 返回扁平化的残差数组
    return np.hstack([x_fit - data[:,0], y_fit - data[:,1]])

# 初始化参数：先用线性拟合得到系数，t初始化为原u
X = np.vstack([np.ones_like(u), u, u**2, u**3]).T
init_x_coeffs, _, _, _ = np.linalg.lstsq(X, x_noisy, rcond=None)
init_y_coeffs, _, _, _ = np.linalg.lstsq(X, y_noisy, rcond=None)
initial_params = np.hstack([init_x_coeffs, init_y_coeffs, u])

# 用Levenberg-Marquardt方法求解非线性最小二乘
fit_result = least_squares(calc_residual, initial_params, args=(data_points,), method='lm')

# 提取最优参数
opt_poly_coeffs = fit_result.x[:8]
opt_t_vals = fit_result.x[8:]

# 生成拟合曲线的密集点用于可视化
t_dense = np.linspace(opt_t_vals.min(), opt_t_vals.max(), 100)
x_fit_dense = opt_poly_coeffs[0] + opt_poly_coeffs[1]*t_dense + opt_poly_coeffs[2]*t_dense**2 + opt_poly_coeffs[3]*t_dense**3
y_fit_dense = opt_poly_coeffs[4] + opt_poly_coeffs[5]*t_dense + opt_poly_coeffs[6]*t_dense**2 + opt_poly_coeffs[7]*t_dense**3

# 可视化结果
plt.figure(figsize=(8, 6))
plt.plot(x_true, y_true, label='真实曲线', color='#1f77b4')
plt.scatter(x_noisy, y_noisy, label='带噪声数据', color='#ff4b5c', s=15)
plt.plot(x_fit_dense, y_fit_dense, label='拟合曲线', color='#ff7f0e', linestyle='--')
plt.legend()
plt.xlabel('x')
plt.ylabel('y')
plt.title('非线性最小二乘拟合（重新参数化）')
plt.show()

print("最优多项式系数：", opt_poly_coeffs.round(4))

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