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?

  • 3
    $\begingroup$ 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. $\endgroup$ – Dr. belisarius Feb 14 '13 at 18:55
  • 1
    $\begingroup$ 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. $\endgroup$ – Oleksandr R. Feb 14 '13 at 18:58
  • 1
    $\begingroup$ @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. $\endgroup$ – Pranas Feb 14 '13 at 19:15
  • 1
    $\begingroup$ 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) $\endgroup$ – Mike Honeychurch Feb 14 '13 at 21:21
  • $\begingroup$ 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 $\endgroup$ – Sjoerd C. de Vries Feb 14 '13 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}}]


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


$\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

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


  • $\begingroup$ 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. $\endgroup$ – TheGrapeBeyond Feb 17 '14 at 15:38
  • 1
    $\begingroup$ @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}}] $\endgroup$ – Jens Feb 20 '14 at 3:35
  • $\begingroup$ 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 $\endgroup$ – step-by-step Dec 13 '18 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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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