There's something that's always been deeply confusing to me about comparing the Jacobian and the Hessian because their nature is very different.

The Hessian shouldn't have been called a matrix.

The Jacobian describes all the first order derivatives of a vector valued function (of multiple inputs), while the Hessian is all the second order derivatives of a scalar valued output function (of multiple inputs). Why doesn't the number of dimensions of the array increase by one as the derivation order increases? It does! The object that fully describes second order derivation of a vector valued function of multiple inputs is actually a 3 dimensionnal tensor. One dimension for the original vector valued output, and one for each derivation order. Mathematicians are afraid of tensors of more than 2 dimensions for some reason and want everything to be a matrix.

In other words, given a function R^n -> R^m:

Order 0: Output value: 1d array of shape (m) (a vector)

Order 1: First order derivative: 2d array of shape (m, n) (Jacobian matrix)

Order 2: Second order derivative: 3d array of shape (m, n, n) (array of Hessian matrices)

It all makes sense!

Talking about "Jacobian and Hessian" matrices as if they are both naturally matrices is highly misleading.

At least in my undergrad multivariate real analysis class, I remember the professor arranging things to strongly suggest that the Hessian should be thought of as ∇⊗∇, and that this was the second term in a higher dimensional Taylor series, so that the third derivative term would be ∇⊗∇⊗∇ etc. Things like tensor products or even quotient spaces weren't assumed knowledge, so it wasn't explicitly covered, but I remember feeling the connection was obvious enough at the time. Then an introductory differential geometry class got into (n,m) tensors. So I'm quite sure mathematicians are fine dealing with tensors. My experience was undergrad engineering math tries to avoid even covectors though, so that will stay well clear of a coherent picture of multi-variable calculus. e.g. my engineering professors would talk of dirac δ as an infinite spike/spooky doesn't-really-exist thing that makes integrals work or whatever. My analysis professor just said δ(f) = f(0) is a linear functional.

You're confusing too many things.

The Hessian is defined as the second order partial derivative of a scalar function. Therefore it will always give you a matrix.

What you're doing with the shape (m,n,n) isn't actually guaranteed at all since the output shape of an arbitrary function can be any tensor and you can apply the Hessian to each scalar value in the tensor to get another arbitrary tensor that has two dimensions more.

It's the Jacobian that is weird, since it is just a vector of gradients and therefore its partial derivative must also be a vector of Hessians.

This doesn't really help with programming, but in physics it's traditional to use up- and down-stairs indices, which makes the distinction you want very clear.

If input x has components xⁿ, and output f(x) components fᵐ, then the Jacobian is ∂ₙfᵐ which has one index upstairs and one downstairs. The derivative has a downstairs index... because x is in the denominator of d/dx, roughly? If x had units seconds, then d/dx has units per second.

Whereas if g(x) is a number, the gradient is ∂ₙg, and the Hessian is ∂ₙ∂ₙ₂g with two downstairs indices. You might call this a (0,2) tensor, while the Jaconian is (1,1). Most of the matrices in ordinary linear algebra are (1,1) tensors.

I’ve been introduced to the Hessian in the context of finding the extrema of functions with multiple variables, where it does not make sense to consider arbitrary output dimensions (what is the minimum of a vector function?). In this context, it is also important to find the definiteness of the underlying quadratic form, which is easier if you treat it as a matrix se you can apply Sylvester’s rule.

I agree it is confusing, because starting with notation will confuse you. I personally don't like the partial derivative-first definition of those concepts, as it all sounds a bit arbitrary.

What made sense to me is to start from the definition of derivative (the best linear approximation in some sense), and then everything else is about how to represent this. vectors, matrices, etc. are all vectors in the appropriate vector space, the derivative is always the same form in a functional form, etc.

E.g. you want the derivative of f(M) ? Just write f(M+h) - f(M), and then look for the terms in h / h^2 / etc. Apply chain rules / etc. for more complicated cases. This is IMO a much better way to learn about this.

As for notation, you use vec/kronecker product for complicated cases: https://janmagnus.nl/papers/JRM093.pdf

Well for me the Hessian is the second order derivative in the special case where the co-domain is of dim 1. It's just very easy to work with...

Mathematicians are afraid of higher order tensors because they are unruly monsters.

There's a whole workshop of useful matrix tools. Decompositions, spectral theory, etc. These tools really break down when you generalize them to k-tensors. Even basic concepts like rank become sticky. (Iirc, the set of 3-tensors of tensor rank ≤k is not even topologically closed in general. Terrifying.) If you hand me some random 5-tensor, it's quite difficult to begin to understand it without somehow turning it into a matrix first by flattening or slicing or whatever.

Don't get me wrong. People work with these things. They do their best. But in general, mathematicians are afraid of higher order tensors. You should be too.