Nonlinear fitting: how do I split the linear and the nonlinear problems?

I am fitting a function to some data I simulated. I managed to get intelligent constraints that help the fit quite a bit, even with a lot of noise.

This is the function and as you can see, c(1) and c(2) are linear, while lam(1), lam(2), lam(3) and lam(4) are nonlinear. I am following the procedure explained here (https://it.mathworks.com/help/optim/ug/nonlinear-data-fitting-problem-based-example.html#NonlinearDataFittingProblemBasedExample-4) to split linear and nonlinear parameters.

% Create a function that computes the value of the response at times t when the parameters are c and lam
diffun = ((c(1)) .* ((1 - exp(-t / lam(1))) .* exp(-t / lam(2))) * (Vm - (-70)) + ...
((c(2)) .* ((1 - exp(-t / lam(3))) .* exp(-t / lam(4))) * -30));

This is the code that I came up with, but for some reason it's not working. To generate the data:

function [EPSC, IPSC, CPSC, t] = generate_current(G_max_chl, G_max_glu, EGlu, EChl, Vm, tau_rise_In, tau_decay_In, tau_rise_Ex, tau_decay_Ex,tmax)
    
    dt = 0.1;                               % time step duration (ms)
    t = 0:dt:tmax-dt;
    
    % Compute compound current
    
    IPSC = ((G_max_chl) .* ((1 - exp(-t / tau_rise_In)) .* exp(-t / tau_decay_In)) * (Vm - EChl));
    
    EPSC = ((G_max_glu) .* ((1 - exp(-t / tau_rise_Ex)) .* exp(-t / tau_decay_Ex)) * (Vm - EGlu));
    
    CPSC = IPSC + EPSC;
    
end

To fit the function:

% Simulated data
[EPSC,IPSC,CPSC,t] = generate_current(80,15,0,-70,-30,0.44,15,0.73,3,120);
ydata = awgn(CPSC,25,'measured'); % Add white noise
% Values
Vm = -30;
% Initial values for fitting
gmc = 40; gmg = 20; tde = 0.2; tdi = 8; tre = 1.56; tri = 3;
% Objective function
c = optimvar('c',2); % Linear parameters
lam = optimvar('lam',4); % Nonlinear parameters
% Bounds
c.LowerBound = [0, 0];
c.UpperBound = [200, 200];
lam.LowerBound = [0.16,7.4,1.1,2.6];
lam.UpperBound = [0.29,8.4,2.3,3.3];
x0.c = [gmc,gmg]; % Starting values
x0.lam = [tri,tdi,tre,tde]; % Starting values
% Create a function that computes the value of the response at times t when the parameters are c and lam
diffun = ((c(1)) .* ((1 - exp(-t / lam(1))) .* exp(-t / lam(2))) * (Vm - (-70)) + ...
((c(2)) .* ((1 - exp(-t / lam(3))) .* exp(-t / lam(4))) * -30));
%Solve the problem using solve starting from initial point x02
x02.lam = x0.lam;
%To do so, first convert the fitvector function to an optimization expression using fcn2optimexpr.
F2 = fcn2optimexpr(@(x) fitvector(x,t,ydata),lam,'OutputSize',[length(t),1]);
% Create a new optimization problem with objective as the sum of squared differences between the converted fitvector function and the data y
ssqprob2 = optimproblem('Objective',sum((F2' - ydata).^2));
[sol2,fval2,exitflag2,output2] = solve(ssqprob2,x02)
% Plot
resp = evaluate(diffun,sol2);
hold on
plot(t,resp)
hold off

The error is:

Solving problem using lsqnonlin.
Error using optim.problemdef.OptimizationProblem/solve
Matrix dimensions must agree.
Error in SplittFit (line 37)
[sol2,fval2,exitflag2,output2] = solve(ssqprob2,x02)
Caused by:
    Failure in initial objective function evaluation. LSQNONLIN cannot continue.

Not sure why.

Answer 1

[Alan Weiss]

You should follow the example more closely. In the example the lambda variables only are declared to be optimization variables; the c variables are not, and are computed by backslash for given values of the lambda variables in the fitvector function.

For your case you will need to update the fitvector function from the example to handle a 4-D lambda vector that differs from your 4-D function because you have some (1-exp(-t/lambda)) terms, not just exp(-t/lambda). You need to write out the linear equations in c and solve those equations in the fitvector function.

Good luck,

Alan Weiss

MATLAB mathematical toolbox documentation

[Samuele Bolotta]

Hi Alan, thanks for the valuable input! I modified this way the code and the fitvector function. I think now almost everything is working correctly. However, I have an error that is actually present also in the mathworks example.

% Simulated data
[EPSC,IPSC,CPSC,t] = generate_current(80,15,0,-70,-30,0.44,15,0.73,3,120);
ydata = awgn(CPSC,25,'measured'); % Add white noise
ydata = ydata';
t = t';

% Values
Vm = -30; DF1 = Vm +70; DF2 = Vm;

% Initial values for fitting
gmc = 40; gmg = 20; tde = 0.2; tdi = 8; tre = 1.56; tri = 3;

% Objective function
c = optimvar('c',2); % Linear parameters
lam = optimvar('lam',4); % Nonlinear parameters

% Bounds
c.LowerBound = [0, 0];
c.UpperBound = [200, 200];
lam.LowerBound = [0.16,7.4,1.1,2.6];
lam.UpperBound = [0.29,8.4,2.3,3.3];

x0.c = [gmc,gmg]; % Starting values
x0.lam = [tri,tdi,tre,tde]; % Starting values

% Create a function that computes the value of the response at times t when the parameters are c and lam
temp1 = 1 - exp(-t / lam(1));
temp2 = exp(-t / lam(2));
temp3 = 1 - exp(-t / lam(3));
temp4 = exp(-t / lam(4));

diffun = (temp1 .* temp2) * DF1 * c(1) + (temp3 .* temp4) * DF2 * c(2);

%Solve the problem using solve starting from initial point x02
%To do so, first convert the fitvector function to an optimization expression using fcn2optimexpr.
F2 = fcn2optimexpr(@(x) fitvector(x,t,ydata),lam,'OutputSize',[length(t),1]);

% Create a new optimization problem with objective as the sum of squared differences between the converted fitvector function and the data y
ssqprob = optimproblem('Objective',sum((F2 - ydata).^2));

[sol,fval,exitflag,output] = solve(ssqprob,x0)

% Plot
resp = evaluate(diffun,sol);
hold on
plot(t,resp)
hold off

And

function yEst = fitvector(lam,xdata,ydata)

temp1 = 1 - exp(-xdata / lam(1));
temp2 = exp(-xdata / lam(2));
temp3 = 1 - exp(-xdata / lam(3));
temp4 = exp(-xdata / lam(4));

A(:,1) = temp1 .* temp2 * 40;
A(:,2) = temp3 .* temp4 * -30;

c = A\ydata; % solve A*c = y for linear parameters c
yEst = A*c; % return the estimated response based on c

However, it still gives me this error:

Error using optim.problemdef.OptimizationExpression/evaluate
Second argument must contain values for variable 'c'.

Error in SplittFit (line 44)
resp = evaluate(diffun,sol);

Looks like the values for c are not being stored correctly.

[Alan Weiss]

I don't see the code for SplittFit. So I cannot see where there is a problem. I can assure you that the "MathWorks example" runs without error. Is it possible that you are running an old version of MATLAB, one that does not support least-squares fitting in the problem-based approach?

Alan Weiss

MATLAB mathematical toolbox documentation

[Yanxin Liu]

I was also studying the Matlab example you pasted link. From my understanding, the equation you want to fit to must have the form y = c(1)*exp(-lam(1)*t) + ... + c(n)*exp(-lam(n)*t). That's say, each columne of matrix A can only represent one nonlinear parameter (lambda). Then the linear parameters c can be derived from backslash by trying different pairs of lambdas. In your problem, the lam1, lam2, also lam3, lam4 are entangled together.

[Alan Weiss]

The response function can have any form. For complicated examples, see Fit ODE, Problem-Based and Fit an Ordinary Differential Equation (ODE), which have objective functions given by the solution of ODEs.

By no means are you required to split the response into a linear and a nonlinear part; when you can, then you can use the trick of fitting the linear part after fitting the nonlinear part.

Alan Weiss

MATLAB mathematical toolbox documentation