2D Convolution as a Doubly Block Circulant Matrix Operating on a Vector
The point is that circular convolution of two 1-D discrete signals can be expressed as the product of a circulant matrix and the vector representation of the other signal. 

The circulant matrix is a toeplitz matrix which is constructed by different circular shifts of a vector in different rows. For example, consider two signls $h[n]$ and $g[n]$, each of length $4$. If we assume $\mathbf{h}=\begin{bmatrix}
h_0 &amp; h_1 &amp; h_2 &amp; h_3
\end{bmatrix}$, then the circulant matrix denoted by $\mathrm{circ}(h)$ is $$\mathrm{circ}(\mathbf{h})=\mathbf{H}=\begin{bmatrix}
h_0 &amp; h_3 &amp; h_2 &amp; h_1\\ 
h_1 &amp; h_0 &amp; h_3 &amp; h_2\\ 
h_2 &amp; h_1 &amp; h_0 &amp; h_3\\ 
h_3 &amp; h_2 &amp; h_1 &amp; h_0
\end{bmatrix}.$$
So we can calculate $h[n]*g[n]$ by evaluating the product $\mathbf{H}\mathbf{g}$, where $\mathbf{g}=\begin{bmatrix}
g_0 &amp; g_1 &amp; g_2 &amp; g_3
\end{bmatrix}^T$. 

This can be extended to 2-D signals. I explain it through the given example.

The size of signals are $2\times2$ and $3\times3$. The size of convolution is $(2+3-1)\times(2+3-1)=4\times4$. So we need to pad zeros to adjust the size of 2-D signal $x_2(m,n)$:
$$\mathbf{A}=\begin{bmatrix}
 0 &amp; 0 &amp; 0 &amp; 0\\ 
0 &amp; -1 &amp; 0 &amp; 0\\ 
-1 &amp; 4 &amp; -1 &amp; 0\\ 
0 &amp; -1 &amp; 0 &amp; 0
\end{bmatrix}$$
Now construct four circulant matrices corresponding to the four rows of $\mathbf{A}$ denoted by $\mathbf{a}_i$ as follows:

$$\begin{align}
\mathbf{X}_0&amp;=\mathrm{circ}(\mathbf{a}_0)=\mathrm{circ}(\begin{bmatrix}
0 &amp; -1 &amp; 0 &amp; 0
\end{bmatrix})\\
\mathbf{X}_1&amp;=\mathrm{circ}(\mathbf{a}_1)=\mathrm{circ}(\begin{bmatrix}
-1 &amp; 4 &amp; -1 &amp; 0
\end{bmatrix})\\
\mathbf{X}_2&amp;=\mathrm{circ}(\mathbf{a}_2)=\mathrm{circ}(\begin{bmatrix}
0 &amp; -1 &amp; 0 &amp; 0
\end{bmatrix})\\
\mathbf{X}_3&amp;=\mathrm{circ}(\mathbf{a}_3)=\mathrm{circ}(\begin{bmatrix}
0 &amp; 0 &amp; 0 &amp; 0
\end{bmatrix})
\end{align}$$
and construct the $16\times16$ doubly circuland matrix $$\mathbf{X}=\begin{bmatrix}
\mathbf{X}_0 &amp; \mathbf{X}_3 &amp; \mathbf{X}_2 &amp; \mathbf{X}_1\\ 
\mathbf{X}_1 &amp; \mathbf{X}_0 &amp; \mathbf{X}_3 &amp; \mathbf{X}_2\\
\mathbf{X}_2 &amp; \mathbf{X}_1 &amp; \mathbf{X}_0 &amp; \mathbf{X}_3\\
\mathbf{X}_3 &amp; \mathbf{X}_2 &amp; \mathbf{X}_1 &amp; \mathbf{X}_0 
\end{bmatrix}.$$
Now we need to construct the vector corresponding to $x_1(m,n)$. Pad zeros from left and bottom to adjust indices considering $m,n$, and pad zeros from top and right to make the size $4\times4$:
$$\mathbf{B}=\begin{bmatrix}
 0 &amp; 0 &amp; 0 &amp; 0\\ 
0 &amp; 3 &amp; 4 &amp; 0\\ 
0 &amp; 1 &amp; 2 &amp; 0\\ 
0 &amp; 0 &amp; 0 &amp; 0
\end{bmatrix},$$
then construct a $16\times 1$ vector $\mathbf{b}=\begin{bmatrix}
 \mathbf{b}_0 \ |&amp; \mathbf{b}_1 \ |&amp; \mathbf{b}_2 \ |&amp;\mathbf{b}_3 
\end{bmatrix}^T$ by concatenating the rows of matrix $\mathbf{B}$. 

Finally, the 2-D circular convolution is $$\boxed{\mathbf{G}=\mathbf{X}\mathbf{b}}$$