scipy.optimize.leastsq fails to fit simple model

Question

I have been trying to use python's scipy.optimize library to estimate the parameters of a model but without success so far. I tried to use scipy.optimize.leastsq which uses the levenberg-marquardt algorithm. Unfortunately, it always fails to find the minimum of my model function even if I set the initial parameter guess very close to the best fit. Actually, it always returns the initial guessed parameters. So, I assume I am doing something wrong. My model is a simple circle but in order to make things even simpler only the radius is the actual parameter, the center of the circle in the data is known and hard-coded. The data is a 10x10 pixel float image with a circle with center 5,5 and radius 4. Actually, the data have been generated using the model I am trying to fit. So, a perfect fit exists. Here is my code:

import math
import numpy
import scipy.optimize

# ========================================================================

g_data_width = 10
g_data_height = 10
g_xo = 5.0
g_yo = 5.0

def evaluate_model01(x,y,r):
    x2 = x*x
    y2 = y*y
    r2 = r*r
    v = 0.0
    if(x2 + y2 <= r2):
        v = 20.0
    return v

def model01(params,data_o):
    data_r = numpy.zeros(g_data_height*g_data_width)
    r = params[0]
    for y in range(g_data_height):
        for x in range(g_data_width):
            xnew = x - g_xo
            ynew = y - g_yo
            data_r[y*g_data_width+x] = math.fabs(data_o[y,x]-evaluate_model01(xnew,ynew,r))
    return data_r

# ========================================================================

g_data_o = numpy.array([[ 0,  0,  0,  0,  0,  0,  0,  0,  0,  0],
                        [ 0,  0,  0,  0,  0, 20,  0,  0,  0,  0],
                        [ 0,  0,  0, 20, 20, 20, 20, 20,  0,  0],
                        [ 0,  0, 20, 20, 20, 20, 20, 20, 20,  0],
                        [ 0,  0, 20, 20, 20, 20, 20, 20, 20,  0],
                        [ 0, 20, 20, 20, 20, 20, 20, 20, 20, 20],
                        [ 0,  0, 20, 20, 20, 20, 20, 20, 20,  0],
                        [ 0,  0, 20, 20, 20, 20, 20, 20, 20,  0],
                        [ 0,  0,  0, 20, 20, 20, 20, 20,  0,  0],
                        [ 0,  0,  0,  0,  0, 20,  0,  0,  0,  0]],dtype=numpy.float32)
g_params = numpy.array([8.0])
print(scipy.optimize.leastsq(model01,g_params,args=(g_data_o),full_output=1))

# ========================================================================

I hard-coded the data in order to remove any data dependencies and allow the code to run out of the box on any machine with scipy installed. What I don't quite understand is what the model01 function is supposed to return. According to the documentation it is supposed to return an array. An array of what? In my code I assumed that I have to return an array of the residual for each data point. Is that correct? My data is a 2D array because it is an image but my residual is a flattened 2D array of residual. Is that OK? Can someone tell me what exactly I am doing wrong? Can someone modify and fix the code? As I mentioned above, the code should run out of the box on any machine with scipy and numpy installed. In case what I want to achieve is not possible with scipy.optimize.leastsq, can you suggest some other library which fits models using the levenberg-marquardt algorithm?

Answer 1

I am afraid you can't solve your particular problem using leastsq. leastsq is a wrap around FORTRAN library minipack, it calls MINPACK’s lmdif and lmder algorithms. Importantly, it is based on the Jacobian and the Hessian of the least squares objective function. Your target function does not have smooth derivative, due to:

if(x2 + y2 <= r2):
    v = 20.0

and

math.fabs(......)

, so leastsq will always return the initial start parameter.

You should try to use some method which doesn't require gradient/derivative, such as Powell’s method fmin_powellor Nelder-Mead fmin.

Regarding what is going on in your model101(). It first makes a 1D array of the size of width*height of g_data, called data_r. Then it iterates over g_data, calculate math.fabs(data_o[y,x]-evaluate_model01(xnew,ynew,r)) for each element and put that value in the 1D array data_r. Finally, it returns date_r.

Your interpretation is correct, model101() returns the flatten 1D residue array for your 2D data.