How do I use the Levenberg-Marquardt algorithm for lsqcurvefit?

I'm trying to run a fit on some data using lsqcurvefit. Here is the basic structure of the code:

opts = optimset('Display', 'off'); 
params0 = [0; 0.60e-4; 1.30e-3];
lb = [0, 0.40e-4, 0.92e-4];
ub = [1, 0.80e-4, 1.76e-4];

[params, sigval] = lsqcurvefit(@sm_signal, params0, double(b_unique),...
                    S_avg_formatted, lb, ub, opts);

function [signal] = sm_signal(params,bvalues)
% Assign parameters from params for use within the function sm_signal.
[f, perp, delta_lambda] = deal(params(1),params(2),params(3));
D = 3.00e-3; 
v = 0.01;

signal = (1-f)*exp(-bvalues*D) + ...
        f*(exp(-bvalues*perp).*sqrt(pi./(4*bvalues*(delta_lambda))).*erf(sqrt(bvalues*(delta_lambda)))) + ...
        v*(perp/delta_lambda);
end

In general, this works well. However, I have a specific data set that I want to fit which has fewer b-values than parameters to be fit. That is to say, size(b_unique) = [2 1] and size(params0) = [3 1]. Trying to fit the data gives the following error:

Error using lsqncommon (line 64)
The Levenberg-Marquardt algorithm does not handle bound constraints and the trust-region-reflective algorithm requires at least as many equations as
variables; aborting.

Error in lsqcurvefit (line 271)
    lsqncommon(funfcn,xCurrent,lb,ub,options,defaultopt,allDefaultOpts,caller,...

I've been reading the documentation and it says that for an underdetermined system, the Levenberg-Marquardt algorithm could be used instead. But whenever I run the code using that algorithm, I receive the same error. What am I doing wrong? Is it the fact that I am imposing bounds on my parameters with lb and ub? If so, how do I fix it?

Thanks,

Warren

[Torsten]

What sense does it make to fit a model with more parameters than measurement data ? None.

[Warren Boschen]

I agree with you, but unfortunately this is the most readily available data at the moment. I'm hopeful that the documentation specifically cites the Levenberg-Marquardt algorithm as a way to handle such a system, but I can't get it to work in the first place.

Do you think there is any way that I can use this data? Or would it require a complete reformulation of the problem?

[Torsten]

It simply doesn't make sense what you try to do. The parameters you will find are almost arbitrary and will only be useful to reproduce the two measurement data. Transfering the model to different conditions - which is the reason to build a mathematical model - will not be possible.

Answer 1

[Matt J]

If the signal is supposed to be smooth, it would make sense to add some sort of roughness penalty term to the reisduals, e.g.,

opts = optimoptions('lsqnonlin','Display', 'off','Algorithm','levenberg-Marquardt'); 
params0 = [0; 0.60e-4; 1.30e-3];
lb = [0, 0.40e-4, 0.92e-4];
ub = [1, 0.80e-4, 1.76e-4];


bvalues=double(b_unique);
fun=@(params) sm_signal(params,bvalues,S_avg_formatted);

[params, sigval] = lsqnonlin(@sm_signal, params0, lb, ub, opts);

function residuals = sm_signal(params,bvalues,S)
% Assign parameters from params for use within the function sm_signal.
[f, perp, delta_lambda] = deal(params(1),params(2),params(3));
D = 3.00e-3; 
v = 0.01;

 signal = (1-f)*exp(-bvalues*D) + ...
         f*(exp(-bvalues*perp).*sqrt(pi./(4*bvalues*(delta_lambda))).*erf(sqrt(bvalues*(delta_lambda)))) + ...
         v*(perp/delta_lambda);

 penalty=0.01*diff(signal);
 
 
 residuals=[signal(:)-S(:) ; penalty(:)];



end

Answer 2

[Matt J]

Here's a silly solution to the problem. The fit will be meaningless, however.

opts = optimset('Display', 'off'); 
params0 = [0; 0.60e-4; 1.30e-3];
lb = [0, 0.40e-4, 0.92e-4];
ub = [1, 0.80e-4, 1.76e-4];


fun=@(p,b) repmat(sm_signal(p,b),1,2);

[params, sigval] = lsqcurvefit(@sm_signal, params0, double(b_unique),...
                    [S_avg_formatted, S_avg_formatted], lb, ub, opts);

[Warren Boschen]

Meaningless in what sense?

[Matt J]

In the sense that the solution is still under-determined. So, you won't be able to make predictions with it, see also Torsten's comment.