Quick Hessian matrix and gradient calculation?

Question

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?

Answer 1

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]

Gradient:

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

Answer 2

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

Answer 3

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}]