# How can the convolution operation be implemented as a matrix multiplication?

How can the convolution operation used by CNNs be implemented as a matrix-vector multiplication? We often think of the convolution operation in CNNs as a kernel that slides across the input. However, rather than sliding this kernel (e.g. using loops), we can perform the convolution operation "in one step" using a matrix-vector multiplication, where the matrix is a circulant matrix containing shifted versions of the kernel (as rows or columns) and the vector is the input.

How exactly can this operation be performed? I am looking for a detailed step-by-step answer that shows how the convolution operation (as usually presented) can be performed using a matrix-vector multiplication.

Is this the usual way the convolution operations are implemented in CNNs?

To show how the convolution (in the context of CNNs) can be viewed as matrix-vector multiplication, let's suppose that we want to apply a $$3 \times 3$$ kernel to a $$4 \times 4$$ input, with no padding and with unit stride.

Here's an illustration of this convolutional layer (where, in blue, we have the input, in dark blue, the kernel, and, in green, the feature map or output of the convolution).

Now, let the kernel be defined as follows

$$\mathbf{W} = \begin{bmatrix} w_{0, 0} & w_{0, 1} & w_{0, 2} \\ w_{1, 0} & w_{1, 1} & w_{1, 2} \\ w_{2, 0} & w_{2, 1} & w_{2, 2} \end{bmatrix} \in \mathbb{R}^{3 \times 3}$$

Similarly, let the input be defined as

$$\mathbf{I} = \begin{bmatrix} i_{0, 0} & i_{0, 1} & i_{0, 2} & i_{0, 3} \\ i_{1, 0} & i_{1, 1} & i_{1, 2} & i_{1, 3} \\ i_{2, 0} & i_{2, 1} & i_{2, 2} & i_{2, 3} \\ i_{3, 0} & i_{3, 1} & i_{3, 2} & i_{3, 3} \\ \end{bmatrix} \in \mathbb{R}^{4 \times 4}$$

Then the convolution above (without padding and with stride 1) can be computed as a matrix-vector multiplication as follows. First, we redefine the kernel $$\mathbf{W}$$ as a sparse matrix $$\mathbf{W}' \in \mathbb{R}^{4 \times 16}$$ (which is a circulant matrix because of its circular nature) as follows.

$${\scriptscriptstyle \mathbf{W}' = \begin{bmatrix} w_{0, 0} & w_{0, 1} & w_{0, 2} & 0 & w_{1, 0} & w_{1, 1} & w_{1, 2} & 0 & w_{2, 0} & w_{2, 1} & w_{2, 2} & 0 & 0 & 0 & 0 & 0 \\ 0 & w_{0, 0} & w_{0, 1} & w_{0, 2} & 0 & w_{1, 0} & w_{1, 1} & w_{1, 2} & 0 & w_{2, 0} & w_{2, 1} & w_{2, 2} & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & w_{0, 0} & w_{0, 1} & w_{0, 2} & 0 & w_{1, 0} & w_{1, 1} & w_{1, 2} & 0 & w_{2, 0} & w_{2, 1} & w_{2, 2} & 0 \\ 0 & 0 & 0 & 0 & 0 & w_{0, 0} & w_{0, 1} & w_{0, 2} & 0 & w_{1, 0} & w_{1, 1} & w_{1, 2} & 0 & w_{2, 0} & w_{2, 1} & w_{2, 2} \end{bmatrix} }$$ Similarly, we reshape the input $$\mathbf{I}$$ as a 16-dimensional vector $$\mathbf{I}' \in \mathbb{R}^{16}$$.

$${\scriptstyle \mathbf{I}' = \begin{bmatrix} i_{0, 0} & i_{0, 1} & i_{0, 2} & i_{0, 3} & i_{1, 0} & i_{1, 1} & i_{1, 2} & i_{1, 3} & i_{2, 0} & i_{2, 1} & i_{2, 2} & i_{2, 3} & i_{3, 0} & i_{3, 1} & i_{3, 2} & i_{3, 3} \end{bmatrix}^T }$$

Then the convolution of $$\mathbf{W}$$ and $$\mathbf{I}$$, that is

$$\mathbf{W} \circledast \mathbf{I} = \mathbf{O} \in \mathbb{R}^{2 \times 2},$$ where $$\circledast$$ is the convolution operator, is equivalently defined as $$\mathbf{W}' \cdot \mathbf{I}' = \mathbf{O}' \in \mathbb{R}^{4},$$ where $$\cdot$$ is the matrix-vector multiplication operator. The produced vector $$\mathbf{O}'$$ can then be reshaped as a $$2 \times 2$$ feature map.

You can easily verify that this representation is correct by multiplying e.g. the 16-dimensional input vector $$\mathbf{I}'$$ with the first row of $$\mathbf{W}'$$ to obtain the top-left entry of the feature map.

$$w_{0, 0} i_{0, 0} + w_{0, 1} i_{0, 1} + w_{0, 2} i_{0, 2} + 0 i_{0, 3} + w_{1, 0} i_{1, 0} + w_{1, 1} i_{1, 1} + w_{1, 2}i_{1, 2} + 0 i_{1, 3} + w_{2, 0} i_{2, 0} + w_{2, 1}i_{2, 1} + w_{2, 2} i_{2, 2} + 0 i_{2, 3} + 0 i_{3, 0} + 0 i_{3, 1} + 0 i_{3, 2} + 0 i_{3, 3} = \\ w_{0, 0} i_{0, 0} + w_{0, 1} i_{0, 1} + w_{0, 2} i_{0, 2} + w_{1, 0} i_{1, 0} + w_{1, 1} i_{1, 1} + w_{1, 2}i_{1, 2} + w_{2, 0} i_{2, 0} + w_{2, 1}i_{2, 1} + w_{2, 2} i_{2, 2} = \\ \mathbf{O}'_{0} \in \mathbb{R} ,$$ which is equivalent to an element-wise multiplication of $$\mathbf{W}$$ with the top-left $$3 \times 3$$ sub-matrix of the input followed by a summation over all elements (i.e. convolution), that is

$$\sum \left( \begin{bmatrix} w_{0, 0} & w_{0, 1} & w_{0, 2} \\ w_{1, 0} & w_{1, 1} & w_{1, 2} \\ w_{2, 0} & w_{2, 1} & w_{2, 2} \end{bmatrix} \odot \begin{bmatrix} i_{0, 0} & i_{0, 1} & i_{0, 2} \\ i_{1, 0} & i_{1, 1} & i_{1, 2} \\ i_{2, 0} & i_{2, 1} & i_{2, 2} \end{bmatrix} \right) = \mathbf{O}_{0, 0} = \mathbf{O}'_{0} \in \mathbb{R},$$ where $$\odot$$ is the element-wise multiplication and $$\sum$$ is the summation over all elements of the resulting matrix.

The advantage of this representation (and computation) is that back-propagation can be computed more easily by just transposing $$\mathbf{W}'$$, i.e. with $$\mathbf{W}'^T$$.

Take also a look at this Github repository that explains how the convolution can be implemented as matrix multiplication.