Optimization problem with objective function that is a convolution of a Gaussian and a step

Hi,

The convolution of a Gaussian and a step function results in something that looks like it models experimental data fairly well. I tried to develop an optimization routine using the convolution of a fixed Gaussian with a step function that has variable parameters, but am getting some error (I suspect a discontinuity), which causes lsqcurvefit to return only the initial guesses.

For example,

x=-10:.1:10; % the independent variable

a=10*exp(-((x)./(0.2)).^2); % the fixed Gaussian

fun = @(p,x)p(3).*conv(a,heaviside(x-(p(1)-p(2)/2)).*heaviside(-x+(p(1)+p(2)/2)),'same') % the objective function - convolution of the fixed Gaussian with a step function that has variable width, x location, and height.

data = fun([5,2,10],x); % some simulated data

data = data + 5*randn(size(x)); % add some noise

p = lsqcurvefit(fun,[1,0.9,1e3],x,data,[0,0,0],[max(x),Inf,Inf]) % the optimization problem

Any ideas about what might be going wrong, and ways to fix it?

Thanks, Lewis

Answer 1

[Matt J]

Instead of using discrete convolution, you need to write an analytical formula for the continuous convolution in your objective function (see the ERF command). If you do not, the least squares cost function, as a function of p(1) and p(2), has a staircase profile (i.e., a gazillion intervals of local minima to get stuck in) like in the figure below

    t= linspace(2, 8,1000);
    f=arrayfun( @(z) norm(fun([z, 2,10],x)-data)^2 , t);
    plot(t,f);