7
$\begingroup$

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?

Improve this question
$\endgroup$
1
  • $\begingroup$ 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". $\endgroup$
    – user6973
    Aug 9, 2012 at 18:45

1 Answer 1

Reset to default
1
$\begingroup$

There is the number theoretic transform where all calculations are done using integers which should completely get round any concerns about numerical stability, at the possible expense of some speed (see e.g. https://ccrma.stanford.edu/~jos/st/Number_Theoretic_Transform.html). However this will only work for you if your numbers are rational and can be converted to integer representation easily.

There is also a whole range of asymptotically slower algorithms known as Toom-Cook multiplication http://en.wikipedia.org/wiki/Toom-Cook_multiplication. However, I am not sure it is true that multiplication via FFT calculations is necessarily numerically unstable. Do you have a reference for this claim?

Improve this answer
$\endgroup$
1
  • $\begingroup$ One interpretation of a DFT is in evaluating a polynomial at N points (these points being the roots of unity). If you interpret each signal being convolved as the coefficients to two different polynomial, and one of these polynomials has its roots concentrated over a small domain, then evaluating it far outside that domain would result in numbers much larger than when evaluating it near the roots, and thus the difference in magnitude may be impossible to capture with any conventional finite number system (e.g. floating point). $\endgroup$
    – Ponkadoodle
    Aug 25, 2016 at 20:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.