5

Discrete convolution can be performed via the Toeplitz matrix, as shown below (Wiki article):

enter image description here

Note that this is not the exact same form as as the general Toeplitz matrix, but it has experienced various shifts and zero-paddings.

Is there a way to achieve this in numpy purely based on roll, hstack etc., i.e. without using any for loops? I have tried all sorts of shifts but I can't really get it in to the form shown above.

ali_m
  • 71,714
  • 23
  • 223
  • 298
Francisco Vargas
  • 653
  • 1
  • 7
  • 22

1 Answers1

10

Yes, you can use scipy.linalg.toeplitz:

import numpy as np
from scipy import linalg

h = np.arange(1, 6)

padding = np.zeros(h.shape[0] - 1, h.dtype)
first_col = np.r_[h, padding]
first_row = np.r_[h[0], padding]

H = linalg.toeplitz(first_col, first_row)

print(repr(H))
# array([[1, 0, 0, 0, 0],
#        [2, 1, 0, 0, 0],
#        [3, 2, 1, 0, 0],
#        [4, 3, 2, 1, 0],
#        [5, 4, 3, 2, 1],
#        [0, 5, 4, 3, 2],
#        [0, 0, 5, 4, 3],
#        [0, 0, 0, 5, 4],
#        [0, 0, 0, 0, 5]])
ali_m
  • 71,714
  • 23
  • 223
  • 298
  • would this work for 2x2 numpy arrays as input signals ? – Francisco Vargas Dec 30 '15 at 20:44
  • No, `scipy.linalg.toeplitz` treats all of its inputs as 1D arrays. I'm not sure what a 2x2 input array would mean in this context, anyway. – ali_m Dec 30 '15 at 20:51
  • convolution between an image and a an image kernel for example. – Francisco Vargas Dec 30 '15 at 20:53
  • which apparently can be achieved in matlab through this method http://stackoverflow.com/questions/16798888/2-d-convolution-as-a-matrix-matrix-multiplication – Francisco Vargas Dec 30 '15 at 20:56
  • 1
    As in that MATLAB question, if your kernel is x-y separable then you could express it as two 1D vectors of weights, then construct two separate Toeplitz matrices as above and compute two separate dot products along different axes of your image array. – ali_m Dec 30 '15 at 20:59