# Numerically stable fast convolution algorithm?

Suppose you have two vectors of real numbers $\langle a_1, \dots, a_n \rangle, \langle b_1, \dots, b_n \rangle$, with $a_i, b_i \geq 0$, and wish to compute the convolution $$c_i = \sum_{j \leq i} a_j b_{i-j}$$

There is an obvious algorithm to compute this in time $O(n^2)$. This obvious algorithm basically involves a sum of many positive summands, so there is no numerical cancellation. Hence the obvious algorithm is numerically stable.

The other obvious algorithm is to use a Fast Fourier transform to compute $\hat a, \hat b$, obtain $\hat c = \hat a \hat b$, and then use an inverse transform to obtain $c$. This algorithm is much faster, like $O(n \log n)$ time. Unfortunately, the Fourier transforms require adding many terms with different complex phases, so there is catastrophic cancellation. So this algorithm is numerically unstable.

Is there any algorithm which is both fast (near-linear time) and numerically stable?

• Do you have a formal definition of numerically stable. $\:$ Your description would suggest $\hspace{0.8 in}$ "whenever the entries are all non-negative, the algorithm won't add reals with opposite signs". – user6973 Aug 9 '12 at 18:45