Mixed directional derivative of second order and the Hessian matrix First order directional derivative in 2 dimensional $xy$-coorinate system is
$$D_{\bf{u}}=u_1\frac{\partial}{\partial x}+u_2\frac{\partial}{\partial y}$$
I thought to go further and analyse second order directional derivatives.
The formula for general 
second order mixed directional derivative is
$$ D^2_{\bf{u}\bf{v}}=D_{\bf{v}}\left[u_1\frac{\partial}{\partial x}+u_2\frac{\partial}{\partial y}\right]=u_1D_{\bf{v}}\frac{\partial}{\partial x}+u_2D_{\bf{v}}\frac{\partial}{\partial y}=$$
$$=u_1(v_1\frac{\partial^2}{\partial x^2}+v_2\frac{\partial^2}{\partial y\partial x})+u_2(v_1\frac{\partial^2}{\partial x\partial y}+v_2\frac{\partial^2}{\partial y^2})=$$
$$=u_1v_1\frac{\partial^2}{\partial x^2}+u_1v_2\frac{\partial^2}{\partial y\partial x}+u_2v_1\frac{\partial^2}{\partial x\partial y}+u_2v_2\frac{\partial^2}{\partial y^2}$$
I see the formula we get has definitely something to do with the following two matrices:
$$U=\left[\matrix{u_1v_1 && u_1v_2\\ u_2v_1 && u_2v_2}\right]$$
$$H=\left[\matrix{\frac{\partial^2}{\partial x^2} && \frac{\partial^2}{\partial x\partial y} \\ \frac{\partial^2}{\partial y\partial x} && \frac{\partial^2}{\partial y^2}}\right]$$
Where the second of them is known as the Hessian operator in the 2-dimensional $xy$-plane. 
Are here any matrix algebra masters, who will find how to express $D_{\bf{u}\bf{v}}$ in terms of $U$ and $H$ ? 
If we had defined an operation like $$\left[\matrix{a_{11} && a_{12} \\ a_{21} && a_{22}}\right]* \left[\matrix{b_{11} && b_{12} \\ b_{21} && b_{22}}\right]=a_{11}b_{11}+a_{12}b_{12}+a_{21}b_{21}+a_{22}b_{22}$$
$$A*B=\sum_{i, j} A_{ij}B_{ij},$$
then the formula for mixed directional 2nd order derivative would look like
$$D_{\bf{u}\bf{v}}=U^T*H$$
Any ideas how to express the $*$ operation using standard operations on matrices?

I'm not sure if this is exactly what you're looking for, but if we write $\mathbf u=\begin{bmatrix} u_1\\u_2\end{bmatrix},\ \mathbf v=\begin{bmatrix} v_1\\v_2\end{bmatrix}$, then the formula becomes
$$D^2_{\bf{u}\bf{v}} = \mathbf v^T\mathbf H\mathbf u$$
Let us write $$\begin{align}
f_{11}&=\frac{\partial^2f}{\partial x^2}\\
f_{12}&=\frac{\partial^2f}{\partial y\partial x}\\
f_{21}&=\frac{\partial^2f}{\partial x\partial y}\\
f_{22}&=\frac{\partial^2f}{\partial y^2}
\end{align}$$
as usual.  If the partials exist and are continuous, we have $f_{12}=f_{21}$ so we can write $$\mathbf H=\begin{bmatrix}f_{11}&f_{12}\\f_{21}&f_{22}\end{bmatrix}$$
and the formula takes the easily-rembered form
$$v_1f_{11}u_1+v_1f_{12}u_2+v_2f_{21}u_1+v_2f_{22}u_2=\begin{bmatrix}v_1&v_2\end{bmatrix}\begin{bmatrix}f_{11}&f_{12}\\f_{21}&f_{22}\end{bmatrix}\begin{bmatrix}u_1\\u_2\end{bmatrix}$$ Since the left-hand side is a scalar, and $\mathbf H$ is symmetric, we can transposes to get$$D^2_{\bf{u}\bf{v}} = \mathbf u^T\mathbf H\mathbf v$$ so we don't even have to remember which side the $\mathbf u$ goes on.
