Gideon Simpson

Frank J. Pinski, Gideon Simpson, Andrew M. Stuart et al.

In this paper we study algorithms to find a Gaussian approximation to a target measure defined on a Hilbert space of functions; the target measure itself is defined via its density with respect to a reference Gaussian measure. We employ the Kullback-Leibler divergence as a distance and find the best Gaussian approximation by minimizing this distance. It then follows that the approximate Gaussian must be equivalent to the Gaussian reference measure, defining a natural function space setting for the underlying calculus of variations problem. We introduce a computational algorithm which is well-adapted to the required minimization, seeking to find the mean as a function, and parameterizing the covariance in two different ways: through low rank perturbations of the reference covariance; and through Schrödinger potential perturbations of the inverse reference covariance. Two applications are shown: to a nonlinear inverse problem in elliptic PDEs, and to a conditioned diffusion process. We also show how the Gaussian approximations we obtain may be used to produce improved pCN-MCMC methods which are not only well-adapted to the high-dimensional setting, but also behave well with respect to small observational noise (resp. small temperatures) in the inverse problem (resp. conditioned diffusion).

Gideon Simpson, Marc Spiegelman

We present a model problem for benchmarking codes that investigate magma migration in the Earth's interior. This system retains the essential features of more sophisticated models, yet has the advantage of possessing solitary wave solutions. The existence of such exact solutions to the nonlinear problem make it an excellent benchmark problem for combinations of solver algorithms. In this work, we explore a novel algorithm for computing high quality approximations of the solitary waves and use them to benchmark a semi-Lagrangian Crank-Nicholson scheme for a finite element discretization of the time dependent problem.

Gideon Simpson, Mitchell Luskin

Parallel replica dynamics is a method for accelerating the computation of processes characterized by a sequence of infrequent events. In this work, the processes are governed by the overdamped Langevin equation. Such processes spend much of their time about the minima of the underlying potential, occasionally transitioning into different basins of attraction. The essential idea of parallel replica dynamics is that the exit time distribution from a given well for a single process can be approximated by the minimum of the exit time distributions of $N$ independent identical processes, each run for only 1/N-th the amount of time. While promising, this leads to a series of numerical analysis questions about the accuracy of the exit distributions. Building upon the recent work in Le Bris et al., we prove a unified error estimate on the exit distributions of the algorithm against an unaccelerated process. Furthermore, we study a dephasing mechanism, and prove that it will successfully complete.

Jeremy L. Marzuola, Gideon Simpson

In this work, we study the spectral properties of matrix Hamiltonians generated by linearizing the nonlinear Schrödinger equation about soliton solutions. By a numerically assisted proof, we show that there are no embedded eigenvalues for the three dimensional cubic equation. Though we focus on a proof of the 3d cubic problem, this work presents a new algorithm for verifying certain spectral properties needed to study soliton stability. Source code for verification of our comptuations, and for further experimentation, are available at http://www.math.toronto.edu/simpson/files/spec_prop_code.tgz.

Brittan A Farmer, Mitchell Luskin, Petr Plecháč et al.

Metastable condensed matter typically fluctuates about local energy minima at the femtosecond time scale before transitioning between local minima after nanoseconds or microseconds. This vast scale separation limits the applicability of classical molecular dynamics methods and has spurned the development of a host of approximate algorithms. One recently proposed method is diffusive molecular dynamics which aims to integrate a system of ordinary differential equations describing the likelihood of occupancy by one of two species, in the case of a binary alloy, while quasistatically evolving the locations of the atoms. While diffusive molecular dynamics has shown to be efficient and provide agreement with observations, it is fundamentally a model, with unclear connections to classical molecular dynamics. In this work, we formulate a spin-diffusion stochastic process and show how it can be connected to diffusive molecular dynamics. The spin-diffusion model couples a classical overdamped Langevin equation to a kinetic Monte Carlo model for exchange amongst the species of a binary alloy. Under suitable assumptions and approximations, spin-diffusion can be shown to lead to diffusive molecular molecular dynamics type models. The key assumptions and approximations include a well defined time scale separation, a choice of spin exchange rates, a low temperature approximation, and a mean field type approximation. We derive several models from different assumptions and show their relationship to diffusive molecular dynamics. Differences and similarities amongst the models are explored in a simple test problem.

Daniel Ginsberg, Gideon Simpson

We consider an equation for a thin-film of fluid on a rotating cylinder and present several new analytical and numerical results on steady state solutions. First, we provide an elementary proof that both weak and classical steady states must be strictly positive so long as the speed of rotation is nonzero. Next, we formulate an iterative spectral algorithm for computing these steady states. Finally, we explore a non-existence inequality for steady state solutions from the recent work of Chugunova, Pugh, & Taranets.

Gideon Simpson, Daniel Watkins

One way of getting insight into non-Gaussian measures, posed on infinite dimensional Hilbert spaces, is to first obtain best fit Gaussian approximations, which are more amenable to numerical approximation. These Gaussians can then be used to accelerate sampling algorithms. This begs the questions of how one should measure optimality and how the optimizers can be obtained. Here, we consider the problem of minimizing the distance with respect to relative entropy. We examine this minimization problem by seeking roots of the first variation of relative entropy, taken with respect to the mean of the Gaussian, leaving the covariance fixed. Adapting a convergence analysis of Robbins-Monro to the infinite dimensional setting, we can justify the application of this algorithm and highlight necessary assumptions to ensure convergence, not only in the context of relative entropy minimization, but other infinite dimensional problems as well. Numerical examples in path space, showing the robustness of this method with respect to dimension, are provided.

Can Huang, Michela Ottobre, Gideon Simpson

We study numerical schemes for Stochastic Partial Differential Equations (SPDEs). We introduce a general method of proof of non-asymptotic uniform in time error bounds on numerical integrators for SPDEs, ensuring the schemes capture both the transient and the long term dynamics faithfully. We then consider SPDEs with non-globally Lipshitz nonlinearities, which include for example the stochastic Allen-Cahn equation and some stochastic advection-diffusion equations. For the case of Allen-Cahn type SPDEs we show that the classic semi-implicit Euler time-discretization can exhibit finite time blow up. This motivates analysing other schemes which do not suffer from this blow-up problem. We consider three numerical schemes for SPDEs with non globally Lipshitz nonlinearity: a fully implicit scheme and two tamed schemes. For these schemes we prove non-asymptotic uniform in time error bounds by leveraging our general criterion, and provide numerical comparisons. While the main emphasis in this paper is on the properties of the time-discretization, the schemes we consider are full space-time discretization of the SPDE.

Rosina O. Weber, Prateek Goel, Shideh Amiri et al.

Previous research in interpretable machine learning (IML) and explainable artificial intelligence (XAI) can be broadly categorized as either focusing on seeking interpretability in the agent's model (i.e., IML) or focusing on the context of the user in addition to the model (i.e., XAI). The former can be categorized as feature or instance attribution. Example- or sample-based methods such as those using or inspired by case-based reasoning (CBR) rely on various approaches to select instances that are not necessarily attributing instances responsible for an agent's decision. Furthermore, existing approaches have focused on interpretability and explainability but fall short when it comes to accountability. Inspired in case-based reasoning principles, this paper introduces a pseudo-metric we call Longitudinal distance and its use to attribute instances to a neural network agent's decision that can be potentially used to build accountable CBR agents.

Aquil D. Jones, Gideon Simpson, William Wilson

Weak turbulence is a phenomenon by which a system generically transfers energy from low to high wave numbers, while persisting for all finite time. It has been conjectured by Bourgain that the 2D defocusing nonlinear Schrödinger equation (NLS) on the torus has this dynamic, and several analytical and numerical studies have worked towards addressing this point. In the process of studying the conjecture, Colliander, Keel, Staffilani, Takaoka, and Tao introduced a "toy model" dynamical system as an approximation of NLS, which has been subsequently studied numerically. In this work, we formulate and examine several numerical schemes for integrating this model equation. The model has two invariants, and our schemes aim to conserve at least one of them. We prove convergence in some cases, and our numerical studies show that the schemes compare favorably to others, such as Trapezoidal Rule and fixed step fourth order Runge-Kutta. The preservation of the invariants is particularly important in the study of weak turbulence as the energy transfer tends to occur on long time scales.

Gideon Simpson, Mitchell Luskin, David J. Srolovitz

Diffusive molecular dynamics is a novel model for materials with atomistic resolution that can reach diffusive time scales. The main ideas of diffusive molecular dynamics are to first minimize an approximate variational Gaussian free energy of the system with respect to the mean atomic coordinates (averaging over many vibrational periods), and to then to perform a diffusive step where atoms and vacancies (or two species in a binary alloy) flow on a diffusive time scale via a master equation. We present a mathematical framework for studying this algorithm based upon relative entropy, or Kullback-Leibler divergence. This adds flexibility in how the algorithm is implemented and interpreted. We then compare our formulation, relying on relative entropy and absolute continuity of measures, to existing formulations. The main difference amongst the equations appears in a model for vacancy diffusion, where additional entropic terms appear in our development.

Derek Olson, Soumitra Shukla, Gideon Simpson et al.

We examine the Petviashvilli method for solving the equation $ ϕ- Δϕ= |ϕ|^{p-1} ϕ$ on a bounded domain $Ω\subset \mathbb{R}^d$ with Dirichlet boundary conditions. We prove a local convergence result, using spectral analysis, akin to the result for the problem on $\mathbb{R}$ by Pelinovsky & Stepanyants, 2004. We also prove a global convergence result by generating a suite of nonlinear inequalities for the iteration sequence, and we show that the sequence has a natural energy that decreases along the sequence.

Andrew Binder, Tony Lelièvre, Gideon Simpson

Metastability is a common obstacle to performing long molecular dynamics simulations. Many numerical methods have been proposed to overcome it. One method is parallel replica dynamics, which relies on the rapid convergence of the underlying stochastic process to a quasi-stationary distribution. Two requirements for applying parallel replica dynamics are knowledge of the time scale on which the process converges to the quasi-stationary distribution and a mechanism for generating samples from this distribution. By combining a Fleming-Viot particle system with convergence diagnostics to simultaneously identify when the process converges while also generating samples, we can address both points. This variation on the algorithm is illustrated with various numerical examples, including those with entropic barriers and the 2D Lennard-Jones cluster of seven atoms.