According to the documentation, when FindMinimum is told to use the method "QuasiNewton" on a unconstrained problem, it uses the BFGS approximation to the Hessian matrix of the function being optimized. How can I recover the Hessian matrix computed by Mathematica using StepMonitor, EvaluationMonitor or any other devious trick?


2 Answers 2


Unfortunately, this requires a lot more devious trickery than I would have preferred. As noted in the documentation at tutorial/UnconstrainedOptimizationQuasiNewtonMethods, the Hessian is not formed directly in the BFGS method, so we have to recover it from the Cholesky factors. However, all of this is done inside the kernel where we cannot access it using ordinary methods, and the only way I can see to obtain these factors is to hook calls to the BLAS function *TRSV (called as LinearAlgebra`BLAS`TRSV in Mathematica) and extract its arguments directly. Needless to say, this is hardly likely to be particularly robust, but it hopefully it will work well enough for your needs.

Let us define:

approximateHessianList[f_, vars_, start_] :=
  With[{fmarg2 = Thread[{vars, start}]},
   Module[{steps, lowerTriangles, upperTriangles},
     trsv : LinearAlgebra`BLAS`TRSV[uplo_, trans_, diag_, a_, x_] /; (
        "N", Sow[LowerTriangularize[a], "LowerTriangle"],
        "T", Sow[UpperTriangularize[a], "UpperTriangle"]
       ]; True
      ) := trsv;
     {steps, lowerTriangles, upperTriangles} =
        f, fmarg2,
        Method -> {"QuasiNewton", "StepMemory" -> Infinity},
        StepMonitor :> Sow[vars, "Step"]
       ], {"Step", "LowerTriangle", "UpperTriangle"}
      ][[2, {1, 2, 3}, 1]];
     Transpose[{steps, MapThread[Dot, {lowerTriangles, upperTriangles}]}]
SetAttributes[approximateHessianList, HoldAll];

We may now write:

approximateHessianList[Cos[x^2 - 3 y] + Sin[x^2 + y^2], {x, y}, {1, 1}]

which gives a list of steps in the BFGS optimization along with the approximate Hessians evaluated at those points. For the sake of brevity (and because it is obviously the most accurate), let us take only the last of these:

{{1.37638, 1.67868}, {{15.1553, 0.986982}, {0.986982, 20.2017}}}

which we may compare to the exact Hessian evaluated at this point:

D[Cos[x^2 - 3 y] + Sin[x^2 + y^2], {{x, y}, 2}] /. {
  x -> 1.376384972443001`, y -> 1.6786760817546214`
{{15.1555, 0.983708}, {0.983708, 20.2718}}

Clearly, it is not a bad approximation.

Now, some caveats. Since the logic of the code inside the kernel is completely opaque, I am not sure whether this business with the upper and lower triangles is really necessary, since the argument a of LinearAlgebra`BLAS`TRSV appears to be the same for successive calls, firstly with trans == "N" and then with trans == "T". However, this represents at most an inefficiency. A more serious problem is that I am not sure whether the first Hessian obtained using this method corresponds to the initial point, $(x,y)=(1,1)$ , or the point at the end of the first step, which here is $(x,y)=(0.811216, 1.68144)$ . While this does not really matter if the approximation has converged, it does influence the interpretation of the earlier approximations substantially, so hopefully someone with access to the kernel code will be able to clarify this aspect of FindMinimum`QuasiNewton's behaviour.

  • $\begingroup$ I'd plus you 10,000 points if I could, Oleksandr. Can you elaborate briefly on the line which implements the hook? I haven't mastered all the ":" and "/;" tricks, and can't quite see what that line does. Brilliant work. $\endgroup$
    – Ian
    Commented Mar 19, 2012 at 22:55
  • 1
    $\begingroup$ This is closely related to the Villegas-Gayley trick, but a little more concise when you don't need another Block (for example, if the definition is already inside an Internal`InheritedBlock). Here we pattern match the function call (with a named pattern), and take advantage of a Condition containing (code; True) so that code executes (once only) each time the hooked function is called. The call proper is passed through the hook by means of the pattern name (here trsv). $\endgroup$ Commented Mar 19, 2012 at 23:21
  • 7
    $\begingroup$ This is insane! $\endgroup$
    – rm -rf
    Commented Mar 20, 2012 at 2:08
  • $\begingroup$ I have to agree with @R.M, the insanity for me was trsv: ... := trsv to attach the hook. $\endgroup$
    – rcollyer
    Commented Mar 20, 2012 at 14:59

Using Experimental`NumericalFunction framework directly (FindMinimum uses it under the hood) it is straightforward to get the numerical approximation of the Hessian:

f = Experimental`CreateNumericalFunction[{x, y}, 
   Cos[x^2 - 3 y] + Sin[x^2 + y^2], {1}, Hessian -> FiniteDifference];

f["Hessian"[{1.376384972443001`, 1.6786760817546214`}]]

{{{15.1555, 0.983708}, {0.983708, 20.2718}}}

  • $\begingroup$ As far as Hessian is concerned, what are the other options apart from FiniteDifference? Thanks! $\endgroup$ Commented Nov 29, 2014 at 9:41
  • 1
    $\begingroup$ As explained in the linked answer, Hessian takes three possible values: Automatic, Symbolic and FiniteDifference with Automatic being equivalent to Symbolic AFAIK. $\endgroup$ Commented Nov 29, 2014 at 9:54
  • $\begingroup$ @Alexey Popkov How is the scale chosen? Are finite differences adaptive? Is the Richardson extrapolation methodology used? $\endgroup$
    – Valerio
    Commented May 18, 2018 at 19:42
  • $\begingroup$ @Valerio Since Experimental`NumericalFunction is undocumented, we can only guess and experiment. You can create a separate question on it. $\endgroup$ Commented May 19, 2018 at 2:32
  • $\begingroup$ @Alexey Popkov I created a new question as you suggested. $\endgroup$
    – Valerio
    Commented May 19, 2018 at 15:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.