Discrete differentiation with the Kronecker product

The Kronecker product has all kinds of interesting mathematical applications, but is also often extremely handy when creating certain kinds of patterned matrices in R or Matlab. Let’s say you are interested in measuring the smoothness of a function ${f(x,y}$ ), along the ${x}$ and ${y}$ direction. One measure is given by ${\int\left(\frac{d}{d^{2}x}f\left(x,y\right)\right)^{2}\mbox{d}x\mbox{d}y}$ and ${\int\left(\frac{d}{d^{2}x}f\left(x,y\right)\right)^{2}\mbox{d}x\mbox{d}y}$ . Let’s say that you sample ${f}$ along a regular grid. The derivatives can be approximated by a discrete differentation operator. If the sampled ${f}$ is stored as a matrix F (image-like storage), then 2nd order discrete differentiation along the columns can be done via:

$\displaystyle \mathbf{D}=\left[\begin{array}{cccccc} -1 & 2 & 1 & 0 & \ldots & \ldots\\ 0 & -1 & 2 & 1 & 0 & \ldots\\ 0 & 0 & -1 & 2 & 1 & \ldots\\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \end{array}\right]$

D is a bunch of shifted copies of the operator ${\left[\begin{array}{ccccc} -1 & 2 & 1 & 0 & \ldots\end{array}\right]}$ . The product ${\mathbf{DF}}$ corresponds to ${\mathbf{F}}$ , differentiated column-wise. Similarly, ${\mathbf{FD}^{t}}$ will be the row-differentiated version of ${\mathbf{F}}$ . What we are interested in computing is the energy of ${\mathbf{DF}}$ and ${\mathbf{FD}^{t}}$ , which amounts to nothing more than taking each element, squaring it and taking the sum.

What if ${\mathbf{F}}$ is not stored as a matrix, but as a vector ${\mathbf{f}}$ ? The operator ${vec()}$ is mathematical notation for taking a matrix and turning it into a column vector by stacking the columns one below the other, i.e. what happens when you do ${F(:)}$ in Matlab, or ${as.numeric(F)}$ in R. We need an operator E that does the same thing to ${\mathbf{f}}$ that D was doing to F, i.e. discrete differentiation along the ${y}$ axis.

Enters the Kronecker product, and this wonderful equality you’ll find for example in the Matrix cookbook:

$\displaystyle vec(\mathbf{AXB})=\left(\mathbf{B}\otimes\mathbf{A}^{t}\right)vec\left(\mathbf{X}\right)$

${\otimes}$ stands for the Kronecker product here.

Taking ${\mathbf{X}=\mathbf{F}}$ , ${\mathbf{A}=\mathbf{F}}$ , and ${\mathbf{B}=\mathbf{I}}$ , we obtain:

$\displaystyle vec(\mathbf{DF})=\left(\mathbf{I}\otimes\mathbf{D}^{t}\right)\mathbf{f}$

${\left(\mathbf{I}\otimes\mathbf{D}^{t}\right)}$ is the operator we were looking for. This is extremely handy because such operators crop up a lot in Gauss-Markov random fields and penalised splines methods, and this gives a very fast way of creating the right precision matrices. More on Kronecker products in this article by C.F. van Loan.

dahtah