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.