I want to compute the discrete convolution of two 3D arrays: $A(i, j, k) \ast B(i, j, k)$
Is there a general way to decompose the array $A$ into a sum of a small number of separable arrays? That is:
$A(i, j, k) = \sum_m X_m(i) \times Y_i(j) \times Z_i(k)$
If not, is there a standard way to approximate $A$ as a sum of a relatively small number of separable arrays? $A$ is typically smooth, with a central maximum that dies off towards the edges.
Motivation: I want to compute many iterations of $A \ast B$, with $A$ constant but $B$ changing. $A$ is quite large (16 x 16 x 16, for example), but $B$ is larger still (64 x 1024 x 1024, for example). Direct convolution is very slow to compute. FFT-based convolution is much faster, but still takes tens of seconds for one FFT of $B$ on my machine, and uses a lot of memory. If $A$ is separable, however, three direct convolutions is much faster and uses less memory:
import time, numpy from scipy.fftpack import fftn, ifftn from scipy.ndimage import gaussian_filter a = numpy.zeros((16, 16, 16), dtype=numpy.float64) b = numpy.zeros((64, 1024, 1024), dtype=numpy.float64) fftn_a = fftn(a, shape=b.shape) start = time.clock() ifftn(fftn_a * fftn(b)).real end = time.clock() print "fft convolution:", end - start start = time.clock() gaussian_filter(b, sigma=3) end = time.clock() print "separable convolution:", end - start
fft convolution: 49 seconds
separable convolution: 9.2 seconds
Alternatively, if someone's got a better way to compute this type of convolution, I'm all ears.