# zbMATH — the first resource for mathematics

Approximation and sampling of multivariate probability distributions in the tensor train decomposition. (English) Zbl 1436.62192
Summary: 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.

##### MSC:
 62H10 Multivariate distribution of statistics 65C05 Monte Carlo methods
##### Software:
ALEA; GitHub; GMRFLib; Rtwalk; SPLIDA; TT Toolbox; t-walk
Full Text:
##### References:
 [1] Atchadé, Yf, An adaptive version for the metropolis adjusted langevin algorithm with a truncated drift, Methodol. Comput. Appl. Probab., 8, 2, 235-254 (2006) · Zbl 1104.65004 [2] Ballani, J.; Grasedyck, L., Hierarchical tensor approximation of output quantities of parameter-dependent PDEs, SIAM/ASA J. Uncertain. Quantif., 3, 1, 852-872 (2015) · Zbl 1327.65010 [3] Brooks, S.; Gelman, A.; Jones, G.; Meng, Xl, Handbook of Markov Chain Monte Carlo (2011), Boca Raton: CRC Press, Boca Raton [4] Christen, J.; Fox, C., A general purpose sampling algorithm for continuous distributions (the t-walk), Bayesian Anal., 5, 2, 263-282 (2010) · Zbl 1330.62007 [5] Devroye, L., Non-Uniform Random Variate Generation (1986), Berlin: Springer, Berlin · Zbl 0593.65005 [6] Dick, J.; Kuo, F.; Sloan, I., High-dimensional integration: the quasi-Monte Carlo way, Acta Numer., 22, 133-288 (2013) · Zbl 1296.65004 [7] Dodwell, T.; Ketelsen, C.; Scheichl, R.; Teckentrup, A., A hierarchical multilevel Markov chain Monte Carlo algorithm with applications to uncertainty quantification in subsurface flow, SIAM/ASA J. Uncertain. Quantif., 3, 1, 1075-1108 (2015) · Zbl 1330.65007 [8] Dolgov, S.; Savostyanov, D., Alternating minimal energy methods for linear systems in higher dimensions, SIAM J. Sci. Comput., 36, 5, A2248-A2271 (2014) · Zbl 1307.65035 [9] Dolgov, S.; Scheichl, R., A hybrid alternating least squares—TT cross algorithm for parametric PDEs, SIAM/ASA J. Uncertain. Quantif., 7, 1, 260-291 (2019) · Zbl 1418.65170 [10] Eigel, M.; Gittelson, C.; Schwab, C.; Zander, E., Adaptive stochastic Galerkin FEM, Comput. Methods Appl., 270, 247-269 (2014) · Zbl 1296.65157 [11] Eigel, M.; Marschall, M.; Schneider, R., Sampling-free Bayesian inversion with adaptive hierarchical tensor representations, Inverse Probl., 34, 3, 035010 (2018) · Zbl 1404.65261 [12] Fox, C., Nicholls, G.: Sampling conductivity images via MCMC. In: The Art and Science of Bayesian Image Analysis, Leeds Annual Statistics Research Workshop, pp. 91-100 (1997) [13] Fox, C.; Norton, R., Fast sampling in a linear-Gaussian inverse problem, SIAM/ASA J. Uncertain. Quantif., 4, 1, 1191-1218 (2016) · Zbl 1398.94081 [14] Fox, C.; Parker, A., Accelerated Gibbs sampling of normal distributions using matrix splittings and polynomials, Bernoulli, 23, 4, 3711-3743 (2017) · Zbl 1457.62085 [15] Fox, C.; Haario, H.; Christen, J.; Damien, P.; Dellaportas, P.; Polson, N.; Stephens, D., Inverse problems, Bayesian Theory and Applications, 619-643 (2013), Oxford: Oxford University Press, Oxford [16] Gilks, W.; Wild, P., Adaptive rejection sampling for Gibbs sampling, Appl. Stat., 41, 337-348 (1992) · Zbl 0825.62407 [17] Golub, Gh; Van Loan, Cf, Matrix Computations (2013), Baltimore: Johns Hopkins University Press, Baltimore [18] Goreinov, S.; Tyrtyshnikov, E.; Zamarashkin, N., A theory of pseudoskeleton approximations, Linear Algebra Appl., 261, 1-3, 1-21 (1997) · Zbl 0877.65021 [19] Goreinov, S.; Oseledets, I.; Savostyanov, D.; Tyrtyshnikov, E.; Zamarashkin, N.; Olshevsky, V.; Tyrtyshnikov, E., How to find a good submatrix, Matrix Methods: Theory, Algorithms, Applications, 247-256 (2010), Singapore: World Scientific, Singapore · Zbl 1215.65078 [20] Gutierrez-Pulido, H.; Aguirre-Torres, V.; Christen, J., A practical method for obtaining prior distributions in reliability, IEEE Trans. Reliab., 54, 2, 262-269 (2005) [21] Haario, H.; Laine, M.; Mira, A.; Saksman, E., DRAM: efficient adaptive MCMC, Stat. Comput., 16, 4, 339-354 (2006) [22] Häggström, O.; Rosenthal, J., On variance conditions for Markov chain CLTs, Electron. Commun. Probab., 12, 454-464 (2007) · Zbl 1191.60082 [23] Hoang, Vh; Schwab, C.; Stuart, Am, Complexity analysis of accelerated MCMC methods for Bayesian inversion, Inverse Probl., 29, 8, 085010 (2013) · Zbl 1288.65004 [24] Hörmann, W.; Leydold, J.; Derflinger, G., Automatic Nonuniform Random Variate Generation (2004), Berlin: Springer, Berlin · Zbl 1038.65002 [25] Johnson, M., Multivariate Statistical Simulation (1987), New York: Wiley, New York · Zbl 0604.62056 [26] Khoromskij, B., Structured rank-$$(r_1,\ldots , r_d)$$ decomposition of function-related operators in $${\mathbb{R}}^d$$, Comput. Methods Appl. Math., 6, 2, 194-220 (2006) · Zbl 1120.65052 [27] Kuo, F.; Scheichl, R.; Schwab, C.; Sloan, I.; Ullmann, E., Multilevel quasi-Monte Carlo methods for lognormal diffusion problems, Math. Comput., 86, 2827-2860 (2017) · Zbl 1368.65005 [28] Liu, J., Metropolized independent sampling with comparisons to rejection sampling and importance sampling, Stat. Comput., 6, 2, 113-119 (1996) [29] Martino, L.; Read, J.; Luengo, D., Independent doubly adaptive rejection Metropolis sampling within Gibbs sampling, IEEE Trans. Signal Process., 63, 12, 3123-3138 (2015) · Zbl 1394.94828 [30] Meeker, W.; Escobar, L., Statsitical Methods for Reliability Data (1998), New York: Wiley, New York [31] Mengersen, Kl; Tweedie, Rl, Rates of convergence of the Hastings and Metropolis algorithms, Ann. Stat., 24, 1, 101-121 (1996) · Zbl 0854.60065 [32] Meyer, R.; Cai, B.; Perron, F., Adaptive rejection Metropolis sampling using Lagrange interpolation polynomials of degree 2, Comput. Stat. Data Anal., 52, 7, 3408-3423 (2008) · Zbl 1452.62099 [33] Mira, A., Ordering and improving the performance of Monte Carlo Markov chains, Stat. Sci., 16, 4, 340-350 (2001) · Zbl 1127.60312 [34] Mira, A., Geyer, C.J.: Ordering Monte Carlo Markov chains. Tech. Rep. 632, Univ. of Minnesota (1999) [35] Niederreiter, H., Quasi-Monte Carlo methods and pseudo-random numbers, Bull. Am. Math. Soc., 84, 6, 957-1041 (1978) · Zbl 0404.65003 [36] Norton, R.; Christen, J.; Fox, C., Sampling hyperparameters in hierarchical models: improving on Gibbs for high-dimensional latent fields and large datasets, Commun. Stat. Simul., 47, 9, 2639-2655 (2018) [37] O’Connor, P.; Kleyner, A., Practical Reliability Engineering (2012), New York: Wiley, New York [38] Oseledets, I., DMRG approach to fast linear algebra in the TT-format, Comput. Methods Appl. Math., 11, 3, 382-393 (2011) · Zbl 1283.15041 [39] Oseledets, I., Tensor-train decomposition, SIAM J. Sci. Comput., 33, 5, 2295-2317 (2011) · Zbl 1232.15018 [40] Oseledets, I., Constructive representation of functions in low-rank tensor formats, Constr. Approx., 37, 1, 1-18 (2013) · Zbl 1282.15021 [41] Oseledets, I.; Tyrtyshnikov, E., TT-cross approximation for multidimensional arrays, Linear Algebra Appl., 432, 1, 70-88 (2010) · Zbl 1183.65040 [42] Oseledets, I., Dolgov, S., Kazeev, V., Savostyanov, D., Lebedeva, O., Zhlobich, P., Mach, T., Song, L.: TT-Toolbox (2011). https://github.com/oseledets/TT-Toolbox [43] Roberts, Go; Rosenthal, Js, Quantitative non-geometric convergence bounds for independence samplers, Methodol. Comput. Appl. Probab., 13, 2, 391-403 (2011) · Zbl 1222.60054 [44] Rosenblatt, M., Remarks on a multivariate transformation, Ann. Math. Stat., 23, 3, 470-472 (1952) · Zbl 0047.13104 [45] Rue, H., Fast sampling of Gaussian Markov random fields, J. R. Stat. Soc. B, 63, 325-338 (2001) · Zbl 0979.62075 [46] Rue, H.; Held, L., Gaussian Markov Random Fields: Theory and Applications (2005), London: Chapman & Hall, London · Zbl 1093.60003 [47] Scheichl, R.; Stuart, A.; Teckentrup, A., Quasi-Monte Carlo and multilevel Monte Carlo methods for computing posterior expectations in elliptic inverse problems, SIAM/ASA J. Uncertain. Quantif., 5, 1, 493-518 (2017) · Zbl 06736512 [48] Schneider, R.; Uschmajew, A., Approximation rates for the hierarchical tensor format in periodic Sobolev spaces, J. Complex., 30, 2, 56-71 (2013) · Zbl 1329.41033 [49] Smith, R.L., Tierney, L.: Exact transition probabilities for the independence Metropolis sampler. Tech. rep., Univ. of North Carolina (1996) [50] Stuart, A., Inverse problems: a Bayesian perspective, Acta Numer., 19, 451-559 (2010) · Zbl 1242.65142 [51] Teckentrup, A.; Scheichl, R.; Giles, M.; Ullmann, E., Further analysis of multilevel Monte Carlo methods for elliptic PDEs with random coefficients, Numer. Math., 125, 3, 569-600 (2013) · Zbl 1306.65009 [52] Tierney, L., A note on Metropolis-Hastings kernels for general state spaces, Ann. Appl. Probab., 8, 1, 1-9 (1998) · Zbl 0935.60053 [53] Tyrtyshnikov, E., Tensor approximations of matrices generated by asymptotically smooth functions, Sbornik Math, 194, 6, 941-954 (2003) · Zbl 1067.65044 [54] Wolff, U., Monte Carlo errors with less errors, Comput. Phys. Commun., 156, 2, 143-153 (2004) · Zbl 1196.65138
This reference list is based on information provided by the publisher or from digital mathematics libraries. Its items are heuristically matched to zbMATH identifiers and may contain data conversion errors. It attempts to reflect the references listed in the original paper as accurately as possible without claiming the completeness or perfect precision of the matching.