Calculate the Hessian of a Vector Function I'm working with optimisation. I am trying to obtain the hessian of a vector function:
$$
\mathbf{F(X) = 0} \quad \text{or} \quad
  \begin{cases}
    f_1(x_1,x_2,\dotsc,x_n) = 0,\\
    f_2(x_1,x_2,\dotsc,x_n) = 0,\\
    \vdots\\
    f_n(x_1,x_2,\dotsc,x_n) = 0,\\
  \end{cases}
$$
I know that the Jacobian for a vector function is calculated as:
$$
\mathbf{J}= \begin{bmatrix}
\frac{\partial f_1}{\partial x_1} & \dots & \frac{\partial f_1}{\partial x_1} \\  \vdots & \ddots & \vdots \\ \frac{\partial f_n}{\partial x_1} & \dots &\frac{\partial f_n}{\partial x_n} 
\end{bmatrix} 
$$
I also know that the hessian for a single function is calculated as:
$$
\mathbf{H}_{f_1}= \begin{bmatrix}
\frac{\partial ^2 f_1}{\partial {x_1}^2} & \frac{\partial ^2 f_1}{\partial {x_1}{x_2}} & \dots & \frac{\partial ^2 f_1}{\partial {x_1}{x_n}}  \\  \frac{\partial ^2 f_1}{\partial {x_2}{x_1}} & \frac{\partial ^2 f_1}{\partial {x_2}^2} & \dots & \frac{\partial ^2 f_1}{\partial {x_2}{x_n}}  \\ \vdots & \vdots & \ddots & \vdots \\ \frac{\partial ^2 f_1}{\partial {x_n}{x_1}} & \frac{\partial ^2 f_1}{\partial {x_n}{x_2}} & \dots & \frac{\partial ^2 f_1}{\partial {x_n}^2}  
\end{bmatrix} 
$$
but I don't have an idea of how does the Hessian for a vector function should look like, neither how to calculate it.
My idea was to calculate the hessian of each function, but I have no idea how to structure the result matrix  
$$
\mathbf{H}_{f_1}, \mathbf{H}_{f_2} , \dots , \mathbf{H}_{f_n}
$$
 A: The hessian of a vector valued function is a 3-tensor, which is simply a trilinear form. So $H \vec{F}=\dfrac{\partial F_i}{\partial x_k\partial x_j}$ and we have that $H \vec{F}(\vec{v},\vec{w},\vec{u})=\sum v_i \vec{u}H F_i \vec{w}$
A: This is something many multivariate calculus courses often miss. Let's recall where the derivative lives. If a function $F\colon V\to W$ where $V$ and $W$ are finite dimensional vector spaces, then $DF\colon V\to L(V,W)$. Using this, let's figure out where the second derivative of a vector function lives. If we have $f\colon \mathbb{R}^{n}\to\mathbb{R}^{n}$ then $DF\colon \mathbb{R}^{n}\to L(\mathbb{R}^{n},\mathbb{R}^{n})$. Thus,  $D^{2}F\colon \mathbb{R}^{n} \to L(\mathbb{R}^{n},L(\mathbb{R}^{n},\mathbb{R}^{n}))$. Thus, $D^{2}F(x)\in L(\mathbb{R}^{n},L(\mathbb{R}^{n},\mathbb{R}^{n}))$.
Thus, the hessian of a vector valued function can be thought of a vector of matrices.  For instance, one can verify that provided $F\in C^{3}$, $D^{2}F(x)\cdot e_{i} = H_{f_{i}}(x)$. 
If you want more education on this matter, I recommend Cartan's Differential Calculus. The book was recently reprinted. 
