1 iteration in fsolve: new x and no more

I want to use fsolve to only use one iteration: get the Jacobian, get x1 and done. However when test the code below, which tries to find the zero of x^2, i see the following:

-iteration 0, 2 function evaluations: it seems f(x0) and 1 for finite difference

-iteration 1, 2 function evaluations: it seems f(x1) and 1 for finite difference

I do not want the last two evaluations. Is this possible?

I mean, this doubles the computation time, but it adds nothing for me. Especially since the Jacobian that is reported is the jacobian at x0 (not at x1). If it would at least give me the jacobian at x1 as output i could calculate the next iteration myself.

The solution would to program a similar routine myself, but I rather use fsolve (i also want to use JacobPattern for example).

clear all;
clc;

fun_res = @(xx)xx.^2;

x0 = 1;

ex_fl = -10;
iter = 0;
while ex_fl < 1 && iter < 100
    iter = iter+1;
    %[xx_it(iter,1),res(iter,1),ex_fl,~,jac(iter,1)] = fsolve(fun_res,x0,optimoptions(@fsolve,'Display','iter','MaxIterations',1));
    [xx_it(iter,1),~,ex_fl] = fsolve(fun_res,x0,optimoptions(@fsolve,'Display','iter','MaxIterations',1));

    x0 = xx_it(iter);   
    
end

[Matt J]

If it would at least give me the jacobian at x1

I don't know why you think it isn't. The documentation states that the Jacobian should that at the "solution".

[Sargondjani]

Actually it's true. I misread the results, I guess. Anyway, i found out that the first step gets really close to the true solution (in my actual program, so this second step doesnt break me much anyway).

Answer 1

[Matt J]

One option is to use fsolve just to get the Jacobian and then generate x1 yourself. This gives you access to JacobPattern and all the other Jacobian calculation features that fsolve offers.

To do so, replace your objective function f(x) with f(x)-f(x0), where x0 is thte initial point. Note that this does not change the Jacobian. Because x0 is a root of the modified objective, fsolve will return x0 and its Jacobian as output in zero iterations. Example:

fun=@(x) x.^2/2;

options=optimoptions(@fsolve,'MaxIter',1);

x0=rand(5,1);
f0=fun(x0);
[x,fval,ef,stats,JacobianNumerical]=fsolve(   @(x)fun(x) - f0 ,x0,options);

Equation solved at initial point.

fsolve completed because the vector of function values at the initial point
is near zero as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

JacobianNumerical,

JacobianNumerical = 5×5

0.3167 0 0 0 0 0 0.5539 0 0 0 0 0 0.7692 0 0 0 0 0 0.0906 0 0 0 0 0 0.5526

JacobianAnalytical=diag(x0)

JacobianAnalytical = 5×5

0.3167 0 0 0 0 0 0.5539 0 0 0 0 0 0.7692 0 0 0 0 0 0.0906 0 0 0 0 0 0.5526

Also, it will do this with no extra Jacobian calculations, as seen from,

stats

stats = struct with fields:
       iterations: 0
        funcCount: 6
        algorithm: 'trust-region-dogleg'
    firstorderopt: 0
          message: '↵Equation solved at initial point.↵↵fsolve completed because the vector of function values at the initial point↵is near zero as measured by the value of the function tolerance, and↵the problem appears regular as measured by the gradient.↵↵<stopping criteria details>↵↵Equation solved. The final point is the initial point.↵The sum of squared function values, r = 0.000000e+00, is less than sqrt(options.FunctionTolerance) = 1.000000e-03.↵The relative norm of the gradient of r, 0.000000e+00, is less than options.OptimalityTolerance = 1.000000e-06.↵↵'

[Sargondjani]

This looks good!! I haven't check it, but since you wrote it, I assume it works :-)

[Matt J]

Note, you should pre-calculate f0=f(x0) to avoid that calculation being repeated in the finite differencing.

f0=fun(x0);

modfun=@(x) fun(x)-f0;

[x,fval,ef,stats,JacobianNumerical]=fsolve( modfun, x0,options);

You can also avoid one extra function evaluation by getting a little fancier:

f0=fun(x0);

[x,fval,ef,stats,JacobianNumerical]=fsolve( @(x) modfun(x,fun,f0,x0), x0,options);

function fval=modfun(x,fun,f0,x0)

  if isequal(x,x0)
      fval=zeros(size(f0));
  else
      fval=fun(x)-f0; 
  end

end

[Sargondjani]

Dear Matt, I have now programmed the Newton–Raphson step myself (using a very simple routine to calculate the Jacobian, since i only need the diagional entries).

This turned out to be 100 times faster than doing it with fsolve. (0.1 second versus 9 seconds, admittedly without using your method, so I could potentially reduce it to 5 seconds, but it seems it will never beat the simple Newton Raphson program).

I should add I didnt change much esle as I pass the function handle to both the numerical jacobian estimator, and fsolve.

Did you ever have similar huge effeciency gains by programming a Newton-Raphson method versus fsolve? I'm curious what causes this enormous difference. I mean, it seems its just overhead from fsolve???

[Matt J]

I'm sure there is a lot of overhead from fsolve. Using it one iteration at a time, just to compute the Jacobian, would re-incur start-up cost every time you run it.

[Sargondjani]

Ok, good to know! Thanks :-)

[Sargondjani]

@Matt J for completeness i justed wanted to say that the speed difference was NOT caused by fsolve. But because the default algorithm of fsolve simply ignores JacobPattern. So apperently it was computing the full Jabobian, even though i specified the pattern. I don't know why it didnt warn me that i specifiy an option which is not avaiable for that algorithm.

Answer 2

[Matt J]

The solution would to program a similar routine myself, but I rather use fsolve (i also want to use JacobPattern for example).

You could also use an existing routine from the File Exchange, like this one,

https://www.mathworks.com/matlabcentral/fileexchange/13490-adaptive-robust-numerical-differentiation

It wouldn't be so difficult to adapt the routine to work with a JacobPattern.

[Sargondjani]

Thanks Matt for helping. And I have a routine to numerically calculate the Jacobian, and maybe it's not that much extra work to implement a JaobianPattern. I had just hoped it would be possible with fsolve.