Reducing 1D Convolution to a Single (Big) Matrix Multiplication

This is perhaps the 3rd time I’ve needed this recipe and it doesnt seem to be readily available on google.  Theano and Tensorflow provide convolution primitives for 1D and 2D, but (correct me if I’m wrong) I think they are generally constrained such that the filter taps you are convolving must be parameters, and not additional tensor values in a big tensor application.   This is unfortunate, and annoying for certain operations, and my work around is to implement my own convolution as a matrix multiplication based on a properly indexed version of the input and tap tensors within an operation.

Anyway, hopefully this snippet will be useful to someone else some day –

The idea here is simply that we can simply use a toeplitz matrix to generate a large 2D matrix (H) which is simply indexes into a 1D input of taps (h).   Multiplying our input (x) by the 2D (H) matrix then simply gives us our convolution output (y).   Its fairly simple but somewhat tedious to set up, an example implementation is shown below for reference.

#!/usr/bin/env python
import numpy as np
from scipy import linalg
from scipy import signal
x = np.array([0,0,1,0,0,2,0,0,0]) # 9
h = np.array([0,1,2,0]) # 4
y = signal.convolve(x, h, mode='same')
print "x", x
print "h", h
print "y(conv):", y
# set up the toeplitz matrix
padding = np.zeros(len(x)-1, h.dtype)
first_col = np.r_[h, padding]
first_row = np.r_[h[0], padding]
H = linalg.toeplitz(first_col, first_row)[1:len(x)+1,:]
print "shape", H.shape, x.shape
y = np.sum(np.multiply(x,H), 1)
print "y(mult):", y
print "**********************"
x = np.array([0,0,1,0,0,2,0,0,0]) # nsamp
x = np.tile(x,[10,1]) # n_ex x n_samp
h = np.array([0,1,2,0]) # n_samp
h = np.tile(h,[10,1]) # n_ex x n_samp
y = np.zeros([x.shape[0], x.shape[1]])
for i in range(0,x.shape[0]):
    y[i,:] = signal.convolve(x[i,:], h[i,:], mode='same')
print "x", x
print "h", h
print "y(conv):", y
# set up the toeplitz matrix
H = np.zeros([ x.shape[0], x.shape[1], x.shape[1] ]) # n_ex x n_samp x n_samp
for i in range(0,x.shape[0]):
    padding = np.zeros(x.shape[1]-1, h.dtype) #
    first_col = np.r_[h[i,:], padding] #
    first_row = np.r_[h[i,0], padding] #
    H[i,:,:] = linalg.toeplitz(first_col, first_row)[1:x.shape[1]+1,:]
print "H shape", H.shape
print H[0,:,:]
x = x.reshape([x.shape[0], 1, x.shape[1]])
x = np.tile(x, [1,x.shape[1],1])
y = np.sum(np.multiply(x,H), 2)
print "y(mult):", y
print "**********************"
h = np.array([0,1,2,3,4,5,6,7,8], dtype='int32')
padding = np.zeros(len(x)-1, h.dtype)
first_col = np.r_[h, padding]
first_row = np.r_[h[0], padding]
H = linalg.toeplitz(first_col, first_row)[1:len(x)+1,:]
print H

 

Leave a Reply

Your email address will not be published. Required fields are marked *