Leslie Greengard

Andreas Klöckner, Alexander Barnett, Leslie Greengard et al.

Integral equation methods for the solution of partial differential equations, when coupled with suitable fast algorithms, yield geometrically flexible, asymptotically optimal and well-conditioned schemes in either interior or exterior domains. The practical application of these methods, however, requires the accurate evaluation of boundary integrals with singular, weakly singular or nearly singular kernels. Historically, these issues have been handled either by low-order product integration rules (computed semi-analytically), by singularity subtraction/cancellation, by kernel regularization and asymptotic analysis, or by the construction of special purpose "generalized Gaussian quadrature" rules. In this paper, we present a systematic, high-order approach that works for any singularity (including hypersingular kernels), based only on the assumption that the field induced by the integral operator is locally smooth when restricted to either the interior or the exterior. Discontinuities in the field across the boundary are permitted. The scheme, denoted QBX (quadrature by expansion), is easy to implement and compatible with fast hierarchical algorithms such as the fast multipole method. We include accuracy tests for a variety of integral operators in two dimensions on smooth and corner domains.

Kenneth L. Ho, Leslie Greengard

We present a fast direct solver for structured linear systems based on multilevel matrix compression. Using the recently developed interpolative decomposition of a low-rank matrix in a recursive manner, we embed an approximation of the original matrix into a larger, but highly structured sparse one that allows fast factorization and application of the inverse. The algorithm extends the Martinsson/Rokhlin method developed for 2D boundary integral equations and proceeds in two phases: a precomputation phase, consisting of matrix compression and factorization, followed by a solution phase to apply the matrix inverse. For boundary integral equations which are not too oscillatory, e.g., based on the Green's functions for the Laplace or low-frequency Helmholtz equations, both phases typically have complexity O(N) in two dimensions, where $N$ is the number of discretization points. In our current implementation, the corresponding costs in three dimensions are $O(N^{3/2})$ and $O(N \log N)$ for precomputation and solution, respectively. Extensive numerical experiments show a speedup of $\sim 100$ for the solution phase over modern fast multipole methods; however, the cost of precomputation remains high. Thus, the solver is particularly suited to problems where large numbers of iterations would be required. Such is the case with ill-conditioned linear systems or when the same system is to be solved with multiple right-hand sides. Our algorithm is implemented in Fortran and freely available.

Felipe Vico, Leslie Greengard, Miguel Ferrando

We introduce a fast algorithm for computing volume potentials - that is, the convolution of a translation invariant, free-space Green's function with a compactly supported source distribution defined on a uniform grid. The algorithm relies on regularizing the Fourier transform of the Green's function by cutting off the interaction in physical space beyond the domain of interest. This permits the straightforward application of trapezoidal quadrature and the standard FFT, with superalgebraic convergence for smooth data. Moreover, the method can be interpreted as employing a Nystrom discretization of the corresponding integral operator, with matrix entries which can be obtained explicitly and rapidly. This is of use in the design of preconditioners or fast direct solvers for a variety of volume integral equations. The method proposed permits the computation of any derivative of the potential, at the cost of an additional FFT.

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.

Eduardo Corona, Leslie Greengard, Manas Rachh et al.

We present a new derivation of a boundary integral equation (BIE) for simulating the three-dimensional dynamics of arbitrarily-shaped rigid particles of genus zero immersed in a Stokes fluid, on which are prescribed forces and torques. Our method is based on a single-layer representation and leads to a simple second-kind integral equation. It avoids the use of auxiliary sources within each particle that play a role in some classical formulations. We use a spectrally accurate quadrature scheme to evaluate the corresponding layer potentials, so that only a small number of spatial discretization points per particle are required. The resulting discrete sums are computed in $\mathcal{O}(n)$ time, where $n$ denotes the number of particles, using the fast multipole method (FMM). The particle positions and orientations are updated by a high-order time-stepping scheme. We illustrate the accuracy, conditioning and scaling of our solvers with several numerical examples.

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.

Andras Pataki, Leslie Greengard

In this paper, we describe a new class of fast solvers for separable elliptic partial differential equations in cylindrical coordinates $(r,θ,z)$ with free-space radiation conditions. By combining integral equation methods in the radial variable $r$ with Fourier methods in $θ$ and $z$, we show that high-order accuracy can be achieved in both the governing potential and its derivatives. A weak singularity arises in the Fourier transform with respect to $z$ that is handled with special purpose quadratures. We show how these solvers can be applied to the evaluation of the Coulomb collision operator in kinetic models of ionized gases.

Jun Lai, Motoki Kobayashi, Leslie Greengard

In this paper, we consider acoustic or electromagnetic scattering in two dimensions from an infinite three-layer medium with thousands of wavelength-size dielectric particles embedded in the middle layer. Such geometries are typical of microstructured composite materials, and the evaluation of the scattered field requires a suitable fast solver for either a single configuration or for a sequence of configurations as part of a design or optimization process. We have developed an algorithm for problems of this type by combining the Sommerfeld integral representation, high order integral equation discretization, the fast multipole method and classical multiple scattering theory. The efficiency of the solver is illustrated with several numerical experiments.

Alex H. Barnett, Leslie Greengard, Tom Hagstrom

We develop a high-order, explicit method for acoustic scattering in three space dimensions based on a combined-field time-domain integral equation. The spatial discretization, of Nyström type, uses Gaussian quadrature on panels combined with a special treatment of the weakly singular kernels arising in near-neighbor interactions. In time, a new class of convolution splines is used in a predictor-corrector algorithm. Experiments on a torus and a perturbed torus are used to explore the stability and accuracy of the proposed scheme. This involved around one thousand solver runs, at up to 8th order and up to around 20,000 spatial unknowns, demonstrating 5-9 digits of accuracy. In addition we show that parameters in the combined field formulation, chosen on the basis of analysis for the sphere and other convex scatterers, work well in these cases.

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.

Yuanxun Bao, Manas Rachh, Eric Keaveny et al.

We present a fluctuating boundary integral method (FBIM) for overdamped Brownian Dynamics (BD) of two-dimensional periodic suspensions of rigid particles of complex shape immersed in a Stokes fluid. We develop a novel approach for generating Brownian displacements that arise in response to the thermal fluctuations in the fluid. Our approach relies on a first-kind boundary integral formulation of a mobility problem in which a random surface velocity is prescribed on the particle surface, with zero mean and covariance proportional to the Green's function for Stokes flow (Stokeslet). This approach yields an algorithm that scales linearly in the number of particles for both deterministic and stochastic dynamics, handles particles of complex shape, achieves high order of accuracy, and can be generalized to three dimensions and other boundary conditions. We show that Brownian displacements generated by our method obey the discrete fluctuation-dissipation balance relation (DFDB). Based on a recently-developed Positively Split Ewald method [A. M. Fiore, F. Balboa Usabiaga, A. Donev and J. W. Swan, J. Chem. Phys., 146, 124116, 2017], near-field contributions to the Brownian displacements are efficiently approximated by iterative methods in real space, while far-field contributions are rapidly generated by fast Fourier-space methods based on fluctuating hydrodynamics. FBIM provides the key ingredient for time integration of the overdamped Langevin equations for Brownian suspensions of rigid particles. We demonstrate that FBIM obeys DFDB by performing equilibrium BD simulations of suspensions of starfish-shaped bodies using a random finite difference temporal integrator.

Carlos Borges, Leslie Greengard

We consider the problem of reconstructing the shape of an impenetrable sound-soft obstacle from scattering measurements. The input data is assumed to be the far-field pattern generated when a plane wave impinges on an unknown obstacle from one or more directions and at one or more frequencies. It is well known that this inverse scattering problem is both ill posed and nonlinear. It is common practice to overcome the ill posedness through the use of a penalty method or Tikhonov regularization. Here, we present a more physical regularization, based simply on restricting the unknown boundary to be band-limited in a suitable sense. To overcome the nonlinearity of the problem, we use a variant of Newton's method. When multiple frequency data is available, we supplement Newton's method with the recursive linearization approach due to Chen. During the course of solving the inverse problem, we need to compute the solution to a large number of forward scattering problems. For this, we use high-order accurate integral equation discretizations, coupled with fast direct solvers when the problem is sufficiently large.

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.

Kenneth L. Ho, Leslie Greengard

We present a fast algorithm for linear least squares problems governed by hierarchically block separable (HBS) matrices. Such matrices are generally dense but data-sparse and can describe many important operators including those derived from asymptotically smooth radial kernels that are not too oscillatory. The algorithm is based on a recursive skeletonization procedure that exposes this sparsity and solves the dense least squares problem as a larger, equality-constrained, sparse one. It relies on a sparse QR factorization coupled with iterative weighted least squares methods. In essence, our scheme consists of a direct component, comprised of matrix compression and factorization, followed by an iterative component to enforce certain equality constraints. At most two iterations are typically required for problems that are not too ill-conditioned. For an $M \times N$ HBS matrix with $M \geq N$ having bounded off-diagonal block rank, the algorithm has optimal $\mathcal{O} (M + N)$ complexity. If the rank increases with the spatial dimension as is common for operators that are singular at the origin, then this becomes $\mathcal{O} (M + N)$ in 1D, $\mathcal{O} (M + N^{3/2})$ in 2D, and $\mathcal{O} (M + N^{2})$ in 3D. We illustrate the performance of the method on both over- and underdetermined systems in a variety of settings, with an emphasis on radial basis function approximation and efficient updating and downdating.

Charles L. Epstein, Zydrunas Gimbutas, Leslie Greengard et al.

A classical problem in electromagnetics concerns the representation of the electric and magnetic fields in the low-frequency or static regime, where topology plays a fundamental role. For multiply connected conductors, at zero frequency the standard boundary conditions on the tangential components of the magnetic field do not uniquely determine the vector potential. We describe a (gauge-invariant) consistency condition that overcomes this non-uniqueness and resolves a longstanding difficulty in inverting the magnetic field integral equation.

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.

Jun Lai, Leslie Greengard, Michael O'Neil

Scattering from large, open cavity structures is of importance in a variety of electromagnetic applications. In this paper, we propose a new well conditioned integral equation for scattering from general open cavities embedded in an infinite, perfectly conducting half-space. The integral representation permits the stable evaluation of both the electric and magnetic field, even in the low-frequency regime, using the continuity equation in a post-processing step. We establish existence and uniqueness results, and demonstrate the performance of the scheme in the cavity-of-revolution case. High-order accuracy is obtained using a Nystrom discretization with generalized Gaussian quadratures.

Lise-Marie Imbert-Gerard, Leslie Greengard

The inversion of the Laplace-Beltrami operator and the computation of the Hodge decomposition of a tangential vector field on smooth surfaces arise as computational tasks in many areas of science, from computer graphics to machine learning to com- putational physics. Here, we present a high-order accurate pseudo-spectral approach, applicable to closed surfaces of genus one in three dimensional space, with a view toward applications in plasma physics and fluid dynamics.

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.

Ryan Denlinger, Zydrunas Gimbutas, Leslie Greengard et al.

We present a fast summation method for lattice sums of the type which arise when solving wave scattering problems with periodic boundary conditions. While there are a variety of effective algorithms in the literature for such calculations, the approach presented here is new and leads to a rigorous analysis of Wood's anomalies. These arise when illuminating a grating at specific combinations of the angle of incidence and the frequency of the wave, for which the lattice sums diverge. They were discovered by Wood in 1902 as singularities in the spectral response. The primary tools in our approach are the Euler-Maclaurin formula and a steepest descent argument. The resulting algorithm has super-algebraic convergence and requires only milliseconds of CPU time.

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.

Alex H. Barnett, Charles L. Epstein, Leslie Greengard et al.

We study the stability properties of explicit marching schemes for second-kind Volterra integral equations that arise when solving boundary value problems for the heat equation by means of potential theory. It is well known that explicit finite difference or finite element schemes for the heat equation are stable only if the time step $Δt$ is of the order $O(Δx^2)$, where $Δx$ is the finest spatial grid spacing. In contrast, for the Dirichlet and Neumann problems on the unit ball in all dimensions $d\ge 1$, we show that the simplest Volterra marching scheme, i.e., the forward Euler scheme, is unconditionally stable. Our proof is based on an explicit spectral radius bound of the marching matrix, leading to an estimate that an $L^2$-norm of the solution to the integral equation is bounded by $c_dT^{d/2}$ times the norm of the right hand side. For the Robin problem on the half space in any dimension, with constant Robin (heat transfer) coefficient $κ$, we exhibit a constant $C$ such that the forward Euler scheme is stable if $Δt < C/κ^2$, independent of any spatial discretization. This relies on new lower bounds on the spectrum of real symmetric Toeplitz matrices defined by convex sequences. Finally, we show that the forward Euler scheme is unconditionally stable for the Dirichlet problem on any smooth convex domain in any dimension, in $L^\infty$-norm.

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.

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.

Manas Rachh, Leslie Greengard

We present new integral representations in two dimensions for the elastance problem in electrostatics and the mobility problem in Stokes flow. These representations lead to resonance-free Fredholm integral equations of the second kind and well conditioned linear systems upon discretization. By coupling our integral equations with high order quadrature and fast multipole acceleration, large-scale problems can be solved with only modest computing resources. We also discuss some applications of these boundary value problems in applied physics.

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.

Alex H. Barnett, Leslie Greengard

In this paper, we consider band-structure calculations governed by the Helmholtz or Maxwell equations in piecewise homogeneous periodic materials. Methods based on boundary integral equations are natural in this context, since they discretize the interface alone and can achieve high order accuracy in complicated geometries. In order to handle the quasi-periodic conditions which are imposed on the unit cell, the free-space Green's function is typically replaced by its quasi-periodic cousin. Unfortunately, the quasi-periodic Green's function diverges for families of parameter values that correspond to resonances of the empty unit cell. Here, we bypass this problem by means of a new integral representation that relies on the free-space Green's function alone, adding auxiliary layer potentials on the boundary of the unit cell itself. An important aspect of our method is that by carefully including a few neighboring images, the densities may be kept smooth and convergence rapid. This framework results in an integral equation of the second kind, avoids spurious resonances, and achieves spectral accuracy. Because of our image structure, inclusions which intersect the unit cell walls may be handled easily and automatically. Our approach is compatible with fast-multipole acceleration, generalizes easily to three dimensions, and avoids the complication of divergent lattice sums.

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).