What's the best way to think about the Hessian? I've always thought about the Hessian like this:
Let $f:\mathbb R^n \to \mathbb R$ be smooth.  Let $g:\mathbb R^n \to \mathbb R^n$ such that $g(x) = \nabla f(x)$.  (I am using the convention that $\nabla f(x)$ is a column vector.)  Then the Hessian of $f$ at $x$ is, by my definition, the $n \times n$ matrix $H(x) = g'(x)$.
However, with this way of looking at the Hessian, I'm not thinking of $H(x)$ as being a quadratic form.  I'm worried that there is a "quadratic form" viewpoint of the Hessian that I am missing.  There is a rule of thumb I've heard that when a matrix (such as the Hessian) is automatically symmetric, then it's often most natural to think of it as defining a quadratic form.  I realize that you can define a quadratic form at $x_0$ by $x \mapsto \langle x, H(x_0) x \rangle$, and that this quadratic form appears in Taylor's formula.  But I still think I'm missing something, because I don't see why it is fundamentally most natural to think of the Hessian as being a quadratic form.
Is it true that there is a quadratic form viewpoint that I'm missing out on?  If so, what is it?  More generally, what do you think is the best way to think about the Hessian?
 A: The Hessian can be treated as a double derivative in the same way that the gradient or Jacobian is treated as a single derivative. The catch is that in multivariable calculus, the first derivative is no longer a scalar but instead is a linear map, while the second derivative is no longer a scalar but instead is a bilinear map. (And so on: the third derivative is a multilinear map accepting three arguments).
Let $f : \mathbb{R}^n \to \mathbb{R}$ be a smooth function, then the first derivative $Df : \mathbb{R}^n \to L(\mathbb{R}^n, \mathbb{R})$ is a smooth function from $\mathbb{R}^n$ to the set $L(\mathbb{R}^n, \mathbb{R})$ of linear maps $\mathbb{R}^n \to \mathbb{R}$. This can look complicated at first, but it really just follows straight from the definition: fixing a point $x \in \mathbb{R}^n$, the derivative must satisfy
$$ f(x + \epsilon) = f(x) + ((Df)(x))(\epsilon) + o(\epsilon),$$
where $o(-)$ is a function which goes to zero fast enough (here, faster than linear).
The derivative $Df$ itself is a smooth function so we can differentiate again to get a map $D^2 f : \mathbb{R}^n \to L(\mathbb{R}^n, L(\mathbb{R}^n, \mathbb{R}))$. The codomain $L(\mathbb{R}^n, L(\mathbb{R}^n, \mathbb{R}))$ is isomorphic to the set of bilinear maps $\operatorname{Bil}(\mathbb{R}^n \times \mathbb{R}^n, \mathbb{R})$. (To see this, just think about the fact that a function which returns another function may as well be a function of two arguments - the "linearised" version is that a linear map returning a linear map may as well just be a multilinear map of two arguments). Therefore after fixing a point $x \in \mathbb{R}^n$, the second derivative $(D^2 f)(x)$ is a bilinear map $\mathbb{R}^n \times \mathbb{R}^n \to \mathbb{R}$, which is the same kind of object that the Hessian is: we write $((D^2 f)(x))(v_1, v_2)$ rather than $v_1^T H v_2$. We then have a kind of Taylor's formula:
$$ f(x + \epsilon) = f(x) + \frac{1}{1!} ((Df)(x))(\epsilon) + \frac{1}{2!} ((D^2 f)(x))(\epsilon, \epsilon) + o(\epsilon),$$
where $o(-)$ goes to zero faster than a quadratic.
I would contend that this is the "right" way of looking at higher derivatives from a theoretical standpoint (perhaps not a computational one), since it generalises even to smooth functions $f : \mathbb{R}^n \to \mathbb{R}^m$ immediately, and does not require a choice of basis or coordinates in order to be defined. The reason that only $Df$ (the Jacobian) and $D^2 f$ (the Hessian) in a first course is that in order to write down $D^n f$ in any other case (or even $D^2 f$ in the $f: \mathbb{R}^n \to \mathbb{R}^m$ case) requires some multilinear algebra / tensors / whatever, which is normally not worth introducing in a first course on vector calculus.
As for what the operator $(D^2 f)(x)$ does intuitively: it models the concavity of a hypersurface at a point. Symmetric bilinear forms can be orthogonally diagonalised, so there is some orthonormal basis $(v_1, \ldots, v_n)$ of $\mathbb{R}^n$ so that the sign of $((D^2 f)(x))(v_i, v_i)$ tells you whether the curve traced by $f(x + tv_i)$ is concave up or down. If the sign is zero, then more work has to be done a-la the second derivative test.
A: One way to understand the Hessian better is to think of it in the language of curves.
Fix $p \in \mathbb{R}^n$. Let $q:T_p\mathbb{R}^n \to \mathbb{R}$  be the quadratic form defined by $q(v) = v^T H(p) v$. $T_p\mathbb{R}^n$ is the tangent space of $\mathbb{R}^n$ at $p$, which is just the set of all vectors starting at $p$. It's basically a copy of $\mathbb{R}^n$ but with $p$ as the origin, and is the space of all possible directions I can go starting from $p$.
$q(v)$ in some sense is the acceleration of the function $f$ in the direction of $v$. Let's see how.
Fix a vector $v\in T_p\mathbb{R}^n$. Let $\gamma:(-1,1) \to \mathbb{R}^n$ be any curve passing through $p$ with velocity $v$. So $\gamma(0) = p, \gamma'(0) = v$.
We can look at the function $f$ restricted to the curve $\gamma$ and study how fast it changes as we pass through $p$ on the curve. So we are now looking at the single variable function $g := f \circ \gamma :(-1,1) \to \mathbb{R}$.
Notice that the speed at which the function changes as we pass through $p$ on the curve is
$$g'(0) = \nabla f \cdot v$$
and the acceleration of the function as we pass through $p$ on the curve is
$$g''(0) = q(v)$$
Now observe the remarkable fact that the right side of both equations does not depend on the curve $\gamma$, but only depends on $p$ and $v$. That is why we can call $g'(0)$ the speed at which $f$ changes as we pass through $p$ in the direction of $v$ (without any mention of $\gamma$). And similarly, we can call $g''(0)$ the acceleration at which $f$ changes as we pass through $p$ in the direction of $v$ (also without any mention of $\gamma$).
A: When I happen to work in analysis, I tend to work in the part that's covered by limits of (possibly multivariable, possibly formal) power series, which you can get quite a lot of mileage from as far as intuition goes. As a basic example, the independence of the order of partial derivatives is obvious in this setting. So, I think of the Hessian in those terms first, and of the rest as generalizations of this setting. Those generalizations are appropriate for when you have less structure, but they aren't the right way to introduce the idea.
Anyway, take a multivariable polynomial $f(x_1, \ldots, x_n)$. We can write it as
$$f(x_1, \ldots, x_n) = f(0, \ldots, 0) + a_1 x_1 + \cdots + a_n x_n + \sum_{i,j=1}^n a_{i,j} x_i x_j + \cdots.$$
The constant bit is clear, and the linear bit is just $\nabla f(0) \cdot (x_1, \ldots, x_n)$. Since $x_i x_j = x_j x_i$, it seems clean to pick $a_{i,j} = a_{j,i}$. The quadratic term is
$$\sum_{i,j=1}^n a_{i,j} x_i x_j = \mathbf{x}^T A \mathbf{x}$$
--oh look, a symmetric quadratic form popped out. Pick off the coefficients with $a_{i,j} = f_{ij}(0)/2 = H_{i,j}/2$ so $A = H/2$ and
$$f(x_1, \ldots, x_n) = f(0) + \nabla f(0) \cdot \mathbf{x}^T + \frac{1}{2} \mathbf{x}^T H \mathbf{x} + \cdots$$
In this way, the Hessian obviously encodes the best quadratic approximation to the local behavior of $f$, after dealing with the constant and linear behavior. Make this precise with the formalism of your choosing--Laithy has provided one such--and bake until golden brown.
Connecting this to your definition, note that your definition of the Hessian is basically the best linear approximation to the local behavior of $\nabla f$. Since $\nabla f$ itself is the best linear approximation to the local behavior of $f$, it's at least intuitively plausible that the best linear approximation of $\nabla f$ gives the best quadratic approximation to the local behavior of $f$. So, I'd say at an intuitive level your definition is explained by my approach.
A: my tutor has notes around the topic, I hope this will be helpful
Mathematics for Intelligent Systems
