How does Matlab compute numerical jacobians?

How do Matlab functions like lsqnonlin compute numerical jacobians, and can those functions be called directly?

And I'm aware of jacobianest and have used it, but I'd like to profile against the code that those Matlab functions use

Answer 1

[Lamont Granquist]

So I got this working on R2021a, which I'll reproduce here since I've seen other people wonder how to make sfdnls or finitedifferences work and it is undocumented and I don't think anyone has posted an answer previously:

function J = jacobian_matlab(x0, tf)
  finDiffOpts.DiffMinChange = 0;
  finDiffOpts.DiffMaxChange = inf;
  finDiffOpts.TypicalX = ones(length(x0),1);
  finDiffOpts.FinDiffType = 'forward';
  finDiffOpts.GradObj = 'off';
  finDiffOpts.GradConstr = 'off';
  finDiffOpts.UseParallel = false;
  finDiffOpts.FinDiffRelStep = sqrt(eps);
  finDiffOpts.Hessian = [];

  sizes.nVar = length(x0);
  sizes.mNonlinIneq = 0;
  sizes.mNonlinEq = 0;
  sizes.xShape = size(x0);

  finDiffFlags.fwdFinDiff = true;
  finDiffFlags.scaleObjConstr = false;
  finDiffFlags.chkComplexObj = false;
  finDiffFlags.isGrad = false;
  finDiffFlags.hasLBs = false(sizes.nVar,1);
  finDiffFlags.hasUBs = false(sizes.nVar,1);
  finDiffFlags.chkFunEval = false;
  finDiffFlags.sparseFinDiff = false;

  funfcn = optimfcnchk(@(x) propagate(x, tf),'fminunc',0,false,false,false)
  f = feval(funfcn{3},x0)
  J = sfdnls(x0,f,[],[],[],funfcn,[],[],finDiffOpts,sizes,finDiffFlags)
end

The results that I got were that it is about 10x faster than jacobianest, and for my purposes has around similar loss of accuracy at the same point (although the determinant drifts in the opposite direction). I found that at 1e-9 tolerances that numerically integrating the jacobian with ode45 produces much more accurate results in about the same runtime (with higher tolerances to ode45 though it starts to take much more time).

The "propagate" function there is just a function handle to my wrapper around an ode45 problem which is the thing that I'm trying to create a jacobian of.

Probably makes more sense in retrospect to test a full problem rather than to microbenchmark the numerical jacobians, but it looks like numerically integrating the jacobian with ode45 may be a more accurate/stable and equally performant approach for me at least to try.

[ And note that if anyone else tries to get this to work in the future and this doesn't work on future versions of Matlab that you're probably on your own, I'm not providing support for figuring this out, you'll need to consult the sources of stuff like fminunc to see how the options/flags/sizes and the construction of the interior wrapper function happens in your version ]

Answer 2

[Walter Roberson]

lsqnonlin, in the case where there are no variables whose upper and lower bounds are equal, and when not configured for marquart, calls sfdnls to estimate the jacobian (the flow might be different if the user supplies jacobian functions.)

This function is not marked as "private" and is not a class member, but is clearly intended as an internal routine.