checkGradients, but the objective function has two inputs: x and xdata?

I'm using lsqcurvefit with the following objective function and Jacobian:

function [f, jacF] = semiCircle(p, Q)

P0 = p(1);
Q0 = p(2);
r = p(3);

f = P0 + sqrt(r^2 - (Q-Q0).^2);

if nargout > 1 % need Jacobian
    jacF = [1, (Q-Q0)./sqrt(r^2-(Q0-Q).^2), r./sqrt(r^2-(Q0-Q).^2)];
end

end

I'd like to use checkGradients to verify if the Jacobian is correct. However, all of the examples in the documentation just have objective functions with one input, the parameters 'x'. Whereas my function semiCircle has two inputs: the parameters 'p' and the xdata 'Q'. Is there a way to use checkGradients for such a function?

Answer 1

[Torsten]

valid = checkGradients(@(p)semiCircle(p, Q),p0)

[Benjamin Pepper]

Many thanks, this seems to work.

From a seperate mistake, I had to modify the Jacobian a bit because there was an error with concatenating the arrays:

function [f, jacF] = semiCircle(p, Q)

P0 = p(1);
Q0 = p(2);
r = p(3);

f = P0 + sqrt(r^2 - (Q-Q0).^2);

if nargout > 1 % need Jacobian
    jacF = zeros(length(Q), length(p));
    
    for i = 1:length(Q)
        jacF(i,:) = [1, (Q(i)-Q0)/sqrt(r^2-(Q0-Q(i))^2), r/sqrt(r^2-(Q0-Q(i))^2)];
    end
    
end

This runs and the checkGradients returns logical 1 (true).

However, the code was running before even though the Jacobian was wrong. Also, the solution lsqcurve gives is the same whether or not 'SpecifyObjectiveGradient' is set to 'true'. This suggests that its somehow ignoring whether the gradient is specified or not. Do you have any insights regarding this?

[Torsten]

If you specified a wrong Jacobian in the sense that you concatenated the arrays incorrectly, MATLAB should have thrown an error.

If you specified a wrong Jacobian in the sense that some elements were not computed correctly, it may happen that nonetheless the code is running and converging.

And it should be the case that whether you specify 'SpecifyObjectiveGradient' to be true or not, the solution MATLAB gives should be the same. What else do you expect in this case ?

[Benjamin Pepper]

When I concatenated the arrays incorrectly, MATLAB did not throw an error for lsqcurvefit. However, it did throw an error for checkGradients.

I expected that the answers would be similar, but maybe some slight differences. For fmincon, I tried with and without a Jacobian, and the answers were different in some cases (specifying the gradient tended to improve the solution).

[Torsten]

When I concatenated the arrays incorrectly, MATLAB did not throw an error for lsqcurvefit. However, it did throw an error for checkGradients.

Does it throw an error if you don't check the gradients, but set 'SpecifyObjectiveGradient' to true ?

If it doesn't, lsqcurvefit does not recognize your analytical gradients.

[Benjamin Pepper]

Does it throw an error if you don't check the gradients, but set 'SpecifyObjectiveGradient' to true ?

No, it seems it does not throw an error for this.

[Torsten]

That's strange and seems to indicate that "lsqcurvefit" does not use your analytical gradients. Can you supply the full code to test ?

[Benjamin Pepper]

Here's the code and attachments:

load('data.mat')

Vb_ll_rms = 690;

% Inital guess
p10 = 5e6;
p20 = 2.5e7;
p30 = 3e7;

p0 = [p10, p20, p30];

[Rls, Xls, Vgls, gradientCheck] = lsqcurvefitNLS(p0, Q, P, Vb_ll_rms)

Rls = 0.0024

Xls = 0.0238

Vgls = 689.9365

gradientCheck = logical
   1

lsqcurvefirNLS calls lsqcurvefit with box and linear constraints, which in turn calls the semiCircle function. Rls, Xls, and Vgls are the estimated parameters, and gradientCheck verifies the objective gradient given.

Here, the code runs fine. However, in my desktop MATLAB it returns the following error:

Error using  - 
Too many output arguments.

Error in lsqcurvefit>@(X)feval(FUN,X,XDATA)-YDATA (line 237)
[xCurrent,Resnorm,FVAL,EXITFLAG,OUTPUT,LAMBDA,JACOB] = cnls(caller, @(X) feval(FUN, X, XDATA) - YDATA,...
    
Error in cnls/funExt (line 355)
[Fvec,JACOBIAN] = feval(FUN, Xk);                    

Error in cnls/objective (line 319)
[f,F] = funExt(xIn);

Error in cnls>@(x)objective(x,true) (line 188)
fungrad = @(x) objective(x, true);

Error in fmincon (line 598)
initVals.g = feval(funfcn{4},X,varargin{:});

Error in cnls (line 215)
[X,~,EXITFLAG,OUTPUT,LAMBDA] = fmincon({fun,fungrad},X(:),...
    
Error in lsqcurvefit (line 237)
[xCurrent,Resnorm,FVAL,EXITFLAG,OUTPUT,LAMBDA,JACOB] = cnls(caller, @(X) feval(FUN, X, XDATA) - YDATA,...
    
Error in lsqcurvefitNLS (line 19)
p = lsqcurvefit(@semiCircle, p0, Q, P, lb, ub, A, b, [], [], [], options);

Caused by:
Failure in initial objective gradient function evaluation. FMINCON cannot continue.
lsqcurvefit cannot continue.

This might be because there are square-root terms in the Jacobian/objective gradient which are maybe giving complex values:

jacF(i,:) = [1, (Q(i)-Q0)/sqrt(r^2-(Q0-Q(i))^2), r/sqrt(r^2-(Q0-Q(i))^2)]

However, the linear contraints are designed to prevent the solution from entering that region. And, the code seems to run here online without any problems.

I could try updating my MATLAB to the 2024b version?

[Torsten]

Could you include the .m-files instead of the .mlx-files ?

[Benjamin Pepper]

I reached my daily uploads limit, so I'll just put the functions here:

function [R, X, Vg, gradientCheck] = lsqcurvefitNLS(p0, Q, P, Vo)

% Box constraints
p1_ub = min(P);
p2_lb = max(Q);
p3_lb = max(Q) - min(Q);

lb = [0, p2_lb, p3_lb];
ub = [p1_ub, inf, inf];

% Linear constraints
A = [0, 1, -1];
b = min(Q); 

gradientCheck = checkGradients(@(p)semiCircle(p,Q),p0);

options = optimoptions('lsqcurvefit','Display','off','SpecifyObjectiveGradient',true);

p = lsqcurvefit(@semiCircle, p0, Q, P, lb, ub, A, b, [], [], [], options);

P0 = p(1);
Q0 = p(2);
r = p(3);

R = P0/(P0^2 + Q0^2)*Vo^2;
X = Q0/(P0^2 + Q0^2)*Vo^2;
Vg = sqrt(r^2/(P0^2 + Q0^2)*Vo^2);

And the other one:

function [f, jacF] = semiCircle(p, Q)

P0 = p(1);
Q0 = p(2);
r = p(3);

f = P0 + sqrt(r^2 - (Q-Q0).^2);

if nargout > 1 % need Jacobian
    jacF = zeros(length(Q), length(p));

    for i = 1:length(Q)
        jacF(i,:) = [1, (Q(i)-Q0)/sqrt(r^2-(Q0-Q(i))^2), r/sqrt(r^2-(Q0-Q(i))^2)];
    end

end