>>> def fft_mul(C1, C2):Takes O(n logn) time. The idea is to right pad each polynomial with enough zeros so that the cyclic convolution becomes a noncyclic convolution. For more speed, pad
... d = len(C1)+len(C2)-1
... c1 = FFT.fft(list(C1) + [0] * (d-len(C1)))
... c2 = FFT.fft(list(C2) + [0] * (d-len(C2)))
... return list(FFT.inverse_fft(c1*c2)[:d].real)
...
>>> fft_mul([1,2,3],[4,5]) # (1+2*x+3*x**2)*(4+5*x)
[4.0, 13.0, 22.0, 15.0]
c1 and c2 so each has power-of-2 length. (Source Code).I posted the algorithm [1] on Wikipedia (alternate link: [2]).
7 comments:
Interesting post,
I wonder if this method could also be applied to finding sums of random variables.
For instance to get this picture(http://web.engr.oregonstate.edu/~bulatov/research/reports/vapnik/plot2.png), I had to get distributions of X1, X1+X2, X1+X2+X3+...+Xl Where X_i are IID variables.
You find this distribution using convolution. If we have P(X=k), then P(X1+X2)=k is
P(X=0)P(X=k)+P(X=1)P(X=k-1)+....
Which is really a convolution. So if a[k] gives P(X=k), then convolve(a,a)[k] will P(X1+X2)=k; def convolve(a,b)=[sum([a[j]*b[i-j] for j in range(0,i+1)]) for i in range(len(a))]
You can repeat the convolution l-1 times to get the sum of X1+...+Xl
The problem is that this gets really slow.
----
In statistics there's method for finding sums called "Method of Moments". They take two-sided Laplace transform of the probability density (which is called the Moment Generating Function http://en.wikipedia.org/wiki/Moment_generating_function),
and then convolutions of the densities become products.
So apparently Fourier transform isn't the only one where convolutions become products...
Is there some general theorem which says that convolutions turn into products for X,Y and Z transforms?
Curiously, in statistics books the inverse of the transform is done by "guess and check"...is the inverse transform really that hard?
Interestingly some proofs of the Central Limit Theorem operate the transformed space. IE, it's hard to see what'll happen if you take X1+X2+X3+....\infty if you viewing it as convolutions (ie, how exactly are you going to do the infinite dimensional integral?!?), but instead if you look at it as infinite product, it's easier. So it turns out that the infinite product converges to the same thing regardless of the thing you are multiplying, which could explain why people were so excited about Central Limit Theorem (http://www-stat.stanford.edu/~susan/courses/s116/node120.html)
You're right about the IID distribution convolutions. You can compute your distribution convolution with the fft_mul() code I gave.
>>> fft_mul([0.6, 0.3, 0.1], [0.8, 0.2])
[0.48, 0.36, 0.14, 0.02]
The reason it works is because the polynomial coefficients of c(x) = a(x)*b(x) are:
c[0] = a[0]*b[0]
c[1] = a[0]*b[1]+a[1]*b[0]
c[2] = a[0]*b[2]+a[1]*b[1]+a[2]*b[0]
...
That is, c[k] is the sum of all a[i] and b[j] such that i + j = k. Which is just P(X1+X2=k).
---
For integral transforms with an arbitrary kernel, you can just follow the proof of the convolution theorem to establish a sufficient condition on the kernel such that a star b = T(a)T(b).
I recall from Math 256 that analytic inverse Laplace transforms are usually impossible if you can't find them in some CRC type handbook. But apparently they can be done numerically.
You can also take polynomial square roots and multiply integers in this way. The NTT does it in exact integer arithmetic, or you can round the floating point results and rely on error bounds proved for IEEE-style floating point arithmetic.
FFT approach for finding convolution is indeed efficient but natural questions arises: what can be said about the error of this approximation. This is of course an approximation due to the representation of the numbers in computer arithmetic (in theory it should be exact). Since computation of the convolution using FFT is fast many you may use very long sequences but what about the round-off errors then? Do you know some sources in which one can find something about the error?
Hi
I tried to add from Numeric import *
but it still complains that it can't find FFT.
what a good method and I wondered the same Yaroslav, after reading this great post I asked myself if this method could also be applied to finding sums of random variables. If it could, it would be awesome!
I remember that once I tried to made a new calculator for times of the employes to Viagra Online labs and the source code it's very simple you don't need to made a lot of lines for that.
Post a Comment