Tony Shardlow

Grigoris Katsiolides, Eike H. Müller, Robert Scheichl et al.

A common way to simulate the transport and spread of pollutants in the atmosphere is via stochastic Lagrangian dispersion models. Mathematically, these models describe turbulent transport processes with stochastic differential equations (SDEs). The computational bottleneck is the Monte Carlo algorithm, which simulates the motion of a large number of model particles in a turbulent velocity field; for each particle, a trajectory is calculated with a numerical timestepping method. Choosing an efficient numerical method is particularly important in operational emergency-response applications, such as tracking radioactive clouds from nuclear accidents or predicting the impact of volcanic ash clouds on international aviation, where accurate and timely predictions are essential. In this paper, we investigate the application of the Multilevel Monte Carlo (MLMC) method to simulate the propagation of particles in a representative one-dimensional dispersion scenario in the atmospheric boundary layer. MLMC can be shown to result in asymptotically superior computational complexity and reduced computational cost when compared to the Standard Monte Carlo (StMC) method, which is currently used in atmospheric dispersion modelling. To reduce the absolute cost of the method also in the non-asymptotic regime, it is equally important to choose the best possible numerical timestepping method on each level. To investigate this, we also compare the standard symplectic Euler method, which is used in many operational models, with two improved timestepping algorithms based on SDE splitting methods.

Peter Kloeden, Tony Shardlow

Mean-field SDEs, also known as McKean-Vlasov equations, are stochastic differential equations where the drift and diffusion depend on the current distribution in addition to the current position. We describe an efficient numerical method for approximating the distribution at time t of the solution to the initial-value problem for one-dimensional mean-field SDEs. The idea is to time march (e.g., using the Euler-Maruyama time-stepping method) an m-point Gauss quadrature rule. With suitable regularity conditions, convergence with first order is proved for Euler-Maruyama time stepping. We also estimate the work needed to achieve a given accuracy in terms of the smoothness of the underlying problem. Numerical experiments are given, which show the effectiveness of this method as well as two second-order time-stepping methods. The methods are also effective for ordinary SDEs in one dimension, as we demonstrate by comparison with the multilevel Monte Carlo method.

Ivan Cheltsov, Federico Cornalba, Clarice Poon et al.

Priors with non-smooth log-densities, such as the l1-prior, are widely used in Bayesian inverse problems for their sparsity-inducing properties. Existing Langevin-based sampling methods typically rely on proximal mappings or smooth approximations, which alter the target distribution. We propose an alternative approach based on a Hadamard product parameterization of the l1-norm, leading to a smooth but nonconvex and non-globally Lipschitz potential whose marginal law exactly recovers the desired posterior. The resulting Hadamard Langevin dynamics (HLD) defines a diffusion process that is analytically distinct from proximal or mirror-type Langevin schemes. Our main contribution is a rigorous well-posedness theory for both the continuous and discrete HLD. We establish existence and uniqueness of strong solutions, geometric ergodicity of the continuous dynamics, and convergence of the discretized scheme as the step size tends to zero. These results provide the first theoretical foundation for sampling from nonconvex, nonsmooth posteriors through overparameterized Langevin dynamics.

Stephen Marsland, Tony Shardlow

Registration of images parameterised by landmarks provides a useful method of describing shape variations by computing the minimum-energy time-dependent deformation field that flows one landmark set to the other. This is sometimes known as the geodesic interpolating spline and can be solved via a Hamiltonian boundary-value problem to give a diffeomorphic registration between images. However, small changes in the positions of the landmarks can produce large changes in the resulting diffeomorphism. We formulate a Langevin equation for looking at small random perturbations of this registration. The Langevin equation and three computationally convenient approximations are introduced and used as prior distributions. A Bayesian framework is then used to compute a posterior distribution for the registration, and also to formulate an average of multiple sets of landmarks.

Tony Shardlow

The Feynman-Kac formula for the exterior-value problem for the fractional Laplacian leads to a walk-outside-spheres algorithm via sampling alpha-stable Levy processes on their exit from maximally inscribed balls and sampling their occupation distribution. Kyprianou, Osojnik, and Shardlow (2017) developed this algorithm, providing a complexity analysis and an implementation, for approximating the solution at a single point in the domain. This paper shows how to efficiently sample the whole field by generating an approximation in L_2(D), for a domain D . The method takes advantage of a hierarchy of triangular meshes and uses the multilevel Monte Carlo method for Hilbert space-valued quantities of interest. We derive complexity bounds in terms of the fractional parameter alpha and demonstrate that the method gives accurate results for two problems with exact solutions. Finally, we show how to couple the method with the variable-accuracy Arnoldi iteration to compute the smallest eigenvalue of the fractional Laplacian. A criteria is derived for the variable accuracy and a comparison is given with analytical results of Dyda (2012).

Teo Deveney, Eike Mueller, Tony Shardlow

We introduce a deep learning accelerated methodology to solve PDE-based Bayesian inverse problems with guaranteed accuracy. This is motivated by the ill-posed problem of inferring a spatio-temporal heat-flux parameter known as the Biot number given temperature data, however the methodology is generalisable to other settings. To accelerate Bayesian inference, we develop a novel training scheme that uses data to adaptively train a neural-network surrogate simulating the parametric forward model. By simultaneously identifying an approximate posterior distribution over the Biot number, and weighting a physics-informed training loss according to this, our approach approximates forward and inverse solution together without any need for external solves. Using a random Chebyshev series, we outline how to approximate a Gaussian process prior, and using the surrogate we apply Hamiltonian Monte Carlo (HMC) to sample from the posterior distribution. We derive convergence of the surrogate posterior to the true posterior distribution in the Hellinger metric as our adaptive loss approaches zero. Additionally, we describe how this surrogate-accelerated HMC approach can be combined with traditional PDE solvers in a delayed-acceptance scheme to a-priori control the posterior accuracy. This overcomes a major limitation of deep learning-based surrogate approaches, which do not achieve guaranteed accuracy a-priori due to their non-convex training. Biot number calculations are involved in turbo-machinery design, which is safety critical and highly regulated, therefore it is important that our results have such mathematical guarantees. Our approach achieves fast mixing in high dimensions whilst retaining the convergence guarantees of a traditional PDE solver, and without the burden of evaluating this solver for proposals that are likely to be rejected. Numerical results are given using real and simulated data.

Tony Shardlow

We study the sticky Cox-Ingersoll-Ross (CIR) process in one dimension, a diffusion on $[0,\infty)$ with a sticky boundary condition at the origin, arising as the marginal process in a sparse Bayesian inference framework based on Hadamard-Langevin dynamics. For the parameter range $δ\in(1,2)$, in which the origin is accessible but not absorbing, we prove well-posedness of the process and uniqueness of its invariant measure, which is a mixture of a point mass at zero and a weighted gamma-type density on the interior. We derive an explicit Green's function for the resolvent in terms of confluent hypergeometric functions, and use this to construct an exact sampler for the invariant measure in the zero-potential case. For a non-trivial potential $G$, we establish existence and uniqueness of the tilted invariant measure via a Girsanov change of measure, and develop two sampling algorithms: a Metropolis-Hastings corrected sampler that targets the invariant measure exactly, and an unadjusted Langevin algorithm (ULA) that is cheaper per step but introduces an $O(h)$ bias. Numerical experiments confirm the predicted behaviour: the Metropolis-Hastings sampler achieves the target invariant measure at all step sizes, while the ULA exhibits the expected $O(h)$ bias.

Allen Paul, Neill Campbell, Tony Shardlow

We derive an algorithm for compression of the currents and varifolds representations of shapes, using ridge leverage score (RLS) sampling, and the theory of Nystrom approximation in Reproducing Kernel Hilbert Spaces. Our method is faster than existing compression techniques and comes with theoretical guarantees on the rate of decay of the compression error as a function of the smoothness of the associated shape representation. The obtained compressions are shown to be useful for accelerating downstream tasks such as nonlinear shape registration in the Large Deformation Diffeomorphic Metric Mapping (LDDMM) framework without loss of quality, even for very high compression ratios. The performance of our algorithm is demonstrated on large-scale shape data from modern geometry processing datasets, and is shown to be fast and scalable with rapid error decay.

Fengpei Wang, Clarice Poon, Tony Shardlow

The use of optimal transport (OT) distances, and in particular entropic-regularised OT distances, is an increasingly popular evaluation metric in many areas of machine learning and data science. Their use has largely been driven by the availability of efficient algorithms such as the Sinkhorn algorithm. One of the drawbacks of the Sinkhorn algorithm for large-scale data processing is that it is a two-phase method, where one first draws a large stream of data from the probability distributions, before applying the Sinkhorn algorithm to the discrete probability measures. More recently, there have been several works developing stochastic versions of Sinkhorn that directly handle continuous streams of data. In this work, we revisit the recently introduced online Sinkhorn algorithm of [Mensch and Peyré, 2020]. Our contributions are twofold: We improve the convergence analysis for the online Sinkhorn algorithm, the new rate that we obtain is faster than the previous rate under certain parameter choices. We also present numerical results to verify the sharpness of our result. Secondly, we propose the compressed online Sinkhorn algorithm which combines measure compression techniques with the online Sinkhorn algorithm. We provide numerical experiments to show practical numerical gains, as well as theoretical guarantees on the efficiency of our approach.

Allen Paul, Neill Campbell, Tony Shardlow

We extend our work for compression of currents and varifolds to a compression algorithm for the embedded normal cycles representation of shape, restricted to the constant normal kernel case, using the Nystrom approximation in Reproducing Kernel Hilbert Spaces (RKHS) and ridge leverage score (RLS) sampling. Our method comes with theoretical guarantees on the compression error decay, and the approximations are shown to be effective for downstream tasks such as nonlinear shape registration in the Large Deformation Diffeomorphic Metric Mapping (LDDMM) framework, even for very high compression ratios. The performance of our algorithm is demonstrated on large-scale shape data from modern geometry processing datasets and is shown to accelerate downstream registration tasks significantly.

Adam Hartshorne, Allen Paul, Tony Shardlow et al.

This work presents a unified framework for the unsupervised prediction of physically plausible interpolations between two 3D articulated shapes and the automatic estimation of dense correspondence between them. Interpolation is modelled as a diffeomorphic transformation using a smooth, time-varying flow field governed by Neural Ordinary Differential Equations (ODEs). This ensures topological consistency and non-intersecting trajectories while accommodating hard constraints, such as volume preservation, and soft constraints, \eg physical priors. Correspondence is recovered using an efficient Varifold formulation, that is effective on high-fidelity surfaces with differing parameterisations. By providing a simple skeleton for the source shape only, we impose physically motivated constraints on the deformation field and resolve symmetric ambiguities. This is achieved without relying on skinning weights or any prior knowledge of the skeleton's target pose configuration. Qualitative and quantitative results demonstrate competitive or superior performance over existing state-of-the-art approaches in both shape correspondence and interpolation tasks across standard datasets.

Allen Paul, George Grammatopoulos, Adwaye Rambojun et al.

Dysplasia is a recognised risk factor for osteoarthritis (OA) of the hip, early diagnosis of dysplasia is important to provide opportunities for surgical interventions aimed at reducing the risk of hip OA. We have developed a pipeline for semi-automated classification of dysplasia using volumetric CT scans of patients' hips and a minimal set of clinically annotated landmarks, combining the framework of the Gaussian Process Latent Variable Model with diffeomorphism to create a statistical shape model, which we termed the Gaussian Process Diffeomorphic Statistical Shape Model (GPDSSM). We used 192 CT scans, 100 for model training and 92 for testing. The GPDSSM effectively distinguishes dysplastic samples from controls while also highlighting regions of the underlying surface that show dysplastic variations. As well as improving classification accuracy compared to angle-based methods (AUC 96.2% vs 91.2%), the GPDSSM can save time for clinicians by removing the need to manually measure angles and interpreting 2D scans for possible markers of dysplasia.

Teo Deveney, Eike Mueller, Tony Shardlow

We investigate a deep learning approach to efficiently perform Bayesian inference in partial differential equation (PDE) and integral equation models over potentially high-dimensional parameter spaces. The contributions of this paper are two-fold; the first is the introduction of a neural network approach to approximating the solutions of Fredholm and Volterra integral equations of the first and second kind. The second is the development of a new, efficient deep learning-based method for Bayesian inversion applied to problems that can be described by PDEs or integral equations. To achieve this we introduce a surrogate model, and demonstrate how this allows efficient sampling from a Bayesian posterior distribution in which the likelihood depends on the solutions of PDEs or integral equations. Our method relies on the direct approximation of parametric solutions by neural networks, without need of traditional numerical solves. This deep learning approach allows the accurate and efficient approximation of parametric solutions in significantly higher dimensions than is possible using classical discretisation schemes. Since the approximated solutions can be cheaply evaluated, the solutions of Bayesian inverse problems over large parameter spaces are efficient using Markov chain Monte Carlo. We demonstrate the performance of our method using two real-world examples; these include Bayesian inference in the PDE and integral equation case for an example from electrochemistry, and Bayesian inference of a function-valued heat-transfer parameter with applications in aviation.

Andreas E. Kyprianou, Ana Osojnik, Tony Shardlow

We consider Monte Carlo methods for simulating solutions to the analogue of the Dirichlet boundary-value problem in which the Laplacian is replaced by the fractional Laplacian and boundary conditions are replaced by conditions on the exterior of the domain. Specifically, we consider the analogue of the so-called `walk-on-spheres` algorithm. In the diffusive setting, this entails sampling the path of Brownian motion as it uniformly exits a sequence of spheres maximally inscribed in the domain. As this algorithm would otherwise never end, it is truncated when the `walk-on-spheres` comes within epsilon > 0 of the boundary. In the setting of the fractional Laplacian, the role of Brownian motion is replaced by an isotropic alpha-stable process with alpha in (0, 2). A significant difference to the Brownian setting is that the stable processes will exit spheres by a jump rather than hitting their boundary. This difference ensures that disconnected domains may be considered and that, unlike the diffusive setting, the algorithm ends after an almost surely finite number of steps.