leastsq vs minimize ! Leastsq wins with an almost 100 times faster fit.

So, after some profiling and testing, I found that the difference in performance for the two fitting routines from scipy.optimize - "minimize" and "leastsq" is huge !!

Minimize seems to be at least 10 times slower than leastsq on the same data. Have a look at the code below for more details.

I used a bi exponential curve which is a sum of two exponents and the model function has 4 parameters to fit. S, f, D_star and D.

All default parameters for fitting were used and the fitting function is the IVIM model function :

S [f e^(-x * D_star) + (1 - f) e^(-x * D)]

Using the time module we got the following results :

Time taken by leastsq : 0.0003180503845214844
Time taken by minimize : 0.011617898941040039

We were in favor of minimize earlier since it allowed setting bounds for the parameters and gave the user different options for the algorithm they wanted to use for fitting such as L-BFGS-B, Truncated Newton etc. but due to the huge difference in performance, we have decided to use leastsq for fitting and use a "hack" to implement the bounds while continuing the use of leastsq. The absence of a bounded leastsq method has been long standing issue in Scipy which was resolved in the version 0.17 by the least_squares function which takes bounds and uses a "Trust Region Reflective algorithm" for fitting. 

This also meant a change the Jacobian which works out as follows :

Jacobian for IVIM function (biexponential)

Shahnawaz Ahmed

Let i be the index which denotes the number of points to fit and j be the number of parameters. In this case i is the number of bvalues. Let xi be the independent parameter (bvalue).
Let β denote the vector of parameters with βj=(S0,f,D∗,D). Thus we have j=0,1,2,3
The IVIM signal is :

S(xi,β)=β0[β1e(−xiβ2)+(1−β1)e(−xiβ3)]

Here :

β0=S0,β1=f,β2=D∗,β3=D

and

xi=b

The residual that we need to find the Jacobian for is :

ri=[yi−S(xi,β)]

The Jacobian is defined as :

−∂ri∂βj

Thus we will have :

Ji,j=−∂S∂βj

The various terms are thus :

∂S∂β0=[β1e(−xiβ2)+(1−β1)e(−xiβ3)]

∂S∂β1=β0[e−xiβ2−e−xiβ3]

∂S∂β2=β0[−xiβ1e−xiβ2]

∂S∂β3=β0[−xi(1−β1)e−xiβ3]

This should give us the Jacobian

More discussion on the issue can be found here : http://stackoverflow.com/questions/6779383/scipy-difference-between-optimize-fmin-and-optimize-leastsq

Code to compare the performance of leastsq and minimize.

Time 

import numpy as np
from scipy.optimize import minimize, leastsq
from time import time


def ivim_function(params, bvals):
    """The Intravoxel incoherent motion (IVIM) model function.

        S(b) = S_0[f*e^{(-b*D\*)} + (1-f)e^{(-b*D)}]

        S_0, f, D\* and D are the IVIM parameters.

    Parameters
    ----------
        params : array
                parameters S0, f, D_star and D of the model

        bvals : array
                bvalues

    References
    ----------
    .. [1] Le Bihan, Denis, et al. "Separation of diffusion
               and perfusion in intravoxel incoherent motion MR
               imaging." Radiology 168.2 (1988): 497-505.
    .. [2] Federau, Christian, et al. "Quantitative measurement
               of brain perfusion with intravoxel incoherent motion
               MR imaging." Radiology 265.3 (2012): 874-881.
    """
    S0, f, D_star, D = params
    S = S0 * (f * np.exp(-bvals * D_star) + (1 - f) * np.exp(-bvals * D))
    return S


def _ivim_error(params, bvals, signal):
    """Error function to be used in fitting the IVIM model
    """
    return (signal - ivim_function(params, bvals))


def sum_sq(params, bvals, signal):
    """Sum of squares of the errors. This function is minimized"""
    return np.sum(_ivim_error(params, bvals, signal)**2)

x0 = np.array([100., 0.20, 0.008, 0.0009])
bvals = np.array([0., 10., 20., 30., 40., 60., 80., 100.,
                  120., 140., 160., 180., 200., 220., 240.,
                  260., 280., 300., 350., 400., 500., 600.,
                  700., 800., 900., 1000.])
data = ivim_function(x0, bvals)

optstart = time()
opt = minimize(sum_sq, x0, args=(bvals, data))
optend = time()
time_taken = optend - optstart
print("Time taken for opt:", time_taken)


lstart = time()
lst = leastsq(_ivim_error,
              x0,
              args=(bvals, data),)
lend = time()
time_taken = lend - lstart
print("Time taken for leastsq :", time_taken)

print('Parameters estimated using minimize :', opt.x)
print('Parameters estimated using leastsq :', lst[0])