Low-rank tensor reconstruction of concentrated densities with application to Bayesian inversion

Transport maps have become a popular mechanic to express complicated probability densities using sample propagation through an optimized push-forward. Beside their broad applicability and well-known success, transport maps suffer from several drawbacks such as numerical inaccuracies induced by the optimization process and the fact that sampling schemes have to be employed when quantities of interest, e.g. moments are to compute. This paper presents a novel method for the accurate functional approximation of probability density functions (PDF) that copes with those issues. By interpreting the pull-back result of a target PDF through an inexact transport map as a perturbed reference density, a subsequent functional representation in a more accessible format allows for efficient and more accurate computation of the desired quantities. We introduce a layer-based approximation of the perturbed reference density in an appropriate coordinate system to split the high-dimensional representation problem into a set of independent approximations for which separately chosen orthonormal basis functions are available. This effectively motivates the notion of h- and p-refinement (i.e. ``mesh size'' and polynomial degree) for the approximation of high-dimensional PDFs. To circumvent the curse of dimensionality and enable sampling-free access to certain quantities of interest, a low-rank reconstruction in the tensor train format is employed via the Variational Monte Carlo method. An a priori convergence analysis of the developed approach is derived in terms of Hellinger distance and the Kullback-Leibler divergence. Applications comprising Bayesian inverse problems and several degrees of concentrated densities illuminate the (superior) convergence in comparison to Monte Carlo and Markov-Chain Monte Carlo methods.

In this work, we propose a new approach based on function space representations with efficient surrogate models in several instances.This is motivated by our previous work on adaptive low-rank approximations of solutions of parametric random PDEs with Adaptive Stochastic Galerkin FEM (ASGFEM, see e.g.Eigel et al. 2017Eigel et al. , 2020) ) and in particular the sampling-free Bayesian inversion presented in Eigel et al. (2018) where the setting of uniform random variables was examined.A generalization to the important case of Gaussian random variables turns out to be non-trivial from a computational point of view due to the difficulties caused by the representation of highly concentrated densities in a compressing tensor format which is required to cope with the high dimensionality of the problem.As a consequence, we develop a discretization approach which takes into account the potentially difficult to represent structure of the probability density at hand by a combination of several transformations and approximations that can be chosen adaptively to counteract the interplay of the employed numerical approximations.With the computed functional representation of the density, the evaluation of moments and statistical quantities of interest can be carried out efficiently and with high accuracy.Additionally, we point out that the generated surrgate can be used for the fast generation of samples from the posterior distribution.
The central idea of the method is to construct a map which transports the target density to some convenient reference density and to employ low-rank regression techniques to obtain a functional representation, for which accurate numerical methods are available.Transport maps for probability densities are a classical topic in mathematics.They are under active research in particular in the area of optimal transport (Villani 2008;Santambrogio 2015) and also have become popular in current machine learning research (Tran et al. 2019;Rezende and Mohamed 2015;Detommaso et al. 2019).A main application we have in mind is Bayesian inversion where, given a prior density and some observations of the forward model, a posterior density should be determined.In this context, the rescaling approaches in Schillings and Schwab (2016) and Schillings et al. (2020) based on the Laplace approximation can be considered as transport maps of a certain (affine) form.More general transport maps have been examined extensively in El Moselhy and Marzouk (2012) and Parno and Marzouk (2018) and other works of the research group.Obtaining a transport map is in general realized by minimizing a certain loss functional, e.g. the Kullback-Leibler distance, between the target and the push-forward of a reference density.This process has been analyzed and improved using iterative maps (Brennan et al. 2020) or multi-scale approaches (Parno et al. 2016).However, the optimization, the loss functional and the chosen model class for the transport map yield only an approximation to an exact transport.We hence suppose that, in general, only an inexact transport is available.Considering the pullback of the target, this can be interpreted as starting from a slightly or severly perturbed reference density.One then has to cope with the degree of the perturbation in subsequent approximation steps to enable an accurate explicit representation of this new reference density.
Finding a suitable approximation relies on concepts from adaptive finite element methods (FEM).In addition to the selection of (local) approximation spaces of a certain degree (in the spirit of "p-refinement"), we introduce a spatial decomposition of the density representation into layers (similar to "h-refinement") around some center of mass of the considered density.This enables to exploit the (assumed) decay behavior of the approximated density.Overall, this "hp-refinement" allows to balance inaccuracies and hence to compensate perturbations of the reference density by putting more computational effort into the discretization part.Consequently, one enjoys the freedom to decide whether more effort should be invested into computing an exact transport map or into a more elaborate discretization of the perturbed reference density.
For eventual computations with the devised (possibly high-dimensional) functional density representation, an efficient representation format is required.In our context, hierarchical tensors and in particular tensor trains (TT) prove to be advantageous, cf.Bachmayr et al. (2016) and Oseledets (2011).These low-rank formats enable to alleviate the curse of dimensionality under suitable conditions and allow for efficient evaluations of high-dimensional objects.For each layer of the discretization we aim to obtain a low-rank tensor representation of the respective perturbed reference density.In certain ideal cases such as transporting to the standard Gaussian density a rank-one representation is sufficient.A discussion on the low-rank approximation of perturbations of Gaussian densities is carried out in Rohrbach et al. (2020).In more general cases, a low-rank representability may be observed numerically.To allow for tensor methods to be applicable, the desired discretization layers have to be tensor domains.Therefore, the underlying perturbed reference density is transformed to an alternative coordinate system which benefits the representation and enables to exploit the regularity and decay behavior of the density.To generate a tensor train representation (coupled with a function basis which is then also called extended or functional TT format Gorodetsky et al. 2015), the Variational Monte Carlo (VMC) method (Eigel et al. 2019b) is employed.It basically is a tensor regression approach based on function samples for which a convergence analysis is available.Notably, depending on the chosen loss functional, it leads to the best approximation in the respective model space.It has previously been examined in the context of random PDEs in Eigel et al. (2019b) as an alternative nonintrusive numerical approach to Stochastic Galerkin FEM in the TT format (Eigel et al. 2017(Eigel et al. , 2020)).The approximation of Eigel et al. (2020) is used in one of the presented examples for Bayesian inversion with the random Darcy equation with lognormal coefficient.We note that surrogate models of the forward model have been used in the context of MCMC e.g. in Li and Marzouk (2014) and tensor representations (obtained by cross approximation) were used in Dolgov et al. (2020) to improve the efficiency of MCMC sampling.
The derivation of our method is supported by an a priori convergence analysis with respect to the Hellinger distance and the Kullback-Leibler divergence.In the analysis, different error sources have to be considered, in particular a layer truncation error depending on decay properties of the density, a low-rank truncation error and model space approximations are introduced.Moreover, the VMC error analysis (Eigel et al. 2019b) comprising statistical estimation and numerical approximation errors is adjusted to be applicable to the devised approach.While not usable for an a posteriori error control in its current initial form, the derived analysis leads the way to more elaborate results for this promising method in future research.
With the constructed functional density surrogate, samplingfree computations of statistical quantities of interest such as moments or marginals become feasible by fast tensor contractions.
While several assumptions have to be satisfied for this method to work most efficiently, the approach is rather general and can be further adapted to the problem at hand.Moreover, it should be emphasized that by constructing a functional representation, structural properties of the density at hand (in particular smoothness, sparsity, low-rank approximability and decay behavior in different parameters) can be exploited in a much more extensive way than what is possible with sampling based methods such as MCMC, leading to more accurate statistical computations and better convergence rates.We note that the perturbed posterior surrogate can be used to efficiently generate samples by rejection sampling or within an MCMC scheme.Since the perturbed transport can be seen as a preconditioner, the sample generation can be based on the perturbed prior.These samples can then be transported to the posterior by the determined push-forward.As a prospective extension, the constructed posterior density could directly be used in a Stochastic Galerkin FEM based on the integral structure, closing the loop of forward and inverse problem, resulting in the inferred forward problem with model data determined by Bayesian inversion from the observed data.
The structure of the paper is as follows.Section 2 is concerned with the representation of probability densities and introduces a relation between a target and a reference density.Such a transport map can be determined numerically by approximation in a chosen class of functions and with an assumed structure, leading to the concept of perturbed ref-erence densities.To counteract the perturbation, a layered truncated discretization is introduced.An efficient low-rank representation of the mappings is described in Sect. 3 where the tensor train format is discussed.In order to obtain this nonintrusively, the Variational Monte Carlo (VMC) tensor reconstruction is reviewed.A priori convergence results with respect to the Hellinger distance and Kullback-Leibler divergence are derived in Sect. 4. For practical purposes, the proposed method is described in terms of an algorithm in Sect. 5. Possible applications we have in mind are examined in Sect.6.In particular, the setting of Bayesian inverse problems is recalled.Moreover, the computation of moments and marginals is scrutinized.Section 7 illustrates the performance of the proposed method.In addition to an examination of the numerical sensitivity of the accuracy with respect to the perturbation of the transport maps, a typical model problem from Uncertainty Quantification (UQ) is depicted, namely the identification of a parametrization for the random Darcy equation with lognormal coefficient given as solution of a stochastic Galerkin FEM.

Density representation
The aim of this section is to introduce the central ideas of the proposed approximation of densities.For this task, two established concepts are reviewed, namely transport maps (El Moselhy and Marzouk 2012;Marzouk et al. 2016), which are closely related to the notion of optimal transport (Villani 2008;Santambrogio 2015), and hierarchical low-rank tensor representations (Oseledets 2011;Hackbusch 2012;Bachmayr et al. 2016).By the combination of these techniques, assuming the access to a suitable transformation, the developed approach yields a functional representation of the density in a format which is suited for computations with high-dimensional functions.In particular, it becomes feasible to accurately handle highly concentrated posterior densities.While transport maps on their own in principle enable the generation of samples of some target distribution, the combination with a functional low-rank representation allows for integral quantities such as (centered) moments to become computable cheaply.Given an approximate transport map, the low-rank representation can be seen as a further approximation step (compensating for the inaccuracy of the used transport) to gain direct access to the target density.
Consider a target measure π with Radon-Nikodym derivative with respect to the Lebesgue measure λ denoted as In the following we assume that point evaluations of f are available up to a multiplicative constant, motivated by the framework of Bayesian posterior density representation with unknown normalization constant.Furthermore, let π 0 be some reference measure exhibiting a Radon-Nikodym derivative with respect to to the Lebesgue measure denoted as f 0 .This is motivated by the prior measure and density in the context of Bayesian inference.
In the upcoming sections we relate π and π 0 with the help of a transport map T .If the exact T is replaced by some T , in Sect.2.2 we alternatively can describe π with the help of a socalled associated auxiliary measure π0 with auxiliary density f0 .As in the Bayesian context, we also may refer to it as a perturbed prior measure with perturbed prior density.This concept introduces a possible workload balancing between an approximation of T and the associated perturbed prior.

Transport maps
The notion of density transport is classical and has become a popular research area with optimal transport, see e.g.Villani (2008) and Santambrogio (2015).From a practical point of view, it has been employed to improve numerical approaches for Bayesian inverse problems for instance in El Moselhy and Marzouk (2012), Brennan et al. (2020) and Dolgov et al. (2020).Similar approaches are discussed in terms of sample transport e.g. for Stein's method (Liu and Wang 2016;Detommaso et al. 2018) or multi-layer maps (Brennan et al. 2020).We review the properties required for our approach in what follows.Note that since our target application is Bayesian inversion, we usually use the terms "prior" and "posterior" instead of the more general "reference" and "target" densities.
Let X := R d and assume that there exists an exact transport map T : X → Y , which is a diffeomorphism1 relating π and π 0 via (2) This change of variables formula allows computations to be carried out in terms of the measure π 0 , which is commonly assumed to be of a simpler structure.In particular, for any function (3) Note that the computation of the right-hand side in (3) may still be a challenging task depending on the actual structure of In "Appendix A" we list several examples of transport maps with an exploitable structure.

Inexact transport and the perturbed prior
In general, the transport map T is unknown or difficult to determine and hence has to be approximated by some T : X → Y , e.g. using a polynomial representation with respect to π 0 (El Moselhy and Marzouk 2012) or with a composition of simple maps in a reduced space such as in Brennan et al. (2020).As a consequence of the approximation, it holds subject to the accuracy of the involved approximation of T .
The approximate transport T relates a measure π0 with density f0 to the target measure π , whereas We henceforth refer to (5) as the auxiliary reference or perturbed prior density.Using this construction, the moment computation reads If one would know f0 , by ( 5) and ( 6) one would also have access to the exact posterior.Equation ( 6) is the starting point of the proposed method by approximating f0 in another coordinate system which is better adapted to the structure of the approximate (perturbed) prior.For this, consider a (fixed) diffeomorphism with Jacobian x → |J Φ ( x)| and define the transformed perturbed prior If ( 8) can be approximated accurately by some function with accuracy determined only by the approximation quality of f h 0 .Thus, (9) enables to shift the complexity of a transport map approximation T ≈ T to a functional approximation f h 0 ≈ f0 in a new coordinate system X and allows for a balancing of both effects.The construction of T and a suitable map in (7) may be used to obtain a convenient transformed perturbed prior given in ( 8).An approximation thereof can be significantly simpler compared to a possibly complicated target density f or the computation of the exact transport T .We give some example scenarios and more details on the choice of Φ in "Appendix B".

Layer based representation
To further refine and motivate the notion of an adapted coordinate system, let L ∈ N and (X ) L =1 be pairwise disjoint domains in X such that is simply connected and compact and define X L+1 := X \ K .Then, for given L ∈ N we may decompose the perturbed prior f0 as where χ denotes the indicator function on X .Moreover, for any tensor set X := × d i=1 X i and diffeomorphism Φ : X →X , 1 ≤ ≤ L +1, we may represent the localized perturbed prior f 0 as a pull-back function where f 0 is a map defined on X as in (8).This layer based coordinate change enables a representation of the density on bounded domains.Even though the remainder layer is unbounded, we assume that K is sufficiently large to cover all probability mass of f0 except for a negligible higher-order error.
Up to this point, the choice of transformation Φ , = 1, . . ., L +1, is fairly general.However, for the further development of the method we assume the following property.
where denotes the Hadamard product of vectors.
Due to the notion of rank 1 stable transformations, the map x → T (Φ( x)) in (9) inherits the rank structure of T , see Sect. 3. Furthermore, since the Jacobian x → |J Φ ( x)| is rank 1, we can construct tensorized orthonormal basis functions which may be used to approximate the transformed perturbed prior in (8).
Remark 1 The described concept can be extended to any rank r ∈ N Jacobian of Φ, i.e.
Motivated by the right-hand side in (9), one may use different approximations of the perturbed transformed prior f0 • Φ in r distinct tensorized spaces, each associated to the rank 1

Layer truncation
This paragraph is devoted to the treatment of the last (remainder or truncation) layer X L+1 introduced in (11) with the aim to suggest some approximation choices.If f0 is represented in the layer format (11), it is convenient to simply extend the function to zero after layer L ∈ N. By this, the remaining (possibly small) probability mass is neglected.Such a procedure is typically employed in numerical applications and does not impose any computational issues since events on the outer truncated domain are usually exponentially unlikely for a truncation value chosen sufficiently large.Nevertheless, in order to present a rigorous treatment, we require properties like absolute continuity, which would be lost by using a cut-off function.Inspired by Schillings et al. (2020) regarding the information limit of unimodal posterior densities,2 we suggest a Gaussian approximation for the last layer L + 1 on the unbounded domain X L+1 , i.e. for some s.p.d.Σ ∈ R d,d and μ ∈ R d we define the density with and f Σ,μ denotes the Gaussian probability density function with mean μ and covariance matrix Σ.
Remark 2 A particular choice for μ and Σ would be the mean and covariance of the normalized version of f0 | K or the MAP and the corresponding square root of the Hessian of f0 .
Note that the constant C < L in ( 16) may exhibit an analytic form whereas computing C > L suffers from the curse of dimensionality and is in general not available.To circumvent this issue and render further use of the representation (15) feasible, we now aim for an approximation model to adequately represent the localized transformed perturbed prior maps f 0 = f 0 • Φ from ( 12).

Approximation model
The computation of high-dimensional integrals and the efficient construction of surrogates is a challenging task with a multitude of approaches.Some of these techniques are sparse grid methods (Chen and Schwab 2016;Garcke and Griebel 2012), collocation (Ernst et al. 2019;Nobile et al. 2008;Foo and Karniadakis 2010) or modern sampling techniques (Gilks et al. 1995;Rudolf and Sprungk 2017;Neal 2001).
A promising class of approximation model relies on the concept of compression in terms of low-rank formats (Grasedyck et al. 2013).Those can be seen as generalization of the singular value decomposition of functions with high-dimensional input.We give a short introduction of the low-rank format considered in Sect.3.1 and present a numerical scheme to construct such an approximation model by the Variational Monte Carlo (VMC) method in Sect.3.2.In order to enhance readability, in the first two sections we use the notation ĝ to underline the connection to the transformed perturbed prior f0 and its localizations, which is part of the total approximation model in Sect.3.3.

Low-rank tensor train format
and consider a map ĝ : X → R. The function ĝ can be represented in the TT format if there exists a rank vector r = (r 1 , . . ., In the forthcoming sections we consider weighted tensorized Lebesgue spaces in which the perturbed prior segments are defined.In particular, for a non-negative weight function This space may be identified by its tensorization We assume that there exists a complete orthonormal basis which is known a priori.For discretization purposes, we introduce the finite dimensional subspaces for i = 1, . . ., d, and n i ∈ N. On these we formulate the extended tensor train format in terms of the coefficient tensors such that every univariate function ĝi ∈ V i,n i can be written as Any function can be expressed in the full tensor format by a high dimensional algebraic tensor G : In contrast to this, the format given by ( 18) and ( 22) admits a linear structure in the dimension.More precisely, the memory complexity of O(max{n 1 , . . ., This observation raises the question of expressibility for certain classes of functions and the existence of a representation rank r where max{r 1 , . . ., r d−1 } stays sufficiently small for practical computations.This issue is e.g.addressed in Schneider and Uschmajew (2014), Bachmayr et al. (2017) and Griebel and Harbrecht (2013) under certain assumptions on the regularity and in Espig et al. (2009), Oseledets (2011), Ballani et al. (2013) and Eigel et al. (2019b) explicit (algorithmic) constructions of the format are discussed even in case that ĝ has no analytic representation.
For later reference we define the finite dimensional lowrank manifold of rank r tensor trains by This is an embedded manifold in the finite full tensor space V Λ from (23).We also require the concept of the algebraic (full) tensor space and the corresponding low-rank form for given r ∈ N d−1 defined by In the following, a method is reviewed that can be used to construct surrogates in the presented extended tensor train format using only samples of the sought high-dimensional function.

Tensor train regression by Variational Monte Carlo
We review the sampling-based VMC method presented in Eigel et al. (2019b) which is employed here to construct TT approximations of the transformed local perturbed priors f 0 as defined in ( 12).The approach generalizes the concept of randomized tensor completion (Eigel et al. 2019a).Its analysis relies on the theory of statistical learning, leading to a priori convergence results.It can also be seen as a generalized tensor least squares technique.In principle, the algorithm for the construction of the surrogate is interchangeable and e.g. an alternative cross-interpolation method for probability densities is presented in Dolgov et al. (2020).
For the VMC framework, consider the model class This space is equipped with the metric Moreover, note that due to the orthonormality of and ι( ĝ, •) is integrable with respect to the weight function w of V( X ) for every ĝ ∈ V( X , c, c) and the cost functional Then, the objective of the method is to find a minimizer Due to the infinite dimensional setting, we confine the minimization problem in (33) to our model class M.This yields the minimization problem Subsequently, instead of J a tractable empirical functional is considered, namely with independent samples { xk } k≤N distributed according to the measure wλ with a (possibly rescaled) weight function w with respect to the Lebesgue measure λ.The corresponding empirical optimization problem reads In our application the loss functional relates to the Kullback-Leibler divergence or the L 2 error for a given function ĝ and target f0 .In order to emphasize this dependence, we may write ι( ĝ, x; f0 ).

The total approximation model
We employ the VMC approach from Sect.3.2 to build an approximation of with samples { xk } N k=1 drawn from the (possibly rescaled) finite measure w λ.Finally, for a given transport T and a choice of transformations Φ : X → X , according to (15) our approximate of the perturbed prior f0 denoted by f Trun,TT 0 is defined by 16) and We refer to Fig. 1 for a visual presentation of the involved objects, approximations and transformations.

Error estimates
We now discuss the accuracy of the actual approximation of the target density f given by Since our approach is based on several components like transport, truncation, low-rank compression and the VMC method, these components are examined separately in the upcoming sections.Our main result is stated in Sect.4.5.

Transport invariant measures of discrepancy
In this section we derive a relation property such that the error of the approximation f Trun,TT 0 to the the perturbed prior transfers directly to the discrepancy between f TT and f .Note that this property is canonical since passing to the image space of some measurable function is fundamental in probability theory.Ideally such a relation is an equivalence of the form Prominent measures of discrepancy for two absolutely continuous Lebesgue probability density functions h 1 and h 2 on some measurable space Z are the squared Hellinger distance and the Kullback-Leibler divergence For the Hellinger distance, the absolute continuity assumption can be dropped from an analytical point of view.We observe that both d Hell and d KL satisfy (41).

Truncation error
Since our approximation scheme relies on the truncation of the density, we introduce a convenient type of decay on the outer layer of the perturbed prior.
Definition 2 (outer polynomial exponential decay) A function f0 : X → R + has outer polynomial exponential decay if there exists a simply connected compact set K ⊂ X with a polynomial π + which is positive on X \ K and some C > 0 such that The error introduced by the Gaussian extension is estimated in the next lemma.
with the incomplete Gamma function Γ .

Proof
The proof follows immediately from the definition of f Trun 0 .
In the case that the perturbed prior is close to a Gaussian standard normal distribution, it holds C ≈ 1.

Low-rank compression error
In this section we discuss the error introduced by compressing a full algebraic tensor into a tensor in low-rank tensor train format.The higher order singular value decomposition (HOSVD) (Oseledets and Tyrtyshnikov 2010) is based on successive unfoldings of the full tensor into matrices, which are orthogonalized and possibly truncated by a singular value decomposition.This algorithm leads to the following result.
Lemma 3 (Theorem 2.2 Oseledets and Tyrtyshnikov 2010) For any g ∈ V Λ and r ∈ R d−1 there exists an extended low-rank tensor train g r ∈ M r such that where σ i is the distance of the i-th unfolding matrix of the coefficient tensor of g in the HOSVD to its best rank r i approximation in the Frobenius norm.
Remark 3 Estimate ( 48) is rather unspecific as the σ i cannot be quantified a priori.In the special case of Gaussian densities we refer to Rohrbach et al. (2020) for an examination of the low-rank representation depending on the covariance structure.When the transport T maps the considered density only "closely" to a standard Gaussian, the results can be applied immediately to our setting and more precise estimates are possible.

VMC error analyis
To examine the VMC convergence in our setting, we recall the analysis of Eigel et al. (2019b) in a slightly more general manner.Analogously to Sect.3.1, we use the notation X and w as a placeholder for any layer X and weight w for = 1, . . ., L.Here we assume that L 2 ( X , w) is continuously embedded in L 1 ( X , w).
Recall the cost functional J from (32) defined by a loss function ι depending on the transformed perturbed prior f0 as in Sect.3.3.As a first step we show compatibility conditions of two specific types of loss functions corresponding to the Kullback-Leibler divergence and the L 2 -norm.
is uniformly bounded and Lipschitz continuous on the model The global Lipschitz continuity of J follows by using ( 50) and with a constant C related to the embedding of ) with coefficient tensors G 1 and G 2 ∈ TT r then by Parseval's identity and the finite dimensionality of M r ( X , c, c) there exists c = which yields the Lipschitz continuity on M r ( X , c, c).
is uniformly bounded and Lipschitz continuous on Due to c < ∞, the Lipschitz property follows as in the proof of Lemma 4 if ĝ1 , ĝ2 ∈ M r ( X , c, c).
Let ĝM and ĝM,N be as in ( 34) and ( 36).The analysis examines different error components with respect to f0 ∈ V( X , 0, c * ) for some 0 < c * < ∞ defined by denoting the VMC, approximation and generalization error, respectively.By a simple splitting, the VMC error can be bounded by the approximation and the generalization error, 3 namely Due to the global Lipschitz property on V( X , c, c) with c > 0 in the setting of ( 49) or c ≥ 0 as in ( 53), the approximation error can be bounded by the best approximation in M. In particular there exists C > 0 such that We note that such a bound by the best approximation in M with respect to the V( X )-norm may not be required when using the Kullback-Leibler divergence if one is interested directly in the best approximation in this divergence.Then the assumption c > 0 can be relaxed in the construction of V( X , c, c) since no global Lipschitz continuity of J in Lemma 4 is necessary.Thus the more natural subspace of V( X , 0, c) of absolutely continuous functions with respect to f0 may be considered instead.
It remains to bound the statistical generalization error E gen .For this the notion of covering numbers is required.Let (Ω, F, P) be an abstract probability space.
Remark 4 (choice of c, c and X ) Due to the layer based representation in ( 11) and (15) on each layer X = Φ −1 (X ) we have the freedom to choose c separately.In particular, assuming that the perturbed prior f0 decays per layer, we can choose c according to the decay and with this control the constant in (50).
3 Neglecting an intractable optimization error.

A priori estimate
In this section we state our main convergence result.
Assumption 1 For a target density f : Y → R + and a transport map T : X → Y , there exists a simply connected compact domain has outer polynomial exponential decay with polynomial π + on X \ K .Consider the symmetric positive definite matrix Σ ∈ R d,d and μ ∈ R d as the covariance and mean for the outer approximation f Σ,μ .Furthermore, let K = L =1 X where X is the image of a rank-1 stable diffeomorphism Φ : X → X such that there exists 0 < c * < ∞ with f 0 ( x) ≤ c * for x ∈ X for = 1, . . ., L.
We can now formulate the main theorem of this section regarding the convergence of the developed approximation with respect to the Hellinger distance and the KL divergence.
Theorem 1 (A priori convergence) Let Assumption 1 hold and let a sequence of sample sizes (N ) L =1 ⊂ N be given.For every = 1, . . ., L, consider bounds 0 < c < c < ∞ and let f TT be defined as in (40).Then there exist constants C, C Σ , C , C ι > 0, = 1, . . ., L, such that for ∈ {KL, Hell} and p Hell = 2 and p KL = 1 (61) Here, E best denotes the error of the best approximation ĝ , * Λ to f 0 in the full truncated space E sing is the low-rank approximation error of the algebraic tensor G : Λ → R associated to ĝ , * Λ and the truncation error E # trun is given by Furthermore, for any ( ) L =1 ⊂ R + the generalization errors E gen can be bounded in probability by Since K = ∪ L =1 X and X are bounded, there exist constants C(X ) > 0, = 1, . . ., L, such that Moreover, by construction The claim follows by application of Lemmas 2, 3 and 6 together with (58).
To show (61) for = KL, note that by Lemma 1 and the construction (38) it holds Using Lemma 2 we can bound the integral over X \ K by the truncation error E trun .Employing the loss function and cost functional of Lemma 4 yields The claim follows by application of Lemmas 3 and 6 together with (58).

Polynomial approximation in weighted L 2 spaces
In order to make the error bound (61) in Theorem 1 more explicit with respect to E best , we consider the case of a smooth density function with analytic extension.The analysis follows the presentation in Babuška et al. (2010) and leads to exponential convergence rates by an iterative interpolation argument based on univariate best approximation bounds by interpolation.An analogous analysis for more general regularity classes is possible but not in the scope of this article.
For a Hilbert space H , a bounded set I ⊂ R and a function f ∈ C(I ; H ) ⊂ L 2 (I , w; H ) with weight w : I → R, let I n : C(I ; H ) → L 2 (I , w; H ) denote the continuous Lagrange interpolation operator.
Assume that f ∈ C(I ; H ) admits an analytic extension in the region of the complex plane Σ(I ; τ ) := {z ∈ C| dist(z, I ) ≤ τ } for some τ > 0.Then, referring to Babuška et al. (2010), Using an iterative argument over d dimensions, a convergence rate for the interpolation of f ∈ C( X ; R) ⊂ L 2 ( X , w; R) can be derived from the 1-dimensional convergence.More specifically, let I Λ : C( X ) → L 2 ( X , w) denote the continuous interpolation operator I Λ := I 1 n 1 • I 2:d n 2 :n d written as composition of a 1-dimensional and a d − 1-dimensional interpolation with continuous Xi , ⊗ d i=2 w i and I 2,...,d n 2 ,...,n d : The second term of the last bound is a d − 1-dimensional interpolation and can hence be bounded uniformly over x1 by a similar iterative argument.We summarize the convergence result for E best in the spirit of (Theorem 4.1 Babuška et al. 2010).

Remark 5
The interpolation argument on polynomial discrete spaces could be expanded to other orthonormal systems such as trigonometric polynomial, admitting well-known Lebesgue constants as in Da Fies and Vianello (2013).
Remark 6 Explicit best approximation bounds for appropriate smooth weights w as in the case of spherical coordinates can be obtained using partial integration techniques as in Mead and Delves (1973).There, the regularity class of f0 uses high-order weighted Sobolev spaces based on derivatives of w as in the case of classical polynomials.

Algorithm
Since a variety of techniques is employed in the proposed density discretization, this section provides an exemplary algorithmic workflow to illustrate the required steps in practical applications (see also Fig. 7 for a sketch of the components of the method).The general method to obtain a representation of the density (1) by its auxiliary reference (5) is summarized in Algorithm 1.Based on this, the computation of possible quantities of interest such as moments (6) or marginals are considered in Sects.6.2.1 and 6.2.2, respectively.In the following we briefly describe the involved algorithmic procedures.
Computing the transformation Obtaining a suitable transport map is a current research topic and examined e.g. in Papamakarios et al. ( 2021), Parno and Marzouk (2018), Tran et al. (2019) and Marzouk et al. (2016).In Sect.2.1, two naive options are introduced.In the numerical applications, we employ an affine transport and also illustrate the capabilities of a quadratic transport in a two-dimensional example.For the affine linear transport we utilize a semi-Newton optimizer to obtain the maximum value of f and an approximation of the Hessian at the optimal value, see Sect.A.1.For the construction of a quadratic transport we rely on the library TransportMaps (Baptista et al. 2015(Baptista et al. -2018)).The task to provide the (possibly inexact) transport map is summarized in the function In the following paragraphs we assume Φ to be the multivariate polar transformation as in "Appendix B.1", defined on the corresponding hyperspherical shells X .We refer to X 1 as the radial dimension and X i as the angular dimensions for 1 < i ≤ d.The computations on each shell X , = 1, . . ., L, are fully decoupled and suitable for parallelization.Note that the proposed method is easily adapted to other transformations Φ .

Generating an orthonormal basis
To obtain suitable finite dimensional subspaces, one has to introduce spanning sets that allow for an efficient computation of e.g.moments (3) and the optimization of the functional (32).Given a fixed dimension vector n ∈ N d for the current X , = 1, . . ., L, and introducing the weight w given by the jacobian of the chosen transformation Φ , the function can be split into three distinct algorithmic parts as follows.
-1st coordinate x1 : The computation of an orthonormal polynomial basis {P 1,α } α with respect to the weight w 1 ( x1 ) = xd−1 1 in the radial dimension by a stabilized Gram-Schmidt method.This is numerically unstable since the involved summations cause cancellation.As a remedy, we define arbitrary precision polynomials with a significant digit length τ mant to represent polynomial coefficients.By this, point evaluations of the orthonormal polynomials and computations of integrals of the form e.g.required for computing moments with polynomial transport, can be realized with high precision.The length τ mant is set to 100 in the numerical examples and the additional run-time is negligible as the respective calculations can be precomputed.-2nd coordinate x2 : Since X 2 = [0, 2π ] and to preserve periodicity, we employ trigonometric polynomials given by Note that here the weight function is constant, i.e. w 2 ( x2 ) ≡ 1, and the defined trigonometric polynomials are orthonormal in L 2 ( X 2 ).coordinate x3 , . . ., xd : On the remaining angular dimensions i = 3, . . ., d, we employ the usual Gram-Schmidt orthogonalization algorithm on [0, π] with weight function w i ( xi ) = sin i ( xi ), based on polynomials.
Fortunately, the basis for dimensions 1 < i ≤ d coincides on every layer = 1, . . ., L. It hence can be computed just once and passed to the individual process handling the current layer.Only the basis in the radial dimension needs to be adjusted to X .The parameter τ GS collects all tolerance parameters for the applied numerical quadrature and the significant digit length τ mant .

Generation of Samples
We utilize Monte Carlo samples to approximate the L 2 integral for the empirical minimization formulation (37).To generate such samples on X , we employ the inverse transform sampling based on the inverse cumulative distribution function of the normalized version of the weight function ω .The generation process of N ∈ N samples is denoted as the function

Reconstruction of a Tensor Train surrogate
The VMC tensor reconstruction of Sect. 3 is summarized in the function The tensor components F ,TT 0,i are associated with the corresponding basis P i to form a rank r extended tensor train as defined in ( 18) and ( 22).The additional parameter τ Recon collects all parameters that determine the VMC algorithm.
The method basically involves the optimization of a loss functional over the set of tensor trains with rank (at most) r .In the presented numerical computations we consider a mean-squared loss and the respective empirical approximation based on a current sample set S .The tensor optimization-based on a rank adaptive, alternating direction fitting (ADF) algorithm-is implemented in the xerus library (Huber andWolf 2014-2017) and wrapped in the ALEA framework (Eigel et al.).Additionally, the machine learning framework PyTorch (Paszke et al. 2017) can be utilized in ALEA to minimize the empirical cost functional from (35) by a wide range of state-of-the-art stochastic optimizers.The latter enable stochastic gradient methods to compute the tensor coefficients as known from machine learning applications.With this setting in mind, the actual meaning of the parameter τ Recon depends on the chosen optimizer.In this article we focus on the ADF implementation and initialize e.g. the starting rank, the number of iteration of the ADF and a target residual norm.

Applications
In the preceding sections the creation of surrogate models of quite generic probability density functions are developed.Based on this, in the following we focus on actual applications where such a representation is beneficial.We start with the framework of Bayesian inverse problems with target density (1) corresponding to the Lebesgue posterior density.Subsequently, we cover the computation of moments and marginals.

Bayesian inversion
This section is devoted to a brief review of the Bayesian paradigm.We assume that the reader is familiar with the concept of statistical inverse problems and hence focus on the general formalism and highlight the notation with the setup of Sect. 2 in mind.We closely follow the presentation in Eigel et al. (2018) and refer to Stuart (2010), Dashti and Stuart (2016) and Kaipio and Somersalo (2006) for a comprehensive overview.
Let Y and Y denote separable Hilbert spaces equipped with inner products •, • H for H ∈ {Y , Y}. Assume there exists a parameter to observation map G : Y → Y such that for some given observation δ ∈ Y the relation δ = G(y) + η holds for some y ∈ Y and noise η ∈ Y.For instance in the Darcy model we take y as parameter determining the permeability of some porous medium and let G describe the observation of water pressure at some location in the domain.By taking η as a random variable with law N (0, C 0 ) for some symmetric positive definite covariance operator C 0 on Y, the inference of y which explains the observation δ becomes a statistical inverse problem.Consequently, the quantities y and δ become random variables over a probability space (Ω, F, P) with values in Y and Y, respectively.In Stuart (2010) mild conditions on the parameter to observation map are derived to show a continuous version of Bayes formula which yields the existence and uniqueness of the Radon-Nikodym derivative of the (posterior) measure π δ of the conditional random variable y|δ with respect to a prior measure π 0 of y.More precisely, by assuming η to be Gaussian and independent with respect to y, both measures π 0 and π δ on Y are related by the Bayesian potential The posterior density is given by Algorithm 1 Tensor surrogate creation of perturbed prior

Input:
Lebesgue target density f : (22) sample size N ∈ N, = 1, . . ., L for level-wise reconstruction tensor train ranks (r 1 , . . ., r L ), r ∈ N d−1 , for = 1, . . ., L (18) Gram-Schmidt tolerance parameter τ GS tensor reconstruction parameter τ Recon Output: Level-wise low-rank approximation of perturbed prior • Generate samples with respect to the weight w • Equip tensor components with basis We assume it exists and denote the normalization constant Z := E π 0 exp (−Ψ (y, δ)) .Note that we interchangeably write y as an element of Y and the corresponding random variable with values in Y .
This simplified version of a Bayesian inverse problem can be cast into the framework of this manuscript by using the notation f for the posterior density (73) and f 0 = dπ 0 for the prior density.

Statistical quantities of interest
A common task is the efficient computation of the expectation of some quantity of interest (QoI) We discuss the special case of moment computations in Sect.6.2.1 and the basis representations of marginals in Sect.6.2.2.In those cases the structure of Q allows for direct computations of the integrals via tensor contractions.For more involved choices of the QoI, we suggest a universal sampling approach by repeated evaluation of the low-rank surrogate.More precisely, by application of the integral transformation we can approximate and replace the integrals over X by Monte Carlo estimates with samples according to the (normalized) weight |J Φ |.
Those samples can be obtained by uniform sampling on the tensor spaces X and the inverse transform approach as mentioned in the paragraph Generating Samples of Sect. 5.Alternatively, efficient MCMC sampling by marginalization can be employed (Weare 2007).

Moment computation
In this section we discuss the computation of for the presented layer-based format with low-rank tensor train approximations.In particular, we are interested in an efficient generation of the moment map Given some transport T : X → Y with an associated per- We fix 1 ≤ ≤ L and assume tensor spaces X , X such that a layer-based splitting can be employed to obtain integrals over X of the form Note that we neglect the remaining unbounded layer X L+1 since for moderate |α| and vol( L =1 X ) sufficiently large, the contribution to the considered moment does not have a significant influence on the overall approximation.Additionally, a rank-1 stable diffeomorphism Φ : X → X is assumed for which there exist univariate functions Φ , j : X j → X with Φ , j = (Φ i, j ) d i=1 and h j : X j → R for every j = 1, . . ., d, such that

Moments under affine transport
be a symmetric positive definite matrix and M = (M i ) d i=1 ∈ R d such that the considered transport map takes the form With the multinomial coefficient for j ∈ N, β ∈ N d 0 for j = |β| given by the computation of moments corresponds to the multinomial theorem as seen in the next lemma.
where the high-dimensional coefficient C H k is given by The statement follows by the multinomial theorem since Generalizing Lemma 8 to multiindices α ∈ N d 0 yields where Exploiting the layer-wise TT representation of f from (38) and using the rank-1 stable map (79), the high-dimensional integral over X simplifies to Note that the right-hand side is composed via decoupled one dimensional integrals only.We point out that while the structure is simplified, the definition of Φ j in (83) a priori results in several integrals (indexed by d k=1 β k ).These integrals, whose number depends on the cardinality of α, have to be computed.In several cases this simplifies further, e.g. when Φ transforms the spherical coordinate system to Cartesian coordinates.

Moment computation using spherical coordinates
In the special case that Φ is the multivariate polar transformation of "Appendix B.1", the number of distinct computation of integrals from (85) reduces significantly.Recall that x1 = ρ, x2:d = θ = (θ 0 , . . ., θ d−2 ) and let The exponential complexity due to the indexing by d k=1 β k reduces to linear complexity in |α|.More precisely, the amount of exponents in ( 86) -( 88) is linear in the dimensions since the sums only depend on |α|, leading to O(|α|d) different integrals that may be precomputed for each tuple (k i−1 , k i ).This exponential complexity in the rank vanishes in the presence of an approximation basis associated with each coordinate dimension as defined in Sect.3.

Computation of marginals
In probability theory and statistics, marginal distributions and especially marginal probability density functions provide insights into an underlying joint density by means of lower dimensional functions that can be visualized.The computa-tion of marginal densities is a frequent problem encountered e.g. in parameter estimation and when using sampling techniques since histograms and corner plots provide easy access to (in general high-dimensional) integral quantities.
In contrast to the Markov Chain Monte Carlo algorithm, the previously presented method of a layer based surrogate for the Lebesgue density function f : Y = R d → R allows for a functional representation and approximation of marginal densities without additional evaluations of f .For y ∈ Y and i = 1, . . ., d, define as the marginalized variable where the i-th component is left out and f (y −i , y i ) := f (y).Then, for given i = 1, . . ., d, the i-th marginal density reads Computing this high-dimensional integral by quadrature or sampling is usually infeasible and the transport map approach as given by (3) fails since the map T : X → Y cannot be used directly in (89).Alternatively, we can represent d f i : R → R in a given orthonormal basis {ϕ j } N ϕ j=1 and consider where β j , j = 1, . . ., N ϕ denotes the L 2 (R) projection coefficient With this, the marginalization can be carried out similarly to the computations in Sect.6.2.1 by replacing f with f T T .A convenient basis is given by monomials since (91) then simplifies to This is the moment corresponding to the multiindex α = Alternatively, indicator functions may be considered in the spirit of histograms.

Stochastic Galerkin finite element method
Consider the case when the target density f : Y → R + with Y = R d should be used in a stochastic Galerkin Finite Element method (SGFEM, see e.g.Eigel et al. 2017Eigel et al. , 2020)).Note that f and its approximation f TT in general is not of product type.If there exists c > 0 and a norm Y on Y such that Y e c y Y f TT (y)dλ(y) < ∞ then there exists a complete orthonormal basis of polynomials (Dunkl and Xu 2014).The moment computation from Sect.6.2.1 can then be used to obtain the polynomials numerically based on a Gram-Schmidt orthonormalization process.This approach is out of the scope of the present work and will be discussed elsewhere.

Samples from the posterior
In case one is interested to draw samples from f , we propose using standard sample generation schemes like rejection sampling or MCMC on the approximated perturbed prior f TT 0 .Such a scheme can be implemented efficiently based on the presented surrogate since fast point evaluations are possible.Once samples are obtained from f TT 0 , mapping them through T yields samples from f TT .Note that we cannot draw samples in the reference space X and map them via the transformations {Φ } to obtain samples from f TT 0 since the transformation is not a transport.

Numerical validation and applications
This section is devoted to a numerical validation of the proposed Algorithm 1 using various types of transformations T while applying the scheme to practical problems.We focus on three example settings.The first consists of an artificial Gaussian posterior density which could be translated to a linear forward model and Gaussian prior assumptions in the Bayesian setting.Second, we study the approximation under inexact transport and conclude as a third setting with an actual Bayesian inversion application governed by the lognormal Darcy flow problem.All considered examples satisfy Assumption 1 and are in fact analytic.This is obvious for the Gaussian and the transformed Gaussian in Sects.7.1 and 7.2.In the log-normal case we refer to Babuška et al. (2007) and Hoang and Schwab (2014).
Remark 7 When dealing with probability density functions, a key requirement for a valid approximation is non-negativity.Some comments are in order since this property cannot be guaranteed when applying usual concepts to obtain tensor train approximiations.Nevertheless, we would like to point out that we did not experience any numerical difficulties when computing the presented experiments and the following two approaches have not been implemented.
Note that the chosen function space of bounded functions V( X , c, c) and the model class , c, c) implies the non-negativity of our approximation in (38).Numerically, those spaces are not immediately accessible in the context of tensor regression.We mention two possible methods that can be implemented to handle the numerical optimization in M: -Constrained optimization: The boundedness in the model class can be translated into constraints on the tensor components of the tensor train model.-Square root trick: One can obtain a tensor train approximation ĝ of f0 , as e.g.applied in Cui and Dolgov (2021).Subsequently squaring the result yields the desired non-negative tensor train approximation ĝ of f0 .However, taking the square of an extended tensor train increases the rank and the number of basis elements in the representation.Therefore, careful adjustments have to be made to ensure that the analysis of Sect. 4 is still valid.Due to the upper bound of f0 in Assumption 1 and the fact that for any ĝ ∈ M it holds for some C > 0 that Then, the error analysis can be carried out as presented with adapted model classes that take care of the increasing rank.

Validation experiment 1: Gaussian density
In this experiment we confirm the theoretical results of Sect. 4 and verify the numerical algorithm.Even though the examined approximation of a Gaussian density is not a challenging task for the proposed algorithm, it can be seen as the most fundamental illustration, revealing the possible rank-1 structure of the perturbed prior under optimal transport.We consider the posterior density determined by a Gaussian density with covariance matrix Σ ∈ R d,d and mean μ ∈ R d given by dπ dλ where C = (2π) − d /2 det Σ − 1 /2 is the normalizing factor of the multivariate Gaussian.We set the covariance operator such that the Gaussian density belongs to uncorrelated random variables, i.e.Σ exhibits a diagonal structure, and it holds for some 0 < σ 1 that Σ = σ 2 I .This Gaussian setting has several benefits when used as validation.On the one hand, we have explicit access to the quantities that are usually of interest in Bayesian inference like the mean, covariance, normalization constant and marginals.On the other hand, the optimal transport to a standard normal density is given by an affine linear function, defined via mean μ and covariance Σ as proposed in Remark 2. We subsequently employ the multivariate polar transformation from Example B.1 and expect a rank-1 structure in the reconstruction of the local approximations of the (perturbed) prior.The remainder of this section considers different choices of σ ∈ R and d ∈ N and highlights the stability of our method under decreasing variance (i.e. with higher density concentration) and increasing dimension.The approximations are compared with their exact counterparts.More specifically, the error of the normalization constant is observed, namely and the relative error of the mean and covariance in the Euclidean and Frobenius norms Additionally the convergence in terms of the Kullback-Leibler divergence (43) and the Hellinger distance (42) are analyzed.Computing the Kullback-Leibler divergence is accomplished by Monte Carlo samples (x i ) N KL i=1 of the posterior (i.e. in this case the multivariate Gaussian posterior) to compute the empirical approximation The index h generically denotes the employed approximation (38).A similar computation can be carried out for the Hellinger distance such that In the numerical experiments, the convergence of these error measures is depicted with respect to the number of calls to the forward model (i.e. the Gaussian posterior density), the discretization of the radial component ρ ∈ [0, ∞) in the polar coordinate system and the number of samples on each layer X , = 1, . . ., L, for fixed L ∈ N.
In Table 1 err Z is depicted for different choices of σ and d.The experiment comprises radial discretizations 0 = ρ 0 < ρ 1 < . . .< ρ L = 10 with L = 19 equidistantly chosen layers and 1000 samples of f 0 on each resulting subdomain X .The generated basis (67) contains polynomials of maximal degree 7 in ρ , = 0, . . ., L, and constant functions in every angular direction.The choice of constant functions relies on the assumption that the perturbed prior that has to be approximated corresponds to the polar transformation of (94), which is a function in ρ only.Additional numerical test show that even much fewer samples and a larger basis lead to the assumed rank-1 structure.It can be observed that the approximation quality of Z is invariant under the choice of σ and fairly robust with the dimension d, which is expected since the transformation is exact and the function to reconstruct is a rank-1 object.
In Fig. 2 we show the convergence behavior of the developed method in terms of the number of calls of the posterior density f .Here, the presented low-rank surrogate is constructed on an increasing number of layers.Taking 100 samples on each layer starting from only one, i.e. from L = 1 up to L = 50 on [0, 10], the VMC tensor reconstruction is carried out.This increase in the number of layers translates directly to an increase in posterior calls.Figure 2 (left) shows the convergence in terms of the Kullback-Leibler divergence and the Hellinger distance.In accordance with Theorem 1, the Kullback-Leibler divergence converges faster.For the computation, N KL = N Hell = 1000 samples of the posterior f are drawn and we show the mean and 90% quantile of 30 repeated evaluations, where each run includes the recomputation of the VMC approximation.For a quantitiative comparison in terms of posterior statistics, in Fig. 2 (right) we show the relative error err Σ for the presented algorithm (TT, blue), a reference Monte Carlo sampling using samples from the posterior f (MC, green) and a Hamiltonian Monte Carlo sampler (HMC, red) as it comes out of the box with Matlab.HMC is a sensible choice as a reference method since the tuning phase of HMC involves the computation of the MAP and Hessian to improve the exploration phase of the sample chains, which is in some sense similar to employed affine linear transport in our algorithm.The HMC tuning phase usually took more than 1000 posterior evaluations, which are not counted in Fig. 2. As a result, we observe fast convergence of our method which takes advantage of the regularity of f in comparison to the sampling techniques, which are inherently limited by the 1/ √ N Monte Carlo convergence.HMC performs as well as MC, which comes from the fact that the estimated MAP and Hessian push the sampler into the posterior region.To emphasize the stability of the results, we repeated the VMC approximation 30 times and again show the mean and 90% quantile.For the sampling methods, we repeated HMC and MC 1000 times each and show the mean of the relative covariance approximation only.Concerning the stagnation of err Σ we suspect a precision problem in the computation, which is confirmed by the small variance.Error of the normalization constant computed via a TT surrogate to Z = 1 Nevertheless, an approximation of around four magnitudes smaller than MC and HMC for the covariance is achieved.

Validation experiment 2: Perturbation of exact transport
In the following experiment we consider a so-called "banana example" as posterior density, see e.g.Marzouk et al. (2016).Let f 0 be the density of a standard normal Gaussian measure and let T Σ be the affine transport of N (0, I ) to the Gaussian measure N (0, Σ).Furthermore, set The exact transport T from N (0, I ) to the curved and concentrated banana distribution with density f is then given by Note that the employed density can be transformed into a Gaussian using a quadratic transport function.For this experiment, we employ transport maps T of varying accuracy for the pull-back of the posterior density to a standard Gaussian.
In particular, we use an approximation T1 (obtained with Baptista et al. 2015Baptista et al. -2018) ) of the optimal affine transport T 1 and the quadratic transport T to build an approximation T given as convex combination For t = 1, the transport map is optimal since it generates the desired reference density.For 0 ≤ t < 1, a perturbed prior density is obtained with perturbation strength determined by t.The impact of the perturbed transport is visualized in Fig. 3.It can be observed that the transformed perturbed prior is not of rank-1 when the transformation is inexact.Furthermore, the difference between the target prior and the perturbed prior is imminent, which implies that e.g. a Laplace approximation to the considered banana density would neglect possible important features of the distribution.In Fig. 4 we illustrate the impact of an inexact transport on the approximation results in terms of err μ and err Σ .For the considered target density, mean and covariance are known analytically and hence no reference sampling has to be carried out.We additionally employ an MCMC sampling to show the improvement due to the additional low-rank reconstruction.For the optimal transport map, one observes that the surrogate reconstruction reduces to the approximation of a rank-1 Gaussian density, which can be done efficiently with few evaluations of f .If the transport is only linear and inaccurate, results comparable to MCMC are achieved.For a (somewhat) more accurate transport, the low-rank reconstruction leads to (already) drastically improved estimates.

Bayesian inversion with log-normal Darcy model
Random partial differential equations (PDEs), i.e.PDEs with correlated random data, play an important role in the popular field of Uncertainty Quantification (UQ).A prominent benchmark example is the (stationary) ground water flow model, also called the Darcy problem, as e.g.examined in Eigel et al. (2014Eigel et al. ( , 2017Eigel et al. ( , 2020)).This linear second order PDE model on some domain D ⊂ R d , d = 1, 2, 3 is determined by a forcing term g ∈ L 2 (D) and the random quantity a(y) ∈ L ∞ (D), which for almost every y ∈ Y models a con-ductivity or permeability coefficient.The physical system is described by − div (a(y)∇q(y)) = g in D, q(y)| ∂ D = 0, (101) and the solution q(y) ∈ V := H 1 0 (D) corresponds to the system response.Pointwise solvability of (101) for almost every y ∈ Y is guaranteed by a Lax-Milgram argument.For details we refer to Schwab and Gittelson (2011).
For the applications in this article, we employ a truncated log-normal coefficient field for some fixed (a k ) d k=1 with a k ∈ L 2 (D) denoting planar Fourier cosine modes and the image of some random variable with law N (0, I ) denoted by y = (y k ) d k=1 ∈ Y .A detailed description and an adaptive Galerkin approach to solve this problem can be found in Eigel et al. (2020).
For the inverse problem, the observation operator is modeled by 144 equidistantly distributed observations in D = [0, 1] 2 of the solution q(y * ) ∈ H 1 0 (D) for some y * ∈ Y = R d , which is drawn from a standard normal distribution.Additionally, the observations are perturbed by a centered as the target density of the measure π on Y according to (1).This target density is now approximated by the developed algorithm.
To obtain the desired relative error quantities as above, we employ reference computations that involve adaptive quadrature for the two dimensional example in Fig. 5 and Markov Chain Monte Carlo integration with 10 6 steps of the chain and a burn-in time of 1000 samples for the experiment in Fig. 6.We point out that the MCMC chains are started in an area defined by the numerically obtained MAP point used for the employed transport map.This is in some sense similar to the HMC algorithm considered in the previous example.For the reconstruction algorithm, an affine linear transport is estimated by Hessian information of the log-likelihood and on every layer we employ 100 samples.As above, this corresponds to taking initially one layer resulting in 100 posterior calls.We increase the number of layers up to L = 25.The respective relative errors are displayed in Fig. 6.
The stagnation of the graphs in Fig. 5 is on the one hand governed by the observation noise and on the other hand explicable by a non-optimal reference solution since the TT approximation yields results equivalent to an adaptive quadrature when taking L = 5 layers of refinement and thus a total of 500 samples.
The improvement of the mean and covariance estimate by the low-rank approach can already be observed for a low sample number.We note that the Monte Carlo estimate does not allow for an adequate computation of the empirical covariance, which therefore is left out of the comparison.

Conclusion
We propose a novel approach to approximate probability densities with high accuracy, combining the notion of transport maps and low-rank functional representations of auxiliary (perturbed) reference densities.Based on a suitable class of transformations, an approximation with respect to a finite tensorized basis can be carried out in extended hierarchical tensor formats.This yields a compressed representation for an efficient computation of statistical quantities of inter-Fig.5 Comparison of the computed reference and the low-rank surrogate of (1) normalization constant (err Z ), and (2) mean (err μ ).For the Darcy setting with d = 2 we observe 144 nodes in the physical domain.The measurements are perturbed by Gaussian noise with deviation η = 1e − 7. We employ an adaptive quadrature in the two dimensional space to obtain the reference quantities.The stagnation of the graphs are due to non-optimal reference solutions.More precisely, the TT approximation yields equivalent results adaptive quadrature when taking 5 nodes of refinement est (e.g.moments or marginals) in a sampling free manner.In this work, the multivariate polar transformation is used as a particular rank-1 stable transformation.The method requires point evaluations of the perturbed reference density (up to a multiplicative constant).The approach can hence be applied to not normalized posterior densities in the context of Bayesian inversion.An a priori convergence analysis is examined and we illustrate the performance of the method via an inverse problem with a log-normal Darcy forward model.A comparison with classical MCMC illustrates the superior convergence in terms of the moment accuracy relative to the number of posterior evaluations.Future research will be concerned with -application: usage of the approximated densities for subsequent computations e.g. with SGFEM, -analysis: Given a function f0 it has to be examined which rank-1 stable transformations Φ lead to a low-rank function f0 • Φ. sin i θ i .
The structure of the transformation and its jacobian leads to the following result.
Proposition 1 The multivariate spherical coordinate transformation is rank-1 stable.
Then the generalized hypercylindrical transformation where Φ 2 : R d → R d denotes the d-dimensional spherical transformation from "Appendix B.1".This construction leads to the following proposition.

Proposition 2
The generalized hypercylindrical coordinate transformation is rank-1 stable.

Fig. 1
Fig. 1 Overview of the presented method sketching the different involved transformations and approximations with references to the respective equations

Definition 3 (
covering number) Let > 0. The covering number ν(M, ) denotes the minimal number of open balls of radius with respect to the metric d V( X ,c,c) needed to cover M. Lemma 6 Let ι be defined as in (49) or (53).Then there exist C 1 , C 2 > 0 only depending on the uniform bound and the Lipschitz constant of M given in Lemmas 4 and 5, respectively, such that for > 0 and N ∈ N denoting the number of samples in the empirical cost functional in (35) it holds with ν denoting the covering number from Definition 3 and δ ( , N ) ≤ 2 exp(−2 2 N /C ι ).Proof We first prove (61) for = Hell.Note that | √ a − √ b| ≤ √ |a − b| for a, b ≥ 0 and with Lemma 1 it holds d 2 Hell

Fig. 4
Fig. 4 Convex combination of affine and quadratic transport for the banana posterior.Affine linear map (t = 0 top left), transport with t = 0.25 (top right), t = 0.5 (bottom left) and exact quadratic map (t = 1, bottom right).Error quantities err μ and err Σ for the employed

Fig. 7 Fig. 8
Fig. 7 Illustration of affine transport: translation, rotation and rescaling Fig. 8 Illustration of quadratic transport: affine properties and bending

Fig. 9
Fig. 9 Illustration of one-dimensional blow-up of a local feature open and bounded sets.A diffeomorphism Φ : X → X is called rank 1 stable if Φ and the Jacobian |J Φ | have rank 1, i.e. there exist univariate functions Φ

Table 1
Numerical approximation of Z in the Gaussian example