fsove is extremely slow

Long story short I am working on a problem where I have 35 equations and 35 unknowns where each unknown can be up to power 3 (e.g. f1=x1^3+x1^2*x2+...+x35^3, ... f35=...). I have tried to use both fsolve and lsqnonlin to solve solve this system, but the run time is extremely long! I'm talking only 2 steps within fsolve over the course of an hour. Is this time frame to be expected for a system like this, or am I implementing something poorly? The basiscs of my implementation are as follows:

xhat_0s=sym('x_0s',[numStateVars,1],'real'); %numStateVars=35
xhat_0s_1=[xhat_0s;1];
xhat_0s_1Mat=xhat_0s_1*xhat_0s_1';

%Code to calculate the arrays Lambda_sr (35x35x36x36) and N_sr (35x36) is
%not shown but takes place here, Lambda_sr and N_sr are dense and 4-D and
%2-D doubles respectively

expandxmat=permute(repmat(xhat_0s_1Mat,[1,1,35,35]),[3,4,1,2]);
finalLambda=sum(sum(Lambda_sr.*expandxmat,4),3);

%An equivalent but slower way to calculate finalLambda is:
        % finalLambda=zeros(numStateVars,numStateVars);
        % for ind1=1:numStateVars+1
        %     for ind2=1:numStateVars+1
        %         finalLambda=finalLambda+Lambda_sr(:,:,ind1,ind2)*xhat_0s_1Mat(ind1,ind2);
        %     end
        % end

finalN=N_sr*xhat_0s_1;

funct=finalLambda*xhat_0s-finalN;
bigfunct=matlabFunction(funct,'Vars',{xhat_0s}); %this might be the bottleneck
residualFunct=@(x) bigfunct(x); %or maybe this is the bottleneck
initialGuess=zeros(numStateVars,1);

options=optimoptions('fsolve','Display','iter',...
    'FunctionTolerance',1e-12,'StepTolerance',1e-16,...
    'OptimalityTolerance',1e-12,'MaxIterations',50,...
    'MaxFunctionEvaluations',1000);
[xhat_0_fsolve,fval,exitFlag_fsolve,~]=fsolve(residualFunct,initialGuess,options)

My hunch is that the bottleneck is how I am setting up the equation residualFunct to be used by fsolve, but I really don't know. Maybe this problem isn't feasible but I feel like it should be doable. I was able to run the same code for 15 unknowns and it took a while to run, but was still doable (maybe 30 minutes for fsolve to find a solution). Any help here would be much appreciated as this is the last test case that I am attempting to run for my thesis.

[Walter Roberson]

For testing, we need several more variables, such as Lambda_sr and N_sr

[dpb]

"For testing, we need several more variables, such as Lambda_sr and N_sr"

"Code to calculate the arrays Lambda_sr (35x35x36x36) and N_sr (35x36) is not shown here..."

Also would have to have the functional code as well.

I would also guess that mixing in 35 symbolic variables with arrays of this size and fsolve is a lethal combination for compute time.

It would also be better to describe the problem from first principles rather than asking folks to debug complex code that isn't working well...the approach may end up the same or similar, but generally speaking, the way to make orders of magnitude improvement in code performance is by modifying the algorithms rather than optimization of existing shown to be slow code---although my first approach would be to remove the Symbolic TB from the problem and go with a purely numerical solution -- after, of course, ensuring the correctness of the system function code first.

[Walter Roberson]

My guess is that you would see performance improvements if you used

bigfunct = matlabFunction(funct, 'file', 'bigfunct.m', 'Vars',{xhat_0s}, 'optimize', true);

However I am uncertain whether optimize true is reliable in R2023b; in older versions it was not reliable.

[Jacob]

Hello, thank you for the responses and apologies for not including the code to calculate Lambda_sr and N_sr. These variables are calulated as a result of many operations and function calls, and so I do not think it would be feasible to include their entire calculation here, but I will try to provide a summary of what is happening.

The overarching code is an orbit determination script that performs a Batch Least Squares Estimation while including second order terms in the state transition matrix (a state transition tensor). The state transition matrix (Phi) and state transition tensor (Psi) are propagated forward in time to all measurements using ode113. The measurement sensitivity matrix at a given time (Hi_tilde) is then mapped to the inital time using the state transition matrix and state transition tensor resulting in the measurement sensitivity matrix at the inital time (Hi). However the equation for this mapping involves xhat_0s which we do not know yet and are trying to solve for. So I keep track of all of the coefficients that would be multiplied by xhat_0s and save them as the variable Hi_sr. This variable Hi_sr is then used to make the variables Lambda_sr and N_sr which are themselves coefficients of the unknown variable xhat_0s. I continually sum all of these coefficient matricies for all measurements (a process that runs fairly quickly) so that I only have to try to use symbolic variables once at the end to solve for my unknown variable xhat_0s (a process that takes an extremely long time, and I am trying to reduce the time of).

%At each measurement the following operations are performed:


%The equation for the mapping of Hi with the troublesome xhat_0s term
%included is:
%Hi=Hi_tilde*Phi+0.5*Hi_tilde*(tensorprod(Psi,xhat_0s,2,1));

%Instead of performing the above calculation with symbolic variables, the
%coefficients of the symbolic xhat_0s are stored as follows:
Hi_sr=zeros(numMeasurementTypes,numStateVars,numStateVars+1);
for place=1:numStateVars
    Hi_sr(:,:,place)=0.5*Hi_tilde*squeeze(Psi(:,place,:));
end
Hi_sr(:,:,end)=Hi_tilde*Phi;

for index1=1:numStateVars+1
    for index2=1:numStateVars+1
        Lambda_sr(:,:,index1,index2)=Lambda_sr(:,:,index1,index2)...
            +Hi_sr(:,:,index1)'*obsData(i,3)*Hi_sr(:,:,index2);
    end
end
for index1=1:numStateVars+1
    N_sr(:,index1)=N_sr(:,index1)+Hi_sr(:,:,index1)'*obsData(i,3)*residual;
end


%Once all measurements have been processed the previously attached code is executed:
expandxmat=permute(repmat(xhat_0s_1Mat,[1,1,35,35]),[3,4,1,2]);
finalLambda=sum(sum(Lambda_sr.*expandxmat,4),3);
finalN=N_sr*xhat_0s_1;
funct=finalLambda*xhat_0s-finalN;
bigfunct=matlabFunction(funct,'Vars',{xhat_0s}); %this might be the bottleneck
residualFunct=@(x) bigfunct(x); %or maybe this is the bottleneck
initialGuess=zeros(numStateVars,1);
options=optimoptions('fsolve','Display','iter',...
    'FunctionTolerance',1e-12,'StepTolerance',1e-16,...
    'OptimalityTolerance',1e-12,'MaxIterations',50,...
    'MaxFunctionEvaluations',1000);
[xhat_0_fsolve,fval,exitFlag_fsolve,~]=fsolve(residualFunct,initialGuess,options)

My initial question was mainly formulated to address whether or not my implementation of matlabFunction and the defining of residualFunct as an anonymous function was efficient for repeated evaluations within fsolve. I do not have a ton of practice with solvers/optimizers like fsolve and lsqnonlin and so I'm not sure if there is a more effecient way to define residualFunct for use within fsolve.

Please let me know your thoughts, and again any help is much appreciated!

[Jacob]

Hello @Walter Roberson thank you for the suggestion, but unfortunately that was one of the fixes I attempted that did not seem to work. When I tried to use matlabFunction to create an optimized file the script seemed to get stuck. After many hours of waiting (at least 8, it ran over night) this file was never generated. Not sure why this was, probably because the function being generated was very complex and was taking a very long time.

[Walter Roberson]

It is true that the optimization phase can take a very long time. In theory the optimization time required rises proportional to the square of the size of the expression.

[dpb]

Given the need for ODE solver, one algorithm alternative could be <Faster Ordinary Differential Equations Solvers>.

Of course, a smaller profiling run to discover where the actual performance bottle neck(s) is(are) first would be a first step in finding out what piece(s) of the code to concentrate on...reducing a 1% area by 50% won't make a noticeable change overall; finding a hot spot that is 70% of the time spent would be something different, indeed.

Answer 1

[Jacob]

After days of testing I was able to rework my code so that it runs extremely quickly! With all the ode113 calls included my new way of doing things was able to go through 10 complete runs in approximately 3 hours. Again that is with the slow ode113 calls included, so the fsolve component is function MUCH quicker.

My problem (as usual) turned out to be the use of symbolic variables. I thought since I only used them once in the definition of residualFunct that they wouldn't be a problem. This wasn't the case however as it appears that they were slowing down fsolve drastically. I was able to resolve this issue by rewriting my code so that residualFunct was defined with exclusively function handles.

I guess I had to learn the hard way for the hundredth time to not use symbolic variables for actual calculations lol

[dpb]

I was pretty sure the symbolics would be getting in the way even without the details of the code... :J>

Glad to hear you were able to recode numerically and solve the problem...thanks for posting the resolution.