PaperlayerView exact revision on arXiv

Source-linked AI summary

A Fast Butterfly Algorithm for the Computation of Fourier Integral Operators

Emmanuel Candes, Laurent Demanet, Lexing Ying

arXiv:0809.0719v1math.NAcs.NA

TL;DR

Direct evaluation of discrete FIOs is costly, and their oscillatory kernels obstruct standard multiscale techniques. The paper uses low-rank structure on paired time-frequency regions with a butterfly algorithm, achieving O(N^2 log N) complexity and substantial speedups on large grids.

  • Problem

    Direct evaluation costs O(N^4), while kernel oscillation is the main obstacle to fast algorithms using standard multiscale techniques.

    Source: §1 · p. 2

  • Method

    The method exploits constant-rank interactions between complementary-level spatial and frequency boxes, constructing approximations with oscillation-prefactored interpolation.

    Source: §2 · p. 6, §3.2 · p. 15

  • Results

    At 4,096 × 4,096 grids, the algorithm achieves nearly two digits of accuracy with a speedup factor exceeding 750 over direct evaluation.

    Source: §5 · p. 21

  • Takeaways & Limitations

    For large grids of practical interest, the algorithm provides substantial speedup over direct evaluation.

    Source: §5 · p. 21

    Tensor-product Chebyshev approximations often have substantially higher separation rank than kernel-optimized interpolative decompositions.

    Source: §4 · p. 18

Abstract

from arXiv · show

This paper is concerned with the fast computation of Fourier integral operators of the general form $\int_{\R^d} e^{2\pi\i \Phi(x,k)} f(k) d k$, where $k$ is a frequency variable, $\Phi(x,k)$ is a phase function obeying a standard homogeneity condition, and $f$ is a given input. This is of interest for such fundamental computations are connected with the problem of finding numerical solutions to wave equations, and also frequently arise in many applications including reflection seismology, curvilinear tomography and others. In two dimensions, when the input and output are sampled on $N \times N$ Cartesian grids, a direct evaluation requires $O(N^4)$ operations, which is often times prohibitively expensive. This paper introduces a novel algorithm running in $O(N^2 \log N)$ time, i. e. with near-optimal computational complexity, and whose overall structure follows that of the butterfly algorithm [Michielssen and Boag, IEEE Trans Antennas Propagat 44 (1996), 1086-1093]. Underlying this algorithm is a mathematical insight concerning the restriction of the kernel $e^{2\pi\i \Phi(x,k)}$ to subsets of the time and frequency domains. Whenever these subsets obey a simple geometric condition, the restricted kernel has approximately low-rank; we propose constructing such low-rank approximations using a special interpolation scheme, which prefactors the oscillatory component, interpolates the remaining nonoscillatory part and, lastly, remodulates the outcome. A byproduct of this scheme is that the whole algorithm is highly efficient in terms of memory requirement. Numerical results demonstrate the performance and illustrate the empirical properties of this algorithm.

1 Introduction

The introduction frames discrete FIO evaluation as expensive because oscillatory kernels obstruct standard multiscale methods, then proposes a butterfly-structured O(N^2 log N) approach based on locally low-rank kernel restrictions.

1 Introduction

  • On N × N grids, direct discrete FIO evaluation costs O(N^4), motivating an O(N^2 log N) algorithm for smooth phases homogeneous of degree 1 in frequency.The input and output are sampled on Cartesian grids.Source: §1 · pp. 1–2

1.1 General strategy

  • When w(A)w(B) ≤ 1/N, the residual phase is O(1/N), making the restricted kernel nonoscillatory and permitting fixed-accuracy low-rank separation.The approximation decouples spatial and frequency variables over the paired boxes.Source: §1.1 · pp. 2–3
  • The proposed interpolation prefactors the oscillation, interpolates the remaining nonoscillatory function, and remodulates the result.Tensor-product Chebyshev interpolation is applied in x when w(A) ≤ 1/N.Source: §1.1 · p. 3
  • The butterfly recursion builds compact expansions by descending the spatial quadtree and ascending the frequency quadtree, reusing coefficients from the previous level.It starts with frequency leaves paired with the spatial root and ends with output values across spatial leaves.Source: §1.1 · p. 4

1.2 Applications

  • FIOs include wave-equation solution operators, while reflection-seismology imaging is an FIO whose O(N^2 log N) evaluation contrasts with O(N^3) Kirchhoff migration.For inhomogeneous wave speeds, the phase is a Hamilton-Jacobi equation solution under mild smoothness assumptions and for sufficiently small times.Source: §1.2 · pp. 4–5

2 The Butterfly Algorithm

The butterfly algorithm exploits approximately constant-rank interactions between paired spatial and frequency boxes, then propagates compact potential expansions through two dyadic trees. Its recursive structure gives O(rε^2 M log M) work, while the paper distinguishes interpolation from earlier collocation or interpolative-decomposition approaches.

2 The Butterfly Algorithm

  • For box pairs satisfying ℓ(A) + ℓ(B) = L, the restricted kernel has approximately constant rank, with ε-separation rank independent of M.The expansions provide compact representations of the potential restricted to A.Source: §2 · pp. 6–8
  • The algorithm initializes at the spatial root and frequency leaves, then traverses the spatial tree top-down and the frequency tree bottom-up, maintaining balanced levels.Paired box widths have product 1/N in the illustrated FIO setting.Source: §2 · pp. 7–8, Figure 2 · p. 9
  • O(rε^2 M log M) operations suffice when each pair's expansion coefficients are computed by a linear transform of size O(rε × rε).There are O(M) box pairs per level and O(log M) levels.Source: §2 · p. 8
  • Interpolative decomposition constructs factors from selected kernel columns, but its quadratic precomputation time and O(rε^2 N log N) storage can hinder large or infrequently repeated applications.The storage burden can become a bottleneck even at moderate sizes when rε^2 is nonnegligible.Source: §2 · pp. 9–10
  • For sparse Fourier transforms supported on curves, prior work uses adaptive quadtrees and an essentially box-independent transform to avoid quadratic-time precomputation and explicit transform storage.That approach concerns the kernel e^(2πi x·p/N).Source: §2 · pp. 10–11
  • The paper's low-rank approximation uses tensor-product Chebyshev grids with interpolation rather than collocation.The passage identifies this as a distinction from the earlier curve-supported sparse Fourier approach.Source: §2 · p. 11

3 Low-rank Approximations

This section studies FIO kernels on uniformly partitioned domains and develops low-rank approximations for spatial-frequency box pairs. Under analyticity and a box-width condition, it establishes a rank bound and uses Chebyshev interpolation to construct the approximations.

3 Low-rank Approximations

  • The section considers points in [0, 1]^2 and uniformly partitions both domains into quadtrees whose finest boxes have sidelength 1/N.The phase function Ψ(x, p) is assumed real-analytic in the joint variables.Source: §3 · p. 11, §3.1 · p. 11

3.1 The low-rank property

  • For boxes with w(A)w(B) = 1/N, the restricted residual kernel admits an ε-accurate separated approximation with rank bounded by a constant multiple of log^4(1/ε).The factors can be chosen as monomials centered at the box centers, with degrees bounded by a constant times log2(1/ε). The full kernel has a separated approximation of the same rank.Source: §3.1 · pp. 11–13

3.2 Interpolation gives good low-rank approximations

  • Tensor-product Chebyshev Lagrange interpolation of the residual kernel in either p or x yields low-rank approximations using qε ≲ log2(1/ε) grid order.The paper forms a special interpolant by prefactoring the oscillation, interpolating the residual, and remodulating the outcome.Source: §3.2 · pp. 14–15

4 Algorithm Description

The algorithm combines low-rank kernel expansions with a butterfly traversal of spatial and frequency quadtrees. It uses tensor-product Chebyshev interpolation and switches representations midway through the recursion.

4 Algorithm Description

  • The butterfly traversal builds spatial and frequency quadtrees, then descends in the spatial tree while ascending in the frequency tree to update box-pair expansions.The quadtrees have leaf boxes of size 1/N × 1/N; initialization constructs coefficients for each frequency leaf against the spatial root.Source: §4 · p. 16
  • At the middle level, the recursion switches from frequency-side equivalent sources to spatial-side potential samples, then continues to the spatial leaves.The frequency-side interpolant is valid through level L/2, while the spatial-side interpolant is valid from level L/2 onward.Source: §4 · pp. 16–17
  • For fixed accuracy ε, the overall cost is O(N^2 log N), with a worst-case polylogarithmic dependence on 1/ε.The stated estimate is O(log^6(1/ε) N^2 log N + log^8(1/ε) N^2); empirically, the dependence grows more moderately, like log(1/ε).Source: §4 · p. 18
  • Fixed Chebyshev grids give explicit expansion formulas without precomputation, and storing only coefficients from two consecutive levels reduces memory to O(rε N^2).The interpolative-decomposition approach instead stores transformations for box pairs, with storage on the order of rε^2 N^2 log N.Source: §4 · p. 18
  • Chebyshev interpolation can produce substantially higher separation rank than optimized interpolative decomposition, although tensor-product structure keeps the added computation from being significant.The Chebyshev approximation exploits residual-kernel smoothness rather than the kernel's finer structure.Source: §4 · p. 18
  • Unlike the cited Fourier-kernel collocation approach, this method interpolates the residual because FIO kernels lack an (A, B)-independent tensor-product decomposition.It retains tensor-product Chebyshev grids while adapting the interpolation variable to the smaller box.Source: §4 · p. 19
  • For variable amplitudes, the method assumes a separated approximation with a number of terms independent of N, runs the butterfly algorithm per term, and combines the results.A vectorized implementation processes the terms simultaneously to avoid redundant kernel evaluations.Source: §4 · p. 19

5 Numerical Results

The numerical tests report near-asymptotic runtime scaling and substantial speedups over direct evaluation, including for a 3D example. They also examine accuracy as interpolation-grid size and problem size vary.

5 Numerical Results

  • The implementation is in C++ and the tests were run on a desktop computer with a 2.8GHz CPU.In the reported runs, the recursion starts at level log2 N − 3 and stops at level 3.Source: §5 · pp. 19–20
  • For the variable-amplitude example, the 1e-7 approximation uses three separated terms, and Table 2 reports results for different N and q combinations.The example uses amplitudes and phase functions specified in (5.3).Source: §5 · pp. 20–21, Table 2 · p. 22
  • At 4,096 × 4,096, nearly two digits of accuracy came with a speedup exceeding 750 over direct evaluation.The empirical running time closely follows O(N^2 log N); doubling N increases runtime by roughly a factor of four.Source: Table 1 · p. 21, §5 · p. 21
  • Increasing the Chebyshev grid size q by a factor of two improves estimated accuracy εa by an average factor of 30, while accuracy decreases only slightly as N increases.The authors attribute the latter behavior to the near-optimal Lebesgue constant of Chebyshev interpolation.Source: §5 · p. 21
  • For N≥512, the butterfly approach clearly outperforms the earlier conical-partitioning method, while that method is faster for N≤256.The two approaches take roughly the same time at N=256 and 512.Source: §5 · p. 22
  • The 3D example shows significant speedup even for moderate N.The tests use a six-dimensional phase function and report results in Table 3.Source: §5 · p. 22, Table 3 · p. 22

6 Conclusions and Discussions

The paper combines geometric low-rank structure with the butterfly algorithm to compute discrete FIOs in O(N^2 log N) time. The conclusion points to related fast transforms and to curvelet or wave-atom approaches as directions for further algorithmic work.

6 Conclusions and Discussions

  • Geometric low-rank structure and the butterfly algorithm yield an accurate discrete-FIO algorithm with O(N^2 log N) complexity.The kernel has approximately low rank when restricted to special time- and frequency-domain subdomains.Source: §6 · p. 23
  • Operationalizing the optimal sparsity of curvelet and wave-atom FIO representations into efficient algorithms remains an open research direction.The paper attributes the curvelet result to Candès and Demanet, notes the same optimality for wave atoms, and connects it to Stein et al.’s second dyadic decomposition.Source: §6 · p. 23
  • The geometric low-rank and butterfly combination has also been used for fast special transforms and Fourier transforms with sparse data or summation constraints.The paper cites O’Neil and Rokhlin for special transforms and Ying [39] [40] for the two Fourier-transform settings.Source: §6 · p. 23
Loading 0809.0719v1…