How to use Fast Fourier Transform to execute convolution of matrix?

Question

I need to add many big 3D arrays (with a shape of 500x500x500) together and want to speed up the process by using multiplication in the Fourier space. The problem is that I don't get the same answer when multiplying in the Fourier space compared to simply adding the matrix.

To test it out, I wrote a minimal example trying to make it work but the answer is not what I expected. Either my math knowledge is wrong or I am not using the function correctly.

Below is the simplest code showing what I am trying to do:

import numpy as np

c = np.asarray(((1,2),(2,3)))
d = np.asarray(((1,4),(1,5)))

print("Transform")
Nc = np.fft.rfft2(c)
Nd = np.fft.rfft2(d)

print("Inverse")
Nnc = np.fft.irfft2(Nc)
Nnd = np.fft.irfft2(Nd)

print("Somme")
S = np.dot(Nc, Nd)
print(np.fft.irfft2(S))

When I print S, I get the result:

[[6, 28],[10,46]]

But from what I understood about the Fourier space, multiplication would mean addition outside of the Fourier space so I should get S = c + d?

Am I doing something wrong using the FFT function or is my assumption that S should equal c plus d wrong?

Answer 1

There is a little misunderstanding here:

Multiplication in Fourier space corresponds to convolution in the spatial domain and not to addition.

There is no way to speed up addition in that way.

Answer 2

If you want to compute c+d through the Fourier domain, you'd have to add the two spectra, not multiply them:

np.fft.irfft2(Nc+Nd) == c+d  # (up to numerical precision)

Of course, this is much slower than simply adding the matrices in the spatial domain.

As @Florian said, it is convolution that can be sped up by multiplying in the spatial domain.