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!

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 }