Alex Barnett

Lin Zhao, Alex Barnett

The drum problem-finding the eigenvalues and eigenfunctions of the Laplacian with Dirichlet boundary condition-has many applications, yet remains challenging for general domains when high accuracy or high frequency is needed. Boundary integral equations are appealing for large-scale problems, yet certain difficulties have limited their use. We introduce two ideas to remedy this: 1) We solve the resulting nonlinear eigenvalue problem using Boyd's method for analytic root-finding applied to the Fredholm determinant. We show that this is many times faster than the usual iterative minimization of a singular value. 2) We fix the problem of spurious exterior resonances via a combined field representation. This also provides the first robust boundary integral eigenvalue method for non-simply-connected domains. We implement the new method in two dimensions using spectrally accurate Nystrom product quadrature. We prove exponential convergence of the determinant at roots for domains with analytic boundary. We demonstrate 13-digit accuracy, and improved efficiency, in a variety of domain shapes including ones with strong exterior resonances.

Abtin Rahimian, Alex Barnett, Denis Zorin

We introduce a quadrature scheme--QBKIX--for the high-order accurate evaluation of layer potentials associated with general elliptic PDEs near to and on the domain boundary. Relying solely on point evaluations of the underlying kernel, our scheme is essentially PDE-independent; in particular, no analytic expansion nor addition theorem is required. Moreover, it applies to boundary integrals with singular, weakly singular, and hypersingular kernels. Our work builds upon Quadrature by Expansion (QBX), which approximates the potential by an analytic expansion in the neighborhood of each expansion center. In contrast, we use a sum of fundamental solutions lying on a ring enclosing the neighborhood, and solve a small dense linear system for their coefficients to match the potential on a smaller concentric ring. We test the new method with Laplace, Helmholtz, Yukawa, Stokes, and Navier (elastostatic) kernels in two dimensions (2D) using adaptive, panel-based boundary quadratures on smooth and corner domains. Advantages of the algorithm include its relative simplicity of implementation, immediate extension to new kernels, dimension-independence (allowing simple generalization to 3D), and compatibility with fast algorithms such as the kernel-independent FMM.

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.

Adrianna Gillman, Alex Barnett

We consider the numerical solution of the scattering of time-harmonic plane waves from an infinite periodic array of reflection or transmission obstacles in a homogeneous background medium, in two dimensions. Boundary integral formulations are ideal since they reduce the problem to $N$ unknowns on the obstacle boundary. However, for complex geometries and/or higher frequencies the resulting dense linear system becomes large, ruling out dense direct methods, and often ill-conditioned (despite being 2nd-kind), rendering fast multipole-based iterative schemes also inefficient. We present an integral equation based solver with O(N) complexity, which handles such ill-conditioning, using recent advances in "fast" direct linear algebra to invert hierarchically the isolated obstacle matrix. This is combined with a recent periodizing scheme that is robust for all incident angles, including Wood's anomalies, based upon the free space Green's function kernel. The resulting solver is extremely efficient when multiple incident angles are needed, as occurs in many applications. Our numerical tests include a complicated obstacle several wavelengths in size, with $N=10^5$ and solution error of $10^{-10}$, where the solver is 66 times faster per incident angle than a fast multipole based iterative solution, and 600 times faster when incident angles are chosen to share Bloch phases.

Alex Barnett, Andrew Hassell, Melissa Tacy

For smooth bounded domains in $\mathbb{R}$, we prove upper and lower $L^2$ bounds on the boundary data of Neumann eigenfunctions, and prove quasi-orthogonality of this boundary data in a spectral window. The bounds are tight in the sense that both are independent of eigenvalue; this is achieved by working with an appropriate norm for boundary functions, which includes a `spectral weight', that is, a function of the boundary Laplacian. This spectral weight is chosen to cancel concentration at the boundary that can happen for `whispering gallery' type eigenfunctions. These bounds are closely related to wave equation estimates due to Tataru. Using this, we bound the distance from an arbitrary Helmholtz parameter $E>0$ to the nearest Neumann eigenvalue, in terms of boundary normal-derivative data of a trial function $u$ solving the Helmholtz equation $(Δ-E)u=0$. This `inclusion bound' improves over previously known bounds by a factor of $E^{5/6}$. It is analogous to a recently improved inclusion bound in the Dirichlet case, due to the first two authors. Finally, we apply our theory to present an improved numerical implementation of the method of particular solutions for computation of Neumann eigenpairs on smooth planar domains. We show that the new inclusion bound improves the relative accuracy in a computed Neumann eigenvalue (around the $42000$th) from 9 digits to 14 digits, with little extra effort.

Chirag Modi, Alex Barnett, Bob Carpenter

The efficiency of Hamiltonian Monte Carlo (HMC) can suffer when sampling a distribution with a wide range of length scales, because the small step sizes needed for stability in high-curvature regions are inefficient elsewhere. To address this we present a delayed rejection variant: if an initial HMC trajectory is rejected, we make one or more subsequent proposals each using a step size geometrically smaller than the last. We extend the standard delayed rejection framework by allowing the probability of a retry to depend on the probability of accepting the previous proposal. We test the scheme in several sampling tasks, including multiscale model distributions such as Neal's funnel, and statistical applications. Delayed rejection enables up to five-fold performance gains over optimally-tuned HMC, as measured by effective sample size per gradient evaluation. Even for simpler distributions, delayed rejection provides increased robustness to step size misspecification. Along the way, we provide an accessible but rigorous review of detailed balance for HMC.

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.

Gary Marple, Alex Barnett, Adrianna Gillman et al.

This paper presents a new boundary integral equation (BIE) method for simulating particulate and multiphase flows through periodic channels of arbitrary smooth shape in two dimensions. The authors consider a particular system---multiple vesicles suspended in a periodic channel of arbitrary shape---to describe the numerical method and test its performance. Rather than relying on the periodic Green's function as classical BIE methods do, the method combines the free-space Green's function with a small auxiliary basis, and imposes periodicity as an extra linear condition. As a result, we can exploit existing free-space solver libraries, quadratures, and fast algorithms, and handle a large number of vesicles in a geometrically complex channel. Spectral accuracy in space is achieved using the periodic trapezoid rule and product quadratures, while a first-order semi-implicit scheme evolves particles by treating the vesicle-channel interactions explicitly. New constraint-correction formulas are introduced that preserve reduced areas of vesicles, independent of the number of time steps taken. By using two types of fast algorithms, (i) the fast multipole method (FMM) for the computation of the vesicle-vesicle and the vesicle-channel hydrodynamic interaction, and (ii) a fast direct solver for the BIE on the fixed channel geometry, the computational cost is reduced to $O(N)$ per time step where $N$ is the spatial discretization size. Moreover, the direct solver inverts the wall BIE operator at $t = 0$, stores its compressed representation and applies it at every time step to evolve the vesicle positions, leading to dramatic cost savings compared to classical approaches. Numerical experiments illustrate that a simulation with $N=128, 000$ can be evolved in less than a minute per time step on a laptop.

Jun Lai, Motoki Kobayashi, Alex Barnett

We present a solver for plane wave scattering from a periodic dielectric grating with a large number $M$ of inclusions lying in each period of its middle layer.Such composite material geometries have a growing role in modern photonic devices and solar cells. The high-order scheme is based on boundary integral equations, and achieves many digits of accuracy with ease. The usual way to periodize the integral equation---via the quasi-periodic Green's function---fails at Wood's anomalies. We instead use the free-space Green's kernel for the near field, add auxiliary basis functions for the far field, and enforce periodicity in an expanded linear system; this is robust for all parameters. Inverting the periodic and layer unknowns, we are left with a square linear system involving only the inclusion scattering coefficients. Preconditioning by the single-inclusion scattering matrix, this is solved iteratively in $O(M)$ time using a fast matrix-vector product. Numerical experiments show that a diffraction grating containing $M=1000$ inclusions per period can be solved to 9-digit accuracy in under 5 minutes on a laptop.