# Quick Hessian matrix and gradient calculation?

I am absolutely new to Mathematica and I actually want to try implementing a little optimization method .

Long story short assuming I have a predefined two-variable function f(x,y) I want to calculate a Hessian matrix and a gradient symbolically. Then I want to be able to quickly plug specific x,y values into them.

How is it done in Mathematica?

• Welcome to the site! Please show some effort. We are always glad to help newcomers here, as far as we aren't taken for a coding service. Commented Feb 14, 2013 at 18:55
• First of all, Mathematica already incorporates many optimization methods--see the documentation pages for FindMinimum (local/gradient optimizer, including Newton's method) and NMinimize (global/derivative-free optimizer). If you really want to code your own implementation (which should not be difficult), see the documentation for the function D. This will give you your Jacobian and Hessian directly. Commented Feb 14, 2013 at 18:58
• @belisarius Thanks for your response. Unfortunately at this moment I'm failing at a very basic level which makes it hard to show the effort. Anyways I'm not asking for a finished implementation. Commented Feb 14, 2013 at 19:15
• Jacobian and Hessian are also given here: tutorial/Differentiation (note that many of these sorts of things are often hard to find unless you already know they are there) Commented Feb 14, 2013 at 21:21
• To be sure: The page @MikeHoneychurch mentions is the one you get if you enter those keywords in the Mathematica Documentation Center (in the Help menu). On the internet that would be: reference.wolfram.com/mathematica/tutorial/Differentiation.html Commented Feb 14, 2013 at 22:07

Your solution was almost correct, except that it should make f an argument of the hessian function and could implement the derivatives in a more compact way. As pointed out by Mike Honeychurch in the above comments, the first place to start would be to look at the documentation on differentiation.

Here is how the derivative operator D can be used to define gradients and hessians:

Clear[f, hessian]


D[f[x, y], {{x, y}}]


$\left\{f^{(1,0)}(x,y),f^{(0,1)}(x,y)\right\}$

Hessian (alternative formulation D[f[x,y],{{x,y},2}]):

D[f[x, y], {{x, y}}, {{x, y}}]


$\left( \begin{array}{cc} f^{(2,0)}(x,y) & f^{(1,1)}(x,y) \\ f^{(1,1)}(x,y) & f^{(0,2)}(x,y) \\ \end{array} \right)$

Now to define the latter as an operator:

hessian[x_, y_] = Function[{f},
D[f, {{x, y}}, {{x, y}}]
];

f[x_, y_] := (x^2 - y)/(x^2 + y^2 + 1);

hessian[x, y][f[x, y]]//FullSimplify


$\left( \begin{array}{cc} -\frac{2 \left(3 x^2-y^2-1\right) \left(y^2+y+1\right)}{\left(x^2+ y^2+1\right)^3} & \frac{2 x \left((2 y+1) x^2-y (y (2 y+3)+2)+1\right)}{\left(x^2+y^2+ 1\right)^3} \\ \frac{2 x \left((2 y+1) x^2-y (y (2 y+3)+2)+1\right)}{\left(x^2+y^2+ 1\right)^3} & -\frac{2 \left(x^4+(1-3 y (y+1)) x^2+y^3-3 y\right)}{\left(x^2+y^2+1\right) ^3} \\ \end{array} \right)$

To plug in specific values for x and y, one approach would be to follow the last result by

%/.{x->2,y->3}


$\left( \begin{array}{cc} -\frac{13}{686} & -\frac{29}{343} \\ -\frac{29}{343} & \frac{53}{686} \\ \end{array} \right)$

Here, the % recalls the result of the previous output, and the /. stands for ReplaceAll.

Edit: generalization to n-th derivative

To generalize the Hessian above, you can get the tensor of n-th derivatives as follows:

D[f[x, y, z], {{x, y, z}, n}]


where n is the order of the derivative. For example, with n=3 you get

Clear[f];
D[f[x, y, z], {{x, y, z}, 3}] // TableForm


• Hello, I share the OP's noobness. I have a question. Can I actually do this on the wolfram alpha website itself? It only shows me one line where I can input stuff... or do I have to buy the full blown Mathematica software? Thanks. Commented Feb 17, 2014 at 15:38
• @TheGrapeBeyond It also works in Alpha if you input the one-line command including the function definition: D[x^2 - y)/(x^2 + y^2 + 1), {{x, y}}, {{x, y}}]
– Jens
Commented Feb 20, 2014 at 3:35
• Quite a long time now, perhaps @Jens code can be reduce a bit by: D[x^2 - y)/(x^2 + y^2 + 1), {{x, y}}, 2}] //FullSimplify Commented Dec 13, 2018 at 15:49

I know this question is a bit old but Wolfram has an official solution on their MathWorld website which I've tested and it works well.

Just use the following function:

HessianH[f_, x_List?VectorQ] := D[f, {x, 2}]


It's used like one would expect:

(* Your function *)
f[x_, y_, z_] := 100*(y - x^2) + (1 - x)^2 + 100 (z - y^2) + (1 - y)^2

(* Just pass whatev vars you want to compute the Hessian of *)
HessianH[f[x, y, z], {x, y, z}] // MatrixForm


Like Samuel, I know this is a bit old, but the documentation for the function Grad gives an even simpler solution, by simple repeated application (and which allows use of curvilinear coordinates etc.):

Grad[Grad[f[x,y,z],{x,y,z}],{x,y,z}]

• In what way is this superior to D[f[x,y,z],{{x,y,z},2}]?
– Alan
Commented Dec 11, 2023 at 14:15