2D matrix convolution with convolve() in R

We could implement the sliding window approach in a for loop. First, pad largemx with zeroes, then take the sum of the window multiplied by the kernel.

> ## zero-pad largemx
> nd <- dim(largemx)
> kd <- dim(kernel)
> pd <- kd/2
> pad <- matrix(0, nd[1] + 2*pd[1] - 1, nd[2] + 2*pd[2] - 1)
> pad[2:(1 + nd[1]), 2:(1 + nd[2])] <- largemx
> pad
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    0    0    0    0    0    0
[2,]    0    1    2    3    4    0
[3,]    0    5    6    7    8    0
[4,]    0    9   10   11   12    0
[5,]    0   13   14   15   16    0
[6,]    0    0    0    0    0    0
> ## sliding window multiplication
> res <- matrix(0, nd[1], nd[2])
> for (i in 1:nd[1]) {
+   for (j in 1:nd[2]) {
+     res[i, j] <- sum(pad[i:(i + 2), j:(j + 2)]*kernel)
+   }
+ }
> res
     [,1] [,2] [,3] [,4]
[1,]   13   22   27   18
[2,]   28   48   56   37
[3,]   48   80   88   57
[4,]   33   58   63   38

Probably faster using Rcpp:

> cpp <- '
+ NumericMatrix conv_2d(NumericMatrix x, NumericMatrix k) {
+   int nr = x.nrow(), nc = x.ncol();
+   int kr = k.nrow(), kc = k.ncol();
+   int pr = kr/2, pc = kc/2;
+   NumericMatrix y(nr, nc);
+   NumericMatrix pad(nr + 2*pr, nc + 2*pc);
+   for (int i = 0; i < nr; i++) {
+     for (int j = 0; j < nc; j++) {
+       pad(i + pr, j + pc) = x(i, j);
+     }
+   }
+   for (int i = 0; i < nr; i++) {
+     for (int j = 0; j < nc; j++) {
+       double acc = 0.0;
+       for (int a = 0; a < kr; a++) {
+         for (int b = 0; b < kc; b++) {
+           acc += pad(i + a, j + b)*k(kr - 1 - a, kc - 1 - b);
+         }
+       }
+       y(i, j) = acc;
+     }
+   }
+   return y;
+ }
+ '
> Rcpp::cppFunction(cpp)
> conv_2d(largemx, kernel)
     [,1] [,2] [,3] [,4]
[1,]   13   22   27   18
[2,]   28   48   56   37
[3,]   48   80   88   57
[4,]   33   58   63   38

The link also mentions a matrix approach which might be even better.

Not sure how this scales but works for this example. Might be packages out there using a similar approach, but sometimes DIY implementing is faster than finding a package :)

Data:

largemx <- matrix(1:16, nrow=4, byrow=TRUE)
kernel <- array(replace(rep(1, 3*3), 5, 0), c(3, 3))

I have found an easy and compact solution how to apply convolve for this case. This is a bit hacky but works for my particular task.

This is the input matrix once again:

largemx <- matrix(1:16, nrow=4, byrow = T)

apply the filter c(1, 1, 1) once to the rows padded with zeroes:

(convoluted_rows <- apply(largemx, 1, \(x)convolve(c(0, x, 0), c(1, 1, 1), type="f")))
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    0    3   11   19   27    0
[2,]    0    6   18   30   42    0
[3,]    0    9   21   33   45    0
[4,]    0    7   15   23   31    0

And then again to the result:

(convoluted_twice <- apply(convoluted_rows, 1, \(x)convolve((0, x, 0), c(1, 1, 1), type="f")))
     [,1] [,2] [,3] [,4]
[1,]   14   24   30   22
[2,]   33   54   63   45
[3,]   57   90   99   69
[4,]   46   72   78   54

It only remains to subtract the original matrix:

(convoluted_twice-largemx)
     [,1] [,2] [,3] [,4]
[1,]   13   22   27   18
[2,]   28   48   56   37
[3,]   48   80   88   57
[4,]   33   58   63   38

How this works? (in case you were wondering)

As it is written in the description to convolve fuction, "Use the Fast Fourier Transform to compute the several kinds of convolutions of two sequences." So it is designed to operate on 1D and not on 2D objects. When you are feeding 2D arrays to convolve, they get flattened, which is unwanted.

In my solution I have started by convolving the rows (by the way, it will work with the columns the same way, but you need to transpose the result) with the kernel c(1, 1, 1). It sums up the rows by the groups of consequent three numbers: Σ[0, 1, 2] = 3, Σ[1, 2, 3] = 6, Σ[2, 3, 4] = 9 etc. Note that apply gives a transposed output, but since it is used twice, that doesn't matter.

In the second run the same kernel is applied to the result. Now the columns are summed up by groups of three. In the end this produces a matrix where each element is a sum of the corresponding 3x3 sumbatrix of pad (see the jay.sf's answer). Because the host number of its cell in the matrix should not be included, we simply subtract the original largemx from the result.