Python实现Newton-Gauss高斯牛顿法拟合二次函数出现奇异矩阵报错求解

你的代码存在两处核心错误，直接导致了奇异矩阵报错和后续结果错误：

• 第一处：雅可比矩阵计算逻辑完全错误
  目标拟合模型为Y = C1 + C2*X + C3*X²，残差定义为r = Y - (C1 + C2*X + C3*X²)，雅可比矩阵的每一列对应残差对一个待求参数的偏导数，正确计算规则如下：
  • 对参数C1的偏导：-1
  • 对参数C2的偏导：-X
  • 对参数C3的偏导：-X²
    你原代码中错误将J[:,1]设为P0[0]、J[:,2]设为P0[2]*X，和模型偏导完全不匹配，构造出的J矩阵不满秩，最终np.dot(J.T, J)为奇异矩阵无法求逆。
• 第二处：预测值计算未适配二次模型
  你得到拟合参数后仍然沿用了教程中分式模型的预测公式pred = C1*X/(C2+X)，没有改为二次函数的计算逻辑，即使迭代成功也会得到错误的拟合结果。

修正后完整代码 #

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

def gauss_newton(X, Y, max_iter=1000, eps=1e-6): 
    P0 = np.array([1.0, 1.0, 1.0]) # 转成numpy数组方便运算
    n_params = len(P0)
    J = np.zeros([len(X), n_params])

    for i in range(max_iter):
        # 修正雅可比矩阵计算
        J[:,0] = -1
        J[:,1] = -X
        J[:,2] = -X**2

        r = Y - (P0[0] + P0[1]*X + P0[2]*X**2)
        # 用伪逆替换普通逆，避免偶尔的奇异问题，鲁棒性更强
        t1 = np.linalg.pinv(np.dot(J.T, J))
        t2 = np.dot(t1, J.T)
        t3 = np.dot(t2, r)

        P1 = P0 - t3
        t4 = np.abs(P1-P0)
        if np.max(t4) <= eps:
            break
        P0 = P1

    return P0[0], P0[1], P0[2]

X = np.asarray([1, 2, 3, 4, 5, 6])
Y = np.asarray([5, 7, 9, 11, 13, 11])
C1, C2, C3 = gauss_newton(X, Y)
# 修正预测值为二次模型
pred = C1 + C2 * X + C3 * X**2

plt.figure(1, figsize=(6,4), dpi=120)
plt.scatter(x=X, y=Y, c='green', marker='o', label='原始数据')
plt.plot(X, pred, '--m', label='拟合二次模型')
plt.legend()
plt.show()

内容的提问来源于stack exchange，提问作者Никита Михалков