Source-linked AI summary
A kernel-independent precorrected-FFT method
Tristan Goodwill, Owen Melia
TL;DR
N-body summation methods have often been limited by their reliance on analytic kernel information. This paper adapts precorrected-FFT using proxy shells and pointwise kernel access, and reports favorable computational qualities and large speedups in integral equation solvers.
Problem
Fast kernel-summation methods rely on analytic properties of the kernel, limiting their direct applicability to general kernels known only through pointwise access.
Source: §1 · p. 1Method
The method uses proxy shells to construct spreading operators and automatically chooses proxy points and spreading grids using pointwise kernel access.
Source: §1 · p. 2, §5 · p. 21Results
The method shares grid-based Ewald's favorable computational qualities and leads to large speedups in integral equation solvers.
Source: §1 · p. 2Takeaways & Limitations
Automatic parameter selection is intended to help when an operator's Fourier transform is complicated or the kernel changes during runtime.
Source: §1 · p. 2The work is limited to scalar-valued kernels; efficient implementations for vector-valued problems remain open.
Source: §5 · p. 21
Abstract
from arXiv · showhide
The precorrected-FFT has been used to rapidly evaluate $N$-body sums arising in electromagnetics and acoustics settings. It has, however, been historically limited to kernels that satisfy a Green's identity. In this work, we demonstrate that, using an improved version of the method of proxy points, the precorrected-FFT can be used to efficiently compute $N$-body sums with a wide variety of kernels. This algorithm requires no analytic information beyond pointwise access to the kernel. Our examples include kernels from non-local PDEs, Matern kernels, anisotropic kernels, and the Green's functions of high order PDEs. We also provide a detailed algorithm for choosing the proxy points and spreading grid, as implemented in the new pcfft package. We show these algorithmic choices result in favorable scaling when compared to other kernel-independent fast summation methods. Finally, we use the algorithm to accelerate a wide variety of integral equations solvers involving non-standard kernels, with applications in surface scattering, geophysics, and microbial community formation.
1 Introduction
Direct evaluation of N-body sums costs O(N^2), while fast methods often depend on analytic kernel properties. This work makes precorrected-FFT kernel-independent, using pointwise kernel access to choose grid and Fourier parameters automatically.
1 Introduction
- Direct evaluation costs O(N^2), and kernels can be singular near source-target pairs while decaying slowly over distance.These features make global approximation difficult while leaving many interactions to account for.Source: §1 · p. 1
- The kernel-independent precorrected-FFT uses pointwise kernel access and equivalent grid charges, with nearby interactions corrected directly.The approach is proposed for a broad class of kernels and general point clouds in two or three dimensions.Source: §1 · p. 1
- Unlike fast multipole methods that rely on analytic kernel structure, the method automatically chooses gridpoints and Fourier modes for a user-specified error tolerance.Grid-based Ewald methods also require a priori kernel information to determine Fourier modes and grid size.Source: §1 · p. 2
- The method shares grid-based Ewald's O(N log(N)) scaling and favorable computational qualities, which lead to large speedups in integral equation solvers.The stated advantages include small constant factors and optimized FFT software.Source: §1 · p. 2
2 Algorithm description
This section describes a proxy-point, FFT-based method for summing kernels without a Green’s identity, including corrections for nearby interactions. It also tests the accuracy of kernel derivatives in a flexural-wave integral equation with clamped plate boundary conditions.
2 Algorithm description
- The algorithm spreads charges to nearby regular-grid points, applies the grid kernel sum through a zero-padded FFT, then uses the adjoint spreading operation.The grid kernel matrix has multilevel Toeplitz structure; sparse corrections restore interactions for pairs that are not well separated.Source: §2 · p. 3
2.1 Corrections for nearby targets
- Nearby interactions are corrected by subtracting the approximate contributions for spreading-bin pairs whose centers lie within 2ρ1 and adding their direct interactions.The correction matrix uses the target points in each bin’s neighborhood and source points in the bin; translation invariance lets the kernel template be precomputed once.Source: §2.1 · pp. 4–5
2.2 Finding the spreading parameters
- Spreading parameters are chosen empirically: resolve the kernel on the first proxy shell to set nproxy, then search nspread and nshell for the requested tolerance.For kernels without a Green’s identity, the method enforces agreement on multiple shells whose reciprocal radii are second-kind Chebyshev nodes; kernels with a Green’s identity use one shell.Source: §2.2 · p. 7
2.3 Proxy point tests
- For kernels without Green’s identities, proxy tests found shell counts grow roughly as |log ϵ| as tolerance shrinks, while nproxy and nspread vary with dimension and kernel frequency.The swapped-dimension tests also show slightly higher interaction rank than the Green’s-function tests.Source: §2.3 · p. 8, Figure 4 · p. 9
- For anisotropic Gaussian and other non-elliptic-PDE kernels, shell counts grew slowly and nproxy logarithmically with tighter tolerance, while the method attained high accuracy.The tested examples include the anisotropic Gaussian and f(r) = sin r/(r + r2 1).Source: §2.3 · p. 9, Figure 6 · p. 10
- Higher powers of the Laplacian required slightly larger nproxy and nspread; a single proxy shell sufficed when the method also matched n −1 radial derivatives for ∆n.The nproxy spike at the smallest ϵ came from finite-difference error in computing radial derivatives.Source: §2.3 · p. 9, Figure 7 · p. 10
2.4 Derivatives of the kernel
- For kernel derivatives, the approximation extends to derivatives when the kernel is bounded outside the bin and admits the low-rank factorization; experiments found errors still scale with tolerance.Derivatives can decrease accuracy, likely because the spread charges are larger.Source: §2.4 · pp. 10–11
- Accuracy scales with the tolerance for all tested derivative combinations, but is not always below the tolerance.The tests concern kernels in the flexural Green’s-function integral equation for the flexural wave equation with clamped plate boundary conditions.Source: §2.4 · p. 11
3 Complexity
The method has O(Npts log Npts) complexity, with FFTs dominating execution. Its application stage is especially fast after precomputation, while total runtime and solver performance depend on setup costs and problem size.
3 Complexity
- O(Npts log Npts) is the execution complexity, with the FFT/IFFT pair as the dominant cost; precomputation has the same complexity.A single application also uses sparse matrices with O(Npts · nnbr) entries.Source: §3 · p. 12
3.1 Scalability
- The 2D and 3D scaling experiments use single-thread execution, while the other examples use multithreading.The compute environment is a MacBook Pro with an Apple M4 Max processor.Source: §3.1 · p. 13
- After precomputation, PCFFT applies the 2D Helmholtz charge vector an order of magnitude faster than the comparison methods, though total runtime includes correction-matrix setup.The comparison uses error tolerance ϵ = 10−5.Source: §3.1 · p. 13, Figure 9 · p. 14
- After precomputation, PCFFT applies the 3D Helmholtz kernel nearly an order of magnitude faster than the other methods, while BBFMM3D could not accurately compress it.For 3D Laplace, optimized kernel-specific codes outperform the kernel-independent methods.Source: §3.1 · p. 14, Figure 11 · p. 15, Figure 10 · p. 15
- At Npts = 13,168, PCFFT is faster than the direct solver, but the direct solver is faster at the two largest sizes as GMRES requires 804 and 916 iterations.At Npts = 13,168, the iterative methods require 526 GMRES iterations.Source: §3.1 · p. 16, Figure 12 · p. 16
4 Examples
The examples apply PCFFT in eigenvalue computations, a surface Helmholtz-Beltrami solve, and lens optimization. They demonstrate its use in integral-equation applications with nonstandard kernels.
4.1 Eigendecomposition of the flexural clamped plate biharmonic problem
- PCFFT accelerates the clamped-plate eigenvalue search on three iceberg geometries, and the authors report solving the nonlinear eigenvalue problems efficiently.The search samples S(k), fits a Chebyshev series for 1/S(k), and identifies its zeros as resonant wavenumbers.Source: §4.1 · pp. 16–17
4.2 Surface Helmholtz
- The surface Helmholtz-Beltrami solve achieves four-digit accuracy on 350,865 points, with 365 GMRES iterations and 48 GB memory use.The discretization uses 7,797 eighth-order patches; each system-matrix application takes 1.50 seconds.Source: §4.2 · p. 18, Figure 14 · p. 18
4.3 Lens optimization
- For small q with (k0, k1) = (4, 6), the proposed representation requires 33 GMRES iterations versus 66 for the simpler representation.The authors report that fewer iterations outweigh the added boundary discretization and second small precorrected-FFT per iteration.Source: §4.3 · p. 20
- The optimized lens reaches intensity 5.1 at x0 = (4, 0), compared with 4.79 for the best constant wavenumber k ∈ [5, 7].The optimization stopped after 132 iterations; the discretization had 20,138 unknowns and took 2,041 seconds overall.Source: §4.3 · p. 20, Figure 16 · p. 21
5 Conclusion
The paper modifies precorrected-FFT for generic kernels using proxy shells and a systematic parameter-selection algorithm that needs no analytic kernel knowledge. The method is currently limited to scalar-valued kernels.
5 Conclusion
- Proxy shells build spreading operators that map point charges to equivalent grid charges, while a systematic algorithm chooses proxy points and spreading grids without analytic kernel knowledge.The authors use this kernel-independent method to accelerate integral-equation approaches for PDE problems.Source: §5 · p. 21
- The presented work is limited to scalar-valued kernels; efficient implementations for vector-valued problems remain open.The authors identify GPU implementations as a direction to explore for offsetting the extra cost.Source: §5 · p. 21