Can lsqnonlin regress morethan one parameter?

Hi everyone,

I have a code that regresses one parameter well (see the code below):

%% Regression Example

%% Parameters
p.KW = 10^-7;   
p.K1 = 10^-6.35;
p.K2 = 10^-10.3;

%% Initial conditions
y0.cB0 = 3.5;        
y0.cT0 = 8;          

%% Known values 
p.cNa = 2*y0.cT0;    
p.cCl = 12;

%% Measurements

t = [0 2 4 6 8 10];
Data.t = [0 2 4 6 8 10];
Data.cB_measured_1 = [3.76883357 2.143052955 -0.589243432 0.689078443 0.286944142 -0.589693715];
Data.cB_measured_2 = [3.283203989 1.397422681 2.329378469 1.642710298 -0.547381948 1.581612167];
Data.cB_measured_3 = [3.862702112 1.194583011 0.897551451 0.155508754 0.065489348 0.808999237];

Data.pH_measured_1 = [7.006085712 6.905004758 6.85705829 6.819587483 6.829728045 6.834313912];
Data.pH_measured_2 = [6.996884305 6.901179764 6.85761217 6.828627943 6.825494374 6.810138731];
Data.pH_measured_3 = [7.000879323 6.879362133 6.839654615 6.823567365 6.793112817 6.832395362];

%% Regressing using lsqnonlin, 
%[p.RegrPara1, p.RegrPara2] = lsqnonlin(@(k) funcLSQ(k, Data, p, y0), 1, 1e-4, 10);
p.RegrPara1 = lsqnonlin(@(k) funcLSQ(k, Data, p, y0), 1, 1e-4, 10);
[cB, cT, pH] = funcSimulate(Data.t, y0.cB0, y0.cT0, p);
disp(['RegrPara1 = ',num2str(p.RegrPara1)]);

%% Plotting the results

subplot(2,1,1)
plot(t, cB, t, cT, t, Data.cB_measured_1,'ko',t, Data.cB_measured_2,'ko',t, Data.cB_measured_3,'ko','LineWidth',1);
xlabel('Time'); ylabel('Concentration'); legend('C_B','C_T','Measured C_B');

subplot(2,1,2)
plot(t, pH, t, Data.pH_measured_1,'ko',t, Data.pH_measured_2,'ko',t, Data.pH_measured_3,'ko','LineWidth',1);
xlabel('Time'); ylabel('pH'); legend('pH','Measured pH');

%-------------------------------------------------------------------------------------------------------------------------------
%% Function used to calculate the error vector (which will be used by lsqnonlin)

function Err = funcLSQ(RegrPara1, Data, p, y0)

p.RegrPara1 = RegrPara1;
% p.RegrPara2 = RegrPara2;

t = [0 2 4 6 8 10];

% Using the guessed value of "k" calculate cB and pH
[cB, ~, pH] = funcSimulate(t, y0.cB0, y0.cT0, p);

% Difference between measured and predicted cB
Err_cB = [ Data.cB_measured_1 - cB; Data.cB_measured_2 - cB; Data.cB_measured_3 - cB];
       
% Difference between measured and predicted pH
Err_pH = [Data.pH_measured_1 - pH; Data.pH_measured_2 - pH; Data.pH_measured_3 - pH];

% Error vector containing both the cB and pH error.
Err = [Err_cB; Err_pH];
end

%% Calcuating state and non-state variables

function [cB, cT, pH] = funcSimulate(t, cB0, cT0, p)

[~, y] = ode23s(@(t,y) funcODE(t, y, p), t, [cB0; cT0]);
cB = y(:,1); 
cT = y(:,2);

pH = zeros(length(t), 1);
for i = 1 : length(t)
    [~, pH(i)] = funcOtherVariables(cB(i), cT(i), p);
end

end

%% ODES

function dy_dt = funcODE(t, y, p)

cB = y(1);
cT = y(2);

r = funcOtherVariables(cB, cT, p);

dcB_dt = -r;
dcT_dt = -r;

dy_dt = [dcB_dt; dcT_dt];

end

%% Other variables
function [r, pH] = funcOtherVariables(cB, cT, p)

% Determines pH by solving the charge balance
pH = fzero(@(pH) funcCharge(pH, cB, cT, p), 7);

% Determines the reaction rate once the pH (and thus cHCO3) is known
cH = 10^-pH;
cHCO3 = ( p.K1/cH / (1 + p.K1/cH + p.K1*p.K2/cH^2) )*cT;
r(1) = p.RegrPara1*cB*cHCO3;
% r(2) = p.RegrPara2*cB*cHCO3;
% r = [r(1) r(2)];
r = r(1);
end

%% Solving algebraic equation using fsolve
function z = funcCharge(pH, cB, cT, p)

cH = 10^-pH;

% Left-hand side of the charge balance
LHS = cB + p.cNa + cH;

% Right-hand side of the charge balance
RHS = p.KW/cH + ( (p.K1/cH + 2*p.K1*p.K2/cH^2) / (1 + p.K1/cH + p.K1*p.K2/cH^2) )*cT + p.cCl;

z = LHS - RHS;  % Should be equal to zero

end

However, when I introduce more paraters, it gives error message (also see code below):

%% Regression Example

%% Parameters
p.KW = 10^-7;   
p.K1 = 10^-6.35;
p.K2 = 10^-10.3;

%% Initial conditions
y0.cB0 = 3.5;        
y0.cT0 = 8;          

%% Known values 
p.cNa = 2*y0.cT0;    
p.cCl = 12;

%% Measurements

t = [0 2 4 6 8 10];
Data.t = [0 2 4 6 8 10];
Data.cB_measured_1 = [3.76883357 2.143052955 -0.589243432 0.689078443 0.286944142 -0.589693715];
Data.cB_measured_2 = [3.283203989 1.397422681 2.329378469 1.642710298 -0.547381948 1.581612167];
Data.cB_measured_3 = [3.862702112 1.194583011 0.897551451 0.155508754 0.065489348 0.808999237];

Data.pH_measured_1 = [7.006085712 6.905004758 6.85705829 6.819587483 6.829728045 6.834313912];
Data.pH_measured_2 = [6.996884305 6.901179764 6.85761217 6.828627943 6.825494374 6.810138731];
Data.pH_measured_3 = [7.000879323 6.879362133 6.839654615 6.823567365 6.793112817 6.832395362];

%% Regressing using lsqnonlin, 
[p.RegrPara1, p.RegrPara2] = lsqnonlin(@(RegrPara1, RegrPara2) funcLSQ(RegrPara1, RegrPara2, Data, p, y0), 1, 1e-4, 10);
[cB, cT, pH] = funcSimulate(Data.t, y0.cB0, y0.cT0, p);
disp(['k = ',num2str(p.k)]);

%% Plotting the results

subplot(2,1,1)
plot(t, cB, t, cT, t, Data.cB_measured_1,'ko',t, Data.cB_measured_2,'ko',t, Data.cB_measured_3,'ko','LineWidth',1);
xlabel('Time'); ylabel('Concentration'); legend('C_B','C_T','Measured C_B');

subplot(2,1,2)
plot(t, pH, t, Data.pH_measured_1,'ko',t, Data.pH_measured_2,'ko',t, Data.pH_measured_3,'ko','LineWidth',1);
xlabel('Time'); ylabel('pH'); legend('pH','Measured pH');

%-------------------------------------------------------------------------------------------------------------------------------
%% Function used to calculate the error vector (which will be used by lsqnonlin)

function Err = funcLSQ(t, RegrPara1,RegrPara2, Data, p, y0)

p.RegrPara1 = RegrPara1;
p.RegrPara2 = RegrPara1;

t = [0 2 4 6 8 10];

% Using the guessed value of "k" calculate cB and pH
[cB, cT, pH] = funcSimulate(t, y0.cB0, y0.cT0, p);

% Difference between measured and predicted cB
Err_cB = [ Data.cB_measured_1 - cB; Data.cB_measured_2 - cB; Data.cB_measured_3 - cB];
       
% Difference between measured and predicted pH
Err_pH = [Data.pH_measured_1 - pH; Data.pH_measured_2 - pH; Data.pH_measured_3 - pH];

% Error vector containing both the cB and pH error.
Err = [Err_cB; Err_pH];
end

%% Calcuating state and non-state variables

function [cB, cT, pH] = funcSimulate(t, cB0, cT0, p)

[~, y] = ode23s(@(t,y) funcODE(t, y, p), t, [cB0; cT0]);
cB = y(:,1); 
cT = y(:,2);

pH = zeros(length(t), 1);
for i = 1 : length(t)
    [~, pH(i)] = funcOtherVariables(cB(i), cT(i), p);
end

end

%% ODES

function dy_dt = funcODE(t, y, p)

cB = y(1);
cT = y(2);

r = funcOtherVariables(cB, cT, p);

dcB_dt = -r;
dcT_dt = -r;

dy_dt = [dcB_dt; dcT_dt];

end

%% Other variables
function [r, pH] = funcOtherVariables(cB, cT, p)

% Determines pH by solving the charge balance
pH = fzero(@(pH) funcCharge(pH, cB, cT, p), 7);

% Determines the reaction rate once the pH (and thus cHCO3) is known
cH = 10^-pH;
cHCO3 = ( p.K1/cH / (1 + p.K1/cH + p.K1*p.K2/cH^2) )*cT;
r(1) = p.RegrPara1*cB*cHCO3;
r(2) = p.RegrPara2*cB*cHCO3;
r = [r(1) r(2)];

end

%% Solving algebraic equation using fsolve
function z = funcCharge(pH, cB, cT, p)

cH = 10^-pH;

% Left-hand side of the charge balance
LHS = cB + p.cNa + cH;

% Right-hand side of the charge balance
RHS = p.KW/cH + ( (p.K1/cH + 2*p.K1*p.K2/cH^2) / (1 + p.K1/cH + p.K1*p.K2/cH^2) )*cT + p.cCl;

z = LHS - RHS;  % Should be equal to zero

end

The error message is :

Not enough input arguments.

Error in RegressionExample>@(RegrPara1,RegrPara2)funcLSQ(RegrPara1,RegrPara2,Data,p,y0) (line 29)
[p.RegrPara1, p.RegrPara2] = lsqnonlin(@(RegrPara1, RegrPara2) funcLSQ(RegrPara1, RegrPara2, Data, p, y0), 1, 1e-4, 10);

Error in lsqnonlin (line 218)
            initVals.F = feval(funfcn{3},xCurrent,varargin{:});

Error in RegressionExample (line 29)
[p.RegrPara1, p.RegrPara2] = lsqnonlin(@(RegrPara1, RegrPara2) funcLSQ(RegrPara1, RegrPara2, Data, p, y0), 1, 1e-4, 10);

Caused by:
    Failure in initial objective function evaluation. LSQNONLIN cannot continue.

Hence my question is, can lsqnonlin regress morethan one parameter or not?

Answer 1

[Matt J]

It is a strange usage of lsqnonlin to solve for only a single unknown parameter. For that, you would normally just use fminsearch or fminbnd. In the more usual case where there is more than a single unknown, lsqnonlin expects your objective function to accept it in the form of a vector:

objective = @(x) funcLSQ(x(1), x(2), Data, p, y0)

[Dursman Mchabe]

Thanks a lot for the answer Matt. The end goal is to use it for more than one unknown.

When I try to use:

Objective = lsqnonlin(@(RegrPara) funcLSQ(RegrPara(1), RegrPara(2), Data, p, y0), 1, 1e-4, 10);

I get the error message,

Error using RegressionExample>funcLSQ
Too many input arguments.

Error in RegressionExample>@(RegrPara)funcLSQ(RegrPara(1),RegrPara(2),Data,p,y0) (line 29)
Objective = lsqnonlin(@(RegrPara) funcLSQ(RegrPara(1), RegrPara(2), Data, p, y0), 1, 1e-4, 10);

Error in lsqnonlin (line 218)
            initVals.F = feval(funfcn{3},xCurrent,varargin{:});

Error in RegressionExample (line 29)
Objective = lsqnonlin(@(RegrPara) funcLSQ(RegrPara(1), RegrPara(2), Data, p, y0), 1, 1e-4, 10);

Caused by:
    Failure in initial objective function evaluation. LSQNONLIN cannot continue.
 
>> 

what could I be doing wrong?

[Matt J]

You have given additional arguments to lsqnonlin (1,1e-4,10) which for some reason are all scalars. They should not be.The second argument to lsqnonlin (your initial guess) also needs to be a vector containing all the unknowns. The same for any upper and lower bound (lb and ub) arguments you wish to give.

[Dursman Mchabe]

Thanks Matt, even when I try this:

Objective = lsqnonlin(@(RegrPara) funcLSQ(RegrPara(1), RegrPara(2), Data, p,[1 1]),[1, 1e-4, 10]);

The error does not go away:

Error using RegressionExample>funcLSQ
Too many input arguments.

Error in RegressionExample>@(RegrPara)funcLSQ(RegrPara(1),RegrPara(2),Data,p,[1,1]) (line 29)
Objective = lsqnonlin(@(RegrPara) funcLSQ(RegrPara(1), RegrPara(2), Data, p,[1 1]),[1, 1e-4, 10]);

Error in lsqnonlin (line 218)
            initVals.F = feval(funfcn{3},xCurrent,varargin{:});

Error in RegressionExample (line 29)
Objective = lsqnonlin(@(RegrPara) funcLSQ(RegrPara(1), RegrPara(2), Data, p,[1 1]),[1, 1e-4, 10]);

Caused by:
    Failure in initial objective function evaluation. LSQNONLIN cannot continue.
 

[Matt J]

A few problems.

First, you have given as your initial guess a 3-vector RegrPara0=[1, 1e-4, 10]. This tells lsqnonlin, among other things, that your objective function will accept 3 unknowns. However, in the expression

@(RegrPara) funcLSQ(RegrPara(1), RegrPara(2), Data, p,[1 1])

it appears that RegrPara contains only 2 unknown elements RegrPara(1) and RegrPara(2).

Additionally, you are passing 5 arguments to funcLSQ, but you have written it to expect 6:

function Err = funcLSQ(t, RegrPara1,RegrPara2, Data, p, y0)

[Dursman Mchabe]

Oops RegrPara0=[1, 1e-4, 10] is a typo...sorry.

I am not sure how many arguments must I pass to funcLSQ.

But when I try this:

function Err = funcLSQ(RegrPara, Data, p, y0)

p.RegrPara1 = RegrPara(1);
p.RegrPara2 = RegrPara(2);

I get this error message:

Error using RegressionExample>funcLSQ
Too many input arguments.

Error in RegressionExample>@(RegrPara)funcLSQ(RegrPara(1),RegrPara(2),Data,p,[1,1]) (line 29)
Objective = lsqnonlin(@(RegrPara) funcLSQ(RegrPara(1), RegrPara(2), Data, p,[1 1]),[1e-4, 10]);

Error in lsqnonlin (line 218)
            initVals.F = feval(funfcn{3},xCurrent,varargin{:});

Error in RegressionExample (line 29)
Objective = lsqnonlin(@(RegrPara) funcLSQ(RegrPara(1), RegrPara(2), Data, p,[1 1]),[1e-4, 10]);

Caused by:
    Failure in initial objective function evaluation. LSQNONLIN cannot continue.

The updated current code now looks like this:

%% Regression Example

%% Parameters
p.KW = 10^-7;   
p.K1 = 10^-6.35;
p.K2 = 10^-10.3;

%% Initial conditions
y0.cB0 = 3.5;        
y0.cT0 = 8;          

%% Known values 
p.cNa = 2*y0.cT0;    
p.cCl = 12;

%% Measurements

t = [0 2 4 6 8 10];
Data.t = [0 2 4 6 8 10];
Data.cB_measured_1 = [3.76883357 2.143052955 -0.589243432 0.689078443 0.286944142 -0.589693715];
Data.cB_measured_2 = [3.283203989 1.397422681 2.329378469 1.642710298 -0.547381948 1.581612167];
Data.cB_measured_3 = [3.862702112 1.194583011 0.897551451 0.155508754 0.065489348 0.808999237];

Data.pH_measured_1 = [7.006085712 6.905004758 6.85705829 6.819587483 6.829728045 6.834313912];
Data.pH_measured_2 = [6.996884305 6.901179764 6.85761217 6.828627943 6.825494374 6.810138731];
Data.pH_measured_3 = [7.000879323 6.879362133 6.839654615 6.823567365 6.793112817 6.832395362];

%% Regressing using lsqnonlin, 
Objective = lsqnonlin(@(RegrPara) funcLSQ(RegrPara(1), RegrPara(2), Data, p,[1 1]),[1e-4, 10]);
[cB, cT, pH] = funcSimulate(Data.t, y0.cB0, y0.cT0, p);
disp(['RegrPara = ',num2str(p.RegrPara)]);

%% Plotting the results

subplot(2,1,1)
plot(t, cB, t, cT, t, Data.cB_measured_1,'ko',t, Data.cB_measured_2,'ko',t, Data.cB_measured_3,'ko','LineWidth',1);
xlabel('Time'); ylabel('Concentration'); legend('C_B','C_T','Measured C_B');

subplot(2,1,2)
plot(t, pH, t, Data.pH_measured_1,'ko',t, Data.pH_measured_2,'ko',t, Data.pH_measured_3,'ko','LineWidth',1);
xlabel('Time'); ylabel('pH'); legend('pH','Measured pH');

%-------------------------------------------------------------------------------------------------------------------------------
%% Function used to calculate the error vector (which will be used by lsqnonlin)

function Err = funcLSQ(RegrPara, Data, p, y0)

p.RegrPara1 = RegrPara(1);
p.RegrPara2 = RegrPara(2);

t = [0 2 4 6 8 10];

% Using the guessed value of "k" calculate cB and pH
[cB, cT, pH] = funcSimulate(t, y0.cB0, y0.cT0, p);

% Difference between measured and predicted cB
Err_cB = [ Data.cB_measured_1 - cB; Data.cB_measured_2 - cB; Data.cB_measured_3 - cB];
       
% Difference between measured and predicted pH
Err_pH = [Data.pH_measured_1 - pH; Data.pH_measured_2 - pH; Data.pH_measured_3 - pH];

% Error vector containing both the cB and pH error.
Err = [Err_cB; Err_pH];
end

%% Calcuating state and non-state variables

function [cB, cT, pH] = funcSimulate(t, cB0, cT0, p)

[~, y] = ode23s(@(t,y) funcODE(t, y, p), t, [cB0; cT0]);
cB = y(:,1); 
cT = y(:,2);

pH = zeros(length(t), 1);
for i = 1 : length(t)
    [~, pH(i)] = funcOtherVariables(cB(i), cT(i), p);
end

end

%% ODES

function dy_dt = funcODE(t, y, p)

cB = y(1);
cT = y(2);

r = funcOtherVariables(cB, cT, p);

dcB_dt = -r;
dcT_dt = -r;

dy_dt = [dcB_dt; dcT_dt];

end

%% Other variables
function [r, pH] = funcOtherVariables(cB, cT, p)

% Determines pH by solving the charge balance
pH = fzero(@(pH) funcCharge(pH, cB, cT, p), 7);

% Determines the reaction rate once the pH (and thus cHCO3) is known
cH = 10^-pH;
cHCO3 = ( p.K1/cH / (1 + p.K1/cH + p.K1*p.K2/cH^2) )*cT;
r(1) = p.RegrPara1*cB*cHCO3;
r(2) = p.RegrPara2*cB*cHCO3;
r = [r(1) r(2)];

end

%% Solving algebraic equation using fsolve
function z = funcCharge(pH, cB, cT, p)

cH = 10^-pH;

% Left-hand side of the charge balance
LHS = cB + p.cNa + cH;

% Right-hand side of the charge balance
RHS = p.KW/cH + ( (p.K1/cH + 2*p.K1*p.K2/cH^2) / (1 + p.K1/cH + p.K1*p.K2/cH^2) )*cT + p.cCl;

z = LHS - RHS;  % Should be equal to zero

end

[Matt J]

Now you have written funcLSQ to accept 4 arguments but you are passing in 5...

[Dursman Mchabe]

I think I am now more confused. I will go through the code thoroughly.

[Matt J]

Well, before you do anything with lsqnonlin, I recommend that you first make sure that funcLSQ is working. When I simply call it with your initial point RegPara=[1e-4,10],

funcLSQ([1e-4,10], Data, p,y0);

I obtain the following errors

Error using odearguments (line 93)
@(T,Y)FUNCODE(T,Y,P) must return a column vector.

Error in ode23s (line 121)
    = odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in test>funcSimulate (line 59)
[~, y] = ode23s(@(t,y) funcODE(t, y, p), t, [cB0; cT0]);

Error in test>funcLSQ (line 43)
[cB, cT, pH] = funcSimulate(t, y0.cB0, y0.cT0, p);

Error in test (line 30)
funcLSQ([1e-4,10], Data, p,y0);

[Dursman Mchabe]

I see, thanks a lot Matt. I only understand the first error. What can I do to make @(T,Y)FUNCODE(T,Y,P) to return a column vector?

[Dursman Mchabe]

Do you think it will be a good idea to abondon the lsqnonlin approach, and go for the lsqcurvefit ?

[Matt J]

The differences between lsqnonlin and lsqcurvefit are very minor. Basically, lsqcurvefit is more convenient if you want to plot a fitted curve after the optimization.