Discrete 3D convolution of matrix valued functions

As stated in the title, I need to make a discrete 3D convolution of two matrix-valued functions. Or, actually, I need to convolve a matrix-valued function with a vector-valued function. That is, I need to compute the discrete version of $(f \ast X)(x)$ for some finite set of $x$'s, where $f:\mathbb{R}^3\to\mathbb{R}^{3 \times 3}$ and $X:\mathbb{R}^3\to\mathbb{R}^3$. Oh, and I'm using MATLAB.

My question is, how do I do this in the fastest way? I am currently using four nested loops, which, needless to say, is quite slow. Is there any other way?

I have thought about doing it as a multiplication in Fourier-space instead. But is there a risk of loosing some information or precision this way, given that $f$ is actually pretty ugly (an odd, slowly decaying function, with a hole cut in the center) and $X$ is random (as in, is actually a stochastic process)?

I'd be pretty surprised if using Fourier techniques didn't help. Moreover, someone has already written a Matlab code for 3D fast convolution using FFT; that author tested it and found that:

Good usage recommendation:
In 1D, this function is faster than CONV for nA, nB > 1000.
In 2D, this function is faster than CONV2 for nA, nB > 20.
In 3D, this function is faster than CONVN for nA, nB > 5

where nA, nB are the lengths of the input vectors A,B. At the very least, you can test that output against your program to compare the two, without worry of wasting time writing a Fourier-based method yourself and having it work poorly.

That said, you may be right to be cautious about using the Fourier transform of your functions. How are these functions given to you? If they're defined at a discrete set of points to begin with, you needn't worry; the discrete Fourier transform is invertible after all.

However, if you're given some analytic expression that defines them in all of 3-space -- for example, $X$ is a Brownian sheet -- then you might be in trouble. In that case, sampling them at a finite number of points and using the discrete Fourier transform may either discard high-frequency components of the fields or alias them to low frequencies. Lack of smoothness is what hinders Fourier analysis.

Finally, if these concerns do prejudice you against using the FFT, you could try using Matlab's parallel programming toolbox to speed up your nested loop. It's also worth considering the order in which you perform those loops. Matlab stores arrays in column-major order, so it's much faster to access A(:,i) than it is to access A(i,:). Nesting your loops one way as opposed to the other could dramatically slow down a program because it has to take big skips around in memory rather than read entries consecutively.

• First of all, thanks for the link! I have Googled, but haven't stumbled upon that yet. I do have an analytic expression for $f$. However, before I do anything else, I evaluate it in a discrete set of points. Naturally, I choose the resolution small enough, for it not to cause too much trouble. I have a version of MATLAB, licensed by my University, and I don't think it includes the Parallel Programming toolbox.:( Besides, doesn't MATLAB use all cores as default? I already index my arrays as suggested, but thanks anyway! May 3 '13 at 9:57
• This is what I experimented with the MATLAB code. >> tic,A = convnfft(rand(300,300,300), ones(5,5,5), 'same');toc which took 8.061082 seconds. The same thing was used for convn, which is a MATLAB function: >> tic,A = convn(rand(300,300,300), ones(5,5,5), 'same');toc. This one took 2.085360 seconds. MATLAB was found to be about 4 times faster than the code which Bruno Luong has developed. mathworks.com/matlabcentral/fileexchange/… Nov 20 '14 at 23:26