{'fraction_non_alphanumeric': 0.038642556353728014, 'fraction_numerical': 0.014243249938072826, 'mean_word_length': 4.525732274842595, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 13, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A transport-based multifidelity preconditioner for Markov chain Monte Carlo

Markov chain Monte Carlo (MCMC) sampling of posterior distributions arising in Bayesian inverse problems is challenging when evaluations of the forward model are computationally expensive. Replacing the forward model with a low-cost, low-fidelity model often significantly reduces computational cost; however, employing a low-fidelity model alone means that the stationary distribution of the MCMC chain is the posterior distribution corresponding to the low-fidelity model, rather than the original posterior distribution corresponding to the high-fidelity model. We propose a multifidelity approach that combines, rather than replaces, the high-fidelity model with a low-fidelity model. First, the low-fidelity model is used to construct a transport map that deterministically couples a reference Gaussian distribution with an approximation of the low-fidelity posterior. Then, the high-fidelity posterior distribution is explored using a non-Gaussian proposal distribution derived from the transport map. This multifidelity"preconditioned"MCMC approach seeks efficient sampling via a proposal that is explicitly tailored to the posterior at hand and that is constructed efficiently with the low-fidelity model. By relying on the low-fidelity model only to construct the proposal distribution, our approach guarantees that the stationary distribution of the MCMC chain is the high-fidelity posterior. In our numerical examples, our multifidelity approach achieves significant speedups compared to single-fidelity MCMC sampling methods.


Introduction
Bayesian inference provides a framework to quantify uncertainties in the solutions of inverse problems [30,58,59]. The Bayesian approach to inverse problems combines observed data, a forward model that maps parameters to observations, a prior distribution on the parameters, and a statistical model for the mismatch between model predictions and observations to define the posterior distribution via Bayes' theorem. The posterior distribution characterizes the parameter values and their uncertainties, given these ingredients. Practically "solving" a Bayesian inverse problem, however, entails exploring the posterior distribution, e.g., computing posterior expectations. A flexible and widely used approach for exploring posterior distributions is to draw samples with Markov chain Monte Carlo (MCMC) methods [60,26]. Using these methods, the forward model typically must be evaluated multiple times at different parameter values for each sample that is drawn, such that MCMC sampling quickly becomes computationally infeasible if each forward model solve is expensive.
In this paper, we propose a multifidelity preconditioner to increase the efficiency of MCMC sampling. Our multifidelity approach exploits low-cost, low-fidelity models to construct a proposal distribution that approximates the posterior distribution at hand, then uses this proposal distribution to perform MCMC sampling of the original (high-fidelity) posterior distribution. The proposal distribution is derived from a transport map that transforms the potentially complex posterior distribution into another distribution from which samples can be drawn more easily; in particular, we seek a map that transforms the posterior into a more Gaussian and more isotropic distribution. See [40,36,44] for an introduction to transport maps in the context of Bayesian inverse problems. The key idea in the present work is to construct an invertible map using low-fidelity models, but to apply it to the high-fidelity posterior. The map then acts as a preconditioner for the high-fidelity posterior, preserving information (due to its invertibility) while enabling MCMC sampling to proceed more efficiently. An alternative but equivalent perspective is that pushing a simple proposal through the inverse of this transport map-for instance, using a Metropolis independence sampler with a standard Gaussian proposal-yields potentially non-Gaussian and tailored proposals that can be efficient for the high-fidelity posterior. The end result is the same: more efficient MCMC sampling, allowing the number of high-fidelity forward model evaluations to be reduced. Since the low-fidelity model is used only for preconditioning (or equivalently, constructing the proposal), the stationary distribution of the chain obtained with our multifidelity approach is guaranteed to be the posterior distribution corresponding to the high-fidelity model.
There is a long tradition of exploiting low-fidelity models to speed up MCMC sampling for Bayesian inverse problems. There is work [65,2,38,37,33,13] that replaces the high-fidelity model with a lowfidelity model (see also the survey [25]); because of this replacement, however, samples are drawn from a different distribution-the posterior induced by the low-fidelity model. Under certain assumptions, it has been shown that the posterior distribution induced by the low-fidelity model converges to the original/highfidelity posterior distribution as the low-fidelity approximation is refined [39,18]. Other approaches adapt low-fidelity models [19,32] over a finite interval of posterior exploration, or quantify the error introduced by sampling the low-fidelity posterior distribution [21,35]. Yet another family of approaches incrementally and infinitely refines approximations of the forward model on-the-fly during MCMC sampling [16,15]; under appropriate conditions, these schemes guarantee that the MCMC chain asymptotically samples the highfidelity posterior distribution.
Instead of replacing the high-fidelity model with low-fidelity models, multifidelity methods combine highand low-fidelity models. The aim is to leverage low-fidelity models for speeding up computations while allowing occasional recourse to the high-fidelity model to establish accuracy guarantees [48]. A variety of multifidelity methods have been developed for uncertainty propagation [8,7,41,23,11,45,47,49,46]; see the survey [48]. For the solution of Bayesian inverse problems, there are multi-stage MCMC methods that aim to reduce the number of high-fidelity model evaluations by first screening proposed moves with lowfidelity models [14,22]. Another line of work builds on hierarchies of low-fidelity models, typically derived from different discretizations of partial differential equations (PDEs) underlying the high-fidelity model, to reduce sampling costs [4,20,31].
An alternative to multi-stage and hierarchical methods for increasing the efficiency of MCMC sampling is to use notions of transport to construct more effective proposal distributions. Effective proposals in MCMC should reflect the local or global geometry of the target distribution [26,9]. In keeping with this idea, Parno et al. [44] use transport maps to precondition MCMC sampling. As described earlier, preconditioning involves constructing transport maps that "Gaussianize" the target distribution, such that it can be sampled more effectively by standard MCMC algorithms; these maps thus encode the geometry of the target. Parno et al. [44] build and refine such transport maps in an online fashion during MCMC sampling, as more and more samples are obtained. This approach can be seen as a form of adaptive MCMC [29,51], wherein a (non-Gaussian) proposal distribution is adapted as the MCMC sampling proceeds. The transport maps in [44] are constructed via the solution of a convex and separable optimization problem, which is simple and fast to obtain numerically. However, this sampling approach faces some of the usual pitfalls of adaptive MCMC. One of these issues is that a certain amount of initial mixing is necessary for the adaptation to be effective, because the online-adapted transport map depends on past samples.
The transport approach of Moselhy et al. [40] instead follows an offline/online decomposition of the computation. In an offline phase, a reference distribution is selected from which independent samples can be drawn cheaply and then a transport map is constructed that pushes forward the reference distribution to the posterior distribution. Then, in the online phase, the transport map is used to transform samples from the reference distribution into samples from the posterior distribution. If a large number of samples are transformed in the online phase, then the one-time high cost of constructing the transport map offline is compensated. The optimization problem for constructing the transport map can employ derivative information from the posterior density, and thus yields accurate maps even if the posterior is concentrated. Furthermore, no samples of the posterior distribution are necessary to solve the optimization problem, in contrast with adaptive MCMC techniques such as [44]. However, evaluations of the objective of the optimization problem entail evaluations of the unnormalized posterior density and thus of the forward model; moreover, this problem is in general not separable across dimensions and not convex (unless the posterior density is log-concave). As a result, map construction in [40] is typically far more computationally demanding than map construction in [44]. Furthermore, since the construction of the transport map involves numerical approximations, the pushforward of the reference distribution by the map in general only approximates the posterior distribution.
We propose a multifidelity approach that combines several of the advantages of the approaches introduced by Parno [44] and Moselhy [40]. We construct transport maps building on the offline/online approach of [40] and mitigate the high computational costs of solving the optimization problem by relying on low-cost, low-fidelity models. Our approach can exploit a wide range of low-fidelity models, including projection-based reduced models [55,53,27,3,12], data-fit interpolation and regression models [24], machine-learning-based models [17,66,54,61], and simplified-physics models [1,42]. Then, in the online phase, the transport map is used to precondition MCMC sampling of the high-fidelity posterior, as in [44]. The corresponding MCMC scheme is ergodic for the high-fidelity posterior. In other words, the Metropolis step corrects errors that otherwise would be introduced by relying on the transport map alone to push forward the reference to the posterior. Thus, we obtain a multifidelity approach that uses low-fidelity models to speed up computations while making occasional recourse to the high-fidelity model to establish convergence, in the sense that the stationary distribution of our MCMC chain is the high-fidelity posterior distribution. In contrast to [44], our MCMC algorithm is not adaptive, because the map is built once offline and then stays fixed during MCMC sampling. In particular, no samples of the high-fidelity posterior distribution are needed to construct the transport maps. Section 2 describes the problem setup and briefly reviews transport maps in the context of Bayesian inverse problems. In Section 3, we introduce our multifidelity preconditioner and discuss it in the context of sampling with the Metropolis-Hastings algorithm. Section 4 demonstrates our multifidelity approach on two examples where we achieve significant speedups compared to using the high-fidelity model alone. Conclusions are given in Section 5.

Preliminaries
Section 2.1 and Section 2.2 define our Bayesian inverse problem setting and describe the Metropolis-Hastings algorithm for MCMC sampling. We refer to, e.g., [30,58,59], for details on Bayesian approaches to inverse problems. Section 2.3 discusses transport maps for coupling probability distributions in the context of Bayesian inverse problems. The problem description is given in Section 2.4.

Bayesian inverse problems
Consider the high-fidelity forward model (parameter-to-observable map) G : P → Y, with parameter θ ∈ P, where P ⊆ R d , and observable y ∈ Y, where Y ⊆ R d ′ . Thus, the parameter θ and the observable y are d-dimensional and d ′ -dimensional vectors, respectively. Consider now observed data y = G(θ * ) + ǫ , corresponding to some parameter value θ * ∈ P. The noise ǫ is assumed to be a realization of a zero-mean Gaussian random variable with covariance matrix Σ ǫ ∈ R d ′ ×d ′ . Define the data-misfit function with the Euclidean norm · 2 . The likelihood function L y : P → R is Combing the prior distribution with density π 0 and the likelihood L via Bayes' theorem gives the posterior density up to a normalizing constant π(θ) ∝ L y (θ)π 0 (θ) .
Above and for the remainder of this paper, we assume that the prior and posterior measures are absolutely continuous with respect to the Lebesgue measure.

The Metropolis-Hastings algorithm
MCMC methods are widely used to sample posterior distributions that arise in Bayesian inverse problems. The Metropolis-Hastings algorithm defines a wide class of MCMC methods, on which we will build in the following. Algorithm 1 describes the Metropolis-Hastings approach. In our Bayesian inverse problem setting, inputs are the likelihood L, the prior density π 0 , a proposal density q, and the number of iterations M ∈ N.
In each iteration, a proposal sample θ ′ is drawn from the proposal distribution, which may depend on the previous sample θ i−1 . Then, the acceptance probability α(θ i−1 , θ ′ ) is computed, which requires evaluating the likelihood L and prior π 0 at the candidate sample θ ′ . The proposal sample θ ′ is accepted with probability α(θ i−1 , θ ′ ) and rejected with probability 1 − α(θ i−1 , θ ′ ). This process is repeated for M iterations and the samples θ 1 , . . . , θ M are returned. The corresponding Markov chain is, by construction, reversible for π and thus has π as a stationary distribution. With some relatively simple additional conditions on π and q, one can show that the chain converges to π, from any starting point; see [50] for a full discussion. In general, the samples produced by MCMC are correlated; this correlation inflates the variance of any expectations estimated with MCMC samples, relative to an expectation estimated with uncorrelated Monte Carlo sample sets of the same size. The efficiency of an MCMC sampler can thus be measured with the effective sample size (ESS) of any sample set it produces, which is inversely proportional to the integrated autocorrelation time of the chain [34]. To define the ESS, consider a function f : P → R that is measurable with respect to the Lebesgue measure, and let us assume we are interested in estimating the expected value with respect to the posterior distribution π. Consider now the Monte Carlo estimator of Draw candidate θ ′ from proposal q(·|θ i−1 )

5:
Compute acceptance probability Set the sample θ i to In other words, n * ≤ n is the number of independent Monte Carlo samples from π that would be required to obtain an estimator with the same variance as E n f . It can be shown that n * = n/τ , where τ is the integrated autocorrelation time associated with the chain {f (θ t )} t . Better MCMC mixing corresponds to smaller τ and larger ESS.

Transport maps
The use of transport maps in the context of Bayesian inference was introduced in [40]. In particular, [40] proposed a variational Bayesian approach involving transport. Rather than using importance sampling or MCMC to characterize the posterior distribution, this approach seeks a transport map that pushes forward a tractable "reference" distribution to the posterior, such that samples drawn from the reference and acted on by the map are distributed according to the posterior. Below we follow [36] to introduce the notion of transport maps.

Definition of transport maps
Let µ π and µ η be two probability measures on R d that are absolutely continuous with respect to the Lebesgue measure. In the following, µ π is the target measure, which corresponds to the posterior distribution in our case, and µ η is the reference measure, which typically is a Gaussian or another distribution from which we can draw independent samples efficiently. The probability density function corresponding to the target measure µ π is the posterior π, and the probability density function corresponding to µ η is denoted by η. A transport map is a function T : R d → R d that pushes forward the reference µ η to the target µ π , which we write as and which means that for any Borel set A ⊆ R d , it holds µ π (A) = µ η (T −1 (A)). Existence of such maps is ensured by the absolute continuity of the reference and the target measure. There may be infinitely many transport maps that push forward a given reference to the target of interest. Uniqueness can be enforced by introducing a cost function that is minimized while imposing the constraint (1). This construction leads to the notion of optimal transport ; see, e.g., [62,63,64]. Instead of introducing a cost function to regularize the problem of finding a transport map, we directly impose structure on the map T as in [40,44,36]. In particular, we will seek lower triangular maps that are monotone increasing. The lower triangular structure of T is as follows where Monotonicity in this context corresponds to the condition that T i is a monotone increasing function of ϑ i , for all i = 1 . . . d. This condition ensures that ∇T 0 and det ∇T ≥ 0; see [36,57] for more detail. Since we assume that the reference and the target measures are absolutely continuous, existence and uniqueness of such a lower-triangular transport map is guaranteed; this map is in fact the Knothe-Rosenblatt rearrangement [6, 10, 52].

Numerical approximations of transport maps
Following [40], we will obtain numerical approximations of the Knothe-Rosenblatt rearrangement by solving an optimization problem. Let T be a finite-dimensional subspace of the space of all smooth lower triangular functions (2) from R d into R d . Then, an approximation T ∈ T of a transport map T can be obtained via numerical optimization over the coefficients of the representation of T in a basis of T . To set up the optimization problem, consider the pullback ( T −1 ) ♯ µ π ≡ T ♯ µ π of µ π through a map T ; the density of this pullback measure can be written as where | det ∇ T (ϑ)| is the absolute value of the determinant of the Jacobian ∇ T (ϑ) of T at ϑ. Note that the functions in T are smooth in the sense that ∇ T exists and is sufficiently regular; see, e.g., [36,57]. Let denote the Kullback-Leibler (KL) divergence of π 1 from π 2 (where π 1 and π 2 in the log term are densities). Then, a solution T * ∈ T of the optimization problem is an approximation of a transport map that pushes forward the reference measure µ η to the target µ π . The constraint ∇ T ≻ 0 means that the Jacobian of T is positive definite. If the approximation space T is sufficiently rich such that D KL (η|| η) = 0, then we have T * ♯ µ π = µ η and T * ♯ µ η = µ π [36]. The KL divergence is not symmetric. The direction of the KL divergence here is chosen such that the expected value is taken with respect to the reference µ η , which is selected so that it can easily be sampled. Furthermore, the objective in (4) can be minimized without knowledge of the normalizing constant of π. To see this, transform the objective D KL (η|| η) into where we used the definition of η in (3). The expectation E η [log η] is independent of the map T and therefore can be ignored when minimizing the objective of (4). Similarly, the normalizing constant of π leads to a constant term in (5) that is independent of T and therefore unnormalized evaluations of π are sufficient to optimize (4).

Problem description
We identify two challenges of directly relying on a solution T * of (4) to solve a Bayesian inverse problem, as proposed in [40]. First, the KL divergence objective (5), which contains an expectation with respect the the reference measure µ η , is typically estimated with a Monte Carlo method because no closed form expression is available for general π. Thus, the (unnormalized version of the) pullback density (3) must be evaluated at a potentially large number of samples for each optimization iteration. Since each evaluation of the pullback entails an evaluation of the (unnormalized) posterior distribution π, and thus of the forward model G, the optimization can become computationally intractable if G is expensive to evaluate. Second, the push forward T * ♯ µ η is only an approximation of the posterior µ π . In particular, approximation errors may follow from the choice of the finite-dimensional approximation space T , the finite number of Monte Carlo samples used to discretize the expectation in the objective, and any other errors in the numerical optimization. While this error can be estimated (see [40,36]), reducing this error to an arbitrarily small threshold-e.g., by enriching T -can be computationally expensive.

Multifidelity preconditioned Metropolis-Hastings
We propose a multifidelity preconditioned Metropolis-Hastings (MFMH). Similar to the approach of Moselhy [40], our multifidelity approach consists of an offline and an online phase. In the offline phase, a transport map is constructed rapidly by using a low-cost, low-fidelity approximation of the high-fidelity model. Then, in the online phase, a proposal distribution is derived from the transport map and used to sample the high-fidelity posterior via a modified Metropolis-Hastings algorithm, following the work of Parno [44]. If the low-fidelity model is an accurate approximation of the high-fidelity model, and if the transport map captures the essential structure of π, then the ESS of this MFMH will be higher and fewer online evaluations of the high-fidelity model will be required to achieve a given accuracy. In any case, any error associated with the transport map is corrected by Metropolization in the online phase; samples are thus drawn (asymptotically) from the posterior distribution corresponding to the high-fidelity model G. In other words, the MFMH approach offers the same convergence guarantees as standard MCMC.

Approximation spaces and numerical optimization
This section provides details on the numerical construction of transport maps from the reference measure to the target measure.

Integrated squared parametrization
Each component function T i , i = 1, . . . , d of the approximation T of the map T defined in (2) is parameterized with the integrated-squared ansatz [5], which is similar to the integrated-exponential ansatz introduced in [57]. Component function T i is parameterized as is a parameter vector and T i , respectively. The integrated-squared parameterization (6) guarantees that the map T is monotone and therefore automatically satisfies the constraint ∇ T ≻ 0 in (4) [57].

Approximation space
Following [44,36], we represent the functions T In the following, the multivariate polynomial functions (8) are simply the monomials. Note that other polynomial families can be used, e.g., Hermite polynomials. (One might also use Hermite functions, as in [57], for better control of tail behavior.) Consider now the sets J i ⊂ N d for i = 1, . . . , d which correspond to the total-degree polynomials of maximal degree ℓ ∈ N. The constraint j k = 0 , ∀k > i in the definition of J i in (9) imposes the lower-triangular structure of the map T as defined in (2). The set J i leads to the definition of the approximation space T i and T where β i,j with j ∈ J i , respectively. Note that we combine β (L) i and β (R) i into the vector β i as in defined in (7). The approximation space T is the product space where ⊕ and ⊗ denote the sum and the product of two spaces.

Numerical solution of optimization problem
Consider the optimization problem (4) and the transformation (5) that shows that it is sufficient to minimize with respect to T . Note that evaluations of the unnormalized version of π are sufficient to minimize (12); see Section 2.3.2. We have dropped the absolute value above since the maps T are guaranteed to be monotone, via our parameterization. We now replace the expected value with its sample-average approximation [56], i.e., a Monte Carlo estimator employing independent draws {ϑ i } n i=1 from the reference distribution η. Making the dependence on the coefficients β 1 , . . . , β d explicit, we obtain the optimization problem min β1,...,β d Our optimization problem (13) is unconstrained because the constraint ∇ T ≻ 0 of (4) is automatically satisfied via the squared-integrated parameterization (6).

Constructing transport maps with low-fidelity models
We propose to approximate the high-fidelity forward model G with a low-cost, low-fidelity model to reduce the computational costs of constructing a transport map via the optimization problem (13). Note that replacing evaluations of G with evaluations of low-fidelity models in (13) will introduce an error that we must correct later; see Section 3.4. Let G : P → Y be a low-fidelity approximation of G that maps the parameter onto an observable. The model G gives rise to a low-fidelity potential function and to the low-fidelity likelihood function L y : P → R L y (θ) = exp − Φ y (θ) .
We now use the low-fidelity posterior density π in our variational construction of the map T ∈ T . Consider the density η that is the pullback of the low-fidelity posterior density π through a map T , We can find T by minimizing D KL (η|| η); this yields the following optimization problem, analogous to (13) but with π replacing π, min β1,...,β d with the coefficients β 1 , . . . , β d defining T * ∈ T .

Constructing deep transport maps
The optimization problem (14) finds a transport map in the approximation space T . Thus, if the approximation space T is chosen too coarse, then the pullback of the target distribution through the map is only a poor approximation of the reference distribution. Instead of choosing richer approximation spaces to find more accurate transport maps, which would lead to a large number of coefficients to be optimized for in (14), the work [40,43] proposes to take compositions of transport maps in coarse approximation spaces. Let T (1) be a transport map that pushes forward the reference µ η onto the posterior µ π and let T (1) ∈ T be a numerical approximation of T (1) derived with optimization problem (14). Since T (1) is an approximation of the map T (1) , we obtain µ π ≈ µ (1) ♯ µ η . Following [44], to account for the discrepancy between µ π and µ (1) π , a second map T (2) ∈ T is constructed to push forward µ (1) π to µ π , i.e., µ π ≈ µ π . The aim is that the map T (2) should capture only a small correction from µ (1) π to µ (2) π , so that the composition T (1,2) = T (2) • T (1) more accurately pushes forward µ η to µ π (in the sense of KL divergence). This process is repeated k times to obtain the "deep" map Note that (15) is only one possible way of constructing deep transport maps. Another possibility is to reverse the order of (15) and compose maps from the right, rather than from the left, by keeping the µ η as the reference and pulling back the target µ π by the current map approximation. This construction is discussed in [57].

MFMH algorithm
We now use the transport map T * to derive an MCMC proposal with the aim of improving sampling efficiency (e.g., reducing integrated autocorrelation time) and at the same time guaranteeing that the stationary distribution of the Markov chain is the posterior distribution corresponding to the high-fidelity model G.
Since the transport map T * is monotone and lower triangular, as enforced by the approximation space T , the inverse map S * = T * −1 is cheap to evaluate by solving d one-dimensional root-finding problems. We refer to [44,36] for details on how to invert lower-triangular, monotone maps. The inverse map S * pushes forward the true posterior µ π onto a distribution that approximates the reference distribution. In other words, S * ♯ µ π = µ η only if the map were exact, which in general it is not; otherwise, S * ♯ µ π is simply closer to µ η (in the sense of KL divergence) than µ π was. If µ η is chosen to be a Gaussian, then the map S * approximately "Gaussianizes" µ π . As shown in [44], and as will be demonstrated with our numerical results in Section 4, MCMC sampling from the approximate reference distribution S * ♯ µ π generally yields higher ESSs for a given computational effort than sampling directly from the posterior distribution. Crucially, because the map S * = T * −1 is invertible, samples from the approximate reference distribution can be exactly (up to machine precision) pushed forward to samples from the posterior distribution via the map T * .
Our MFMH algorithm is summarized in Algorithm 2 and follows the preconditioned Metropolis-Hastings algorithm introduced in [44]. The MFMH algorithm has the same inputs as the single-fidelity Metropolis-Hastings in Algorithm 1, except that additionally the maps T * and S * are required. The current state θ i−1 is mapped with S * onto the approximate reference to obtain ϑ i−1 . Then, a candidate sample is drawn from the proposal distribution based on ϑ i−1 , and that candidate sample is mapped back with T * to a sample of the posterior π.
Whether the candidate sample is accepted or rejected is based on the probability α that is determined using the high-fidelity posterior density π; the maps, derived from the low-fidelity posterior, only affect the Map state θ i−1 onto reference with S * (θ i−1 ) = ϑ i−1

7:
Compute acceptance probability Set the sample θ i to end for 10: return θ 1 , . . . , θ M 11: end procedure proposal distribution. Thus it is guaranteed that the stationary distribution of the chain is the posterior distribution π. Note that there is considerable flexibility in the choice of the reference-space proposal distribution q. Below we will mostly choose q(·|ϑ i−1 ) to be a d-dimensional Gaussian, independent of ϑ i−1 . Thus our MCMC algorithm reduces to a Metropolis independence sampler [50], whether viewed on the reference space or on the target space. Given the monotonicity of the transport maps and the full support of the Gaussian reference on R d , the proposal distribution on the target space is guaranteed to dominate the posterior distribution, as required for Metropolis independence sampling to converge.

Numerical results
This section demonstrates our multifidelity approach on two examples. All runtime measurements were performed on compute nodes with Intel Xeon E5-1620 and 64GB RAM on a single core using a Matlab implementation.

Diffusion equation with reaction term
We first consider a model with diffusion and a nonlinear reaction term, where our goal is to infer the parameters of the reaction term.
We discretize (16) with finite differences on a grid with equidistant grid points and mesh width h > 0. The corresponding system of nonlinear equations is solved with Newton's method and inexact line search based on the Armijo condition. The model G h : P → Y derived with mesh width h maps from P into Y = R 12 . The components of the observable y correspond to the values of the approximated solution function at the spatial coordinates [0.25i, 0.2j] T ∈ Ω with i = 1, 2, 3 and j = 1, 2, 3, 4.
A low-fidelity model G h of G h is derived via projection-based model reduction [3]. Solutions of (16) for parameters on an 100 × 100 equidistant grid in [−π/2, −π/2] × [1, 5] ⊂ P and mesh width h are computed and a 20-dimensional reduced space with proper orthogonal decomposition is constructed. The operators corresponding to the high-fidelity G h are projected via Galerkin projection onto the reduced space and the low-fidelity model G h is obtained.

Setup of inverse problems
We set θ * = [0.5, 2] T and consider the data y = G H (θ * ) + ǫ, where H = 1/64 and ǫ adds Gaussian noise with zero mean and variance 0.0026, which corresponds to 0.1% noise with respect to the Euclidean norm of the output G H (θ * ). We then use the model G h with h = 1/32 and the corresponding low-fidelity model G h for inference. Note that the data y is computed with a discretization of the PDE (16) with mesh width H = 1/64, while we use a mesh width h = 1/32 for inference. The prior distribution is a Gaussian distribution with mean [π/4, 1.2] T and covariance matrix The low-fidelity G h is about 80 times faster to evaluate than the high-fidelity model G h for h = 1/32. We construct a transport map from the reference Gaussian distribution to the posterior corresponding to the low-fidelity model G h . The reference Gaussian distribution has zero mean and standard deviation 0.1. Note that instead we could have scaled and centered the posterior distribution to have zero mean and marginal variances of one and then used a standard Gaussian as reference distribution.
The optimization is performed with Matlab's fmincon optimizer, where the tolerance (TolX) is set to 10 −3 and where we use n = 250 samples of the reference distribution to approximate the expected value in the objective function, see Section 3.1.3. We compose two maps following Section 3.3. The approximation space for the first map T (1) corresponds to first-order polynomials (linear) and the approximation space of the second map T (2) corresponds to second-order polynomials (quadratic). The starting point for the optimization is the identity map. The transport map is then used to precondition Metropolis-Hastings as shown in Algorithm 2. The proposal is the reference distribution, i.e., we obtain an independence sampler with a proposal that is independent of the current state of the chain. Because we expect that the transport map approximately pushes the reference onto the posterior, it is reasonable to consider an independence sampler with the reference distribution as proposal. We discard every other sample, which means that we perform M = 2m iterations if we want m samples.

Euler Bernoulli beam problem
We now infer the effective stiffness of an Euler Bernoulli beam, for which a model is available on GitHub 1 .
where u : Ω → R is the vertical deflection of the beam, f : Ω → R is the load, andΩ = (0, L). The effective stiffness of the beam is given by E : Ω → R and describes beam geometry and material properties. The beam is in cantilever configuration, where the left boundary is fixed and the right boundary is free, i.e., the boundary conditions are The length of the beam is L = 1 in the following. The PDE (17) is discretized with finite differences on a mesh of N = 601 equidistant grid points in Ω. The same effective stiffness E as available in the GitHub repository 1 is used by interpolating on the grid points. Let x obs , . . . , x obs be equidistant points of the N = 601 grid points. The observation y ∈ R 41 is the displacement u at the 41 points x (1) obs , . . . , x (41) obs polluted with zero-mean Gaussian noise with variance 10 −4 .

Models for the Euler Bernoulli problem
We now derive the high-fidelity forward model. Consider the function I : R × Ω → R defined as such that there is a smooth transition from 0 to 1 at α. Define further k = 3 and let α 1 , . . . , α k+1 be the k + 1 equidistant points in Ω. Let R + = {z ∈ R : z > 0} and consider the parameter for i = 2, . . . , k and E 1 (x, θ 1 ) = θ 1 . Given a parameter θ, the function E k is a smooth approximation of the piecewise constant function k i=1 θ i I (αi,αi+1] , where I (αi,αi+1] is the indicator function of the interval (α i , α i+1 ] ⊂ R. The high-fidelity forward model G maps a parameter θ ∈ R k + onto the displacement u with effective stiffness E k at the observation points x (1) obs , . . . , x (41) obs . The same discretization as described in Section 4.2.1 is used. The low-fidelity model G is a spline interpolant of G on a logarithmically spaced grid in the domain [0.5, 4] 3 . The interpolant is obtained with the griddedInterpolant procedure available in Matlab. Extrapolation is turned off. The low-fidelity model G is about 1, 100 times faster to evaluate than the high-fidelity model G.

Setup of inverse problem
The observation y is obtained as described in Section 4.2.1. The prior is a log-normal distribution with mean [1, 1, 1] T and covariance matrix We construct a transport map T using the low-fidelity model G. The transport map T is based on quadratic polynomials and approximately maps the reference Gaussian distribution with mean [1, 1, 1] T and covariance with diagonal entries 0.1 onto the low-fidelity posterior distribution. Note that a standard Gaussian can be used as reference distribution if the posterior distribution is centered to zero mean and scaled to have marginal variances of one. This effectively would perform a linear transformation. We optimize for the coefficients of T with fmincon available in Matlab with tolerance set to 10 −8 (TolX option) and with 500 samples from the reference distribution. The initial point of the optimization corresponds to the identity map. We then  use the transport map in Algorithm 2 with a random-walk Metropolis algorithm on the reference space, i.e., with a local Gaussian proposal ϑ ′ ∼ N (ϑ i−1 , Σ ′ ). The covariance matrix Σ ′ is diagonal. To select the diagonal of Σ ′ , we ran MFMH for proposals with variance in {10 −4 , 5 × 10 −4 , 10 −3 , 5 × 10 −3 , 10 −2 , 2 × 10 −2 , 3 × 10 −2 , 4 × 10 −2 , 5 × 10 −2 , 10 −1 , 5 × 10 −1 } and then selected the variance that leads to the highest ESS, cf. Section 4.1.2. The corresponding results are reported in the following. Additionally, we consider an independence MFMH sampler where the reference distribution together with the transport map serves as a proposal distribution that is independent of the previous sample. The rest of the setup is the same as in Section 4.1.2. Figure 3a visualizes the posterior density function corresponding to the high-fidelity model. The plot indicates that the posterior distribution is a non-Gaussian distribution. We use the same procedure as in Section 4.1.2 to draw samples from the posterior with DRAM. The samples with the highest ESS are shown in Figure 3b.

Results
Let us now consider our MFMH approach. First, we consider the independence MFMH sampler that uses a proposal that is independent of the previous sample, see Figure 4a. The samples are in agreement with the samples drawn with DRAM from the high-fidelity posterior distribution, cf. Figure 3b. Figure 4b plots the samples mapped with the transport map onto an approximation of the reference Gaussian distribution, which demonstrates that the transport map T * captures well the posterior distribution. The ESSs of m ∈ {10 5 , 5 × 10 5 , 10 6 , 5 × 10 6 } samples drawn with DRAM and the independence MFMH sampler are compared in Figure 5. Our independence MFMH sampler achieves a higher ESS than DRAM in this example. Figure 5 also shows that a higher ESS compared to DRAM is achieved only if sufficiently many samples are drawn because otherwise the offline costs of constructing the transport maps are not compensated. The MFMH sampler with a local random walk proposal leads to a higher ESS than DRAM too, but the improvement is smaller compared to the independence MFMH sampler. Additionally, an even larger number of samples is necessary to compensate the offline costs of constructing the transport maps. Note that the MFMH sampler with a local random walk proposal and the DRAM sampler both use a burn-in of 10 4 samples as in Section 4.1.2.

Conclusions
Our MFMH approach leverages low-fidelity models to precondition MCMC sampling, with the aim of improving MCMC mixing while guaranteeing that the stationary distribution of the chain is the posterior distribution corresponding to the high-fidelity model. In a one-time expensive step, a low-fidelity model is used to construct a transport map that approximately maps an easy-to-sample reference distribution to the posterior distribution corresponding to the low-fidelity model. In the second step, the transport map is used to precondition the posterior distribution corresponding to the high-fidelity model during sampling with Metropolis-Hastings. Since the low-fidelity model is used for preconditioning only, the stationary distribution of the chain obtained in the second step is the posterior distribution corresponding to the high-fidelity model. Our MFMH approach achieves significant speedups compared to single-fidelity sampling with DRAM in our numerical examples.