Getting the Jacobian with FindRoot output

Questions
1

I would like to extract from FindRoot the Jacobian matrix

Is there a command for this purpose that uses directly the calculation of FindRoot to get the Jacobian using Reap and Sow?

Thanks for your time.

2

I think this isn’t possible, but I may be wrong.

Example see here:

 functions = {x^2 + y^2 - 4, E^x + y - 1};
 ResourceFunction["JacobianMatrix", ResourceSystemBase -> 
 "https://www.wolframcloud.com/obj/resourcesystem/api/1.0"][functions, {x,y}]//MatrixForm

or:

vars = {x, y};
functions = {x^2 + y^2 - 4, E^x + y - 1};
jacobianUsingOuter = Outer[D, functions, vars];
MatrixForm[jacobianUsingOuter]

Regards M.I.

3

I doubt it can be done. As an alternative, you can supply your own Jacobian and monitor it. EvaluationMonitor can be used to monitor when the Jacobian is computed; however, for some reason, WRI did not supply a hook to get the value that is computed. I don’t understand why not.

Below are some variant approaches. If you have a symbolic system and the Jacobian can be computed with D[], then it’s straightforward. This probably doubles the number of times the Jacobian is evaluated though. If the objective function or functions are produced by a numerical process, then probably FindRoot[] numerically approximates the Jacobian.

(* symbolic input, jacobian *)
jac = D[{10 (-x^2 + y), 1 - x}, {{x, y}}];
Reap@FindRoot[{10 (-x^2 + y), 1 - x}, {x, -1.2`}, {y, 1.`},
  Jacobian -> {jac, EvaluationMonitor :> Sow[{x, y, jac}]}]

(* numeric objective function; by-hand central difference derivatives *)
ClearAll[obj, jac];
obj[x_?NumericQ, y_?NumericQ] := {10 (-x^2 + y), (1 - x)};
jac[x_?NumericQ, y_?NumericQ] := With[{h = 2.^-14},
  Transpose@{
   obj[x + h, y] - obj[x - h, y],
   obj[x, y + h] - obj[x, y - h]
  }/(2 h)];
Reap@FindRoot[obj[x, y], {x, -1.2`}, {y, 1.`},
  Jacobian -> {jac[x, y], EvaluationMonitor :> Sow[{x, y, jac[x, y]}]}]

(* numeric objective function; 4th-order FiniteDifferenceDerivative[] *)
ClearAll[obj, jac, dx];
dx[f_, x_, h_] = NDSolve`FiniteDifferenceDerivative[
  Derivative[1],
  x + Range[-2, 2] h,
  f /@ (x + Range[-2, 2] h)][[3]];
obj[x_?NumericQ, y_?NumericQ] := {10 (-x^2 + y), (1 - x)};
jac[x_?NumericQ, y_?NumericQ] := With[{h = 2.^-14}, 
  Transpose@{dx[obj[#, y] &, x, h], dx[obj[x, #] &, y, h]}];
Reap@FindRoot[obj[x, y], {x, -1.2`}, {y, 1.`},
  Jacobian -> {jac[x, y], EvaluationMonitor :> Sow[{x, y, jac[x, y]}]}]

(* numeric vector objective function; built-in numerical derivative *)
ClearAll[obj, jac, obj1, obj2];
obj[x_?NumericQ, y_?NumericQ] := {10 (-x^2 + y), (1 - x)};
obj1[x_?NumericQ, y_?NumericQ] := Indexed[obj[x, y], 1];
obj2[x_?NumericQ, y_?NumericQ] := Indexed[obj[x, y], 2];
jac[x_?NumericQ, y_?NumericQ] = D[{obj1[x, y], obj2[x, y]}, {{x, y}}];
Reap@FindRoot[obj[x, y], {x, -1.2`}, {y, 1.`},
  Jacobian -> {jac[x, y], EvaluationMonitor :> Sow[{x, y, jac[x, y]}]}]

If the Jacobian is expensive to compute and there are not too many evaluations, you might save time memoizing it:

jac[x_?NumericQ, y_?NumericQ] := jac[x, y] = ...

This stores and reuses computed values.

I, probably others, too, wish the numerical-function case was easier to deal with, especially the last one if obj[x, y] returns a vector that cannot be computed componentwise. I don’t know a simpler way to do it. Maybe someone else knows how to deal with it. There are special cases that could be handled, I suppose.

4

It would help to have a specific example so readers will know what sorts of equations are being considered. Offhand all I can suggest is to add tag variables and equations to recover the needed partial derivatives. E.g. add equation
df1dx3==D[f1,x3]
and add df1dx3 to the variable list