Hessian of a $\mathbb{R}^{n \times m} \to \mathbb{R}$ function I have a function $f : \mathbb{R}^{n \times m} \to \mathbb{R}$ whose input is denoted by $\bf U$. I know that $\nabla_{\mathbf{U}}f \in \mathbb{R}^{n \times m}$. I would like to compute $\nabla_{\mathbf{U}}^2 f$, but am currently stuck. For example, I would like to compute a derivative of $\mathbf{A}\mathbf{U}\mathbf{B}$ with respect to $\mathbf{U}$, where $\mathbf{A} \in \mathbb{R}^{n \times n}$ and $\mathbf{B} \in \mathbb{R}^{m \times m}$.
I tried googling some relevant materials (such as The Matrix Cookbook), but was not able to find useful information. Please provide me with some enlightenment to tackle this issue. Thank you!
 A: The gradient of a matrix-valued function with respect to a matrix variable is a 4th order tensor. There are two ways to work with such problems. 
The first is to introduce a special isotropic 4th order tensor, whose components can be expressed in terms of Kronecker deltas
$$\eqalign{
 {\mathcal H}_{ijkl} &= {\delta}_{ik} {\delta}_{jl} \cr
}$$
The utility of this tensor is that it allows you to rearrange matrix products, e.g. in index notation
$$\eqalign{
F_{mj}
 &= A_{mi}{\mathcal H}_{ijkl}B_{ln}^T\,\,U_{kn} \cr
 &= A_{mi}{\delta}_{ik} {\delta}_{jl}B_{nl}\,U_{kn} \cr
 &= A_{mi}U_{in}B_{nj} \cr
}$$
or in matrix notation
$$\eqalign{
F &= AUB &= A{\mathcal H}B^T:U \cr
dF &= A\,dU\,B &= A{\mathcal H}B^T:dU \cr
\frac{\partial F}{\partial U} &= A{\mathcal H}B^T \cr
}$$
The other approach is to use the vec-operation to flatten matrices into vectors. 
$$\eqalign{
{\rm vec}(dF) &= {\rm vec}(A\,dU\,B) \cr
 &= (B^T\otimes A)\,{\rm vec}(dU) \cr
\frac{\partial{\,\rm vec}(F)}{\partial{\,\rm vec}(U)} &= B^T\otimes A\cr
}$$
