Approximation and sampling of multivariate probability distributions in the tensor train decomposition

Abstract

General multivariate distributions are notoriously expensive to sample from, particularly the high-dimensional posterior distributions in PDE-constrained inverse problems. This paper develops a sampler for arbitrary continuous multivariate distributions that is based on low-rank surrogates in the tensor-train format, a methodology that has been exploited for many years for scalable, high-dimensional density function approximation in quantum physics and chemistry. We build upon recent developments of the cross approximation algorithms in linear algebra to construct a tensor-train approximation to the target probability density function using a small number of function evaluations. For sufficiently smooth distributions the storage required for accurate tensor-train approximations is moderate, scaling linearly with dimension. In turn, the structure of the tensor-train surrogate allows sampling by an efficient conditional distribution method since marginal distributions are computable with linear complexity in dimension. Expected values of non-smooth quantities of interest, with respect to the surrogate distribution, can be estimated using transformed independent uniformly-random seeds that provide Monte Carlo quadrature, or transformed points from a quasi-Monte Carlo lattice to give more efficient quasi-Monte Carlo quadrature. Unbiased estimates may be calculated by correcting the transformed random seeds using a Metropolis-Hastings accept/reject step, while the quasi-Monte Carlo quadrature may be corrected either by a control-variate strategy, or by importance weighting. We show that the error in the tensor-train approximation propagates linearly into the Metropolis-Hastings rejection rate and the integrated autocorrelation time of the resulting Markov chain; thus the integrated autocorrelation time may be made arbitrarily close to 1, implying that, asymptotic in sample size, the cost per effectively independent sample is one target density evaluation plus the cheap tensor-train surrogate proposal that has linear cost with dimension. These methods are demonstrated in three computed examples: fitting failure time of shock absorbers; a PDE-constrained inverse diffusion problem; and sampling from the Rosenbrock distribution. The delayed rejection adaptive Metropolis (DRAM) algorithm is used as a benchmark. In all computed examples, the importance-weight corrected quasi-Monte Carlo quadrature performs best, and is more efficient than DRAM by orders of magnitude across a wide range of approximation accuracies and sample sizes. Indeed, all the methods developed here significantly outperform DRAM in all computed examples.

Solutions of the Multivariate Inverse Frobenius–Perron Problem

“…To overcome this cost, Dolgov et al [ 25 ] precomputed an approximation of in a compressed tensor train representation that allows the fast computation of integrals in ( 17 ) and subsequent simulation of the inverse Rosenblatt transformation from the conditionals in ( 17 ) and showed that computational cost scales linearly with dimension d . Practical examples presented in [ 25 ], in dimension , demonstrate that operation by the forward and inverse Rosenblatt transformations is computationally feasible for multivariate problems with no special structure.…”

“…These induce pseudo-random and quasi-Monte Carlo sequences, respectively, on the space X via the inverse Rosenblatt transformation [ 29 ]. Both these schemes were demonstrated in practical high-dimensional settings in [ 25 ].…”

“…The RHS of Equation ( 20 ) in dimension is exactly the standard computational route for implementing inverse cumulative transformation sampling from , since computational uniform pseudo-random number generators perform a deterministic iteration on to implement a uniform map [ 29 ]. For , the RHS of Equation ( 20 ) is the conditional distribution method that generalizes the inverse cumulative transformation method, as mentioned above [ 22 , 23 , 24 , 25 ]. Thus, Lemma 1 shows that the standard computational implementation of both the inverse cumulative transformation method in and the conditional distribution method in is equivalent to implementing a solution to the IFFP.…”

“…Then, if , it follows that ; i.e., x is distributed as the desired target distribution [ 22 ]. This is the basis of the conditional distribution method for generating multi-variate random variables, which generalizes the inverse cumulative transformation method for univariate distributions [ 22 , 23 , 24 , 25 ]. These results may also be established by substituting R or into the the FP Equation ( 2 ), noting that there is a single inverse image and that the Jacobian determinant of R equals the target PDF .…”

Deep Composition of Tensor-Trains Using Squared Inverse Rosenblatt Transports

“…The recently developed transport map idea, e.g. [4,17,39,43,50], offers new insights for this task by identifying a measurable mapping, T : U → X , such that the pushforward of μ, denoted by T μ, is a close approximation to ν π . Then, the mapping T can be used to either accelerate classical sampling methods such as MCMC or to improve the efficiency of importance sampling.…”