Determining convergence rate of fsolve - Matlab

Question

Suppose I'm solving a system of nonlinear equations. A simple example would be:

function example
    x0 = [15; -2];
    options = optimoptions('fsolve','Display','iter','TolFun',eps,'TolX',eps);
    [x,fval,exitflag,output] = fsolve(@P1a,x0,options);
end

function f1 = P1a(x)
    f1 = [x(1)+x(2)*(x(2)*(5-x(2))-2)- 13; x(1)+x(2)*(x(2)*(1+x(2))-14)-29];
end

How do I determine the rate of convergence? 'Display','iter' shows me the norm at each step, but I can't find a way to extract these values. (For this particular example, I believe fsolve does not converge to the right solution, but rather to a local minimum. That is not the issue, however. I just want to find a way to estimate the convergence rate.)

Answer 1

You can get plenty out of fsolve. However, you'll need to do some work. Read up on the 'OutputFcn' option and writing output functions for Matlab's Optimization methods. This is vary similar to the option of the same name used by Matlab's ODE solvers. Here's an example that replicates the 'Display','iter' option-value for fsolve (for the default 'trust-region-dogleg' algorithm specifically):

function stop = outfun(x,optimValues,state)
% See private/trustnleqn
stop = false;

switch state
    case 'init'
        header = sprintf(['\n                                         Norm of      First-order   Trust-region\n',...
                          ' Iteration  Func-count     f(x)          step         optimality    radius']);
        disp(header);
    case 'iter'
        iter = optimValues.iteration;               % Iteration
        numFevals = optimValues.funccount;          % Func-count
        F = optimValues.fval;                       % f(x)
        normd = optimValues.stepsize;               % Norm of step
        normgradinf = optimValues.firstorderopt;    % First-order optimality
        Delta = optimValues.trustregionradius;      % Trust-region radius

        if iter > 0
            formatstr = ' %5.0f      %5.0f   %13.6g  %13.6g   %12.3g    %12.3g';
            iterOutput = sprintf(formatstr,iter,numFevals,F'*F,normd,normgradinf,Delta);
        else
            formatstr0 = ' %5.0f      %5.0f   %13.6g                  %12.3g    %12.3g';
            iterOutput = sprintf(formatstr0,iter,numFevals,F'*F,normgradinf,Delta);
        end
        disp(iterOutput);
    case 'done'

    otherwise

end

You can then call this via:

function example
P1a=@(x)[x(1)+x(2)*(x(2)*(5-x(2))-2)- 13; x(1)+x(2)*(x(2)*(1+x(2))-14)-29];
x0 = [15; -2];
opts = optimoptions('fsolve','Display','off','OutputFcn',@outfun,'TolFun',eps,'TolX',eps);
[x,fval,exitflag,output] = fsolve(P1a,x0,opts);

This still just prints to the Command Window. From here it's a matter of creating an output function that can write data to an array, file, or other data structure. Here's how you might do that with a global variable (in general, not a good idea):

function stop = outfun2(x,optimValues,state)
stop = false;
global out;   % Global variable, define in main function too

switch state
    case 'init'
        out = [];
    case 'iter'
        iter = optimValues.iteration;               % Iteration
        numFevals = optimValues.funccount;          % Func-count
        F = optimValues.fval;                       % f(x)
        normd = optimValues.stepsize;               % Norm of step
        normgradinf = optimValues.firstorderopt;    % First-order optimality
        Delta = optimValues.trustregionradius;      % Trust-region radius

        out = [out;iter numFevals F'*F normd normgradinf Delta];
    case 'done'

    otherwise

end

Then just declare global out; in your main function before calling fsolve. You could also accomplish this by making your output function a nested function, in which case the out array would be shared with the outer main function.

The second output function example performs dynamic memory allocation instead of reallocating the entire out array. There's no way around this because both we and the algorithm don't know how many iterations it will take to converge. However, for a few hundred iterations, dynamic memory allocation will be plenty fast.

I'll leave "determining the rate of convergence" to you now that you have the tools in hand...