Autotuning Hamiltonian Monte Carlo for efficient generalized nullspace exploration

We propose methods to efficiently explore the generalized nullspace of (non-linear) inverse problems, defined as the set of plausible models that explain observations within some misfit tolerance. Owing to the random nature of observational errors, the generalized nullspace is an inherently probabilistic entity, described by a joint probability density of tolerance values and model parameters. Our exploration methods rest on the construction of artificial Hamiltonian systems, where models are treated as high-dimensional particles moving along a trajectory through model space. In the special case where the distribution of misfit tolerances is Gaussian, the methods are identical to standard Hamiltonian Monte Carlo, revealing that its apparently meaningless momentum variable plays the intuitive role of a directional tolerance. Its direction points from the current towards a new acceptable model, and its magnitude is the corresponding misfit increase. We address the fundamental problem of producing independent plausible models within a high-dimensional generalized nullspace by autotuning the mass matrix of the Hamiltonian system. The approach rests on a factorized and sequentially preconditioned version of the L-BFGS method, which produces local Hessian approximations for use as a near-optimal mass matrix. An adaptive time stepping algorithm for the numerical solution of Hamilton's equations ensures both stability and reasonable acceptance rates of the generalized nullspace sampler. In addition to the basic method, we propose variations of it, where autotuning focuses either on the diagonal elements of the mass matrix or on the macroscopic (long-range) properties of the generalized nullspace distribution. We quantify the performance of our methods in a series of numerical experiments, involving analytical, high-dimensional, multimodal test functions. These are designed to mimic realistic inverse problems, where sensitivity to different model parameters varies widely, and where parameters tend to be correlated. The tests indicate that the effective sample size may increase by orders of magnitude when autotuning is used. Finally, we present a proof of principle of generalized nullspace exploration in viscoelastic full-waveform inversion. In this context, we demonstrate (1) the quantification of inter- and intraparameter trade-offs, (2) the flexibility to change model parametrization a posteriori, for instance, to adapt averaging length scales, (3) the ability to perform dehomogenization to retrieve plausible subwavelength models and (4) the extraction of a manageable number of alternative models, potentially located in distinct local minima of the misfit functional.