{'fraction_non_alphanumeric': 0.050336774142122676, 'fraction_numerical': 0.012727007027752315, 'mean_word_length': 3.9996014877789587, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 5, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 32, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Certified dimension reduction in nonlinear Bayesian inverse problems

We propose a dimension reduction technique for Bayesian inverse problems with nonlinear forward operators, non-Gaussian priors, and non-Gaussian observation noise. The likelihood function is approximated by a ridge function, i.e., a map which depends non-trivially only on a few linear combinations of the parameters. We build this ridge approximation by minimizing an upper bound on the Kullback--Leibler divergence between the posterior distribution and its approximation. This bound, obtained via logarithmic Sobolev inequalities, allows one to certify the error of the posterior approximation. Computing the bound requires computing the second moment matrix of the gradient of the log-likelihood function. In practice, a sample-based approximation of the upper bound is then required. We provide an analysis that enables control of the posterior approximation error due to this sampling. Numerical and theoretical comparisons with existing methods illustrate the benefits of the proposed methodology.


Introduction
Solving Bayesian inverse problems [38,61] is a challenging task in many domains of application, due to the complexity of the posterior distribution. One of the primary sources of complexity is the dimension of the parameters to be inferred, which is often high or in principle infinitee.g., when the posterior is a distribution over functions or their finite-dimensional discretization. High dimensionality presents difficulties for posterior sampling: care is required to design sampling algorithms that mix effectively while remaining robust under refinement of the discretization. High dimensionality also raises significant hurdles to the use of model reduction or approximation schemes (e.g., [10,14,19,23,48,57]) that attempt to reduce the cost of likelihood or forward model evaluations.
Successful strategies for high-dimensional Bayesian inverse problems typically exploit the presence of some low-dimensional structure. A common structure in inverse problems is that the posterior is a low-dimensional update of the prior, in the sense that change from prior to posterior is most prominent on a low-dimensional subspace of the parameter space. Put another way, the likelihood is influential, relative to the prior, only on a low-dimensional subspace (we will make these intuitive notions more precise later). Sources of such structure include the smoothing properties of the forward operator and limitations on the number or accuracy of the observations. In the linear-Gaussian case, this structure is already well understood. For instance, low-rank approximations of the (prior-preconditioned) Hessian of the log-likelihood have been used in [27] to approximate the posterior covariance. In [60], this approach is shown to yield optimal approximations of the posterior covariance and of the posterior mean. A heuristic extension of this approach to nonlinear forward models, known as the likelihood-informed subspace (LIS) method, is proposed in [21] and shown to perform well in many applications [6,24,42]. A similar idea underlies the active subspace method applied to Bayesian problems [16]. Identifying the subspace on which changes from the prior to posterior are most prominent also provides a foundation for infinite-dimensional MCMC algorithms that "split" the parameter space, such as the dimension-independent likelihoodinformed (DILI) MCMC samplers of [20]. In [22], LIS-based dimension reduction is shown to be an effective prelude to model order reduction, where the forward model is replaced with a computationally inexpensive approximation. Yet none of these dimension reduction approaches, outside of the linear-Gaussian setting of [60], comes with a rigorous error analysis. The subject of this paper is a new dimension reduction method for the nonlinear/non-Gaussian setting that is, in contrast, certified -in the sense that the error induced by the approximation is endowed with a computable upper bound.
Other forms of dimension reduction for Bayesian inverse problems have also been proposed, besides the "update" form of dimension reduction described above. For instance, in [45,49], a truncated Karhunen-Loève decomposition of the prior distribution is used to reduce the parameter dimension of the entire inverse problem. This approach exploits only the low-dimensional structure of the prior distribution, however, and does not take advantage of structure in the likelihood function or forward model. In [46], a greedy algorithm is used to identify a parameter subspace capable of reproducing the forward model. This approach, in contrast, does not take advantage of the prior correlation structure, and moreover can remove directions that are uninformed by the likelihood but that still retain large variation under the true posterior. More recent results on the intrinsic dimension of linear Bayesian inverse problems [1] reinforce the idea that the update from prior to posterior should instead be a central object, e.g., when characterizing the performance of importance sampling schemes.

Contribution
In this paper we propose a methodology to detect and exploit the low-dimensional structure of the update from prior to posterior. Our approach addresses Bayesian inverse problems with nonlinear forward operators, non-Gaussian priors, and non-Gaussian observation noise: the posterior ν is a measure on R d given by dν ∝ f dµ where µ is the prior measure and f is any unnormalized likelihood function defined on supp(µ). Our approach consists in approximating x → f (x) by a ridge function, i.e., a function of the form x → h(Ax) which depends non-trivially only on r d linear combinations of the parameters Ax, with A ∈ R r×d . More precisely, we seek a controlled approximation such that the Kullback-Leibler (KL) divergence D KL (ν||ν r ) from the resulting posterior approximation dν r (x) ∝ h(Ax)dµ(x) to the posterior ν is below some user-defined threshold. As our main contribution, we derive in Theorem 2.10 and Corollary 2.11 an upper bound on the KL divergence. Minimizing this bound permits us to construct the ridge approximation so that the bound falls below the threshold.
This bound is obtained via logarithmic Sobolev inequalities, an important class of inequalities in measure theory [31,32,51] with many implications on concentration of measure phenomena [9,43]. Using logarithmic Sobolev inequalities requires some assumptions on the prior distribution and on the regularity of the likelihood function. In particular, we need the gradient of the log-likelihood to be square-integrable over the posterior distribution. A similar methodology has been proposed in [65] to reduce the input dimension of multivariate functions. In that paper, a bound on the function approximation error in L 2 norm is obtained via Poincaré inequalities, another class of Sobolev inequality.
In the proposed dimension reduction method, the informed subspace is constructed as the dominant eigenspace of the matrix Once the informed subspace is identified, the posterior approximation is constructed by computing a conditional expectation of the likelihood function. We prove that the resulting divergence D KL (ν||ν r ) is bounded by the sum of the d − r smallest eigenvalues of H: a quickly decaying spectrum in H thus reveals the low-dimensional structure in the posterior.
The proposed method requires computing (i) the matrix H and (ii) a conditional expectation of the likelihood function. In practice, we approximate both quantities with Monte Carlo estimates, requiring samples from the posterior for (i) and samples from the prior for (ii). Provided the sample sizes are sufficiently large, we prove that the resulting random posterior approximation is a quasioptimal approximation compared to the theoretical posterior approximation obtained by exact integration. This quasi-optimality result is given in expectation for (ii) and in high probability for (i). In particular, we show that the number of posterior samples to approximate (i) should scale in proportion to the rank of the matrix to be estimated, which can be much smaller than the ambient dimension. Finally, even if the method only requires a limited number of posterior samples, it can be difficult to obtain the posterior samples in (i). We thus propose several alternatives for computing (i) that do not require sampling from the exact posterior. These alternatives include an iterative procedure that builds a sequence of low-dimensional posterior approximations in order to obtain an accurate final estimate.
The outline of the paper is as follows. In Section 2 we introduce the theoretical tools to derive and to minimize the error bound. We demonstrate the benefits of this method on an analytical example. In Section 3 we propose algorithms for the numerical construction of the lowdimensional posterior approximation. We give a theoretical analysis of the convergence of these sample-based estimators. In order to provide some context for our developments, we show in Section 4 how the proposed methodology compares with existing dimension reduction methods, including a comparison with a Karhunen-Loève decomposition. Finally, Section 5 illustrates the benefits of our approach on two numerical examples.

Dimension reduction for the approximation of highdimensional distributions
Let µ be a probability distribution defined over the Borel sets B(R d ) of R d . Given a measurable function f : supp(µ) → R + such that f dµ < ∞, let ν be the probability distribution such that dν dµ ∝ f.
In the context of Bayesian inverse problems, one can view µ as the prior distribution and ν as the posterior distribution, while f represents, up to a multiplicative constant, the likelihood function.
We consider the problem of approximating the posterior ν by a probability distribution ν r such that where P r : R d → R d is a linear projector with rank r and g : Im(P r ) → R + is a Borel function called the profile function. Throughout this paper we identify the projector with its matrix representation P r ∈ R d×d . Notice that P r is not restricted to be orthogonal: it can be any matrix P r ∈ R d×d which satisfies P 2 r = P r and rank(P r ) = r, but not necessarily P T r = P r . Any vector x ∈ R d can be uniquely decomposed as where I d ∈ R d×d denotes the identity matrix. If r d, the approximation of ν by ν r consists essentially in replacing the high-dimensional likelihood f by a function of fewer variables. Indeed, x → g(P r x) = g(x r ) depends only on the variable x r ∈ Im(P r ) ∼ = R r and is constant along Ker(P r ) ∼ = R d−r . Yet ν r cannot itself be considered a low-dimensional distribution, since its support could be the same as that of ν.
We use the Kullback-Leibler divergence D KL ( · || · ) to measure the dissimilarity between probability distributions. It is defined as for any probability distributions ν 1 and ν 2 such that ν 1 is absolutely continuous with respect to ν 2 , and D KL (ν 1 ||ν 2 ) = ∞ otherwise. Given a prescribed tolerance ε ≥ 0, our goal is to build an approximation ν r of ν of the form (1) such that Of course if P r is the identity matrix and g = f , the distribution ν r is exactly the posterior ν, so that (3) is trivially satisfied. In that case, the rank of P r = I d is d and there is no dimension reduction. The goal of this section is to give sufficient conditions under which we can build an approximation ν r such that (3) holds with a rank r = r(ε) that is significantly smaller than d.
Remark 2.1. Instead of approximating ν by ν r as given by (1), we could have considered an approximation ν r of the form with h : R r → R + a Borel function and A r ∈ R r×d a matrix of full row rank. Functions of the form x → h(A r x) are a particular kind of ridge function [53]. They are used in various domains to approximate multivariate functions; see [13,28,63], for example. However, any function of the form h • A r can be written as g • P r and vice-versa, where g : Im(P r ) → R + is some Borel function and P r is some rank-r projector. More precisely, for any r ≤ d we have This means that approximating ν by ν r or by ν r are essentially the same problem: we do not gain approximation power in using ν r instead of ν r .

Optimal profile function
To begin our analysis, we assume that the projector P r is given (fixed); we will address the problem of constructing P r later in Section 2.2. In this subsection, we characterize the optimal profile function g * for a given P r , defined as a minimizer of g → D KL (ν||ν r ) over the set of positive and measurable functions. This function g * yields an optimal probability measure ν * r defined as dν * r dµ ∝ g * • P r . We will also show that, for any probability measure ν r with dνr dµ ∝ g • P r , the Kullback-Leibler divergence can then be decomposed as D KL (ν||ν r ) = D KL (ν||ν * r ) + D KL (ν * r ||ν r ). Let us denote by σ(P r ) the σ-algebra generated by P r . It is defined by σ( The following lemma corresponds to the Doob-Dynkin lemma; see, for example, Lemma 1.13 in [39]. It states that, for any projector P r , the set of all functions of the form g • P r for some measurable function g is exactly the set of all σ(P r )-measurable functions.
Lemma 2.2. Let P r ∈ R d×d be a projector. Given a measurable function g defined on R d , the function h = g • P r is σ(P r )-measurable. Conversely, given a σ(P r )-measurable function h, there exists a measurable function g defined on R d such that h = g • P r .
We denote by E µ (f |σ(P r )) : R d → R the conditional expectation of f given σ(P r ) under the distribution µ. It is the unique σ(P r )-measurable function that satisfies for any σ(P r )-measurable function h : R d → R. We consider the probability distribution ν * r defined by dν * This distribution is well defined because E µ (f |σ(P r )) is a positive function (the conditional expectation preserves the sign) and because E µ (f |σ(P r )) dµ = f dµ < ∞ (by letting h = 1 in (4)).
By the definition of the Kullback-Leibler divergence (2), we can write where Z 1 and Z 2 are two normalizing constants defined by Z 1 = f dµ = E µ (f |σ(P r ))dµ and To go from (5) to (6), we used relation (4) with h = log , which is a σ(P r )-measurable function. By (7) we have for any measurable function g. From this inequality we deduce that any function g * which satisfies g * • P r = E µ (f |σ(P r )) is a minimizer of g → D KL (ν||ν r ). By Lemma 2.2 such a function g * exists because E µ (f |σ(P r )) is a σ(P r )-measurable function.
Remark 2.3. The conditional expectation E µ (f |σ(P r )) is known to be the best approximation of a function f with respect to the L 2 µ -norm, meaning that it minimizes h → (f − h) 2 dµ among all σ(P r )-measurable functions h. Here we showed that it is also optimal with respect to the Kullback-Leibler divergence. As pointed out in [5], E µ (f |σ(P r )) is also an optimal approximation of f with respect to the class of expected Bregman divergences, which includes the Kullback-Leibler divergence and the L 2 µ -norm distance. The following proposition recalls the well-known analytical expression for the conditional expectation. For the sake of completeness, the proof is given in Section 7.1.
Proposition 2.4 (Explicit expression for the conditional expectation). Let µ be a probability distribution on R d which admits a Lebesgue density ρ. Given a rank-r projector P r , we denote by U ⊥ ∈ R d×(d−r) a matrix whose columns form a basis of Ker(P r ). Let p ⊥ ( · |P r x) be the conditional probability density on R d−r defined by for all ξ ⊥ ∈ R d−r and any x ∈ R d , with the convention p ⊥ (ξ ⊥ |P r x) = 0 whenever the denominator of (8) is zero. Then, for any measurable function f , the conditional expectation E µ (f |σ(P r )) can be written as We conclude this section with a remarkable property of the conditional expectation E µ (f |σ(P r )). Let Q r be any projector whose kernel is the same as that of P r . Then the relations P r = P r •Q r and Q r = Q r • P r hold; see, for instance, the proof of Proposition 2.2 in [65]. Then Lemma 2.2 ensures that any σ(P r )-measurable function is σ(Q r )-measurable and vice-versa. In other words, being σ(P r )-measurable or σ(Q r )-measurable are equivalent such that, by definition of the conditional expectation (4), we have In other words, the conditional expectation E µ (f |σ(P r )) is invariant with respect to the image of P r -so that the error D KL (ν||ν * r ), seen as a function of the projector P r , is actually only a function of Ker(P r ). In the context of Bayesian inference, this property means that when the optimal profile function g * is used in (1), the important feature of P r to be discovered is its kernel, and not its image. Thus reducing the dimension of a Bayesian inverse problem consists in identifying the subspace Ker(P r ) on which the data are not informative, meaning the subspace where the data do not modify the prior knowledge one has about the parameter.

A controlled approximation
In this section we show how to build a projector P r with a sufficiently large rank so that ν * r , defined by dν * r dµ ∝ E µ (f |σ(P r )), satisfies D KL (ν||ν * r ) ≤ ε for some prescribed tolerance ε. We argue that under some conditions on the likelihood f , the rank of P r can be significantly smaller than d.
We make the following assumption on the prior distribution µ. Here, the notation A B means that the matrix A − B is positive semidefinite.
Assumption 2.5. The distribution µ is supported on a convex set supp(µ) ⊆ R d and admits a Lebesgue density ρ such that ρ ∝ exp(−V − Ψ), where V and Ψ are two functions defined on supp(µ) which satisfy the following two properties.
• V is twice continuously differentiable and there exists a symmetric positive definite matrix Γ ∈ R d×d such that for all x ∈ supp(µ) we have • Ψ is a bounded function over supp(µ) such that for some κ ≥ 1.
Let us make some comments on this assumption. First consider the case κ = 1, which means that Ψ is a constant function so that ρ ∝ exp(−V −Ψ) ∝ exp(−V ). Assumption (10) corresponds to a strong convexity property of the function V . Intuitively, it means that V is "at least quadratically convex." Now, when κ > 1, Assumption (11) means that c exp(−V (x)) ≤ ρ(x) ≤ C exp(−V (x)) holds for all x ∈ supp(µ), where c > 0 and C < ∞ are two constants that satisfy C/c ≤ κ < ∞. We note that if Assumption 2.5 is satisfied, then κ and Γ are not unique. Indeed, if κ and Γ satisfy (10) and (11), then any κ ≥ κ and symmetric positive definite matrix Γ Γ also satisfy (10) and (11). In the results that follow, tighter bounds will be obtained by having the smallest κ and the largest Γ, but this optimality can be difficult to verify in general.
Assumption 2.5 provides sufficient conditions for the prior distribution µ to satisfy the following logarithmic Sobolev inequality: for any continuously differentiable function h : supp(µ) → R, where · Γ −1 denotes the norm on R d such that Logarithmic Sobolev inequalities are a fundamental tool in the analysis of Markov semigroups [4] and in concentration inequalities; see [9]. The inequality (12) results from the Bakry-Émery theorem [3,8,51], which uses Assumption (10), and from the Holley-Stroock perturbation lemma [36], which uses Assumption (11). More precisely, Proposition 3.1 in [8] states that in the case κ = 1 (i.e., Ψ is constant), Assumption (10) is sufficient to have (12). Then, following the lines of the original proof in [36] (see also the proof of Theorem 1.9 in [30], for instance), we have that (12) still holds when Ψ is such that κ > 1.
The logarithmic Sobolev inequality (12) permits control of the KL divergence from prior to posterior, D KL (ν||µ), provided that f is sufficiently regular. Indeed, by letting h = f /( f dµ) in (12), we obtain In order to obtain a similar bound for the KL divergence D KL (ν||ν * r ), we need the following subspace logarithmic Sobolev inequalities. Theorem 2.10 (Subspace logarithmic Sobolev inequality). Let µ be a probability distribution on R d that satisfies Assumption 2.5 for some Γ ∈ R d×d and κ ≥ 1. Then the relation holds for any projector P r ∈ R d×d and for any continuously differentiable function h : supp(µ) → R such that Proof. The proof is given in Section 7.2.
Corollary 2.11. Let µ be a distribution which satisfies the subspace logarithmic Sobolev inequality (13) for some Γ ∈ R d×d and κ ≥ 1, and let ν be such that dν dµ ∝ f for some continuously differentiable function f such that ∇ log f 2 Γ −1 dµ < ∞. Then for any projector P r ∈ R d×d , the distribution ν * r such that dν * r dµ ∝ E µ (f |σ(P r )) satisfies where and Proof. The proof consists in rewriting (13) with h = (f /Z) 1/2 , where Z = E µ (f |σ(P r )) dµ = f dµ. We have ∇h = 1 2 (f /Z) 1/2 ∇ log f so that inequality (13) becomes holds. This concludes the proof. Corollary 2.11 provides an upper bound for D KL (ν||ν * r ) that is proportional to the reconstruction error R Γ (P r , H) = ∇ log f − P T r ∇ log f 2 Γ −1 dν, see (17). Minimizing this bound over P r corresponds to finding the projector P r so that P T r ∇ log f best approximates ∇ log f in the mean squared sense: this corresponds to principal component analysis of ∇ log f ; see [54]. The following proposition gives a closed form expression for a minimizer of P r → R Γ (P r , H) over the set of the rank-r projectors. It corresponds to Proposition 2.6 in [65] where Σ is replaced by Γ −1 . The proof relies on the Eckart-Young theorem, which provides the optimality properties of the singular value decomposition.
Proposition 2.12. Let Γ ∈ R d×d be a symmetric positive definite matrix and H ∈ R d×d be a symmetric positive semidefinite matrix. Denote by (λ i , v i ) ∈ R ≥0 × R d the i-th generalized eigenpair of the matrix pencil (H, Γ), meaning Hv i = λ i Γv i with λ i ≥ λ i+1 and v i Γ = 1, where · Γ = (·) T Γ(·). For any r ≤ d we have min Pr∈R d×d rank-r projector Furthermore, a solution to (18) is the following Γ-orthogonal projector (i.e., satisfying Corollary 2.11 and Proposition 2.12 ensure that, provided µ satisfies the subspace logarithmic Sobolev inequality (13) and provided P r is defined as in (19), the approximation ν * r of ν defined by This relation holds for any r ≤ d. Then for any ε ≥ 0, the choice r = r(ε) := min{r : κ Observe that a strong decay in the generalized eigenvalues of H ensures that r(ε) d. In particular, if H is rank deficient, we have λ i = 0 for all i > rank(H) so that r ≥ rank(H) implies D KL (ν||ν * r ) = 0. The spectrum of the matrix H reveals the low effective dimensionality of the posterior distribution: a strong decay in the generalized spectrum of H, or certainly a rank deficiency of H, ensures that there exists an approximation ν * r of ν such that D KL (ν||ν * r ) ≤ ε with small r(ε) d. Remark 2.13 (Prior-based dimension reduction). For any projector P r we can write In the case of a Gaussian prior, µ = N (m, Σ), we have Γ = Σ −1 (see Example 2.6) and so Proposition 2.12 ensures that the minimizer of P r → R I d (P T r , Γ −1 ) = R I d (P T r , Σ) is the orthogonal projector onto the leading eigenspace of the prior covariance Σ. This projector corresponds to the truncated Karhunen-Loève decomposition; see Section 4.1 for a detailed discussion.

An illustrative example
To illustrate how sharp the bound given by Corollary 2.11 can be, we consider a simple example for which the Kullback-Leibler divergence D KL (ν||ν * r ) is computable in closed form. This allows a comparison of the error D KL (ν||ν * r ) and its upper bound κ 2 R Γ (P r , H).
Assume that the prior µ is the standard Gaussian distribution N (0, I d ) and let the likelihood function be given by , so that µ satisfies Assumption 2.5 with Γ = I d and κ = 1. The posterior ν defined by dν dµ ∝ f is also Gaussian with zero mean and covariance Σ = (I d + A) −1 . In this setting, the matrix H defined by (16) can be written as follows Consider the generalized eigenvalue problem Hv i = λ i Γv i which, since Γ = I d , is simply the standard eigenvalue problem Hv i = λ i v i . Notice that H is a rational function in A, so that H and A have the same eigenvectors and their eigenvalues satisfy the relation According to Proposition 2.12, a projector minimizing P r → R Γ (P r , H) over the set of rank-r projectors is To analyze the sharpness of this inequality, we now compute D KL (ν||ν * r ). Using Proposition 2.4, we can express the conditional expectation E µ (f |σ(P r )) as follows: is Gaussian with zero mean and covariance Σ r . The Kullback-Leibler divergence from ν * r = N (0, Σ r ) to ν = N (0, Σ) admits the following closed form expression [41]: To continue the calculation, one needs the eigenvalues of Σ −1 r Σ. Let U be the orthogonal matrix containing the eigenvectors of A and let D = diag(α 1 , . . . , α d ) so that A = U DU T . By construction, From this relation, we deduce that the i-th eigenvalue of Σ −1 r Σ is 1 when i ≤ r and (1 + α i ) −1 otherwise. Substituting these eigenvalues into (22) yields We now analyze the deficit in inequality (20). With a Taylor expansion as α i goes to zero for all i ≥ r, we can write If the function f is nearly constant (α i ≈ 0) along the subspace span{v r+1 , . . . , v d }, then the upper bound (20) is close to 2 D KL (ν||ν * r ). In this particular example, the projector obtained by minimizing the upper bound in fact yields the optimal approximation of ν in Kullback-Leibler divergence, for any given rank r. This can be shown by Theorem 2.3 in [60]. Of course, minimizing the upper bound does not produce the optimal projector in general.

Building the approximation
In this section we propose and analyze algorithms for the numerical construction of a low-rank projector P r and of a profile function g, such that the distribution ν r given by dνr dµ ∝ g • P r is a controlled approximation of the posterior distribution ν. Recall that in the previous section we obtained the following decomposition of the Kullback-Leibler divergence: where dν * r dµ ∝ E µ (f |σ(P r )) gives the optimal (but intractable) probability distribution ν * r we can obtain for a given P r . The term D KL (ν||ν * r ) thus measures the error between the target ν and its best approximation for a given projector. Using the material presented in Section 2.2, Section 3.1 below shows how to build P r so that the error D KL (ν||ν * r )-or, more precisely, an upper bound for this error-is arbitrarily small. The second term D KL (ν * r ||ν r ) can be interpreted as a distance between g • P r and E µ (f |σ(P r )), so that the construction of ν r consists essentially in approximating a conditional expectation. In Section 3.2, we propose an approximation ν r = ν r that relies on a sample estimate F r of the conditional expectation E µ (f |σ(P r )). In the end, the computational strategy we adopt here is to: (i) construct P r to control the first term; and (ii) given such a projector, build the function g so that the second term is arbitrarily close to zero-or, at least, of the same order of magnitude as the first term.

Construction of the projector
Let us recall some of the results of Section 2.2. Under Assumption 2.5 and provided that f is sufficiently regular, Corollary 2.11 provides an upper bound for the first term D KL (ν||ν * r ) in the decomposition (23) by means of the reconstruction error R Γ (P r , This bound holds for any projector P r . We denote by P * r a rank-r projector which minimizes the reconstruction error, meaning By Proposition 2.12, P * r can be obtained by means of the generalized eigendecomposition of H. In practice, however, the matrix H may be difficult to compute exactly because it requires computing a high-dimensional integral when d 1. Instead, we consider the rank-r projector P r that is a minimizer of the approximate reconstruction error R Γ (P r , H), i.e., where H is a Monte Carlo approximation of H defined by Here, for the purpose of theoretical analysis, we assume that X 1 , . . . , X K are independent draws from the posterior distribution ν. Since drawing independent samples from the posterior is computationally infeasible for most problems, we will present in Section 3.3 an iterative algorithm that can be used to estimate H in practice. The minimization problem in (25) and (26) corresponds to the principal component analysis (PCA) of the random vector ∇ log f (X), X ∼ ν, and its solution follows directly from the solution of a generalized eigenvalue problem; see [40]. An important question is how large should K be in order to control the reconstruction error R Γ ( P r , H)? We refer to [7,54] for recent progress in this direction. The following proposition gives a new sufficient condition on K so that P r is a quasioptimal solution to (24). This result relies on concentration properties of sub-Gaussian random vectors; see [9,64]. The proof is given in Section 7.3.
Here · ψ 2 denotes the sub-Gaussian norm, meaning ξ ψ 2 = sup k≥1 k −1/2 E(|ξ| k ) 1/k for any real-valued random variable ξ. Let H be a Ksample Monte Carlo approximation of H. For any 0 < δ, η < 1, the condition for some absolute (numerical) constant Ω, ensures that with probability at least 1 − η, the following relation holds for any projector P r : Assumption (27) means that ∇ log f (X) is a sub-Gaussian vector. Intuitively it says that the tails of the distribution of ∇ log f (X) are at most Gaussian. Together with the independence of X 1 , . . . , X K , the sub-Gaussian property is an essential ingredient in the proof of Proposition 3.1. We give now examples for which (27) is satisfied.
u T Hu. This shows that relation (27) holds true for all u ∈ R d .

Example 3.3 (Unbounded gradients)
. We now consider an example where ∇ log f (X) satisfies (27) without being bounded. As in Section 2.3, let µ = N (0, I d ) be the standard normal prior and let f : Then Z is sub-Gaussian and Z ψ 2 ≤ Lσ Z holds for some absolute constant L; see, for instance, Example 5.8 in [64]. This means that (27) holds for all u ∈ R d with a constant L which does not depend on A.
Proposition 3.1 gives a sufficient condition for (29) to hold with high probability. This relation yields the following quasi-optimality result: Then, even though P r is not a minimizer of P r → R Γ (P r , H), the reconstruction error R Γ ( P r , H) is no greater than an arbitrary multiplicative constant 1+δ 1−δ ≥ 1 times the minimum of the reconstruction error R Γ (P * r , H). In particular R Γ (P * r , H) = 0 implies R Γ ( P r , H) = 0. Together with Corollary 2.11 and Proposition 2.12, the above quasi-optimality result allows us to deduce that, with probability at least 1 − η, we have where dν * r dµ ∝ E µ (f |σ(P r )) with P r = P r . This relation ensures that, up to a multiplicative constant 1+δ 1−δ , the bound on the error D KL (ν||ν * r ) decays with r at the same rate as if we had used the projector P r = P * r that minimizes the true reconstruction error. We note that condition (28) requires K to be at least proportional to the rank of H. If H is full rank and if d 1, then K ought to be large in order to satisfy (28). In practice, though, we observe that suitable projectors can be obtained by choosing K proportional to the rank r of P r (which is PCA's "rule of thumb"). Here there is a challenge in finding weaker conditions on K that nonetheless allow one to obtain a quasi-optimality result. As in [54], one should try to exploit some properties of H such as a rapid spectral decay or a large spectral gap.

Approximation of the conditional expectation
In this section we assume that a rank-r projector P r ∈ R d×d is given. We consider the problem of approximating the conditional expectation E µ (f |σ(P r )). Recall that, using the notation of Proposition 2.4, we have For any x ∈ R d , an approximation of E µ (f |σ(P r ))(x) can be obtained via a Monte Carlo estimate of the form 1 In general, the law of Ξ depends on x and we should redraw the samples Ξ 1 , . . . , Ξ M for each different value of x, which can be computationally demanding.
Instead, we will consider an approximation of E µ (f |σ(P r )) of the form where Y 1 , . . . , Y M are independent samples from the prior µ. The form of F r is motivated by situations where X r and X ⊥ are in fact independent under the prior, i.e., where p ⊥ (·|P r x) = p ⊥ (·). This occurs, for example, in the case of a Gaussian prior with covariance Γ −1 , since the projector P r is by construction Γ-orthogonal. In this case, F r is an unbiased estimator of E µ (f |σ(P r )). From a computational perspective, the samples Y i will be drawn once and reused for each evaluation of x → F r (x). We define ν r to be the random probability distribution defined by Notwithstanding the motivation above, the laws of (I d − P r )Y i and of U ⊥ Ξ i are in general different, so that F r (x) is in general a biased estimate of E µ (f |σ(P r ))(x). In this case we cannot hope for D KL (ν * r || ν r ) to go to zero with M . In order to analyze this bias, let us introduce the probability distribution µ such that for any Borel function h, where X and Y are independent random variables with distribution µ. Equation (35) is equivalent to saying that µ is the probability distribution of the random vector P r X + (I d − P r )Y . The following proposition is proven in Section 7.4.
By Proposition 3.4 we can write E( F r ) = E µ (f |σ(P r )) which allows us to interpret F r as a Monte Carlo estimate of E µ (f |σ(P r )). The following proposition gives a bound for the expectation of the error D KL (ν * r || ν r ). The proof is given in Section 7.5. Proposition 3.5. Let ν r be the distribution such that dν r dµ ∝ E µ (f |σ(P r )) and assume that f (x) > 0 for µ-a.e. x. We have Proposition 3.5 shows that, up to a third order term, the expectation of the error D KL (ν * r || ν r ) is bounded by the sum of two terms. The first term is the Kullback-Leibler divergence from ν r to ν * r and corresponds to the bias E( F r ) = E µ (f |σ(P r )). The second term can be interpreted as a measure of the variance of F r . Under some assumptions on the distribution µ, the following proposition provides an upper bound for those two terms.
Proposition 3.6. In addition to the assumptions of Proposition 3.5, assume that µ admits a Lebesgue density ρ ∝ exp(−V −Ψ), where Ψ is a bounded function such that exp(sup Ψ−inf Ψ) ≤ κ and where V : x → 1 2 x − m 2 Γ for some m ∈ R d and some symmetric positive-definite matrix Γ ∈ R d×d . Then, for any · Γ -orthogonal projector P r , we have and where The proof is given is Section 7.6. Proposition 3.6 requires µ to be a bounded perturbation of a Gaussian distribution N (m, Γ −1 ). This is a stronger assumption than Assumption 2.5. Neglecting the third order term in (36), Proposition 3.6 allows one to bound E(D KL (ν * r || ν r )) by (C 1 + C 2 M )E Γ (P r , f ) where C 1 and C 2 are two constants which depend only on κ. If µ is Gaussian (κ = 1) then C 1 = 0 and C 2 = 1/2. In that case, E(D KL (ν * r || ν r )) goes to zero with M . This is not surprising because µ = µ holds whenever µ is Gaussian, 1 so that ν * r = ν r and hence D KL (ν * r ||ν r ) = 0. In the general case κ = 1, (37) and (38) show that both the variance and the bias of F r are no greater than a constant independent of P r times E Γ (P r , f ). Note that the quantity E Γ (P r , f ) differs from the reconstruction error R Γ (P r , H) only by the term f E µ (f |σ(P r )) −1 ; compare (15) and (39). Provided that this term is uniformly bounded from above, we can write The above relation shows that the error E(D KL (ν * r || ν r )) can be controlled by the reconstruction error R Γ (P r , H), provided the supremum of f E µ (f |σ(P r )) −1 is finite. Then a small number of samples M will be enough to guarantee that D KL (ν * r || ν r ) is, in expectation, of the same order of magnitude as the first term D KL (ν||ν * r ) in the decomposition (23), which was our initial motivation. We now give an example for which the term sup (f /E µ (f |σ(P r ))) decreases to one with the rank of the projector. Example 3.7. As in Section 2.3, we again consider µ = N (0, I d ) and f : x → exp(− 1 2 x T Ax) for some symmetric matrix A 0. By applying Proposition 2.12, the projector that minimizes where v i is the i-th eigenvector of A and let Σ = (I d + A) −1 and Σ r = (I d +P r AP r ) −1 . Using the closed-form expression (21) for E µ (f |σ(P r )) and since A−P r AP r = where α i ≥ 0 is the i-th eigenvalue of A. This shows that, in this example, the supremum of f E µ (f |σ(P r )) −1 goes monotonically to one with the rank r of P r .
Remark 3.8. In the way they are presented, the error analyses of the projector and of the conditional expectation (i.e., the two terms on the right-hand side of (23)) are not unified: the first term is bounded in high probability (31), while the second is controlled in expectation (via Proposition 3.5 and subsequent discussion). One could unify these results, for instance by using a Markov inequality to control the probability that D KL (ν * r || ν r ) exceeds a certain value. However, we believe this step is not essential to the purpose of this section, which is to show that the associated sample approximations are feasible and sound.

Ideal algorithm
Algorithm 1 can be used to construct an approximation ν r of the posterior distribution ν. It assumes that we can draw samples from the posterior distribution ν. Since this is typically not possible in practice, this algorithm is called "ideal."

Algorithm 1 Ideal algorithm
Require: Likelihood function f , prior distribution µ, error threshold ε, sample sizes K and M . 1: Draw K independent samples X 1 , . . . , X K from dν ∝ f dµ 2: Compute ∇ log f (X k ) for k = 1, . . . , K 3: Assemble the matrix H = 1 : Compute a rank-r projector P r with the smallest rank such that R Γ (P r , H) ≤ ε 5: Draw M samples Y 1 , . . . , Y M from µ 6: Return the approximate distribution ν r defined by

Construction using approximate measures
Obtaining samples from the posterior distribution ν at step 1 of Algorithm 1 can be difficult in practice; indeed, this challenge is one of the motivations for the dimension reduction approach presented here. To alleviate this difficulty, we can construct H using samples from another distribution ν that can be directly simulated. By the definition (16) of H, we have so that, thanks to Corollary 2.11, we can write for any projector P r . Assuming ν satisfies sup dν d ν < +∞ (a common assumption in importance sampling [52,55]), the above relation suggests that minimizing the approximate reconstruction error P r → R Γ (P r , H) can yield suitable projectors. In particular if H is rank deficient, then R Γ (P r , H) = 0 and thus D KL (ν||ν * r ) = 0 for a suitable projector P r with r = rank( H). In general, however, the constant sup dν d ν will be unknown in practice, and we will not be able to use (40) as a bound for the error D KL (ν||ν * r ). By drawing samples from ν rather than from ν at step 1 of Algorithm 1, the matrix H computed at step 2 is a Monte Carlo approximation of H. Proposition 3.1 still applies when replacing ν by ν. Then for K sufficiently large (K = O(rank( H)), for instance), it holds with high probability that any rank-r projector that minimizes P r → R Γ (P r , H) will be a quasi-optimal solution to the minimization problem of P r → R Γ (P r , H); see the discussion of Section 3.1. Here we list two common choices of ν.
1. Laplace approximation. The Laplace approximation [59] constructs a Gaussian distribution ν = N ( m, Σ), where m is the mode of ν and the matrix Σ −1 is obtained from Hessian of the negative log density of ν evaluated at m. Note that if ν is Gaussian then ν is exactly ν. Even though Laplace's method is a simple way to obtain a Gaussian approximation of ν, there is no guarantee that sup dν d ν is finite in general.
2. Prior distribution. Sampling from the prior distribution µ is usually tractable. With the choice ν = µ, we have dν d ν = dν dµ = f f dµ . For most applications, the likelihood function f is bounded so that sup dν d ν < ∞. Note that this choice has been considered in [16].

Iterative construction
As suggested by [22], rather than limiting ourselves to a fixed approximation H, we can approximate the true H (16) using a sequential importance sampling framework. Let us consider a sequence of posterior approximations ν Notice that the same samples Y i are used for every 1 ≤ l ≤ L. The idea is to use ν so that is a self-normalized importance sampling estimator of H. Having computed H (l) , the (l + 1)-th projector P (l+1) r is defined as a projector with minimal rank such that the approximate reconstruction error R Γ (P (l+1) r , H (l) ) is below some prescribed tolerance. This iterative construction is detailed in Algorithm 2. Notice that at the first iteration, the importance weights w (0) k are set to one. We make this choice in order to avoid the potential degeneracy (i.e., large variance of the weights) that might occur when µ is a poor approximation of ν. Also, our implementation includes a constraint on the rank of the projector, so that it cannot exceed a user-defined maximum rank r max . By doing so, we avoid any explosion of the rank in the earlier stages of the algorithm, i.e., when a poor posterior approximation might yield a crude approximation of H. We emphasize that this algorithm only involves sampling from the low-dimensional posterior approximations { ν (l) r }.

Alternative approaches to dimension reduction 4.1 Karhunen-Loève-based dimension reduction
The Karhunen-Loève decomposition is a simple and powerful tool for reducing the dimension of a given random vector X ∈ R d . Letting m = E(X), this method exploits the fact that X − m may take values mostly on a low-dimensional subspace of R d , so that X − m can be well approximated by P r (X − m) for some low-rank projector P r ∈ R d×d . The standard approach is to seek P r such that the mean squared error E( (X − m) − P r (X − m) 2 2 ) is below some prescribed tolerance. We can write E( (X − m) − P r (X − m) 2 2 ) = trace((I d − P r )Σ(I d − P T r )), where Σ = E((X − m)(X − m) T ) is the covariance matrix of X. Using Proposition 2.12, we have that the orthogonal projector onto the r-dimensional leading eigenspace of Σ is a minimizer of the mean squared error over the set of rank-r projectors. Furthermore we have min Pr∈R d×d , rank-r projector where σ 2 i is the i-th eigenvalue of Σ. This relation shows that a strong decay in the spectrum of Σ ensures the mean squared error can be arbitrarily small for some r d. The eigenvectors of Σ are called the Karhunen-Loève modes of X and m + P r (X − m) corresponds to the truncated Karhunen-Loève decomposition of X.
This methodology can be used to approximate the prior distribution µ. Assuming X ∼ µ has mean m = E(X) and covariance matrix Σ, let µ r be the distribution of m + P r (X − m), where P r is the rank-r orthogonal projector onto the dominant eigenspace of the prior covariance matrix Σ. As proposed in [45,49,44], we can then introduce the approximate posteriorν r such that dν r dµ r ∝ f. Compute ∇ log f (X Draw K samples X 1 , . . . , X K from ν (l) r using any MCMC algorithm 8: Assemble the matrix

11:
Compute the lowest rank r ε such that R Γ (P, H (l) ) ≤ ε for some rank-r ε projector P

13:
Define the approximate distribution ν This Karhunen-Loève-based dimension reduction and the dimension reduction method primarily considered in this paper yield different approximation formats for the posterior measure: the latter considers approximate measures that are fully supported on supp(µ) (which is convex by Assumption 2.5), while the former seeks an approximate measure that is supported on supp(µ r ) = supp(µ) ∩ {m + Im(P r )}. The difference in support makes the two approximations hard to compare since, by construction, the divergence D KL (ν||ν r ) is infinite. This is hardly surprising, as Karhunen-Loève dimension reduction does not try to minimize D KL (ν||ν r ), but rather focuses on the mean squared error E( (X − m) − P r (X − m) 2 ), with X ∼ µ. From a computational perspective, Karhunen-Loève dimension reduction relies only on the prior measure µ: there is no need to compute expectations over the posterior measure ν or to evaluate gradients of the likelihood function, and hence the resulting approximation is easy to compute. Yet these computational advantages highlight some fundamental limitations of the method. The efficiency of the reduction strategy is limited to cases where there is a sharp decay in the eigenvalues of the prior covariance matrix Σ; the method does not exploit any low-dimensional structure that the likelihood function might have.
Remark 4.1. It is also possible to apply the Karhunen-Loève dimension reduction method to the posterior measure ν. In that case, we seek an approximation of ν defined as the distribution of m + P r (X − m) for some projector P r , where X ∼ ν follows the posterior distribution and m = E(X) is the posterior mean. As before, the projector can be defined as the minimizer of P r → E( (X − m) − P r (X − m) 2 2 ) for X ∼ ν, which turns out to be an orthogonal projector onto the leading eigenspace of the posterior covariance. This approach is no longer an a priori method since it requires computing an integral over the posterior distribution.

Likelihood-informed subspace
The likelihood-informed subspace (LIS) reduction method [21] leverages optimality results available for the Bayesian linear-Gaussian model [60] to propose a structure-exploiting approximation of the posterior distribution. The method assumes the prior to be Gaussian µ = N (m, Σ) and the likelihood function f , up to a multiplicative constant, to have the following form: Here G is a suitably regular forward operator, y is the observed data vector, and Σ obs is the covariance matrix of the observational noise which is assumed to be additive and Gaussian. The LIS reduction approximates the posterior ν by ν LIS r , defined as Here P r is a projector that is self-adjoint with respect to the inner product induced by Σ −1 and whose range is spanned by the leading generalized eigenvectors of the matrix pencil (H LIS , The matrix H LIS is the expectation over the posterior ν of the Gauss-Newton Hessian of the log-likelihood ∇G T Σ −1 obs ∇G. The work of [22] generalizes the construction of the matrix H LIS to reference measures other than the posterior-for example, replacing ν with the prior measure or the Laplace approximation to the posterior. The LIS methodology differs from the methodology considered in this paper in the following ways. First, the likelihood function is approximated by x → f (P r x + (I d − P r )m), rather than by the conditional expectation E µ (f |σ(P r )). This choice leads to a suboptimal approximation with respect to the Kullback-Leibler divergence; see Section 2.1. Notice however that this approximation choice is quite similar to the Monte Carlo approximation of the conditional expectation presented in Section 3.2. If the sample size satisfies M = 1, the approximation of the conditional expectation defined by (33) is F r : x → f (P r x + (I d − P r )Y ), where Y is a sample drawn from the prior distribution µ.
The other important difference is in the definition of the projector. Recall that in Section 2.2, P r is defined as the projector onto the leading eigenspace of the pencil (H, Σ −1 ), where H is defined in (16). The matrices H and H LIS are in general different, and so are the resulting projectors. The projector P r in the present paper is defined as the minimizer of an upper bound on the Kullback-Leibler divergence between ν and its approximation (see Propositon 2.12 in Section 2.2). In contrast, the projector introduced in [21] is only justified by analogy with optimality results developed in the linear-Gaussian case [60]. As a consequence, the LIS projector does not come with an error analysis on the resulting posterior approximation, while the strategy presented in this paper does. Yet both methodologies can perform remarkably well in applications, and can be even comparable; cf. Section 5.

Active subspace for Bayesian inverse problems
The active subspace (AS) method [15,58] is a dimension reduction technique which addresses the approximation (in the L 2 sense) of a real-valued function by means of a ridge function. These directions span the so-called active subspace [58]. In [16], this methodology is used to approximate the log-likelihood function in a Bayesian inverse problem. Denoting the prior covariance by Σ, the posterior distribution ν is approximated by where P r is defined as the Σ −1 -orthogonal projector onto the leading generalized eigenspace of the matrix pencil (H AS , Σ −1 ), with Active subspace reduction differs from the methodology introduced in this paper in the following aspects. First, the conditional expectation in (42) applies to the log-likelihood function rather than to the likelihood function itself. This choice is motivated by the fact that E µ (log f |σ(P r )) yields an optimal approximation of log f in the set {g • P r , g : R d → R} with respect to the L 2 µ -norm. However, the function exp E µ (log f |σ(P r )) is not optimal with respect to the Kullback-Leibler divergence; see Section 2.1.
Another difference is in the choice of projector. Comparing (43) with (16), we see that the integral in (43) is taken over the prior µ rather than the posterior ν, and thus the matrices H AS and H are in general different. Notice that H AS corresponds to the matrix H defined in (40) with the choice ν = µ.
Finally, we mention that the active subspace method comes with an upper bound on the Hellinger distance between ν and its approximation; see Theorem 3.1 in [16]. This analysis relies on a Poincaré inequality rather than on a logarithmic Sobolev inequality. Moreover, the Hellinger bound in [16] contains unknown (and uncontrolled) constants. In this paper, we provide sufficient conditions on the prior µ to control the constants associated with the logarithmic Sobolev inequality and hence with our Kullback-Leibler error bound. From a theoretical point of view, it is challenging to relate the two bounds, and we cannot make a definitive statement about their relative merit: comparing upper bounds of two different metrics is not really informative. Instead we will compare the two methodologies by means of numerical experiments in Section 5.

Numerical illustration
We use two Bayesian inverse problems to numerically demonstrate various theoretical aspects of the proposed dimension reduction method. In both examples we assume that the prior distribution is Gaussian, µ = N (m, Σ). This choice of prior distribution satisfies Assumption 2.5 with Γ = Σ −1 and κ = 1. In both examples, we suppose that the data y correspond to predictions of a nonlinear forward model x → G(x) that are corrupted by measurement noise, where the latter is additive and normally distributed with zero mean and covariance Σ obs . Thus we have a likelihood function, up to a multiplicative constant, of the form Computer codes used to produce these numerical results, along with instructions, are available at https://github.com/tcui001/certified.

Example 1: Atmospheric remote sensing
Our first example is a realistic atmospheric remote sensing problem, where satellite observations from the Global Ozone MOnitoring System (GOMOS) are used to estimate the concentration profiles of various gases in the atmosphere.

Problem setup
The GOMOS instrument repeatedly measures light intensities at different wavelengths λ and different altitudes "alt." The light transmissions T λ,alt are modeled using Beer's law: where the integral is taken along the ray path, and z(ζ) is the height in the atmosphere of a point ζ on the path. The curvature of the earth is taken into account; see the illustration in Figure 1. The quantity a gas λ (z), known from laboratory measurements, is called the cross-section. It is a measure of how much a gas absorbs light of a given wavelength λ at a given height z. In this model there are N gas gases, and for each we would like to infer the density profile κ gas : z → κ gas (z). Each of the N gas profiles are modeled as independent random processes with log-normal prior distributions. That is, log κ gas (·) follows a Gaussian distribution N (m gas , Σ gas ) where Σ gas is the covariance operator associated with the squared exponential kernel These priors are chosen to promote smooth gas density profiles with large variations. We consider a discretization of the vertical axis z into N alt layers with piecewise constant densities. We denote by x ∈ R N alt Ngas the vector containing the logarithms of the (unknown) N gas gas densities in the N alt layers, i.e., where z 1 , . . . , z N alt denote the N alt altitudes and vec(·) denotes the vectorization operator. Because κ gas is log-normally distributed, the prior distribution µ = N (m, Σ) of the parameter x is Gaussian, and its prior mean m and its prior covariance Σ are derived from the discretization of m gas and Σ gas .
The data y is a noisy measurement of the light transmissions T λ,alt at N λ different wavelengths and at N alt different altitudes. It is a vector y ∈ R N λ N alt given by y = G(x) + ε obs where the forward model G(x) results from the discretization of Beer's law (45), and where ε obs ∼ N (0, Σ obs ) is centered Gaussian noise with known observation covariance Σ obs . With this model, the likelihood function f takes the form of (44).
Here we adopt the same model setup and synthetic data set used in [21,22]. The atmosphere is discretized into N alt = 50 layers and with N gas = 4 profiles to infer so the total parameter dimension is d = N alt N gas = 200. We have observations at N λ = 1416 wavelengths, and thus the dimension of the data is N alt N λ = 70800. We refer the readers to [62,33] for a further description of the model setup, the data set, and the Bayesian treatment of this inverse problem.  Figure 1: The principle of the GOMOS measurement. The atmosphere is represented as spherical layers around the Earth. Note that the thickness of the layers is much larger relative to the Earth in this figure than in reality. The figure is adopted from [33], with the permission of the authors.

Comparison of approximations
We compare the posterior approximations obtained by the proposed method (see Sections 2 and 3) to those obtained by the existing methods described in Section 4. The comparison is summarized in Figures 2 and 3.
Except for the prior-based method described in Remark 2.13, which builds a projector only from the prior covariance matrix Σ, the projector P r is built with a generalized eigendecomposition of where ρ is either the prior µ, the posterior ν, or the Laplace approximation of the posterior ν.  [20,Proposal 3.4] is used as a state-of-the-art method for obtaining posterior samples. We first run the DILI algorithm for 10 5 iterations to adaptively estimate the matrix H (∇G) ν and the associated reduced dimensional subspace required internally within the algorithm. Then, we continue to run the DILI algorithm with the previously computed subspace for another 10 6 iterations to generate posterior samples. This choice is not essential; we could equivalently have used many other algorithms to generate a reference posterior sample set-for instance, an adaptive Metropolis method [11,34], a delayed-acceptance method [12,17,18], Hamiltonian Monte Carlo [26,50] and improved versions thereof [29,35], or sequential Monte Carlo samplers [25].
Once we have estimated the matrix H for ρ ∈ {µ, ν, Laplace(ν)}, we compute the generalized eigendecomposition of the estimate (i.e., ( H , Σ −1 )) and assemble the rank-r projector P r onto the corresponding leading eigenspace for r ∈ {10, 15, . . . , 50}. Then we approximate the conditional expectation E µ (f |σ(P r )) by with either M = 10 prior samples (yellow line with squares in Figure 2) or with M = 1 prior sample (solid grey lines). We also consider the deterministic approximation x → f ( P r x+(I d − P r )m), which corresponds to using M = 1 sample that is fixed to the prior mean, Y 1 = m (red line with crosses). Finally, for each of the projectors P r obtained above, we also calculate and plot (blue line with circles) the reconstruction error 1 2 R Γ ( P r , H ref ) (15), where H ref is an independent Monte Carlo estimate of H = (∇ log f )(∇ log f ) T dν, again computed using 10 6 posterior samples generated by the DILI algorithm. As described in Corollary 2.11, the reconstruction error is an upper bound for the KL divergence D KL (ν||ν * r ) achieved by any given projector. Once both the projector and the approximation of the conditional expectation are defined, the KL divergence D KL (ν|| ν r ) is approximated with an independent set of 10 6 samples drawn from the posterior using the DILI algorithm [20]. This approximation is assumed to be sufficiently accurate for the purpose of our experiments.
Among the different ways of approximating the conditional expectation, the sample-based approach with M = 10 outperforms both the one-sample approach (M = 1) and the deterministic approach (Y 1 = m). As shown in the theoretical analysis of Section 3.2, the error due to approximating the conditional expectation is, in expectation, proportional to the reconstruction error; this is why we observe that the error decays quickly with the rank of the projector even when M is small. Note that the cost for evaluating the approximation of the conditional expectation is proportional to M : it is exactly M times the cost of evaluating the exact likelihood function f , which might be expensive when M is large. Therefore we would prefer using the deterministic prior mean approach (Y 1 = m) which, in addition to having a low evaluation cost, has an accuracy comparable to that of the sample-based approach with M = 10 and seems better than the one-sample approach (M = 1). Figure 2 also shows that the reconstruction error provides a good estimate of the KL divergence, a quantity which is rarely available in practice. In theory, the reconstruction error gives a bound on the KL divergence for any projector, provided that the conditional expectation is computed exactly-which is not what we do here. In practice, however, we observe that the reconstruction error is in general a fairly good error indicator even when using a sample approximation of the conditional expectation. This might not be true when sample sizes are very small; see, for instance, the cases ρ = ν and M = 1 in Figure 2.
In Figure 3 we compare the performance of the different definitions of the projector. Here we only consider the KL divergence of posterior approximations obtained using the sample-based approximation of the conditional expectation with M = 10 (left) and the prior mean approximation Y 1 = m (right). The projector obtained from H (∇ log f ) ν achieves the best accuracy overall. The LIS method, obtained from H (∇ log G) ν also performs reasonably well. We observe a performance gap between these two projectors and the other five: the latter lead to significantly less accurate posterior approximations, at any given rank.

Demonstration of Algorithm 2
Now we demonstrate the iterative procedure given in Algorithm 2 for constructing projectors and posterior approximations. In this exercise, we approximate the conditional expectation using the prior mean option (Y 1 = m). For each iteration l, the rank of the projector, the reconstruction error R Γ (P  Figure 4. The left column of Figure 4 illustrates the iterative procedure where the error threshold ε = 10 −2 is used to determine the rank of the projector in each iteration. Here we also set the maximum rank of the projector to be r max = 40. The right column of Figure 4 shows a fixed-rank variant of the iterative procedure. In this variant, the rank of the projector is held constant at r = 30. In both cases, we use the adaptive MALA algorithm [2] to draw samples from the r-dimensional distribution ν (l) r at each iteration. This choice is incidental; any other standard MCMC algorithm could have been used instead.
The iterative procedure appears to be effective in approximating the posterior. When the rank of the projector is dynamically adjusted (left column of Figure 4), the iterative procedure achieves the desired level of accuracy within the first two iterations; then the rank of the projector, the reconstruction error, and the KL divergence are stabilized in later iterations. When the rank of the projector is fixed (right column of Figure 4), the KL divergence decays significantly in the first two iterations, and then is stabilized in later iterations.

Impact of sample size
We now analyze the impact of the sample size K in the Monte Carlo approximation H of H; see (26). As discussed in Section 3.1, Proposition 3.1 ensures that for sufficiently large K, the reconstruction error is sufficiently well approximated so that the resulting P r in (25) is quasi-optimal with high probability. We now illustrate numerically the impact of the sample size K via the following criteria: 1. Is R Γ ( P r , H) close to the minimum of the true reconstruction error R Γ (P * r , H)? Since our primary goal is to minimize a bound on the KL divergence, it makes sense to measure the quality of P r in terms of its ability to minimize the reconstruction error. Thus we do not look at distances between P r and the minimizer P * r (such as the operator norm P r − P * r or the Frobenius norm P r − P * r F ) because these distances are not directly related to the problem of posterior approximation.
2. Can we use R Γ ( P r , H) as an error indicator for the KL divergence? Since the quantity R Γ ( P r , H) is the only one that is accessible in practice, we would like to know whether it is safe to use as an error estimator.
Numerical results are summarized in the left plot of Figure 5 for different sample sizes K ∈ {50, 100, 200, 500, 1000, 2000, 10 6 } generated by the DILI algorithm. To assess the variability of these results due to the randomness in H, we also repeat the previous experiment 10 times for K ∈ {50, 200, 1000, 10 6 } (right plot of Figure 5). Concerning the first criterion, we observe that the quality of P r depends not only on the sample size K, but also on the rank r of the projector. Indeed, R Γ ( P r , H ref ) (solid colored lines in Figure  5) is closest to the minimum of the reconstruction error (black line) for small r and large K. This suggests that the sample size should be chosen larger when r is chosen to be large. Note that this fact is not revealed by the theoretical analysis of Section 3.1.
Concerning the second criterion, we observe that the approximate reconstruction error R Γ ( P r , H) (dashed colored lines) is always smaller than the reference reconstruction error R Γ ( P r , H ref ) (solid colored lines). This relationship is particularly apparent when the sample size K is small and the rank r is large. This means that one should be careful when using R Γ ( P r , H) in place of the exact reconstruction error, as it tends to underestimate the error.
Overall, in this example, both criteria are validated for sufficiently large sample sizes. With K ≥ 1000, the resulting projector P r has a reconstruction error comparable to that of the projector P * r , and the approximate reconstruction error R Γ ( P r , H) provides a fairly good estimate for R Γ ( P r , H).

Example 2: Elliptic PDE
Our second example is an inverse problem aiming at estimating the spatially inhomogeneous coefficient of an elliptic PDE, adopted from [22]. In physical terms, our problem setup corresponds to inferring the transmissivity field of a two-dimensional groundwater aquifer from partial observations of the stationary drawdown field of the water table, measured from well bores.  The discretized transitivity field T (ζ) is endowed with a log-normal prior distribution, i.e.,

Problem setup
where the prior mean is set to log(1000 [m/day]) and the inverse of the covariance matrix Σ −1 is defined through the discretization of an Laplace-like stochastic partial differential equation [47], where W(ζ) is white noise. In this example, we set the stationary, anisotropic correlation tensor K to and put κ = 50. The "true" transmissivity field is a realization from the prior distribution. The true transmissivity field, sources/sinks, simulated drawdown field, and synthetic data are shown in Figure 6. Partial observations of the pressure field are collected at 13 sensors whose locations are depicted by black dots in Figure 6(c), with additive error e ∼ N (0, σ 2 I 13 ). The standard deviation σ of the measurement noise is prescribed so that the observations have signal-to-noise ratio 20. The noisy data are shown in Figure 6(d).
In this example, the finite element discretization of (48) uses 120 × 40 bilinear elements to represent the drawndown field p(ζ), while the transmissivity field T (ζ) is modeled as piecewise constant for each element. This yields the discretized forward model G : R 4800 → R 13 , where the parameter x is of dimension d = 4800.

Numerical results
We proceed with the same comparison as in the GOMOS example. We compare the KL divergence D KL (ν|| ν r ) for different posterior approximations ν r obtained using the various projectors given in Sections 3 and 4. The projectors used in this comparison are built in the same way as in the GOMOS example. Also we consider the same three options (M = 1, M = 10, and Y 1 = m) for approximating the conditional expectation within ν r . Again, to compute the KL divergence D KL (ν|| ν r ), we draw 10 6 samples from the posterior using the DILI algorithm, which follows the same setup as that of the previous example. The comparison is summarized in Figures 7 and 8. In Figure 7 we observe that the approximation of the conditional expectation with M = 10 samples outperforms the two other options, M = 1 and Y 1 = m. This holds true regardless of how the projector is constructed. We also notice that the prior-mean approximation Y 1 = m in general performs better than the one-sample approximation M = 1. Similar to the GOMOS example, the reconstruction error R Γ (P r , H ref ) provides an effective error indicator for posterior approximations that employ the estimate of the conditional expectation with M = 10 for any choice of projector P r . Also, the KL divergence decays with r at the same rate as the reconstruction error, independently of the choice of approximation scheme for the conditional expectation. with ρ = ν outperforms all the other methods. In other wordsand as observed in the previous example-the method proposed in Section 2 provides more effective dimension reduction than the alternatives. Note also that the projectors obtained from H (∇G) ρ with ρ = ν and with ρ = Laplace(ν) perform better than the projector obtained from H (∇ log f ) ρ with ρ = Laplace(ν), indicating that efficient dimension reduction can also be obtained from gradient of the forward model itself. All other projectors have a rather large accuracy gap relative to the abovementioned projectors. In particular, choosing ρ = µ in this example seems to be much less effective than averaging over the posterior distribution or its Laplace approximation.

Conclusion
We have addressed the problem of reducing the dimension of a Bayesian inverse problem, in the nonlinear/non-Gaussian setting. A Bayesian inverse problem has a low intrinsic dimension when the update from the prior distribution to the posterior distribution is essentially low-dimensional, meaning that the data only inform a few directions in the parameter space. We proposed a methodology that reveals and exploits such structure by seeking an approximation of the likelihood as a ridge function, i.e., a function that varies only on a low-dimensional subspace of its input space. To obtain this approximation, we first identified the optimal profile function of the ridge approximation, meaning the profile that minimizes the Kullback-Leibler divergence from the posterior approximation to the exact posterior. Then, using logarithmic Sobolev inequalities, we derived an upper bound for the remaining error. The bound admits a simple form and can be easily minimized: this is how the informed directions are discovered. Moreover, our methodology provides a computable upper bound for the approximation error in the posterior distribution, measured in terms of Kullback-Leibler divergence, a quantity that is otherwise difficult to compute in practice.
Our method is fundamentally gradient-based and requires the second moment matrix of the gradient of log-likelihood function. Computing this matrix can be challenging in practice because it entails integrating over the posterior distribution. We thus propose several sample-based approximation schemes, including an iterative algorithm which employs only low-dimensional posterior approximations. One direction which remains to be explored is the optimal balance between the computational effort used to reduce the dimension of the problem and the computational effort of exploring the final dimension-reduced posterior distribution.
The use of logarithmic Sobolev inequalities requires some (rather strong) assumptions on the prior distribution, e.g., being Gaussian, being a Gaussian mixture, or being a bounded perturbation of a Gaussian. One open question left for future work is how to weaken the assumption on the prior so that the resulting theory applies to priors with tails that are heavier than Gaussian tails. Another natural improvement of the methodology is the extension to infinite-dimensional parameter spaces.
Analytical and numerical examples demonstrate good performance of the proposed method. In particular, we show that it outperforms other state-of-the-art gradient-based dimension reduction schemes. We also note that effective MCMC algorithms for large-scale Bayesian inverse problems, such as the DILI sampler of [20], fundamentally rely on the low-dimensional structure of the posterior distribution. In its original form, [20] uses the LIS method to discover this low-dimensional structure, a method which is less efficient than the one we propose here. Incorporating our new developments in the DILI algorithm (or in any other MCMC algorithm that exploits dimension reduction) could thus yield better sampling performance. Making this precise-e.g., understanding how the quality of the subspace Im(P r ) affects MCMC mixing rates-is another important topic for future work. We denote by U r ∈ R d×r a matrix whose columns form a basis of Im(P r ) and we let U ∈ R d×d be defined as the horizontal concatenation of U r and U ⊥ . We have U ξ = U r ξ r + U ⊥ ξ ⊥ for any ξ ∈ R d , where ξ r ∈ R r and ξ ⊥ ∈ R d−r are the vectors containing respectively the r first and the d − r last components of ξ.
By definition of p ⊥ (ξ ⊥ |P r x) and since P r U r = U r , this function satisfies for any ξ r ∈ R r . Let h be σ(P r )-measurable function. By Lemma 2.2, h can be written as the composition of a function with P r , so that h(U r ξ r + U ⊥ ξ ⊥ ) = h(U r ξ r ) holds for any ξ r ∈ R r and ξ ⊥ ∈ R d−r . Thus we have where |U | denotes the determinant of U . This shows that F satisfies (4). Since it is a σ(P r )measurable function (as the composition of a function with P r ), then it is the conditional expectation E µ (f |σ(P r )).

Proof of Theorem 2.10
In this proof, let K = supp(µ). Assumption 2.5 implies dµ dx = ρ ∝ e −(V +Ψ) where V : K → R and Ψ : K → R satisfy respectively (10) and (11). Let P r ∈ R d×d be a rank-r projector and let U ⊥ ∈ R d×(d−r) be a matrix whose columns form a basis of Ker(P r ). For any x ∈ K, the conditional density ξ ⊥ → p ⊥ (ξ ⊥ |P r x) defined by (8) and Ψ ⊥ (ξ ⊥ ) = Ψ(P r x + U ⊥ ξ ⊥ ) for any ξ ⊥ ∈ K ⊥ (x). Because K is convex, the set K ⊥ (x) is also convex. The function V ⊥ inherits the convexity property of V : using the chain rule, we have for any ξ ⊥ ∈ K ⊥ (x). Also, Ψ ⊥ is a bounded function which satisfies We conclude that the distribution on the convex set K ⊥ (x) with Lebesgue density p ⊥ (·|P r x) satisfies Assumption 2.5. Then it satisfies the logarithmic Sobolev inequality (12), meaning that holds for any function g with sufficient regularity, where the norm · * is defined by v 2 Inequality (49) holds in particular when the function g is defined by g(ξ ⊥ ) = h(P r x + U ⊥ ξ ⊥ ) for any ξ ⊥ ∈ K ⊥ (x), where h is any function on R d with sufficient regularity. By the chain rule we have ∇g(ξ ⊥ ) = U T ⊥ ∇h(P r x + U ⊥ ξ ⊥ ) so that (49) can be written as Note that the denominator in the logarithm is the conditional expectation E µ (h 2 |σ(P r )) evaluated at x, see Proposition 2.4. Let U r ∈ R d×r be a matrix whose columns form a basis of Im(P r ), and let ξ r ∈ K r (x) = {ξ r ∈ R r : U r ξ r + (I d − P r )x ∈ K} ⊂ R r . Replacing x by U r ξ r in the above inequality and using the definition of the conditional density (8), we obtain By integrating over ξ r we obtain To conclude the proof, it remains to show that U T ⊥ ∇h 2 where · = (·) T (·) denotes the Euclidean norm of R d and Γ 1/2 denotes a symmetric positivedefinite square root of Γ. This concludes the proof.

Proof of Proposition 3.1
The proof of Proposition 3.1 requires concentration properties of sub-Gaussian random vectors. We will need the following lemma, which is essentially Theorem 5.39 in [64].
Lemma 7.1. Let Z be an isotropic sub-Gaussian random vector in R n , meaning that E(ZZ T ) = I n , and let L be a constant such that w T Z ψ 2 := sup p≥1 p −1/2 E(|w T Z| p ) 1/p ≤ L w , for any w ∈ R n , where · denotes the Euclidean norm of R n . Let Σ = 1 K K k=1 Z k Z T k , where Z 1 , . . . , Z K are K independent copies of Z. Then for every t ≥ 0, with probability at least 1 − 2 exp(−ct 2 ), we have where Ω = log(9)/c 1 . Then (28) implies (53) so that (52) holds with probability greater than 1 − η, which concludes the proof of Proposition 3.1.

Proof of Proposition 3.4
Let F : x → E(f (P r x + (I d − P r )Y )) where Y ∼ µ and let h be any σ(P r )-measurable function. By Doob-Dynkin's lemma, see Lemma 2.2, h can be written as a composition of a function with P r , so that h(P r x + (I d − P r )y) = h(x) holds for all x, y ∈ R d . F also satisfies the same property. Using these observations and using the definition of µ , we can write F h dµ = F (P r x + (I d − P r )y) h(P r x + (I d − P r )y) µ(dx)µ(dy) = F (x)h(x) µ(dx) = f (P r x + (I d − P r )y)µ(dy) h(x) µ(dx) = f (P r x + (I d − P r )y) h(P r x + (I d − P r )y) µ(dx)µ(dy) = f h dµ , which shows that F is the conditional expectation E µ (f |σ(P r )).

Proof of Proposition 3.6
To prove this proposition, we will need the following lemmas and corollaries.
As first observed in [56], a distribution which satisfies the logarithmic Sobolev inequality also satisfies the Poincaré inequality; see for instance Proposition 5.1.3 in [4]. The next proposition shows that the subspace logarithmic Sobolev inequality (13) implies the subspace Poincaré inequality.
We now have all the material to prove Proposition 3.6.