Multiscale Butterfly Algorithm

• The multiscale butterfly algorithm is a hierarchical method that exploits complementary low-rank structures to approximate high-dimensional oscillatory kernels with near-linear complexity.
• It uses corona-based decomposition to partition the frequency domain, effectively isolating singularities and enabling efficient recursive factorization.
• Recent extensions apply tensor decompositions to achieve linear scaling, reducing memory and computation in large-scale 3D and high-frequency problems.

The multiscale butterfly algorithm is a class of hierarchical, data-sparse algorithms designed for efficient approximation and application of high-dimensional oscillatory integral operators and kernel matrices, especially those arising in the discretization of Fourier integral operators (FIOs), Green’s functions, and related transforms. The multiscale butterfly approach achieves near-linear or subquadratic complexity for problems that would otherwise require $\mathcal{O}(N^{2d})$ operations, where $N$ is the number of degrees of freedom per dimension and $d$ the dimensionality. The algorithm exploits special low-rank structure, notably the complementary low-rank property, present in off-diagonal blocks away from kernel singularities, and uses a recursive multiscale decomposition of the frequency or spatial domain to localize and factorize interactions efficiently (Li et al., 2015, Li et al., 2014, Kielstra et al., 2024).

1. Mathematical Foundation: Complementary Low-Rank Structure and Corona Decomposition

The foundational principle of the multiscale butterfly algorithm is the complementary low-rank property. Let $K = (K(x,\xi))_{x \in X, \,\xi \in \Omega}$ be an $N \times N$ kernel matrix with $X$ and $\Omega$ discrete sets in $\mathbb{R}^d$. The complementary low-rank property states that for two quadtrees $T_X$ and $T_\Omega$ of depth $N$0, each submatrix $N$1, indexed by blocks $N$2 (at level $N$3 in $N$4) and $N$5 (at level $N$6 in $N$7), has numerical rank $N$8 independent of $N$9 (with polynomial dependence on $d$0 for accuracy $d$1):

$d$2

This property is leveraged only away from singularities (e.g., $d$3 in FIOs). For domains with singularities, the frequency (or analogous) domain is decomposed into a central low-frequency region and a sequence of “coronas” (dyadic annuli or axis-aligned cubes) $d$4, achieving low-rank structure within each corona (Li et al., 2015, Li et al., 2014).

2. Algorithmic Structure and Hierarchical Factorization

The algorithm builds two hierarchical trees (typically quadtrees in 2D, octrees in 3D) over the spatial and frequency domains. At the “middle” scale (level $d$5, the kernel matrix is block-diagonalized:

$d$6

where $d$7 and $d$8 are block-diagonal matrices with blocks of size $d$9, and $K = (K(x,\xi))_{x \in X, \,\xi \in \Omega}$0 is block-sparse. Recursive compression is then performed up and down the tree structures:

$K = (K(x,\xi))_{x \in X, \,\xi \in \Omega}$1

Iteratively, this forms an overall factorization:

$K = (K(x,\xi))_{x \in X, \,\xi \in \Omega}$2

Each factor is data-sparse with $K = (K(x,\xi))_{x \in X, \,\xi \in \Omega}$3 non-zeros, and the total number of factors is $K = (K(x,\xi))_{x \in X, \,\xi \in \Omega}$4 (Li et al., 2015).

In the multiscale butterfly algorithm, this process is repeated independently for each corona. The central small region is treated directly, and each corona’s contribution is handled using the above recursive butterfly factorization (Li et al., 2014).

3. Extension to Fourier Integral Operators and Corona-Based Partitioning

For Fourier integral operators (FIOs), where the kernel $K = (K(x,\xi))_{x \in X, \,\xi \in \Omega}$5 exhibits singular behavior near $K = (K(x,\xi))_{x \in X, \,\xi \in \Omega}$6, the frequency domain is decomposed as

$K = (K(x,\xi))_{x \in X, \,\xi \in \Omega}$7

Here, $K = (K(x,\xi))_{x \in X, \,\xi \in \Omega}$8 is the central cube, and $K = (K(x,\xi))_{x \in X, \,\xi \in \Omega}$9 are successive dyadic coronas. On each $N \times N$0, the phase is sufficiently regular to admit a rank-$N \times N$1 separated form over appropriately chosen spatial/frequency boxes (see Theorem 3.1 in (Li et al., 2014)). The full operator action is then written as:

$N \times N$2

with $N \times N$3 restriction operators. Each $N \times N$4 is then factorized independently using the butterfly scheme (Li et al., 2014, Li et al., 2015).

This cartesian corona-based multiscale decomposition allows the butterfly algorithm to efficiently bypass the singular origin and to maintain low application and factorization costs.

4. Construction Methods and Implementation Considerations

Two principal construction methods are employed for the kernel factorizations at the middle levels:

• Black-box matvec (fast apply): Requires access to routines that compute $N \times N$5 and $N \times N$6 at $N \times N$7 cost. Block SVDs at the middle level are built via randomized algorithms using Gaussian test matrices, with costs scaling as $N \times N$8 (Li et al., 2015).
• Entrywise sampling (random sampling): Assumes access to individual $N \times N$9 in $X$0. Randomized sampling and pivoted QR/ID are used to compute the basis matrices, with costs $X$1 (Li et al., 2015).

Recursive upward/downward compression of the hierarchical block factors is applied after the middle-level SVDs. Chebyshev interpolation in both the spatial and frequency directions is typically employed for numerical stability and error control (Börm et al., 2017, Li et al., 2014).

5. Complexity, Accuracy, and Numerical Performance

The multiscale butterfly algorithm achieves, for $X$2 degrees of freedom:

Step	Complexity (Sampling)	Complexity (Fast Matvec)
Factorization	$X$3	$X$4
Application	$X$5	same
Memory	$X$6	same (can be reduced to $X$7 with interleaving)

Accuracy depends on the prescribed butterfly rank $X$8, which scales polynomially in $X$9. For example, for $\Omega$0, $\Omega$1, direct and factorized applications are accelerated by factors of $\Omega$2 at errors $\Omega$3 (Li et al., 2015, Li et al., 2014).

The corona-based multiscale version exhibits substantially reduced prefactors and eliminates the need for polar (coordinate) transformations. Empirical results indicate speedups of $\Omega$4 to $\Omega$5 over polar-butterfly and direct FIO evaluation in both 2D and 3D (Li et al., 2014).

6. High-Dimensional and Tensor Extensions

Recent algorithmic developments generalize the multiscale butterfly structure to tensors, yielding linear $\Omega$6 scaling for $\Omega$7-mode tensors in $\Omega$8 dimensions (Kielstra et al., 2024). The tensor butterfly algorithm partitions modes along each dimension and applies mode-wise interpolative decomposition (Tucker-ID), maintaining small ranks and achieving superior constants relative to quantized tensor train (QTT) and matrix-based butterfly algorithms. In large-scale 3D Helmholtz problems (up to $\Omega$9 per side), this approach enables $\mathbb{R}^d$0 speedup and $\mathbb{R}^d$1 memory reduction versus matrix butterfly (Kielstra et al., 2024).

This tensor extension is the first to achieve linear (in $\mathbb{R}^d$2) factorization and contraction costs for high-dimensional oscillatory integral operators, outperforming existing FFT, matrix butterfly, and QTT-based schemes in large and high-frequency regimes (Kielstra et al., 2024).

7. Analytical Stability and Theoretical Guarantees

Refined stability analysis demonstrates that the use of Chebyshev polynomial interpolation and carefully controlled nested reinterpolation keeps operator norms and error amplification under control. The required interpolation degree $\mathbb{R}^d$3 scales as $\mathbb{R}^d$4 in practical regimes, which further reduces storage and computational complexity. For analytic kernels (e.g., those with phases and amplitudes analytic in neighborhoods of the real domain), the overall error exhibits exponential decay in $\mathbb{R}^d$5 at all hierarchical levels (Börm et al., 2017). The stability of the nested butterfly representation guarantees that complexity remains near-optimal, and that catastrophic error growth is avoided (Börm et al., 2017).

---

References:

(Li et al., 2014) "A Multiscale Butterfly Algorithm for Multidimensional Fourier Integral Operators" (Li et al., 2015) "Multidimensional Butterfly Factorization" (Börm et al., 2017) "An analysis of a butterfly algorithm" (Kielstra et al., 2024) "A Linear-complexity Tensor Butterfly Algorithm for Compressing High-dimensional Oscillatory Integral Operators"