---
title: Multiscale Butterfly Algorithm
url: https://www.emergentmind.com/topics/multiscale-butterfly-algorithm
type: topic
---

# Multiscale Butterfly Algorithm

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 [1509.07925][1411.7418][2411.03029].

## 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 $L = O(\log N)$, each submatrix $K^{\ell}_{i,j}$, indexed by blocks $A^\ell_i$ (at level $\ell$ in $T_X$) and $B^{L-\ell}_j$ (at level $L-\ell$ in $T_\Omega$), has numerical rank $r$ independent of $N$ (with polynomial dependence on $\log(1/\varepsilon)$ for accuracy $\varepsilon$):
$$
K^{\ell}_{i,j} = (K(x,\xi))_{x \in A^\ell_i,\, \xi \in B^{L-\ell}_j} \approx U^\ell_{i,j} S^\ell_{i,j} (V^\ell_{j,i})^*
$$
This property is leveraged only away from singularities (e.g., $\xi=0$ 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) $\Omega_t$, achieving low-rank structure within each corona [1509.07925][1411.7418].

## 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 $h = L/2)$, the kernel matrix is block-diagonalized:
$$
K \approx U^h M^h (V^h)^*
$$
where $U^h$ and $V^h$ are block-diagonal matrices with blocks of size $\sim N/2^h$, and $M^h$ is block-sparse. Recursive compression is then performed up and down the tree structures:
$$
U^h \approx U^{h+1} G^h, \qquad (V^h)^* \approx (H^h)^* (V^{h+1})^*
$$
Iteratively, this forms an overall factorization:
$$
K \approx U^L G^{L-1} \cdots G^h M^h (H^h)^* \cdots (H^{L-1})^* (V^L)^*
$$
Each factor is data-sparse with $\mathcal{O}(N)$ non-zeros, and the total number of factors is $\mathcal{O}(\log N)$ [1509.07925].

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 [1411.7418].

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

For Fourier integral operators (FIOs), where the kernel $K(x,\xi) = a(x,\xi)e^{2\pi i \Phi(x,\xi)}$ exhibits singular behavior near $\xi=0$, the frequency domain is decomposed as
$$
\Omega = \Omega_C \cup \bigcup_{t=0}^{L-s} \Omega_t,\quad \Omega_t = \{\xi : n/2^{t+2} < \max{|\xi_1|,|\xi_2|} \le n/2^{t+1}\}
$$
Here, $\Omega_C$ is the central cube, and $\Omega_t$ are successive dyadic coronas. On each $\Omega_t$, the phase is sufficiently regular to admit a rank-$r$ separated form over appropriately chosen spatial/frequency boxes (see Theorem 3.1 in [1411.7418]). The full operator action is then written as:
$$
K = K_C R_C + \sum_{t=0}^{L-s} K_t R_t
$$
with $R_t$ restriction operators. Each $K_t$ is then factorized independently using the butterfly scheme [1411.7418][1509.07925].

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 $K g$ and $K^* g$ at $\mathcal{O}(N \log N)$ cost. Block SVDs at the middle level are built via randomized algorithms using Gaussian test matrices, with costs scaling as $\mathcal{O}(N^{3/2} \log N)$ [1509.07925].
- **Entrywise sampling (random sampling):** Assumes access to individual $K(x,\xi)$ in $\mathcal{O}(1)$. Randomized sampling and pivoted QR/ID are used to compute the basis matrices, with costs $\mathcal{O}(N^{3/2} r^2)$ [1509.07925].

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 [1703.01941][1411.7418].

## 5. Complexity, Accuracy, and Numerical Performance

The multiscale butterfly algorithm achieves, for $N = n^d$ degrees of freedom:

| Step              | Complexity (Sampling)   | Complexity (Fast Matvec) |
|-------------------|------------------------|--------------------------|
| Factorization     | $\mathcal{O}(N^{3/2})$ | $\mathcal{O}(N^{3/2} \log N)$ |
| Application       | $\mathcal{O}(N \log N)$| same                     |
| Memory            | $\mathcal{O}(N^{3/2})$ | same (can be reduced to $\mathcal{O}(N \log N)$ with interleaving) |

Accuracy depends on the prescribed butterfly rank $r$, which scales polynomially in $\log(1/\varepsilon)$. For example, for $n=256$, $r=14$, direct and factorized applications are accelerated by factors of $10^2$ at errors $\sim 2 \times 10^{-3}$ [1509.07925][1411.7418].

The corona-based multiscale version exhibits substantially reduced prefactors and eliminates the need for polar (coordinate) transformations. Empirical results indicate speedups of $6\times$ to $100\times$ over polar-butterfly and direct FIO evaluation in both 2D and 3D [1411.7418].

## 6. High-Dimensional and Tensor Extensions

Recent algorithmic developments generalize the multiscale butterfly structure to tensors, yielding linear $\mathcal{O}(n^d)$ scaling for $2d$-mode tensors in $d$ dimensions [2411.03029]. 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 $512\lambda$ per side), this approach enables $200\times$ speedup and $30\times$ memory reduction versus matrix butterfly [2411.03029].

This tensor extension is the first to achieve linear (in $n^d$) 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 [2411.03029].

## 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 $m$ scales as $\mathcal{O}(\log \log N)$ 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 $m$ at all hierarchical levels [1703.01941]. The stability of the nested butterfly representation guarantees that complexity remains near-optimal, and that catastrophic error growth is avoided [1703.01941].


---

**References:**  
[1411.7418] "A Multiscale Butterfly Algorithm for Multidimensional Fourier Integral Operators"  
[1509.07925] "Multidimensional Butterfly Factorization"  
[1703.01941] "An analysis of a butterfly algorithm"  
[2411.03029] "A Linear-complexity Tensor Butterfly Algorithm for Compressing High-dimensional Oscillatory Integral Operators"

Source: https://www.emergentmind.com/topics/multiscale-butterfly-algorithm