# Coupled input-output dimension reduction

Qiao Chen, Élise Arnaud, Ricardo Baptista, Olivier Zahm

► To cite this version:

Qiao Chen, Élise Arnaud, Ricardo Baptista, Olivier Zahm. Coupled input-output dimension reduction: Application to goal-oriented bayesian experimental design and global sensitivity analysis. 2024. hal-04614897

HAL Id: hal-04614897

https://hal.science/hal-04614897

Preprint submitted on 19 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
COUPLED INPUT-OUTPUT DIMENSION REDUCTION:
APPLICATION TO GOAL-ORIENTED BAYESIAN EXPERIMENTAL
DESIGN AND GLOBAL SENSITIVITY ANALYSIS

QIAO CHEN*, ÉLISE ARNAUD*, RICARDO BAPTISTA†, AND OLIVIER ZAHM*

**Abstract.** We introduce a new method to jointly reduce the dimension of the input and output space of a high-dimensional function. Choosing a reduced input subspace influences which output subspace is relevant and vice versa. Conventional methods focus on reducing either the input or output space, even though both are often reduced simultaneously in practice. Our coupled approach naturally supports goal-oriented dimension reduction, where either an input or output quantity of interest is prescribed. We consider, in particular, goal-oriented sensor placement and goal-oriented sensitivity analysis, which can be viewed as dimension reduction where the most important output or, respectively, input *components* are chosen. Both applications present difficult combinatorial optimization problems with expensive objectives such as the expected information gain and Sobol' indices. By optimizing gradient-based bounds, we can determine the most informative sensors and most sensitive parameters as the largest diagonal entries of some diagnostic matrices, thus bypassing the combinatorial optimization and objective evaluation.

**Key words.** dimension reduction, Poincaré inequality, Cramér-Rao inequality, derivative-based, goal-oriented, Bayesian optimal experimental design, global sensitivity analysis

**MSC codes.** 65D40, 62F15, 62K05

**1. Introduction.** High dimensionality poses a challenge in many fields of applied mathematics, seriously limiting the effectiveness of established solution methods. Tasks like approximation, sampling, and optimization become exponentially harder as the size of the respective exploration spaces explodes with growing dimensions. And yet, the dimensionality of computational problems continues to rise due to improved processing capacities that allow for higher resolution simulations and an unprecedented accumulation of data. In response to this challenge, dimension reduction methods seek to leverage the underlying low-dimensional structures present in many problems.

Classical dimension reduction methods target either the input or the output space of a model. For example, truncated Karhunen-Loève expansion [37] is often applied to reduce input parameters, while the reduced basis method [5], proper orthogonal decomposition [48, 39] and principal component analysis (PCA) [28] are used on the output states. Derivative-based reduction methods like ours have been mostly proposed for the parameter space like the active subspace [17, 53] and likelihood informed subspace [18, 54] methods. Although input and output reduction methods are often applied simultaneously, the two spaces are commonly treated separately (see e.g. [4, 15, 41]). However, choosing an input subspace inevitably impacts the relevant output subspace and vice versa. Taking this interplay into account can significantly reduce computational cost by allowing lower-dimensional approximations with the same level of accuracy [4].

Few works consider the coupling between the input and output space reduction, though the coupling is often only algorithmic and one-directional. For example, [36] uses a greedy reduced basis method for output state reduction and re-purposes the

<footnote>* Université Grenoble Alpes, Inria Grenoble, France (elise.arnaud@univ-grenoble-alpes.fr, qiao.chen@inria.fr, olivier.zahm@inria.fr).</footnote>

<footnote>† Computing and Mathematical Sciences, California Institute of Technology, Pasadena, CA (rsb@caltech.edu).</footnote>
COUPLED INPUT-OUTPUT DIMENSION REDUCTION:
APPLICATION TO GOAL-ORIENTED BAYESIAN EXPERIMENTAL
DESIGN AND GLOBAL SENSITIVITY ANALYSIS

QIAO CHEN*, ÉLISE ARNAUD*, RICARDO BAPTISTA†, AND OLIVIER ZAHM*

**Abstract.** We introduce a new method to jointly reduce the dimension of the input and output space of a high-dimensional function. Choosing a reduced input subspace influences which output subspace is relevant and vice versa. Conventional methods focus on reducing either the input or output space, even though both are often reduced simultaneously in practice. Our coupled approach naturally supports goal-oriented dimension reduction, where either an input or output quantity of interest is prescribed. We consider, in particular, goal-oriented sensor placement and goal-oriented sensitivity analysis, which can be viewed as dimension reduction where the most important output or, respectively, input *components* are chosen. Both applications present difficult combinatorial optimization problems with expensive objectives such as the expected information gain and Sobol' indices. By optimizing gradient-based bounds, we can determine the most informative sensors and most sensitive parameters as the largest diagonal entries of some diagnostic matrices, thus bypassing the combinatorial optimization and objective evaluation.

**Key words.** dimension reduction, Poincaré inequality, Cramér-Rao inequality, derivative-based, goal-oriented, Bayesian optimal experimental design, global sensitivity analysis

**MSC codes.** 65D40, 62F15, 62K05

**1. Introduction.** High dimensionality poses a challenge in many fields of applied mathematics, seriously limiting the effectiveness of established solution methods. Tasks like approximation, sampling, and optimization become exponentially harder as the size of the respective exploration spaces explodes with growing dimensions. And yet, the dimensionality of computational problems continues to rise due to improved processing capacities that allow for higher resolution simulations and an unprecedented accumulation of data. In response to this challenge, dimension reduction methods seek to leverage the underlying low-dimensional structures present in many problems.

Classical dimension reduction methods target either the input or the output space of a model. For example, truncated Karhunen-Loève expansion [37] is often applied to reduce input parameters, while the reduced basis method [5], proper orthogonal decomposition [48, 39] and principal component analysis (PCA) [28] are used on the output states. Derivative-based reduction methods like ours have been mostly proposed for the parameter space like the active subspace [17, 53] and likelihood informed subspace [18, 54] methods. Although input and output reduction methods are often applied simultaneously, the two spaces are commonly treated separately (see e.g. [4, 15, 41]). However, choosing an input subspace inevitably impacts the relevant output subspace and vice versa. Taking this interplay into account can significantly reduce computational cost by allowing lower-dimensional approximations with the same level of accuracy [4].

Few works consider the coupling between the input and output space reduction, though the coupling is often only algorithmic and one-directional. For example, [36] uses a greedy reduced basis method for output state reduction and re-purposes the

<footnote>* Université Grenoble Alpes, Inria Grenoble, France (elise.arnaud@univ-grenoble-alpes.fr, qiao.chen@inria.fr, olivier.zahm@inria.fr).</footnote>

<footnote>† Computing and Mathematical Sciences, California Institute of Technology, Pasadena, CA (rsb@caltech.edu).</footnote>
snapshot parameters to span a reduced parameter space. Conversely, [19] employs likelihood informed parameter subspaces, and then reduces the output using PCA on snapshots sampled from a reduced parameter distribution. The operator learning framework of [38] optimizes the output reduction (seen here as a linear decoder parametrized by neural networks) after fixing a reduced input, e.g., using Karhunen-Loève or active subspaces; see also the discussion in [31]. [12] derive error bounds for joint input-output dimension reduction which permits one to balance the contribution of the two reductions to meet an overall level of accuracy.

The coupling between input and output space is also the motivation behind goal-oriented methods. These approaches tailor the dimension reduction of the input or output space to specific quantities of interest in the opposite space to enhance computational efficiency and significance of the results. For example, optimal sensor placement can be viewed as dimension reduction that selects the most informative output components. The goal-oriented optimal sensor placement problem [21, 9, 51] seeks optimality with respect to a prescribed lower dimensional parameter of interest (instead of the whole parameter). On the other hand, sensitivity analysis can be viewed as input dimension reduction where parameter components that cause the most (or least) output variations are identified. Goal-oriented sensitivity analysis considers sensitivities with respect to a specified output of interest.

**1.1. Contribution and Outline.** We develop a coupled dimension reduction method that accounts for the interdependence between the input and output space of a non-linear function $G: \mathbb{R}^d \to \mathbb{R}^m$. As a consequence of the coupling, our method naturally supports goal-oriented dimension reduction. Using the Poincaré and a Cramér-Rao-like inequality, we establish gradient-based upper and lower bounds for several joint and goal-oriented dimension reduction objectives. These bounds provide us with computable error estimates and an easy optimization algorithm based on the eigen-decomposition of two diagnostic matrices

$$H_X(V_s) = \mathbb{E}[\nabla G(X)^\top V_s V_s^\top \nabla G(X)],$$

$$H_Y(U_r) = \mathbb{E}[\nabla G(X) U_r U_r^\top \nabla G(X)^\top],$$

for the input and output space, respectively. Here, $\nabla G(X) \in \mathbb{R}^{m \times d}$ denotes the Jacobian of $G$, while $U_r \in \mathbb{R}^{d \times r}$ and $V_s \in \mathbb{R}^{m \times s}$ are orthogonal matrices spanning the reduced input and output space, respectively. We note that each diagnostic matrix is a function of a specified subspace in the opposing space. In the goal-oriented case where a fixed $U_r$ prescribes a parameter of interest $U_r^\top X$, the optimal $V_s^*(U_r)$ consists of the dominant eigenvectors of $H_Y(U_r)$. Equivalently, for a given output of interest $V_s^\top G(X)$, the optimal $U_r^*(V_s)$ consists of the dominant eigenvectors of $H_X(V_s)$. To obtain coupled subspaces $(U_r^*, V_s^*)$, we propose an alternating eigendecomposition of both diagnostic matrices.

In the goal-oriented optimal sensor placement scenario, we seek $V_\tau(U_r) = [e_\tau^m]$ consisting of canonical basis vectors $e_i^m \in \mathbb{R}^m$ indexed by $\tau \subseteq \{1, \dots, m\}$. By maximizing a derivative-based lower bound of the expected information gain, we obtain the optimal $\tau^*$ as the index set of the largest diagonal entries of $H_Y(U_r)$. Similarly, in the goal-oriented sensitivity analysis scenario, we seek $U_\tau(V_r) = [e_\tau^d]$ consisting of canonical basis vectors $e_i^d \in \mathbb{R}^d$. The solution $\tau^* \subseteq \{1, \dots, d\}$ that optimizes a derivative-based bound of the total Sobol' index is given as the largest diagonal entries of $H_X(V_s)$. In both applications, we circumvent the expensive evaluation of the objective functions and the combinatorial optimization problem by maximizing our derivative-based bounds.
snapshot parameters to span a reduced parameter space. Conversely, [19] employs likelihood informed parameter subspaces, and then reduces the output using PCA on snapshots sampled from a reduced parameter distribution. The operator learning framework of [38] optimizes the output reduction (seen here as a linear decoder parametrized by neural networks) after fixing a reduced input, e.g., using Karhunen-Loève or active subspaces; see also the discussion in [31]. [12] derive error bounds for joint input-output dimension reduction which permits one to balance the contribution of the two reductions to meet an overall level of accuracy.

The coupling between input and output space is also the motivation behind goal-oriented methods. These approaches tailor the dimension reduction of the input or output space to specific quantities of interest in the opposite space to enhance computational efficiency and significance of the results. For example, optimal sensor placement can be viewed as dimension reduction that selects the most informative output components. The goal-oriented optimal sensor placement problem [21, 9, 51] seeks optimality with respect to a prescribed lower dimensional parameter of interest (instead of the whole parameter). On the other hand, sensitivity analysis can be viewed as input dimension reduction where parameter components that cause the most (or least) output variations are identified. Goal-oriented sensitivity analysis considers sensitivities with respect to a specified output of interest.

**1.1. Contribution and Outline.** We develop a coupled dimension reduction method that accounts for the interdependence between the input and output space of a non-linear function $G: \mathbb{R}^d \to \mathbb{R}^m$. As a consequence of the coupling, our method naturally supports goal-oriented dimension reduction. Using the Poincaré and a Cramér-Rao-like inequality, we establish gradient-based upper and lower bounds for several joint and goal-oriented dimension reduction objectives. These bounds provide us with computable error estimates and an easy optimization algorithm based on the eigen-decomposition of two diagnostic matrices

$$H_X(V_s) = \mathbb{E}[\nabla G(X)^\top V_s V_s^\top \nabla G(X)],$$

$$H_Y(U_r) = \mathbb{E}[\nabla G(X) U_r U_r^\top \nabla G(X)^\top],$$

for the input and output space, respectively. Here, $\nabla G(X) \in \mathbb{R}^{m \times d}$ denotes the Jacobian of $G$, while $U_r \in \mathbb{R}^{d \times r}$ and $V_s \in \mathbb{R}^{m \times s}$ are orthogonal matrices spanning the reduced input and output space, respectively. We note that each diagnostic matrix is a function of a specified subspace in the opposing space. In the goal-oriented case where a fixed $U_r$ prescribes a parameter of interest $U_r^\top X$, the optimal $V_s^*(U_r)$ consists of the dominant eigenvectors of $H_Y(U_r)$. Equivalently, for a given output of interest $V_s^\top G(X)$, the optimal $U_r^*(V_s)$ consists of the dominant eigenvectors of $H_X(V_s)$. To obtain coupled subspaces $(U_r^*, V_s^*)$, we propose an alternating eigendecomposition of both diagnostic matrices.

In the goal-oriented optimal sensor placement scenario, we seek $V_\tau(U_r) = [e_\tau^m]$ consisting of canonical basis vectors $e_i^m \in \mathbb{R}^m$ indexed by $\tau \subseteq \{1, \dots, m\}$. By maximizing a derivative-based lower bound of the expected information gain, we obtain the optimal $\tau^*$ as the index set of the largest diagonal entries of $H_Y(U_r)$. Similarly, in the goal-oriented sensitivity analysis scenario, we seek $U_\tau(V_r) = [e_\tau^d]$ consisting of canonical basis vectors $e_i^d \in \mathbb{R}^d$. The solution $\tau^* \subseteq \{1, \dots, d\}$ that optimizes a derivative-based bound of the total Sobol' index is given as the largest diagonal entries of $H_X(V_s)$. In both applications, we circumvent the expensive evaluation of the objective functions and the combinatorial optimization problem by maximizing our derivative-based bounds.
The remainder of this article is organized as follows. We begin in Section 2 by introducing several coupled and goal-oriented dimension reduction objectives. Section 3 contains our main results. We derive derivative-based upper and lower bounds for the input-output dimension-reduction objectives, as well as an algorithm that computes the optimal coupled subspaces with respect to the bounds. In Section 4, we consider two goal-oriented applications. In particular, we look at goal-oriented Bayesian optimal experimental design in Subsection 4.1 and goal-oriented global sensitivity analysis in Subsection 4.2. We finally demonstrate our theoretical results on numerical examples in Section 5.

**2. Problem Setting and Objectives.** Let $G: \mathbb{R}^d \to \mathbb{R}^m$ be a deterministic non-linear function of a random input parameter vector $X$. We reduce the input vector $X$ and output $Y = G(X)$ by retaining only some low-dimensional projections

$$X_r = U_r^\top X,$$

$$Y_s = V_s^\top Y,$$

where $U_r \in \mathbb{R}^{d \times r}$, $V_s \in \mathbb{R}^{m \times s}$ are two matrices with $r \ll d$ and $s \ll m$ orthonormal columns. The identification of $U_r$ and $V_s$ depends on the objective at hand. In the following, we introduce four potential objectives. The first two consider general function approximation and Bayesian inference problems, while the second two consider goal-oriented applications of our framework.

In reduced order modelling (ROM), $G$ is typically a parametrized partial differential equation (PDE) model taking parameter $X$ and returning the (discretized) solution $G(X)$ of the PDE. ROM aims to construct a fast-to-evaluate approximation $\tilde{G}$ that can be used in place of $G$ for real-time or multiple-query problems such as control, optimization or sampling. A possible formalization of the problem is to find matrices $U_r$ and $V_s$ which minimize the $L^2$-error

$$(O1) \qquad \min_{U_r, V_s} \mathbb{E}[\|G(X) - \tilde{G}(X)\|^2],$$

$$(2.1) \qquad \text{where} \quad \tilde{G}(x) = V_s \tilde{g}(U_r^\top x) + \tilde{G}_\perp,$$

for some optimal low-dimensional function $\tilde{g}: \mathbb{R}^r \to \mathbb{R}^s$ and vector $\tilde{G}_\perp \in \mathbb{R}^m$ to be specified later in Proposition 3.3. If the minimum value of the objective (O1) is close to 0, then $G \approx \tilde{G}$ is essentially constant along the $(d-r)$-dimensional subspace $\text{Ker}(U_r^\top)$ and the variation in its output is confined in the $m$-dimensional affine subspace $\text{Im}(V_s) + \tilde{G}_\perp$.

In Bayesian inference, the aim is to compute statistics for the posterior distribution $\pi_{X|Y}$. Existing knowledge on $X$, represented by the prior distribution $\pi_X$, is thereby updated with new information from a noisy observation $Y = G(X) + \eta$, where $\eta \sim N(0, \sigma^2 I_m)$ denotes some additive Gaussian noise with $I_m \in \mathbb{R}^{m \times m}$ being the identity matrix. While one may consider more general likelihood functions, we will restrict ourselves to the setting of additive and Gaussian observational noise in this work. By this definition, the likelihood $\pi_{Y|X}$ is Gaussian so that Bayes theorem yields the posterior density, which is known up to a normalization constant, as

$$\pi_{X|Y}(x|y) \propto \exp\left(-\frac{1}{2\sigma^2} \|y - G(x)\|^2\right) \pi_X(x).$$

In the case of a high-dimensional $X$, standard Markov-chain Monte Carlo (MCMC) methods need an excessive amount of iterations (and thus of model evaluations) to
The remainder of this article is organized as follows. We begin in Section 2 by introducing several coupled and goal-oriented dimension reduction objectives. Section 3 contains our main results. We derive derivative-based upper and lower bounds for the input-output dimension-reduction objectives, as well as an algorithm that computes the optimal coupled subspaces with respect to the bounds. In Section 4, we consider two goal-oriented applications. In particular, we look at goal-oriented Bayesian optimal experimental design in Subsection 4.1 and goal-oriented global sensitivity analysis in Subsection 4.2. We finally demonstrate our theoretical results on numerical examples in Section 5.

**2. Problem Setting and Objectives.** Let $G: \mathbb{R}^d \to \mathbb{R}^m$ be a deterministic non-linear function of a random input parameter vector $X$. We reduce the input vector $X$ and output $Y = G(X)$ by retaining only some low-dimensional projections

$$X_r = U_r^\top X,$$

$$Y_s = V_s^\top Y,$$

where $U_r \in \mathbb{R}^{d \times r}$, $V_s \in \mathbb{R}^{m \times s}$ are two matrices with $r \ll d$ and $s \ll m$ orthonormal columns. The identification of $U_r$ and $V_s$ depends on the objective at hand. In the following, we introduce four potential objectives. The first two consider general function approximation and Bayesian inference problems, while the second two consider goal-oriented applications of our framework.

In reduced order modelling (ROM), $G$ is typically a parametrized partial differential equation (PDE) model taking parameter $X$ and returning the (discretized) solution $G(X)$ of the PDE. ROM aims to construct a fast-to-evaluate approximation $\tilde{G}$ that can be used in place of $G$ for real-time or multiple-query problems such as control, optimization or sampling. A possible formalization of the problem is to find matrices $U_r$ and $V_s$ which minimize the $L^2$-error

$$(O1) \qquad \min_{U_r, V_s} \mathbb{E}[\|G(X) - \tilde{G}(X)\|^2],$$

$$(2.1) \qquad \text{where} \quad \tilde{G}(x) = V_s \tilde{g}(U_r^\top x) + \tilde{G}_\perp,$$

for some optimal low-dimensional function $\tilde{g}: \mathbb{R}^r \to \mathbb{R}^s$ and vector $\tilde{G}_\perp \in \mathbb{R}^m$ to be specified later in Proposition 3.3. If the minimum value of the objective (O1) is close to 0, then $G \approx \tilde{G}$ is essentially constant along the $(d-r)$-dimensional subspace Ker($U_r^\top$) and the variation in its output is confined in the $m$-dimensional affine subspace $\text{Im}(V_s) + \tilde{G}_\perp$.

In Bayesian inference, the aim is to compute statistics for the posterior distribution $\pi_{X|Y}$. Existing knowledge on $X$, represented by the prior distribution $\pi_X$, is thereby updated with new information from a noisy observation $Y = G(X) + \eta$, where $\eta \sim \mathcal{N}(0, \sigma^2 I_m)$ denotes some additive Gaussian noise with $I_m \in \mathbb{R}^{m \times m}$ being the identity matrix. While one may consider more general likelihood functions, we will restrict ourselves to the setting of additive and Gaussian observational noise in this work. By this definition, the likelihood $\pi_{Y|X}$ is Gaussian so that Bayes theorem yields the posterior density, which is known up to a normalization constant, as

$$\pi_{X|Y}(x|y) \propto \exp\left(-\frac{1}{2\sigma^2} \|y - G(x)\|^2\right) \pi_X(x).$$

In the case of a high-dimensional $X$, standard Markov-chain Monte Carlo (MCMC) methods need an excessive amount of iterations (and thus of model evaluations) to
accurately sample the posterior. Reducing the dimension of X aims at identifying the components $X_r$ which are most *informed* by the data Y in order to run the inference algorithm, e.g., MCMC on a low-dimensional subspace [18]. Reducing the dimension of Y permits one to only consider the most *informative* data $Y_s$, typically yielding smaller autocorrelation between posterior samples and further reducing the number of required MCMC iterations [12]. We formalize this problem as finding $U_r$ and $V_s$ that minimize the posterior error with respect to the expected Kullback-Leibler (KL) divergence

$$ (O2) \qquad \min_{U_r, V_s} \mathbb{E}_Y [\mathbb{D}_{\text{KL}} (\pi_{X|Y} \| \tilde{\pi}_{X|Y})], $$

$$ (2.2) \qquad \text{where } \tilde{\pi}_{X|Y}(x|y) \propto \pi_{Y_s|X_r}(y_s|x_r)\pi_X(x), $$

with $\mathbb{D}_{\text{KL}}(\nu\|\mu) = \int \log(\nu/\mu) d\nu$. For (O2) to be close to 0, the parameter X needs to be essentially independent of the observations in the subspace $\text{Ker}(V_s)$, and observations of Y should primarily provide information within the parameter subspace $\text{Im}(U_r)$. To sample from the approximate density $\tilde{\pi}_{X|Y}$, we can independently sample $X_r \in \mathbb{R}^r$, which depends on the projected observation (see algorithms to sample the reduced posterior in [12, Section 6]), and its orthogonal complement $X_{\perp} \in \mathbb{R}^{d-r}$, which only depends on the prior.

Bayesian optimal experimental design (BOED) aims to select a set of sensors $\tau \subseteq \{1, \dots, m\}$ with $|\tau| = s$ that maximizes the information gained about the parameter that is to be inferred. Measurements from the sensors $\tau$ are gathered in the vector $Y_\tau = V_\tau^\top Y$, where $V_\tau = [e_\tau^m] \in \mathbb{R}^{m \times s}$ contains the canonical basis vectors $e_\tau^m \in \mathbb{R}^m$ indexed by $\tau$. The information content of the experimental design $V_\tau$ can be quantified by the expected information gain (EIG), which is defined as the KL-divergence between posterior and prior in expectation over the data $\Phi(V_\tau) = \mathbb{E}_Y [\mathbb{D}_{\text{KL}}(\pi_{X|Y_\tau} \| \pi_X)]$. When the inference problem focuses on a specific parameter of interest $X_r = U_r^\top X$, the sensor placement can be improved by maximizing the *goal-oriented EIG*

$$ (O3) \qquad \max_{V_\tau} \Phi(V_\tau|U_r) = \mathbb{E}_Y [\mathbb{D}_{\text{KL}} (\pi_{X_r|Y_\tau} \| \pi_{X_r})], $$

where $\pi_{X_r|Y_\tau}$ and $\pi_{X_r}$ are the corresponding marginalized densities. BOED is notoriously challenging as each evaluation of the EIG requires a costly nested Monte-Carlo estimator [47]. Existing methods often resort to linear Gaussian models for which the EIG can be computed explicitly [3]. Nevertheless, it remains a combinatorial optimization problem that typically only allows for sub-optimal greedy solutions [34, 16].

Global sensitivity analysis (GSA) aims to identify which set $\tau \subseteq \{1, \dots, d\}$ of $r = |\tau|$ input variables causes the least output variations in order to fix them to a nominal value and simplify the model. A widely used sensitivity index is the total Sobol' index $S_\tau^{\text{tot}} = 1 - \frac{\text{Tr}(\text{Cov}(G(X)|X_{-\tau}))}{\text{Tr}(\text{Cov}(G(X)))}$ with $X_{-\tau} = U_{-\tau}^\top X \in \mathbb{R}^{d-r}$ where $U_{-\tau} = [e_{-\tau}^d] \in \mathbb{R}^{d \times (d-r)}$ contains the canonical basis vectors $e_{-\tau}^d \in \mathbb{R}^d$ indexed by the complement of $\tau$ [24]. A smaller $S_\tau^{\text{tot}}$ implies that the output is less sensitive to the parameters indexed by $\tau$. When there are particular output directions of interest $V_s^\top G(X)$, it is more suitable to consider minimizing the *goal-oriented total Sobol' index*

$$ (O4) \qquad \min_{U_\tau} S_\tau^{\text{tot}}(U_\tau|V_s) = 1 - \frac{\text{Tr}(\text{Cov}(V_s^\top G(X)|U_{-\tau}^\top X))}{\text{Tr}(\text{Cov}(V_s^\top G(X)))}, $$

where $U_\tau = [e_\tau^d] \in \mathbb{R}^{d \times r}$. GSA faces similar challenges as BOED. First, the objective function is expensive to evaluate and the development of fast computational methods is an ongoing research topic [23, 43]. Secondly, identifying the optimal set that
accurately sample the posterior. Reducing the dimension of X aims at identifying the components $X_r$ which are most *informed* by the data Y in order to run the inference algorithm, e.g., MCMC on a low-dimensional subspace [18]. Reducing the dimension of Y permits one to only consider the most *informative* data $Y_s$, typically yielding smaller autocorrelation between posterior samples and further reducing the number of required MCMC iterations [12]. We formalize this problem as finding $U_r$ and $V_s$ that minimize the posterior error with respect to the expected Kullback-Leibler (KL) divergence

$$ (O2) \qquad \min_{U_r, V_s} \mathbb{E}_Y [\mathbb{D}_{\text{KL}} (\pi_{X|Y} \| \tilde{\pi}_{X|Y})], $$

$$ (2.2) \qquad \text{where } \tilde{\pi}_{X|Y}(x|y) \propto \pi_{Y_s|X_r}(y_s|x_r)\pi_X(x), $$

with $\mathbb{D}_{\text{KL}}(\nu\|\mu) = \int \log(\nu/\mu) d\nu$. For (O2) to be close to 0, the parameter X needs to be essentially independent of the observations in the subspace $\text{Ker}(V_s)$, and observations of Y should primarily provide information within the parameter subspace $\text{Im}(U_r)$. To sample from the approximate density $\tilde{\pi}_{X|Y}$, we can independently sample $X_r \in \mathbb{R}^r$, which depends on the projected observation (see algorithms to sample the reduced posterior in [12, Section 6]), and its orthogonal complement $X_{\perp} \in \mathbb{R}^{d-r}$, which only depends on the prior.

Bayesian optimal experimental design (BOED) aims to select a set of sensors $\tau \subseteq \{1, \dots, m\}$ with $|\tau| = s$ that maximizes the information gained about the parameter that is to be inferred. Measurements from the sensors $\tau$ are gathered in the vector $Y_{\tau} = V_{\tau}^{\top} Y$, where $V_{\tau} = [e_{\tau}^{m}] \in \mathbb{R}^{m \times s}$ contains the canonical basis vectors $e_{i}^{m} \in \mathbb{R}^{m}$ indexed by $\tau$. The information content of the experimental design $V_{\tau}$ can be quantified by the expected information gain (EIG), which is defined as the KL-divergence between posterior and prior in expectation over the data $\Phi(V_{\tau}) = \mathbb{E}_Y [\mathbb{D}_{\text{KL}}(\pi_{X|Y_{\tau}} \| \pi_X)]$. When the inference problem focuses on a specific parameter of interest $X_r = U_r^{\top} X$, the sensor placement can be improved by maximizing the *goal-oriented EIG*

$$ (O3) \qquad \max_{V_{\tau}} \Phi(V_{\tau}|U_r) = \mathbb{E}_Y [\mathbb{D}_{\text{KL}} (\pi_{X_r|Y_{\tau}} \| \pi_{X_r})], $$

where $\pi_{X_r|Y_{\tau}}$ and $\pi_{X_r}$ are the corresponding marginalized densities. BOED is notoriously challenging as each evaluation of the EIG requires a costly nested Monte-Carlo estimator [47]. Existing methods often resort to linear Gaussian models for which the EIG can be computed explicitly [3]. Nevertheless, it remains a combinatorial optimization problem that typically only allows for sub-optimal greedy solutions [34, 16].

Global sensitivity analysis (GSA) aims to identify which set $\tau \subseteq \{1, \dots, d\}$ of $r = |\tau|$ input variables causes the least output variations in order to fix them to a nominal value and simplify the model. A widely used sensitivity index is the total Sobol' index $S_{\tau}^{\text{tot}} = 1 - \frac{\text{Tr}(\text{Cov}(G(X)|X_{-\tau}))}{\text{Tr}(\text{Cov}(G(X)))}$ with $X_{-\tau} = U_{-\tau}^{\top} X \in \mathbb{R}^{d-r}$ where $U_{-\tau} = [e_{-\tau}^{d}] \in \mathbb{R}^{d \times (d-r)}$ contains the canonical basis vectors $e_{i}^{d} \in \mathbb{R}^{d}$ indexed by the complement of $\tau$ [24]. A smaller $S_{\tau}^{\text{tot}}$ implies that the output is less sensitive to the parameters indexed by $\tau$. When there are particular output directions of interest $V_s^{\top} G(X)$, it is more suitable to consider minimizing the *goal-oriented total Sobol' index*

$$ (O4) \qquad \min_{U_{\tau}} S_{\tau}^{\text{tot}}(U_{\tau}|V_s) = 1 - \frac{\text{Tr}(\text{Cov}(V_s^{\top} G(X)|U_{-\tau}^{\top} X))}{\text{Tr}(\text{Cov}(V_s^{\top} G(X)))}, $$

where $U_{\tau} = [e_{\tau}^{d}] \in \mathbb{R}^{d \times r}$. GSA faces similar challenges as BOED. First, the objective function is expensive to evaluate and the development of fast computational methods is an ongoing research topic [23, 43]. Secondly, identifying the optimal set that
minimizes the total Sobol' index poses a combinatorial problem. In practice, only single-element sets $\tau = \{i\}$ are considered and Sobol' indices are computed for all $i = 1, \dots, d$ to rank them by importance.

**3. Coupled Input-Output Dimension Reduction.** The following Lemmas allow us to address both the posterior approximation problem (O2) as well as the goal-oriented BOED problem (O3) by solving the $L^2$-approximation problem of the forward model (O1). The proofs can be found in Appendices A.1 and A.2.

**LEMMA 3.1.** With the above notations, for any $U_r, V_s$ and $\tilde{G}$ as in (2.1), we have

$$ (3.1) \qquad \mathbb{E}_Y [\mathbb{D}_{KL}(\pi_{X|Y} \| \tilde{\pi}_{X|Y})] \le \frac{1}{2\sigma^2} \mathbb{E}[\|G(X) - \tilde{G}(X)\|^2]. $$

**LEMMA 3.2.** For any $U_r, V_s$ and $\tilde{G}$ as in (2.1), we have

$$ (3.2) \qquad \Phi(V_s|U_r) \ge \mathbb{E}_Y [\mathbb{D}_{KL}(\pi_{X|Y} \| \pi_X)] - \frac{1}{2\sigma^2} \mathbb{E}[\|G(X) - \tilde{G}(X)\|^2]. $$

Inequality (3.1) shows that minimizing the $L^2$-error of $\tilde{G}$ corresponds to minimizing an upper bound of the posterior error of $\tilde{\pi}_{X|Y}$. At the same time, minimizing the $L^2$-error also maximizes a lower bound of the goal-oriented EIG $\Phi(V_s|U_r)$ according to (3.2). For this reason, we will focus on the $L^2$-approximation of $\tilde{G}$ in this section. The link to the goal-oriented GSA problem (O4) will be addressed later in Proposition 4.1. For any fixed pair of matrices $(U_r, V_s)$, the following proposition gives analytical expressions for the optimal function $\tilde{G}$ in (2.1). The proof is given in Appendix B.

**PROPOSITION 3.3.** Let $X$ be a random vector taking values in $\mathbb{R}^d$ and let $G: \mathbb{R}^d \to \mathbb{R}^m$ be such that $\mathbb{E}[\|G(X)\|^2] < \infty$. For any matrices $U_r \in \mathbb{R}^{d \times r}$ and $V_s \in \mathbb{R}^{m \times s}$ with orthogonal columns, we define $G^*: \mathbb{R}^d \to \mathbb{R}^m$ as

$$ G^*(x) = V_s g^*(U_r^\top x) + G_\perp^*, $$

where $g^*(x_r) = V_s^\top \mathbb{E}[G(X)|U_r^\top X = x_r]$ and $G_\perp^* = (I_m - V_s V_s^\top) \mathbb{E}[G(X)]$. Then $G^*$ minimizes $\tilde{G} \mapsto \mathbb{E}[\|G(X) - \tilde{G}(X)\|^2]$ over the set of square integrable functions in the form of (2.1). Furthermore, it holds that

$$ (3.3) \qquad \mathbb{E}[\|G(X) - G^*(X)\|^2] = \mathrm{Tr}(\mathrm{Cov}(G(X))) - \mathrm{Tr}(V_s^\top \mathrm{Cov}(\mathbb{E}[G(X)|U_r^\top X]) V_s), $$
where $\mathrm{Cov}(\mathbb{E}[G(X)|U_r^\top X]) \in \mathbb{R}^{m \times m}$ denotes the covariance of $\mathbb{E}[G(X)|U_r^\top X]$.

Proposition 3.3 allows us to obtain optimal $U_r$ and $V_s$ by minimizing the $L^2$-error of the forward model as given by (3.3), the minimum value for the objective in (O1). After that, one can approximate the low-dimensional function $g^*$ using state-of-the-art approximation methods and compute $G_\perp^*$ to obtain a reduced order model of $G$. We focus on reducing the inputs and/or outputs by minimizing (3.3) over $U_r$ and/or $V_s$, which is equivalent to

$$ (3.4) \qquad \max \mathrm{Tr} (V_s^\top \mathrm{Cov} (\mathbb{E}[G(X)|U_r^\top X]) V_s). $$

Let us look at some computational aspects of (3.4). Consider the goal-oriented case with fixed $U_r$. The matrix $V_s^*(U_r)$ that maximizes (3.4) spans the $s$-dimensional subspace which most effectively captures the output variations caused by parameter of interest $X_r = U_r^\top X$. It can be computed as the matrix containing the $s$ dominant eigenvectors of $\mathrm{Cov}(\mathbb{E}[G(X)|U_r^\top X])$, see Lemma 3.9 below. Without input dimension reduction $r = d$, this would correspond to the PCA solution which consists of
minimizes the total Sobol' index poses a combinatorial problem. In practice, only single-element sets $\tau = \{i\}$ are considered and Sobol' indices are computed for all $i = 1, \dots, d$ to rank them by importance.

**3. Coupled Input-Output Dimension Reduction.** The following Lemmas allow us to address both the posterior approximation problem (O2) as well as the goal-oriented BOED problem (O3) by solving the $L^2$-approximation problem of the forward model (O1). The proofs can be found in Appendices A.1 and A.2.

**LEMMA 3.1.** With the above notations, for any $U_r, V_s$ and $\tilde{G}$ as in (2.1), we have

$$ (3.1) \qquad \mathbb{E}_Y [\mathbb{D}_{KL}(\pi_{X|Y} \| \tilde{\pi}_{X|Y})] \le \frac{1}{2\sigma^2} \mathbb{E}[\|G(X) - \tilde{G}(X)\|^2]. $$

**LEMMA 3.2.** For any $U_r, V_s$ and $\tilde{G}$ as in (2.1), we have

$$ (3.2) \qquad \Phi(V_s|U_r) \ge \mathbb{E}_Y [\mathbb{D}_{KL}(\pi_{X|Y} \| \pi_X)] - \frac{1}{2\sigma^2} \mathbb{E}[\|G(X) - \tilde{G}(X)\|^2]. $$

Inequality (3.1) shows that minimizing the $L^2$-error of $\tilde{G}$ corresponds to minimizing an upper bound of the posterior error of $\tilde{\pi}_{X|Y}$. At the same time, minimizing the $L^2$-error also maximizes a lower bound of the goal-oriented EIG $\Phi(V_s|U_r)$ according to (3.2). For this reason, we will focus on the $L^2$-approximation of $\tilde{G}$ in this section. The link to the goal-oriented GSA problem (O4) will be addressed later in Proposition 4.1. For any fixed pair of matrices $(U_r, V_s)$, the following proposition gives analytical expressions for the optimal function $\tilde{G}$ in (2.1). The proof is given in Appendix B.

**PROPOSITION 3.3.** Let $X$ be a random vector taking values in $\mathbb{R}^d$ and let $G: \mathbb{R}^d \to \mathbb{R}^m$ be such that $\mathbb{E}[\|G(X)\|^2] < \infty$. For any matrices $U_r \in \mathbb{R}^{d \times r}$ and $V_s \in \mathbb{R}^{m \times s}$ with orthogonal columns, we define $G^*: \mathbb{R}^d \to \mathbb{R}^m$ as

$$ G^*(x) = V_s g^*(U_r^\top x) + G_\perp^*, $$

where $g^*(x_r) = V_s^\top \mathbb{E}[G(X)|U_r^\top X = x_r]$ and $G_\perp^* = (I_m - V_s V_s^\top) \mathbb{E}[G(X)]$. Then $G^*$ minimizes $\tilde{G} \mapsto \mathbb{E}[\|G(X) - \tilde{G}(X)\|^2]$ over the set of square integrable functions in the form of (2.1). Furthermore, it holds that

$$ (3.3) \qquad \mathbb{E}[\|G(X) - G^*(X)\|^2] = \mathrm{Tr}(\mathrm{Cov}(G(X))) - \mathrm{Tr}(V_s^\top \mathrm{Cov}(\mathbb{E}[G(X)|U_r^\top X]) V_s), $$
where $\mathrm{Cov}(\mathbb{E}[G(X)|U_r^\top X]) \in \mathbb{R}^{m \times m}$ denotes the covariance of $\mathbb{E}[G(X)|U_r^\top X]$.

Proposition 3.3 allows us to obtain optimal $U_r$ and $V_s$ by minimizing the $L^2$-error of the forward model as given by (3.3), the minimum value for the objective in (O1). After that, one can approximate the low-dimensional function $g^*$ using state-of-the-art approximation methods and compute $G_\perp^*$ to obtain a reduced order model of $G$. We focus on reducing the inputs and/or outputs by minimizing (3.3) over $U_r$ and/or $V_s$, which is equivalent to

$$ (3.4) \qquad \max \mathrm{Tr} (V_s^\top \mathrm{Cov}(\mathbb{E}[G(X)|U_r^\top X]) V_s). $$

Let us look at some computational aspects of (3.4). Consider the goal-oriented case with fixed $U_r$. The matrix $V_s^*(U_r)$ that maximizes (3.4) spans the $s$-dimensional subspace which most effectively captures the output variations caused by parameter of interest $X_r = U_r^\top X$. It can be computed as the matrix containing the $s$ dominant eigenvectors of $\mathrm{Cov}(\mathbb{E}[G(X)|U_r^\top X])$, see Lemma 3.9 below. Without input dimension reduction $r = d$, this would correspond to the PCA solution which consists of
the dominant eigenvectors of Cov(G(X)) [28]. Replacing G(X) by its conditional expectation E[G(X)|U_r^T X] is a natural modification to couple the output dimension reduction with a prescribed input subspace U_r. However, computing the covariance of the conditional expectation is challenging. It poses a common problem in global sensitivity analysis and requires dedicated algorithms such as the pick-freeze method [49]. The alternative goal-oriented problem with fixed V_s gives us a matrix U_r^*(V_s) that extracts the inputs which best explain the projected output V_r^T G(X). This problem is even more difficult to solve as there exists no closed form expression of U_r^*(V_s) in general. The same computational difficulties affect the problem of finding the optimal coupled matrix pair (U_r^*, V_s^*) that maximizes (3.4). In the following section, we derive a gradient-based bound that provides a feasible optimization procedure to identify the reduced subspaces.

**3.1. Derivative-based Bounds.** In this subsection, we derive derivative-based bounds for the L²-error (3.3) which permit more manageable optimization problems compared to (3.4). To derive the bounds, we invoke the Poincaré inequality [42] and a Cramér-Rao-like inequality [44] for the input random variable X.

**DEFINITION 3.4 (Poincaré and Cramér-Rao inequality).** Given a random vector X taking values in $\mathbb{R}^d$, we denote the smallest and the largest constants by $C(X) \ge 0$ and $c(X) \ge 0$, respectively, such that

$$ (3.5) \qquad c(X) \|\mathbb{E}[\nabla f(X)]\|^2 \le \mathbb{E}[(f(X) - \mathbb{E}[f(X)])^2] \le C(X) \mathbb{E}[\|\nabla f(X)\|^2], $$

holds for any continuously differentiable function $f: \mathbb{R}^d \to \mathbb{R}$. If $C(X) < \infty$, the right inequality is called the Poincaré inequality with constant $C(X)$ and, if $c(X) > 0$, the left inequality is called the Cramér-Rao inequality with constant $c(X)$.

It is well known that a Gaussian random vector $X \sim N(\mu, \Sigma)$ on $\mathbb{R}^d$ satisfies (3.5) with $c(X) = \lambda_{\min}(\Sigma)$ and $C(X) = \lambda_{\max}(\Sigma)$, see Lemma D.1 and [6]. More details regarding the Poincaré inequality can be found in [10, 54] and references therein. For the Cramér-Rao inequality see the discussion in Appendix D. We defer to Proposition 3.6 for sufficient conditions so that $c(X) > 0$ and $C(X) < \infty$ hold. We now state our main result. The proof can be found in Appendix C.

**THEOREM 3.5.** Given a random vector X taking values in $\mathbb{R}^d$, we let $\bar{C}(X) \ge 0$ and $\bar{c}(X) \ge 0$ be the smallest and the largest constant, respectively, such that

$$ (3.6) \qquad C(X_{\perp} | X_r = x_r) \le \bar{C}(X), $$

$$ (3.7) \qquad c(X_{\perp} | X_r = x_r) \ge \bar{c}(X), $$

holds for all $x_r \in \mathbb{R}^r$ and all $U_r \in \mathbb{R}^{d \times r}$ with orthonormal columns, where $X_r = U_r^T X$ and $X_{\perp} = U_{\perp}^T X$ with $U_{\perp} \in \mathbb{R}^{d \times (d-r)}$ being an orthogonal complement to $U_r$. In other words, we require every conditional distribution of $X_{\perp}$ to satisfy the Poincaré and Cramér-Rao inequality.

Let $G: \mathbb{R}^d \to \mathbb{R}^m$ be a continuously differentiable function such that $E[\|\nabla G(X)\|_F^2] < \infty$, where $\nabla G(X) \in \mathbb{R}^{m \times d}$ denotes the Jacobian of G and $\|\cdot\|_F$ is the Frobenius norm. Then, for any $U_r \in \mathbb{R}^{d \times r}$ and $V_s \in \mathbb{R}^{m \times s}$ with orthonormal columns, we have

$$ (3.8) \qquad \mathbb{E}[\|G(X) - G^*(X)\|^2] \le \bar{C}(X) (\mathbb{E}[\|\nabla G(X)\|_F^2] - \mathbb{E}[\|V_s^\top \nabla G(X) U_r\|_F^2]), $$

$$ (3.9) \qquad \mathbb{E}[\|G(X) - G^*(X)\|^2] \ge \bar{c}(X) (\|\mathbb{E}[\nabla G(X)]\|_F^2 - \|V_s^\top \mathbb{E}[\nabla G(X)] U_r\|_F^2), $$

with $G^*$ defined as in Proposition 3.3.
the dominant eigenvectors of Cov(G(X)) [28]. Replacing G(X) by its conditional expectation E[G(X)|U_r^T X] is a natural modification to couple the output dimension reduction with a prescribed input subspace U_r. However, computing the covariance of the conditional expectation is challenging. It poses a common problem in global sensitivity analysis and requires dedicated algorithms such as the pick-freeze method [49]. The alternative goal-oriented problem with fixed V_s gives us a matrix U_r^*(V_s) that extracts the inputs which best explain the projected output V_r^T G(X). This problem is even more difficult to solve as there exists no closed form expression of U_r^*(V_s) in general. The same computational difficulties affect the problem of finding the optimal coupled matrix pair (U_r^*, V_s^*) that maximizes (3.4). In the following section, we derive a gradient-based bound that provides a feasible optimization procedure to identify the reduced subspaces.

**3.1. Derivative-based Bounds.** In this subsection, we derive derivative-based bounds for the L²-error (3.3) which permit more manageable optimization problems compared to (3.4). To derive the bounds, we invoke the Poincaré inequality [42] and a Cramér-Rao-like inequality [44] for the input random variable X.

**DEFINITION 3.4 (Poincaré and Cramér-Rao inequality).** Given a random vector X taking values in $\mathbb{R}^d$, we denote the smallest and the largest constants by $C(X) \ge 0$ and $c(X) \ge 0$, respectively, such that

$$ (3.5) \qquad c(X) \|\mathbb{E}[\nabla f(X)]\|^2 \le \mathbb{E}[(f(X) - \mathbb{E}[f(X)])^2] \le C(X) \mathbb{E}[\|\nabla f(X)\|^2], $$

holds for any continuously differentiable function $f: \mathbb{R}^d \to \mathbb{R}$. If $C(X) < \infty$, the right inequality is called the Poincaré inequality with constant $C(X)$ and, if $c(X) > 0$, the left inequality is called the Cramér-Rao inequality with constant $c(X)$.

It is well known that a Gaussian random vector $X \sim N(\mu, \Sigma)$ on $\mathbb{R}^d$ satisfies (3.5) with $c(X) = \lambda_{\min}(\Sigma)$ and $C(X) = \lambda_{\max}(\Sigma)$, see Lemma D.1 and [6]. More details regarding the Poincaré inequality can be found in [10, 54] and references therein. For the Cramér-Rao inequality see the discussion in Appendix D. We defer to Proposition 3.6 for sufficient conditions so that $c(X) > 0$ and $C(X) < \infty$ hold. We now state our main result. The proof can be found in Appendix C.

**THEOREM 3.5.** Given a random vector X taking values in $\mathbb{R}^d$, we let $\bar{C}(X) \ge 0$ and $\bar{c}(X) \ge 0$ be the smallest and the largest constant, respectively, such that

$$ (3.6) \qquad C(X_{\perp} | X_r = x_r) \le \bar{C}(X), $$

$$ (3.7) \qquad c(X_{\perp} | X_r = x_r) \ge \bar{c}(X), $$

holds for all $x_r \in \mathbb{R}^r$ and all $U_r \in \mathbb{R}^{d \times r}$ with orthonormal columns, where $X_r = U_r^T X$ and $X_{\perp} = U_{\perp}^T X$ with $U_{\perp} \in \mathbb{R}^{d \times (d-r)}$ being an orthogonal complement to $U_r$. In other words, we require every conditional distribution of $X_{\perp}$ to satisfy the Poincaré and Cramér-Rao inequality.

Let $G: \mathbb{R}^d \to \mathbb{R}^m$ be a continuously differentiable function such that $E[\|\nabla G(X)\|_F^2] < \infty$, where $\nabla G(X) \in \mathbb{R}^{m \times d}$ denotes the Jacobian of G and $\|\cdot\|_F$ is the Frobenius norm. Then, for any $U_r \in \mathbb{R}^{d \times r}$ and $V_s \in \mathbb{R}^{m \times s}$ with orthonormal columns, we have

$$ (3.8) \qquad \mathbb{E}[\|G(X) - G^*(X)\|^2] \le \bar{C}(X) (\mathbb{E}[\|\nabla G(X)\|_F^2] - \mathbb{E}[\|V_s^\top \nabla G(X) U_r\|_F^2]), $$

$$ (3.9) \qquad \mathbb{E}[\|G(X) - G^*(X)\|^2] \ge \bar{c}(X) (\|\mathbb{E}[\nabla G(X)]\|_F^2 - \|V_s^\top \mathbb{E}[\nabla G(X)] U_r\|_F^2), $$

with $G^*$ defined as in Proposition 3.3.
Instead of solving the computationally expensive optimization problem (3.4), we
can now optimize the upper bound (3.8) of the L²-error to find the subspaces, meaning

$$
(3.10) \qquad \max \mathbb{E} \left[ \| V_s^\top \nabla G(X) U_r \|_F^2 \right].
$$

We will show in the next subsection that maximizing over either $U_r$ or $V_s$ yields closed-form solutions that can be computed as dominant eigenvectors of certain diagnostic matrices. There is, however, no closed-form solution for the joint maximization over $U_r$ and $V_s$. Nonetheless, since the Stiefel manifold $\text{St}(p, q) = \{A \in \mathbb{R}^{p \times q} \mid A^\top A = I_q\}$ is compact [1] and the objective function is continuous with respect to $U_r$ and $V_s$, the joint maximization problem has a solution. We propose a natural heuristic for its computation in the next subsection. After computing $U_r$ and/or $V_s$, the lower-bound (3.9) can be evaluated with no additional cost and allows us to certify the error made by optimizing the upper bound.

The following proposition gives sufficient condition on the density $\pi_X$ which en-
sures $\bar{c}(X) > 0$ and $\bar{C}(X) < \infty$.

PROPOSITION 3.6. The following holds:

1. Assume that $\text{supp}(\pi_X)$ is convex and that there exists $\rho > 0$ so that

$$
(3.11) \qquad - \nabla^2 \log \pi_X(x) \geq \rho I_d
$$

holds for any $x \in \mathbb{R}^d$. Then $\bar{C}(X) \le 1/\rho$.

2. Assume that $\text{supp}(\pi_X) = \mathbb{R}^d$ and that there exists $\eta > 0$ so that

$$
(3.12) \qquad - \nabla^2 \log \pi_X(x) \leq \eta I_d
$$

holds for any $x \in \mathbb{R}^d$. Then $\bar{c}(X) \ge 1/\eta$.

Proof. For the proof of 1 see [54], and for 2 see Lemma D.2.

<figure><img src="image_1.png" /></figure>

While assumption (3.11) is a classical strict log-concavity assumption on $\pi_X$, the assumption (3.12) means that the density $\pi_X$ has to lie above some scaled Gaussian density (see Remark D.3).

Remark 3.7 (Preconditioning). Using Proposition 3.6, one can show for $X \sim \mathcal{N}(\mu, \Sigma)$ that $\bar{C}(X) = \lambda_{\max}(\Sigma)$ and $\bar{c}(X) = \lambda_{\min}(\Sigma)$. Using a preconditioned vector $\bar{X} = \Sigma^{-1/2}(X - \mu) \sim \mathcal{N}(0, I_d)$, we get sharper bounds for (3.8) and (3.9) with $\bar{C}(\bar{X}) = \bar{c}(\bar{X}) = 1$. Applying Theorem 3.5 on $\bar{X} \mapsto G(\Sigma^{1/2}\bar{X} + \mu)$ yields

$$
\mathbb{E} \left[ \| G(X) - G^*(X) \|^2 \right] \le \mathbb{E} \left[ \| \nabla G(X) \Sigma^{1/2} \|^2_F \right] - \mathbb{E} \left[ \| \bar{V}_s^\top \nabla G(X) \Sigma^{1/2} \bar{U}_r \|^2_F \right],
$$

where we use $G^*(X) = V_s^\top \mathbb{E}[G(X)|X_r] + (I_m - V_s V_s^\top) \mathbb{E}[G(X)]$ with reduced pa-
rameter $X_r = (\bar{U}_r^\top \Sigma^{-1/2})(X - \mu)$. In practice, even for non-Gaussian $X$, it can be
beneficial to apply the change of variable $\bar{X} = \text{Cov}(X)^{-1/2}(X - \mathbb{E}[X])$ and compute
the dimension-reduction with the isotropic random variable $\bar{X}$ instead of $X$.

Remark 3.8 (Equality and Closed-form Solution for Affine Models). For affine models of the form $G(x) = a + Mx$ with $a \in \mathbb{R}^m, M \in \mathbb{R}^{m \times d}$, the Jacobian $\nabla G(x) = M$ is constant. If we further assume $x \sim \mathcal{N}(0, I_d)$, so that $\bar{C}(x) = \bar{c}(x) = 1$, the upper and lower bound coincide with the L²-error

$$
\mathbb{E} \left[ \| G(X) - G^*(X) \|^2 \right] = \| M \|_F^2 - \| V_s^\top M U_r \|_F^2.
$$
Instead of solving the computationally expensive optimization problem (3.4), we
can now optimize the upper bound (3.8) of the L²-error to find the subspaces, meaning

$$
(3.10) \qquad \max \mathbb{E} \left[ \| V_s^\top \nabla G(X) U_r \|_F^2 \right].
$$

We will show in the next subsection that maximizing over either $U_r$ or $V_s$ yields closed-form solutions that can be computed as dominant eigenvectors of certain diagnostic matrices. There is, however, no closed-form solution for the joint maximization over $U_r$ and $V_s$. Nonetheless, since the Stiefel manifold $\text{St}(p, q) = \{A \in \mathbb{R}^{p \times q} \mid A^\top A = I_q\}$ is compact [1] and the objective function is continuous with respect to $U_r$ and $V_s$, the joint maximization problem has a solution. We propose a natural heuristic for its computation in the next subsection. After computing $U_r$ and/or $V_s$, the lower-bound (3.9) can be evaluated with no additional cost and allows us to certify the error made by optimizing the upper bound.

The following proposition gives sufficient condition on the density $\pi_X$ which en-
sures $\bar{c}(X) > 0$ and $\bar{C}(X) < \infty$.

PROPOSITION 3.6. The following holds:

1. Assume that $\text{supp}(\pi_X)$ is convex and that there exists $\rho > 0$ so that

$$
(3.11) \qquad - \nabla^2 \log \pi_X(x) \geq \rho I_d
$$

holds for any $x \in \mathbb{R}^d$. Then $\bar{C}(X) \le 1/\rho$.

2. Assume that $\text{supp}(\pi_X) = \mathbb{R}^d$ and that there exists $\eta > 0$ so that

$$
(3.12) \qquad - \nabla^2 \log \pi_X(x) \leq \eta I_d
$$

holds for any $x \in \mathbb{R}^d$. Then $\bar{c}(X) \ge 1/\eta$.

Proof. For the proof of 1 see [54], and for 2 see Lemma D.2.

<figure><img src="image_2.png" /></figure>

While assumption (3.11) is a classical strict log-concavity assumption on $\pi_X$, the assumption (3.12) means that the density $\pi_X$ has to lie above some scaled Gaussian density (see Remark D.3).

Remark 3.7 (Preconditioning). Using Proposition 3.6, one can show for $X \sim \mathcal{N}(\mu, \Sigma)$ that $\bar{C}(X) = \lambda_{\max}(\Sigma)$ and $\bar{c}(X) = \lambda_{\min}(\Sigma)$. Using a preconditioned vector $\bar{X} = \Sigma^{-1/2}(X - \mu) \sim \mathcal{N}(0, I_d)$, we get sharper bounds for (3.8) and (3.9) with $\bar{C}(\bar{X}) = \bar{c}(\bar{X}) = 1$. Applying Theorem 3.5 on $\bar{X} \mapsto G(\Sigma^{1/2}\bar{X} + \mu)$ yields

$$
\mathbb{E} \left[ \| G(X) - G^*(X) \|^2 \right] \le \mathbb{E} \left[ \| \nabla G(X) \Sigma^{1/2} \|^2_F \right] - \mathbb{E} \left[ \| \bar{V}_s^\top \nabla G(X) \Sigma^{1/2} \bar{U}_r \|^2_F \right],
$$

where we use $G^*(X) = V_s^\top \mathbb{E}[G(X)|X_r] + (I_m - V_s V_s^\top) \mathbb{E}[G(X)]$ with reduced pa-
rameter $X_r = (\bar{U}_r^\top \Sigma^{-1/2})(X - \mu)$. In practice, even for non-Gaussian $X$, it can be
beneficial to apply the change of variable $\bar{X} = \text{Cov}(X)^{-1/2}(X - \mathbb{E}[X])$ and compute
the dimension-reduction with the isotropic random variable $\bar{X}$ instead of $X$.

Remark 3.8 (Equality and Closed-form Solution for Affine Models). For affine models of the form $G(x) = a + Mx$ with $a \in \mathbb{R}^m, M \in \mathbb{R}^{m \times d}$, the Jacobian $\nabla G(x) = M$ is constant. If we further assume $x \sim \mathcal{N}(0, I_d)$, so that $\bar{C}(x) = \bar{c}(x) = 1$, the upper and lower bound coincide with the L²-error

$$
\mathbb{E} \left[ \| G(X) - G^*(X) \|^2 \right] = \| M \|_F^2 - \| V_s^\top M U_r \|_F^2.
$$
Thus, the solution to (3.10) is actually the same as the one of (3.4). By the Eckart-Young-Mirsky theorem, the optimal $U_r$ and $V_s$ minimizing the right hand side have columns consisting of the dominant right and left eigenvectors of $M$, respectively, meaning that $U_r^* = [u_1 \dots u_r]$ and $V_s^* = [v_1 \dots v_s]$ where

$$M = \sum_{i=1}^{\min(m,d)} \lambda_i v_i u_i^\top$$

is the singular value decomposition of $M$.

**3.2. Minimizing the Error Bound.** We define the diagnostic matrices

$$(3.13) \qquad H_X(V_s) = \mathbb{E}[\nabla G(X)^\top V_s V_s^\top \nabla G(X)] \in \mathbb{R}^{d \times d},$$

$$(3.14) \qquad H_Y(U_r) = \mathbb{E}[\nabla G(X) U_r U_r^\top \nabla G(X)^\top] \in \mathbb{R}^{m \times m}.$$

Then the objective function in (3.10) can be written as

$$(3.15) \qquad \mathbb{E}[\|V_s^\top \nabla G(X) U_r\|_F^2] = \mathrm{Tr}(V_s^\top H_Y(U_r) V_s) = \mathrm{Tr}(U_r^\top H_X(V_s) U_r).$$

To find the optimal subspaces that minimize the objective above, we recall the following property on Hermitian matrices (see e.g. Corollary 4.3.39 [27]).

**LEMMA 3.9 (Variational characterization of eigenvalues of Hermitian matrices).**
Let $H \in \mathbb{R}^{d \times d}$ be a symmetric positive definite matrix with eigenpairs $(\lambda_i, w_i)$ meaning $Hw_i = \lambda_i w_i$, where $\lambda_{i+1} \le \lambda_i$ and $\|w_i\|_2 = 1$ for $i = 1, \dots, d$. Then for any $r < d$ we have

$$\max_{\substack{W_r \in \mathbb{R}^{d \times r} \\ W_r^\top W_r = I_r}} \mathrm{Tr}(W_r^\top H W_r) = \sum_{i=1}^r \lambda_i,$$

where a solution is given by $W_r = [w_1, \dots, w_r]$.

Applying Lemma 3.9 to the second term $\mathrm{Tr}(V_s^\top H_Y(U_r) V_s)$ in (3.15), we obtain the solution $V_s^*(U_r)$ to the goal-oriented problem (3.10) with fixed $U_r$. The solution requires computing the eigendecomposition $H_Y(U_r) = \sum_{i=1}^m \lambda_i^Y(U_r) v_i v_i^\top$ and assembling

$$(3.16) \qquad V_s^*(U_r) = [v_1, \dots, v_s].$$

In a similar fashion, the third term in (3.15) suggests that the solution $U_r^*(V_s)$ to the goal-oriented problem (3.10) with fixed $V_s$ can be obtained by computing the eigendecomposition $H_X(V_s) = \sum_{i=1}^d \lambda_i^X(V_s) u_i u_i^\top$ and assembling

$$(3.17) \qquad U_r^*(V_s) = [u_1, \dots, u_r].$$

By (3.8), we also get the following error bounds for the two goal-oriented problems

$$\mathbb{E}[\|G(X) - G^*(X)\|^2] \le \bar{C}(X) \left( \mathbb{E}[\|\nabla G(X)\|_F^2] - \sum_{i=1}^s \lambda_i^Y(U_r) \right),$$

$$\mathbb{E}[\|G(X) - G^*(X)\|^2] \le \bar{C}(X) \left( \mathbb{E}[\|\nabla G(X)\|_F^2] - \sum_{i=1}^r \lambda_i^X(V_s) \right).$$
Thus, the solution to (3.10) is actually the same as the one of (3.4). By the Eckart-Young-Mirsky theorem, the optimal $U_r$ and $V_s$ minimizing the right hand side have columns consisting of the dominant right and left eigenvectors of $M$, respectively, meaning that $U_r^* = [u_1 \dots u_r]$ and $V_s^* = [v_1 \dots v_s]$ where

$$M = \sum_{i=1}^{\min(m,d)} \lambda_i v_i u_i^{\top}$$

is the singular value decomposition of $M$.

### 3.2. Minimizing the Error Bound. We define the diagnostic matrices

$$(3.13) \qquad H_X(V_s) = \mathbb{E}[\nabla G(X)^\top V_s V_s^\top \nabla G(X)] \in \mathbb{R}^{d \times d},$$

$$(3.14) \qquad H_Y(U_r) = \mathbb{E}[\nabla G(X) U_r U_r^\top \nabla G(X)^\top] \in \mathbb{R}^{m \times m}.$$

Then the objective function in (3.10) can be written as

$$(3.15) \qquad \mathbb{E}[\|V_s^\top \nabla G(X) U_r\|_F^2] = \mathrm{Tr}(V_s^\top H_Y(U_r) V_s) = \mathrm{Tr}(U_r^\top H_X(V_s) U_r).$$

To find the optimal subspaces that minimize the objective above, we recall the following property on Hermitian matrices (see e.g. Corollary 4.3.39 [27]).

**LEMMA 3.9 (Variational characterization of eigenvalues of Hermitian matrices).**
Let $H \in \mathbb{R}^{d \times d}$ be a symmetric positive definite matrix with eigenpairs $(\lambda_i, w_i)$ meaning $Hw_i = \lambda_i w_i$, where $\lambda_{i+1} \le \lambda_i$ and $\|w_i\|_2 = 1$ for $i = 1, \dots, d$. Then for any $r < d$ we have

$$\max_{\substack{W_r \in \mathbb{R}^{d \times r} \\ W_r^\top W_r = I_r}} \mathrm{Tr}(W_r^\top H W_r) = \sum_{i=1}^r \lambda_i,$$

where a solution is given by $W_r = [w_1, \dots, w_r]$.

Applying Lemma 3.9 to the second term $\mathrm{Tr}(V_s^\top H_Y(U_r) V_s)$ in (3.15), we obtain the solution $V_s^*(U_r)$ to the goal-oriented problem (3.10) with fixed $U_r$. The solution requires computing the eigendecomposition $H_Y(U_r) = \sum_{i=1}^m \lambda_i^Y(U_r) v_i v_i^\top$ and assembling

$$(3.16) \qquad V_s^*(U_r) = [v_1, \dots, v_s].$$

In a similar fashion, the third term in (3.15) suggests that the solution $U_r^*(V_s)$ to the goal-oriented problem (3.10) with fixed $V_s$ can be obtained by computing the eigendecomposition $H_X(V_s) = \sum_{i=1}^d \lambda_i^X(V_s) u_i u_i^\top$ and assembling

$$(3.17) \qquad U_r^*(V_s) = [u_1, \dots, u_r].$$

By (3.8), we also get the following error bounds for the two goal-oriented problems

$$\mathbb{E}[\|G(X) - G^*(X)\|^2] \le \bar{C}(X) \left( \mathbb{E}[\|\nabla G(X)\|_F^2] - \sum_{i=1}^s \lambda_i^Y(U_r) \right),$$

$$\mathbb{E}[\|G(X) - G^*(X)\|^2] \le \bar{C}(X) \left( \mathbb{E}[\|\nabla G(X)\|_F^2] - \sum_{i=1}^r \lambda_i^X(V_s) \right).$$
For the joint minimization over $(U_r, V_s)$ for the coupled problem (3.10), there is no closed form solution in general. To approximate the solution, we propose the iterative alternating eigendecomposition scheme

$$ (3.18) \qquad \begin{aligned} U_r^{k+1} &= U_r^*(V_s^k), \\ V_s^{k+1} &= V_s^*(U_r^{k+1}), \end{aligned} $$

where $U_r^*(\cdot)$ and $V_s^*(\cdot)$ are defined as in (3.17) and (3.16). If $\{(U_r^k, V_s^k)\}_{k \ge 1}$ converges towards a limit point $(U_r^*, V_s^*)$ then, by construction, this point will be a stationary point of the objective function. In practice, the diagnostic matrices can be computed using Monte Carlo estimators and the operational algorithm is formalized in Algorithm 3.1.

**Algorithm 3.1 Alternating eigendecomposition for coupled input and output dimension reduction.**

**Input:** $s, r \in \mathbb{N}$; $V_s^0 \in \mathbb{R}^{m \times s}$; $G_i = \nabla G(x^{(i)})$ where $x^{(i)} \stackrel{i.i.d.}{\sim} \pi_X$

1: Set $k = 0$

2: **while** $k \le$ maximal iteration **do**

3: Assemble $\hat{H}_X(V_s^k) = \frac{1}{M} \sum_{i=1}^M G_i^\top V_s^k V_s^{k\top} G_i$

4: Compute SVD of $\hat{H}_X(V_s^k) = U\Lambda U^\top$

5: Set $U_r^k = [u_1, \dots, u_r]$

6: Assemble $\hat{H}_Y(U_r^k) = \frac{1}{M} \sum_{i=1}^M G_i U_r^k U_r^{k\top} G_i^\top$

7: Compute SVD of $\hat{H}_Y(U_r^k) = V\Lambda V^\top$

8: Set $V_s^{k+1} = [v_1, \dots, v_s]$

9: $k = k + 1$

10: **end while**

11: **return** $U_r, V_s$

**Remark 3.10 (Implicit diagnostic matrices).** Algorithm 3.1 requires computing eigencompositions of the matrices $\hat{H}_X(V_s)$ and $\hat{H}_Y(U_r)$ which, when $d \gg 1$ or $m \gg 1$, can be numerically costly to assemble and to store. To remedy this, one can compute the eigencompositions using iterative algorithms (e.g., Krylov methods) which only require access to matrix-vector product operations of the form

$$ v \mapsto \hat{H}_Y(U_r)v = \frac{1}{M} \sum_{i=1}^{M} \nabla G(X^{(i)}) (U_r U_r^\top (\nabla G(X^{(i)})^\top v)), $$

$$ u \mapsto \hat{H}_X(V_s)u = \frac{1}{M} \sum_{i=1}^{M} \nabla G(X^{(i)})^\top (V_s V_s^\top (\nabla G(X^{(i)})u)). $$

Computing these matrix-vector products requires only evaluating the linear-tangent operator $u \mapsto \nabla G(X^{(i)})u$ and its adjoint $v \mapsto \nabla G(X^{(i)})^\top v$, which are readily available for many applications or can be implemented efficiently, e.g., using finite difference and automatic differentiation.

**Remark 3.11.** In a recent paper [12], the authors proposed a joint dimension reduction method for the parameter and data space based on similar diagnostic matrices. Using an information theoretic approach, they derive an upper bound for the
For the joint minimization over $(U_r, V_s)$ for the coupled problem (3.10), there is no closed form solution in general. To approximate the solution, we propose the iterative alternating eigendecomposition scheme

$$ (3.18) \qquad \begin{aligned} U_r^{k+1} &= U_r^*(V_s^k), \\ V_s^{k+1} &= V_s^*(U_r^{k+1}), \end{aligned} $$

where $U_r^*(\cdot)$ and $V_s^*(\cdot)$ are defined as in (3.17) and (3.16). If $\{(U_r^k, V_s^k)\}_{k \ge 1}$ converges towards a limit point $(U_r^*, V_s^*)$ then, by construction, this point will be a stationary point of the objective function. In practice, the diagnostic matrices can be computed using Monte Carlo estimators and the operational algorithm is formalized in Algorithm 3.1.

**Algorithm 3.1 Alternating eigendecomposition for coupled input and output dimension reduction.**

Input: $s, r \in \mathbb{N}$; $V_s^0 \in \mathbb{R}^{m \times s}$; $G_i = \nabla G(x^{(i)})$ where $x^{(i)} \stackrel{i.i.d.}{\sim} \pi_X$

1: Set $k = 0$

2: **while** $k \le$ maximal iteration **do**

3: Assemble $\hat{H}_X(V_s^k) = \frac{1}{M} \sum_{i=1}^M G_i^\top V_s^k V_s^{k\top} G_i$

4: Compute SVD of $\hat{H}_X(V_s^k) = U\Lambda U^\top$

5: Set $U_r^k = [u_1, \dots, u_r]$

6: Assemble $\hat{H}_Y(U_r^k) = \frac{1}{M} \sum_{i=1}^M G_i U_r^k U_r^{k\top} G_i^\top$

7: Compute SVD of $\hat{H}_Y(U_r^k) = V\Lambda V^\top$

8: Set $V_s^{k+1} = [v_1, \dots, v_s]$

9: $k = k + 1$

10: **end while**

11: **return** $U_r, V_s$

Remark 3.10 (Implicit diagnostic matrices). Algorithm 3.1 requires computing eigencompositions of the matrices $\hat{H}_X(V_s)$ and $\hat{H}_Y(U_r)$ which, when $d \gg 1$ or $m \gg 1$, can be numerically costly to assemble and to store. To remedy this, one can compute the eigencompositions using iterative algorithms (e.g., Krylov methods) which only require access to matrix-vector product operations of the form

$$ v \mapsto \hat{H}_Y(U_r)v = \frac{1}{M} \sum_{i=1}^{M} \nabla G(X^{(i)}) (U_r U_r^\top (\nabla G(X^{(i)})^\top v)), $$

$$ u \mapsto \hat{H}_X(V_s)u = \frac{1}{M} \sum_{i=1}^{M} \nabla G(X^{(i)})^\top (V_s V_s^\top (\nabla G(X^{(i)}) u)). $$

Computing these matrix-vector products requires only evaluating the linear-tangent operator $u \mapsto \nabla G(X^{(i)})u$ and its adjoint $v \mapsto \nabla G(X^{(i)})^\top v$, which are readily available for many applications or can be implemented efficiently, e.g., using finite difference and automatic differentiation.

Remark 3.11. In a recent paper [12], the authors proposed a joint dimension reduction method for the parameter and data space based on similar diagnostic matrices. Using an information theoretic approach, they derive an upper bound for the
posterior error $\mathbb{E}_Y[\mathbb{D}_{KL}(\pi_{X|Y} \Vert \tilde{\pi}_{X|Y})]$ which is minimised by $U_r^*$ and $V_s^*$ consisting of the dominant eigenvectors of

$$H_X = \mathbb{E}[\nabla G(X)^\top \nabla G(X)],$$

$$H_Y = \mathbb{E}[\nabla G(X) \nabla G(X)^\top],$$

respectively. Their results correspond to our goal-oriented solutions with $V_s^* = V_s^*(I_d)$ and $U_r^* = U_r^*(I_m)$. Contrary to our error bound, their results depend on the logarithmic Sobolev constant of the joint distribution $\pi_{XY}$ whose existence is generally difficult to establish.

# 4. Application to Goal-oriented Coordinate Selection.

## 4.1. Goal-Oriented Bayesian Optimal Experimental Design (BOED).

Consider the problem of finding an optimal subset of $s \le m$ sensors out of $m$ candidate locations that maximizes the information gained about the parameter $X$ that is to be inferred. The vector $Y = G(X) + \eta$ gathers all possible measurements at the $m$ candidate locations. Data collected at select sensor locations $\tau \subseteq \{1, \dots, m\}$, $|\tau| = s$ are denoted by $Y_\tau = V_\tau^\top Y$, where $V_\tau = [e_\tau^m] \in \mathbb{R}^{m \times s}$ contains the canonical basis vectors $e_\tau^m \in \mathbb{R}^m$ indexed by $\tau$. Later, we will propose an alternative approach for experimental design where the reduced data is characterized by a matrix $V_s \in \mathbb{R}^{m \times s}$ with orthonormal columns. We will explain the benefits of considering this relaxed design space later and how to associate sensor locations to the design matrix $V_s$. For any $V_\star$, which can be a coordinate selection design $V_\star = V_\tau$ or a relaxed orthogonal design $V_\star = V_s$, we have the posterior density of $X$ knowing $Y_\star = V_\star^\top Y$ given as

$$(4.1) \qquad \pi_{X|V_\star^\top Y}(x|y_\star) \propto \exp\left(-\frac{1}{2\sigma^2} \|y_\star - V_\star^\top G(x)\|^2\right) \pi_X(x).$$

To quantify the information content of an experimental design $V_\star$, we use the expected information gain (EIG) defined as the KL-divergence between posterior and prior distribution in expectation over the data

$$(4.2) \qquad \Phi(V_\star) = \mathbb{E}_Y[\mathbb{D}_{KL}(\pi_{X|V_\star^\top Y}(\cdot | V_\star^\top Y) \Vert \pi_X)].$$

The divergence reflects the extent to which the posterior distribution differs from the prior due to observed data and, thereby, represents the information gain of the design $V_\star$. The EIG is equivalent to the mutual information between $X$ and $Y$ given as $I(X, Y) = \mathbb{D}_{KL}(\pi_{X,Y} \Vert \pi_X \pi_Y)$, i.e., the KL-divergence from the product of marginal densities to their joint density. This perspective can be useful for models with inaccessible likelihood functions [33]. For a linear Gaussian inference problem, optimizing the EIG corresponds to minimizing the determinant of the posterior covariance matrix, also referred to as the D-optimality criterion. It can be interpreted as minimizing the volume of the uncertainty ellipsoid around the maximum a-posteriori estimator of the parameter $X$ [2].

Estimating the EIG poses an even greater challenge than solving the Bayesian inference problem as it requires solving the inverse problem itself as a sub-problem. The EIG consists of an expectation over both the posterior and data distributions, necessitating an expensive nested Monte Carlo estimator [47, 14]. Efficient EIG estimation is an active research area, and various approaches exist. Some rely on constructing tractable approximations to the posterior distribution using measure transport [32, 11] or variational inference [22]. Others employ different types of surrogate models [30, 52, 7]. A widely used strategy is to resort to Gaussian linear models for which
posterior error $E_Y[D_{KL}(\pi_{X|Y} || \tilde{\pi}_{X|Y})]$ which is minimised by $U_r^*$ and $V_s^*$ consisting of the dominant eigenvectors of

$$H_X = E[\nabla G(X)^{\top} \nabla G(X)],$$

$$H_Y = E[\nabla G(X) \nabla G(X)^{\top}],$$

respectively. Their results correspond to our goal-oriented solutions with $V_s^* = V_s^*(I_d)$ and $U_r^* = U_r^*(I_m)$. Contrary to our error bound, their results depend on the logarithmic Sobolev constant of the joint distribution $\pi_{XY}$ whose existence is generally difficult to establish.

# 4. Application to Goal-oriented Coordinate Selection.

## 4.1. Goal-Oriented Bayesian Optimal Experimental Design (BOED).

Consider the problem of finding an optimal subset of $s \le m$ sensors out of $m$ candidate locations that maximizes the information gained about the parameter $X$ that is to be inferred. The vector $Y = G(X) + \eta$ gathers all possible measurements at the $m$ candidate locations. Data collected at select sensor locations $\tau \subseteq \{1, \dots, m\}$, $|\tau| = s$ are denoted by $Y_\tau = V_\tau^\top Y$, where $V_\tau = [e_\tau^m] \in \mathbb{R}^{m \times s}$ contains the canonical basis vectors $e_\tau^m \in \mathbb{R}^m$ indexed by $\tau$. Later, we will propose an alternative approach for experimental design where the reduced data is characterized by a matrix $V_s \in \mathbb{R}^{m \times s}$ with orthonormal columns. We will explain the benefits of considering this relaxed design space later and how to associate sensor locations to the design matrix $V_s$. For any $V_\star$, which can be a coordinate selection design $V_\star = V_\tau$ or a relaxed orthogonal design $V_\star = V_s$, we have the posterior density of $X$ knowing $Y_\star = V_\star^\top Y$ given as

$$(4.1) \qquad \pi_{X|V_\star^\top Y}(x|y_\star) \propto \exp\left(-\frac{1}{2\sigma^2} \|y_\star - V_\star^\top G(x)\|^2\right) \pi_X(x).$$

To quantify the information content of an experimental design $V_\star$, we use the expected information gain (EIG) defined as the KL-divergence between posterior and prior distribution in expectation over the data

$$(4.2) \qquad \Phi(V_\star) = E_Y[D_{KL}(\pi_{X|V_\star^\top Y}(\cdot | V_\star^\top Y) || \pi_X)].$$

The divergence reflects the extent to which the posterior distribution differs from the prior due to observed data and, thereby, represents the information gain of the design $V_\star$. The EIG is equivalent to the mutual information between $X$ and $Y$ given as $I(X, Y) = D_{KL}(\pi_{X,Y} || \pi_X \pi_Y)$, i.e., the KL-divergence from the product of marginal densities to their joint density. This perspective can be useful for models with inaccessible likelihood functions [33]. For a linear Gaussian inference problem, optimizing the EIG corresponds to minimizing the determinant of the posterior covariance matrix, also referred to as the D-optimality criterion. It can be interpreted as minimizing the volume of the uncertainty ellipsoid around the maximum a-posteriori estimator of the parameter $X$ [2].

Estimating the EIG poses an even greater challenge than solving the Bayesian inference problem as it requires solving the inverse problem itself as a sub-problem. The EIG consists of an expectation over both the posterior and data distributions, necessitating an expensive nested Monte Carlo estimator [47, 14]. Efficient EIG estimation is an active research area, and various approaches exist. Some rely on constructing tractable approximations to the posterior distribution using measure transport [32, 11] or variational inference [22]. Others employ different types of surrogate models [30, 52, 7]. A widely used strategy is to resort to Gaussian linear models for which
a closed-form expression of the EIG exists [3, 51, 9]. In addition to the challenging EIG estimation, identifying the optimal set of sensor placements requires solving a combinatorial optimization problem with $\binom{m}{s}$ possible combinations. This NP-hard problem is intractable even for moderate dimensions. Hence, sub-optimal solutions are commonly identified using greedy methods [34, 16, 51] or by solving continuous relaxations of the problem with sparsity constraints [26, 9].

We consider, in particular, goal-oriented BOEDs, where the final target is to infer the parameter in a lower-dimensional subspace of interest. For example, if the PDE model describes a spatio-temporal field with $X$ as the initial state, the initial state in a localized region of the domain could be a parameter of interest. Tailoring the experimental design to maximize the information gained with respect to the parameter of interest allows for more efficient and accurate inference from experimental data. For a given parameter of interest $X_r = U_r^\top X$, we seek the optimal design $V_\star$ that maximizes the goal-oriented EIG (see (O3))

$$ (4.3) \qquad \max_{V_\star} \Phi(V_\star|U_r) = \mathbb{E}_Y [\mathbb{D}_{\text{KL}} (\pi_{X_r|V_\star^\top Y}(\cdot|V_\star^\top Y)\|\pi_{X_r})], $$

where $\pi_{X_r|V_\star^\top Y}$ and $\pi_{X_r}$ are the marginalized densities. The goal-oriented BOED problem faces the same obstacles as the non-goal-oriented one, namely the expensive EIG estimation and combinatorial optimization problem. We avoid both challenges by maximizing a lower bound for the EIG. Although our method requires the assumption of a Gaussian likelihood, it allows for non-linear forward models.

Combining Lemma 3.2 with Theorem 3.5 permits us to optimize a lower bound of the goal-oriented EIG by maximizing $\text{Tr}(V_\star^\top H_Y(U_r)V_\star)$. In the coordinate selection case with $V_\star = V_\tau$, this reduces to solving the problem

$$ (4.4) \qquad \max_{\substack{V_\tau=[e_\tau]\in\mathbb{R}^{m\times s} \\ \tau\subset\{1,\dots,m\}, |\tau|=s}} \text{Tr}(V_\tau^\top H_Y(U_r)V_\tau) = \sum_{i\in\tau^*} H_Y(U_r)_{ii}, $$

where the solution $\tau^* \subseteq \{1, \dots, m\}$ is the index set of the $s$ largest diagonal entries of $H_Y(U_r)$. Not only do we bypass computing the EIG, but our solution $\tau^*$ can be computed without combinatorial optimization.

For orthogonal designs $V_\star = V_s$, we obtain with Lemma 3.9 the solution

$$ (4.5) \qquad \max_{\substack{V_s \in \mathbb{R}^{m \times s} \\ V_s^\top V_s = I_m}} \text{Tr}(V_s^\top H_Y(U_r) V_s) = \sum_{i=1}^s \lambda_i^Y(U_r), $$

with $V_s^*(U_r) = [v_1, \dots, v_s]$ and $(\lambda_i^Y, v_i)$ being the $i$-th largest eigenpair of $H_Y(U_r) = \sum_{i=1}^m \lambda_i^Y(U_r)v_i v_i^\top$ with $\lambda_{i+1}^Y(U_r) \le \lambda_i^Y(U_r)$. Given that (4.5) is a relaxation of (4.4), we expect $V_s^*(U_r)$ to yield better results compared to $V_{\tau^*}$, meaning $\sum_{i=1}^s \lambda_i^Y(U_r) \ge \sum_{i \in \tau^*} H_Y(U_r)_{ii}$. To associate sensor locations with the optimal subspace $V_s^*(U_r)$, we propose to employ the empirical interpolation method (EIM) [13, 40].

The EIM takes as input a set of basis vectors $v_1, \dots, v_s \in \mathbb{R}^m$ and returns a set of indices $\{i_1, \dots, i_s\} \subseteq \{1, \dots, m\}$ corresponding to the empirically-best interpolation points for the given basis. Combining the EIM with dimension reduction methods for sensor placement has been considered, for example, in [25, 8]. By using the EIM on the columns of $V_s^*(U_r)$, we can thus associate sensor locations $\tilde{\tau} = \{i_1, \dots, i_s\}$ to the orthogonal design matrix $V_s^*(U_r)$. Compared to $\tau^*$, the EIM-based sensors $\tilde{\tau}$ are less optimal with respect to the goal-oriented EIG lower bound (4.4). However, our
a closed-form expression of the EIG exists [3, 51, 9]. In addition to the challenging EIG estimation, identifying the optimal set of sensor placements requires solving a combinatorial optimization problem with $\binom{m}{s}$ possible combinations. This NP-hard problem is intractable even for moderate dimensions. Hence, sub-optimal solutions are commonly identified using greedy methods [34, 16, 51] or by solving continuous relaxations of the problem with sparsity constraints [26, 9].

We consider, in particular, goal-oriented BOEDs, where the final target is to infer the parameter in a lower-dimensional subspace of interest. For example, if the PDE model describes a spatio-temporal field with $X$ as the initial state, the initial state in a localized region of the domain could be a parameter of interest. Tailoring the experimental design to maximize the information gained with respect to the parameter of interest allows for more efficient and accurate inference from experimental data. For a given parameter of interest $X_r = U_r^\top X$, we seek the optimal design $V_\star$ that maximizes the goal-oriented EIG (see (O3))

$$ (4.3) \qquad \max_{V_\star} \Phi(V_\star|U_r) = \mathbb{E}_Y \left[ \mathbb{D}_{\text{KL}} \left( \pi_{X_r|V_\star^\top Y}(\cdot|V_\star^\top Y) \| \pi_{X_r} \right) \right], $$

where $\pi_{X_r|V_\star^\top Y}$ and $\pi_{X_r}$ are the marginalized densities. The goal-oriented BOED problem faces the same obstacles as the non-goal-oriented one, namely the expensive EIG estimation and combinatorial optimization problem. We avoid both challenges by maximizing a lower bound for the EIG. Although our method requires the assumption of a Gaussian likelihood, it allows for non-linear forward models.

Combining Lemma 3.2 with Theorem 3.5 permits us to optimize a lower bound of the goal-oriented EIG by maximizing $\text{Tr}(V_\star^\top H_Y(U_r)V_\star)$. In the coordinate selection case with $V_\star = V_\tau$, this reduces to solving the problem

$$ (4.4) \qquad \max_{\substack{V_\tau = [e_\tau] \in \mathbb{R}^{m \times s} \\ \tau \subset \{1, \dots, m\}, \ |\tau| = s}} \text{Tr} \left( V_\tau^\top H_Y(U_r) V_\tau \right) = \sum_{i \in \tau^*} H_Y(U_r)_{ii}, $$

where the solution $\tau^* \subseteq \{1, \dots, m\}$ is the index set of the $s$ largest diagonal entries of $H_Y(U_r)$. Not only do we bypass computing the EIG, but our solution $\tau^*$ can be computed without combinatorial optimization.

For orthogonal designs $V_\star = V_s$, we obtain with Lemma 3.9 the solution

$$ (4.5) \qquad \max_{\substack{V_s \in \mathbb{R}^{m \times s} \\ V_s^\top V_s = I_m}} \text{Tr} \left( V_s^\top H_Y(U_r) V_s \right) = \sum_{i=1}^s \lambda_i^Y(U_r), $$

with $V_s^*(U_r) = [v_1, \dots, v_s]$ and $(\lambda_i^Y, v_i)$ being the $i$-th largest eigenpair of $H_Y(U_r) = \sum_{i=1}^m \lambda_i^Y(U_r)v_i v_i^\top$ with $\lambda_{i+1}^Y(U_r) \le \lambda_i^Y(U_r)$. Given that (4.5) is a relaxation of (4.4), we expect $V_s^*(U_r)$ to yield better results compared to $V_{\tau^*}$, meaning $\sum_{i=1}^s \lambda_i^Y(U_r) \ge \sum_{i \in \tau^*} H_Y(U_r)_{ii}$. To associate sensor locations with the optimal subspace $V_s^*(U_r)$, we propose to employ the empirical interpolation method (EIM) [13, 40].

The EIM takes as input a set of basis vectors $v_1, \dots, v_s \in \mathbb{R}^m$ and returns a set of indices $\{i_1, \dots, i_s\} \subseteq \{1, \dots, m\}$ corresponding to the empirically-best interpolation points for the given basis. Combining the EIM with dimension reduction methods for sensor placement has been considered, for example, in [25, 8]. By using the EIM on the columns of $V_s^*(U_r)$, we can thus associate sensor locations $\tilde{\tau} = \{i_1, \dots, i_s\}$ to the orthogonal design matrix $V_s^*(U_r)$. Compared to $\tau^*$, the EIM-based sensors $\tilde{\tau}$ are less optimal with respect to the goal-oriented EIG lower bound (4.4). However, our
numerical experiments show that the EIM-based sensors often slightly outperform $\tau^*$ in terms of the true goal-oriented EIG rather than the lower bound. A reason could be that the sensors indexed by $\tilde{\tau}$ tend to be less localized than $\tau^*$. In the numerical experiments, we also apply the EIM on PCA bases to obtain sensor placements $\tilde{\tau}^{\text{PCA}}$, which we will refer to as PCA-EIM sensors.

**4.2. Goal-Oriented Sensitivity Analysis (GSA).** Global sensitivity analysis aims to assign an importance value to each set of input variables, reflecting its contribution to the output variations. It is often used for factor prioritization and/or fixing in complex models. Factor prioritization identifies the subset $\tau \subseteq \{1, \dots, d\}$, $|\tau| = r$ of input parameters causing the most output variation so that the parameters $X_{\tau} = (X_i)_{i \in \tau}$ can be prioritized e.g., in uncertainty reduction. Factor fixing, on the other hand, identifies the input set $\tau$ with the least influence on the output. Thus, $X_{\tau}$ can be fixed at a nominal value to reduce the complexity of, e.g., parameter estimation problems. For a scalar-valued square-integrable function $f: \mathbb{R}^d \to \mathbb{R}$, commonly used sensitivity indices are the closed Sobol' indices and the total Sobol' indices defined as

$$ (4.6) \qquad S_{\tau}^{\text{cl}} = \frac{\text{Var}(\mathbb{E}[f(X)|X_{\tau}])}{\text{Var}(f(X))} \quad \text{and} \quad S_{\tau}^{\text{tot}} = 1 - \frac{\text{Var}(\mathbb{E}[f(X)|X_{-\tau}])}{\text{Var}(f(X))}, $$

respectively, where $X_{-\tau} = (X_i)_{i \notin \tau}$ and the input variables $X_i$ for $i = 1, \dots, d$ are assumed to be independent. The closed Sobol' index quantifies the portion of output variance attributed to $X_{-\tau}$. A higher value implies a stronger influence of $X_{-\tau}$ on the output. The total Sobol' index encompasses output variation due to $X_{-\tau}$ along with variations caused by interactions of $X_{-\tau}$ with any of the remaining input variables. As a result, it holds that $S_{-\tau}^{\text{cl}} \le S_{-\tau}^{\text{tot}}$ and $1 = S_{-\tau}^{\text{tot}} + S_{-\tau}^{\text{cl}}$. A lower value in total Sobol' index suggests that the parameters indexed by $\tau$ have little direct nor indirect impact on the output.

We consider, in particular, goal-oriented GSA where some outputs of interest $V_s^\top G(X)$ are prescribed. We use definitions of Sobol' indices for vector-valued functions in [24] and generalize them to the goal-oriented case as follows

$$ (4.7) \qquad \begin{aligned} S^{\text{cl}}(U_{\tau}|V_s) &= \frac{\text{Tr}(\text{Cov}(\mathbb{E}[V_s^{\top}G(X)|U_{\tau}^{\top}X]))}{\text{Tr}(\text{Cov}(V_s^{\top}G(X)))}, \\ S^{\text{tot}}(U_{\tau}|V_s) &= 1 - \frac{\text{Tr}(\text{Cov}(V_s^{\top}G(X)|U_{-\tau}^{\top}X))}{\text{Tr}(\text{Cov}(V_s^{\top}G(X)))}. \end{aligned} $$

For $V_s \in \mathbb{R}^{m \times 1}$ being a vector, the provided indices correspond to the ones for a scalar-valued output (4.6). Using the Poincaré and Cramér-Rao inequality, we now derive derivative-based bounds for the goal-oriented Sobol' indices involving similar derivative-based quantities as in Theorem 3.5.

**PROPOSITION 4.1.** With the above notations, for any matrices $U_r \in \mathbb{R}^{d \times r}$ and $V_s \in \mathbb{R}^{m \times s}$ with orthogonal columns, we have

$$ (4.8) \qquad 1 - \bar{C}(X) \frac{\text{Tr}(U_{\perp}^{\top} H_X(V_s) U_{\perp})}{\text{Tr}(\text{Cov}(V_s^{\top} G(X)))} \le S^{\text{cl}}(U_r|V_s) \le 1 - \bar{c}(X) \frac{\|V_s^{\top} \mathbb{E}[\nabla G(X)] U_{\perp}\|_F^2}{\text{Tr}(\text{Cov}(V_s^{\top} G(X)))}, $$

$$ (4.9) \qquad \bar{c}(X) \frac{\|V_s^{\top} \mathbb{E}[\nabla G(X)] U_r\|_F^2}{\text{Tr}(\text{Cov}(V_s^{\top} G(X)))} \le S^{\text{tot}}(U_r|V_s) \le \bar{C}(X) \frac{\text{Tr}(U_r^{\top} H_X(V_s) U_r)}{\text{Tr}(\text{Cov}(V_s^{\top} G(X)))}, $$

where $U_{\perp} \in \mathbb{R}^{d \times (d-r)}$ is the orthogonal complement to $U_r$.
numerical experiments show that the EIM-based sensors often slightly outperform $\tau^*$ in terms of the true goal-oriented EIG rather than the lower bound. A reason could be that the sensors indexed by $\tilde{\tau}$ tend to be less localized than $\tau^*$. In the numerical experiments, we also apply the EIM on PCA bases to obtain sensor placements $\tilde{\tau}^{\text{PCA}}$, which we will refer to as PCA-EIM sensors.

**4.2. Goal-Oriented Sensitivity Analysis (GSA).** Global sensitivity analysis aims to assign an importance value to each set of input variables, reflecting its contribution to the output variations. It is often used for factor prioritization and/or fixing in complex models. Factor prioritization identifies the subset $\tau \subseteq \{1, \dots, d\}$, $|\tau| = r$ of input parameters causing the most output variation so that the parameters $X_{\tau} = (X_i)_{i \in \tau}$ can be prioritized e.g., in uncertainty reduction. Factor fixing, on the other hand, identifies the input set $\tau$ with the least influence on the output. Thus, $X_{\tau}$ can be fixed at a nominal value to reduce the complexity of, e.g., parameter estimation problems. For a scalar-valued square-integrable function $f: \mathbb{R}^d \to \mathbb{R}$, commonly used sensitivity indices are the closed Sobol' indices and the total Sobol' indices defined as

$$ (4.6) \qquad S_{\tau}^{\text{cl}} = \frac{\text{Var}(\mathbb{E}[f(X)|X_{\tau}])}{\text{Var}(f(X))} \quad \text{and} \quad S_{\tau}^{\text{tot}} = 1 - \frac{\text{Var}(\mathbb{E}[f(X)|X_{-\tau}])}{\text{Var}(f(X))}, $$

respectively, where $X_{-\tau} = (X_i)_{i \notin \tau}$ and the input variables $X_i$ for $i = 1, \dots, d$ are assumed to be independent. The closed Sobol' index quantifies the portion of output variance attributed to $X_{-\tau}$. A higher value implies a stronger influence of $X_{-\tau}$ on the output. The total Sobol' index encompasses output variation due to $X_{-\tau}$ along with variations caused by interactions of $X_{-\tau}$ with any of the remaining input variables. As a result, it holds that $S_{\tau}^{\text{cl}} \le S_{\tau}^{\text{tot}}$ and $1 = S_{\tau}^{\text{tot}} + S_{-\tau}^{\text{cl}}$. A lower value in total Sobol' index suggests that the parameters indexed by $\tau$ have little direct nor indirect impact on the output.

We consider, in particular, goal-oriented GSA where some outputs of interest $V_s^\top G(X)$ are prescribed. We use definitions of Sobol' indices for vector-valued functions in [24] and generalize them to the goal-oriented case as follows

$$ (4.7) \qquad \begin{aligned} S^{\text{cl}}(U_{\tau}|V_s) &= \frac{\text{Tr}(\text{Cov}(\mathbb{E}[V_s^{\top}G(X)|U_{\tau}^{\top}X]))}{\text{Tr}(\text{Cov}(V_s^{\top}G(X)))}, \\ S^{\text{tot}}(U_{\tau}|V_s) &= 1 - \frac{\text{Tr}(\text{Cov}(V_s^{\top}G(X)|U_{-\tau}^{\top}X))}{\text{Tr}(\text{Cov}(V_s^{\top}G(X)))}. \end{aligned} $$

For $V_s \in \mathbb{R}^{m \times 1}$ being a vector, the provided indices correspond to the ones for a scalar-valued output (4.6). Using the Poincaré and Cramér-Rao inequality, we now derive derivative-based bounds for the goal-oriented Sobol' indices involving similar derivative-based quantities as in Theorem 3.5.

**PROPOSITION 4.1.** With the above notations, for any matrices $U_r \in \mathbb{R}^{d \times r}$ and $V_s \in \mathbb{R}^{m \times s}$ with orthogonal columns, we have

$$ (4.8) \qquad 1 - \bar{C}(X) \frac{\text{Tr}(U_{\perp}^{\top} H_X(V_s) U_{\perp})}{\text{Tr}(\text{Cov}(V_s^{\top} G(X)))} \le S^{\text{cl}}(U_r|V_s) \le 1 - \bar{c}(X) \frac{\|V_s^{\top} \mathbb{E}[\nabla G(X)] U_{\perp}\|_F^2}{\text{Tr}(\text{Cov}(V_s^{\top} G(X)))}, $$

$$ (4.9) \qquad \bar{c}(X) \frac{\|V_s^{\top} \mathbb{E}[\nabla G(X)] U_r\|_F^2}{\text{Tr}(\text{Cov}(V_s^{\top} G(X)))} \le S^{\text{tot}}(U_r|V_s) \le \bar{C}(X) \frac{\text{Tr}(U_r^{\top} H_X(V_s) U_r)}{\text{Tr}(\text{Cov}(V_s^{\top} G(X)))}, $$

where $U_{\perp} \in \mathbb{R}^{d \times (d-r)}$ is the orthogonal complement to $U_r$.
The proof can be found in Appendix E. Similarly to the BOED problem, GSA struggles with an objective function that is expensive to evaluate combined with a difficult combinatorial optimization problem to identify an optimal index set $\tau$ (see discussion in Section 2). On the other hand, optimizing our derivative-based bound instead of the Sobol' indices provides us with a closed-form solution for the optimal index set. That is,

$$ (4.10) \qquad \min_{\substack{U_{-\tau_{\mathrm{cl}}} \in \mathbb{R}^{d \times (d-r)} \\ \tau_{\mathrm{cl}} \subset \{1, \dots, d\}, \, |\tau_{\mathrm{cl}}|=r}} \mathrm{Tr} \left( U_{-\tau_{\mathrm{cl}}}^{\top} H_X(V_s) U_{-\tau_{\mathrm{cl}}} \right) = \sum_{i \notin \tau_{\mathrm{cl}}^{*}} H_X(V_s)_{ii}, $$

where the index set $\tau_{\mathrm{cl}}^*$ maximizing the lower bound for the closed Sobol' index is given by the indices corresponding to the largest diagonal entries of $H_X(V_s)$. Correspondingly, we have

$$ (4.11) \qquad \min_{\substack{U_{\tau_{\mathrm{tot}}} \in \mathbb{R}^{d \times r} \\ \tau_{\mathrm{tot}} \subset \{1, \dots, d\}, \, |\tau_{\mathrm{tot}}|=r}} \mathrm{Tr} \left( U_{\tau_{\mathrm{tot}}}^{\top} H_X(V_s) U_{\tau_{\mathrm{tot}}} \right) = \sum_{i \in \tau_{\mathrm{tot}}^{*}} H_X(V_s)_{ii}, $$

where the index set $\tau_{\mathrm{tot}}^*$ minimizing the upper bound for the total Sobol' index is given by the indices corresponding to the smallest diagonal entries of $H_X(V_s)$. Note that we can handle index sets of any size for the same computational cost as the case of single-element index sets, given that only the diagonal elements of $H_X(V_s)$ need to be computed in either setting. Moreover, evaluating the respective upper and lower bound of the closed and total Sobol' index allows for complete certification of the induced error.

Remark 4.2. Our method generalises the derivative-based global sensitivity measures (DGSM) introduced in [50] which are defined as

$$ \nu_i = \mathbb{E} \left[ \left( \frac{\partial f}{\partial x_i}(X) \right)^2 \right] $$

for a scalar function $f: \mathbb{R}^d \to \mathbb{R}$. Setting $V_s \in \mathbb{R}^{m \times 1}$ to be a vector, the diagonal elements $H_X(V_s)_{ii}$ correspond to the DGSM $\nu_i$ of the scalar function $V_s^\top G(X)$. Our approach additionally enables sensitivity measures for goal-oriented vector-valued outputs and is not restricted to single-element sets $\tau = \{i\}$. The authors in [35] also derived the following DGSM-based bounds for the total Sobol' index where $X_i \sim \mathcal{N}(\mu_i, \sigma_i)$

$$ \frac{\sigma_i^4}{(\mu_i^2 + \sigma_i^2)\bar{V}} \omega_i^2 \le S_i^{\mathrm{tot}} \le \frac{\sigma_i^2}{\bar{V}} \nu_i, $$

where $\omega_i = \mathbb{E}[\frac{\partial f(X)}{\partial x_i}]$ and $\bar{V} = \mathrm{Var}(f(X))$. For normally distributed $X_i$, we have $\bar{c}(X) = \bar{C}(X) = \sigma_i^2$ (see Remark 3.7) and thus obtain a slightly improved lower bound $\frac{\sigma_i^2}{\bar{V}} \omega_i^2 \le S_i^{\mathrm{tot}} \le \frac{\sigma_i^2}{\bar{V}} \nu_i$. Note that lower-bounds for the total Sobol' index are generally receiving growing attention in GSA [46, 45].

**5. Numerical Experiments.** The following numerical experiments are computed with code provided at https://github.com/qchen95/Coupled-input-output-DR.

**5.1. Conditioned Diffusion.** Here, we consider a model for the movement of a diffusive particle acting under a double-well potential in molecular dynamics. The goal is to infer the driving force $dB_t$ applied to the particle given observations of its
The proof can be found in Appendix E. Similarly to the BOED problem, GSA struggles with an objective function that is expensive to evaluate combined with a difficult combinatorial optimization problem to identify an optimal index set $\tau$ (see discussion in Section 2). On the other hand, optimizing our derivative-based bound instead of the Sobol' indices provides us with a closed-form solution for the optimal index set. That is,

$$ (4.10) \qquad \min_{\substack{U_{-\tau_{\mathrm{cl}}} \in \mathbb{R}^{d \times (d-r)} \\ \tau_{\mathrm{cl}} \subset \{1, \dots, d\}, \ |\tau_{\mathrm{cl}}|=r}} \mathrm{Tr} \left( U_{-\tau_{\mathrm{cl}}}^{\top} H_X(V_s) U_{-\tau_{\mathrm{cl}}} \right) = \sum_{i \notin \tau_{\mathrm{cl}}^{*}} H_X(V_s)_{ii}, $$

where the index set $\tau_{\mathrm{cl}}^*$ maximizing the lower bound for the closed Sobol' index is given by the indices corresponding to the largest diagonal entries of $H_X(V_s)$. Correspondingly, we have

$$ (4.11) \qquad \min_{\substack{U_{\tau_{\mathrm{tot}}} \in \mathbb{R}^{d \times r} \\ \tau_{\mathrm{tot}} \subset \{1, \dots, d\}, \ |\tau_{\mathrm{tot}}|=r}} \mathrm{Tr} \left( U_{\tau_{\mathrm{tot}}}^{\top} H_X(V_s) U_{\tau_{\mathrm{tot}}} \right) = \sum_{i \in \tau_{\mathrm{tot}}^{*}} H_X(V_s)_{ii}, $$

where the index set $\tau_{\mathrm{tot}}^*$ minimizing the upper bound for the total Sobol' index is given by the indices corresponding to the smallest diagonal entries of $H_X(V_s)$. Note that we can handle index sets of any size for the same computational cost as the case of single-element index sets, given that only the diagonal elements of $H_X(V_s)$ need to be computed in either setting. Moreover, evaluating the respective upper and lower bound of the closed and total Sobol' index allows for complete certification of the induced error.

**Remark 4.2.** Our method generalises the derivative-based global sensitivity measures (DGSM) introduced in [50] which are defined as

$$ \nu_i = \mathbb{E} \left[ \left( \frac{\partial f}{\partial x_i}(X) \right)^2 \right] $$

for a scalar function $f: \mathbb{R}^d \to \mathbb{R}$. Setting $V_s \in \mathbb{R}^{m \times 1}$ to be a vector, the diagonal elements $H_X(V_s)_{ii}$ correspond to the DGSM $\nu_i$ of the scalar function $V_s^\top G(X)$. Our approach additionally enables sensitivity measures for goal-oriented vector-valued outputs and is not restricted to single-element sets $\tau = \{i\}$. The authors in [35] also derived the following DGSM-based bounds for the total Sobol' index where $X_i \sim \mathcal{N}(\mu_i, \sigma_i)$

$$ \frac{\sigma_i^4}{(\mu_i^2 + \sigma_i^2)\bar{V}} \omega_i^2 \leq S_i^{\mathrm{tot}} \leq \frac{\sigma_i^2}{\bar{V}} \nu_i, $$

where $\omega_i = \mathbb{E}[\frac{\partial f(X)}{\partial x_i}]$ and $\bar{V} = \mathrm{Var}(f(X))$. For normally distributed $X_i$, we have $\bar{c}(X) = \bar{C}(X) = \sigma_i^2$ (see Remark 3.7) and thus obtain a slightly improved lower bound $\frac{\sigma_i^2}{\bar{V}} \omega_i^2 \leq S_i^{\mathrm{tot}} \leq \frac{\sigma_i^2}{\bar{V}} \nu_i$. Note that lower-bounds for the total Sobol' index are generally receiving growing attention in GSA [46, 45].

**5. Numerical Experiments.** The following numerical experiments are computed with code provided at https://github.com/qchen95/Coupled-input-output-DR.

**5.1. Conditioned Diffusion.** Here, we consider a model for the movement of a diffusive particle acting under a double-well potential in molecular dynamics. The goal is to infer the driving force $dB_t$ applied to the particle given observations of its
path $u: [0, T] \to \mathbb{R}$ with final time $T = 1$. The particle movement is described by the solution of the stochastic differential equation

$$du_t = f(u_t) dt + dB_t, \quad u_0 = 0,$$

with non-linear drift function $f(u) = u(1 - u^2)/(1 + u^2)$ and $B$ denoting a Wiener process. We discretize the differential equation using an Euler-Maruyama scheme with time step $\Delta t = 10^{-2}$ resulting in 100 uniformly-spaced time steps. Noisy observations $Y$ are taken at each time step so that the forward model $G: \mathbb{R}^d \to \mathbb{R}^m$ is defined as the mapping from the incremental driving force $X$ at $d = 100$ times to the particle position $u$ at $m = 100$ times. We further assume a standard normal prior $N(0, I_d)$ for the driving force and i.i.d. additive observational noise with distribution $N(0, 0.1)$. Figures 1a and 1b illustrate some realizations of the particle path $u$.

In this problem, we aim to identify optimal subspaces and indices of the inputs and outputs of $G$ that reduce the dimensions of the Bayesian inference problem for inferring $X$ from $Y$. To compute the diagnostic matrices $H_X(V_s)$ (3.13) and $H_Y(U_r)$ (3.14), we use $M = 10000$ prior samples and set $r = s = 10$, if not stated otherwise. The alternating eigendecomposition algorithm is run for a fixed number of 10 iterations and initiated from a random $V_s \in \mathbb{R}^{m \times s}$ with orthonormal columns when no goal is specified. We recall that the optimal subspace is given without alternation when only solving for one of $U_r$ or $V_s$.

<figure><img src="image_2.png" /><figcaption>Fig. 1: Illustration of the inputs and outputs of interest (i.e., goals) within the specified orange time intervals, and their corresponding goal-oriented solutions for the vectors defining the reduced output and input spaces, respectively. Notice that the solutions correctly reflect the time causality of the problem.</figcaption></figure>
path $u: [0, T] \to \mathbb{R}$ with final time $T = 1$. The particle movement is described by the solution of the stochastic differential equation

$$du_t = f(u_t) dt + dB_t, \quad u_0 = 0,$$

with non-linear drift function $f(u) = u(1 - u^2)/(1 + u^2)$ and $B$ denoting a Wiener process. We discretize the differential equation using an Euler-Maruyama scheme with time step $\Delta t = 10^{-2}$ resulting in 100 uniformly-spaced time steps. Noisy observations $Y$ are taken at each time step so that the forward model $G: \mathbb{R}^d \to \mathbb{R}^m$ is defined as the mapping from the incremental driving force $X$ at $d = 100$ times to the particle position $u$ at $m = 100$ times. We further assume a standard normal prior $N(0, I_d)$ for the driving force and i.i.d. additive observational noise with distribution $N(0, 0.1)$. Figures 1a and 1b illustrate some realizations of the particle path $u$.

In this problem, we aim to identify optimal subspaces and indices of the inputs and outputs of $G$ that reduce the dimensions of the Bayesian inference problem for inferring $X$ from $Y$. To compute the diagnostic matrices $H_X(V_s)$ (3.13) and $H_Y(U_r)$ (3.14), we use $M = 10000$ prior samples and set $r = s = 10$, if not stated otherwise. The alternating eigendecomposition algorithm is run for a fixed number of 10 iterations and initiated from a random $V_s \in \mathbb{R}^{m \times s}$ with orthonormal columns when no goal is specified. We recall that the optimal subspace is given without alternation when only solving for one of $U_r$ or $V_s$.

<figure><img src="image_3.png" /><figcaption>Fig. 1: Illustration of the inputs and outputs of interest (i.e., goals) within the specified orange time intervals, and their corresponding goal-oriented solutions for the vectors defining the reduced output and input spaces, respectively. Notice that the solutions correctly reflect the time causality of the problem.</figcaption></figure>
In Figure 1, we observe solutions of the goal-oriented problems. The goals $U_r$ and $V_s$, illustrated in Figure 1a and Figure 1b, project onto the driving force and, respectively, the particle path observations within the specified orange time intervals. For the goal $U_r$, the resulting observation space bases $V_s^*(U_r)$ (3.16) in Figure 1c are zero up to the interval of interest. This behavior can be explained as the driving force within the prescribed time interval does not influence the particle position at earlier times. Similarly, the parameter space bases $U_r^*(V_s)$ (3.17) in Figure 1d derived from the goal $V_s$ do not extend beyond the observed time interval. An explanation is that observations of the particle's path do not provide information about the driving force at future times. Thus, our goal-oriented solutions appear to capture the time causality within the inference problem correctly.

<figure><img src="image_3.png" /><figcaption>Fig. 2: Goal-oriented total Sobol' indices for each index set as compared with the derived derivative-based bounds. The bounds tightly track the true value for the Sobol' indices. The bounds also take into account the time causality of the problem.</figcaption></figure>

In Figure 2, we present the goal-oriented total Sobol' indices for the output of interest $V_s^\top G(X)$ illustrated in Figure 1d. For each coordinate $i = 1, \dots, d$ of the input parameter, i.e., the driving force at time $i \cdot \Delta t$, we compute $S^{\text{tot}}(U_{\{i\}|V_s})$ along with the derivative-based bounds in (4.9). We present single-element input sets $\tau = \{i\}$, $1 \le i \le d$ for easier illustration. However, we recall that the derived bounds can be computed for arbitrary parameter sets by summing up the elementwise values according to (4.11). Both the upper and lower bounds in Figure 2 tightly track the total Sobol' indices. The bounds also capture the problem's temporal causality, i.e., driving forces at later times exert no influence on the model's output variations during the time interval of interest.

Figure 3 shows solutions of the coupled dimension reduction problem (top) in comparison with dimension reduction by applying PCA separately on the parameter and data space (bottom) [28]. The coupled modes of the parameter space focus on approximating the driving force at initial times, reflecting the strong influence of the force at these times on the overall particle path. Contrarily, the modes of the data space focus on approximating the particle path at final times. This phenomenon reflects the fact that the final position of the particle reveals more about the overall driving force. The coupled eigenvectors of both spaces exhibit notable disparities from those derived through PCA. The PCA modes of the input space reflect the white noise behavior of the driving force. The PCA modes of the observation space are more global as compared to the coupled ones and moreover they do not exhibit awareness of the inverse problem.

Next, we consider the convergence of the alternating eigendecomposition algo-
In Figure 1, we observe solutions of the goal-oriented problems. The goals $U_r$ and $V_s$, illustrated in Figure 1a and Figure 1b, project onto the driving force and, respectively, the particle path observations within the specified orange time intervals. For the goal $U_r$, the resulting observation space bases $V_s^*(U_r)$ (3.16) in Figure 1c are zero up to the interval of interest. This behavior can be explained as the driving force within the prescribed time interval does not influence the particle position at earlier times. Similarly, the parameter space bases $U_r^*(V_s)$ (3.17) in Figure 1d derived from the goal $V_s$ do not extend beyond the observed time interval. An explanation is that observations of the particle's path do not provide information about the driving force at future times. Thus, our goal-oriented solutions appear to capture the time causality within the inference problem correctly.

<figure><img src="image_4.png" /><figcaption>Fig. 2: Goal-oriented total Sobol' indices for each index set as compared with the derived derivative-based bounds. The bounds tightly track the true value for the Sobol' indices. The bounds also take into account the time causality of the problem.</figcaption></figure>

In Figure 2, we present the goal-oriented total Sobol' indices for the output of interest $V_s^\top G(X)$ illustrated in Figure 1d. For each coordinate $i = 1, \dots, d$ of the input parameter, i.e., the driving force at time $i \cdot \Delta t$, we compute $S^{\text{tot}}(U_{\{i\}|V_s})$ along with the derivative-based bounds in (4.9). We present single-element input sets $\tau = \{i\}$, $1 \le i \le d$ for easier illustration. However, we recall that the derived bounds can be computed for arbitrary parameter sets by summing up the elementwise values according to (4.11). Both the upper and lower bounds in Figure 2 tightly track the total Sobol' indices. The bounds also capture the problem's temporal causality, i.e., driving forces at later times exert no influence on the model's output variations during the time interval of interest.

Figure 3 shows solutions of the coupled dimension reduction problem (top) in comparison with dimension reduction by applying PCA separately on the parameter and data space (bottom) [28]. The coupled modes of the parameter space focus on approximating the driving force at initial times, reflecting the strong influence of the force at these times on the overall particle path. Contrarily, the modes of the data space focus on approximating the particle path at final times. This phenomenon reflects the fact that the final position of the particle reveals more about the overall driving force. The coupled eigenvectors of both spaces exhibit notable disparities from those derived through PCA. The PCA modes of the input space reflect the white noise behavior of the driving force. The PCA modes of the observation space are more global as compared to the coupled ones and moreover they do not exhibit awareness of the inverse problem.

Next, we consider the convergence of the alternating eigendecomposition algo-
<figure><img src="image_4.png" /><figcaption>Fig. 3: Comparison of subspaces from coupled and PCA-based dimension reduction.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4: Convergence of alternating eigendecomposition for 10 random initializations.</figcaption></figure>

Figure 4 depicts the relative $L^2$-error (3.3) as well as the relative upper (3.8) and lower bound (3.9) over the iterations of the alternating eigendecomposition for 10 random initializations of the matrices. We observe that minimizing the upper error bound indeed decreases the true $L^2$-error. Moreover, the alternating eigendecomposition converges after just three iterations independent of the initiation.

Figure 5 shows the convergence of the relative $L^2$-error and upper bound with increasing subspace dimensions $r$ and $s$. We compare the errors of our coupled dimension reduction method with other dimension reduction methods in the data and parameter spaces such as PCA, CCA [29] and the joint dimension reduction method [12] discussed in Remark 3.11. Our coupled method has very similar error and upper
<figure><img src="image_5.png" /><figcaption>Fig. 3: Comparison of subspaces from coupled and PCA-based dimension reduction.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 4: Convergence of alternating eigendecomposition for 10 random initializations.</figcaption></figure>

Figure 4 depicts the relative $L^2$-error (3.3) as well as the relative upper (3.8) and lower bound (3.9) over the iterations of the alternating eigendecomposition for 10 random initializations of the matrices. We observe that minimizing the upper error bound indeed decreases the true $L^2$-error. Moreover, the alternating eigendecomposition converges after just three iterations independent of the initiation.

Figure 5 shows the convergence of the relative $L^2$-error and upper bound with increasing subspace dimensions $r$ and $s$. We compare the errors of our coupled dimension reduction method with other dimension reduction methods in the data and parameter spaces such as PCA, CCA [29] and the joint dimension reduction method [12] discussed in Remark 3.11. Our coupled method has very similar error and upper
<figure><img src="image_6.png" /><figcaption>Fig. 5: Comparison of model error and upper bound convergence across different dimension reduction methods as subspace dimensions $r, s$ increase.</figcaption></figure>

bound values as the joint dimension reduction method. Both decrease much faster with increasing subspace dimensions and do not get stuck on a plateau like CCA.

**5.2. Burgers' Equation.** In this subsection, we demonstrate our goal-oriented BOED framework on the one-dimensional non-linear Burgers' equation

$$ \frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} - D \frac{\partial u}{\partial x^2} = 0, $$

with $t \in [0, 1]$, $x \in [0, 1]$, $D = 0.001$ and periodic boundaries. Our parameters are the initial conditions $u_0(x)$. We assume a normal prior distribution for $u_0(x)$ with mean $\mu(x) = \frac{1}{\sqrt{2\pi}} \exp(-\frac{1}{2}\|0.5 - x\|^2)$ and Matérn covariance with correlation length $l = 0.1$ and smoothness $\nu = 2.5$. Using a uniform space discretization our parameters are $(u_0(x_1), \dots, u_0(x_d)) \in \mathbb{R}^d$ with $d = 100$. As data $y \in \mathbb{R}^m$, we use state observations at the final time $T = 0.1$ with i.i.d. measurement noise $\eta \sim N(0, 0.01 I_m)$ so that $m = 100$. Therefore, $G: \mathbb{R}^d \to \mathbb{R}^m$ maps the discretized initial condition $u_0(x)$ to the discretized solution $u(x, T)$ at the final time. The left subplot in Figure 6 shows the prior mean $\mu(x)$ together with several prior samples and also visualizes each of the goals $U_r^{(i)}$, which project the initial conditions $X$ onto the grid points $i, \dots, i+r$ marked by an orange bar. The right subplot in Figure 6 illustrates the space-time model solutions $u(x, t)$ with the prior mean as initial condition $u_0(x) = \mu(x)$. We consider $r = 5$ and the sensor placement for $s = 10$ sensors. To compute the diagnostic matrices, we use $M = 1000$ prior samples, and the goal-oriented true EIG (4.3) is computed via the double-loop importance sampling estimator from [14] with 30 inner and $M$ outer loop samples.

In Figure 7, we present the goal-oriented EIG for different goals $U_r^{(i)}$ and sensor placement methods. For each goal, we select sensors based on three different methods: the optimal sensor placement $\tau^*$ in (4.4), the goal-oriented EIM-based sensors $\tilde{\tau}$ and the PCA-EIM sensors $\tilde{\tau}^{\text{PCA}}$. The violin plots also show the EIG values for 100 random sensor placements. Both sensors $\tau^*$ and $\tilde{\tau}$ almost always outperform the random sensors, i.e., yielding larger EIG. Moreover, the EIM-based sensors $\tilde{\tau}$ sometimes slightly outperform $\tau^*$. On the other hand, the PCA-EIM sensors $\tilde{\tau}^{\text{PCA}}$ are often only as good as average random selections, especially for goals $U_r^{(i)}$ with lower indices $i < 45$. This phenomenon can be explained by looking at the actual sensor positions
<figure><img src="image_7.png" /><figcaption>Fig. 5: Comparison of model error and upper bound convergence across different dimension reduction methods as subspace dimensions r, s increase.</figcaption></figure>

bound values as the joint dimension reduction method. Both decrease much faster with increasing subspace dimensions and do not get stuck on a plateau like CCA.

**5.2. Burgers' Equation.** In this subsection, we demonstrate our goal-oriented BOED framework on the one-dimensional non-linear Burgers' equation

$$ \frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} - D \frac{\partial u}{\partial x^2} = 0, $$

with $t \in [0, 1]$, $x \in [0, 1]$, $D = 0.001$ and periodic boundaries. Our parameters are the initial conditions $u_0(x)$. We assume a normal prior distribution for $u_0(x)$ with mean $\mu(x) = \frac{1}{\sqrt{2\pi}} \exp(-\frac{1}{2}\|0.5 - x\|^2)$ and Matérn covariance with correlation length $l = 0.1$ and smoothness $\nu = 2.5$. Using a uniform space discretization our parameters are $(u_0(x_1), \dots, u_0(x_d)) \in \mathbb{R}^d$ with $d = 100$. As data $y \in \mathbb{R}^m$, we use state observations at the final time $T = 0.1$ with i.i.d. measurement noise $\eta \sim \mathcal{N}(0, 0.01 I_m)$ so that $m = 100$. Therefore, $G: \mathbb{R}^d \to \mathbb{R}^m$ maps the discretized initial condition $u_0(x)$ to the discretized solution $u(x, T)$ at the final time. The left subplot in Figure 6 shows the prior mean $\mu(x)$ together with several prior samples and also visualizes each of the goals $U_r^{(i)}$, which project the initial conditions $X$ onto the grid points $i, \dots, i+r$ marked by an orange bar. The right subplot in Figure 6 illustrates the space-time model solutions $u(x, t)$ with the prior mean as initial condition $u_0(x) = \mu(x)$. We consider $r = 5$ and the sensor placement for $s = 10$ sensors. To compute the diagnostic matrices, we use $M = 1000$ prior samples, and the goal-oriented true EIG (4.3) is computed via the double-loop importance sampling estimator from [14] with 30 inner and $M$ outer loop samples.

In Figure 7, we present the goal-oriented EIG for different goals $U_r^{(i)}$ and sensor placement methods. For each goal, we select sensors based on three different methods: the optimal sensor placement $\tau^*$ in (4.4), the goal-oriented EIM-based sensors $\tilde{\tau}$ and the PCA-EIM sensors $\tilde{\tau}^{\text{PCA}}$. The violin plots also show the EIG values for 100 random sensor placements. Both sensors $\tau^*$ and $\tilde{\tau}$ almost always outperform the random sensors, i.e., yielding larger EIG. Moreover, the EIM-based sensors $\tilde{\tau}$ sometimes slightly outperform $\tau^*$. On the other hand, the PCA-EIM sensors $\tilde{\tau}^{\text{PCA}}$ are often only as good as average random selections, especially for goals $U_r^{(i)}$ with lower indices $i < 45$. This phenomenon can be explained by looking at the actual sensor positions
<figure><img src="image_7.png" /><figcaption>Fig. 6: Left: Illustration of the prior mean $\mu(x)$ (solid line) and samples (transparent lines) of the initial condition parameter with orange bars marking the quantities of interest $U_r^{(i)}$. Right: Solutions of Burgers' equation with the prior mean as initial condition, i.e., $u_0(x) = \mu(x)$.</figcaption></figure>

in Figure 8. Since the solution has a wavefront that is moving to the right, the final state has strong features on the right side of the domain. As a result, the PCA-EIM sensors are concentrated on the right side of the domain, regardless of the goal. Thus, the PCA-EIM sensors seem to be more informative for parameters of interests $U_r^{(i)}$ that also lie on the right side of the domain (large $i$) and less so for parameters of interest located on the left side of the domain (small $i$). Figure 8 also shows that the EIM-based indices $\tilde{\tau}$ are often more spread than the optimal sensors $\tau^*$, which could explain their better performance.

<figure><img src="image_8.png" /><figcaption>Fig. 7: Goal-oriented EIG values of different goals $U_r^{(i)}$ and sensor placements. Violin plots represent 100 random sensor placements.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Fig. 6: Left: Illustration of the prior mean $\mu(x)$ (solid line) and samples (transparent lines) of the initial condition parameter with orange bars marking the quantities of interest $U_r^{(i)}$. Right: Solutions of Burgers' equation with the prior mean as initial condition, i.e., $u_0(x) = \mu(x)$.</figcaption></figure>

in Figure 8. Since the solution has a wavefront that is moving to the right, the final state has strong features on the right side of the domain. As a result, the PCA-EIM sensors are concentrated on the right side of the domain, regardless of the goal. Thus, the PCA-EIM sensors seem to be more informative for parameters of interests $U_r^{(i)}$ that also lie on the right side of the domain (large $i$) and less so for parameters of interest located on the left side of the domain (small $i$). Figure 8 also shows that the EIM-based indices $\tilde{\tau}$ are often more spread than the optimal sensors $\tau^*$, which could explain their better performance.

<figure><img src="image_9.png" /><figcaption>Fig. 7: Goal-oriented EIG values of different goals $U_r^{(i)}$ and sensor placements. Violin plots represent 100 random sensor placements.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Fig. 8: Goal-oriented modes $V_s^*(U_r)$ (denoted by lines) and sensor placements (denoted by markers at select indices below each plot) for four different parameters of interest $U_r^\top X$ marked by orange bars.</figcaption></figure>

**6. Conclusions.** We derived a coupled input-output dimension reduction method that is easily computed using aN alternating eigendecomposition of two diagnostic matrices. The reduced spaces are identified given only access to matrix-vector products of the model's gradient. Our method supports goal-oriented dimension reduction when some quantity of interest is prescribed in the input or the output space. We demonstrate our dimension reduction framework on Bayesian optimal experimental design and global sensitivity analysis, which can be seen as selecting the most important components of each variable instead of subspaces. In both applications, we not only bypass the expensive evaluation of the objective function for each task but also a combinatorial optimization problem by computing the relevant components from the largest diagonal entries of the diagnostic matrices.

Next, we explain a few possible directions for further developments. First, the convergence of the alternating eigendecomposition algorithm and the (adaptive) choice of the ranks $r$ and $s$ are important open questions. Second, we plan to examine quantities of interest described by some nonlinear functions $\psi_r: \mathbb{R}^d \to \mathbb{R}^r$ for the input and/or $\phi_s: \mathbb{R}^m \to \mathbb{R}^s$ for the output. By the chain rule, the diagnostic matrices become

$$ 
\begin{aligned}
H_X(\phi_s) &= \mathbb{E}[\nabla G(X)^\top \nabla \phi_s(G(X))^\top \nabla \phi_s(G(X)) \nabla G(X)], \\
H_Y(\psi_r) &= \mathbb{E}[\nabla G(X) \nabla \psi_r(X)^\top \nabla \psi_r(X) \nabla G(X)^\top].
\end{aligned}
$$

The error analysis associated with such diagnostics is a natural question for future work. Third, our coupled dimension reduction framework can also be extended to operator learning [38, 31] where $G$ is considered as operator $G: \mathcal{U} \to \mathcal{V}$ acting on function
<figure><img src="image_10.png" /><figcaption>Fig. 8: Goal-oriented modes $V_s^*(U_r)$ (denoted by lines) and sensor placements (denoted by markers at select indices below each plot) for four different parameters of interest $U_r^\top X$ marked by orange bars.</figcaption></figure>

**6. Conclusions.** We derived a coupled input-output dimension reduction method that is easily computed using aN alternating eigendecomposition of two diagnostic matrices. The reduced spaces are identified given only access to matrix-vector products of the model's gradient. Our method supports goal-oriented dimension reduction when some quantity of interest is prescribed in the input or the output space. We demonstrate our dimension reduction framework on Bayesian optimal experimental design and global sensitivity analysis, which can be seen as selecting the most important components of each variable instead of subspaces. In both applications, we not only bypass the expensive evaluation of the objective function for each task but also a combinatorial optimization problem by computing the relevant components from the largest diagonal entries of the diagnostic matrices.

Next, we explain a few possible directions for further developments. First, the convergence of the alternating eigendecomposition algorithm and the (adaptive) choice of the ranks $r$ and $s$ are important open questions. Second, we plan to examine quantities of interest described by some nonlinear functions $\psi_r: \mathbb{R}^d \to \mathbb{R}^r$ for the input and/or $\phi_s: \mathbb{R}^m \to \mathbb{R}^s$ for the output. By the chain rule, the diagnostic matrices become

$$H_X(\phi_s) = \mathbb{E}[\nabla G(X)^\top \nabla \phi_s(G(X))^\top \nabla \phi_s(G(X)) \nabla G(X)], \\ H_Y(\psi_r) = \mathbb{E}[\nabla G(X) \nabla \psi_r(X)^\top \nabla \psi_r(X) \nabla G(X)^\top].$$

The error analysis associated with such diagnostics is a natural question for future work. Third, our coupled dimension reduction framework can also be extended to operator learning [38, 31] where $G$ is considered as operator $G: \mathcal{U} \to \mathcal{V}$ acting on function
spaces $\mathcal{U}, \mathcal{V}$, typically separable Hilbert spaces. This extension is possible since the Poincaré and Cramér-Rao inequalities are known to hold in infinite dimensions.

**Appendix A. Linking the $L^2$-error of $\tilde{G}$.** Recall Proposition 1 from [12].

**PROPOSITION A.1.** Set $\tilde{\pi}_{X|Y} = \pi_{X_r|Y_s} \pi_{X_\perp|X_r}$, where $X_\perp = U_\perp^\top X$ is the orthogonal complement of $X_r$ with $U_\perp \in \mathbb{R}^{d \times (d-r)}$ spanning $\text{Ker}(U_\perp^\top)$. Then, we have

$$ \mathbb{E}_Y [\mathbb{D}_{\text{KL}} (\pi_{X|Y} \| \tilde{\pi}_{X|Y})] \le \mathbb{E}_Y [\mathbb{D}_{\text{KL}} (\pi_{X|Y} \| \hat{\pi}_{X|Y})] $$

for any approximate posterior of the form $\hat{\pi}_{X|Y}(x|y) = f_1(x_r, y_s)f_2(x_\perp, x_r)$ with non-negative functions $f_1, f_2$.

**A.1. Proof of Lemma 3.1.** First note that the normalized form of $\tilde{\pi}_{X|Y} \propto \pi_{Y_s|X_r} \pi_X$ in (2.2) can be written as

$$ \tilde{\pi}_{X|Y} = \frac{\pi_{Y_s|X_r} \pi_X}{\int \pi_{Y_s|X_r} \pi_X dx} = \frac{\pi_{Y_s|X_r} \pi_X}{\pi_{Y_s}} = \frac{\pi_{Y_s, X_r} \pi_X}{\pi_{X_r} \pi_{Y_s}} = \pi_{X_r|Y_s} \pi_{X_\perp|X_r} $$

Thus, it corresponds to the definition in Proposition A.1 and we can use it to obtain

$$ (A.1) \qquad \mathbb{E}_Y [\mathbb{D}_{\text{KL}} (\pi_{X|Y} \| \tilde{\pi}_{X|Y})] \le \mathbb{E}_Y [\mathbb{D}_{\text{KL}} (\pi_{X|Y} \| \hat{\pi}_{X|Y})] $$

where $\hat{\pi}_{X|Y}(x|y) \propto \exp(-\frac{1}{2\sigma^2}\|\tilde{G}(x) - y\|^2)\pi_X(x)$ satisfies the definition in Proposition A.1 with $f_1(x_r, y_s) = \frac{1}{Z} \exp(-\frac{1}{2\sigma^2}\|\tilde{g}(x_r) + \tilde{G}_\perp - y_s\|^2)$ and $f_2(x_\perp, x_r) = \pi_X(x)$. Here, $Z \in \mathbb{R}$ represents a normalization constant independent of $x$ and $y$. By Proposition 4.1 in [20] we have the following bound for any posterior $\pi_{X|Y}$ with Gaussian likelihood and its approximation $\hat{\pi}_{X|Y}$ with an approximate forward model

$$ (A.2) \qquad \mathbb{E}_Y [\mathbb{D}_{\text{KL}} (\pi_{X|Y} \| \hat{\pi}_{X|Y})] \le \frac{1}{2\sigma^2} \mathbb{E}_Y [\|G(x) - \tilde{G}(x)\|^2]. $$

Combining (A.1) and (A.2) finishes the proof.

**A.2. Proof of Lemma 3.2.** We have from the previous proof that $\tilde{\pi}_{X|Y} = \pi_{X_r|Y_s} \pi_{X_\perp|X_r}$. We get for the goal-oriented EIG

$$ \begin{aligned} \Phi(V_s|U_r) &= \mathbb{E}_Y [\mathbb{D}_{\text{KL}} (\pi_{X_r|Y_s} \| \pi_{X_r})] = \int \log \left( \frac{\pi_{X_r|Y_s}}{\pi_{X_r}} \right) d\pi_{XY} \\ &= \int \log \left( \frac{\tilde{\pi}_{X|Y}}{\pi_{X_r} \pi_{X_\perp|X_r}} \right) d\pi_{XY} \\ &= \int \log \left( \frac{\pi_{Y|X}}{\pi_Y} \right) d\pi_{XY} - \int \log \left( \frac{\pi_{X|Y}}{\tilde{\pi}_{X|Y}} \right) d\pi_{XY} \\ &= \int \left( \int \log \left( \frac{\pi_{X|Y}}{\pi_X} \right) d\pi_{X|Y} \right) d\pi_Y - \int \left( \log \left( \frac{\pi_{X|Y}}{\tilde{\pi}_{X|Y}} \right) d\pi_{X|Y} \right) d\pi_Y \\ &= \mathbb{E}_Y [\mathbb{D}_{\text{KL}} (\pi_{X|Y} \| \pi_X)] - \mathbb{E}_Y [\mathbb{D}_{\text{KL}} (\pi_{X|Y} \| \tilde{\pi}_{X|Y})]. \end{aligned} $$

Applying Lemma 3.1 on the last term concludes the proof.

**Appendix B. Proof of Proposition 3.3.** Let us rewrite

$$ \begin{aligned} \mathbb{E}[\|G(X) - \tilde{G}(X)\|^2] &= \mathbb{E}[\|(G(X) - G^*(X)) + (G^*(X) - \tilde{G}(X))\|^2] \\ &= \mathbb{E}[\|G(X) - G^*(X)\|^2] + \mathbb{E}[\|G^*(X) - \tilde{G}(X)\|^2] \\ &\quad + 2\mathbb{E}[(G(X) - G^*(X))^{\top}(G^*(X) - \tilde{G}(X))] \end{aligned} $$
spaces $\mathcal{U}, \mathcal{V}$, typically separable Hilbert spaces. This extension is possible since the Poincaré and Cramér-Rao inequalities are known to hold in infinite dimensions.

**Appendix A. Linking the $L^2$-error of $\tilde{G}$.** Recall Proposition 1 from [12].

**PROPOSITION A.1.** Set $\tilde{\pi}_{X|Y} = \pi_{X_r|Y_s} \pi_{X_{\perp}|X_r}$, where $X_{\perp} = U_{\perp}^{\top} X$ is the orthogonal complement of $X_r$ with $U_{\perp} \in \mathbb{R}^{d \times (d-r)}$ spanning $\text{Ker}(U_r^{\top})$. Then, we have

$$ \mathbb{E}_Y [\mathbb{D}_{\text{KL}} (\pi_{X|Y} \| \tilde{\pi}_{X|Y})] \leq \mathbb{E}_Y [\mathbb{D}_{\text{KL}} (\pi_{X|Y} \| \hat{\pi}_{X|Y})] $$

for any approximate posterior of the form $\hat{\pi}_{X|Y}(x|y) = f_1(x_r, y_s)f_2(x_{\perp}, x_r)$ with non-negative functions $f_1, f_2$.

**A.1. Proof of Lemma 3.1.** First note that the normalized form of $\tilde{\pi}_{X|Y} \propto \pi_{Y_s|X_r} \pi_X$ in (2.2) can be written as

$$ \tilde{\pi}_{X|Y} = \frac{\pi_{Y_s|X_r} \pi_X}{\int \pi_{Y_s|X_r} \pi_X dx} = \frac{\pi_{Y_s|X_r} \pi_X}{\pi_{Y_s}} = \frac{\pi_{Y_s, X_r} \pi_X}{\pi_{X_r} \pi_{Y_s}} = \pi_{X_r|Y_s} \pi_{X_{\perp}|X_r} $$

Thus, it corresponds to the definition in Proposition A.1 and we can use it to obtain

$$ (A.1) \qquad \mathbb{E}_Y [\mathbb{D}_{\text{KL}} (\pi_{X|Y} \| \tilde{\pi}_{X|Y})] \leq \mathbb{E}_Y [\mathbb{D}_{\text{KL}} (\pi_{X|Y} \| \hat{\pi}_{X|Y})] $$

where $\hat{\pi}_{X|Y}(x|y) \propto \exp(-\frac{1}{2\sigma^2}\|\tilde{G}(x) - y\|^2)\pi_X(x)$ satisfies the definition in Proposition A.1 with $f_1(x_r, y_s) = \frac{1}{Z} \exp(-\frac{1}{2\sigma^2}\|\tilde{g}(x_r) + \tilde{G}_{\perp} - y_s\|^2)$ and $f_2(x_{\perp}, x_r) = \pi_X(x)$. Here, $Z \in \mathbb{R}$ represents a normalization constant independent of $x$ and $y$. By Proposition 4.1 in [20] we have the following bound for any posterior $\pi_{X|Y}$ with Gaussian likelihood and its approximation $\hat{\pi}_{X|Y}$ with an approximate forward model

$$ (A.2) \qquad \mathbb{E}_Y [\mathbb{D}_{\text{KL}} (\pi_{X|Y} \| \hat{\pi}_{X|Y})] \leq \frac{1}{2\sigma^2} \mathbb{E}_Y [\|G(x) - \tilde{G}(x)\|^2]. $$

Combining (A.1) and (A.2) finishes the proof.

**A.2. Proof of Lemma 3.2.** We have from the previous proof that $\tilde{\pi}_{X|Y} = \pi_{X_r|Y_s} \pi_{X_{\perp}|X_r}$. We get for the goal-oriented EIG

$$ \begin{aligned} \Phi(V_s|U_r) &= \mathbb{E}_Y [\mathbb{D}_{\text{KL}} (\pi_{X_r|Y_s} \| \pi_{X_r})] = \int \log \left( \frac{\pi_{X_r|Y_s}}{\pi_{X_r}} \right) d\pi_{XY} \\ &= \int \log \left( \frac{\tilde{\pi}_{X|Y}}{\pi_{X_r} \pi_{X_{\perp}|X_r}} \right) d\pi_{XY} \\ &= \int \log \left( \frac{\pi_{Y|X}}{\pi_Y} \right) d\pi_{XY} - \int \log \left( \frac{\pi_{X|Y}}{\tilde{\pi}_{X|Y}} \right) d\pi_{XY} \\ &= \int \left( \int \log \left( \frac{\pi_{X|Y}}{\pi_X} \right) d\pi_{X|Y} \right) d\pi_Y - \int \left( \log \left( \frac{\pi_{X|Y}}{\tilde{\pi}_{X|Y}} \right) d\pi_{X|Y} \right) d\pi_Y \\ &= \mathbb{E}_Y [\mathbb{D}_{\text{KL}} (\pi_{X|Y} \| \pi_X)] - \mathbb{E}_Y [\mathbb{D}_{\text{KL}} (\pi_{X|Y} \| \tilde{\pi}_{X|Y})]. \end{aligned} $$

Applying Lemma 3.1 on the last term concludes the proof.

**Appendix B. Proof of Proposition 3.3.** Let us rewrite

$$ \begin{aligned} \mathbb{E}[\|G(X) - \tilde{G}(X)\|^2] &= \mathbb{E}[\|(G(X) - G^*(X)) + (G^*(X) - \tilde{G}(X))\|^2] \\ &= \mathbb{E}[\|G(X) - G^*(X)\|^2] + \mathbb{E}[\|G^*(X) - \tilde{G}(X)\|^2] \\ &\quad + 2\mathbb{E}[(G(X) - G^*(X))^{\top}(G^*(X) - \tilde{G}(X))] \end{aligned} $$
We show that the last term vanishes so that for any $\tilde{G}$ it holds $\mathbb{E}[\|G(X) - \tilde{G}(X)\|^2] \ge \mathbb{E}[\|G(X) - G^*(X)\|^2]$. This proves that $G^*$ is a minimizer of $\tilde{G} \mapsto \mathbb{E}[\|G(X) - \tilde{G}(X)\|^2]$. Let $V_{\perp} \in \mathbb{R}^{m \times (m-s)}$ be the orthogonal complement of $V_s$. Without loss of generality, we assume that $\tilde{G}_{\perp} \in \text{Ker}(V_s)$ so that $V_{\perp}V_{\perp}^{\top}\tilde{G}_{\perp} = \tilde{G}_{\perp}$. By construction, we have

$$ 
\begin{aligned}
G(X) - G^*(X) &= V_s(V_s^\top G(X) - g^*(U_r^\top X)) + V_\perp V_\perp^\top (G(X) - \mathbb{E}[G(X)]), \\
G^*(X) - \tilde{G}(X) &= V_s(g^*(U_r^\top X) - \tilde{g}(U_r^\top X)) + V_\perp V_\perp^\top (G_\perp^* - \tilde{G}_\perp).
\end{aligned}
$$

Thus, we get that

$$ 
\begin{aligned}
\mathbb{E}[(G(X) - G^*(X))^\top (G^*(X) - \tilde{G}(X))] \\
&= \mathbb{E}[(V_s^\top G(X) - g^*(U_r^\top X))^\top (g^*(U_r^\top X) - \tilde{g}(U_r^\top X))] \\
&\quad + \underbrace{\mathbb{E}[(G(X) - \mathbb{E}[G(X)])^\top V_\perp V_\perp^\top (G_\perp^* - \tilde{G}_\perp)]}_{=0},
\end{aligned}
$$

where the last term vanishes by taking the expectation of $G(X) - \mathbb{E}[G(X)]$. Next, we exploit the fundamental property of the conditional expectation $\mathbb{E}[G(X)f(U_r^\top X)] = \mathbb{E}[\mathbb{E}[G(X)|U_r^\top X]f(U_r^\top X)]$ which holds for any function $f: \mathbb{R}^r \to \mathbb{R}$. Considering the definition of $g^*$ and setting $f(U_r^\top X) = V_s(g^*(U_r^\top X) - \tilde{g}(U_r^\top X))$, we obtain

$$ \mathbb{E}[(G(X) - G^*(X))^\top (G^*(X) - \tilde{G}(X))] = \mathbb{E}[(G(X) - \mathbb{E}[G(X)|U_r^\top X])^\top f(U_r^\top X)] = 0. $$

To prove (3.3), we again use the fundamental property of the conditional expectation to write

$$ 
\begin{aligned}
\mathbb{E}[\|G(X) - G^*(X)\|^2] &= \mathbb{E}[\|G(X) - V_s V_s^\top \mathbb{E}[G(X)|U_r^\top X] - (I_m - V_s V_s^\top) \mathbb{E}[G(X)]\|^2] \\
&= \mathbb{E}[\|(G(X) - \mathbb{E}[G(X)]) - V_s V_s^\top (\mathbb{E}[G(X)|U_r^\top X] - \mathbb{E}[G(X)])\|^2] \\
&= \mathbb{E}[\|G(X) - \mathbb{E}[G(X)]\|^2] - \mathbb{E}[\|V_s^\top (\mathbb{E}[G(X)|U_r^\top X] - \mathbb{E}[G(X)])\|^2] \\
&= \text{Tr}(\text{Cov}(G(X))) - \text{Tr}(V_s^\top \text{Cov}(\mathbb{E}[G(X)|U_r^\top X]) V_s),
\end{aligned}
$$

which concludes the proof.

**Appendix C. Proof of Theorem 3.5.** Let $V_{\perp} \in \mathbb{R}^{m \times (m-s)}$ be the orthogonal complement of $V_s$. We have

$$ 
\begin{aligned}
\mathbb{E}[\|G(X) - G^*(X)\|^2] &= \mathbb{E}[\|G(X) - V_s V_s^\top \mathbb{E}[G(X)|U_r^\top X] - V_\perp V_\perp^\top \mathbb{E}[G(X)]\|^2] \\
&= \mathbb{E}[\|V_\perp V_\perp^\top (G(X) - \mathbb{E}[G(X)]) - V_s V_s^\top (\mathbb{E}[G(X)|U_r^\top X] - G(X))\|^2] \\
(C.1) \quad &= \mathbb{E}[\|V_\perp V_\perp^\top (G(X) - \mathbb{E}[G(X)])\|^2] + \mathbb{E}[\|V_s V_s^\top (\mathbb{E}[G(X)|U_r^\top X] - G(X))\|^2]
\end{aligned}
$$
We show that the last term vanishes so that for any $\tilde{G}$ it holds $\mathbb{E}[\|G(X) - \tilde{G}(X)\|^2] \ge \mathbb{E}[\|G(X) - G^*(X)\|^2]$. This proves that $G^*$ is a minimizer of $\tilde{G} \mapsto \mathbb{E}[\|G(X) - \tilde{G}(X)\|^2]$. Let $V_{\perp} \in \mathbb{R}^{m \times (m-s)}$ be the orthogonal complement of $V_s$. Without loss of generality, we assume that $\tilde{G}_{\perp} \in \text{Ker}(V_s)$ so that $V_{\perp}V_{\perp}^{\top}\tilde{G}_{\perp} = \tilde{G}_{\perp}$. By construction, we have

$$ 
\begin{aligned}
G(X) - G^*(X) &= V_s(V_s^\top G(X) - g^*(U_r^\top X)) + V_\perp V_\perp^\top (G(X) - \mathbb{E}[G(X)]), \\
G^*(X) - \tilde{G}(X) &= V_s(g^*(U_r^\top X) - \tilde{g}(U_r^\top X)) + V_\perp V_\perp^\top (G_\perp^* - \tilde{G}_\perp).
\end{aligned}
$$

Thus, we get that

$$ 
\begin{aligned}
\mathbb{E}[(G(X) - G^*(X))^\top (G^*(X) - \tilde{G}(X))] \\
&= \mathbb{E}[(V_s^\top G(X) - g^*(U_r^\top X))^\top (g^*(U_r^\top X) - \tilde{g}(U_r^\top X))] \\
&\quad + \underbrace{\mathbb{E}[(G(X) - \mathbb{E}[G(X)])^\top V_\perp V_\perp^\top (G_\perp^* - \tilde{G}_\perp)]}_{=0},
\end{aligned}
$$

where the last term vanishes by taking the expectation of $G(X) - \mathbb{E}[G(X)]$. Next, we exploit the fundamental property of the conditional expectation $\mathbb{E}[G(X)f(U_r^\top X)] = \mathbb{E}[\mathbb{E}[G(X)|U_r^\top X]f(U_r^\top X)]$ which holds for any function $f: \mathbb{R}^r \to \mathbb{R}$. Considering the definition of $g^*$ and setting $f(U_r^\top X) = V_s(g^*(U_r^\top X) - \tilde{g}(U_r^\top X))$, we obtain

$$ \mathbb{E}[(G(X) - G^*(X))^\top (G^*(X) - \tilde{G}(X))] = \mathbb{E}[(G(X) - \mathbb{E}[G(X)|U_r^\top X])^\top f(U_r^\top X)] = 0. $$

To prove (3.3), we again use the fundamental property of the conditional expectation to write

$$ 
\begin{aligned}
\mathbb{E}[\|G(X) - G^*(X)\|^2] &= \mathbb{E}[\|G(X) - V_s V_s^\top \mathbb{E}[G(X)|U_r^\top X] - (I_m - V_s V_s^\top) \mathbb{E}[G(X)]\|^2] \\
&= \mathbb{E}[\|(G(X) - \mathbb{E}[G(X)]) - V_s V_s^\top (\mathbb{E}[G(X)|U_r^\top X] - \mathbb{E}[G(X)])\|^2] \\
&= \mathbb{E}[\|G(X) - \mathbb{E}[G(X)]\|^2] - \mathbb{E}[\|V_s^\top (\mathbb{E}[G(X)|U_r^\top X] - \mathbb{E}[G(X)])\|^2] \\
&= \text{Tr}(\text{Cov}(G(X))) - \text{Tr}(V_s^\top \text{Cov}(\mathbb{E}[G(X)|U_r^\top X]) V_s),
\end{aligned}
$$

which concludes the proof.

**Appendix C. Proof of Theorem 3.5.** Let $V_{\perp} \in \mathbb{R}^{m \times (m-s)}$ be the orthogonal complement of $V_s$. We have

$$ 
\begin{aligned}
\mathbb{E}[\|G(X) - G^*(X)\|^2] &= \mathbb{E}[\|G(X) - V_s V_s^\top \mathbb{E}[G(X)|U_r^\top X] - V_\perp V_\perp^\top \mathbb{E}[G(X)]\|^2] \\
&= \mathbb{E}[\|V_\perp V_\perp^\top (G(X) - \mathbb{E}[G(X)]) - V_s V_s^\top (\mathbb{E}[G(X)|U_r^\top X] - G(X))\|^2] \\
(C.1) \quad &= \mathbb{E}[\|V_\perp V_\perp^\top (G(X) - \mathbb{E}[G(X)])\|^2] + \mathbb{E}[\|V_s V_s^\top (\mathbb{E}[G(X)|U_r^\top X] - G(X))\|^2]
\end{aligned}
$$
Using the Poincaré inequality (3.5) for the component $f_i(X) = [V_⊥V_⊥^T G(X)]_i$ yields

$$
\begin{align*}
\mathbb{E} \left[ \| V_{\perp} V_{\perp}^{\top} (G(X) - \mathbb{E}[G(X)]) \|^2 \right] &= \sum_{i=1}^{m} \mathbb{E} \left[ (f_i(X) - \mathbb{E}[f_i(X)])^2 \right] \\
&\leq \sum_{i=1}^{m} \mathcal{C}(X) \mathbb{E} \left[ \| \nabla f_i(X) \|^2 \right] \\
&= \mathcal{C}(X) \mathbb{E} \left[ \|(I_m - V_s V_s^{\top}) \nabla G(X)\|_F^2 \right] \\
(C.2) \qquad &\qquad = \mathcal{C}(X) \left( \mathbb{E} \left[ \| \nabla G(X) \|_F^2 \right] - \mathbb{E} \left[ \| V_s^{\top} \nabla G(X) \|_F^2 \right] \right)
\end{align*}
$$

Moreover, the subspace Poincaré inequality (3.6) for $f_i(X) = [V_s V_s^\top G(X)]_i$ yields

$$
\begin{align*}
\mathbb{E} \left[ \| V_s V_s^\top (\mathbb{E}[G(X)|U_r^\top X] - G(X)) \|^2 \right] &= \sum_{i=1}^m \mathbb{E} \left[ \mathbb{E} \left[ (f_i(X) - \mathbb{E}[f_i(X)|U_r^\top X])^2 \right] \middle| U_r^\top X \right] \\
&\le \sum_{i=1}^m \mathbb{E} \left[ \mathcal{C}(X_{\perp}|X_r) \mathbb{E} \left[ \|(I_d - U_r U_r^\top) \nabla f_i(X)\|^2 \middle| U_r^\top X \right] \right] \\
&\le \sum_{i=1}^m \bar{\mathcal{C}}(X) \mathbb{E} \left[ \|(I_d - U_r U_r^\top) \nabla f_i(X)\|^2 \right] \\
&= \bar{\mathcal{C}}(X) \mathbb{E} \left[ \| V_s V_s^\top \nabla G(X) (I_d - U_r U_r^\top) \|_F^2 \right] \\
(C.3) \qquad &\qquad = \bar{\mathcal{C}}(X) \left( \mathbb{E} \left[ \| V_s^\top \nabla G(X) \|_F^2 \right] - \mathbb{E} \left[ \| V_s^\top \nabla G(X) U_r \|_F^2 \right] \right)
\end{align*}
$$

Finally, since $\mathcal{C}(X) \le \bar{\mathcal{C}}(X)$, injecting (C.2) and (C.3) in (C.1) yields the upper bound (3.8). The lower bound (3.9) is obtained in the same way by permuting the expectation and the norm in the above calculation. This concludes the proof.

**Appendix D. Cramér-Rao-like Inequality.**

LEMMA D.1. Let $X$ be a random variable on $\mathbb{R}^d$ with Lebesgue density $\pi_X$ such that $\text{supp}(\pi_X) = \mathbb{R}^d$. Let $f: \mathbb{R}^d \to \mathbb{R}$ be a differentiable function. Then it holds

$$
\begin{align*}
\text{(D.1)} \qquad & c(X) \|\mathbb{E}[\nabla f(X)]\|^2 \le \mathbb{E}[(f(X) - \mathbb{E}[f(X)])^2], \\
& \text{where } c(X) = \lambda_{\max}\{\mathbb{E}[\nabla \log \pi_X(X) \nabla \log \pi_X(X)^{\top}]\}^{-1}.
\end{align*}
$$

*Proof.* By partial integration, it holds for any differentiable function $h: \mathbb{R}^d \to \mathbb{R}$

$$
\begin{align*}
\mathbb{E}[\nabla h(X)] &= \int \nabla h(x) \pi_X(x) \, dx = - \int h(x) \nabla \pi_X(x) \, dx \\
(D.2) \qquad &\qquad = - \int h(x) \nabla \log \pi_X(x) \pi_X(x) \, dx = - \mathbb{E}[h(X) \nabla \log \pi_X(X)]
\end{align*}
$$

In the second to last equality we used the chain rule $\nabla \log h(x) = \frac{\nabla h(x)}{h(x)}$. Let $\alpha \in \mathbb{R}^d$ be such that $\|\alpha\| = 1$. Applying (D.2) to $h(x) = f(x) - \mathbb{E}[f(X)]$ and using the Cauchy-Schwarz inequality yields

$$
\begin{align*}
|\mathbb{E}[\nabla f(X)]^\top \alpha|^2 &= |\mathbb{E}[\nabla h(X)]^\top \alpha|^2 = \mathbb{E}[h(X) (\nabla \log \pi_X(X)^\top \alpha)]^2 \\
&\le \mathbb{E}[h(X)^2] \mathbb{E}[\nabla \log \pi_X(X)^\top \alpha]^2 \\
&= \operatorname{Var}(f(X)) \alpha^\top \mathbb{E}[\nabla \log \pi_X(X) \nabla \log \pi_X(X)^\top] \alpha \\
&\le \operatorname{Var}(f(X)) \lambda_{\max} (\mathbb{E}[\nabla \log \pi_X(X) \nabla \log \pi_X(X)^\top])
\end{align*}
$$
Using the Poincaré inequality (3.5) for the component $f_i(X) = [V_⊥V_⊥^T G(X)]_i$ yields

$$
\begin{align*}
\mathbb{E} \left[ \| V_{\perp} V_{\perp}^{\top} (G(X) - \mathbb{E}[G(X)]) \|^2 \right] &= \sum_{i=1}^{m} \mathbb{E} \left[ (f_i(X) - \mathbb{E}[f_i(X)])^2 \right] \\
&\le \sum_{i=1}^{m} \mathcal{C}(X) \mathbb{E} \left[ \| \nabla f_i(X) \|^2 \right] \\
&= \mathcal{C}(X) \mathbb{E} \left[ \|(I_m - V_s V_s^{\top}) \nabla G(X)\|_F^2 \right] \\
(C.2) \qquad \qquad \qquad \qquad \qquad \qquad \qquad &= \mathcal{C}(X) \left( \mathbb{E} \left[ \| \nabla G(X) \|_F^2 \right] - \mathbb{E} \left[ \| V_s^{\top} \nabla G(X) \|_F^2 \right] \right)
\end{align*}
$$

Moreover, the subspace Poincaré inequality (3.6) for $f_i(X) = [V_s V_s^\top G(X)]_i$ yields

$$
\begin{align*}
\mathbb{E} \left[ \| V_s V_s^\top (\mathbb{E}[G(X)|U_r^\top X] - G(X)) \|^2 \right] &= \sum_{i=1}^m \mathbb{E} \left[ \mathbb{E} \left[ (f_i(X) - \mathbb{E}[f_i(X)|U_r^\top X])^2 \right] \middle| U_r^\top X \right] \\
&\le \sum_{i=1}^m \mathbb{E} \left[ \mathcal{C}(X_{\perp}|X_r) \mathbb{E} \left[ \|(I_d - U_r U_r^\top) \nabla f_i(X)\|^2 \middle| U_r^\top X \right] \right] \\
&\le \sum_{i=1}^m \bar{\mathcal{C}}(X) \mathbb{E} \left[ \|(I_d - U_r U_r^\top) \nabla f_i(X)\|^2 \right] \\
&= \bar{\mathcal{C}}(X) \mathbb{E} \left[ \| V_s V_s^\top \nabla G(X) (I_d - U_r U_r^\top) \|_F^2 \right] \\
(C.3) \qquad \qquad \qquad \qquad \qquad \qquad \qquad &= \bar{\mathcal{C}}(X) \left( \mathbb{E} \left[ \| V_s^\top \nabla G(X) \|_F^2 \right] - \mathbb{E} \left[ \| V_s^\top \nabla G(X) U_r \|_F^2 \right] \right)
\end{align*}
$$

Finally, since $\mathcal{C}(X) \le \bar{\mathcal{C}}(X)$, injecting (C.2) and (C.3) in (C.1) yields the upper bound (3.8). The lower bound (3.9) is obtained in the same way by permuting the expectation and the norm in the above calculation. This concludes the proof.

**Appendix D. Cramér-Rao-like Inequality.**

LEMMA D.1. Let $X$ be a random variable on $\mathbb{R}^d$ with Lebesgue density $\pi_X$ such that $\text{supp}(\pi_X) = \mathbb{R}^d$. Let $f: \mathbb{R}^d \to \mathbb{R}$ be a differentiable function. Then it holds

$$
\begin{align*}
\text{(D.1)} \qquad & c(X) \|\mathbb{E}[\nabla f(X)]\|^2 \le \mathbb{E}[(f(X) - \mathbb{E}[f(X)])^2], \\
& \text{where } c(X) = \lambda_{\max}\{\mathbb{E}[\nabla \log \pi_X(X) \nabla \log \pi_X(X)^{\top}]\}^{-1}.
\end{align*}
$$

*Proof.* By partial integration, it holds for any differentiable function $h: \mathbb{R}^d \to \mathbb{R}$

$$
\begin{align*}
\mathbb{E}[\nabla h(X)] &= \int \nabla h(x) \pi_X(x) \, dx = - \int h(x) \nabla \pi_X(x) \, dx \\
\text{(D.2)} \qquad \qquad \qquad \qquad \qquad \qquad \qquad &= - \int h(x) \nabla \log \pi_X(x) \pi_X(x) \, dx = - \mathbb{E}[h(X) \nabla \log \pi_X(X)]
\end{align*}
$$

In the second to last equality we used the chain rule $\nabla \log h(x) = \frac{\nabla h(x)}{h(x)}$. Let $\alpha \in \mathbb{R}^d$ be such that $\|\alpha\| = 1$. Applying (D.2) to $h(x) = f(x) - \mathbb{E}[f(X)]$ and using the Cauchy-Schwarz inequality yields

$$
\begin{align*}
|\mathbb{E}[\nabla f(X)]^\top \alpha|^2 &= |\mathbb{E}[\nabla h(X)]^\top \alpha|^2 = \mathbb{E} \left[ h(X) (\nabla \log \pi_X(X)^\top \alpha) \right]^2 \\
&\le \mathbb{E} [h(X)^2] \mathbb{E} \left[ \nabla \log \pi_X(X)^\top \alpha \right]^2 \\
&= \operatorname{Var}(f(X)) \alpha^\top \mathbb{E} \left[ \nabla \log \pi_X(X) \nabla \log \pi_X(X)^\top \alpha \right] \\
&\le \operatorname{Var}(f(X)) \lambda_{\max} \left( \mathbb{E} \left[ \nabla \log \pi_X(X) \nabla \log \pi_X(X)^\top \right] \right)
\end{align*}
$$
Thus, we obtain inequality (D.1). □

We call (D.1) the Cramér-Rao-like inequality since it poses a similar lower bound to the variance involving the inverse Fisher information matrix in $c(X)$ [44]. The following Lemma establishes sufficient conditions for every conditional distribution of $X_⊥$ to satisfy the Cramér-Rao inequality, i.e., $\bar{c}(X) \ge 0$.

LEMMA D.2 (see Proposition 3.6). Let $X$ be a random variable on $\mathbb{R}^d$ with Lebesgue density $\pi_X$ such that $\text{supp}(\pi_X) = \mathbb{R}^d$ and there exists $\eta < \infty$ so that

$$ (D.3) \qquad -\nabla^2 \log \pi_X(x) \le \eta I_d. $$

Then there exists a positive constant $\bar{c}(X) \ge 1/\eta$ such that

$$ (D.4) \qquad \bar{c}(X) \|(I_d - U_r U_r^\top) \mathbb{E}[\nabla f(X)]\|^2 \le \mathbb{E}[(f(X) - \mathbb{E}[f(X)|U_r^\top X])^2], $$

holds for any differentiable function $f: \mathbb{R}^d \to \mathbb{R}$ with $\mathbb{E}[f(X)^2] < \infty$ and $\mathbb{E}[\|\nabla f(X)\|^2] < \infty$ and for any matrix $U_r \in \mathbb{R}^{d \times r}$ with orthogonal columns.

Proof. Let $U_⊥ \in \mathbb{R}^{d \times (d-r)}$ be orthogonal to $U_r$ so that $X = U_r X_r + U_⊥ X_⊥$. Further, let $\nabla_⊥$ denote the derivative with respect to the argument $X_⊥$. We write $\mathbb{E}_⊥$ and $\text{Var}_⊥$ for the expectation and variance with respect to the conditional density $\pi_{X_⊥|X_r}$. From the proof of Lemma D.1, we have for any differentiable function $h: \mathbb{R}^d \to \mathbb{R}$

$$ \mathbb{E}_⊥[\nabla_⊥h(X_⊥, x_r)] = -\mathbb{E}_⊥[h(X_⊥, x_r) \nabla_⊥ \log \pi_{X_⊥|X_r}(X_⊥|x_r)] $$

Let $\alpha \in \mathbb{R}^d$ such that $\|\alpha\| = 1$. Similar to the proof of Lemma D.1, we set $h(x_⊥, x_r) = f(x_⊥, x_r) - \mathbb{E}_⊥[f(X_⊥, x_r)]$ to be the centered version of $f$ and obtain

$$ (D.5) \qquad \begin{aligned} |\mathbb{E}_⊥[\nabla_⊥ f(X_⊥, x_r)]^\top \alpha|^2 &= |\mathbb{E}_⊥[\nabla_⊥ h(X_⊥, x_r)]^\top \alpha|^2 \\ &\le \mathbb{E}_⊥[h(X_⊥, x_r)^2] \mathbb{E}_⊥\left[\left(\nabla_⊥ \log \pi_{X_⊥|X_r}(X_⊥|x_r)^\top \alpha\right)^2\right] \\ &= \text{Var}_⊥(f(X_⊥, x_r)) \alpha^\top \mathcal{I}_⊥(x_r) \alpha, \end{aligned} $$

where

$$ \begin{aligned} \mathcal{I}_⊥(x_r) &= \int \nabla_⊥ \log \pi_{X_⊥|X_r}(x_⊥|x_r) \nabla_⊥ \log \pi_{X_⊥|X_r}(x_⊥|x_r)^\top \mathrm{d}\pi_{X_⊥|X_r} \\ &= \int \nabla_⊥ \pi_{X_⊥|X_r}(x_⊥|x_r) \nabla_⊥ \log \pi_{X_⊥|X_r}(x_⊥|x_r)^\top \mathrm{d}x_⊥ \\ &= - \int \nabla_⊥^2 \log \pi_X(x_⊥|x_r) \pi_{X_⊥|X_r}(x_⊥|x_r) \mathrm{d}x_⊥ \\ &= -U_⊥^\top \int \nabla_⊥^2 \log \pi_X(x) \mathrm{d}\pi_{X_⊥|X_r} U_⊥. \end{aligned} $$

The second line uses the chain rule for $\nabla \log$ while the third uses partial integration. In the last line, we apply the relationship $\log \pi_{X_⊥|X_r} = \log \pi_X - \log \pi_{X_r}$ as well as the chain rule for $\nabla_⊥$. Due to the assumption $-\nabla_⊥^2 \log \pi_X(x) \le \eta I_d$, we have

$$ \mathcal{I}_⊥(x_r) \le \eta I_{d-r} $$

for any $U_r$. Thus, Equation (D.5) yields

$$ |\mathbb{E}_⊥[\nabla_⊥ f(X_⊥, x_r)]^\top \alpha|^2 \le \eta \text{Var}_⊥(f(X_⊥, x_r)) $$
Thus, we obtain inequality (D.1). □

We call (D.1) the Cramér-Rao-like inequality since it poses a similar lower bound to the variance involving the inverse Fisher information matrix in $c(X)$ [44]. The following Lemma establishes sufficient conditions for every conditional distribution of $X_⊥$ to satisfy the Cramér-Rao inequality, i.e., $\bar{c}(X) \ge 0$.

LEMMA D.2 (see Proposition 3.6). Let $X$ be a random variable on $\mathbb{R}^d$ with Lebesgue density $\pi_X$ such that $\text{supp}(\pi_X) = \mathbb{R}^d$ and there exists $\eta < \infty$ so that

$$ (D.3) \qquad -\nabla^2 \log \pi_X(x) \le \eta I_d. $$

Then there exists a positive constant $\bar{c}(X) \ge 1/\eta$ such that

$$ (D.4) \qquad \bar{c}(X) \|(I_d - U_r U_r^\top) \mathbb{E}[\nabla f(X)]\|^2 \le \mathbb{E}[(f(X) - \mathbb{E}[f(X)|U_r^\top X])^2], $$

holds for any differentiable function $f: \mathbb{R}^d \to \mathbb{R}$ with $\mathbb{E}[f(X)^2] < \infty$ and $\mathbb{E}[\|\nabla f(X)\|^2] < \infty$ and for any matrix $U_r \in \mathbb{R}^{d \times r}$ with orthogonal columns.

Proof. Let $U_⊥ \in \mathbb{R}^{d \times (d-r)}$ be orthogonal to $U_r$ so that $X = U_r X_r + U_⊥ X_⊥$. Further, let $\nabla_⊥$ denote the derivative with respect to the argument $X_⊥$. We write $\mathbb{E}_⊥$ and $\text{Var}_⊥$ for the expectation and variance with respect to the conditional density $\pi_{X_⊥|X_r}$. From the proof of Lemma D.1, we have for any differentiable function $h: \mathbb{R}^d \to \mathbb{R}$

$$ \mathbb{E}_⊥[\nabla_⊥ h(X_⊥, x_r)] = -\mathbb{E}_⊥[h(X_⊥, x_r) \nabla_⊥ \log \pi_{X_⊥|X_r}(X_⊥|x_r)] $$

Let $\alpha \in \mathbb{R}^d$ such that $\|\alpha\| = 1$. Similar to the proof of Lemma D.1, we set $h(x_⊥, x_r) = f(x_⊥, x_r) - \mathbb{E}_⊥[f(X_⊥, x_r)]$ to be the centered version of $f$ and obtain

$$ (D.5) \qquad \begin{aligned} |\mathbb{E}_⊥[\nabla_⊥ f(X_⊥, x_r)]^\top \alpha|^2 &= |\mathbb{E}_⊥[\nabla_⊥ h(X_⊥, x_r)]^\top \alpha|^2 \\ &\le \mathbb{E}_⊥[h(X_⊥, x_r)^2] \mathbb{E}_⊥\left[\left(\nabla_⊥ \log \pi_{X_⊥|X_r}(X_⊥|x_r)^\top \alpha\right)^2\right] \\ &= \text{Var}_⊥(f(X_⊥, x_r)) \alpha^\top \mathcal{I}_⊥(x_r) \alpha, \end{aligned} $$

where

$$ \begin{aligned} \mathcal{I}_⊥(x_r) &= \int \nabla_⊥ \log \pi_{X_⊥|X_r}(x_⊥|x_r) \nabla_⊥ \log \pi_{X_⊥|X_r}(x_⊥|x_r)^\top \mathrm{d}\pi_{X_⊥|X_r} \\ &= \int \nabla_⊥ \pi_{X_⊥|X_r}(x_⊥|x_r) \nabla_⊥ \log \pi_{X_⊥|X_r}(x_⊥|x_r)^\top \mathrm{d}x_⊥ \\ &= - \int \nabla_⊥^2 \log \pi_X(x_⊥|x_r) \pi_{X_⊥|X_r}(x_⊥|x_r) \mathrm{d}x_⊥ \\ &= -U_⊥^\top \int \nabla_⊥^2 \log \pi_X(x) \mathrm{d}\pi_{X_⊥|X_r} U_⊥. \end{aligned} $$

The second line uses the chain rule for $\nabla \log$ while the third uses partial integration. In the last line, we apply the relationship $\log \pi_{X_⊥|X_r} = \log \pi_X - \log \pi_{X_r}$ as well as the chain rule for $\nabla_⊥$. Due to the assumption $-\nabla_⊥^2 \log \pi_X(x) \le \eta I_d$, we have

$$ \mathcal{I}_⊥(x_r) \le \eta I_{d-r} $$

for any $U_r$. Thus, Equation (D.5) yields

$$ |\mathbb{E}_⊥[\nabla_⊥ f(X_⊥, x_r)]^\top \alpha|^2 \le \eta \text{Var}_⊥(f(X_⊥, x_r)) $$
so that taking the expectation over the marginal $\pi_{X_r}$ and applying Jensen's inequality on the left hand side, we get

$$ \left| \alpha^{\top} U_{\perp}^{\top} \int \nabla f(x) d\pi_X \right|^2 \leq \eta \mathbb{E} \left[ \left( f(X) - \mathbb{E}[f(X) | U_r^{\top} X] \right)^2 \right] $$

Taking the supremum over $\|\alpha\| = 1$ finally yields

$$ \frac{1}{\eta} \left\| (I_d - U_r U_r^{\top}) \mathbb{E}[\nabla f(X)] \right\|^2 \leq \mathbb{E} \left[ (f(X) - \mathbb{E}[f(X) | U_r^{\top} X])^2 \right] $$

for any $U_r$.

□

Remark D.3. Assumption (D.3) is equivalent to saying that the function $h(x) := \frac{\eta}{2}\|x\|^2 + \log \pi_X(x)$ is convex. This implies that $h(x) \geq h(y) + \nabla h(y)^{\top}(x - y)$ for all $x, y \in \mathbb{R}^d$ and, by letting $y = 0$, we deduce that

$$ \pi_X(x) \geq C \exp \left( -\frac{\eta}{2} \|x - \eta^{-1} \nabla \log \pi_X(0)\|^2 \right), $$

holds for all $x \in \mathbb{R}^d$, where $C > 0$ is a constant independent of $x$. This means that if (3.12) holds, then the density $\pi_X$ is necessarily above some rescaled Gaussian density.

## Appendix E. Proof of Proposition 4.1.

On one side, the law of total variance enables us to write

$$ \begin{aligned} S^{\text{cl}}(U_r|V_s) &= 1 - \frac{\mathbb{E}[\|V_s^{\top} G(X) - \mathbb{E}[V_s^{\top} G(X) | U_r^{\top} X]\|^2]}{\text{Tr}(\text{Cov}(V_s^{\top} G(X)))} \\ S^{\text{tot}}(U_r|V_s) &= \frac{\mathbb{E}[\|V_s^{\top} G(X) - \mathbb{E}[V_s^{\top} G(X) | U_{\perp}^{\top} X]\|^2]}{\text{Tr}(\text{Cov}(V_s^{\top} G(X)))}. \end{aligned} $$

On the other side, applying Theorem 3.5 with $V_s^{\top} G(X)$ instead of $G(X)$ yields

$$ \begin{aligned} & \mathbb{E}[\|V_s^{\top} G(X) - \mathbb{E}[V_s^{\top} G(X) | U_r^{\top} X]\|^2] \\ & \leq \bar{c}(X) (\mathbb{E}[\|V_s^{\top} \nabla G(X)\|_F^2] - \mathbb{E}[\|V_s^{\top} \nabla G(X) U_r\|_F^2]) \\ & = \bar{c}(X) \text{Tr}(U_{\perp}^{\top} H_X(V_s) U_{\perp}) \\ & \mathbb{E}[\|V_s^{\top} G(X) - \mathbb{E}[V_s^{\top} G(X) | U_r^{\top} X]\|^2] \\ & \geq \bar{c}(X) (\|V_s^{\top} \mathbb{E}[\nabla G(X)]\|_F^2 - \|V_s^{\top} \mathbb{E}[\nabla G(X)] U_r\|_F^2) \\ & = \bar{c}(X) \|V_s^{\top} \mathbb{E}[\nabla G(X)] U_{\perp}\|_F^2 \end{aligned} $$

By combining these relations we directly obtain (4.8) and, by replacing $U_r$ with $U_{\perp}$ in the above inequalities, we obtain (4.9). This concludes the proof.

## REFERENCES

[1] P.-A. ABSIL, R. MAHONY, AND R. SEPULCHRE, *Optimization algorithms on matrix manifolds*, in Optimization Algorithms on Matrix Manifolds, Princeton University Press, 2009.

[2] A. ALEXANDERIAN, *Optimal experimental design for infinite-dimensional bayesian inverse problems governed by pdes: A review*, Inverse Problems, 37 (2021), p. 043001.

[3] A. ALEXANDERIAN, N. PETRA, G. STADLER, AND O. GHATTAS, *A fast and scalable method for a-optimal design of experiments for infinite-dimensional bayesian nonlinear inverse problems*, SIAM Journal on Scientific Computing, 38 (2016), pp. A243–A272.
so that taking the expectation over the marginal $\pi_{X_r}$ and applying Jensen's inequality on the left hand side, we get

$$ \left| \alpha^{\top} U_{\perp}^{\top} \int \nabla f(x) d\pi_X \right|^2 \leq \eta \mathbb{E} \left[ \left( f(X) - \mathbb{E}[f(X) | U_r^{\top} X] \right)^2 \right] $$

Taking the supremum over $\|\alpha\| = 1$ finally yields

$$ \frac{1}{\eta} \left\| (I_d - U_r U_r^{\top}) \mathbb{E}[\nabla f(X)] \right\|^2 \leq \mathbb{E} \left[ (f(X) - \mathbb{E}[f(X) | U_r^{\top} X])^2 \right] $$

for any $U_r$.

□

Remark D.3. Assumption (D.3) is equivalent to saying that the function $h(x) := \frac{\eta}{2}\|x\|^2 + \log \pi_X(x)$ is convex. This implies that $h(x) \geq h(y) + \nabla h(y)^{\top}(x - y)$ for all $x, y \in \mathbb{R}^d$ and, by letting $y = 0$, we deduce that

$$ \pi_X(x) \geq C \exp \left( -\frac{\eta}{2} \|x - \eta^{-1} \nabla \log \pi_X(0)\|^2 \right), $$

holds for all $x \in \mathbb{R}^d$, where $C > 0$ is a constant independent of $x$. This means that if (3.12) holds, then the density $\pi_X$ is necessarily above some rescaled Gaussian density.

## Appendix E. Proof of Proposition 4.1.

On one side, the law of total variance enables us to write

$$ \begin{aligned} S^{\text{cl}}(U_r|V_s) &= 1 - \frac{\mathbb{E}[\|V_s^{\top} G(X) - \mathbb{E}[V_s^{\top} G(X)|U_r^{\top} X]\|^2]}{\text{Tr}(\text{Cov}(V_s^{\top} G(X)))} \\ S^{\text{tot}}(U_r|V_s) &= \frac{\mathbb{E}[\|V_s^{\top} G(X) - \mathbb{E}[V_s^{\top} G(X)|U_{\perp}^{\top} X]\|^2]}{\text{Tr}(\text{Cov}(V_s^{\top} G(X)))}. \end{aligned} $$

On the other side, applying Theorem 3.5 with $V_s^{\top} G(X)$ instead of $G(X)$ yields

$$ \begin{aligned} & \mathbb{E}[\|V_s^{\top} G(X) - \mathbb{E}[V_s^{\top} G(X)|U_r^{\top} X]\|^2] \\ & \leq \bar{c}(X) (\mathbb{E}[\|V_s^{\top} \nabla G(X)\|_F^2] - \mathbb{E}[\|V_s^{\top} \nabla G(X) U_r\|_F^2]) \\ & = \bar{c}(X) \text{Tr}(U_{\perp}^{\top} H_X(V_s) U_{\perp}) \\ & \mathbb{E}[\|V_s^{\top} G(X) - \mathbb{E}[V_s^{\top} G(X)|U_r^{\top} X]\|^2] \\ & \geq \bar{c}(X) (\|V_s^{\top} \mathbb{E}[\nabla G(X)]\|_F^2 - \|V_s^{\top} \mathbb{E}[\nabla G(X)] U_r\|_F^2) \\ & = \bar{c}(X) \|V_s^{\top} \mathbb{E}[\nabla G(X)] U_{\perp}\|_F^2 \end{aligned} $$

By combining these relations we directly obtain (4.8) and, by replacing $U_r$ with $U_{\perp}$ in the above inequalities, we obtain (4.9). This concludes the proof.

## REFERENCES

[1] P.-A. ABSIL, R. MAHONY, AND R. SEPULCHRE, *Optimization algorithms on matrix manifolds*, in Optimization Algorithms on Matrix Manifolds, Princeton University Press, 2009.

[2] A. ALEXANDERIAN, *Optimal experimental design for infinite-dimensional bayesian inverse problems governed by pdes: A review*, Inverse Problems, 37 (2021), p. 043001.

[3] A. ALEXANDERIAN, N. PETRA, G. STADLER, AND O. GHATTAS, *A fast and scalable method for a-optimal design of experiments for infinite-dimensional bayesian nonlinear inverse problems*, SIAM Journal on Scientific Computing, 38 (2016), pp. A243–A272.
[4] A. ALEXANDERIAN, W. REESE, R. C. SMITH, AND M. YU, Model input and output dimension reduction using karhunen-loève expansions with application to biotransport, ASCE-ASME J Risk and Uncert in Engrg Sys Part B Mech Engrg, 5 (2019).

[5] B. O. ALMROTH, P. STERN, AND F. A. BROGAN, Automatic choice of global shape functions in structural analysis, Aiaa Journal, 16 (1978), pp. 525–528.

[6] C. ANÈ, S. BLACHÈRE, D. CHAFAÏ, P. FOUGÈRES, I. GENTIL, F. MALRIEU, C. ROBERTO, AND G. SCHEFFER, *Sur les inégalités de Sobolev logarithmiques*, vol. 10, Société mathématique de France Paris, 2000.

[7] N. ARETZ, P. CHEN, D. DEGEN, AND K. VEROY, *A greedy sensor selection algorithm for hyper-parameterized linear bayesian inverse problems*, arXiv preprint arXiv:2301.12019, (2023).

[8] J.-P. ARGAUD, B. BOURIQUET, F. DE CASO, H. GONG, Y. MADAY, AND O. MULA, *Sensor placement in nuclear reactors based on the generalized empirical interpolation method*, Journal of Computational Physics, 363 (2018), pp. 354–370.

[9] A. ATTIA, A. ALEXANDERIAN, AND A. K. SAIBABA, *Goal-oriented optimal design of experiments for large-scale bayesian linear inverse problems*, Inverse Problems, 34 (2018), p. 095009.

[10] D. BAKRY, I. GENTIL, M. LEDOUX, ET AL., *Analysis and geometry of Markov diffusion operators*, vol. 103, Springer, 2014.

[11] R. BAPTISTA, L. CAO, J. CHEN, O. GHATTAS, F. LI, Y. M. MARZOUK, AND J. T. ODEN, *Bayesian model calibration for block copolymer self-assembly: Likelihood-free inference and expected information gain computation via measure transport*, Journal of Computational Physics, (2024), p. 112844.

[12] R. BAPTISTA, Y. MARZOUK, AND O. ZAHM, *Gradient-based data and parameter dimension reduction for bayesian models: an information theoretic perspective*, arXiv preprint arXiv:2207.08670, (2022).

[13] M. BARRAULT, Y. MADAY, N. C. NGUYEN, AND A. T. PATERA, *An 'empirical interpolation' method: application to efficient reduced-basis discretization of partial differential equations*, Comptes Rendus Mathematique, 339 (2004), pp. 667–672.

[14] J. BECK, B. M. DIA, L. F. ESPATH, Q. LONG, AND R. TEMPONE, *Fast bayesian experimental design: Laplace-based importance sampling for the expected information gain*, Computer Methods in Applied Mechanics and Engineering, 334 (2018), pp. 523–553.

[15] K. BHATTACHARYA, B. HOSSEINI, N. B. KOVACHKI, AND A. M. STUART, *Model reduction and neural networks for parametric pdes*, The SMAI journal of computational mathematics, 7 (2021), pp. 121–157.

[16] P. BINEV, A. COHEN, O. MULA, AND J. NICHOLS, *Greedy algorithms for optimal measurements selection in state estimation using reduced models*, SIAM/ASA Journal on Uncertainty Quantification, 6 (2018), pp. 1101–1126.

[17] P. G. CONSTANTINE, E. DOW, AND Q. WANG, *Active subspace methods in theory and practice: applications to kriging surfaces*, SIAM Journal on Scientific Computing, 36 (2014), pp. A1500–A1524.

[18] T. CUI, J. MARTIN, Y. M. MARZOUK, A. SOLONEN, AND A. SPANTINI, *Likelihood-informed dimension reduction for nonlinear inverse problems*, Inverse Problems, 30 (2014), p. 114015.

[19] T. CUI, Y. MARZOUK, AND K. WILLCOX, *Scalable posterior approximations for large-scale bayesian inverse problems via likelihood-informed parameter and state reduction*, Journal of Computational Physics, 315 (2016), pp. 363–387.

[20] T. CUI AND O. ZAHM, *Data-free likelihood-informed dimension reduction of bayesian inverse problems*, Inverse Problems, 37 (2021), p. 045009.

[21] C. FENG AND Y. M. MARZOUK, *A layered multiple importance sampling scheme for focused optimal bayesian experimental design*, arXiv preprint arXiv:1903.11187, (2019).

[22] A. FOSTER, M. JANKOWIAK, E. BINGHAM, P. HORSFALL, Y. W. TEH, T. RAINFORTH, AND N. GOODMAN, *Variational bayesian optimal experimental design*, Advances in Neural Information Processing Systems, 32 (2019).

[23] F. GAMBOA, P. GREMAUD, T. KLEIN, AND A. LAGNOUX, *Global sensitivity analysis: A novel generation of mighty estimators based on rank statistics*, Bernoulli, 28 (2022).

[24] F. GAMBOA, A. JANON, T. KLEIN, AND A. LAGNOUX, *Sensitivity analysis for multidimensional and functional outputs*, (2014).

[25] L. GIRALDI, O. P. LE MAÎTRE, I. HOTEIT, AND O. M. KNIO, *Optimal projection of observations in a bayesian setting*, Computational Statistics & Data Analysis, 124 (2018), pp. 252–276.

[26] E. HABER, L. HORESH, AND L. TENORIO, *Numerical methods for experimental design of large-scale linear ill-posed inverse problems*, Inverse Problems, 24 (2008), p. 055012.

[27] R. A. HORN AND C. R. JOHNSON, *Matrix analysis*, Cambridge university press, 2012.

[28] H. HOTELLING, *Analysis of a complex of statistical variables into principal components.*, Journal of educational psychology, 24 (1933), p. 417.
[4] A. ALEXANDERIAN, W. REESE, R. C. SMITH, AND M. YU, Model input and output dimension reduction using karhunen-loève expansions with application to biotransport, ASCE-ASME J Risk and Uncert in Engrg Sys Part B Mech Engrg, 5 (2019).

[5] B. O. ALMROTH, P. STERN, AND F. A. BROGAN, Automatic choice of global shape functions in structural analysis, Aiaa Journal, 16 (1978), pp. 525–528.

[6] C. ANÈ, S. BLACHÈRE, D. CHAFAÏ, P. FOUGÈRES, I. GENTIL, F. MALRIEU, C. ROBERTO, AND G. SCHEFFER, *Sur les inégalités de Sobolev logarithmiques*, vol. 10, Société mathématique de France Paris, 2000.

[7] N. ARETZ, P. CHEN, D. DEGEN, AND K. VEROY, *A greedy sensor selection algorithm for hyper-parameterized linear bayesian inverse problems*, arXiv preprint arXiv:2301.12019, (2023).

[8] J.-P. ARGAUD, B. BOURIQUET, F. DE CASO, H. GONG, Y. MADAY, AND O. MULA, *Sensor placement in nuclear reactors based on the generalized empirical interpolation method*, Journal of Computational Physics, 363 (2018), pp. 354–370.

[9] A. ATTIA, A. ALEXANDERIAN, AND A. K. SAIBABA, *Goal-oriented optimal design of experiments for large-scale bayesian linear inverse problems*, Inverse Problems, 34 (2018), p. 095009.

[10] D. BAKRY, I. GENTIL, M. LEDOUX, ET AL., *Analysis and geometry of Markov diffusion operators*, vol. 103, Springer, 2014.

[11] R. BAPTISTA, L. CAO, J. CHEN, O. GHATTAS, F. LI, Y. M. MARZOUK, AND J. T. ODEN, *Bayesian model calibration for block copolymer self-assembly: Likelihood-free inference and expected information gain computation via measure transport*, Journal of Computational Physics, (2024), p. 112844.

[12] R. BAPTISTA, Y. MARZOUK, AND O. ZAHM, *Gradient-based data and parameter dimension reduction for bayesian models: an information theoretic perspective*, arXiv preprint arXiv:2207.08670, (2022).

[13] M. BARRAULT, Y. MADAY, N. C. NGUYEN, AND A. T. PATERA, *An 'empirical interpolation' method: application to efficient reduced-basis discretization of partial differential equations*, Comptes Rendus Mathematique, 339 (2004), pp. 667–672.

[14] J. BECK, B. M. DIA, L. F. ESPATH, Q. LONG, AND R. TEMPONE, *Fast bayesian experimental design: Laplace-based importance sampling for the expected information gain*, Computer Methods in Applied Mechanics and Engineering, 334 (2018), pp. 523–553.

[15] K. BHATTACHARYA, B. HOSSEINI, N. B. KOVACHKI, AND A. M. STUART, *Model reduction and neural networks for parametric pdes*, The SMAI journal of computational mathematics, 7 (2021), pp. 121–157.

[16] P. BINEV, A. COHEN, O. MULA, AND J. NICHOLS, *Greedy algorithms for optimal measurements selection in state estimation using reduced models*, SIAM/ASA Journal on Uncertainty Quantification, 6 (2018), pp. 1101–1126.

[17] P. G. CONSTANTINE, E. DOW, AND Q. WANG, *Active subspace methods in theory and practice: applications to kriging surfaces*, SIAM Journal on Scientific Computing, 36 (2014), pp. A1500–A1524.

[18] T. CUI, J. MARTIN, Y. M. MARZOUK, A. SOLONEN, AND A. SPANTINI, *Likelihood-informed dimension reduction for nonlinear inverse problems*, Inverse Problems, 30 (2014), p. 114015.

[19] T. CUI, Y. MARZOUK, AND K. WILLCOX, *Scalable posterior approximations for large-scale bayesian inverse problems via likelihood-informed parameter and state reduction*, Journal of Computational Physics, 315 (2016), pp. 363–387.

[20] T. CUI AND O. ZAHM, *Data-free likelihood-informed dimension reduction of bayesian inverse problems*, Inverse Problems, 37 (2021), p. 045009.

[21] C. FENG AND Y. M. MARZOUK, *A layered multiple importance sampling scheme for focused optimal bayesian experimental design*, arXiv preprint arXiv:1903.11187, (2019).

[22] A. FOSTER, M. JANKOWIAK, E. BINGHAM, P. HORSFALL, Y. W. TEH, T. RAINFORTH, AND N. GOODMAN, *Variational bayesian optimal experimental design*, Advances in Neural Information Processing Systems, 32 (2019).

[23] F. GAMBOA, P. GREMAUD, T. KLEIN, AND A. LAGNOUX, *Global sensitivity analysis: A novel generation of mighty estimators based on rank statistics*, Bernoulli, 28 (2022).

[24] F. GAMBOA, A. JANON, T. KLEIN, AND A. LAGNOUX, *Sensitivity analysis for multidimensional and functional outputs*, (2014).

[25] L. GIRALDI, O. P. LE MAÎTRE, I. HOTEIT, AND O. M. KNIO, *Optimal projection of observations in a bayesian setting*, Computational Statistics & Data Analysis, 124 (2018), pp. 252–276.

[26] E. HABER, L. HORESH, AND L. TENORIO, *Numerical methods for experimental design of large-scale linear ill-posed inverse problems*, Inverse Problems, 24 (2008), p. 055012.

[27] R. A. HORN AND C. R. JOHNSON, *Matrix analysis*, Cambridge university press, 2012.

[28] H. HOTELLING, *Analysis of a complex of statistical variables into principal components.*, Journal of educational psychology, 24 (1933), p. 417.
[29] H. HOTELLING, *Relations between two sets of variates*, in *Breakthroughs in statistics: methodology and distribution*, Springer, 1992, pp. 162–190.

[30] X. HUAN AND Y. M. MARZOUK, *Simulation-based optimal bayesian experimental design for nonlinear systems*, Journal of Computational Physics, 232 (2013), pp. 288–317.

[31] N. B. KOVACHKI, S. LANTHALER, AND A. M. STUART, *Operator learning: Algorithms and analysis*, arXiv preprint arXiv:2402.15715, (2024).

[32] K. KOVAL, R. HERZOG, AND R. SCHEICHL, *Tractable optimal experimental design using transport maps*, arXiv preprint arXiv:2401.07971, (2024).

[33] A. KRASKOV, H. STÖGBAUER, AND P. GRASSBERGER, *Estimating mutual information*, Physical review E, 69 (2004), p. 066138.

[34] A. KRAUSE, A. SINGH, AND C. GUESTRIN, *Near-optimal sensor placements in gaussian processes: Theory, efficient algorithms and empirical studies.*, Journal of Machine Learning Research, 9 (2008).

[35] S. KUCHERENKO AND B. IOOSS, *Derivative based global sensitivity measures*, arXiv preprint arXiv:1412.2619, (2014).

[36] C. LIEBERMAN, K. WILLCOX, AND O. GHATTAS, *Parameter and state model reduction for large-scale statistical inverse problems*, SIAM Journal on Scientific Computing, 32 (2010), pp. 2523–2542.

[37] M. LOÈVE AND M. LOÈVE, *Elementary probability theory*, Springer, 1977.

[38] L. LU, P. JIN, G. PANG, Z. ZHANG, AND G. E. KARNIADAKIS, *Learning nonlinear operators via deeponet based on the universal approximation theorem of operators*, Nature machine intelligence, 3 (2021), pp. 218–229.

[39] J. L. LUMLEY, *The structure of inhomogeneous turbulent flows*, Atmospheric turbulence and radio wave propagation, (1967), pp. 166–178.

[40] Y. MADAY, N. C. NGUYEN, A. T. PATERA, AND G. S. PAU, *A general, multipurpose interpolation procedure: the magic points*, (2007).

[41] T. O'LEARY-ROSEBERRY, P. CHEN, U. VILLA, AND O. GHATTAS, *Derivative-informed neural operator: an efficient framework for high-dimensional parametric derivative learning*, Journal of Computational Physics, 496 (2024), p. 112555.

[42] H. POINCARÉ, *Sur les équations aux dérivées partielles de la physique mathématique*, American Journal of Mathematics, (1890), pp. 211–294.

[43] E. QIAN, B. PEHERSTORFER, D. O'MALLEY, V. V. VESSELINOV, AND K. WILLCOX, *Multifidelity monte carlo estimation of variance and sensitivity indices*, SIAM/ASA Journal on Uncertainty Quantification, 6 (2018), pp. 683–706.

[44] C. RADHAKRISHNA RAO, *Information and accuracy attainable in the estimation of statistical parameters*, Bulletin of the Calcutta Mathematical Society, 37 (1945), pp. 81–91.

[45] O. ROUSTANT, F. GAMBOA, AND B. IOOSS, *Sensitivity analysis and generalized chaos expansions. lower bounds for sobol indices*, arXiv preprint arXiv:1906.09883, (2019).

[46] O. ROUSTANT, F. GAMBOA, AND B. IOOSS, *Parseval inequalities and lower bounds for variance-based sensitivity indices*, (2020).

[47] K. J. RYAN, *Estimating expected information gains for experimental designs with application to the random fatigue-limit model*, Journal of Computational and Graphical Statistics, 12 (2003), pp. 585–603.

[48] L. SIROVICH, *Turbulence and the dynamics of coherent structures. i. coherent structures*, Quarterly of applied mathematics, 45 (1987), pp. 561–571.

[49] I. SOBOL, *Sensitivity estimates for nonlinear mathematical models*, Math. Model. Comput. Exp., 1 (1993), p. 407.

[50] I. SOBOL' AND S. KUCHERENKO, *Derivative based global sensitivity measures and their link with global sensitivity indices*, Mathematics and Computers in Simulation, 79 (2009), pp. 3009–3017.

[51] K. WU, P. CHEN, AND O. GHATTAS, *An efficient method for goal-oriented linear bayesian optimal experimental design: Application to optimal sensor placement*, arXiv preprint arXiv:2102.06627, (2021).

[52] K. WU, T. O'LEARY-ROSEBERRY, P. CHEN, AND O. GHATTAS, *Large-scale bayesian optimal experimental design with derivative-informed projected neural network*, Journal of Scientific Computing, 95 (2023), p. 30.

[53] O. ZAHM, P. G. CONSTANTINE, C. PRIEUR, AND Y. M. MARZOUK, *Gradient-based dimension reduction of multivariate vector-valued functions*, SIAM Journal on Scientific Computing, 42 (2020), pp. A534–A558.

[54] O. ZAHM, T. CUI, K. LAW, A. SPANTINI, AND Y. MARZOUK, *Certified dimension reduction in nonlinear bayesian inverse problems*, Mathematics of Computation, 91 (2022), pp. 1789–1835.
[29] H. HOTELLING, *Relations between two sets of variates*, in Breakthroughs in statistics: methodology and distribution, Springer, 1992, pp. 162–190.

[30] X. HUAN AND Y. M. MARZOUK, *Simulation-based optimal bayesian experimental design for nonlinear systems*, Journal of Computational Physics, 232 (2013), pp. 288–317.

[31] N. B. KOVACHKI, S. LANTHALER, AND A. M. STUART, *Operator learning: Algorithms and analysis*, arXiv preprint arXiv:2402.15715, (2024).

[32] K. KOVAL, R. HERZOG, AND R. SCHEICHL, *Tractable optimal experimental design using transport maps*, arXiv preprint arXiv:2401.07971, (2024).

[33] A. KRASKOV, H. STÖGBAUER, AND P. GRASSBERGER, *Estimating mutual information*, Physical review E, 69 (2004), p. 066138.

[34] A. KRAUSE, A. SINGH, AND C. GUESTRIN, *Near-optimal sensor placements in gaussian processes: Theory, efficient algorithms and empirical studies.*, Journal of Machine Learning Research, 9 (2008).

[35] S. KUCHERENKO AND B. IOOSS, *Derivative based global sensitivity measures*, arXiv preprint arXiv:1412.2619, (2014).

[36] C. LIEBERMAN, K. WILLCOX, AND O. GHATTAS, *Parameter and state model reduction for large-scale statistical inverse problems*, SIAM Journal on Scientific Computing, 32 (2010), pp. 2523–2542.

[37] M. LOÈVE AND M. LOÈVE, *Elementary probability theory*, Springer, 1977.

[38] L. LU, P. JIN, G. PANG, Z. ZHANG, AND G. E. KARNIADAKIS, *Learning nonlinear operators via deeponet based on the universal approximation theorem of operators*, Nature machine intelligence, 3 (2021), pp. 218–229.

[39] J. L. LUMLEY, *The structure of inhomogeneous turbulent flows*, Atmospheric turbulence and radio wave propagation, (1967), pp. 166–178.

[40] Y. MADAY, N. C. NGUYEN, A. T. PATERA, AND G. S. PAU, *A general, multipurpose interpolation procedure: the magic points*, (2007).

[41] T. O'LEARY-ROSEBERRY, P. CHEN, U. VILLA, AND O. GHATTAS, *Derivative-informed neural operator: an efficient framework for high-dimensional parametric derivative learning*, Journal of Computational Physics, 496 (2024), p. 112555.

[42] H. POINCARÉ, *Sur les équations aux dérivées partielles de la physique mathématique*, American Journal of Mathematics, (1890), pp. 211–294.

[43] E. QIAN, B. PEHERSTORFER, D. O'MALLEY, V. V. VESSELINOV, AND K. WILLCOX, *Multifidelity monte carlo estimation of variance and sensitivity indices*, SIAM/ASA Journal on Uncertainty Quantification, 6 (2018), pp. 683–706.

[44] C. RADHAKRISHNA RAO, *Information and accuracy attainable in the estimation of statistical parameters*, Bulletin of the Calcutta Mathematical Society, 37 (1945), pp. 81–91.

[45] O. ROUSTANT, F. GAMBOA, AND B. IOOSS, *Sensitivity analysis and generalized chaos expansions. lower bounds for sobol indices*, arXiv preprint arXiv:1906.09883, (2019).

[46] O. ROUSTANT, F. GAMBOA, AND B. IOOSS, *Parseval inequalities and lower bounds for variance-based sensitivity indices*, (2020).

[47] K. J. RYAN, *Estimating expected information gains for experimental designs with application to the random fatigue-limit model*, Journal of Computational and Graphical Statistics, 12 (2003), pp. 585–603.

[48] L. SIROVICH, *Turbulence and the dynamics of coherent structures. i. coherent structures*, Quarterly of applied mathematics, 45 (1987), pp. 561–571.

[49] I. SOBOL, *Sensitivity estimates for nonlinear mathematical models*, Math. Model. Comput. Exp., 1 (1993), p. 407.

[50] I. SOBOL' AND S. KUCHERENKO, *Derivative based global sensitivity measures and their link with global sensitivity indices*, Mathematics and Computers in Simulation, 79 (2009), pp. 3009–3017.

[51] K. WU, P. CHEN, AND O. GHATTAS, *An efficient method for goal-oriented linear bayesian optimal experimental design: Application to optimal sensor placement*, arXiv preprint arXiv:2102.06627, (2021).

[52] K. WU, T. O'LEARY-ROSEBERRY, P. CHEN, AND O. GHATTAS, *Large-scale bayesian optimal experimental design with derivative-informed projected neural network*, Journal of Scientific Computing, 95 (2023), p. 30.

[53] O. ZAHM, P. G. CONSTANTINE, C. PRIEUR, AND Y. M. MARZOUK, *Gradient-based dimension reduction of multivariate vector-valued functions*, SIAM Journal on Scientific Computing, 42 (2020), pp. A534–A558.

[54] O. ZAHM, T. CUI, K. LAW, A. SPANTINI, AND Y. MARZOUK, *Certified dimension reduction in nonlinear bayesian inverse problems*, Mathematics of Computation, 91 (2022), pp. 1789–1835.