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.