Multi-index Sequential Monte Carlo Ratio Estimators for Bayesian Inverse problems

• Published in: Foundations of computational mathematics Vol. 24; no. 4; pp. 1249 - 1304
• Main Authors: Jasra, Ajay, Law, Kody J. H., Walton, Neil, Yang, Shangda
• Format: Journal Article
• Language: English
• Published: New York Springer US 01.08.2024; Springer Nature B.V
• Subjects: Applications of Mathematics; Bayesian analysis; Complexity; Computer Science; Context; Economics; Elliptic differential equations; Estimators; Fast Fourier transformations; Fourier transforms; Gaussian process; Inverse problems; Linear and Multilinear Algebras; Machine learning; Math Applications in Computer Science; Mathematics; Mathematics and Statistics; Matrix Theory; Normalizing (statistics); Numerical Analysis; Parabolic differential equations; Partial differential equations; Statistical inference; Statistical methods; 65C40; Bayesian inverse problems; Sequential Monte Carlo; 65C05; Multi-index Monte Carlo; 62F15; 62M30; 60J22
• ISSN: 1615-3375; 1615-3383
• DOI: 10.1007/s10208-023-09612-z

We consider the problem of estimating expectations with respect to a target distribution with an unknown normalising constant, and where even the un-normalised target needs to be approximated at finite resolution. This setting is ubiquitous across science and engineering applications, for example in the context of Bayesian inference where a physics-based model governed by an intractable partial differential equation (PDE) appears in the likelihood. A multi-index sequential Monte Carlo (MISMC) method is used to construct ratio estimators which provably enjoy the complexity improvements of multi-index Monte Carlo (MIMC) as well as the efficiency of sequential Monte Carlo (SMC) for inference. In particular, the proposed method provably achieves the canonical complexity of MSE - 1 , while single-level methods require MSE - ξ for ξ > 1 . This is illustrated on examples of Bayesian inverse problems with an elliptic PDE forward model in 1 and 2 spatial dimensions, where ξ = 5 / 4 and ξ = 3 / 2 , respectively. It is also illustrated on more challenging log-Gaussian process models, where single-level complexity is approximately ξ = 9 / 4 and multilevel Monte Carlo (or MIMC with an inappropriate index set) gives ξ = 5 / 4 + ω , for any ω > 0 , whereas our method is again canonical. We also provide novel theoretical verification of the product-form convergence results which MIMC requires for Gaussian processes built in spaces of mixed regularity defined in the spectral domain, which facilitates acceleration with fast Fourier transform methods via a cumulant embedding strategy, and may be of independent interest in the context of spatial statistics and machine learning.