高斯-牛顿法求解最小二乘问题-python

拟合曲线

其中 a, b, c 为曲线的参数，w 为高斯噪声，满足 w ∼ (0, σ2 )。

问题本质就是求解下面的最小二乘问题以估计曲线参数：

拟合结果

参数估计结果，与设置的a=1, b=2, c=1还是比较接近的

直接上代码：

import numpy as np
import matplotlib.pyplot as plt

# 定义真实参数值
ar, br, cr = 1.0, 2.0, 1.0

# 定义估计参数值
ae, be, ce = 2.0, -1.0, 5.0

# 数据点数量
N = 100

# 噪声标准差
w_sigma = 1.0

# 生成数据
rng = np.random.default_rng()
x_data = np.linspace(0, 1, N)
y_data = np.exp(ar * x_data**2 + br * x_data + cr) + rng.normal(0, w_sigma, N)

# 高斯牛顿法迭代
iterations = 100
inv_sigma = 1.0 / w_sigma

# 开始迭代
for _ in range(iterations):
    H = np.zeros((3, 3))
    b = np.zeros(3)
    cost = 0

    for i in range(N):
        xi, yi = x_data[i], y_data[i]
        error = yi - np.exp(ae * xi**2 + be * xi + ce)
        J = np.array([
            -xi**2 * np.exp(ae * xi**2 + be * xi + ce),
            -xi * np.exp(ae * xi**2 + be * xi + ce),
            -np.exp(ae * xi**2 + be * xi + ce)
        ])

        H += inv_sigma**2 * np.outer(J, J)
        b += -inv_sigma**2 * error * J

        cost += error**2

    # 求解线性方程 Hx=b
    dx = np.linalg.solve(H, b)

    # 更新参数
    ae += dx[0]
    be += dx[1]
    ce += dx[2]

    print(f"Total cost: {cost}, Update: {dx}, Estimated params: {ae}, {be}, {ce}")

# 绘制拟合曲线和原始数据
x_range = np.linspace(0, 1, 100)
y_fit = np.exp(ae * x_range**2 + be * x_range + ce)

plt.scatter(x_data, y_data, label='Original data')
plt.plot(x_range, y_fit, label='Fitted curve', color='red')
plt.legend()
plt.xlabel('x')
plt.ylabel('y')
plt.title('Nonlinear Least Squares Fit with Gauss-Newton')
plt.show()

运行结果