multithreading optimization in python

Question

I am using lmfit to do some optimization and it is tediously slow. I have a large image and I am basically running least square minimization at every pixel. It seems to me that perhaps it is ideal for multithreading or other types of optimizations as it is painfully slow at the moment.

So, my optimization code is as follows:

I define the objective function as:

def objective(params, x, data):
    s0 = params['S0']
    t1 = params['T1']
    m = s0 * (1.0 - np.exp(-x / t1))
    return m - data

So, I am trying to minimise the difference between the model and observation. I think lmfit ensures that the absolute value is minimized but I am not sure and needs to check.

The main loop is as follows:

The parameters I need to estimate are as initialized as follows with some initial values:

p = Parameters()
p.add('S0', value=1.0)
p.add('T1', value=3.0)

final_data = <numpy image>
tis = np.asarray([1.0, 2.0, 3.0])

for i in range(final_data.shape[0]):
    print "Processing pixel: ", i
    minner = Minimizer(objective, params=p, fcn_args=np.asarray(tis),
                       final_data[i, :]),
                       nan_policy='propagate')
    result = minner.minimize(method='least_squares')
    s0_data[i] = result.params['S0']
    t1_data[i] = result.params['T1']

This works fine but it is tediously slow. I was trying to figure out how to do multithreading in python and got thoroughly confused about posts regarding GIL locking and that multithreading in python does not really exist.

My question is: 1: Can this be scaled with multithreading in python easily? 2: Is there any other optimizations that I can try?

Answer 1

As the comments suggest, multithreading here won't be very fruitful. Basically, any single fit with lmfit or scipy ends up with a single-threaded fortran routine calling your python objective function repeatedly, and using those results to generate the next step. Trying to use multithreading means that the python objective function and parameters have to be managed among the threads -- the fortran code is not geared for this, and the calculation is not really I/O bound anyway.

Multiprocessing in order to use multiple cores is a better approach. But trying to use multiprocessing for a single fit is not as trivial as it sounds, as the objective function and parameters have to be pickle-able. For your simple example, this should work, but the approach can break as the problem uses more complex objects. The dill package can help with that.

But also: there is an even easier solution for your problem, as it is naturally parallelized. Just to do a separate fit per pixel, each in their own process. You can use multiprocessing for that, or you could even break the problem up into N separate real processes run N different scripts each fitting 1/N of the pixels.