FFT for distance-limited sum of pairwise products

Question

I have a sequence of fixed length $N$ input arrays ($x_0$, $x_1$, ..., $x_{n-1}$), and my goal is to compute output arrays as such:

$$d_{00}[l] = \sum_{i = 0}^{N - l} x_0[i]x_0[i + l]$$ $$d_{01}[l] = \sum_{i = 0}^{N - l} x_0[i]x_1[i + l]$$ $$d_{02}[l] = \sum_{i = 0}^{N - l} x_0[i]x_2[i + l]$$ $$...$$ $$d_{10}[l] = \sum_{i = 0}^{N - 1} x_1[i]x_0[i + l]$$ $$...$$ With $l$ having a maximum value that scales linearly in $N$. That is to say, $d_{ij}$ will contain, in it's $k^{th}$ entry, the sum of all pairs in $x_i$ and $x_j$ that are distance $k$ away.

My understanding - Clearly, since each $d_{ij}$ output requires $l \propto N$ entries, each containing $N$ multiplies, the time complexity of producing each $d_{ij}$ output array is $O(N^2)$. These pairwise multiplies hidden behind a summation seem optimal for a Fourier Transform to improve the time complexity to $O(N \log N)$. However, I understand the standard Fourier Transform convolution does not allow for limiting to constant distance, instead allowing easy sum of pairwise indices in the form $\sum_{i=0}^{c} x[i] y[c - i]$.

How can I use this standard convolution, or any variation of the Fourier transform algorithm to make this problem computationally tractable?

Answer 1

Your sequences are all correlations between $x_i$ and $x_j$, where you assume that $x_i[n]=0$ for $n<0$ or $n\geq N$. Correlation is very related to convolution and it has a similar expression in frequency domain:

$$ c[n]=\sum_l x_1[l]x_2[l+n] \Leftrightarrow C[k]=X_1^*[k]X_2[k]$$

(capital letters are Fourier transforms of the small ones).

Hence, you can efficiently calculate the correlation using this theorem.

One caveat is that you are actually only looking for linear correlation, but using DFT you calculate circular correlation. So, the solution is to zero-pad the sequences. Have a look at the following Python code:

N = 5
x1 = np.random.randn(N)
x2 = np.random.randn(N)

x1z = np.hstack([x1, np.zeros(N)])
x2z = np.hstack([x2, np.zeros(N)])

# calculate according to your expression
d12 = []
for l in range(N):
    d12.append(sum(x1[i]*x2[i+l] for i in range(N-l)))

# calculate using zero-padding of the sequences with no 
# dependence of the range of i
d12z = []
for l in range(N):
    d12z.append(sum(x1z[i]*x2z[i+l] for i in range(N)))

# calculate using circular correlation, using the zero-padded versions
d12c = []
for l in range(N):
    d12c.append(sum(x1z[i]*x2z[(i+l)%len(x2z)] for i in range(len(x1z))))
# check they all match   
print ("x1:", x1)
print ("x2:", x2)
print ("d12_1:", d12)
print ("d12_2:", d12z)
print ("d12_3:", d12c)

# now, exploit the correlation theorem to calculate the same stuff:
C = (np.fft.ifft(np.fft.fft(x1z)*np.fft.fft(x2z).conj()).real)
# take only the first N values, as the others correspond to parts 
# where the sequences do not fully overlap.
print ("FFT:", C[-np.arange(N)])  

A sample output could then be:

x1: [ 1.11998979 -0.81341718  1.47117346  0.59699824 -0.9233909 ]
x2: [-0.31884542 -0.2214878  -0.02678967 -0.25943789  1.36355875]
d12_1: [-1.6303355800547898, 0.20609114482882795, 2.1870585244152387, -1.3997099032157734, 1.5271718773447336]
d12_2: [-1.6303355800547898, 0.20609114482882795, 2.1870585244152387, -1.3997099032157734, 1.5271718773447336]
d12_3: [-1.6303355800547898, 0.20609114482882795, 2.1870585244152387, -1.3997099032157734, 1.5271718773447336]
FFT: [-1.63033558  0.20609114  2.18705852 -1.3997099   1.52717188]