A Butterfly-Accelerated Manifold Harmonic Transform

Problem Statement
The Fast Fourier Transform (FFT) and Nonuniform FFT (NUFFT) have revolutionized spectral analysis on simple, periodic domains like the circle and sphere by accelerating discrete trigonometric sums from $O(n^2)$ to $O(n \log n)$. These algorithms rely on the specific algebraic structure of complex exponentials, which are the eigenfunctions of the Laplacian on flat, periodic domains. However, for general compact manifolds $M$, spectral analysis requires the use of Laplace-Beltrami eigenfunctions (manifold harmonics). Unlike the Fourier case, these eigenfunctions $\phi_k$ generally lack closed-form expressions, and no general fast algorithm exists for evaluating linear combinations of these functions at arbitrary points. Existing methods for specific geometries (e.g., spheres) rely on symmetries and analytical expansions not available for arbitrary surfaces. Consequently, computing the Manifold Harmonic Transform (MHT), defined as $f(x_j) = \sum_{k=1}^m c_k \phi_k(x_j)$, typically requires $O(nm)$ operations and storage, which becomes prohibitive for large-scale problems in computer graphics, statistics, and machine learning.

Methodology
The authors propose a Butterfly-Accelerated Manifold Harmonic Transform (BF-MHT). This approach generalizes the hierarchical algebraic structure of the FFT to arbitrary operators by exploiting the "complementary low-rank" property of oscillatory kernels.

1. Butterfly Factorization: The method constructs a low-rank factorization of the dense MHT matrix $\Phi \in \mathbb{C}^{n \times m}$ (where $\Phi_{jk} = \phi_k(x_j)$). It utilizes two binary trees: a "space tree" $T_x$ on the spatial indices (points on the manifold) and a "frequency tree" $T_\lambda$ on the eigenvalue indices. By traversing these trees simultaneously (moving down the space tree and up the frequency tree), the algorithm identifies submatrices that are approximately low-rank. These blocks are compressed using transfer matrices and basis vectors, reducing the storage and application cost.
2. Fiedler Trees: A critical component is the construction of the space tree $T_x$. Standard spatial partitioning (e.g., octrees) often fails on manifolds because points close in Euclidean distance may be far in geodesic distance, leading to poor compressibility. The authors employ Fiedler trees, which recursively partition the manifold based on the nodal domains of the second eigenfunction ($\phi_2$) of the Laplace-Beltrami operator. This approach respects the intrinsic geometry of the manifold, ensuring that subdomains have approximately equal surface areas and minimal connection, thereby maximizing the low-rank property of the resulting matrix blocks.
3. Streaming Algorithm: To address memory constraints where the full dense matrix $\Phi$ cannot fit in RAM, the authors implement a streaming factorization. By traversing the frequency tree in post-order, the algorithm computes and compresses eigenfunctions in small batches, discarding intermediate data immediately. This ensures memory usage remains bounded by the size of the final compressed representation.

Key Contributions

• Generalization of FFT: The paper extends fast spectral transforms from simple periodic domains to arbitrary compact manifolds by utilizing Laplace-Beltrami eigenfunctions.
• Discretization-Agnostic Framework: The BF-MHT is applicable regardless of the underlying discretization method (finite elements, spectral collocation, radial basis functions, or mesh-free point clouds).
• Theoretical Complexity Analysis: The authors provide a rigorous analysis for the special case of the flat torus ($M = [-\pi, \pi]^2$). They prove that the $\epsilon$-rank of the annulus-to-disk Fourier kernel scales as $O(bR)$, leading to an asymptotic storage and application complexity of $O(n + m^{5/3})$ for the BF-MHT.
• Empirical Scaling: Numerical experiments on various geometries (flat torus, deformed torus, dragon mesh, Lake Superior) demonstrate that for moderate problem sizes, the empirical scaling is often closer to $O(n + m^{3/2})$, offering significant memory and computational savings over the $O(nm)$ dense approach.

Results
The paper presents numerical experiments verifying the algorithm's performance:

• Memory Reduction: On a dragon mesh with $n \approx 460,000$ degrees of freedom and $m \approx 5,600$ eigenfunctions, the dense matrix requires $\sim 19$ GB of RAM. The BF-MHT reduces this to 1.33 GB (for $\epsilon=10^{-6}$) and 533 MB (for $\epsilon=10^{-3}$), achieving compression factors of $\sim 14$ and $\sim 37$, respectively.
• Scalability: Experiments on the flat torus and deformed torus show a clear $\sim n^{3/2}$ scaling for memory usage, consistent with the empirical $O(n + m^{3/2})$ complexity.
• Applications: The algorithm is successfully applied to:
  • Spectral Geometry Processing: Filtering and modifying mesh coordinates via inverse MHT and spectral filtering.
  • Gaussian Random Fields (GRFs): Rapidly sampling from GRFs on manifolds by applying the transform to random coefficients.
  • Laplacian Eigenmaps: Compressing eigenvectors of weighted graph Laplacians for dimensionality reduction on noisy point clouds.

Significance and Claims
The authors claim that the BF-MHT provides a practical, fast algorithm for computing linear combinations of manifold harmonics, a task previously limited to $O(nm)$ costs. They emphasize that while the theoretical bound for the flat torus is $O(n + m^{5/3})$, the empirical performance on general manifolds suggests even better scaling ($O(n + m^{3/2})$) for practical problem sizes. The work bridges the gap between spectral analysis on simple domains and complex geometries, enabling efficient spectral processing in computer graphics, uncertainty quantification, and machine learning. The authors note that the precomputation of eigenfunctions remains a bottleneck ($O(nm)$), but the BF-MHT significantly reduces the cost of storing and applying the transform once these functions are available. They acknowledge that the discrepancy between the theoretical $O(m^{5/3})$ and empirical $O(m^{3/2})$ scaling is an open question for future work.