Selectively set parameters to fit to 0 in lsqnonlin

Hi everyone,

I was wondering if there's a more elegant and less computing-time consuming way to implement the following problem.

I'm performing a lsqnonlin fit of some objective function with let's say 4 parameters to some data. I now want to try to set single or combinations of parameters to 0, fit again and compare the results.

What I do in a run script is define a matrix describing which parameters should be set to 0 (0 in matrix) and which ones should be present (1 in matrix). The rows determine the parameter combinations and the columns the presence of the parameter. So, for example:

 comb = [0,0,0,1; 0,1,0,1; 1,1,1,1];
 for i = 1:size(comb,1)
   combRow = comb(i, :);
   [pars, chi2] = fitFcn(xdata, ydata, p0, combRow)
 end

My function doing one fit looks like this:

function [pars, chi2] = fitFcn(xdata, ydata, p0, combRow)

Now I call lsqnonlin by

 [pars, chi2] = lsqnonlin(@objFcn, p0, [], [], [], xdata, ydata);

And then the objective function is defined as:

 function error = objFcn(p, x, y)
    q(combRow == 0) = 0; % set non-existent parameters to 0
    q(combRow ~= 0) = p(1:nnz(combRow)); % generate parameter vector with just as many elements as non-zeros in the specific row of combRow

    % some fitting function
    fit = q(1)*x + q(2)*x^2 + q(3)*x^3 + q(4)*x^4;

    error = fit - y;
 end

Ok, so that's my algorithm. Keep in mind that I just generated a small example, because the real thing is hell-a-lot more complex. That's why I'm hoping that someone has an idea how to implement in a way that it's computed much faster. Because when I profile the objective function (in my real code) the lines

 q(combRow == 0) = 0;
 q(combRow ~= 0) = p(1:nnz(combRow));

need almost 20% of the time and that is just too much.

I'll appreciate any help and/or hint.

Cheers, Pascal

Answer 1

[Matt Tearle]

Have you tried using the bounds to fix the parameters to be zero? You could set the lower bound to 0 and the upper bound to realmin (or something else tiny). Example:

% Make some data
c = rand(3,1)
x = linspace(0,pi)';
y = polyval(c,x).*(1+0.2*randn(size(x)));
plot(x,y,'o')
% Make objective function
f = @(c,x) polyval(c,x);
% Solve with all parameters
lsqcurvefit(f,rand(3,1),x,y)
% Solve with one parameter fixed to 0
lsqcurvefit(f,rand(3,1),x,y,[-Inf,0,-Inf],[Inf,realmin,Inf])

[Pascal Schulthess]

I've tried this. But I set the bounds both to zeros which isn't accepted by the solver. I didn't know about realmin but I'll give it a shot and report back. Thanks a lot

Edit: It's not really faster. After some trial runs I get that my method always needs around 11 secs while your suggestion finishes after about 15 secs.

[Matt Tearle]

Oh well. Plan B...

Do any of the coefficients you want to set to zero actually show up nonlinearly in the model? It seems like you're trying to "turn off" terms in the model, in which case I would assume these are linear combinations.

You say the real thing is complex -- can you give some idea of the sizes of things? The number of parameters, how many you typically are turning off, how many different model combinations you're trying, etc...