Tiangang Cui
email: tiangang.cui@monash.edu

Olivier Zahm
email: olivier.zahm@inria.fr

Data-Free Likelihood-Informed Dimension Reduction of Bayesian Inverse Problems

Identifying a low-dimensional informed parameter subspace offers a viable path to alleviating the dimensionality challenge in the sampled-based solution to large-scale Bayesian inverse problems. This paper introduces a novel gradient-based dimension reduction method in which the informed subspace does not depend on the data. This permits online-offline computational strategy where the expensive low-dimensional structure of the problem is detected in an offline phase, meaning before observing the data. This strategy is particularly relevant for multiple inversion problems as the same informed subspace can be reused. The proposed approach allows to control the approximation error (in expectation over the data) of the posterior distribution. We also present sampling strategies which exploit the informed subspace to draw efficiently samples from the exact posterior distribution. The method is successfully illustrated on two numerical examples: a PDE-based inverse problem with a Gaussian process prior and a tomography problem with Poisson data and a Besov-B 2 11 prior.

Introduction

The Bayesian approach to inverse problems builds a probabilistic representation of the parameter of interest conditioned on observed data. Denoting the parameter and data by x ∈ R d and y ∈ R m , respectively, the solution to the inverse problem is encapsulated in the posterior distribution, which has the probability density function (pdf)

π y pos (x) = 1 π data (y) L y (x)π pr (x), (1)
where π pr (x) denotes the prior density, L y (x) is the likelihood function, and π data (y) is the marginal density of the data y that can be expressed as

π data (y) = R d L y (x)π pr (x) dx. ( 2
)
This way, one can encode the posterior into summary statistics, for example, moments, quantiles, or probabilities of some events of interest [START_REF] Kaipio | Statistical and computational inverse problems[END_REF][START_REF] Stuart | Inverse problems: a Bayesian perspective[END_REF][START_REF] Tarantola | Inverse problem theory and methods for model parameter estimation[END_REF], to provide parameter inference and associated uncertainty quantification.

In practice, computing these summary statistics requires dedicated methods to efficiently characterize the posterior distribution. Markov chain Monte Carlo (MCMC) methods [START_REF] Brooks | Handbook of Markov Chain Monte Carlo[END_REF][START_REF] Jun | Monte Carlo strategies in scientific computing[END_REF], originating with the Metropolis-Hastings algorithm [START_REF] Hastings | Monte Carlo sampling using Markov chains and their applications[END_REF][START_REF] Metropolis | Equation of state calculations by fast computing machines[END_REF], and sequential Monte Carlo methods [START_REF] Del Moral | Sequential Monte Carlo samplers[END_REF] have been developed as workhorses in this context. However, many inverse problems have high-dimensional or infinite-dimensional parameter space, which present a significant hurdle to the applicability of MCMC, SMC, and other related sampling methods in general.

The efficiency of these sampling methods, measured by the required number of posterior density evaluations, may deteriorate with the dimension of the parameter space, see [START_REF] Roberts | Weak convergence and optimal scaling of random walk Metropolis algorithms[END_REF][START_REF] Gareth | Optimal scaling for various metropolis-Hastings algorithms[END_REF] and references therein. Even with the rather strong log-concave assumption, start-of-the-art MCMC methods can still be sensitive to the dimension of the problem, see for instance [START_REF] Arnak S Dalalyan | Theoretical guarantees for approximate sampling from smooth and log-concave densities[END_REF][START_REF] Durmus | High-dimensional Bayesian inference via the unadjusted langevin algorithm[END_REF][START_REF] Dwivedi | Log-concave sampling: Metropolis-Hastings algorithms are fast[END_REF].

One promising way to alleviating the challenge of dimensionality is to exploit the effectively low-dimensional structures of the posterior distribution. Such low-dimensional structures can be used to construct certified low-dimensional approximations of the posterior distribution [START_REF] Spantini | Optimal low-rank approximations of Bayesian linear inverse problems[END_REF][START_REF] Zahm | Certified dimension reduction in nonlinear Bayesian inverse problems[END_REF] and efficient MCMC proposals that are robust in the parameter dimension [START_REF] Beskos | Geometric MCMC for infinite-dimensional inverse problems[END_REF][START_REF] Beskos | Multilevel sequential Monte Carlo with dimension-independent likelihood-informed proposals[END_REF][START_REF] Cui | Dimension-independent likelihood-informed MCMC[END_REF][START_REF] Cui | Multilevel dimension-independent likelihoodinformed MCMC for large-scale inverse problems[END_REF][START_REF] Lan | Adaptive dimension reduction to accelerate infinite-dimensional geometric Markov chain Monte Carlo[END_REF]. There exists several ways to detect low-dimensional structures. A widely accepted method is to utilize the regularity of the prior, in which the dominant eigenvectors of the prior covariance operator [START_REF] Youssef | Dimensionality reduction and polynomial chaos acceleration of Bayesian inference in inverse problems[END_REF] can be used to define such a low-dimensional subspace. This prior-based dimension reduction also plays a key role in the analysis of high-dimensional integration methods such as [START_REF] Ivan G Graham | Quasi-Monte Carlo finite element methods for elliptic PDEs with lognormal random coefficients[END_REF][START_REF] Ivan G Graham | Quasi-Monte Carlo methods for elliptic PDEs with random coefficients and applications[END_REF]. In addition to the prior regularity, the limited accuracy of the observations and the ill-posed nature of the forward model often allow one to express the posterior as a low-dimensional update from the prior. Methods such as the active subspace (AS) [START_REF] Paul G Constantine | Accelerating Markov chain Monte Carlo with active subspaces[END_REF][START_REF] Cortesi | Forward and backward uncertainty quantification with active subspaces: application to hypersonic flows around a cylinder[END_REF] and the likelihood-informed subspace (LIS) [START_REF] Cui | Likelihood-informed dimension reduction for nonlinear inverse problems[END_REF][START_REF] Cui | A unified performance analysis of likelihood-informed subspace methods[END_REF][START_REF] Zahm | Certified dimension reduction in nonlinear Bayesian inverse problems[END_REF] utilize gradients of the forward model and/or of the likelihood function in order to better identify the low-dimensional structure of the problem. We refer to [START_REF] Cui | A unified performance analysis of likelihood-informed subspace methods[END_REF][START_REF] Zahm | Certified dimension reduction in nonlinear Bayesian inverse problems[END_REF] for an overview and a comparison of the existing methods.

The success of AS and LIS relies on the computation of the gradient or the Hessian of the log-likelihood function. Since the likelihood function depends on the observed data, the resulting subspaces need to be reconstructed each time a new data set is observed. This can add a significant computational burden to the solution of inverse problems. In this paper, we present a new datafree strategy for constructing the informed subspace in which the computationally costly subspace construction can be performed in an offline phase, meaning before observing any data sets. In the subsequent online phase, the data set is observed and the precomputed informed subspace is utilized to accelerate the inversion process. This computational strategy is particularly relevant for real-time systems such as medical imaging where multiple inversions are needed.

The rest of the paper is organized as follows. To begin, we introduce the problem setting in Section 2. In Section 3, we present a new data-free likelihood-informed approach to construct the subspace. Denoting the Fisher information matrix of the likelihood function by I(x) = (∇ log L y (x))(∇ log L y (x)) L y (x)dy, this approach defines the informed subspace as the rankr dominant eigenspace of the matrix

H = I(x) π pr (x)dx, (3)
with r d. This definition makes no particular assumption on the likelihood function, so it can be applied to a wide range of measurement processes, e.g., Gaussian likelihood and Poisson likelihood. It also does not involve any particular data set y, and hence can be constructed offline. Given the informed subspace, we approximate the posterior density π y pos (x) by π y pos (x) = π y pos (x r )π pr (x ⊥ |x r ), [START_REF] Bebendorf | A note on the Poincaré inequality for convex domains[END_REF] where x r and x ⊥ denote respectively the informed and the non-informed components of x. We prove that the expected Kullback-Leibler (KL) divergence of the full posterior from its approximation is bounded as

E[D KL (π Y pos || π Y pos )] ≤ κ i>r λ i (H), (5)
where the expectation is taken over the data Y ∼ π data (y), κ being the subspace Poincaré constant of the prior [START_REF] Zahm | Gradient-based dimension reduction of multivariate vector-valued functions[END_REF][START_REF] Zahm | Certified dimension reduction in nonlinear Bayesian inverse problems[END_REF] and λ i (H) the i-th largest eigenvalue of H. This way, a problem with a fast decay in the spectrum of H yields an accurate low-dimensional posterior approximation in expectation over the data.

In Section 4, we restrict the analysis to Gaussian likelihood. In this case, we show that the vector-valued extension [START_REF] Zahm | Gradient-based dimension reduction of multivariate vector-valued functions[END_REF] of the AS method [START_REF] Paul G Constantine | Active subspace methods in theory and practice: applications to kriging surfaces[END_REF], which reduces parameter dimensions via approximating forward models, also leads to the same data-free informed subspace as that obtained using [START_REF] Andrieu | Establishing some order amongst exact approximations of MCMCs[END_REF]. We can further show that, although the likelihood-informed approach and AS employ different approximations to the posterior density, the resulting approximations share the same structure as shown in [START_REF] Bebendorf | A note on the Poincaré inequality for convex domains[END_REF] and follow the same error bound as in [START_REF] Beskos | Geometric MCMC for infinite-dimensional inverse problems[END_REF].

As suggested by [START_REF] Bebendorf | A note on the Poincaré inequality for convex domains[END_REF], the factorized form of the approximate posterior densities allows for dimension-robust sampling. One can explore the low-dimensional intractable parameter reduced posterior π y pos (x r ) using methods such as MCMC, followed by direct sampling of the high-dimensional but tractable conditional prior π pr (x ⊥ |x r ). This strategy has been previously investigated, see [START_REF] Cui | Likelihood-informed dimension reduction for nonlinear inverse problems[END_REF][START_REF] Zahm | Certified dimension reduction in nonlinear Bayesian inverse problems[END_REF] and references therein. We provide a brief summary to this existing sampling strategy in Section 5. Despite the accelerated sampling offered by the informed subspace, the resulting inference results are subject to the dimension truncation error that is bounded in [START_REF] Beskos | Geometric MCMC for infinite-dimensional inverse problems[END_REF]. In Section 6, by integrating the pseudo-marginal approach [START_REF] Andrieu | The pseudo-marginal approach for efficient Monte Carlo computations[END_REF][START_REF] Andrieu | Convergence properties of pseudo-marginal Markov chain Monte Carlo algorithms[END_REF] and the surrogate transition approach [START_REF] Andrés | Markov chain Monte Carlo using an approximation[END_REF][START_REF] Jun | Monte Carlo strategies in scientific computing[END_REF][START_REF] Jun | Sequential Monte Carlo methods for dynamic systems[END_REF] into the abovementioned sampling strategy, we present new exact inference algorithms that can enjoy the same subspace acceleration while target on the full posterior. Our exact inference algorithms only require minor modifications to the sampling strategy of [START_REF] Cui | Likelihood-informed dimension reduction for nonlinear inverse problems[END_REF][START_REF] Zahm | Certified dimension reduction in nonlinear Bayesian inverse problems[END_REF].

While our dimension reduction method readily apply for Gaussian priors, its application to non-Gaussian priors might not be straightforward. In Section 7, we show how to use the propose method for problems with Besov priors [START_REF] Dashti | Besov priors for Bayesian inverse problems[END_REF][START_REF] Ville Kolehmainen | Sparsity-promoting Bayesian inversion[END_REF][START_REF] Lassas | Discretization-invariant Bayesian inversion and Besov space priors[END_REF] which are commonly used in image reconstruction problems.

We demonstrate the accuracy of the proposed data-free LIS and the efficiency of new sampling strategies on a range of problems. These include the identification of the diffusion coefficient of a two-dimensional elliptic partial differential equation (PDE) with a Gaussian prior in Section 8 and Positron emission tomography (PET) with Poisson data and a Besov prior in Section 9.

Problem setting

For high-dimensional ill-posed inverse problems, the data are often informative only along a few directions in the parameter space. To detect and exploit this low-dimensional structure, we introduce a projector P r ∈ R d×d of rank r d such that Im(P r ) is the informed subspace and Ker(P r ) the non-informed one. This splits the parameter space as

R d = Im(P r ) ⊕ Ker(P r ),
where the subspaces Im(P r ) and Ker(P r ) are not necessarily orthogonal unless P r is orthogonal. The fact that the data are only informative in Im(P r ) means there exists an approximation to the posterior density π y pos (x) ∝ L y (x)π pr (x) under the form π y pos (x) ∝ L y (P r x)π pr (x), [START_REF] Beskos | Multilevel sequential Monte Carlo with dimension-independent likelihood-informed proposals[END_REF] in which the likelihood function x → L y (x) is replaced by a ridge function x → L y (P r x). A ridge function [START_REF] Pinkus | Ridge functions[END_REF] is a function which is constant on a subspace, here Ker(P r ). Let x r = P r x and x ⊥ = (I d -P r )x be the components of x in Im(P r ) and Ker(P r ), respectively. We have the parameter decomposition

x = x r + x ⊥ .
Using a slight abuse of notation, we factorize the prior density as π pr (x) = π pr (x r )π pr (x ⊥ |x r ), where

π pr (x r ) = Ker(P r ) π pr (x r + x ⊥ )dx ⊥ and π pr (x ⊥ |x r ) = π pr (x r + x ⊥ )/π pr (x r )
denote the marginal prior and the conditional prior. The approximate posterior [START_REF] Beskos | Multilevel sequential Monte Carlo with dimension-independent likelihood-informed proposals[END_REF] writes

π y pos (x r + x ⊥ ) ∝ L y (x r )π pr (x r ) π y pos (x r ) π pr (x ⊥ |x r ).
This factorization shows that, under the approximate posterior density, the Bayesian update is effective on the informed subspace Im(P r ) (first term π y pos (x r )), while the non-informed subspace Ker(P r ) is characterized by the prior (second term π pr (x ⊥ |x r )). This property will be exploited later on to design efficient sampling strategies for exploring both the approximate posterior and the full posterior.

The challenge of dimension reduction is to construct both the low-rank projector P r and the ridge approximation L y such that the KL divergence of the full posterior from its approximation

D KL (π y pos || π y pos ) = log π y pos (x) π y pos (x)
π y pos (x)dx, can be controlled. In this work, we specifically focus on constructing a projector P r which is independent on the data y and which allows to bound D KL (π y pos || π y pos ).

Dimension reduction via optimal parameter-reduced likelihood

In this section, we first briefly review the optimal parameter-reduced likelihood and the datadependent LIS proposed in [START_REF] Zahm | Certified dimension reduction in nonlinear Bayesian inverse problems[END_REF], and then we will introduce the data-free LIS.

Optimal parameter-reduced likelihood using a given projector

As shown in Section 2.1 of [START_REF] Zahm | Certified dimension reduction in nonlinear Bayesian inverse problems[END_REF], for a given data set y and a given projector P r , the parameterreduced likelihood function

L * ,y (x r ) = Ker(P r ) L y (x r + x ⊥ )π pr (x ⊥ |x r ) dx ⊥ , (7)
is an optimal approximation in the sense that it minimizes L y → D KL (π y pos || π y pos ). We denote by π * ,y pos (x) ∝ L * ,y (P r x)π pr (x), the resulting approximate posterior density. The marginal density π * ,y pos (x r ) = Ker(P r ) π * ,y pos (x r + x ⊥ )dx ⊥ can be expressed as π * ,y pos (x r ) ∝ L * ,y (x r )π pr (x r )

∝

Ker(P r ) L y (x r + x ⊥ )π pr (x r + x ⊥ )dx ⊥ (1)
∝ π y pos (x r ), [START_REF] Billaud-Friess | A tensor approximation method based on ideal minimal residual formulations for the solution of high-dimensional problems[END_REF] for all x r ∈ Im(P r ), where π y pos (x r ) = Ker(P r ) π y pos (x r + x ⊥ )dx ⊥ is the marginal density of the full posterior. Thus, for any projector P r and any data y, the approximate posterior π * ,y pos and the full posterior π y pos have the same marginal density on Im(P r ). In summary we have

π y pos (x) = π y pos (x r )π y pos (x ⊥ |x r ), π * ,y pos (x) (8)
= π y pos (x r )π pr (x ⊥ |x r ), which shows that the optimal approximation π * ,y pos (x) to π y pos (x) replaces the conditional posterior π y pos (x ⊥ |x r ) with the conditional prior π pr (x ⊥ |x r ).

Data-dependent dimension reduction

We denote by P r = P y r a projector built by a data-dependent approach. Ideally, we would like to build P y r that minimizes D KL (π y pos ||π * ,y pos ) over the manifold of rank-r projectors. However, this non-convex minimization problem can be challenge to solve. Instead, the strategy proposed in [START_REF] Zahm | Certified dimension reduction in nonlinear Bayesian inverse problems[END_REF] minimizes an upper bound of the KL divergence obtained by logarithmic Sobolev inequalities, in which the following assumption on the prior density is adopted.

Assumption 3.1 (Subspace logarithmic Sobolev inequality

). There exists a symmetric positive definite matrix Γ ∈ R d×d and a scalar κ > 0 such that for any projector P r ∈ R d×d and for any continuously differentiable function h : R d → R the inequality

R d h(x) 2 log h(x) 2 h P r (x) 2 π pr (x)dx ≤ 2κ R d (I d -P r ) ∇h(x) 2 Γ -1 π pr (x)dx,
holds, where h 2 P r is the conditional expectation of h 2 given by h P r (x) 2 = Ker(P r ) h(P [START_REF] Zahm | Certified dimension reduction in nonlinear Bayesian inverse problems[END_REF] gives sufficient conditions on the prior density such that Assumption 3.1 holds. In particular, any Gaussian prior π pr = N (m pr , Σ pr ) with mean m pr ∈ R d and non-singular covariance matrix Σ pr ∈ R d×d satisfies Assumption 3.1 with κ = 1 and Γ = Σ -1 pr . As shown in [54, example 2], any Gaussian mixture also satisfies this assumption, but with a constant κ which might not be accessible in practice. We refer to [START_REF] Guionnet | Lectures on logarithmic Sobolev inequalities[END_REF][START_REF] Ledoux | Logarithmic Sobolev inequalities for unbounded spin systems revisited[END_REF] for nicely written introductions to logarithmic Sobolev inequalities and examples of distributions which satisfy it. Proposition 3.2. Suppose π pr satisfies Assumption 3.1 and the likelihood function L y is continuously differentiable. Then, for any projector P r ∈ R d×d , the posterior approximation π * ,y pos (x) ∝ L * ,y (x r )π pr (x r ) π pr (x ⊥ |x r ) induced by the optimal parameter-reduced likelihood as in [START_REF] Beskos | MCMC methods for diffusion bridges[END_REF] satisfies

r x + x ⊥ ) 2 π pr (x ⊥ |P r x) dx ⊥ . Here the norm • Γ -1 is defined by v 2 Γ -1 = v Γ -1 v for any v ∈ R d . Theorem 1 in
D KL (π y pos ||π * ,y pos ) ≤ κ 2 trace Γ -1 (I d -P r )H(y)(I d -P r ) , (9)
where the matrix H(y) ∈ R d×d is defined by

H(y) = R d ∇ log L y (x))(∇ log L y (x) π y pos (x)dx. ( 10
)
Proof. See the proof of Corollary 1 in [START_REF] Zahm | Certified dimension reduction in nonlinear Bayesian inverse problems[END_REF].

The resulting projector P y r yields an approximate posterior density π * ,y pos that satisfies

D KL (π y pos ||π * ,y pos ) ≤ κ 2 d i=r+1 λ y i .
The above relation can be used to choose the rank r = rank(P y r ) to guarantee that the D KL (π y pos ||π * ,y pos ) is bounded below some user-defined tolerance. A rapid decay in the spectrum (λ y 1 , λ y 2 , . . .) ensures that one can choose a rank r that is much lower than the original dimension d. Note that the projector P y r may not be unique, unless there exists a spectral gap λ y r > λ y r+1 which ensures the r-dimensional dominant eigenspace of (H(y), Γ) is unique. Remark 3.3 (Coordinate selection). The projector defined in [START_REF] Andrés | Markov chain Monte Carlo using an approximation[END_REF] is, in general, not aligned with the canonical coordinates. However, in some parametrizations-for example, different components of x represent physical quantities of different nature-we may prefer coordinate selection than subspace identification to make the dimension reduction more interpretable. Denoting the i-th canonical basis vector of R d by e i , we let P y r = i∈I e i e i , be the projector of rank r = #I, which extracts the components of x indexed by the index set I ⊂ {1, . . . , d} such that

Im(P y r ) = span{x i : i ∈ I}, Ker(P y r ) = span{x i : i / ∈ I}.
Using such a projector, the bound (9) becomes

D KL (π y pos ||π * ,y pos ) ≤ κ 2 i / ∈I (Γ -1 ) ii H(y) ii ,
which suggests to define the index set I that selects the r largest values of (Γ -1 ) ii H(y) ii .

Because of the dependency on the data set y, the projector P y r must be built after a data set has been observed, see Algorithm 1. For scenarios where one wants to solve multiple inverse problems with multiple data sets, the matrix H(y) and the resulting projector have to be reconstructed for each data set. This can be a computationally challenging task. In addition, H(y) is defined as an expectation over the high-dimensional posterior distribution, which further raises the computational burden.

Algorithm 1: Data-dependent dimension reduction. Requires: π pr satisfying Assumption 3.1, tolerance ε > 0 and maximal rank r max Online phase: given the data y do:

Compute the matrix H(y) using [START_REF] Carvalho | Handling sparsity via the horseshoe[END_REF].

Compute the generalized eigendecomposition H(y)v y i = λ y i Γv y i . Find the smallest r such that κ 2 d i=r+1 λ y i ≤ ε. If r ≥ r max , set r = r max . Assemble the projector P y r using [START_REF] Andrés | Markov chain Monte Carlo using an approximation[END_REF]. Define the conditional expectation L * ,y (x r ) defined in [START_REF] Beskos | MCMC methods for diffusion bridges[END_REF]. Return: Approximate posterior π * ,y pos (x) ∝ L * ,y (P y r x)π pr (x).

Data-free dimension reduction

To overcome the abovementioned computational burden of recomputing the data-dependent projector for every new data set, we present a new data-free dimension reduction method. The key idea is to control the KL divergence in expectation over the marginal density of data. We introduce an m-dimensional random vector

Y ∼ π data (y),
where π data is the marginal density of data defined in [START_REF] Andrieu | Convergence properties of pseudo-marginal Markov chain Monte Carlo algorithms[END_REF]. Note that the observed data y corresponds to a particular realization of Y . For a given projector P r independent on the data, replacing y with Y in [START_REF] Brooks | Handbook of Markov Chain Monte Carlo[END_REF] and taking the expectation over Y yields

E D KL (π Y pos ||π * ,Y pos ) ≤ κ 2 trace Γ -1 (I d -P r )E[H(Y )](I d -P r ) . (12)
Here, the approximate posterior π * ,Y pos depends on Y via the optimal likelihood L * ,Y . Similar to the data-dependent case, the leading generalized eigenvectors of the matrix pair (E[H(Y )], Γ) can be used to obtain a projector that minimizes the error bound. However, in this case, the matrix E[H(Y )] is the expectation of H(y) over the marginal density of data, and thus it is independent of observed data. andλ y i ≥ λ y j for all i ≤ j. The data-free projector P r that minimizes the right-hand side of ( 12) is given by Im(P r ) = span{v 1 , . . . , v r },

Let (λ i , v i ) ∈ R ≥0 × R d denotes the i-th eigenpair of (E[H(Y )], Γ) such that E[H(Y )]v i = λ i Γv i , with v i Γv j = δ i,j
Ker(P r ) = span{v r+1 , . . . , v d }. (13)
When using this projector for defining the approximate posterior π * ,Y pos , the expectation of the KL divergence D KL (π Y pos ||π * ,Y pos ) can be controlled as

E D KL (π Y pos ||π * ,Y pos ) ≤ κ 2 d i=r+1 λ i . ( 14
)
Remark 3.4 (Bound in high probability). Inequality [START_REF] Cortesi | Forward and backward uncertainty quantification with active subspaces: application to hypersonic flows around a cylinder[END_REF] gives a bound on D KL (π Y pos ||π * ,Y pos ) in expectation. In order to obtain a bound in high probability, let us use the Markov inequality

P{D KL (π Y pos ||π * ,Y pos ) ≤ ε} ≥ 1 -ε -1 E[D KL (π Y pos ||π * ,Y
pos )] for some ε > 0. Thus, for a given 0 < η ≤ 1, the condition κ

2 d i=r+1 λ i ≤ ε η is sufficient to ensure that D KL (π Y pos ||π * ,Y pos ) ≤ ε,
holds with a probability greater than 1 -η. Remark 3.5 (Coordinate selection). Similarly to Remark 3.3, instead of defining P r as in [START_REF] Paul G Constantine | Accelerating Markov chain Monte Carlo with active subspaces[END_REF], we can define a coordinate-aligned projector P r = i∈I e i e i by selecting an index set I corresponding to the r largest values of (Γ -1 ) ii E[H(Y )] ii . Now we show that the matrix E[H(Y )] admits a simple expression in terms of the Fisher information matrix associated with the likelihood function. This leads to a computationally convenient way to construct the data-free projector. Recall that the likelihood L y (x), seen as a function of y, is the pdf of the data y conditionned on the parameter x ∈ R d . The Fisher information matrix associated with this family of pdf is

I(x) = R m ∇ log L y (x))(∇ log L y (x) L y (x) dy. ( 15
)
We can write

E[H(Y )] = R m H(y)π data (y)dy (10) = R m R d ∇ log L y (x))(∇ log L y (x) π y pos (x)dx π data (y)dy (1) = R m ×R d ∇ log L y (x))(∇ log L y (x) L y (x)π pr (x) R d L y (x )π pr (x )dx π data (y)dxdy (2) = R m ×R d ∇ log L y (x))(∇ log L y (x) L y (x)π pr (x)dxdy (15) = R d I(x)π pr (x)dx, (16)
which shows that the matrix E[H(Y )] is the expectation of the Fisher information matrix over the prior. This expression does not involve any expectation over the posterior density, which is a major advantage compared to the expression [START_REF] Carvalho | Handling sparsity via the horseshoe[END_REF] of the data-dependent matrix H(y). The methodology presented here is summarized in Algorithm 2.

Algorithm 2: Data-free dimension reduction Requires: π pr satisfying Assumption 3.1, Fisher information matrix I(x) of L y , tolerance ε > 0, and maximal rank r max Offline phase

Compute the matrix

H I = R d I(x)π pr (x)dx. Compute the generalized eigendecomposition H I v i = λ i Γv i . Find the smallest r such that κ 2 d i=r+1 λ i ≤ ε. If r ≥ r max , set r = r max .
Return: Projector P r defined by [START_REF] Paul G Constantine | Accelerating Markov chain Monte Carlo with active subspaces[END_REF] Online phase: given the data y do:

Define L * ,y as the conditional expectation defined in [START_REF] Beskos | MCMC methods for diffusion bridges[END_REF]. Return: Approximate posterior π * ,y pos (x) ∝ L * ,y (P r x)π pr (x)

Example 3.6 (Gaussian likelihood). Consider the parameter-to-data map is represented by a smooth forward model G : R d → R m and corrupted by an additive Gaussian noise ξ obs ∼ N (0, Σ obs ) with non-singular covariance matrix Σ obs ∈ R m×m , i.e., y = G(x) + ξ obs , where ξ obs ∼ N (0, Σ obs ).

The likelihood function takes the form

L y (x) = Z -1 exp(-1 2 G(x) -y 2 Σ - 1 obs
), where Z = (2π) m det(Σ obs ) is a normalizing constant. The Slepian-Bangs formula gives an explicit expression for the Fisher information matrix

I(x) = ∇G(x) Σ -1
obs ∇G(x), where ∇G(x) ∈ R m×d denotes the Jacobian of the forward model G(x). By relation [START_REF] Cui | Multilevel dimension-independent likelihoodinformed MCMC for large-scale inverse problems[END_REF] we obtain

E[H(Y )] = R d ∇G(x) Σ -1 obs ∇G(x) π pr (x)dx. (17)
A similar matrix was considered in [START_REF] Cui | Likelihood-informed dimension reduction for nonlinear inverse problems[END_REF] in the context of data-dependent dimension reduction. The major difference with [START_REF] Cui | Dimension-independent likelihood-informed MCMC[END_REF] is that, in [START_REF] Cui | Likelihood-informed dimension reduction for nonlinear inverse problems[END_REF], the expectation is taken over the posterior density rather than over the prior.

Dimension reduction via parameter-reduced forward model

In the previous Section 3, the detection of the data-free informed subspace is based on an approximation of the likelihood function. In this section, we present an alternative strategy which, under Gaussian likelihood assumption, consist in approximating the forward model instead of the likelihood itself. This approach is similar to the vector-valued extension of the AS method [START_REF] Zahm | Gradient-based dimension reduction of multivariate vector-valued functions[END_REF] and still yields error bounds for the expected KL divergence.

As in Example 3.6, let us start with a Gaussian likelihood of the form

L y (x) = 1 Z exp - 1 2 G(x) -y 2 Σ -1 obs , (18)
where x → G(x) is a continuously differentiable forward model, Σ obs ∈ R m×m is a non-singular covariance matrix and Z = (2π) m det(Σ obs ) a normalizing constant. Our goal is to build a low-dimensional approximation to the likelihood [START_REF] Cui | Likelihood-informed dimension reduction for nonlinear inverse problems[END_REF] by replacing the forward model with a ridge approximation x → G(P r x). That is, we look for a likelihood approximation of the form

L y (P r x) = 1 Z exp - 1 2 G(P r x) -y 2 Σ -1 obs , (19)
where P r is a low-rank projector and where G is some parameter-reduced function defined over Ker(P r ). In general, this approximate likelihood [START_REF] Cui | A unified performance analysis of likelihood-informed subspace methods[END_REF] is different than the previous one L * ,y , see [START_REF] Beskos | MCMC methods for diffusion bridges[END_REF], and therefore L y might not be optimal with respect to the KL divergence as discussed in Section 3.1. The following proposition will guide the construction of the approximate forward model.

Proposition 4.1. Consider the posterior density π y pos (x) ∝ L y (x)π pr (x) with a Gaussian likelihood as in [START_REF] Cui | Likelihood-informed dimension reduction for nonlinear inverse problems[END_REF]. For any approximate forward model G : R d → R m , the resulting approximate likelihood

L y (x) = 1 Z exp(-1 2 G(x) -y 2 Σ - 1 obs
) defines an approximate posterior density π y pos (x) ∝ L y (x)π pr (x) such that

E D KL (π Y pos || π Y pos ) + D KL (π data || π data ) = 1 2 R d G(x) -G(x) 2 Σ -1 obs π pr (x)dx.
Here the expectation is taken over Y ∼ π data and π data (y) = R d L y (x)π pr (x) dx is the approximate marginal density of data.

Proof. See Appendix A.

Using an approximate forward model in the form of G(x) = G(P r x), Proposition 4.1 ensures that the approximate posterior π y pos (x) ∝ L y (P r x)π pr (x) with L y (P r x) as in ( 19) satisfies

E[D KL (π Y pos || π Y pos )] ≤ 1 2 R d G(x) -G(P r x) 2 Σ -1 obs π pr (x)dx, (20)
This suggests to construct a ridge approximation G(P r x) to G(x) in the L 2 π pr sense. To accomplish this, we follow the methodology proposed in [START_REF] Zahm | Gradient-based dimension reduction of multivariate vector-valued functions[END_REF] for the approximation of multivariate function using gradient information. First, for any projector P r , the optimal function G * that minimizes the right-hand side of ( 20) is the conditional expectation

G * (x r ) = Ker(P r ) G(x r + x ⊥ )π pr (x ⊥ |x r ) dx ⊥ . (21)
Then, similarly to Assumption 3.1, we assume that π pr satisfies the following subspace Poincaré inequality.

Assumption 4.2 (Subspace Poincaré inequality

). There exists a symmetric positive definite matrix Γ ∈ R d×d and a scalar κ > 0 such that for any projector P r ∈ R d×d and for any continuously differentiable function h : R d → R the inequality

R d h(x) -h P r (x) 2 π pr (x)dx ≤ κ R d (I d -P r ) ∇h(x) 2 Γ -1 π pr (x)dx,
holds, where h P r is the conditional expectation of h defined by h P r (x) = Ker(P r ) h(P r x +

x ⊥ )π pr (x ⊥ |P r x)dx ⊥ .
Assumption 4.2 is weaker than Assumption 3.1, in the sense that any distribution which satisfies the subspace logarithmic Sobolev inequality automatically satisfies the subspace Poincaré inequality with the same κ and the same Γ, see for instance [START_REF] Zahm | Certified dimension reduction in nonlinear Bayesian inverse problems[END_REF]Corollary 2]. We refer to the recent contributions [START_REF] Bebendorf | A note on the Poincaré inequality for convex domains[END_REF][START_REF] Teixeira Parente | Generalized bounds for active subspaces[END_REF][START_REF] Roustant | Poincaré inequalities on intervals-application to sensitivity analysis[END_REF] for examples of probability distribution which satisfy (subspace) Poincaré inequality. As for the logarithmic-Sobolev constant, the Poincaré constant is hard to compute in practice, except the case of Gaussian prior. Using similar arguments as in the proof of Proposition 2.5 in [START_REF] Zahm | Gradient-based dimension reduction of multivariate vector-valued functions[END_REF], Assumption 4.2 allows to write

R d G(x) -G * (P r x) 2 Γ -1 obs π pr (x)dx ≤ κ trace Γ -1 (I d -P r )H G (I d -P r ) , (22)
holds for any projector P r , where the matrix H G ∈ R d×d is defined by

H G = R d ∇G(x) Σ -1 obs ∇G(x) π pr (x)dx, (23)
with ∇G(x) the Jacobian matrix of G(x) given by

∇G(x) =     ∂G 1 ∂x 1 (x) . . . ∂G 1 ∂x d (x) . . . . . . . . . ∂G m ∂x 1 (x) . . . ∂G m ∂x d (x)     .
Again, the projector P G r that minimizes the right-hand side of ( 22) can be constructed via the generalized eigenvalue problem

H G v G i = λ G i Γv G i : Im(P G r ) = span{v G 1 , . . . , v G r }, Ker(P G r ) = span{v G r+1 , . . . , v G d }. (24)
Using this projector to construct the approximate forward model G * in ( 21) and the approximate likelihood as in [START_REF] Cui | A unified performance analysis of likelihood-informed subspace methods[END_REF], Proposition 4.1 and the inequality in ( 22) yield

E D KL (π Y pos || π Y pos ) ≤ κ 2 d i=r+1 λ G i . ( 25
)
The methodology is summarized in Algorithm 3.

The matrix H G used in this case takes the same form as the matrix E[H(Y )] in Section 3.3 with the Gaussian likelihood (cf. Example 3.6), and hence results in the same data-free projector. However, the resulting approximate likelihood functions are not the same. Indeed in Section 3.3, the optimal approximate likelihood L * ,y is given as the conditional expectation of the likelihood function (cf. ( 7)), whereas here, L y is defined by the conditional expectation of the forward model G * (cf. ( 21)). Using either the parameter-reduced likelihood in [START_REF] Beskos | MCMC methods for diffusion bridges[END_REF] or the parameter-reduced forward model in [START_REF] Dashti | Besov priors for Bayesian inverse problems[END_REF] results in the same parameter truncation error bound in terms of expected KL divergence.

Algorithm 3: Data-free dimension reduction via forward model approximation Requires: π pr satisfying Assumption 3.1, Jacobian ∇G(x), tolerance ε > 0 and maximal rank r = r max .

Offline phase

Compute the matrix H G defined in [START_REF] Doucet | Efficient implementation of Markov chain Monte Carlo when using an unbiased likelihood estimator[END_REF] Compute the generalized eigendecomposition

H G v G i = λ G i Γv G i Find the smallest r such that κ 2 d i=r+1 λ i ≤ ε. If r ≥ r max , set r = r max Assemble the projector P G
r defined in (24) Define G * as the conditional expectation [START_REF] Dashti | Besov priors for Bayesian inverse problems[END_REF] Return: Approximate forward model x → G * (P G r x) Online phase: given the data y do:

Assemble L y (P G r x) as in ( 19) Return: Approximate posterior π y pos (x) ∝ L y (P r x)π pr (x)

Remark 4.3. Despite the similarity between the approximate likelihood functions given in ( 7) and ( 19), these two approaches offer different computational characteristics. Given the data-free informed subspace, the optimal parameter-reduced forward model x r → G * (x r ) can be further replaced by a surrogate model x r → G ROM (x r ) constructed in the offline phase. The surrogate model can be obtained using tensor methods [START_REF] Billaud-Friess | A tensor approximation method based on ideal minimal residual formulations for the solution of high-dimensional problems[END_REF][START_REF] Nouy | Low-Rank Tensor Methods for Model Order Reduction[END_REF], the reduced basis method [START_REF] Anthony T Patera | Reduced basis approximation and a posteriori error estimation for parametrized partial differential equations[END_REF], polynomial techniques [START_REF] Le | Spectral methods for uncertainty quantification: with applications to computational fluid dynamics[END_REF], etc., just to cite a few. All these approximation techniques do not scale well with the apparent parameter dimensions d, and thus parameter reduction can greatly improve the scalability of surrogate models.

In contrast, the conditional expectation of the likelihood function in ( 7) cannot be replaced with offline surrogate models because of the data-dependency of the likelihood.

Sampling the approximate posterior

Given a data-free informed subspace, the approximate posterior density has the factorized form

π y pos (x) ∝ π y post (x r )π pr (x ⊥ |x r ), (26)
with either π y post (x r ) = π y post (x r ) in the optimal parameter-reduced likelihood approach of Section 3, or with π y post (

x r ) = L y (x r )π y pr (x r ), L y (x r ) ∝ exp(-1 2 G * (x r )-y 2 Σ - 1 obs
) in the optimal parameterreduced forward model approach of Section 4. The factorization (26) naturally suggests a dimension robust way to sampling the approximate posterior. The sampling method consists in first drawing samples x (1) r , x (2) r , ..., x (K) r from the low-dimensional density π y pos (x r ) using either MCMC or SMC method. Then, for each sample x (j) r , we simulate a conditional prior sample x

(j) ⊥ from π pr (x ⊥ |x (j) r ). In the end, x (j) = x (j) r + x (j)
⊥ are samples from the approximate posterior π y pos (x). We emphasis here that the key is to be able to sample from the conditional prior π pr (x ⊥ |x r ). This task is rather easy for Gaussian priors. We show in Section 7 how to sample from π pr (x ⊥ |x r ) for non-Gaussian priors with a particular structure that can be exploited. Remark 5.1. If the end goal is to compute expectation of some function h over of the approximate posterior, the factorization (26) leads to

R d h(x) π pos (x)dx = Im(P r ) Ker(P r ) h(x r + x ⊥ )π pr (x ⊥ |x r )dx ⊥ π pos (x r )dx r ≈ 1 K K j=1 Ker(P r ) h(x (j) r + x ⊥ )π pr (x ⊥ |x (j) r )dx ⊥ ,
where x (1) r , ..., x (K) r are samples from the approximate marginal posterior π y pos (x r ). This way, if the expectation over the conditional prior π pr (x ⊥ |x (j) r ) can be carried out analytically, one can can simply avoid using conditional prior samples. Alternatively, the K conditional expectations h(x (j) r + x ⊥ )π pr (x ⊥ |x (j) r )dx ⊥ can also be approximated via other accurate quadrature rule for π pr (x ⊥ |x (j) r ). Either way, we assume that integration with respect to the conditional prior is tractable.

In Algorithm 4 we provide the details of an MCMC-based sampling procedure in which the approximate likelihood (defined by either optimal parameter-reduced likelihood or optimal parameter-reduced forward model) can be obtained as sample averages over the conditional prior π pr (x ⊥ |x r ). To make these approximations generally applicable, we replace the conditional prior with the marginal prior π pr (x ⊥ ) in computing those conditional expectations in the Equations ( 27) and ( 28) in Algorithm 4. Note that the typical class of inverse problems equipped with a Gaussian prior π pr = N (m pr , Σ pr ) is a special case. Since the projector P r is orthogonal with respect to Σ -1 pr , the marginal prior π pr (x ⊥ ) coincides with the conditional prior π pr (x ⊥ |x r ).

A remaining question is how to choose the sample size N for computing the conditional expectations in ( 27) and [START_REF] Guionnet | Lectures on logarithmic Sobolev inequalities[END_REF]. The following heuristic is developed based on the optimal parameterreduced forward model. Consider the exact parameter-reduced forward model G * (P r x) = G * (x r ) and its sample-averaged approximation G N (P r x) = G N (x r ). The sample-averaged approximation defines an approximate posterior density

π y pos (x) ∝ exp - 1 2 G N (P r x) -y 2 Σ -1 obs π pr (x),
that satisfies

E[D KL (π Y pos || π Y pos )] (20)
≤ E 1 2 R d G(x) -G N (P r x) 2 Σ -1 obs π pr (x)dx = 1 2 E R d G(x) -G * (P r x) 2 Σ -1 obs + G * (P r x) -G N (P r x) 2 Σ -1 obs π pr (x)dx = 1 + 1 N 1 2 R d G(x) -G * (P r x) 2 Σ -1 obs π pr (x)dx. ( 30
)
Here, the expectation is taken jointly over the data Y and the sample {x

(i) ⊥ } N i=1 .
The above inequality directly follows from Proposition 4.1 and the fact that G * (P r x) is the conditional expectation of G(x) over Ker(P r ). We refer to Theorem 3.2 in [START_REF] Paul G Constantine | Active subspace methods in theory and practice: applications to kriging surfaces[END_REF] for more details on this derivation. Inequality [START_REF] Heikkinen | Statistical inversion theory in x-ray tomography[END_REF] implies that the random approximate posterior π y pos (x) can be used in place of π y pos (x), as the bounds on the expected Kullback-Leibler divergence in [START_REF] Arnak S Dalalyan | Theoretical guarantees for approximate sampling from smooth and log-concave densities[END_REF] and (30) are comparable. In addition, this suggests that the sample size N in [START_REF] Guionnet | Lectures on logarithmic Sobolev inequalities[END_REF] does not have to be large.

Algorithm 4: MCMC-based approach for sampling the approximate posterior. Requires: A projector P r , a sample size N for approximating the likelihood, a total posterior sample size K, and a proposal density q(•|x r ) on Im(P r ).

Sample-averaged likelihood approximation

Draw N i.i.d. samples x

⊥ , . . . , x

(N ) ⊥ from the marginal π pr (x ⊥ ) if optimal parameter-reduced likelihood is used then

L y N (x r ) = 1 N N i=1 L y (x r + x (i) ⊥ ) (27) for any x r ∈ Im(P r ) end if optimal parameter-reduced forward model is used then L y N (x r ) ∝ exp - 1 2 G N (x r ) -y 2 Σ -1 obs , G N (x r ) = 1 N N i=1 G(x r + x (i) ⊥ ) (28)
for any x r ∈ Im(P r ) end Return: a sample-averaged approximate likelihood function L y N (x r ) Subspace MCMC sampling for j = 1, 2, ..., K do Given the Markov chain state X (j-1)

r = x r , propose a candidate x † r ∼ q(•|x r ) Evaluate the approximate likelihood function L y N (x † r ) Compute the acceptance probability α(x † r |x r ) = min 1 , L y N (x † r )π pr (x † r ) q(x r |x † r ) L y N (x r )π pr (x r ) q(x † r |x r ) . (29)
With probability α(x † r |x r ), accept x † r by setting X (j) r = x † r , otherwise reject x † r by setting X (j) r = x r . end Return: a Markov chain X (1) r , X (2) r , ..., X (K) r with invariant density π y pos (x r ) Approximate posterior sampling for j = 1, 2, ..., K do Given the state X (j) r = x (j) r , draw a conditional prior sample x

⊥ ∼ π pr (•|x (j) r ) Compute the i-th approximate posterior sample x (j) = x (j) r + x (j) ⊥ end Return: approximate marginal posterior samples x (1) , x (2) , ..., x (K) Even with N = 1, ( 20) and ( 30) differs only by a factor of 2. For the optimal parameter-reduced likelihood function, it is not obvious how to obtain a similar bound for the sampled-averaged conditional expectation in [START_REF] Ivan G Graham | Quasi-Monte Carlo methods for elliptic PDEs with random coefficients and applications[END_REF], see for instance the result [54, Proposition 5]. In this case, we adopt the identity (30) as a heuristic.

Sampling from the exact posterior

In this section, we present new strategies for sampling the exact posterior by adding minor modifications to Algorithm 4.

6.1. Pseudo-marginal for the optimal parameter-reduced likelihood For the optimal parameter-reduced likelihood approach, Algorithm 4 replaces the optimal likelihood L * ,y (x r ) with the sample-average L y N (x r ) defined by ( 27) using frozen (fixed) samples {x

(i) ⊥ } N i=1
. This way, Algorithm 4 produces samples from an estimation to the posterior approximation π * ,y pos (x) = π pos (x r )π pr (x ⊥ |x r ). In this section, we first show that replacing the frozen samples with freshly drawing samples {x

(i) ⊥ } N
i=1 at each MCMC iteration yields a pseudo-marginal MCMC [START_REF] Andrieu | The pseudo-marginal approach for efficient Monte Carlo computations[END_REF] which samples exactly from π * ,y pos (x). In addition, we also show that an appropriate recycling of the data generated by this modified algorithm allows obtaining samples from the exact posterior π y pos (x) itself. We propose to modify Algorithm 4 by replacing the acceptance rate α N (x † r |x r ) in ( 29) with

α N (x † r |x r ) = min    1 , π pr (x † r ) 1 N N i=1 L y (x † r + x †(i) ⊥ ) q(x r |x † r ) π pr (x r ) 1 N N i=1 L y (x r + x (i) ⊥ ) q(x † r |x r )    . (31)
Here, {x

⊥ } N i=1 = {x †(i) ⊥ } N i=1
where frozen, the new acceptance rate (31) requires to redraw fresh samples at each proposal candidate x † r . This is summarized in Algorithm 5.

Algorithm 5: Pseudo-marginal MCMC for sampling the exact marginal posterior. Requires: A projector P r , a sample size N for approximating the likelihood, a total posterior sample size K, and a proposal density q(•|x r ) on Im(P r ). for j = 1, 2, ..., K do Given the previous state of the Markov chain X (j-1) r = x r and the associated set of conditional prior samples {X

(j-1,i) ⊥ } N i=1 = {x (i) ⊥ } N i=1 Propose a candidate x † r ∼ q(•|x r ) Draw N independent samples x †(1) ⊥ , ..., x †(N ) ⊥ ∼ π pr (x ⊥ |x † r ) Compute the acceptance probability α(x † r |x r ) as in (31) With probability α(x † r |x r ), accept X (j) r = x † r and {X (j,i) ⊥ } N i=1 = {x ( †,i) ⊥ } N i=1 .
Otherwise reject and set X (j) r = x r and {X (j,i)

⊥ } N i=1 = {x (i) ⊥ } N i=1 . end Return: the Markov chain {(X (j) r , {X (j,i) ⊥ } N i=1 )} K j=1
In the next proposition we apply the analysis of pseudo-marginal MCMC [START_REF] Andrieu | The pseudo-marginal approach for efficient Monte Carlo computations[END_REF] to show that π y pos (x r ) is the invariant density of the Markov chain constructed by Algorithm 5. The key step is to interpret Algorithm 5 as a classical Metropolis-Hastings algorithm that operates on the product space Im(P r ) × Ker(P r ) N . Proposition 6.1. Algorithm 5 constructs an ergodic Markov chain {(X (j) r , {X (j,i)

⊥ } N i=1
)} j≥1 on the product space Im(P r ) × Ker(P r ) N with invariant density π y,N tar (x r , {x

(i) ⊥ } N i=1 ) ∝ π pr (x r ) N i=1 L y (x r + x (i) ⊥ ) N j=1 π pr (x (j) ⊥ |x r ). ( 32
)
The marginal of this target density satisfies π y,N tar (x r ) = π y pos (x r ) so that the sequence {X (j) r } N j=1 is an ergodic Markov chain with the invariant density π y pos (x r ).

Proof. See Appendix B.

Remark 6.2 (Choosing N in Algorithm 5). The statistical performance of pseudo-marginal methods depends on the variance of the sample-averaged estimate 1

N N i=1 L y (x r +x (i)
⊥ ). This variance being inversely proportional to the sample size N , a larger N may result in better statistical efficiency of the MCMC chain. However, the computational cost per MCMC iteration increases linearly with N , while the improvement of the statistical efficiency will not follow the same rate. We refer the readers to [START_REF] Andrieu | Establishing some order amongst exact approximations of MCMCs[END_REF][START_REF] Doucet | Efficient implementation of Markov chain Monte Carlo when using an unbiased likelihood estimator[END_REF] for a detailed discussion on this topic and only provide an interpretation as follows. With N → ∞, the Markov chain constructed by the pseudo-marginal MCMC converges to that of an idealized standard MCMC, which has the acceptance probability defined by the same proposal density and the exact evaluation of L * ,y (x r ). This way, even with a very large N , the statistical efficiency of the pseudo-marginal MCMC cannot be improved further beyond that of the idealized standard MCMC. As suggested by [START_REF] Doucet | Efficient implementation of Markov chain Monte Carlo when using an unbiased likelihood estimator[END_REF], the standard deviation of the logarithm of the parameter-reduced likelihood estimate, var[log L y N ]

1 2 , can be used to monitor the quality of the sample-averaged estimator.

It is remarkable to observe that, for N = 1, the target density (32) becomes the true posterior π y,1 tar (x r , x ⊥ ) = π y pos (x r + x ⊥ ). This means that Algorithm 5 actually produces samples x = x r + x ⊥ from π y pos (x). For N > 1, we propose to recycle the Markov chain {X

(1,i)

⊥ } N i=1 , . . . , {X (K,i) ⊥ } N i=1
produced by Algorithm 5 in order to generate samples from the exact posterior π y pos (x). This procedure is summarized in Algorithm 6 and a justification is provided in the following proposition.

Proposition 6.3. Let {(X (j) r , {X (j,i)
⊥ } N i=1 )} j≥1 be a Markov chain generated by Algorithm 5. For any j ≥ 1 we randomly select X (j)

⊥ ∈ {X (j,i) ⊥ } N i=1 according to the probability P X (j) ⊥ = X (j,k) ⊥ X (j) r , {X (j,i) ⊥ } N i=1 = L y (X (j) r + X (j,k) ⊥ ) N i=1 L y (X (j) r + X (j,i) ⊥ ) , 1 ≤ k ≤ N, (33)
and we let X (j) = X (j) r + X

⊥ . Then {X (j) } j≥1 is a Markov chain with the exact posterior π y pos (x) as invariant density.

Proof. See Appendix C

Algorithm 6: Recycling the Markov chain generated by Algorithm 5 to generate exact posterior samples Requires: MCMC chain generated {(X (j) r , {X (j,i)

⊥ } N i=1 )} K j=1 by Algorithm 5 for j = 1, 2, ..., K do Subsample X (j) ⊥ ∈ {X (j,i)
⊥ } N i=1 according to the probability (33) Assemble X (j) = X (j) r + X (j) ⊥ end Return: the Markov chain {X (j) } K j=1 with invariant density π y pos (x)

Delayed acceptance for the optimal parameter-reduced forward model

For the optimal parameter-reduced forward model, the marginal density of the resulting approximate posterior does not coincide with that of the exact posterior in general. However, we can still modify the approximate inference algorithm 4 using the delayed acceptance technique [START_REF] Andrés | Markov chain Monte Carlo using an approximation[END_REF][START_REF] Jun | Monte Carlo strategies in scientific computing[END_REF][START_REF] Jun | Sequential Monte Carlo methods for dynamic systems[END_REF] to explore the exact posterior. The delayed acceptance modifies Algorithm 4 by adding a second stage acceptance rejection within each MCMC iteration. Here we consider the sampleaveraged likelihood L y N (x r ) defined by either ( 27) (the optimal parameter-reduced likelihood) or (28) (the optimal parameter-reduced forward model), where the marginal prior sample set {x

(i) ⊥ } N i=i is prescribed.
The following Proposition and Algorithm 7 detail this modification. Proposition 6.4. Suppose we have a proposal distribution q(•|x r ) defined in the parameter reduced subspace Im(P r ). We consider the following two stage Metropolis-Hastings method. In the first stage, we draw a proposal candidate x † r ∼ q(•|x r ). Then, with the probability

α(x † r |x r ) = min 1 , L y N (x † r )π pr (x † r ) q(x r |x † r ) L y N (x r )π pr (x r ) q(x † r |x r ) , (34)
we move the proposal candidate x † r to the next stage. In the second stage, we draw a proposal candidate π pr (x † ⊥ |x † r ) in the complement subspace Ker(P r ) and then accept the pair of proposal candidates (x † r , x † ⊥ ) with the probability

β(x † r , x † ⊥ |x r , x ⊥ ) = min 1, L y (x † r + x † ⊥ ) L y N (x r ) L y (x r + x ⊥ ) L y N (x † r ) . (35)
Then, the above procedure constructs an ergodic Markov chain with the full posterior π y pos (x) as the invariant density.

Proof. This result can be derived from the standard delayed acceptance [START_REF] Andrés | Markov chain Monte Carlo using an approximation[END_REF]. For completeness, we provide the proof in Appendix D. Remark 6.5. It worth to note that the delayed acceptance also opens the door to further accelerate the exact inference using surrogate models instead of the original forward model. The approximate likelihood L y N (x r ) is deterministic and dimension reduced, which makes it possible to further approximate L y N (x r ) using computationally fast surrogate models. In this case, the same delayed acceptance MCMC (Algorithm 7) can still produce ergodic Markov chains that converge to the Algorithm 7: Delayed acceptance MCMC for sampling the exact posterior. Requires: A projector P r , a sample-averaged likelihood approximation defined in Algorithm 4, a total sample size K, and a proposal density q(•|x r ) on Im(P r ). for j = 1, 2, ..., K do Given the Markov chain state X (j-1) = x r + x ⊥ , propose a candidate x † r ∼ q(•|x r ) Compute the parameter-reduced likelihood L y N (x † r ) using either using ( 27) or ( 28) With probability α(x † r |x r ) in ( 34) move x † r to the next stage as follows Propose a candidate

x † ⊥ ∼ π pr (•|x † r ) Compute the full likelihood L y (x † r + x † ⊥ ) With probability β(x † r , x † ⊥ |x r , x ⊥ ) in (35) accept (x † r , x † ), otherwise reject (x † r , x † ) Otherwise reject x † r end Accept: set X j = x † r + x † ⊥ Reject: set X j = X (j-1)
Return: a Markov chain X (1) , X (2) , ..., X (K) with the invariant density π y pos (x) full posterior π y pos (x). In contrast, the pseudo-marginal method requires an unbiased Monte Carlo estimate of the exact marginal posterior π y pos (x r ) at every iteration, which is not straightforward to accelerate using surrogate models.

Non-Gaussian priors

The dimension reduction techniques presented in Sections 3 and 4 require one to evaluate the marginal prior density π pr (x r ) = Ker(P r ) π pr (x r + x ⊥ )dx ⊥ and draw samples from the conditional prior π pr (x ⊥ |x r ) = π pr (x r + x ⊥ )/π pr (x r ). While these tasks are readily doable for Gaussian distributions, it might not be the case in general. In this section, we use Besov priors as an example to present strategies that can extend the proposed dimension reduction methods to some non-Gaussian priors.

Besov priors

Besov measure [START_REF] Dashti | Besov priors for Bayesian inverse problems[END_REF][START_REF] Ville Kolehmainen | Sparsity-promoting Bayesian inversion[END_REF][START_REF] Lassas | Discretization-invariant Bayesian inversion and Besov space priors[END_REF] naturally appears in image reconstruction problems in which the detection of edges and interfaces is important. Following [START_REF] Ville Kolehmainen | Sparsity-promoting Bayesian inversion[END_REF][START_REF] Lassas | Discretization-invariant Bayesian inversion and Besov space priors[END_REF], we construct Besov priors using wavelet functions and consider functions on the one-dimensional torus T = (0, 1]. Starting with a suitable compactly supported mother wavelet function ψ * ∈ L 2 (T), we can define an orthogonal basis

ψ j,k (s) = 2 j 2 ψ * (2 j s -k), j, k ∈ N ≥0 , k ∈ [0, 2 j -1].
This way, given a smoothness parameter r > 0 and integrability parameters 1 ≤ p, q ≤ ∞, a function f : s → f (s) in the Besov space B r pq (T) can be written as

f (s) = c 0 + ∞ j=0 2 j -1 k=0 2 -j(r+ 1 2 -1 p ) b j,k ψ j,k (s), (36)
and satisfies

f B r pq := |c 0 | q + ∞ j=0 2 j -1 k=0 |b j,k | p q p 1 q < ∞.
In a Bayesian setting, we can set p = q and define the Besov-B r pp prior with the pdf ‡

π pr (f ) ∝ exp -γ f p B r pp , (37)
where γ > 0 is a scale parameter. One can easily generalize the above definition of Besov priors to multidimensional settings by taking tensor products of the one-dimensional basis and associated coefficients.

We can discretize the Besov prior by truncating the infinite sum in [START_REF] Ledoux | Logarithmic Sobolev inequalities for unbounded spin systems revisited[END_REF] to the first D terms. This way, collecting all the coefficients into a parameter vector x = (c 0 , b 0,0 , b 0,1 , ..., b 1,0 , b 1,1 , ...) ∈ R d , where d = 2 D+1 , the discretized Besov-B r pp prior can be equivalently expressed as a productform distribution over the parameter x with the pdf

π pr (x) = d i=1 π (i) pr (x i ) with π (i) pr (x i ) ∝ exp (-γ|x i | p ) . (38)

Dimension reduction via coordinate selection

In general, we do not have closed form expressions for both the marginal π pr (x r ) and the conditional π pr (x ⊥ |x r ), unless the projector P r is aligned with the canonical basis. This leads to the construction of reduced subspace by selecting a subset of canonical basis. As discussed in Remarks 3.3 and 3.5, this task can be achieved by identifying an index set I ⊂ {1, . . . , d} with cardinality r such that I contains the indices of the r largest values of i

→ (Γ -1 ) ii E[H(Y )] ii in the data-dependent case or those of i → (Γ -1 ) ii E[H(Y )]
ii in the data-free case. This leads to the projector P r = i∈I e i e i , where {e 1 , . . . , e d } is the canonical basis of R d . Thus, the product-form of (38) yields the marginal prior and the conditional prior

π pr (x r ) = i∈I π (i) pr (x i ) and π pr (x ⊥ |x r ) = i / ∈I π (i) pr (x i ),
respectively. In this formulation, evaluating the marginal prior density and drawing samples from the conditional prior become straightforward tasks. Remark 7.1. For 1 ≤ q < 2, the tails of π (i) pr (x i ) defined in [START_REF] Jun | Monte Carlo strategies in scientific computing[END_REF] are heavier than Gaussian tails, and hence Assumptions 3.1 and 4.2 may not be satisfied. Nonetheless, one can still numerically apply the proposed dimension reduction methods without having the error bounds in ( 14) and [START_REF] Dwivedi | Log-concave sampling: Metropolis-Hastings algorithms are fast[END_REF]. In this case, we set Γ to be the identity matrix in accordance with the fact that the prior components are independent and identically distributed. Remark 7.2 (Other sparsity-inducing prior). There exist other shrinkage priors similar to Besov priors, in which the random function is expressed as a weighted linear combination of basis functions and the associated random weights follow other type of heavy tail distributions. For example, the horseshoe prior and the Student's t prior. See [START_REF] Carvalho | Handling sparsity via the horseshoe[END_REF] for further discussions and references therein. The coordinate selection technique introduced here may also be applicable to those shrinkage priors. ‡ This pdf is used for demonstrating the intuition rather than a rigorous characterization, as it is defined with respect to the (non-existent) infinite-dimensional Lebesgue measure. However, the finite-dimensional discretization of the Besov measure, which is used in numerical simulations, has a pdf in this form with respect to Lebesgue measure.

Dimension reduction via prior normalization

Alternatively, we consider the case where the prior can be defined as the pushforward of the standard Gaussian measure with pdf µ(x) ∝ exp(-

1 2 x 2 2 ) under a C 1 -diffeomorphism T : R d → R d , which takes the form π pr (x) = T µ(x). (39)
In other words, π pr (x) is the pdf of the random vector X = T (Z) where Z ∼ µ(z). For the Besov-B r pp prior defined in [START_REF] Jun | Monte Carlo strategies in scientific computing[END_REF], the diffeomorphism T has a diagonal form

T (z) = (T 1 (z 1 ), . . . , T d (z d )) with T i (z i ) = Φ -1 i (Ψ(z i ))
, where Φ i (•) is the cumulative density function (cdf) of π (i) pr (x i ) defined in [START_REF] Jun | Monte Carlo strategies in scientific computing[END_REF] and Ψ(•) is the cdf of the standard Gaussian. We provide details of the cdf Φ(•) in Appendix E.

The invertibility of T allows us to reparametrize the Bayesian inverse problem in terms of the variable z = T -1 (x), which is endowed with the Gaussian prior µ. With this change of variable, the likelihood function becomes z → L y (T (z)), and thus the matrix H(y) used to reduce the dimension of z should be

H z (y) = R d ∇T (z) ∇ log L y (T (z)))(∇ log L y (T (z)) ∇T (z) µ(z)dz,
in the data-dependent case and E[H z (Y )] in the data-dependent case. For the optimal parameterreduced forward model in the Gaussian likelihood case (cf. Section 4), the forward model x → G(x) is replaced by z → G(T (z)). This way, the matrix H G should be replaced by

H G z = R d ∇T (z) ∇G(T (z)) Σ -1 obs ∇G(T (z))∇T (z) µ(z)dz.
Using either of these matrices, we obtain a projector P r to reduce the dimension in the variable z = z r + z ⊥ , where z r = P r z and z ⊥ = (I d -P r )z. In term of the original variable x, the dimension reduction method allows one to identify x r = T (P r T -1 (x)) with the observed data, while x ⊥ = T ((I d -P r )T -1 (x)) is informed by the prior only. Since x r and x ⊥ are nonlinear with respect to x, the resulting method can be interpreted as a nonlinear dimension reduction method.

Example 1: elliptic PDE

We first validate our methods using an inverse problem of identifying the coefficient of a twodimensional elliptic PDE from point observations of its solution.

Problem setup

Consider the problem domain Ω = [0, 1]×[0, 1], with boundary ∂Ω. We denote the spatial coordinate by s = (s 1 , s 2 ) ∈ Ω. We model the steady state potential solution field p(s) for a given conductivity field κ(s) and forcing function f (s) using the Poisson's equation

-∇ • (κ(s)∇p(s)) = f (s), s ∈ Ω. ( 40
) Let ∂Ω n = {s ∈ ∂Ω | s 2 = 0} ∪ {s ∈ ∂Ω | s 2 = 1}
f (s, t) = c exp - 1 2r 2 s -a 2 -exp - 1 2r 2 s -b 2 , ∀t ≥ 0,
with r = 0.05, which is the superposition of two Gaussian-shaped sink/source terms centered at a = (0.5, 0.5) and b = (2.5, 0.5), scaled by a constant c = 6 × 10 -4 . The conductivity field κ(s) is endowed with a log-normal prior. That is, letting x(s) = log κ(s), the Gaussian process prior for x(s) is defined by the stochastic PDE (see [START_REF] Lindgren | An explicit link between gaussian fields and gaussian Markov random fields: the stochastic partial differential equation approach[END_REF] and references therein):

-x(s) + γx(s) = W(s), s ∈ Ω, ( 41
)
where is the Laplace operator and W(s) is the white noise process. We impose a no-flux boundary condition on the above SPDE and set γ = 10. Equations ( 40) and ( 41) are solved using the finite element method with bilinear basis functions. A mesh with 80×80 elements is used in this example. This leads to n = 6400 dimensional discretised parameters. We generate three "true" conductivity fields from the prior distribution and use them to simulate synthetic observed data sets. The true conductivity fields and the simulated potential fields are shown in Figure 1. Observations of the potential fields are measured at the m = 36 discrete locations shown as black dots in Figure 1. We set the standard derivation of the observation noise to σ = 0.0415, which corresponds to a signal-to-noise ratio of about 20.

Figure 2: Elliptic PDE example. KL divergence of the full posterior densities from the approximate posterior densities defined by projectors with various ranks. Each column represents posteriors conditioned on a given data set. Top row: approximate posteriors are defined by the optimal parameter-reduced likelihood. Bottom: approximate posteriors are defined by the optimal parameter-reduced forward model.

Low-dimensional posterior approximations

We first compare the approximate posterior densities defined by the data-free dimension reduction with that of the data-dependent dimension reduction and that of the truncated Karhunen-Loéve expansion of the prior. We build five sets of projectors: the data-free projectors as detailed in Section 3.3, three sets of data-dependent projectors (see Section 3.2) that correspond to three synthetic data sets, and projectors defined by leading eigenvectors of the prior covariance (i.e. the truncated Karhunen-Loéve, referred to as the "prior-based projector" from hereinafter). For each of the data sets, the corresponding data-dependent projectors are constructed using the adaptive MCMC algorithm of [START_REF] Zahm | Certified dimension reduction in nonlinear Bayesian inverse problems[END_REF]. Each set consists of projectors with ranks r = 2 2 , 2 4 , ..., 2 8 . For each projector, we compute the KL divergences of the full posteriors from the approximated posterior densities defined by the optimal parameter-reduced likelihood [START_REF] Beskos | MCMC methods for diffusion bridges[END_REF]. The results are shown in the top row of Figure 2. Using the same set of projectors, we also compare the KL divergences of the full posteriors from the resulting approximated posterior densities defined by the optimal parameterreduced forward model [START_REF] Dashti | Besov priors for Bayesian inverse problems[END_REF]. The results are shown in the bottom row of Figure 2.

In these experiments, we estimate the KL divergence using Monte Carlo integration with N posterior samples, which yields

D KL (π y pos πy pos ) ≈ 1 N N i=1 log L y (x (i) )-log Ly (x (i) ) +log 1 N N i=1 exp log Ly (x (i) )-log L y (x (i) ) ,
where the second sample average accounts for the ratio between normalizing constants. For approximations that are close to the full posterior, using a reasonable number of (independent) posterior samples, e.g., 10 5 used here, make the standard deviations of the estimated KL divergence insignificant compared with the mean estimates in our numerical examples. We observe that the optimal parameter-reduced likelihood and the optimal parameter-reduced forward model result in approximate posteriors with similar accuracy. For sufficiently large ranks (r ≥ 8), the most accurate approximate posterior densities are obtained by the data-dependent projectors of the corresponding data set, followed by those obtained by the data-free projectors. We also observe that, for each data set, the data-dependent projectors constructed using other data sets result in less accurate approximations. By allowing a marginal loss of accuracy compared to the data-dependent construction, the data-free construction bypasses the computationally costly online dimension reduction process for every new data set. Compared with the prior-based dimension reduction, which is also an offline method, the data-free construction offers significantly more accurate approximations in this example.

For each of the data sets, we also compare the errors of the approximate posterior densities with the bounds defined in ( 14) and [START_REF] Dwivedi | Log-concave sampling: Metropolis-Hastings algorithms are fast[END_REF]. Note that the right hand sides of ( 14) and ( 25) are the same up to the constant κ in this example. We plot the errors and the bounds (with κ = 1 for Gaussian prior) in Figure 3, in which all approximate posterior densities are defined by the datafree projectors. In this example, we observe that the errors of the approximate posterior densities follow the same trend as their corresponding error bounds. Note that both ( 14) and [START_REF] Dwivedi | Log-concave sampling: Metropolis-Hastings algorithms are fast[END_REF] give upper bounds on the expected KL divergence, and thus they may not bound the KL divergence for a realization of the data. 14) and ( 25) (with κ = 1) are compared to KL divergence of the full posterior densities from the approximate posterior densities defined by the data-free projectors with various ranks. Left: approximate posteriors are defined by the optimal parameter-reduced likelihood. Right: approximate posteriors are defined by the optimal parameterreduced forward model.

Subspace accelerated sampling

Table 1: Acronyms of various inference algorithms used in numerical comparisons.

OL

approximate inference using Algorithm 4 and the optimal parameterreduced likelihood function in (27) PM exact inference using the pseudo-marginal method (Algorithms 5 and 6) OF approximate inference using Algorithm 4 and the optimal parameterreduced forward model in (28) DA exact inference using the delayed acceptance algorithm (Algorithm 7) with the approximation defined by the parameter-reduced forward model in (28) H-MALA exact inference using the Hessian preconditioned Langevin MCMC [START_REF] Petra | A computational framework for infinitedimensional Bayesian inverse problems, part ii: Stochastic newton MCMC with application to ice sheet flow inverse problems[END_REF] PCN exact inference using the preconditioned Crank-Nicolson MCMC [START_REF] Beskos | MCMC methods for diffusion bridges[END_REF][START_REF] Simon L Cotter | MCMC methods for functions: modifying old algorithms to make them faster[END_REF] We demonstrate the sampling performance of various approximate and exact inference algorithms introduced in Sections 5 and 6 using the posterior density conditioned on the first data set. All the methods used in the comparison and their acronyms are summarized in Table 1.

We use the Hessian-preconditioned Metropolis-Adjusted Langevin Algorithm (H-MALA) and the preconditioned Crank-Nicolson (PCN) MCMC as reference MCMC methods for sampling the full-dimensional posterior. Since H-MALA uses the low-rank decomposition of the Hessian matrix of the logarithm of the posterior density computed at the maximum a posteriori point to precondition MCMC, it can also be viewed as a data-dependent subspace-accelerated method. We refer to [START_REF] Cui | Dimension-independent likelihood-informed MCMC[END_REF][START_REF] Petra | A computational framework for infinitedimensional Bayesian inverse problems, part ii: Stochastic newton MCMC with application to ice sheet flow inverse problems[END_REF] for a detailed discussion. In order to make a fair comparison with H-MALA, the MCMC algorithm we use on our data-free informed subspace is based on a Langevin proposal preconditioned by the same Hessian matrix used by H-MALA projected onto the data-free informed subspaces.

In Figure 4, the contours of the marginal posterior densities (marginalized onto the first two data-free LIS basis vectors) produced by approximate inference methods (with r = 16 and r = 48) are compared with those produced by their exact inference modifications (with r = 48). We can observe that the results produced by approximate inference methods approach those of their modifications as the rank of informed subspace increases.

To measure the efficiency of various MCMC methods, we use the average integrated autocorrelation times (IACTs) of parameters

τ = 1 d d i=1 IACT(x i ),
where IACT(x i ) is the IACT of the i-th component of x. See [START_REF] Jun | Monte Carlo strategies in scientific computing[END_REF]Section 12.7] for the definition of IACT. The data-free projectors with different ranks r and two different sample sizes N = 2 and N = 5 are used in this experiment. Here H-MALA and PCN are used as base cases to benchmark those MCMC methods accelerated by the informed subspace. All the methods are simulated for 5 × 10 5 iterations and repeated 10 times to report the mean and the standard deviation of τ . The initial state of all the simulations are randomly selected from a pre-computed Markov chain of posterior samples to avoid burn-in. The results are reported in Table 2.

For the approximate inference methods (OL and OF), the average IACTs consistently increase with the rank of the projectors, as the sampling performance of the Langevin proposal is expected Compared to the OL method, the PM method (the exact inference counterpart for OL) has a different behavior. Here we recall that the sample-averaged parameter-reduced likelihood, L y N , in the PM method is a random estimator, whereas L y N in the OL method is deterministic because of the usage of prescribed samples. The standard deviation of the logarithm of L y N in Table 2 confirms that low-rank projectors have rather large Monte Carlo errors as the approximation accuracy is controlled by the rank truncation (cf. ( 14)). The exactness of the PM method comes at the cost of Monte Carlo error, which is controlled by the sample size N and the rank of the projector. We observe that increasing either the rank or the sample size can narrow the gap between the IACTs produced by PM and its OL counterpart. This experiment clearly suggests that PM needs to balance the sample size N and the rank of the projector to achieve the optimal performance.

Compared to the OF method, the DA method (the exact inference counterpart for OF) produces the largest IACTs among all subspace inference methods. This result is not surprising, as the second stage acceptance/rejection of DA necessarily deteriorates the statistical performance [START_REF] Andrés | Markov chain Monte Carlo using an approximation[END_REF]. In Table 2, we observe that the second stage acceptance rates, E[β], increase with more accurate approximations obtained with higher projector ranks and larger sample sizes. As the result, the gaps between the IACTs produced by OF and DA are smaller for higher projector ranks and larger sample sizes.

Overall, approximate inference methods have better statistical performance compared to other methods in this example (cf. Table 2) and can obtain reasonably accurate results as shown in Figures 2 andFigure 4. With the additional cost that comes in the form of either Monte Carlo error (PM) or the second stage acceptance/rejection (DA), the exact inference modifications produce Markov chains with larger IACTs. Among all the exact inference methods, PM produces smaller IACTs compared with the full-dimensional H-MALA, PCN, and DA.

It is worth to mention that each iteration of the subspace MCMC method needs N number of forward model simulations, whereas H-MALA requires only one forward model simulation per iteration. In this example, approximate inference methods (OL and OF) with N = 2 still outperforms H-MALA in terms of IACTs per model evaluation. Exact inference methods, however, need more model evaluations than H-MALA to obtain the same number of effective samples (we will show in the subsection another example where H-MALA is outperformed by PM and DA). Notice that the forward model evaluations in each iteration can be embarrassingly parallelized: with parallel computing resources available, the subspace MCMC methods can still be more efficient than H-MALA in terms of the effective sample size per wall-clock time.

Example 2: PET with Poisson data

The second example is a two dimensional PET imaging problem, where we aim to reconstruct the density of the object from integer-valued Poisson observed data. We use here a Besov prior for which we access the coordinate selection technique and the prior normalization method presented in Sections 7.2 and 7.3.

Problem setup

In PET imaging, the goal is to identify an object of interest located inside a domain Ω subjected to gamma rays. The rays travel through Ω from multiple sources and the detectors count the number of incident photons (thus the data are integer-valued), see Figure 5a. The object of interest is described by its density of mass which is represented by s → exp(f (s)), where f : Ω → R follows a Besov-B 2 11 prior with the Haar wavelet, see Section 7.1. The change of intensity of a gamma ray along the path, i (s), s ∈ Ω, can be modelled using Beer's law:

I d,i = I s,i exp - i (s) exp f (s) ds , (42)
where I d,i ∈ R ≥0 and I s,i ∈ R ≥0 are the intensities at the detector and at the source, respectively. We assume that all the gamma ray sources have the same intensity, I s,i = I s for i = 1, ..., m.

In this example, the domain Ω is discretized into a regular grid with d cells and the logarithm of the density is assumed to be piecewise constant. This yields the discretized parameter x ∈ R d . The line integrals in [START_REF] Nouy | Low-Rank Tensor Methods for Model Order Reduction[END_REF] are approximated by

i (s) exp(f (s))ds ≈ n j=1 A ij exp(f j ),
where A ij ∈ R ≥0 is the length of the intersection between line i and cell j, and exp(f j ) is the discretized density in cell j. By discretizing the wavelet basis on the same grid and following the parametrization discussed in Section 7, we can write

f = Bx,
where B ∈ R d×d consists of discretized basis functions and x consists of associated coefficients. In this setting, x follows a product-form Laplace distribution given by [START_REF] Jun | Monte Carlo strategies in scientific computing[END_REF] with p = 1 and with the scale parameter arbitrarily set to γ = 1. Suppose we have a total of m number of gamma ray paths and the corresponding matrix A ∈ R m×d , the forward model G : R d → R m is defined as

G(x) = I s exp(-A exp(Bx)). ( 43
)
We consider a PET setup shown in Figure 5a We denote the observed data by y ∈ N m where each element y i is associated with the ith gamma ray in the model. For the i-th gamma ray, recall that the intensity at the detector is computed by G i (x) for some input parameter x, and then the probability mass function of observing the counting data Y i = y i is given by

P(Y i = y i |x) = G i (x) y i exp(-G i (x)) y i ! .
Suppose we can observe the counting data at all the detectors and assume the measurement processes are independent, we can write the likelihood function with the complete data as

L y (x) = m i=1 P(Y i = y i |x) = m i=1 G i (x) y i exp(-G i (x)) y i ! . ( 44
)
As shown in Appendix F, the Fisher information matrix of the above likelihood function takes the form

I(x) = ∇G(x) M (x)∇G(x), (45)
where M (x) is a diagonal matrix with M ii (x) = G i (x) -1 along its diagonal. We generate three "true" density functions from the prior distribution and use them to simulate synthetic data sets. The true density functions and the simulated data are shown in Figure 5b.

Numerical results using coordinate selection

We first present the results obtained by applying the coordinate selection method (cf. Section 7.2). Similar to the first example, here we will first compare the accuracy of approximate posterior densities defined by various approaches and projectors, and then benchmark the performance of MCMC methods. We adopt the same setup and acronyms as in Example 1 and Table 1.

For the approximate posterior densities, five sets of projectors built from selected coordinates, including the data-free projectors, three sets of data-dependent projectors (corresponding to three data sets), and the prior-based truncated wavelet basis, are considered. Each set consists of projectors with ranks r = 2 3 , 2 4 , ..., 2 7 . The KL divergences of the full posteriors from the approximated posterior densities defined by the optimal parameter-reduced likelihood [START_REF] Beskos | MCMC methods for diffusion bridges[END_REF] are shown in the top row of Figure 6, while those of the optimal parameter-reduced forward model [START_REF] Dashti | Besov priors for Bayesian inverse problems[END_REF] are shown in the bottom row of Figure 6.

In this example, we observe similar results as the results of the elliptic PDE example. The optimal parameter-reduced likelihood and the optimal parameter-reduced forward model result in approximate posteriors with similar accuracy. The most accurate approximate posterior densities are obtained by the data-dependent projectors of the corresponding data set, followed by those obtained by the data-free projectors. For each data set, the data-dependent projectors constructed using other data sets result in less accurate approximations in general. However, the accuracy gaps between the data-free projectors and the data-dependent projectors (using other data sets) are not as significant as the elliptic PDE example. This can be caused by either the coordinate selection method or the rather large data size in this example. Compared with the prior-based dimension reduction, which is also an offline method, the data-free construction offers significantly more accurate approximations in this example. Overall, the data-free dimension reduction provides reasonably accurate posterior approximations for the Poisson observation process considered here. Although it remains an open question if the bounds in ( 14) and ( 25) can be applied for Besov priors, we provide a comparison of the errors of the approximate posterior densities (defined by the data-free projectors) with the bounds. The results are shown in Figure 7 with κ being replaced by 1. Interesting, we still observe that the errors of the approximate posterior densities follow the same trend as their corresponding bounds.

We then compare the performance of various subspace driven inference methods. In Figure 8, the contours of the the marginal posterior densities produced by approximate inference methods (with r = 16 and r = 48) are compared with those produced by their exact inference modifications (with r = 48). In this example, we observe that the contours produced by approximate inference methods are visually similar to those of exact inference methods. In addition, with increasing ranks, the contours produced by approximate inference methods approach those of the exact inference methods.

Figure 8: Same as Figure 4, but for PET and data set 1. Here r = 16, 48 are used in OL and OF, and r = 48 is used in PM and DA. Here (x 1 , x 2 ) represent the first two coordinates selected by the data-free method. For the approximate inference methods (OL and OF), the average IACTs consistently increase with the rank of the projectors, as the sampling performance of the Langevin proposal is expected to decay with underlying the parameter dimension. Both OL and OF produce significantly smaller IACTs compared with the full-dimensional PCN and H-MALA method. Compared to the OL method, the PM method, has a slightly higher IACTs in this example. This rather mild loss of performance (compared with the elliptic PDE example) is justified by the rather small values of L y N (with N = 2, 5) in Table 3. Compared to the OF method, the DA method, again produces the largest IACTs among all subspace inference methods. However, the loss of performance here is not as severe as the elliptic PDE example, this is also justified by the improved second stage acceptance rates, E[β].

Overall, approximate inference methods have better statistical performance compared to other methods in this example and can obtain reasonably accurate results as shown in Figures 6 and Figure 8. With improved approximation errors, the exact inference methods also produces Markov chains with better mixing. Among all the exact inference methods, PM produces significantly smaller IACTs compared with other methods.

Numerical results using prior normalization

Then, we present the results obtained by applying the prior normalization method (cf. Section 7.3). The KL divergences of the full posteriors from the approximated posterior densities are shown in Figure 9. Here the result of prior-based dimension reduction is not presented, as the prior in the transformed space an identity covariance matrix. We observe similar results as those obtained by the coordinate selection. We notice that the accuracy gaps between the data-free projectors and the data-dependent projectors (built using other data sets) are more significant compared with those obtained by the coordinate selection. The comparison of the errors of the approximate posterior densities (defined by the data-free projectors) with the bounds in ( 14) and ( 25) are provided in Figure 10. Here we have κ = 1 because the transformed coordinate is endowed with a Gaussian prior. We observe that the errors of the approximate posterior densities follow the same trend as their corresponding bounds. The IACTs of various MCMC methods are reported in Table 4. Again, the efficiency of subspace MCMC methods defined by the prior normalization is very close to that defined by the coordinate selection. Overall, both the coordinate selection and the prior normalization can be applied in this example to obtain accurate reduced-dimensional posterior approximations and derive efficient subspace MCMC methods. approach identifies the computationally costly subspace construction in an offline phase. Our data-free dimension reduction is certified in the sense that its development is directly guided by factorizable posterior approximations and associated error bounds. The factorizable posterior approximations naturally offer dimension robust sampling methods for exploring the approximate posterior densities. More interestingly, by adding minor modifications to those approximate inference algorithms, we further develop exact inference methods using the pseudo-marginal approach and the delayed acceptance approach. The resulting exact inference methods also scale well with parameter dimensionality, as the backbone of those methods is based on the dimension robust approximate inference methods. We also demonstrate the efficiency of our datafree dimensional reduction and various inference methods on two inverse problems involving a two-dimensional elliptic PDE with a Gaussian process prior and a PET problem with Poisson data and a Besov-B 2 11 prior. For the last equality we used the fact that y → L y (x) is a pdf so that Consider a Metropolis-Hastings algorithm which targets the pdf π y,N tar defined by (32) using the following proposal density

q x † r , {x †(i) ⊥ } N i=1 x r , {x (i)
⊥ } N i=1 = q(x † r |x r ) N i=1 π pr (x †(i) ⊥ |x † r ), (B.1)
where q(x † r |x r ) is the same proposal density as the one used at step 1 of Algorithm [START_REF] Beskos | Geometric MCMC for infinite-dimensional inverse problems[END_REF]. The acceptance probability of this Metropolis-Hastings algorithm is given by

α x † r , {x †(i) ⊥ } N i=1 x r , {x (i)
⊥ } N i=1 = min    1, π y,N tar (x † r , {x †(i) ⊥ } N i=1 ) q x r , {x (i)
⊥ } N i=1 x † r , {x †(i) ⊥ } N i=1 π y,N tar (x r , {x (i)
⊥ } N i=1 ) q x † r , {x †(i) ⊥ } N i=1 x r , {x (i)
⊥ } N i=1    = min    1, π pr (x † r ) N i=1 L y (x † r + x †(i) ⊥ ) q(x r |x † r ) π pr (x r ) N i=1 L y (x r + x (i) ⊥ ) q(x † r |x r )    ,
which is precisely α N (x † r |x r ) defined in [START_REF] Kaipio | Statistical and computational inverse problems[END_REF]. Note that the first two steps of Algorithm 5 consists of drawing a sample (x † r , {x †(i) ⊥ } N i=1 ) from the proposal (B.1). This way, Algorithm 5 can be interpreted as a MCMC algorithm which targets π y,N tar . It remains to show that the marginal distribution

  denote the top and bottom boundaries, and ∂Ω d = {s ∈ ∂Ω | s 1 = 0} ∪ {s ∈ ∂Ω | s 1 = 1} denote the left and right boundaries. We impose the mixed boundary condition: p(s) = 0, ∀s ∈ ∂Ω d , and (κ(s)∇p(s)) • n(s) = 0, ∀x ∈ ∂Ω n , and let the forcing function take the form

Figure 1 :

 1 Figure 1: Setup of three test cases for the elliptic PDE example. The observation locations are shown as dots. Each column represent a test case, in which the top row shows the true conductivity fields and the bottom row shows the corresponding potential field.

Figure 3 :

 3 Figure 3: Elliptic PDE example. The bounds in (14) and (25) (with κ = 1) are compared to KL divergence of the full posterior densities from the approximate posterior densities defined by the data-free projectors with various ranks. Left: approximate posteriors are defined by the optimal parameter-reduced likelihood. Right: approximate posteriors are defined by the optimal parameterreduced forward model.

Figure 4 :r 5 r

 45 Figure 4: Elliptic PDE example. Contours of the marginal posterior densities computed by various inference algorithms using the data-free projector. (a): OL (with r = 16 and r = 48) and PM (with r = 48). (b): OF (r = 16 and r = 48) and DA (with r = 48). Here (x 1 , x 2 ) represent the directions spanned by the first two data-free LIS basis vectors.

  : the problem domain Ω = [-10, 10] 2 is discretised into a d = 64 × 64 regular grid, five radiation sources with intensity I s = 10 are positioned with equal spaces on one side of a circle, spanning a 90 • angle, and each source sends a fan of 30 gamma rays that are measured by detectors. This leads to m = 150 observations. The model setup is based on the code of [30].

( a )

 a Discretised domain of interest Ω = [-10, 10] 2 (mesh), radiation sources (red dots), and detectors (blue dots). (b) Top row: three density functions generated from the prior. Bottom row: corresponding intensity function (solids lines) and measured counting data sets (dots).

Figure 5 :

 5 Figure 5: The PET setup and three test cases.

Figure 6 :

 6 Figure 6: Same as Figure 2, but for the PET test case and where we used coordinate selection.

Figure 7 :

 7 Figure 7: Same as Figure 3, but for the PET test case and where we used coordinate selection.

Figure 9 :

 9 Figure 9: Same as Figure 2, but for PET. The prior normalization is used in dimension reduction.

Figure 10 :

 10 Figure 10: Same as Figure 3, but for PET. The prior normalization is used in dimension reduction.

D

  KL (π y pos || π y pos ) = x) Σ -1 obs (G(x) -y) π y pos (x)dx (A.1) where e(x) = G(x) -G * (x) is independent on y. Next we replace y by Y ∼ π data and we take the expectation over Y . The first term in the above expression becomes E log π data (Y ) π data (Y ) = R m log π data (y) π data (y) π data (y)dy = -D KL (π data || π data ). (A.2)Next, by definition of π y pos (x), we have E

Rm 4 )

 4 L y (x)dy = 1. Using the same arguments, we haveE R d e(x) Σ -1 obs G(x) -Y π Y pos (x)dx = R d ×R m e(x) Σ -1 obs G(x) -y L y (x)π pr (x)dxdy = R d e(x) Σ -1 obs G(x)π pr (x)dx -The last equality is obtained by noting that the expectation of the data knowing the parameter x is R m yL y (x)dy = G(x). Combining (A.2) (A.3) and (A.4), we obtainE D KL (π Y pos || π Y pos ) (A.1) = -D KL (π data || π data ) + 1 2 R d G(x) -G(x) 2 Σ -1 obs π pr (x)dx,which concludes the proof.Appendix B. Proof of Proposition 6.1

Table 3 :

 3 Same as Table2, but for PET and data set 1. The coordinate selection is used in dimension reduction. Ranks of the subspaces are chosen to be 16, 32, and 48.

		OL	IACT PM	√	var[log L y N ] OF	IACT DA	E[β]	IACT HMALA PCN
	N = 2 N = 5	r = 16 33.2±1.7 85.1±2.7 1.54±0.02 33.9±1.1 214±44 0.18±0.06 r = 32 40.0±1.8 54.1±3.1 0.61±.007 41.0±2.2 87.8±6.5 0.55±0.01 r = 48 45.3±1.2 49.4±2.6 0.45±.002 46.0±2.2 73.5±5.8 0.66±0.01 r = 16 31.4±1.9 60.0±6.2 0.93±.006 31.8±1.4 220±65 0.22±0.04 r = 32 40.8±2.5 47.6±2.5 0.39±.004 42.8±2.4 88.0±6.8 = 48 46.5±1.4 0.29±.001 69.5±4.0 0.67±0.01	95.9±3.3 387±79
		We use the average IACTs of the density function, τ = 1 d	d i=1 IACT(f i ), to measure
	efficiency of various MCMC methods. The results are reported

Table 3 .

 3 Here both PCN H-MALA are implemented to sample the posterior the transformed coordinate equipped with a Gaussian prior Section 7.3).

Table 4 :

 4 Same as Table 2, but for PET and data set 1. The prior normalization is used in dimension reduction. Ranks of the subspaces are chosen to be 16, 32, and 48. = 32 42.8±2.0 55.1±2.9 0.64±.006 41.2±1.6 86.3±6.2 0.55±0.01 r = 48 45.0±2.4 51.8±2.0 0.46±.005 44.3±2.2 74.6±7.4 0.65±0.01 N = 5 r = 16 35.1±1.9 54.1±4.1 0.88±0.01 32.8±2.8 151±21 0.26±0.02 r = 32 45.0±1.7 49.0±1.9 0.41±.003 42.0±2.6 83.1±5.1 0.55±0.01 r = 48 45.9±2.9 46.3±2.2 0.29±.003 44.4±0.8 70.6±3.7 0.66±0.01 10. Conclusion We present a new data-free strategy for reducing the dimensionality of large-scale statistical inverse problems. Compared to existing gradient-based dimension reduction technique, this new

		OL	IACT PM	√	var[log L y N ] OF	IACT DA	E[β]	IACT HMALA PCN
	N = 2	r = 16 35.8±1.7 81.4±6.0 1.48±0.02 33.5±1.8 168±23 0.25±0.02 r
								95.9±3.3 387±79

Acknowledgments

T. Cui and O. Zahm would like to acknowledge support from the INRIA associate team UNQUESTIONABLE. T. Cui also acknowledges support from the Australian Research Council.

π data (y) . By definition of L y (x) and L y (x) we have π y,N tar (x r ) is the marginal posterior π y pos (x r ). We can write π y,N tar (x r ) =

⊥ , ..., x

⊥ ...dx

which concludes the proof.

Appendix C. Proof of Proposition 6.3

Recall that {(X (j) r , {X

)} j≥1 admits π y,N tar (32) as the invariant density, see Proposition 6.1. It remains to prove that {X (j) } j≥1 admits π y pos (x) as the invariant density. For a given state (X (j) r , {X (j,i)

⊥ } N i=1 is selected with respect to the probability

Thus, we need to prove that the pdf π(x) where x = x r +x ⊥ is the posterior density π(x) = π y pos (x). We can write

⊥ ...dx

⊥ ...dx

where π(x ⊥ |x r , {x

⊥ } N i=1 ). By construction we have

where δ x (k) ⊥ denotes the Dirac mass function at point x (k)

⊥ . We can write

⊥ ...dx

which concludes the proof.

Appendix D. Proof of Proposition 6.4

To show the result of Proposition 6.4, we first interpret the first stage acceptance/rejection and the conditional prior sampling π pr (x † ⊥ |x † r ) as a joint proposal acting in the full parameter space Im(P r )⊕Ker(P r ). The proposal q(•, |x r ) and the acceptance probability α(x † r |x r ) defines an effective proposal distribution

where δ x r (•) denotes the Dirac delta and the term in the bracket represents the probability of a proposal candidate being rejected. Then, we can define a joint proposal distribution

for the MH to sample the full posterior. Following the exactly same derivation in [START_REF] Andrés | Markov chain Monte Carlo using an approximation[END_REF], one can show that accepting

defines a Markov transition kernel with the full posterior π y pos (x r + x ⊥ ) as its invariant distribution. Since the above acceptance probability is only used in the case where the first stage proposal candidate x † r is accepted, i.e., x † r = x r , we do not need to consider the Dirac delta term. This way, the above acceptance probability can be simplified as

.

,

into the above equation, we obtain

, which is identical to the second stage acceptance probability in [START_REF] Le | Spectral methods for uncertainty quantification: with applications to computational fluid dynamics[END_REF]. Thus, the result follows.

Appendix E. Cumulative density function of p(x) ∝ exp(-γ|x| p )

Given the pdf p(x) = 1 c γ,p exp(-γ|x| p ), x ∈ R, we want to find its normalizing constant c γ,p and cdf. Using symmetry, the normalizing constant takes the form

and the cdf can be expressed as

We first introduce the change-of-variable z = γt p so that

This yields

where Γ lower (•, •) is the lower incomplete gamma function. Following a similar derivation, we obtain c γ,p = 2 p -1 γ -1 p Γ(p -1 ), where Γ(•) is the Gamma function. This way, we have the cdf

There are two notable special cases. The Gaussian distribution N (0, σ 2 ) can be specified using γ = (2σ 2 ) -1 and p = 2, in which the cdf can be equivalently expressed using the error function. The Laplace distribution can be specified using p = 1, so that the cdf yields a simpler (but equivalent) expression in the form of

Appendix F. Derivation of Fisher information matrices

Here we derive the Fisher information matrix for the Poisson likelihood case. Recall the Fisher information matrix

Defining the predata η = G(x), we can express the gradient of the likelihood function as

where

This way, the Fisher information matrix can be rewritten as where M (x) is a diagonal matrix with M ii (x) = G i (x) -1 along its diagonal.