# Methods for fast approximation of convolution

What are the state of the art methods for fast 2D convolution approximation?
I'm familiar with SVD based multiplication and cross approximation approaches, but would be thankful to get additional references.

The kernel size is typically ~3x3 - 11x11, not separable, with input matrices from sizes ~200x200 - 1000x1000.

Thanks!

• Do you are a repeated object in the convolution, i.e. a kernel that is used everytime? – Laurent Duval Mar 10 '16 at 16:22
• Is the kernel that you're convolving with the same size as the thing it is being convolved with or smaller? How big are they? – Brian Borchers Mar 11 '16 at 4:57
• @LaurentDuval - I don't use repeated objects. Kernel sizes vary from 3x3 to 11x11, and input matrices are ~200x200-1000x1000. Thanks. – rursw1 Mar 13 '16 at 7:25
• @BrianBorchers - no, the kernel is much smaller than the input. Thanks. – rursw1 Mar 13 '16 at 7:25
• One might use the ideas of compressed sensing the sense of very sparse cloud of Dirichlet samples. I wonder if there is a sample-at-a-time approach that allows stopping once the resultant change is below a threshold. That might be fast. – EngrStudent Mar 17 '16 at 16:36

Fast filter approximations have been studied for a long time, especially to implement IIR filters, like Gaussians and their derivatives. You may want to reuse such concepts, with keywords like integral image, summed-area tables, box filter, recursive filtering. You can start from a recent review:

and dig other cited papers.

I am thinking about approximating an impulse response as a sum of boxes, on top of the others. The next one could be of interest "We have presented a method for automatically approximating an arbitrary 2-D filter by a box filter"

I am not sure though you can gain a lot in number of operations, at least for separable filters.

Disclaimer:
I'm new(er) to 'R' and so I don't know the sparse FFT library (SE,MIT,Berk) off the top of my head. This stuff has been published for ~2 years. I would be surprised if it doesn't exist - it should and would be valuable addition to the language if it doesn't. If I were highly motivated (not at the moment), I would make it myself.

The question is:
You are wanting to convolute a, presumably zero-padded, 3x3 to 11x11 matrix with a 200x200 to 1000x1000 matrix very quickly.

Disclaimer:
I'm not "hip" on the state of the art on this, sorry.

My comment from above:
One might use the ideas of compressed sensing the sense of very sparse cloud of Dirichlet samples. I wonder if there is a sample-at-a-time approach that allows stopping once the resultant change is below a threshold. That might be fast.

What I mean by compressed sensing from sparse Dirichlet samples is this:
Uniformly randomly sample single pixels across the large picture and one-at-a-time add to your convolution result.

Let's say that you are looking for "pac-man" in the screen image.

Here is pac-man

Here is the game-screen image:

How do you find the "pac-man" in the image using convolution?

Here is rev0 code:

library(png)             #image load
library(pracma)          #matlab like matrix utils
library(microbenchmark)  #for very tight timing
library(spatstat)        #for gaussian blur

#raw data
im1 <- readPNG(source = "./data/pacman1.png",native=F) #pac man
im2 <- readPNG(source = "./data/pac_game.png")         #game board

#convert to grayscale (aka 2d matrix)
im1.g <- rot90(0.2126*im1[,,1] + 0.7152*im1[,,2] + 0.0722*im1[,,3],k=-1)
im2.g <- rot90(0.2126*im2[,,1] + 0.7152*im2[,,2] + 0.0722*im2[,,3],k=-1)

#normalize
im1.g <- im1.g/max(im1.g)
im2.g <- im2.g/max(im2.g)

#housekeeping on original
rm(im1,im2)

# #check images
# image(im1.g,col = gray.colors(256))
# image(im2.g,col = gray.colors(256))

#classic method for image registration using convolution

#make low intensity values negative,
# (improves curvature for subpixel approaches.)
im1.g <- im1.g - 0.1
im2.g <- im2.g - 0.1

im1b <- 0 * im2.g
temp <- size(im1.g)

im1b[1:temp[1],1:temp[2]] <- im1.g

#convolute
fim1 <- fft(im1b)
fim2 <- fft(im2.g)
im12 <- Re(fft(fim1*fim2,inverse=T))

#ground negatives
im12[which(im12<0)] <- 0

#scale to height
im12 <- im12/max(im12)

#plot the result
image(im12, col = terrain.colors(256), new = TRUE)

drape.plot(1:nrow(im12), 1:ncol(im12), (im12), border = NA, theta = 25, phi = 55, )


The resulting 2d image is:

The drape plot is:

While I don't exactly like what is going on, I am reasonably confident that the one of the highest 6 peaks happens at the actual "packman". Personally I would prefer to use the edges with a gaussian blur so that I am not confusing "Speedy" with Pac-Man, but you are caring about compute speed, not necessarily registration.

So I can wrap the "convolution" part in benchmark and get compute time

# benchmark convolute
mybench <- microbenchmark({
fim1 <- fft(im1b);
fim2 <- fft(im2.g);
im12 <- Re(fft(fim1*fim2,inverse=T));},
times=100)
#display value
print(mybench)


The displayed results are:

Unit: milliseconds
...
min       lq     mean   median       uq      max neval
9.33394 9.665923 10.00031 9.816518 9.940719 14.81459   100


Now there are 64512 pixels in our image. For larger images we could have millions of pixels. Performing a convolution there can be really expensive.

Let's uniformly randomly sub-select 10% of the pixels, and perform the convolution on them, with a sparse fft transform.

#make empty
im2b <- matrix(0, nrow = nrow(im2.g), ncol = ncol(im2.g))

#uniformly randomly sample 10% of pixels
idx2b <- sample(x = 1:length(im2.g), size = floor(length(im2.g)/10),replace = FALSE)
im2b[idx2b] <- im2.g[idx2b]

fim1 <- fft(im1b);
fim2 <- fft(im2b);
im12 <- Re(fft(fim1 * fim2, inverse = T));


With only 10% of pixels, the registration property of the convolution actually improves.

The true location retains its high value in the convolution even when much of the data is removed.

--> INSERT UPDATED BENCHMARK from SPARSE FFT HERE <-- note: the following is work in progress

We can consider each pixel of the sub-sampled 'im2.g' as its own image, the sum of which comprise the full background.

What we then look for is a way to make successive convolutions. The first way may be more expensive, but then we want to find a way to make it less expensive. I know symbolically that we can treat it as a dirac-delta. The forward transform is going to be very cheap. I know that a spatial tight variance function has a wide variance in wavenumber.

• @EngrStudent - could you please explain the "one at a time" approach? I don't find it in your R code. - Thanks. – NewB Apr 7 '16 at 11:50
• @NewB - Several years ago (2007?) there was some excitement and rumors that Gilbert Strang had "cracked the code" on a very fast matrix inverse but it had been classified and he wasn't allowed to publish. There was one obscure comment attributed to him that said a variation on the Woodbury matrix identity for 2x2 diagonally dominant matrices should be investigated. Woodbury can be constructed to update the inverse efficiently for a few changes in a matrix. Let me work on formulation for the one-at-a-time here. (it is a to-do item) – EngrStudent Apr 8 '16 at 16:40
• @EngrStudent - thank you. Very interesting approach! – NewB Apr 10 '16 at 7:10
• @EngrStudent - maybe this is a naïve question, or I missed something. By parsifying the input (uniquely randomly subsample the input image) you mean that this is the measurement stage? i.e. this is the matrix which you can recover the original picture from? Or this is just a way to make sure that the signal is sparse and the compressed sensing idiom comes later? Because I'm pretty sure that for downsampling both input image and kernel you should preserve registration quality in a faster way than sparse FFT. What am I missing? Thanks! – NewB May 10 '16 at 8:34
• @NewB - A photo on an iPhone 6s is 12 megapixels. If you take two pictures and want to perform "image registration" then each transform and each error computation is operating on 12 million pixels. It turns out that to "get in the neighborhood" you don't need all 12 million of them. A million really is huge. You can get away with a few thousand, but they have to be "well behaved". The uniform random helps with that. The cost of a few thousand in the transform and a few thousand in the cost function is thousands of times smaller than 12 million. This comes before sparse FFT. – EngrStudent May 10 '16 at 10:58