最小二乘拟合正弦幂级数

Question

我正在尝试拟合形式的函数：

其中 A 和 B 是固定常数。在 scipy 中，我通常（并且我认为是相当规范的）处理此类问题的方法如下：

def func(t, coefs):
    phase = np.poly1d(coefs)(t)
    return A * np.cos(phase) + B

def fit(time, data, guess_coefs): 
    residuals = lambda p: func(time, p) - data
    fit_coefs = scipy.optimize.leastsq(residuals, guess_coefs) 
    return fit_coefs

这行得通，但我想提供一个分析雅可比行列式来提高收敛性。因此：

def jacobian(t, coefs):
    phase = np.poly1d(coefs, t)
    the_jacobian = []
    for i in np.arange(len(coefs)):
        the_jac.append(-A*np.sin(phase)*(t**i))
    return the_jac

def fit(time, data, guess_coefs):
    residuals = lambda p: func(time, p) - data
    jac = lambda p: jacobian(time, p)
    fit_coefs = scipy.optimize.leastsq(residuals, guess_coefs, 
                                       Dfun=jac, col_deriv=True)

即使是 2 个或更少的订单，这也不起作用。使用 optimize.check_gradient() 进行快速检查也不会产生积极的结果。

我几乎可以肯定 Jacobian 和代码是正确的（尽管请纠正我）并且问题更根本：Jacobian 中的 t**i 项会导致溢出错误。这不是函数本身的问题，因为这里的单项式项乘以它们的系数，系数非常小。

我的问题是：

1. 我在上面所做的事情在代码方面有问题吗？
2. 还是有其他问题？
3. 如果我的假设是正确的，有没有办法对拟合函数进行预处理，以便雅可比行列式表现更好？也许我可以拟合数据和时间的对数，或者其他什么。

谢谢！

编辑：忘记了原始函数形式中的正方形

Answer 1

poly1D 函数首先具有最高系数，而您的 jacobian 函数首先假定最低系数。如果在jacobian 中你做出了返回语句return the_jac[::-1]（并且还修复了更明显的拼写错误）你的函数将通过optimize.check_gradient() 并在leastsq() 中正常工作。

您关于数值稳定性的进一步问题在这里也是有道理的。如果您有较大的 t 值和大量系数，则很容易出现数值精度问题：例如，在 32 位精度下，如果您的多项式计算结果超过 10^8，则正弦的相位将完全丢失在尾数中，正弦评估的结果基本上是垃圾！

您可以通过在多项式中使用 模幂运算 来解决这个问题：基本上，您在多项式的每个项中关心的是(a_k t ** k) % p，其中p = 2 * np.pi 是正弦曲线的周期。您可以使用 (a_k * (t % (p / a_k)) ** k) % p 为大型 t 计算此模块化指数以获得更高的精度；为了准确k 以及，事情变得有点复杂。请参阅 this answer 了解有关此类问题的精彩讨论。