Convolution and the Fast Fourier Transform Tarik Guelzim

Convolution and the Fast Fourier Transform Tarik Guelzim

Vectors • Given two vectors • a=(a 0, a 1, . . . , an-1) • b=(b 0, b 1, . . . , bn-1) • we can compute the sum a+b • or the inner product, producing a. b • we can also produce the convolution a*b

Convolution • The convolution of two vectors of length n is a vector with 2 n-1 coordinates where: ai b i (i, j): i+j=k i, j<n In other words, a*b = (a 0 b 0, a 0 b 1 + a 1 b 0, a 0 b 2+a 1 b 1+a 2 b 0, . . . , an-2 bn-1+ an-2 bn-2, an-1 bn-1).

Simply speaking. . . • This definition is simply: a 0 b 0 a 1 b 0 a 2 b 0. . . an-1 b 0 a 0 b 1 a 1 b 1 a 2 b 1. . . an-1 b 1 . . . a 0 bn-1 a 1 bn-1 a 2 bn-1. . . an-1 bn-1 we compute the coordinates in the convolution vector by summing along the diagonals.

simply it is. . . • This definition is simply: a 0 b 0 a 1 b 0 a 2 b 0. . . an-1 b 0 a 0 b 1 a 1 b 1 a 2 b 1. . . an-1 b 1 a 0 b 2. . . a 0 bn-1 a 1 bn-1 a 2 bn-1. . . an-1 bn-1 a*b = (a 0 b 0, a 0 b 1 + a 1 b 0, a 0 b 2+a 1 b 1+a 2 b 0, . . . , an-2 bn-1+ an-2 bn-2, an-1 bn-1).

Polynomial Multiplication A(x) = a 0 + a 1 x + a 2 x 2 +. . . + am-1 xm-1 B(x) = b 0 + b 1 x + b 2 x 2 +. . . + bm-1 xm-1 we can represent them by vector coefficients where: a = (a 0, a 1, . . , am-1) and b = (b 0, b 1, . . , bm-1) so now given A(x) and B(x) consider C(x) such that C(x) = A(x)B(x) • the coefficients of C(x) are c = (c 0, c 1, . . , c 2 n-2) • we can compute the coef. of xk • • • Ck = aibj (i, j): i+j=k

Example • A(x) = 1 + 2 x + 3 x 2 and B(x) = 2 + 3 x +4 x 2 • A = (1, 2, 3) and B = (2, 3, 4) 1*2 2*2 3*2 1*3 2*3 3*3 1*4 2*4 3*4 = 2 4 6 3 6 9 A(x) * B(x) = 2+7+16 x 2+17 x 3+12 x 4 C = (2, 7, 16, 17, 12) 4 8 12

Computing the convolution • Computing the convolution with this method takes O(n 2) because for every pair (i, j) we calculate the product aibi • Can we do better that this quadratic running time?

Algorithm design • Again suppose we are given two vectors: • a = (a 0, a 1, . . , am-1) and b = (b 0, b 1, . . , bm-1) • what we want is to compute the product of C(x) = A(x)B(x) in O(n log n). • once done, we can read off the desired answer directly from the coef. of C(x).

• Rather that multiplying A and B symbolically we treat them as functions of x. • step 1: we choose 2 n values (C has 2 n-1 coordinates) x 1, x 2, . . . , x 2 n and evaluate A(xj) and b(xj) for each j = 1, 2, . . . , 2 n • step 2: we compute C(Xj) simply by A(xj)*b(xj) • step 3: we can recover C from its values on x 1, x 2, . . , x 2 n by using polynomial interpolation.

Polynomial interpolation • any polynomial of degree d can be reconstructed from its values on any set of d+1 or more points. • we can observe that since A and B each have a degree at most n-1, the product C has degree at most 2 n-2 • [a(x) = 0+0 x+x 2 and b(x) = 0+0 x+x 2 • a*b = x 4 ]

The problem is. . . • step 1: we choose 2 n values (C has 2 n-1 coordinates) x 1, x 2, . . . , x 2 n and evaluate A(xj) and b(xj) for each j = 1, 2, . . . , 2 n • step 2: we compute C(Xj) simply by A(xj)*b(xj) O(n) • step 3: we can recover C from its values on x 1, x 2, . . , x 2 n by using polynomial interpolation. since we are trying to perform 2 n evaluations we fall back into O(n 2) again

The Trick • the idea is to find a set of 2 n values x 1, x 2, . . . , x 2 n that are intimately related in some way such that the work in evaluating A and B on all of them can be shared across different evaluations.

Complex Roots of Unity • recall that a complex number x+yi can also be written in polar coordinate as: reφi • where eπi = -1 and e 2πi = 1 • for a positive integer k, the polynomial equation xk = 1 has k distinct complex roots and they are • ωj, k= e 2πji/k ( for j = 0, 1, 2, . . . , k-1) • (e 2πji/k )k = e 2πji = (e 2πi)j = 1 • each of these numbers is distinct, so they are all of the roots.

• The representation of polynomial P by its values on the (d +1)st roots of unity is referred to as the discrete Fourier transform (DFT).

Polynomial Evaluation • our goal is to design an algorithm that evaluates A on each of the (2 n)th roots of unity recursively. • we need to satisfy the recurrence relation 2 T(n/2)+O(n).

Divide and Conquer we need to define two polynomials Aeven and Aodd Aeven(x)=a 0+a 2 x+a 4 x 2+. . . +an-2 x(n-2)/2 Aodd(x)=a 1+a 3 x+a 5 x 2+. . . +an-1 x(n-1)/2 A(x) = Aeven(x 2)+x. Aodd(x 2) Aeven and Aodd both have half the input and only n roots of unity to evaluate so the total time is • 2 T(n/2) • • •

• To satisfy the recurrence relation we need to evaluate A(x) = Aeven(x 2)+x. Aodd(x 2) in O(n). • given ωj, 2 n= e 2πji/2 n the quantity ω2 j, 2 n= (e 2πji/2 n )2 = e 2πji/n and this is an nth root of unity so when we try to compute A(x) = Aeven(x 2)+ x. Aodd(x 2) the right side has already been performed in a recursive step. • thus, A(x) for 2 n can be done in O(n). with this we satisfied the O(n log n) bound.

Reconstructing C • we can do the same to evaluate B on the (2 n)th roots of unity and to compute C(ωj, n) = • • A(ωj, 2 n)B(ωj, 2 n) consider a polynomial : c(x) = Σ 2 n-1 c xs s s=0 we want to reconstruct it from its values C(ωj, 2 n) we need to define D(x) = 2 n-1 Σ d xs s=0 Where ds = C(ωj, 2 n) s

Theorem 2 n-1 s=0 For any polynomial 2 n-1 • D(ωj, 2 n) = Σ s=0 we have that – c(x) =Σ C(ωs, 2 n) xs Cs c sxs = 1/2 n D(ω2 n-s, 2 n) After evaluating the polynomial D(x) at the (2 n)th roots of unity give us the coefficient of the polynomial C(x) in reverse order All the evaluation of D can be done in O (n log n) using the same divide and conquer approach we showed previously.

What we proved is. . . • Using the FFT to determine the product polynomial C(x), we can compute the convolution of the original vectors a and b in o(n log n) instead of O(n 2)

Applications of Convolution • Combining histograms is one application of Convolution. • For example if we have two histograms of men and women and their salaries and we want to get another histogram from the previous to of all pairs (man, woman) that have a combined income of k we in a sense add all the terms i+j=k. • Most of the applications of convolution are in the signal processing domain where it is used to convert data from the time domain to the frequency domain.