Least squares fit to sinusoidal power series

Question

I'm trying to fit a function of form:

where A and B are fixed constants. In scipy, my usual (and I think reasonably canonical) approach to such problems goes like so:

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

This works OK, but I'd like to supply an analytic Jacobian to improve convergence. Thus:

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)

This doesn't work as well, even at orders of 2 or less. Neither does a quick check with optimize.check_gradient() give positive results.

I'm almost certain the Jacobian and the code are correct (although please correct me) and that the problem is more fundamental: the t**i term in the Jacobian causes overflow errors. This isn't a problem in the function itself, because here the monomial terms are multiplied by their coefficients, which are very small.

My questions are:

1. Is there something wrong code-wise with what I've done above?
2. Or is there some other problem?
3. If my hypothesis is correct, is there a way to precondition the fit-function so the Jacobian is better-behaved? Maybe I could fit to logarithms of the data and time, or something.

Thanks!

Edit: forgot the square in the original functional form

Answer 1

The poly1D function has the highest coefficient first, while your jacobian function assumes the lowest coefficient first. If in jacobian you make the return statement return the_jac[::-1] (and also fix the more obvious typos) your functions will pass optimize.check_gradient() and work correctly in leastsq().

Your further question about numerical stability is also warranted here. If you have large values of t and large numbers of coefficients, you could very easily have numerical accuracy issues: for example, in 32-bit precision, if your polynomial evaluates to over about 10^8, the phase of the sinusoid will be entirely lost in the mantissa and the result of the sine evaluation will be basically garbage!

You could address this by using modular exponentiation within the polynomial: essentially all you care about in each term of the polynomial is (a_k t ** k) % p, where p = 2 * np.pi is the period of the sinusoid. You can compute this modular exponential to higher precision for large t with (a_k * (t % (p / a_k)) ** k) % p; for accuracy at large k as well, things get a bit more complicated. See this answer for a nice discussion of such things.