Publications

We consider when a sparse nonnegative matrix S can be recovered, via an elementwise nonlinearity, from a real-valued matrix L of significantly lower rank. Of particular interest is the setting where the positive elements of S encode the similarities of nearby points on a low dimensional manifold. The recovery can then be posed as a problem in manifold learning— in this case, how to learn a norm-preserving and neighborhood-preserving mapping of high dimensional inputs into a lower dimensional space. We describe an algorithm for this problem based on a generalized low-rank decomposition of sparse matrices. This decomposition has the interesting property that it can be encoded by a neural network with one layer of rectified linear units; since the algorithm discovers this encoding, it can be viewed as a layerwise primitive for deep learning. The algorithm regards the inputs xi and xj as similar whenever the cosine of the angle between them exceeds some threshold τ ∈ (0,1). Given this threshold, the algorithm attempts to discover a mapping xi → yi by matching the elements of two sparse matrices; in particular, it seeks a mapping for which S = max(0,L), where Sij = max(0,xi·xj−τ∥xi∥∥xj∥) and Lij = yi·yj−τ∥yi∥∥yj∥. We apply the algorithm to data sets where vector magnitudes and small cosine distances have interpretable meanings (e.g., the brightness of an image, the similarity to other words). On these data sets, the algorithm discovers much lower dimensional representations that preserve these meanings.

We consider the analysis of singular waveguides separating insulating phases in two-space dimensions. The insulating domains are modeled by a massive Schrödinger equation and the singular waveguide by appropriate jump conditions along the one-dimensional interface separating the insulators. We present an integral formulation of the problem and analyze its mathematical properties. We also implement a fast multipole and sweeping-accelerated iterative algorithm for solving the integral equations, and demonstrate numerically the fast convergence of this method. Several numerical examples of solutions and scattering effects illustrate our theory.

We consider the inverse acoustic obstacle problem for sound-soft star-shaped obstacles in two dimensions wherein the boundary of the obstacle is determined from measurements of the scattered field at a collection of receivers outside the object. One of the standard approaches for solving this problem is to reformulate it as an optimization problem: finding the boundary of the domain that minimizes the $L^2$ distance between computed values of the scattered field and the given measurement data. The optimization problem is computationally challenging since the local set of convexity shrinks with increasing frequency and results in an increasing number of local minima in the vicinity of the true solution. In many practical experimental settings, low frequency measurements are unavailable due to limitations of the experimental setup or the sensors used for measurement. Thus, obtaining a good initial guess for the optimization problem plays a vital role in this environment.
We present a neural network warm-start approach for solving the inverse scattering problem, where an initial guess for the optimization problem is obtained using a trained neural network. We demonstrate the effectiveness of our method with several numerical examples. For high frequency problems, this approach outperforms traditional iterative methods such as Gauss-Newton initialized without any prior (i.e., initialized using a unit circle), or initialized using the solution of a direct method such as the linear sampling method. The algorithm remains robust to noise in the scattered field measurements and also converges to the true solution for limited aperture data. However, the number of training samples required to train the neural network scales exponentially in frequency and the complexity of the obstacles considered. We conclude with a discussion of this phenomenon and potential directions for future research.

We introduce an efficient numerical method for second order linear ODEs whose solution may vary between highly oscillatory and slowly changing over the solution interval. In oscillatory regions the solution is generated via a nonoscillatory phase function that obeys the nonlinear Riccati equation. We propose a defect-correction iteration that gives an asymptotic series for such a phase function; this is numerically approximated on a Chebyshev grid with a small number of nodes. For analytic coefficients we prove that each iteration, up to a certain maximum number, reduces the residual by a factor of order of the local frequency. The algorithm adapts both the step size and the choice of method, switching to a conventional spectral collocation method away from oscillatory regions. In numerical experiments we find that our proposal outperforms other state-of-the-art oscillatory solvers, most significantly at low-to-intermediate frequencies and at low tolerances, where it may use up to 106 times fewer function evaluations. Even in high frequency regimes, our implementation is on average 10 times faster than other specialized solvers.

We develop efficient and accurate sum-of-exponential (SOE) approximations for the Gaussian using rational approximation of the exponential function on the negative real axis. Six digit accuracy can be obtained with eight terms and ten digit accuracy can be obtained with twelve terms. This representation is of potential interest in approximation theory but we focus here on its use in accelerating the fast Gauss transform (FGT) in one and two dimensions. The one-dimensional scheme is particularly straightforward and easy to implement, requiring only twenty-four lines of MATLAB code. The two-dimensional version requires some care with data structures, but is significantly more efficient than existing FGTs. Following a detailed presentation of the theoretical foundations, we demonstrate the performance of the fast transforms with several numerical experiments.

We discuss the challenges of motivating, constructing, and quantising a canonically-normalised inflationary perturbation in spatially curved universes. We show that this has historically proved challenging due to the interaction of non-adiabaticity with spatial curvature. We propose a novel curvature perturbation which is canonically normalised, unique up to a single scalar parameter. This corrected quantisation has potentially observational consequences via modifications to the primordial power spectrum at large angular scales, as well as theoretical implications for quantisation procedures in curved cosmologies filled with a scalar field.

We describe a new, adaptive solver for the two-dimensional Poisson equation in complicated geometries. Using classical potential theory, we represent the solution as the sum of a volume potential and a double layer potential. Rather than evaluating the volume potential over the given domain, we first extend the source data to a geometrically simpler region with high order accuracy. This allows us to accelerate the evaluation of the volume potential using an efficient, geometry-unaware fast multipole-based algorithm. To impose the desired boundary condition, it remains only to solve the Laplace equation with suitably modified boundary data. This is accomplished with existing fast and accurate boundary integral methods. The novelty of our solver is the scheme used for creating the source extension, assuming it is provided on an adaptive quad-tree. For leaf boxes intersected by the boundary, we construct a universal "stencil" and require that the data be provided at the subset of those points that lie within the domain interior. This universality permits us to precompute and store an interpolation matrix which is used to extrapolate the source data to an extended set of leaf nodes with full tensor-product grids on each. We demonstrate the method's speed, robustness and high-order convergence with several examples, including domains with piecewise smooth boundaries.

Rapid advances in deep learning have brought not only myriad powerful neural networks, but also breakthroughs that benefit established scientific research. In particular, automatic differentiation (AD) tools and computational accelerators like GPUs have facilitated forward modeling of the Universe with differentiable simulations. Current differentiable cosmological simulations are limited by memory, thus are subject to a trade-off between time and space/mass resolution. They typically integrate for only tens of time steps, unlike the standard non-differentiable simulations. We present a new approach free of such constraints, using the adjoint method and reverse time integration. It enables larger and more accurate forward modeling, and will improve gradient based optimization and inference. We implement it in a particle-mesh (PM) N-body library pmwd (particle-mesh with derivatives). Based on the powerful AD system JAX, pmwd is fully differentiable, and is highly performant on GPUs.

The formation of the large-scale structure, the evolution and distribution of galaxies, quasars, and dark matter on cosmological scales, requires numerical simulations. Differentiable simulations provide gradients of the cosmological parameters, that can accelerate the extraction of physical information from statistical analyses of observational data. The deep learning revolution has brought not only myriad powerful neural networks, but also breakthroughs including automatic differentiation (AD) tools and computational accelerators like GPUs, facilitating forward modeling of the Universe with differentiable simulations. Because AD needs to save the whole forward evolution history to backpropagate gradients, current differentiable cosmological simulations are limited by memory. Using the adjoint method, with reverse time integration to reconstruct the evolution history, we develop a differentiable cosmological particle-mesh (PM) simulation library pmwd (particle-mesh with derivatives) with a low memory cost. Based on the powerful AD library JAX, pmwd is fully differentiable, and is highly performant on GPUs.

Recent developments in Markov chain Monte Carlo (MCMC) algorithms now allow us to run thousands of chains in parallel almost as quickly as a single chain, using hardware accelerators such as GPUs. We explore the benefits of running many chains, with an emphasis on achieving a target precision in as short a time as possible. One expected advantage is that, while each chain still needs to forget its initial point during a warmup phase, the subsequent sampling phase can be almost arbitrarily short. To determine if the resulting short chains are reliable, we need to assess how close the Markov chains are to convergence to their stationary distribution. The R̂ statistic is a general purpose and popular convergence diagnostic but unfortunately can require a long sampling phase to work well. We present a nested design to overcome this challenge and a generalization called nested R̂ . This new diagnostic works under conditions similar to R̂ and completes the MCMC workflow for many GPU-friendly samplers. In addition, the proposed nesting provides theoretical insights into the utility of R̂ , in both classical and short-chains regimes.