1 FFT on Multipoles (FFTM) FFT accelerated SLFMM Lecture 23 FFTM (Fast Fourier Transform on Multipoles) • FFTM is, in fact, SLFMM speeded up with the FFT, which brings its complexity to O(NlogN). E. T. Ong, K. M. Lim, K. H. Lee and H. P. Lee: A fast algorithm for three-dimensional potential fields calculation: fast Fourier transform on multipoles Journal of Computational Physics, Volume 192, Issue 1, 20 November 2003, Pages 244-2612 Idea of a Single Level FMM Sources Sources Evaluation Points Evaluation Points Standard algorithm SLFMM N M N M Total number of operations: O(NM) Total number of operations: O(N+M+KL) K groups L groups Spatial Domains E E 1 E 3 E 2 1 E 3 E 2 n n n I1(n) = n I2(n) = {Neighbors(n)} n I3(n) = {All boxes}\ I2(n) Boxes with these numbers belong to these spatial domains F1(n)(y) F2(n)(y) F3(n)(y) Potentials due to sources in these spatial domains3 Definition of potentials SLFMM Algorithm Step 1. Generate S-expansion coefficients for each box For n  NonEmptySource Get xc(n) , the center of the box; loop over all non-empty source boxes For xi  E1(n) Get B (xi, xc(n)) , the S-expansion coefficients near the center of the box; C(n) = C(n) + uiB (xi, xc(n)); C(n) = 0; End; End; loop over all sources in the box Implementation can be different! All we need is to get C(n).4 SLFMM Algorithm Step 2. (S|R)-translate expansion coefficients For n  NonEmptyEvaluation Get xc(n) , the center of the box; loop over all non-empty evaluation boxes For m  I3(n) Get xc(m) , the center of the box; D(n) = D(n) + (S|R)(xc(n)- xc(m)) C(m) ; D(n) = 0; End; End; loop over all non-empty source boxes outside the neighborhood of the n-th box Implementation can be different! All we need is to get D(n). S|R-translation5 SLFMM Algorithm Step 3. Final Summation For n  NonEmptyEvaluation Get xc(n) , the center of the box; For xi  E2(n) vj = D(n) R(yj - xc(n)) ; End; End; loop over all evaluation points in the box Implementation can be different! All we need is to get vj loop over all boxes containing evaluation points For yj  E1(n) vj = vj +F(yj , xi); End; loop over all sources in the neighborhood of the n-th box FFTM6 Convolution and FFT This can be extended to d-dimensions • So, for any kernel, which can be approximated by Fourier series O(NlogN) algorithm exist to compute convolution.7 Comparison of the FFTM and FMM • FFTM has better error bounds; • FFTM can be faster for relatively uniform distributions; • FFTM does not require hierarchical data structures; • But… for non-uniform distributions FMM, perhaps, is faster…8 From Ong, et al. (2003). Comparison with a non-efficient FMM Non Uniform FFT Lecture 239 NUFFT • Uniform and Nonuniform Discrete Fourier Transforms – DFT, FFT, NUFFT, IDFT, IFFT, INUFFT – Statement of the problems • INUFFT – Data Structure – Kernel Decomposition – “Middleman” part of the algorithm – FMM part of the algorithm – Error bounds – Complexity and Optimization • Some Numerical Results – Error Analysis – Performance Band-limited functions bandwidth 2p-periodic function (complex) Fourier coefficients (complex) real some people like functions10 Inverse and Forward Discrete Fourier Transform 2p f(x) x x1 x0 x2 xk xN-1 f0 f1 f2 fk fN-1 xk = 2pk/N, fk = f(xk), k = 0, … , N-1 IDFT: {cn}→{fk} Straightforward: O(N2) IFFT/FFT: O(NlogN) DFT: {fk} → {cn} (Cooley & Tukey, 1965) Fast Fourier Transform • Divide and conquer type O(N logN) algorithm; • Exact in exact arithmetic; • Revolutionized computations: – Signal processing; – Scientific computing; – Large data sets handling/compression; – Many other applications… • Requires equispaced (uniform) data11 Fast Fourier Transform • FFT Requires equispaced (uniform) data. • In many applications data is not uniform: – Jitter; – Non-uniform grids for computations; – Events that occur randomly – Etc. • An algorithm of complexity O(NlogN) to perform the forward and inverse Fourier Transform for such cases is highly desired. Publications on the NUFFT • Dutt & Rokhlin (1993): – Interpolation using Gaussians; – Forward problem O(NlogN); – Inverse problem: iteratively O(NiterNlogN). • Dutt & Rokhlin (1995): – FMM for kernel cot(t/2) and ln(sin(t/2)); – Forward problem O(NlogN); – Inverse problem: O(NlogN). • Beylkin (1995): – Interpolation using multiresolution analysis (MRA). • Liu & Nguyen (1997,1998): – Least square approximations; – Forward problem O(NlogN); – Inverse problem: iteratively, O(NiterNlogN). OUR WORK IS CLOSE TO THIS ONE, WHILE FACTORIZATION AND DATA STRUCTURE ARE DIFFERENT12 Statement of the Problem Inverse Nonuniform Discrete Fourier Transform gj = f(yj), j = 0, … , N-1 INUDFT: {cn}→{gj} Straightforward: O(N2) Looking for: O(NlogN) (INUFFT) 0 2p f(x) x y1 y0 y2 yj yN-1 g0 g1 g2 gj gN-113 INUDFT: Reduction to Interpolation Problem INUFFT: {cn} → {gj} {fk} IFFT FFIA = “Fast Fourier Interpolation Algorithm” O(NlogN) O(NlogN) equispaced samples FFIA Fast Oscillations Modulation14 Forward Nonuniform Discrete Fourier Transform gj = f(yj), j = 0, … , N-1 NUDFT: {gj} → {cn} Straightforward matrix inversion: O(N3) Dutt & Rokhlin (1995): Straightforward: O(N2) Looking for: O(NlogN) (NUFFT) 0 2p f(x) x y1 y0 y2 yj yN-1 g0 g1 g2 gj gN-1 NUDFT: Reduction to Inverse Interpolation Problem NUFFT: {gj} → {cn} {fk} FFT IFFIA O(NlogN) O(NlogN) equispaced samples15 Inverse Interpolation Problem Computation of Products16 FFIA and IFFIA can be reduced to computation of sums: Here xk and yj can be arbitrary (non-uniform) We can perform fast summation using the Fast Multipole Methods Solution of the Problem17 Kernel Decomposition P0P1P2P3xkyjP2P1P3P0P0P3P2t = yj-xH(t) = cot(t/2)HtP0P1P2P3xkyjP0P1P2P3xkyjP2P1P3P0P0P3P2t = yj-xH(t) = cot(t/2)HtW1 W2 Kernel Decomposition (2) B2m – Bernoulli numbers18 Bernoulli Numbers Generating Function: Relation to the Riemann Zeta Function: Bounds: Rational. First 60 numbers can be found in Abramowitz and Stegun (1964): Factorization of the regular part of the Kernel19 Error bounds for the regular part Middleman for the regular part Solution Since we have polynomial factorization of the regular part for the entire domain, we can apply the “Middleman” method (the fastest possible) to compute this part. Complexity: O(qN): For the forward interpolation problem (FFIA): q~log4(1/e) (does not depend on