optimization of a problem using arrayfun

Dear all!

I am currently facing a problem optimizing a large equation made up of many sum-terms.

let's say, for the sake of simplicity, that the basic function has some gaussian-like character

f_i(p)=p(1).*exp(-(p(2)-xdata).^2/p(3))

and i is an index running from 1 to something like 20. Since the function looks the same for every set of parameters p (size(p)=20,3) it would be nice to "recycle" the basic formula for every set of parameters, and I want to use arrayfun

(f)=fun1(p)
f(p)=p(1).*exp(-(p(2)-xdata).^2./p(3))
end

[sse]=fun2(p)
X=arrayfun(@(k) fun1(params(k,:)), 1:size(p,1), 'UniformOutput',false);
f=sum(cell2mat(X),2);
sse=((f-ydata).^2);
end

this basically works because I can plot the result and I see all my nice gaussian-like peaks.

However when I pass this to an optimization

x=lsqnonlin(@fun2, p);

it is not working (or the start condition is too poor, but from visual inspection the start point is reasonably close to the measured data, see the attached file).

is there some conceptual error in my approach?

Any help is appreciated!

Yours, Chris

Answer 1

[Matt J]

There is no conceptual error in your approach, that I can see. The code you have posted does not run, however, so is obviously different from the version you are actually using yourself. This makes it difficult to comment on the estimation failure you are seeing.

One potential problem with your code (as posted) is that you are returning (f-ydata).^2 whereas lsqnonlin expects fun2 to return just (f-ydata). [Side note: using lsqcurvefit would be more appropriate here, since then you only have to return f.] This causes lsqnonlin to minimize the 4th order error sum((f-ydata).^4), which might be making things ill-determined.

I also notice from your plots that your function has units of 1e-13. You should probably normalize it to something on the order of 1, so that the tolerance parameters of lsqnonlin don't trigger an early stop.

You should also call lsqnonlin or lsqcurvefit with more output arguments so that you know under what conditions it exited,

   [x,resnorm,residual,exitflag,output] = lsqnonlin(...)