Scipy.optimize.leastsq returns the best guessing options

Question

I want to fit the Lorentzian peak to the dataset x and y, the data is fine. Other programs like OriginLab are perfect, but I wanted to automate the fitting with python, so I have the code below, which is based on http://mesa.ac.nz/?page_id=1800

The problem is that scipy.optimize.leastsq returns, as best fit, the same initial guess parameters that I passed in that I did practically nothing. Here is the code.

#x, y are the arrays with the x,y axes respectively #defining funcitons def lorentzian(x,p): return p[2]*(p[0]**2)/(( x - (p[1]) )**2 + p[0]**2) def residuals(p,y,x): err = y - lorentzian(x,p) return err p = [0.055, wv[midIdx], y[midIdx-minIdx]] pbest = leastsq(residuals, p, args=(y, x), full_output=1) best_parameters = pbest[0] print p print pbest 

p are the initial guesses, and best_parameters are the least-squares returned best-fit parameters, but they are always the same.

this is what full_output = 1 is returned (long numeric arrays have been shortened, but still represent)

  [0.055, 855.50732, 1327.0] (array([ 5.50000000e-02, 8.55507324e+02, 1.32700000e+03]), None, {'qtf':array([ 62.05192947, 69.98033905, 57.90628052]), 'nfev': 4, 'fjac': array([[-0., 0., 0., 0., 0., 0., 0.,], [ 0., -0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.], [ 0., 0., -0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]]), 'fvec': array([ 62.05192947, 69.98033905, 53.41218567, 45.49879837, 49.58242035, 36.66483688, 34.74443436, 50.82238007, 34.89669037]), 'ipvt': array([1, 2, 3])}, 'The cosine of the angle between func(x) and any column of the\n Jacobian is at most 0.000000 in absolute value', 4) 

Can anyone understand what is wrong?

Answer 1

A quick Google search indicates a problem with single data precision (your other programs almost certainly upgrade to double precision too, although this is clearly a problem for scipy too, see also this error report ). If you look at your result full_output=1 , you will see that the Jacobian is approximately equal to zero. Therefore, explicitly providing the Jacobians can help (although even then you might want to take off, because the minimum accuracy of the relative error that you can get with single precision is very limited).

Solution: the simplest and numerically best solution (of course, giving a real Jacobian is also a bonus) is simply converting the x and y data to double precision (for example, x = x.astype(np.float64) ).

I would not suggest this, but you can also fix it by setting the epsfcn keyword argument (as well as the admission keyword arguments possibly) manually, something like epsfcn=np.finfo(np.float32).eps . It seems that the problem is fixed, but (since most calculations are performed with scalars, and scalars do not lead to overflow in your calculation), the calculations are performed in float32, and the loss accuracy seems rather large, at least when not providing Dfunc.