Leslie Greengard

Jiuyang Liang, Libin Lu, Alex Barnett et al.

The evaluation of long-range Coulomb interactions is a significant cost in molecular dynamics (MD), even when using Particle Mesh Ewald (PME) or Particle-Particle-Particle-Mesh (PPPM) methods, which rely on Ewald splitting and the fast Fourier transform to achieve near-linear scaling. We introduce ESP -- Ewald summation with prolate spheroidal wave functions (PSWFs) -- which leads to a more efficient Fourier representation and a reduction in the required grid size, global communication, and particle-grid operations, without loss of accuracy. We have integrated the ESP method into two widely-used open-source MD packages, LAMMPS and GROMACS, enabling rapid comparison and adoption. Relative to PME/PPPM baselines at error tolerances $10^{-3}$ to $10^{-4}$, ESP gives roughly a $3$-fold acceleration of electrostatic interactions, and a $2.5$-fold speed-up in the MD simulation when using about $10^3$ compute cores. At high accuracy ($10^{-5}$), these increase to $10$-fold for the far-field electrostatics and $5$-fold for MD simulation. Furthermore, we show that the accelerated codes have improved strong scaling with core count, and validate them in realistic long-time biological and material simulations. ESP thus offers a practical, drop-in path to reduce the time-to-solution and energy footprint of MD workflows.

Charles L. Epstein, Leslie Greengard, Michael O'Neil

In this paper, we develop a new integral representation for the solution of the time harmonic Maxwell equations in media with piecewise constant dielectric permittivity and magnetic permeability in R^3. This representation leads to a coupled system of Fredholm integral equations of the second kind for four scalar densities supported on the material interface. Like the classical Muller equation, it has no spurious resonances. Unlike the classical approach, however, the representation does not suffer from low frequency breakdown. We illustrate the performance of the method with numerical examples.

Charles L. Epstein, Leslie Greengard

In this paper, we develop a new representation for outgoing solutions to the time harmonic Maxwell equations in unbounded domains in $\bbR^3.$ This representation leads to a Fredholm integral equation of the second kind for solving the problem of scattering from a perfect conductor, which does not suffer from spurious resonances or low frequency breakdown, although it requires the inversion of the scalar surface Laplacian on the domain boundary. In the course of our analysis, we give a new proof of the existence of non-trivial families of time harmonic solutions with vanishing normal components that arise when the boundary of the domain is not simply connected. We refer to these as $k$-Neumann fields, since they generalize, to non-zero wave numbers, the classical harmonic Neumann fields. The existence of $k$-harmonic fields was established earlier by Kress.

Zydrunas Gimbutas, Leslie Greengard

A variety of problems in device and materials design require the rapid forward modeling of Maxwell's equations in complex micro-structured materials. By combining high-order accurate integral equation methods with classical multiple scattering theory, we have created an effective simulation tool for materials consisting of an isotropic background in which are dispersed a large number of micro- or nano-scale metallic or dielectric inclusions.

Charles L. Epstein, Leslie Greengard, Andreas Klöckner

In a recently developed quadrature method (quadrature by expansion or QBX), it was demonstrated that weakly singular or singular layer potentials can be evaluated rapidly and accurately on surface by making use of local expansions about carefully chosen off-surface points. In this paper, we derive estimates for the rate of convergence of these local expansions, providing the analytic foundation for the QBX method. The estimates may also be of mathematical interest, particularly for microlocal or asymptotic analysis in potential theory.

Michael O'Neil, Leslie Greengard, Andras Pataki

A classical problem in acoustic (and electromagnetic) scattering concerns the evaluation of the Green's function for the Helmholtz equation subject to impedance boundary conditions on a half-space. The two principal approaches used for representing this Green's function are the Sommerfeld integral and the (closely related) method of complex images. The former is extremely efficient when the source is at some distance from the half-space boundary, but involves an unwieldy range of integration as the source gets closer and closer. Complex image-based methods, on the other hand, can be quite efficient when the source is close to the boundary, but they do not easily permit the use of the superposition principle since the selection of complex image locations depends on both the source and the target. We have developed a new, hybrid representation which uses a finite number of real images (dependent only on the source location) coupled with a rapidly converging Sommerfeld-like integral. While our method applies in both two and three dimensions, we restrict the detailed analysis and numerical experiments here to the two-dimensional case.

Andras Pataki, Antoine J. Cerfon, Jeffrey P. Freidberg et al.

We present a new fast solver to calculate fixed-boundary plasma equilibria in toroidally axisymmetric geometries. By combining conformal mapping with Fourier and integral equation methods on the unit disk, we show that high-order accuracy can be achieved for the solution of the equilibrium equation and its first and second derivatives. Smooth arbitrary plasma cross-sections as well as arbitrary pressure and poloidal current profiles are used as initial data for the solver. Equilibria with large Shafranov shifts can be computed without difficulty. Spectral convergence is demonstrated by comparing the numerical solution with a known exact analytic solution. A fusion-relevant example of an equilibrium with a pressure pedestal is also presented.

Charles L. Epstein, Leslie Greengard, Michael O'Neil

The generalized Debye source representation of time-harmonic electromagnetic fields yields well-conditioned second-kind integral equations for a variety of boundary value problems, including the problems of scattering from perfect electric conductors and dielectric bodies. Furthermore, these representations, and resulting integral equations, are fully stable in the static limit as $ω\to 0$ in multiply connected geometries. In this paper, we present the first high-order accurate solver based on this representation for bodies of revolution. The resulting solver uses a Nyström discretization of a one-dimensional generating curve and high-order integral equation methods for applying and inverting surface differentials. The accuracy and speed of the solvers are demonstrated in several numerical examples.

Jun Wang, Leslie Greengard

A variety of problems in computational physics and engineering require the convolution of the heat kernel (a Gaussian) with either discrete sources, densities supported on boundaries, or continuous volume distributions. We present a unified fast Gauss transform for this purpose in two dimensions, making use of an adaptive quad-tree discretization on a unit square which is assumed to contain all sources. Our implementation permits either free-space or periodic boundary conditions to be imposed, and is efficient for any choice of variance in the Gaussian.

Lise-Marie Imbert-Gerard, Felipe Vico, Leslie Greengard et al.

We present a collection of well-conditioned integral equation methods for the solution of electrostatic, acoustic or electromagnetic scattering problems involving anisotropic, inhomogeneous media. In the electromagnetic case, our approach involves a minor modification of a classical formulation. In the electrostatic or acoustic setting, we introduce a new vector partial differential equation, from which the desired solution is easily obtained. It is the vector equation for which we derive a well-conditioned integral equation. In addition to providing a unified framework for these solvers, we illustrate their performance using iterative solution methods coupled with the FFT-based technique of [1] to discretize and apply the relevant integral operators.

Jun Wang, Leslie Greengard

We present a hybrid asymptotic/numerical method for the accurate computation of single and double layer heat potentials in two dimensions. It has been shown in previous work that simple quadrature schemes suffer from a phenomenon called "geometrically-induced stiffness," meaning that formally high-order accurate methods require excessively small time steps before the rapid convergence rate is observed. This can be overcome by analytic integration in time, requiring the evaluation of a collection of spatial boundary integral operators with non-physical, weakly singular kernels. In our hybrid scheme, we combine a local asymptotic approximation with the evaluation of a few boundary integral operators involving only Gaussian kernels, which are easily accelerated by a new version of the fast Gauss transform. This new scheme is robust, avoids geometrically-induced stiffness, and is easy to use in the presence of moving geometries. Its extension to three dimensions is natural and straightforward, and should permit layer heat potentials to become flexible and powerful tools for modeling diffusion processes.

Leslie Greengard, Shidong Jiang, Jun Wang

Two fundamental difficulties are encountered in the numerical evaluation of time-dependent layer potentials. One is the quadratic cost of history dependence, which has been successfully addressed by splitting the potentials into two parts - a local part that contains the most recent contributions and a history part that contains the contributions from all earlier times. The history part is smooth, easily discretized using high-order quadratures, and straightforward to compute using a variety of fast algorithms. The local part, however, involves complicated singularities in the underlying Green's function. Existing methods, based on exchanging the order of integration in space and time, are able to achieve high order accuracy, but are limited to the case of stationary boundaries. Here, we present a new quadrature method that leaves the order of integration unchanged, making use of a change of variables that converts the singular integrals with respect to time into smooth ones. We have also derived asymptotic formulas for the local part that lead to fast and accurate hybrid schemes, extending earlier work for scalar heat potentials and applicable to moving boundaries. The performance of the overall scheme is demonstrated via numerical examples.

Nour G. Al Hassanieh, Leslie Greengard, Alex H. Barnett

We present a fast algorithm for evaluating the (non-smooth) solution of the free-space two-dimensional (2D) scalar wave equation with many point sources, each with a high-frequency band-limited time signature. Such an algorithm is key to an efficient time-domain scattering solver using spatially-discretized hyperbolic layer potentials. Given $M$ sources/targets and $N_t$ time steps, direct evaluation costs $O(M^2N_t^2)$, due to the history dependence. We develop a quasi-linear scaling algorithm that splits the solution at a given time into (a) a non-smooth time-local part, (b) a (smooth) near history involving sources up to ${\mathcal O}(1)$ domain traversal times into the past, plus (c) a (very smooth) far history comprising all waves emitted before the near history. The local part is computed directly via high-order quadrature. A naive spatial Fourier transform for (b) plus (c) would be both slowly converging and arbitrarily oscillatory as time progresses. Yet in (b) the oscillations are controlled, so we use the recent truncated windowed Fourier projection (TK-WFP) method to give rapid convergence. For (c) -- present due to the weak Huygens' principle -- we exploit a new large-time sum-of-exponentials approximation of the free-space wave kernel. Numerical examples with up to a million sources and targets, a domain of $300\times 300$ wavelengths, and 6-digit accuracy, show an acceleration of five orders of magnitude relative to direct evaluation.

Leslie Greengard, Shidong Jiang

We present a new integral representation for the unsteady, incompressible Stokes or Navier-Stokes equations, based on a linear combination of heat and harmonic potentials. For velocity boundary conditions, this leads to a coupled system of integral equations: one for the normal component of velocity and one for the tangential components. Each individual equation is well-condtioned, and we show that using them in predictor-corrector fashion, combined with spectral deferred correction, leads to high-order accuracy solvers. The fundamental unknowns in the mixed potential representation are densities supported on the boundary of the domain. We refer to one as the vortex source, the other as the pressure source and the coupled system as the combined source integral equation.

Alex Barnett, Leslie Greengard, Andras Pataki et al.

Determining the three-dimensional structure of proteins and protein complexes at atomic resolution is a fundamental task in structural biology. Over the last decade, remarkable progress has been made using "single particle" cryo-electron microscopy (cryo-EM) for this purpose. In cryo-EM, hundreds of thousands of two-dimensional images are obtained of individual copies of the same particle, each held in a thin sheet of ice at some unknown orientation. Each image corresponds to the noisy projection of the particle's electron-scattering density. The reconstruction of a high-resolution image from this data is typically formulated as a nonlinear, non-convex optimization problem for unknowns which encode the angular pose and lateral offset of each particle. Since there are hundreds of thousands of such parameters, this leads to a very CPU-intensive task---limiting both the number of particle images which can be processed and the number of independent reconstructions which can be carried out for the purpose of statistical validation. Here, we propose a deterministic method for high-resolution reconstruction that operates in an ab initio manner---that is, without the need for an initial guess. It requires a predictable and relatively modest amount of computational effort, by marching out radially in the Fourier domain from low to high frequency, increasing the resolution by a fixed increment at each step.

Carlos Borges, Adrianna Gillman, Leslie Greengard

We describe a fast, stable algorithm for the solution of the inverse acoustic scattering problem in two dimensions. Given full aperture far field measurements of the scattered field for multiple angles of incidence, we use Chen's method of recursive linearization to reconstruct an unknown sound speed at resolutions of thousands of square wavelengths in a fully nonlinear regime. Despite the fact that the underlying optimization problem is formally ill-posed and non-convex, recursive linearization requires only the solution of a sequence of linear least squares problems at successively higher frequencies. By seeking a suitably band-limited approximation of the sound speed profile, each least squares calculation is well-conditioned and involves the solution of a large number of forward scattering problems, for which we employ a recently developed, spectrally accurate, fast direct solver. For the largest problems considered, involving 19,600 unknowns, approximately one million partial differential equations were solved, requiring approximately two days to compute using a parallel MATLAB implementation on a multi-core workstation.

Alex H. Barnett, Jeremy F. Magland, Leslie F. Greengard

We describe a suite of validation metrics that assess the credibility of a given automatic spike sorting algorithm applied to a given electrophysiological recording, when ground-truth is unavailable. By rerunning the spike sorter two or more times, the metrics measure stability under various perturbations consistent with variations in the data itself, making no assumptions about the noise model, nor about the internal workings of the sorting algorithm. Such stability is a prerequisite for reproducibility of results. We illustrate the metrics on standard sorting algorithms for both in vivo and ex vivo recordings. We believe that such metrics could reduce the significant human labor currently spent on validation, and should form an essential part of large-scale automated spike sorting and systematic benchmarking of algorithms.

Jun Lai, Leslie Greengard, Michael O'Neil

A variety of problems in acoustic and electromagnetic scattering require the evaluation of impedance or layered media Green's functions. Given a point source located in an unbounded half-space or an infinitely extended layer, Sommerfeld and others showed that Fourier analysis combined with contour integration provides a systematic and broadly effective approach, leading to what is generally referred to as the Sommerfeld integral representation. When either the source or target is at some distance from an infinite boundary, the number of degrees of freedom needed to resolve the scattering response is very modest. When both are near an interface, however, the Sommerfeld integral involves a very large range of integration and its direct application becomes unwieldy. Historically, three schemes have been employed to overcome this difficulty: the method of images, contour deformation, and asymptotic methods of various kinds. None of these methods make use of classical layer potentials in physical space, despite their advantages in terms of adaptive resolution and high-order accuracy. The reason for this is simple: layer potentials are impractical in layered media or half-space geometries since they require the discretization of an infinite boundary. In this paper, we propose a hybrid method which combines layer potentials (physical-space) on a finite portion of the interface together with a Sommerfeld-type (Fourier) correction. We prove that our method is efficient and rapidly convergent for arbitrarily located sources and targets, and show that the scheme is particularly effective when solving scattering problems for objects which are close to the half-space boundary or even embedded across a layered media interface.

Sivaram Ambikasaran, Carlos Borges, Lise-Marie Imbert-Gerard et al.

We present a fast direct solver for two dimensional scattering problems, where an incident wave impinges on a penetrable medium with compact support. We represent the scattered field using a volume potential whose kernel is the outgoing Green's function for the exterior domain. Inserting this representation into the governing partial differential equation, we obtain an integral equation of the Lippmann-Schwinger type. The principal contribution here is the development of an automatically adaptive, high-order accurate discretization based on a quad tree data structure which provides rapid access to arbitrary elements of the discretized system matrix. This permits the straightforward application of state-of-the-art algorithms for constructing compressed versions of the solution operator. These solvers typically require $O(N^{3/2})$ work, where $N$ denotes the number of degrees of freedom. We demonstrate the performance of the method for a variety of problems in both the low and high frequency regimes.

Charles L. Epstein, Leslie Greengard, Thomas Hagstrom

We give a principled approach for the selection of a boundary integral, retarded potential representation for the solution of scattering problems for the wave equation in an exterior domain.

Sivaram Ambikasaran, Daniel Foreman-Mackey, Leslie Greengard et al.

A number of problems in probability and statistics can be addressed using the multivariate normal (Gaussian) distribution. In the one-dimensional case, computing the probability for a given mean and variance simply requires the evaluation of the corresponding Gaussian density. In the $n$-dimensional setting, however, it requires the inversion of an $n \times n$ covariance matrix, $C$, as well as the evaluation of its determinant, $\det(C)$. In many cases, such as regression using Gaussian processes, the covariance matrix is of the form $C = σ^2 I + K$, where $K$ is computed using a specified covariance kernel which depends on the data and additional parameters (hyperparameters). The matrix $C$ is typically dense, causing standard direct methods for inversion and determinant evaluation to require $\mathcal O(n^3)$ work. This cost is prohibitive for large-scale modeling. Here, we show that for the most commonly used covariance functions, the matrix $C$ can be hierarchically factored into a product of block low-rank updates of the identity matrix, yielding an $\mathcal O (n\log^2 n) $ algorithm for inversion. More importantly, we show that this factorization enables the evaluation of the determinant $\det(C)$, permitting the direct calculation of probabilities in high dimensions under fairly broad assumptions on the kernel defining $K$. Our fast algorithm brings many problems in marginalization and the adaptation of hyperparameters within practical reach using a single CPU core. The combination of nearly optimal scaling in terms of problem size with high-performance computing resources will permit the modeling of previously intractable problems. We illustrate the performance of the scheme on standard covariance kernels.

Leslie Greengard, Chris Stucchio

A recurring task in image processing, approximation theory, and the numerical solution of partial differential equations is to reconstruct a piecewise-smooth real-valued function f(x) in multiple dimensions from its truncated Fourier transform (its truncated spectrum). An essential step is edge detection for which a variety of one-dimensional schemes have been developed over the last few decades. Most higher-dimensional edge detection algorithms consist of applying one-dimensional detectors in each component direction in order to recover the locations in R^N where f(x) is singular (the singular support). In this paper, we present a multidimensional algorithm which identifies the wavefront of a function from spectral data. The wavefront of f(x) is the set of points $(x,k) \in R^N \times (S^{N-1} / \{\pm 1\})$ which encode both the location of the singular points of a function and the orientation of the singularities. (Here $S^{N-1}$ denotes the unit sphere in N dimensions.) More precisely, k is the direction of the normal line to the curve or surface of discontinuity at x. Note that the singular support is simply the projection of the wavefront onto its x-component. In one dimension, the wavefront is a subset of R^1 and coincides with the singular support. In higher dimensions, geometry comes into play and they are distinct. We discuss the advantages of wavefront reconstruction and indicate how it can be used for segmentation in magnetic resonance imaging (MRI).