Trying to supply fsolve with system's Jacobian but receiving error ''check for incorrect argument data type"

Hello,

I am trying to solve a system of approximately 600 non-linear equations in MATLAB with the help of fsolve. While I managed to get the programme I am writing to provide a solution to said system, I am currently attempting to supply the system's Jacobian to fsolve, in order to improve efficiency as well as accuracy. However, so far I have been prevented from doing so by the following error:

Check for incorrect argument data type or missing argument in call to function 'jacobian'.

Error in fun_name (line 509)
    JAC = jacobian(xss,x);

Error in
File_name>@(x)fun_name(x,inputs) (line 712)
    sol = @(x)fun_name(x,inputs)

Error in fsolve (line 269)
            [fuser,JAC] = feval(funfcn{3},x,varargin{:});

Error in File_name(line 731)
[xss(ITER,:),x(ITER,:),exitflag,output,jacobian] = fsolve(sol,x0,options);

Here is the relevant (simplified) code I used:

clear;
close all;

%%% CALIBRATION, ETC. %%%

options = optimoptions('fsolve','Display','iter','Algorithm','Levenberg-Marquardt','MaxFunEval',1.0e+08,...
    'SpecifyObjectiveGradient',true,'Jacobian','on');

maxiter = 3; % could be anything

%%% Section initialising endogenous variables and any variables that become
%%% updated so that solution values from each iteration can be saved

for ITER=1:maxiter
    % Vector of initialised variables
    if ITER == 1
        x0 = [list of initialised endogenous variables]; % about 600
    elseif ITER > 1
        x0 = xss(ITER-1,:); % Use the current solution as initial guess for the next period
    end
    
    sol = @(x)fun_name(x,inputs); % function handle in separate file [this is the (line 712) mentioned in the error]

    a = sym('EV',[1 593]); % EV = Endogenous Variables; create a symbolic vector
    FUN = sol(a); 
    [xss(ITER,:),x(ITER,:),exitflag,output,jacobian] = fsolve(sol,x0,options); % [This is the (line 731) mentioned in the error]

    % Some more (unrelated) code
    
end

The function containing the equations and the udpate mechanism is then as follows:

function [xss,JAC] = fun_name(x,inputs) 
% A. Endogenous variables. For instance:
A = x(1:a_given_length);                CT = a_given_length;
B = x(CT+1:CT+another_length);          CT = CT+another_length;
.
.
.
%%% B. Section with updating procedures for certain variables %%%

% Reshape any matrices that have been vectorised at beginning of function

%%% C. Equations. For instance:

xss1 = A.*B; xss{1} = xss1;
.
.
.
xss = horzcat(xss{:});

% Supplying the system's Jacobian
if nargout(@fun_name) > 1
    JAC = jacobian(xss,x); % This is the (line 509) mentioned in the error
end

% some more code, including e.g. 
assignin('base','A',A);
% to save updated values from each iteration

end

Having tried to debug the programme, I think I know what the issue is: while the function 'jacobian' requires symbolic inputs, due to some mistake I am making in the code the function ends up being called at JAC with 'double'-type data as inputs. I checked this by running

class(x)

and

class(xss)

before it, and I saw both are 'double'. Having placed a breakpoint at JAC = jacobian(xss,x); the first time it is reached, however, I also noticed that my attempt at creating symbolic variables through a and FUN does create a symbolic JAC matrix of the required dimensions, but when fsolve comes around both xss and x are already both of type 'double' again. How can I get fsolve to solve with its Jacobian in this settings? And, on a related note, is it worth doing so in my case? On my laptop it takes approximately 3 minutes to solve three iterations, so if it takes a long time to create the symbolic variables then I might not see enough gains anyway. And would there be a 'visible' improvement in terms of solution accuracy, in your opinion? Thank you all for the help!

Answer 1

[Matt J]

It is not clear (from what you posted), why you are using symbolic variable manipulation at all. fsolve is a numerical solver, not a symbolic solver.

If you need the Symbolic Math Toolbox to derive the analytical expression for some part of the Jacobian calculation, that's fine. However, you should not be repeating the symbolic analysis in every iteration of fsolve. The normal and efficient workflow is to derive any symbolic formulas that you need once and once only, and to do it prior to the optimization. Before running hte optimization, you would use matlabFunction() to convert your symbolic functions to an ordinary numerical function which the objective function code fun_name() can then use.

[Mark]

Thank you for your answer. I am not looking to use symbolic variables if they are not needed, but I thought they were required as data type for the arguments of the 'jacobian' function: the programme does not throw the error above when xss and x are symbolic, only when they are double. Does this mean that I need to move the section computing the Jacobian from within the function to the main file, above the for loop? But I am defining the function handle within the loop as the iteration number itself is an input to fun_name() -- so how would I get it to recognize the variables as required before the loop?

[Matt J]

but I thought they were required as data type for the arguments of the 'jacobian' function:

But you don't need to use the jacobian() function if you know how to write code that calculates the Jacobian of fun_name at a given x already.

If you don't know how to write such code, the jacobian() command from the Symbolic Math Toolbox can be helpful, but jacobian() is merely a code-writing tool. You can use it get the expression for the derivative of something when you don't know how to perform the necessary calculus yourself, but that is all. Once you have the symbolic expression for a derivative, you still have to put it in a form that fsolve can use, i.e., a code routine that will take in numbers, apply the expression for the derivative, and spit out numbers. However, translating a symbolic expression into an actual numeric function is something that matlabFunction() can help automate. Thus, as a simplified example, you could do something like this:

syms x


for ITER=1:3
    
    funcSYM=(x-ITER)^2;       %a symbolic expresssion, depending on ITER
    jacSYM=jacobian(funcSYM); %the expression's symbolic Jacobian
    
    inputs.func=matlabFunction(funcSYM); %Convert to actual numeric functions
    inputs.jac=matlabFunction(jacSYM)
    
    sol=@(x) fun_name(x,inputs);

end

inputs = struct with fields:
    func: @(x)(x-1.0).^2
     jac: @(x)x.*2.0-2.0

inputs = struct with fields:
    func: @(x)(x-2.0).^2
     jac: @(x)x.*2.0-4.0

inputs = struct with fields:
    func: @(x)(x-3.0).^2
     jac: @(x)x.*2.0-6.0

Now, fun_name has access to both inputs.func() and its derivative inputs.jac in its calculations.

Does this mean that I need to move the section computing the Jacobian from within the function to the main file, above the for loop?

Not above the for-loop. Inside the for-loop, but outside of fun_name, like in my example above.

[Mark]

Thank you for the suggestion! To make sure I understand: do I need to place the system's equations where your funcSYM is? If that is the case, when I tried doing so (see code below) the conversion from sym to numeric through matlabFunction took a very long time (1 hr 4 min on my laptop). This could be because I am working with approximately 600 endogenous variables, which is a lot. As far as I can tell, however, writing the Jacobian manually is not really a feasible option in my case, given the huge number of both equations and variables... Is there a way to shorten the conversion time with the method above, or to supply the Jacobian in some other way? My goal is to speed up the solver by providing the Jacobian, so of course a conversion time of 1 hr 4 min would not help there. Or was I wrong in placing all the equations where funcSYM is? Schematically, what I did is this:

syms x
for ITER=1:3
    sol=@(x) fun_name(x,inputs);
    funcSYM = sol(x);
    jacSYM = jacobian(funcSYM,x);
    JAC = matlabFunction(JAC); % This is the line which took 1 hr 4 min
end

As you can see I placed the function handle sol above the symbolic function and the Jacobian, because that way I could call sol to provide me with the model's equations directly. Thank you for the help!

[Matt J]

Nothing can be said really without seeing the full contents of fun_name.

[Mark]

Ah -- I see. Unfortunately, I cannot share the full code, as it is for an ongoing project, but I will try to reproduce the main architecture here:

function [xss] = fun_name(x, list of inputs: e.g. a, b, c, etc.)
                                                       
% A. Endogenous variables

A = x(1:length_of_A);                   CT = length_of_A;
B = x(CT+1:CT+length_of_B);             CT = CT+length_of_B;
.
.   % List of all endogenous variables with their vectorised size: ca. 600 variables,
.   % grouped into 53 vectors (e.g. A, B...)
.
X = x(CT+1:CT+length_of_X);             CT = CT+length_of_X;
Y = x(CT+1:CT+length_of_Y);

% B. Specify updating mechanism, as the function is being called within a for-loop
if ITER > 1
% Sample updating
    a_update = a(ITER-1)*(1-b)+N(ITER-1,c~=0); % Some inputs, e.g. c, are needed for
    % logical indexing
    a(ITER,:) = a_update;
end

% C. Some shocks
if ITER == 2    
    m(ITER) = m(ITER-1)+210;
    % other shocks like the one above, at different time periods
end

% D. EQUATIONS

% i. Reshape any vectorised matrices 
P = reshape(P,[rows_of_P columns_of_P]);
% etc.

% ii. Model equations

xss1 = N(ITER,c~=0).*C.*(1+a(ITER,:)).*(1+b); xss{1} = xss1;
.
.   % Repeat the same procedure for all equations (I have 59 such equations, often of
.   % varying length)
.
xss59 = W-sum(N(ITER,c~=0).*L,2).'-R; xss{59} = xss59;

xss = horzcat(xss{:});

% Report any exogenous variables that have been updated back to the main
% workspace, so that they can be saved and used in the next iteration
assignin('base','a',a);
assignin('base','m',m);
% etc.

end

This is the code I have, but let me know if you would like me to be more specific (and apologies for any inconsistencies in the code when it comes to e.g. naming, as I had to simplify and omit irrelevant information).

Let me know if this is enough for you to be able to help me figure out how to supply the Jacobian correctly, using symbolic variables if needed. Thank you for your time!

[Matt J]

There's nothing really solid that I can say while only seeing part of the code. From what is shown, it looks like all the operation done with yor unknown variables are linear. If so, then this is really a system of linear equations. Instead of fsolve, you should really be putting it in matrix-vector form A*x=b and solving it as x=A\b.

[Mark]

I see -- as I said, unfortunately I cannot share more of the code. Thank you anyways: your suggestions were useful in making the programme run, at least!