# Extract Jacobian from FindRoot?

I have to solve a system of fixed point equations and then calculate the eigenvalues of the Jacobian at the fixed point. There are about 50 equations with 50 variables and they contain a lot of numerical integrals, it would be really hard for me to give an explicit working example. I have the system of equations in the form

eqnSys={expr1,expr2,expr3,..,expr50};


the point is, it is a list of expressions which has to equal zero in the end and it is not defined as a function, but I have the variables saved in

eqnVars={x1,x2,..,x50};


I have very good initial guess for the root finding:

eqnGuess={{x1,1},{x2,2},..,{x50,50}};


The basic Newton method of FindRoot breaks down immediately, complaining about singular Jacobian. So this

FindRoot[eqnSys,eqnGuess]


does not work. I used the secant method for a long time to find the roots:

eqnGuessSec={{x1,0.9,1.1},{x2,1.8,2.2},..,{x50,45,55}};
FindRoot[eqnSys,eqnGuessSec]


Recently I came across the AffineCovariantNewton method for FindRoot which works like a charm and outperforms the secant method in time by a factor of 4:

FindRoot[eqnSys,eqnGuess,Method -> {"AffineCovariantNewton"}]


Judging from the evaluation monitor, it has several Jacobian evaluations:

FindRoot[eqnSys,eqnGuess,Method -> {"AffineCovariantNewton"},Jacobian -> {Automatic, EvaluationMonitor :> Print["J evaluated here"]}]


My question is: It would be really lucrative for me to be able to save the Jacobian directly from FindRoot. Is it possible to extract the Jacobian matrix constructed by FindRoot? I'm thinking about something like

Reap@FindRoot[eqnSys,eqnGuess,Method -> {"AffineCovariantNewton"},Jacobian -> {Automatic,Sow[jacobian]}]


I'm only interested in the purely numerical matrix and not a symbolic one. Bonus question: What is the most efficient way to turn the system of equations to a function? So something like

FeqnSys[x1_,x2_,...,x50_]:=eqnSys


Edit: I have implemented a very simple version of the problem. Increasing UTrunc increaseas the number of equations (but then need additional initial conditions). I basically need the object with the name ineedthisguy. I hoped that it can be obtained without this analytical differention, because for a real problem I can only generate the full matrix in chunks due to memory limitation.

d = 3;
WorPrec = 16;
\[Alpha] = 1;

UTrunc = 6;

Z[r_] := 0
W[r_] := 0

U[r_] := Sum[
ToExpression["u" <> ToString[n]]/n! (r - \[Kappa])^n, {n, 2,
UTrunc}];
\[Omega][r_] := U'[r] + 2 r U''[r]

MasterKernel1[d_, n1_, \[Omega]_?NumericQ, w_?NumericQ] :=
MasterKernel1[d, n1, \[Omega],
w] = -2 \[Alpha] NIntegrate[
E^-y y^(-1 + d/
2) (1 + y) (y + w y^2 + E^-y \[Alpha] + \[Omega])^-n1, {y,
0, \[Infinity]},
Method -> {Automatic, "SymbolicProcessing" -> False},
WorkingPrecision -> WorPrec]

Derivative[1][MasterL[n_, d_]][\[Rho]_] :=
Derivative[1][
MasterL[n,
d]][\[Rho]] = -n (MasterL[n + 1, d][\[Rho]] \[Omega]'[\[Rho]] +
MasterL[pa][n + 1, d + 2][\[Rho]] Z'[\[Rho]] +
MasterL[pa][n + 1, d + 4][\[Rho]] W'[\[Rho]])

MasterL[n_, d_][\[Kappa]] :=
MasterL[n, d][\[Kappa]] =
MasterKernel1[d, n, 2 \[Kappa] u2, W[\[Kappa]]]

BetaU[r_] := -d U[r] + (d - 2) r U'[r] -
1/(4 \[Pi]^2) MasterL[1, d][r]

dExpr[f_, betafunc_, n_] := D[k D[f[r], k] == betafunc[r], {r, n}]

GenBeta[f_, betafunc_, min_, max_] := Block[{expr, result, tmpres},
expr = dExpr[f, betafunc, min];
result = {(expr /. r -> \[Kappa])};
Do[
expr = D[expr, r];
tmpres = Block[{r = \[Kappa]}, expr];
result = Join[result, {tmpres}];
, {i, min + 1, max}
];
Return[result];
];

listU = GenBeta[U, BetaU, 1, UTrunc];

FPEqn = ((Flatten@(List @@@ Flatten[listU]))[[2 ;; ;; 2]]);

varTrf = {g_[n_] :> ToExpression[ ToString[g] <> ToString[n]]};
varList = Flatten[{\[Kappa], Table[u[i], {i, 2, UTrunc}]}];

iniGuess =
Rationalize[
List @@@ {\[Kappa] -> 0.0417487541261041756617205337339609668612.,
u2 -> 6.1458403749048580482270685737667568587812.,
u3 -> 60.0491811653211896544374917466544653009612.,
u4 -> 390.901060703305764622212.,
u5 -> -3513.611214090298842396512.,
u6 -> -93676.707982735664990099912.}, 0];
(*real solution:
{\[Kappa]\[Rule]0.0726928522670547,u2\[Rule]4.570711765672155,u3\
\[Rule]28.871831592476088,u4\[Rule]134.9966784017132,u5\[Rule]-371.\
15673934569224,u6\[Rule]-14195.11815231752}
*)

fOPT = ExperimentalOptimizeExpression[FPEqn,
"OptimizationLevel" -> 2]; (*im not sure if this helps*)

fpLocator[initial__] :=
FindRoot[fOPT // First, List @@@ initial,
Method -> {"AffineCovariantNewton"}, WorkingPrecision -> WorPrec,
StepMonitor :> {Print[initial[[All, 1]]]} ] ;

sol = fpLocator[iniGuess]

\!$$\*SuperscriptBox[\(MasterKernel1$$,
TagBox[
RowBox[{"(",
RowBox[{"0", ",", "0", ",", "1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[d_, n1_, \[Omega]_?NumericQ,
w0_?NumericQ] :=
\!$$\*SuperscriptBox[\(MasterKernel1$$,
TagBox[
RowBox[{"(",
RowBox[{"0", ",", "0", ",", "1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[d, n1, \[Omega],
w0] = -n1 MasterKernel1[d, n1 + 1, \[Omega], w0]

\!$$\*SuperscriptBox[\(MasterKernel1$$,
TagBox[
RowBox[{"(",
RowBox[{"0", ",", "0", ",", "0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[d_, n1_, \[Omega]_?NumericQ,
w0_?NumericQ] :=
\!$$\*SuperscriptBox[\(MasterKernel1$$,
TagBox[
RowBox[{"(",
RowBox[{"0", ",", "0", ",", "0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[d, n1, \[Omega],
w0] = -n1 MasterKernel1[d + 4, n1 + 1, \[Omega], w0]

ineedthisguy = Eigenvalues[D[FPEqn, {varList /. varTrf}] /. sol]

• D[eqnSys /. Equal -> Subtract, {eqnVars}]? If not, it must depend on eqnSys, which you haven't shared. The Jacobian might be evaluated numerically, and then there is no single matrix; the matrix is updated at each step (unless it accepts the "UpdateJacobian", in which case you can specify how often it is updated). Sep 3 '20 at 13:08
• I can generate symbolically the Jacobian, but it takes a lot of time both generating it and then substituting the values into it. I only need the matrix with only numerical entries/numerically evaluated matrix. Sep 3 '20 at 13:19
• Do you want each matrix at each step, or just one matrix at some particular point in the domain? Sep 3 '20 at 13:23
• I only want the last one, which corresponds to the solution itself or close to it (last few digits does not really matter of the solution). Sep 3 '20 at 13:26
• You could try complex-step differentiation. I worked up an example before you posted yours, and it relied on the system being made up of Listable expressions; you could adapt it. Here's key code: jac = Im[ sys /. Thread[ vars -> x0 + $MachineEpsilon*I*IdentityMatrix[Length@vars] ] ]/$MachineEpsilon, but this was faster for me: jac = Hold@Evaluate@vars /. Hold[v_] :> Block[v, v = x0 + $MachineEpsilon*I*IdentityMatrix[Length@vars]; Im@sys/$MachineEpsilon ] Sep 3 '20 at 17:38

Here is a way to extract the Jacobian. The idea is to re-write the linear solver to Sow the Jacobian and then to Reap that:

Setup a simple problem:

f[X_] := Block[{x, y}, {x, y} = X; {Exp[x - 2] - y, y^2 - x}]
vars = {x, y};
start = {1, 1};
callFindRoot[f_, vars_, start_, opts___] :=
vars /. FindRoot[f[vars], Evaluate[Transpose[{vars, start}]], opts]


Re-write the LinearSolver to Sow the Jacobian:

MyLinearSolver = (Sow[#]; LinearSolve[##]) &;


Call FindRoot

Reap[callFindRoot[f, vars, start, Method -> {"AffineCovariantNewton"
, "LinearSolver" -> {MyLinearSolver}
}]]

(* {{0.019026,
0.137935}, {{{{0.367879, -1.}, {-1., 2.}}, {{0.295112, -1.}, {-1.,
1.7796}}, {{0.191064, -1.}, {-1.,
1.28831}}, {{0.096665, -1.}, {-1., 0.121763}}}}} *)


I'll add this example to the documentation.

To turn the equations into a function you could use something along the lines of:

cf = With[{vars = vars, fun = f[vars]},
Compile[{{X, _Real, 1}},
Block[vars,
vars = X;
fun
]
]
]

• Awesome, thank you! I'll accept your answer in the evening (morning here atm). Maybe if you have time, could you include an answer to my bonus question? Sep 4 '20 at 6:45
• Is it possible to something like this with FindMinimum to sow the Hessian matrices that were used? Sep 4 '20 at 8:22
• @SjoerdSmit, not that I know. But you can create a NumericalFunction and extract the Hessian from there. I am sure FindMinimum` uses NumericalFunctions - though I have never tried not looked. Sep 4 '20 at 8:47
• @Sjoerd, have you already seen this? Sep 12 '20 at 5:59