lsqnonlin not changing some parameter values

I am trying to fit some data with a rect function convolved with a gaussian kernel (1D) using lsqnonlin. The parameters I am trying to optimize are: rect function height, width and centre location, as well as the gaussian FWHM. When I run this, the rect height and gaussian FWHM are optimized, but the rect width and centre location are not being changed. I have normalized the parameters so their values and search ranges are ~same order magnitude, this makes no difference. The exit message is:

Optimization terminated: first-order optimality less than OPTIONS.TolFun,
 and no negative/zero curvature detected in trust region model.

Here is my code:

function[x,x0,xtruth] = fit_somedata

r = 1:17;
c = 10;
mw = 4;
ma = 1.2;
FWHM = 4;
inner_bkgd = 0;
outer_bkgd = 0;
mm_pp_rect = 0.25;
[r_rect,rect] = mk_rect(inner_bkgd,ma,outer_bkgd,c,mw,r,mm_pp_rect);

mm_pp_kernel = mm_pp_rect;
[r_kernel,gauss_kernel] = get_gauss_kernel1D(FWHM,mm_pp_kernel);
n_kernel = length(gauss_kernel);
n_hw_kernel = (n_kernel-1)/2;

conv_data = conv(gauss_kernel,rect);
conv_data = conv_data(n_hw_kernel+1:length(conv_data)-n_hw_kernel);
r_conv = r_rect;
data = interp1(r_conv,conv_data,r);

xnames = {'ma';'mw';'c';'FWHM'};
xunits = {'MMI';'units';'units';'units'};
xtruth = [ma    mw    c     FWHM];
x0 = [ma+0.1    mw+1    c-1     FWHM+1];
lb = [0         1       1       1];
ub = [2*ma      10      17      12];

options = optimset('Display','iter');

[x,resnorm,residual,exitflag,output,lambda,jacobian] = lsqnonlin(@(x) ...
    costfunction(x,inner_bkgd,outer_bkgd,r,data,mm_pp_rect),...
    x0,lb,ub,options);

end

function[cost,predicted_data] = costfunction(x,inner_bkgd,outer_bkgd,r,data,mm_pp_rect)

ma = x(1);
mw = x(2);
c = x(3);
FWHM = x(4);

[r_rect,rect] = mk_rect(inner_bkgd,ma,outer_bkgd,c,mw,r,mm_pp_rect);

mm_pp_kernel = mm_pp_rect;
[r_kernel,gauss_kernel] = get_gauss_kernel1D(FWHM,mm_pp_kernel);
n_kernel = length(gauss_kernel);
n_hw_kernel = (n_kernel-1)/2;

conv_data = conv(gauss_kernel,rect);
conv_data = conv_data(n_hw_kernel+1:length(conv_data)-n_hw_kernel);
r_conv = r_rect;
predicted_data = interp1(r_conv,conv_data,r);

cost = data - predicted_data;
end

function[r_rect_mm,rect] = mk_rect(inner_bkgd,ma,outer_bkgd,c_mm,mw_mm,r_data_mm,mm_pp_rect)

% length of rect to r & l of r_data_mm = 0:
ll_rect_mm = - min(r_data_mm);
lr_rect_mm = max(r_data_mm);
nl_rect = ceil(ll_rect_mm/mm_pp_rect);
nr_rect = ceil(lr_rect_mm/mm_pp_rect);
i_rect = -nl_rect:1:nr_rect;
r_rect_mm = i_rect*mm_pp_rect;

% boundaries bla/ma/bkgda:
lb_ma_mm = c_mm - mw_mm/2;
rb_ma_mm = c_mm + mw_mm/2;

% initialize & set rect:
rect = zeros(1,length(r_rect_mm));
rect(r_rect_mm < lb_ma_mm) = inner_bkgd;
rect(r_rect_mm > rb_ma_mm) = outer_bkgd;
i_rect_myo = ((r_rect_mm >= lb_ma_mm) + (r_rect_mm <= rb_ma_mm))>1;
rect(i_rect_myo) = ma;
end

function[r_kernel_mm,gauss_kernel,s_mm] = get_gauss_kernel1D(FWHM_mm, mm_pp_kernel)

% want kernel size allowing values to fall off to ~0 at edges
% for gaussian, take this as 3stds from the centre (origin)
% mm_pp_kernel: the sample frequency of the kernel in mm
% odd kernel length guaranteed, centred, this ensures symmetry

s_mm = FWHM_mm/(2*sqrt(2*log(2)));
hw_kernel_mm = 3*s_mm;
n_hw_kernel = ceil(hw_kernel_mm/mm_pp_kernel);
r_kernel_mm = [-n_hw_kernel:1:n_hw_kernel]*mm_pp_kernel; % guarantees odd b/c 0 in middle
c_kernel_mm = 0;
gauss_kernel = exp(-(r_kernel_mm - c_kernel_mm).^2/(2*s_mm.^2));
gauss_kernel = gauss_kernel/sum(gauss_kernel);
end

Answer 1

[Matt J]

It looks very similar to what happened here,

http://www.mathworks.com/matlabcentral/answers/122216-optimization-problem-with-objective-function-that-is-a-convolution-of-a-gaussian-and-a-step

Also, interp1 using the default linear interpolation is a non-differentiable operation. LSQNONLIN assumes that your objective function is differentiable.