The standard algorithm for discrete convolution of two vectors $x\in \mathbb{R}^{n}$ and $y \in \mathbb{R}^{m}$ is (in essence) a FFT of the two input vectors, multiplication of the two elementwise, and then an inverse FFT of the result. In practice there are complications, for example if $m \ne n$ then we need might need some padding of one of the inputs, and padding of the outputs.
From the following graph, you can see that the length of the zero padding affects the results, especially at the boundary. I have searched through the literature, but I could not find anything on the error analysis of the fast convolution algorithm (like Higham's analysis of the FFT in Accuracy and Stability). What error level should I expect for the fast convolution algorithm, and should the padding have this large of an affect?
For reference, the code to generate the figure is here:
#!/usr/bin/env python3
import numpy as np
from numpy import fft as np_fft
from matplotlib import pyplot as plt
def convolve_by_fft(a, b, len_after_padding=None):
# scipy.fftconvolve doesn't allow control of the padding length,
# so make this function mimics what is done in `scipy.fftconvolve`
# i.e., ifft(fft(a) * fft(b))
fft, ifft = np_fft.rfftn, np_fft.irfftn
result_len = a.size + b.size - 1
if len_after_padding is None:
len_after_padding = result_len
assert len_after_padding >= result_len
sp1 = fft(a, (len_after_padding,))
sp2 = fft(b, (len_after_padding,))
ret = ifft(sp1 * sp2, (len_after_padding,))
return ret[:result_len]
def next_power_of_2(x):
return 1 << (x - 1).bit_length()
rng = np.random.default_rng(seed=10000)
a, b = rng.random((800,)), rng.random((800,))
result_direct = np.convolve(a, b) # direct
result_fft_no_padding = convolve_by_fft(a, b) # fft no padding
len1 = next_power_of_2(a.size + b.size - 1)
result_fft_padding_len1 = convolve_by_fft(a, b, len1) # fft padding to next power of 2
len2 = 2 * len1
result_fft_padding_len2 = convolve_by_fft(
a, b, len2
) # fft padding to next next power of 2
# verify the error of last element
print(f"last element golden: {a[-1] * b[-1]}")
print(
f"last element direct: {result_direct[-1]}, error = {(result_direct[-1] - a[-1] * b[-1]) / (a[-1] * b[-1])}"
)
print(
f"last element fft_no_padding: {result_fft_no_padding[-1]}, error = {(result_fft_no_padding[-1] - a[-1] * b[-1]) / (a[-1] * b[-1])}"
)
print(
f"last element fft_padding_len1: {result_fft_padding_len1[-1]}, error = {(result_fft_padding_len1[-1] - a[-1] * b[-1]) / (a[-1] * b[-1])}"
)
print(
f"last element fft_padding_len2: {result_fft_padding_len2[-1]}, error = {(result_fft_padding_len2[-1] - a[-1] * b[-1]) / (a[-1] * b[-1])}"
)
# plot the error compare to the direct method
plt.plot(
(result_fft_no_padding - result_direct) / result_direct,
label="error(fft_no_padding)",
)
plt.plot(
(result_fft_padding_len1 - result_direct) / result_direct,
label="error(fft_padding_len1)",
)
plt.plot(
(result_fft_padding_len2 - result_direct) / result_direct,
label="error(fft_padding_len2)",
)
plt.legend()
plt.show()
```
