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.