{'fraction_non_alphanumeric': 0.042975369458128077, 'fraction_numerical': 0.011152709359605911, 'mean_word_length': 4.273379052369077, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 16, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Efficient D-optimal design of experiments for infinite-dimensional Bayesian linear inverse problems

We develop a computational framework for D-optimal experimental design for PDE-based Bayesian linear inverse problems with infinite-dimensional parameters. We follow a formulation of the experimental design problem that remains valid in the infinite-dimensional limit. The optimal design is obtained by solving an optimization problem that involves repeated evaluation of the log-determinant of high-dimensional operators along with their derivatives. Forming and manipulating these operators is computationally prohibitive for large-scale problems. Our methods exploit the low-rank structure in the inverse problem in three different ways, yielding efficient algorithms. Our main approach is to use randomized estimators for computing the D-optimal criterion, its derivative, as well as the Kullback--Leibler divergence from posterior to prior. Two other alternatives are proposed based on a low-rank approximation of the prior-preconditioned data misfit Hessian, and a fixed low-rank approximation of the prior-preconditioned forward operator. Detailed error analysis is provided for each of the methods, and their effectiveness is demonstrated on a model sensor placement problem for initial state reconstruction in a time-dependent advection-diffusion equation in two space dimensions.

1. Introduction. Mathematical models of physical systems enable predictions of future outcomes. These mathematical models-typically described by a system of partial differential equations (PDEs)-often contain parameters such as initial conditions, boundary data, or coefficient functions that are unspecified and have to be estimated. This estimation requires solving inverse problems where one uses experimental data and the model to infer unknown model parameters. Experimental design is the process of specifying how measurement data, used in parameter inversion, are collected. On the one hand, the quantity and quality of data can critically impact the accuracy of parameter reconstruction. On the other hand, acquisition of experimental data is an expensive process that requires the deployment of scarce resources. Hence, it is of great practical importance to optimize the design of experiments so as to maximize the information gained in the parameter estimation or reconstruction.
Optimal experimental design (OED) controls experimental conditions by optimizing certain design criteria subject to physical or budgetary constraints. In large-scale applications, using commonly used experimental design criteria requires repeated evaluations of functionals (e.g., trace or determinant) of high dimensional and expensiveto-apply operators-applications of these operators to vectors involve expensive PDE solves. This is a fundamental challenge in OED for large-scale PDE-based Bayesian inverse problems for which we develop novel algorithms.
In our target applications, the Bayesian inverse problem seeks to infer an infinitedimensional parameter using experimental measurements, which are collected at a set of sensor locations. The OED problem aims to find an optimal placement of sensors that optimize the statistical quality of the inferred parameter. We focus on D-optimal design of sensor networks at which observational data are collected. That is, the sensor placement problem as a D-optimal design problem. Section 4 is devoted to developing efficient algorithms for computing the KL divergence from posterior to prior, and for evaluating the D-optimal objective function and its gradient; also, in that section, we formulate the optimization problem for finding D-optimal sensor placements. In section 5, we outline our model sensor placement problem for initial state inversion in a time-dependent advection-diffusion equation. Section 6 contains our numerical results, where we illustrate effectiveness of our methods. Finally, in section 7, we provide some concluding remarks.
2. Background. In this section, we provide the requisite background material required for the formulation and numerical computation of optimal experimental designs for infinite-dimensional Bayesian inverse problems.

Bayesian linear inverse problems.
Here, we describe the main ingredients of a Bayesian linear inverse problem, and our assumptions about them. The inference parameter is denoted by θ and is an element of an infinite-dimensional real separable Hilbert space V. We consider a Gaussian measure µ pr = N (θ pr , C pr ) as the prior distribution law of θ. The prior mean θ pr is assumed to be a sufficiently regular element of V, and C pr : V → V a strictly positive self-adjoint trace-class operator. Following [12,18,40], we consider prior covariance operators of the form C = A −2 , where A is a Laplacian-like operator (in the sense defined in [40]). This ensures that, in two and three space dimensions, C is trace-class.
We consider observations taken at n s measurement points, which henceforth we refer to as sensor locations, and at n t discrete points in time. We denote the vector of experimental (sensor) data by y ∈ R ny , where n y = n s n t . We assume an additive Gaussian noise model, y = F θ + η, η ∼ N (0, Γ noise ) .
Here Γ noise ∈ R ny×ny is the noise covariance matrix, and F : V → R ny is a linear parameter-to-observable map, whose evaluation involves solving the governing PDEs, and applying an observation operator that extracts solution values at the sensor locations and at observation times. Due to our choice of the noise model, the likelihood probability density function (pdf), π like (y|θ), is given by The solution of a Bayesian inverse problem is the posterior measure, µ y post , which describes the distribution law of the inference parameter θ conditioned on the observed data. The Bayes Theorem describes the relation between the prior measure, the data likelihood, and the posterior measure. In the infinite-dimensional Hilbert space settings, the Bayes Formula is given by [40], dµ y post dµ pr ∝ π like (y | θ).
Here, dµ y post dµpr is the Radon-Nikodym derivative [48] of µ y post with respect to the µ pr . In the present work, the parameter-to-observable map is linear and is assumed to be continuous. Under these assumptions on the parameter-to-observable map, the noise model, and the prior, it is known [40,Section 6.4] that the solution µ y post of the Bayesian linear inverse problem is a Gaussian measure µ y post = N (θ MAP , C post ) whose mean and covariance operator are given by Note that in this case the posterior mean coincides with the MAP estimator.
2.2. Discretization of the Bayesian inverse problems. Following the setup in [3,12], we state the finite-dimensional Bayesian inverse problem so that it is consistent with the underlying inference problem with V = L 2 (D) as its parameter space.
Finite-dimensional Hilbert space setup. We use a finite-element discretization of the Bayesian inverse problem, and hence work in the n-dimensional inner product space V n , which is R n equipped the mass-weighted inner product, x, y M = Mx, y ; here ·, · denotes the Euclidean inner product, M is the finite-element mass matrix, and n is the dimension of the discretized inversion parameter; see [3,12].
We denote by L(V n ), and L(R k ) the spaces of linear operators on V n , and R k , respectively. We also define the spaces of the linear transformations from V n to R k , and from R k to V n by L(V n , R k ), and L(R k , V n ), respectively. Here R k is assumed to be equipped with the Euclidean inner product. Next, we give explicit formulas for the adjoints of elements of these spaces, which will be useful in the rest of the paper: The superscript T denotes matrix transpose. Also, we use the notations L sym (R n ) and L sym (V n ) for the subspaces of self-adjoint operators in L(R n ) and L(V n ), respectively.
Remark 2.1. It is straightforward to note that for A ∈ L sym (V n ), Thus, the similarity transform A → M 1/2 AM −1/2 gives a matrix that is self-adjoint with respect to the Euclidean inner product; this is useful in practical computations.
The discretized inverse problem. We denote by θ ∈ V n the discretized parameter that we wish to infer, and as before, let y ∈ R ny be a vector containing sensor measurements, i.e., the experimental data. As discussed in [12], the discretized inverse problem is stated on the space V n . The prior density is given by Following [12], we define A = M −1 L, with where K is the finite-element stiffness matrix, and α and β are positive constants. Note that A ∈ L sym (V n ) and Γ prior = A −2 . The posterior measure is also Gaussian with mean and covariance operator given by [12,40], where F ∈ L(V n , R ny ) is the discretized parameter-to-observable map (forward operator). Note that (7) is the discretized version of (1). As is well known, the posterior mean θ y post is the unique global minimizer of the regularized data misfit functional The Hessian operator of this functional is noise F is the data misfit Hessian [12,42]. Note also that

The operator
prior ∈ L sym (V n ) is referred to as the prior-preconditioned data misfit Hessian [12] and will be of relevance in the sequel. Following Remark 2.1, we also define which is a symmetric matrix and is convenient to use in computations. The rank of H H H m is min{n y , n}. For some applications, even though the rank can be large, the eigenvalues of H H H m decay rapidly. This is further discussed in [19]. We also define the prior-preconditioned forward operator Remark 2.2. Preconditioning the forward operator by the prior covariance operator is expected to lead to faster decay of singular values of the preconditioned operator, when using a smoothing Γ prior . This can be explained using the multiplicative singular value inequalities [10, page 188], The first inequality shows that preconditioning by Γ 1/2 prior does not degrade the decay of singular values of F, and the second one shows that a smoothing Γ prior for which singular values σ i (Γ 1/2 prior ) decay faster than that of F can improve the spectral decay of F. The latter was observed numerically in many cases; see e.g. [3,19].
Remark 2.3. The expression for H H H m given above involves square root of M. While this can be computed for some problems, a numerically viable alternative is to compute a Cholesky factorization M = LL T , where L is lower triangular with positive diagonal entries. Then, it is simple to show that the transformation A → L T AL −T , maps elements of L sym (V n ) to L sym (R n ). This alternate formulation can be used to redefine H H H m as an element of L sym (R n ). Computing explicit factorizations may be prohibitively expensive for large-scale problems in three-dimensional physical domains. In this case, the action of M 1/2 (and its inverse) on vectors can be computed using polynomial approaches [14,15] or rational approaches [8,23].

Randomized log-determinant estimators.
Here we describe the randomized log-determinant estimator developed in our previous work [38]. These estimators will be used for efficient estimation of log det(I + H H H m ). Recall that H H H m is symmetric positive semidefinite and has rapidly decaying eigenvalues-a problem structure that is exploited by the randomized algorithms discussed here. Moreover, H H H m is a dense operator that is too computationally expensive to build elementwise, but matrix-vector products involving H H H m can be efficiently computed. To compute the log-determinant of I +H H H m , we first use randomized subspace iteration to estimate the dominant subspace. Algorithm 1 is an idealized version of randomized subspace iteration used for this purpose. As the starting guess, we sample random matrix Ω with ≥ k columns, with i.i.d. entries from the standard normal distribution. In practice, = k + p, with p a small (e.g., p = 10) oversampling parameter. We then apply q steps of the power iterations to this random matrix using H H H m to obtain Y. A thin QR decomposition of Y produces a matrix Q with orthonormal columns. The output of Algorithm 1 is the × restriction of H H H m to Span(Q), i.e., T = Q T H H H m Q. In practice, the number of power iterations is taken to be q = 1, or q = 2. Specifically, for all the numerical tests in the present work, q = 1 was found to be sufficient.
Remark 2.4. Throughout this paper, we assume that the random starting guess Ω is a standard Gaussian random matrix. An extension to other random starting guesses, such as Rademacher, can be readily derived, following the approach in [38].

Algorithm 1 Randomized subspace iteration (idealized version)
Input: Hermitian positive semi-definite matrix H H H m ∈ R n×n with target rank k, number of subspace iterations q ≥ 1, starting guess Ω ∈ R n× with k ≤ ≤ n columns.
With the output of Algorithm 1, we can approximate the objective function as log det(I + H H H m ) ≈ log det(I + T).
We also mention that this algorithm can be used to compute an (approximate) lowrank approximation of H H H m as follws: compute the eigenvalue decomposition of T = U T Λ T U T T , and compute U = QU T . Then we have the low-rank approximation For the error analysis, we will need the following definitions. Since H H H m is symmetric positive semidefinite, its eigenvalues are nonnegative and can be arranged in descending order. Consider the eigendecomposition H H H m = UΛU T ; let us suppose that Λ 1 ∈ R k×k contains the dominant eigenvalues and Λ 2 ∈ R (n−k)×(n−k) , where k is the target rank. Partition the eigendecomposition as (10) H We assume that Λ 1 is nonsingular, and that there is a gap between the eigenvalues λ k and λ k+1 . The size of the gap is inversely proportional to As noted before, the rank of H H H m , which we denote by K, satisfies K = min{n y , n}. 6 3. D-optimal design of experiments for PDE based inverse problems. In this section we examine the information theoretic connections of the D-optimality criterion, and formulate the sensor placement problem as an OED problem.
3.1. Kullback-Leibler divergence. The Kullback-Liebler (KL) divergence [30] is a "measure" of how one probability distribution deviates from another. It is important to note that it is not a true metric on the set of probability measures, since it is not symmetric and does not satisfy the triangle inequality [41]. Despite these shortcomings, the KL divergence is widely used since it has many favorable properties. In the context of Bayesian inference, the KL divergence measures the information gain between the prior and the posterior distributions. We motivate the discussion by recalling the form of the KL divergence from posterior to prior distributions in the finite-dimensional case. To avoid introducing new notation, we continue to denote the (finite-dimensional) posterior and prior measures by (11) µ pr = N (0, Γ prior ) , µ y post = N θ y post , Γ post , respectively. Here, for convenience, and without loss of generality, we assumed the prior mean is zero. The following expression for D KL µ y post µ pr is well known [41, where we have used the notation, x, y Γ −1 We note that the expression (12) is not meaningful in the infinite-dimensional case, and it is not immediately clear what the infinite-dimensional limit will be. In [2], this issue is investigated in detail in the infinite-dimensional Hilbert space setting. Here we provide a derivation of an expression for the KL divergence, within the context of the discretized problem, that remains meaningful in infinite dimensions. The log det term in (12) may be simplified as: Since H H H m is the discretization of a trace-class operator, det(H H H m + I) is the Fredholm determinant, and is well-defined. Next, we consider the term −n + tr (Γ −1 prior Γ post ). Since n = tr (I) and the trace operator is linear, where in the penultimate step we used the relation Γ −1 post = H m + Γ −1 prior . Notice that the argument of the trace in the final expression is in fact a trace-class operator in 7 the infinite-dimensional case and has a well-defined trace. Combining (13) and (14) we rewrite (12) as follows: While this equation is equivalent to (12) in finite dimensions, it will be used henceforth, because it is well defined in the infinite dimensional limit.
3.2. Expected information gain and the D-optimal criterion. The following result provides an explicit expression for the expected information gain leading to the D-optimal criterion. Theorem 1. Let µ pr = N (0, Γ prior ) and µ y post = N θ y post , Γ post be the prior and posterior measures corresponding to a Bayesian linear inverse problem with additive Gaussian noise model. Then, Proof. See Appendix A.1.
Several important observations are worth mentioning here. First, this theorem provides the expression for the D-optimal criterion in terms of the expected KL divergence from the posterior to the prior distribution; this expression remains meaningful in the infinite-dimensional setting. Moreover, (16) reveals an important problem structure in the D-optimal criterion: the operator H H H m admits a low-rank representation in a wide class of inverse problems. This low-rank structure is exploited in our algorithms for numerical solution of the D-optimal experimental design problem. Finally, the KL divergence (15) depends on the data, through the term θ y post . For the linear-Gaussian problem, which we consider in this paper, the D-optimality criterion is independent of the data since the expectation of the KL divergence is taken over the data and the prior.
In the context of this paper, an optimal experimental design is one that maximizes the D-optimal criterion, or alternatively, maximizes the expected information gain. Moreover, in practical computations, we can use H H H m in replace of H H H m in the D-optimal criterion (16) and the expression for the KL divergence (15).
3.3. The D-optimal criterion for sensor placement problems. For the sensor placement problem, we consider the following setup [3,5,20,22,45]. A set of points x i , i = 1, . . . , n s , representing the candidate sensor locations, is identified. Each sensor location is associated with a non-negative weight w i ∈ R encoding its relative importance. The experimental design is then fully specified by the vector For sensor placemement, binary weight vectors are desired. If w i equals one, then a sensor will be placed at x i ; otherwise, no sensors will be placed at that location. In this context, nonbinary weights are difficult to interpret and implement. From a computational standpoint, however, the experimental design problem with binary weights is a combinatorial problem and intractable even for a modest number of sensors. Therefore, as in [3,5], we consider a relaxation of the problem with weights w i ∈ [0, 1], i = 1, . . . , n s . Subsequently, binary weights are obtained using sparsifying penalty functions; see section 4.3.
We next describe how the design w is incorporated in the D-optimal criterion. The weight vector w enters the Bayesian inverse problem through the data likelihood. 8 Subsequently, following the setup in [3], the data misfit Hessian takes the form, where W is an n s n t × n s n t diagonal matrix given by with e j the jth unit vector in R ns . We consider the case of uncorrelated observations 1 , and hence, Γ noise is a diagonal matrix, indicates the noise level at each sensor. For convenience, we introduce the notation, With the above definitions in place, the D-optimal objective becomes with H H H m (w) defined according to (8), Using the expression for derivative of log-determinant of matrix-valued functions [45,Theorem B.17], and letting ∂ j be shorthand for ∂ ∂wj , the derivative of the above function is given by, (9). The last expression follows since ∂ j W σ = E σ j . 4. Numerical algorithms. In this section, we present our proposed numerical algorithms for efficient computation of KL divergence (section 4.1), and the D-optimal criterion and its derivative (section 4.2). The optimization problem for D-optimal sensor placements is formulated in section 4.3.

Fast estimation of KL divergence.
To estimate the KL-divergence, we propose two methods-based on the the truncated spectral representation of H H H m or using a randomized estimator.
Truncated spectral decomposition. We can estimate the KL divergence, using only the first k exact eigenvalues of H H H m . Let Λ 1 and Λ 2 be as in (10). We propose the following estimator: Since the posterior is Gaussian, the MAP estimator coincides with the conditional mean, and θ y post is obtained by solving (7a). Therefore, θ y post , θ y post Γ −1 prior is assumed to be available and need not be estimated using the eigenvalue decomposition. The error in the KL-divergence, defined as follows The proof is a straightforward consequence of the properties of the trace and the determinant and is omitted.
Randomized approach. Computing the exact eigenvalues of H H H m can be computationally challenging. Based on the approach in section 2.3, we propose a randomized estimator for the KL divergence, which is cheaper to compute. Given the output T of Algorithm 1, our estimator for the KL divergence is It is clear that we only need to compute the determinant and the trace of an × matrix instead of an n × n matrix. This has significant computational benefits when n. Theorem 2 presents the error of the KL decomposition with the expectation taken over the random starting guess Ω. A similar result can be derived for concentration, or tail bounds of the distribution, but is omitted.
Theorem 2. Let T be computed using Algorithm 1 with Gaussian starting guess Ω ∈ R n×(k+p) with p ≥ 2 and D KL µ y post µ pr be computed using (24). Then, we have the following estimate for the expected value of the error The result can be further simplified with the observation that for x ≥ 0, log(1 + x) ≤ x. Therefore, the upper bound in Theorem 2 can be simplified to The importance of this result is that the expected error incurred in the KL divergence is no worse than the error in the approximation of tr (H H H m ). However, it is worth comparing the result of (23) with Theorem 2. The error in Theorem 2 is higher; however, the estimator (24) is much easier to compute.

Efficient computation of OED objective and gradient.
To efficiently estimate the OED objective function and gradient, we propose three different methods: truncated spectral decomposition approach (section 4.2.1); randomized approach (section 4.2.2); and fixed (frozen) low-rank approximation of F F F (section 4.2.3); then, in section 4.2.4, we discuss the computational costs and pros and cons of each approach.

Spectral decomposition approach.
We show how to compute the objective function J(w) and its derivatives ∂ j J(w) using the dominant eigenpairs of H H H m . Using the spectral decomposition of H H H m in section 2.3, we can compute By Sherman-Morrison-Woodbury formula, we have ( Substituting this expression into the gradient (19), we have, for j = 1, . . . , n s where Z j is as in (20).
the constants {z j } ns j=1 are independent of the weights and can be precomputed; see section 4.2.4 for more details. To summarize, the gradient can be computed as The estimators for the D-optimal objective and its gradient, denoted, respectively, by J eig (w) and ∂J eig (w) only retain k ≤ K eigenpairs 2 . In what follows, we refer to this as the Eig-k approach, where k is the target rank of the approximation. The steps for the computation of the D-optimal objective and gradient using the Eig-k approach are detailed in Algorithm 2.
Input: Constants {z j } ns j=1 , target rank k, and design w. Output: OED objective J eig (w) and gradient ∇J eig (w).
1: Solve an eigenvalue problem for eigenpairs Error analysis. It should be clear that if k = K, the objective function and gradient evaluation is exact. For the case of k < K, we record the following result summarizing the approximation errors in objective function and gradient.
Theorem 3. Let J eig and ∇J eig (w) respectively be the approximation to the OED objective and gradient computed by Algorithm 2, with a rank k truncation. Then, for every w ∈ R ns + , with Z j defined in (20), and for j = 1, . . . , n s 2 Here and henceforth, we denote the approximate quantities with a symbol.
Here λ i are eigenvalues of H H H m and Λ 2 is defined in (10).
The interpretation of this theorem is as follows: if the truncated eigenvalues, {λ j } K j=k+1 are very small the errors are negligible. In practical computations, we choose an appropriate target rank k based on the spectral decay of H H H m .
Theorem 3 gives the elementwise error in the gradient. If the error in the 2−norm of the gradient is desired, it can be obtained from Note also that, due to boundedness of the prior-preconditioned forward operator, { Z j } ns j=1 remain bounded with respect to mesh refinements. 4.2.2. Randomized approach. The matrix H H H m is symmetric positive semidefinite and has rapidly decaying eigenvalues (and is exactly of rank K = min{n, n s n t }). Under these conditions, the randomized estimator described in section 2.3 is likely to be accurate. Hence, we use it to compute estimate the log-determinant log det(I+H H H m ).
Let T be the output of of Algorithm 1 applied to H H H m . The objective function can be approximated as (32) J(w) ≈ J rand (w) ≡ log det(I + T).
Next, we consider approximation of the derivative. Following the approach in section 2.3, we can approximate H H H m ≈ QTQ T = U Λ U T . Then, using the Woodbury matrix identity, we have (1 +λ 1 ), . . . ,λ /(1 +λ )) ∈ R × , whereλ i are the eigenvalues of T. This approximation can be used to estimate the gradient ∂ j J(w) as follows. Write The constants z j are as before and can be precomputed in the same manner. The steps for the computation of the D-optimal objective and gradient using the randomized approach are detailed in Algorithm 3. Error analysis. The error in the objective function approximation, when using the randomized approach, is quantified by the following result: Theorem 4. Let T be computed using Algorithm 1 with Gaussian random starting guess Ω ∈ R n×(k+p) with p ≥ 2 andĴ(w) be computed using (32). Then, we have the following estimate for the expected value of the error (35) E Ω |J(w) − J rand (w)| ≤ log det(I + Λ 2 ) + log det( where C ge is defined in Theorem 2.
1: Apply Algorithm 1 with = k + p and q to obtain T ∈ R × and Q ∈ R n× . 2: Evaluate J rand (w) = log det(I + T).
Next, we present a result that quantifies the error in the gradient using the randomized approach.
Theorem 5. Let T be computed using Algorithm 1 with starting guess Ω ∈ R n×(k+p) with p ≥ 2 and the approximate derivative ∂ j J rand (w) defined as (34). The expected error in the approximation to ∂ j J(w) is where C ge is defined in Theorem 2, Z j is defined in (20) and j = 1, . . . , n s .
The interpretation of Theorem 5 is similar to Theorems 2 and 4. The estimator is accurate if the eigenvalues contained in Λ 2 are small relative to Λ 1 . It is worth mentioning that ∇J rand (w) is not the exact derivative of the approximate objective function J rand (w) but its accuracy is comparable to the exact gradient. This is also illustrated numerically in section 6.1.
Theorem 5 gives the component-wise error in the gradient. Under the same assumptions, it can be shown that (36) E This follows from the vector norm inequality · 2 ≤ · 1 and the linearity of expectations. Compared to (31), this result is clearly suboptimal. It is not clear, if the bounds can be further tightened.

Frozen low-rank approximation of F F F.
The dominant cost in the evaluation of the objective function and the derivative is the applications of F F F to vectors. This operator is typically low-rank for many applications. Since the operator F F F is independent of the weights w, a low-rank approximation to F F F can be precomputed at the start of the optimization routine. We refer to this as the "Frozen low-rank approach." The idea of precomputing a low-rank SVD was used in [3,22], for computing A-optimal designs. Here we show how this can be used efficiently for D-optimal experimental design, and provide a result quantifying the approximation error. 13 The matrix F F F is approximated by the rank-k thin SVD F F F ≡ U Σ V * , where U ∈ R nsnt×k has orthonormal columns, Σ ∈ R k×k and V ∈ R n×k has columns that are orthonormal with respect to the mass-weighted inner product. With this low-rank approximation, the approximate objective function is defined to be Evaluating J froz (w) still requires computation of the log-determinant. Fortunately, the low-rank approximation can be used to evaluate this function efficiently. Note that computing the objective function using the Frozen low-rank approach only requires the evaluation of the determinant of an × matrix, instead of an n × n matrix. The computation of the derivative, denoted by ∇J froz (w), is similar to the discussion in the exact (and randomized) case and is omitted. We derive the error in the Frozen low-rank approximation when F F F is the best rank−k approximation to F F F.
Theorem 7. Let F F F = F F F k be the rank−k truncated singular value decomposition to F F F. Then, assuming w i ∈ [0, 1], i = 1, . . . , n s , The interpretation of this theorem is that the Frozen low-rank approach is accurate so long as the discarded singular values of F F F are small. Moreover, this result is notable in the sense that the error is independent of the weight vector.
The precomputation of F F F is the dominant cost of this algorithm and requires O(2k) PDE solves. The main computational benefit of this approach is that during the optimization iterations no PDE solves are necessary, making this algorithm very fast in practice. The disadvantages are that the errors in the objective function and gradient computation cannot be controlled since F F F is fixed; moreover, storing the low-rank approximation can be expensive. The other major disadvantage is that this approach is specific to our formulation of the optimal sensor placement problem. The assumption that F F F is independent of the design parameters may not hold in general. The Eig-k and the randomized approaches do not make this assumption and are thus more generally applicable. Despite its shortcomings, the Frozen low-rank approach can be viable and very efficient for some problems. Finally, we note that a low-rank approximation of F F F can also be used for fast computation of the KL divergence (15).

Computational costs.
Here we discuss precomputation of the constants {z j } ns j=1 defined in (27), and the overall computational cost of our proposed methods.
Precomputing {z j }. The identity matrix can be expressed as the sum of outer products of vectors {v m } nt m=1 which form the standard canonical basis for R nt , i.e., Then with e j and E j defined in (17) Using the cyclic property of the trace and with f jm ≡ F F FM −1/2 (v m ⊗ e j ) Computing each z j requires applying F F FM −1/2 to n t matrices v m ⊗ e j ; however, since v m ⊗ e j only has n t nonzero columns, precomputing {z j } ns j=1 costs n s n t PDE solves. Summary of Computational Costs. Here we discuss the computational cost of objective and gradient evaluation using our proposed approaches. For the spectral approach, the objective function evaluation requires O(k) matrix-vector products (matvecs) involving the transformed forward operator F F F, O(k) matvecs involving its adjoint F F F * ; derivative evaluation requires k PDE solves and an additional computational cost of O(k 2 n) flops. The computational cost of the randomized approach is similar. The computational cost of our proposed methods is summarized in Table 1.

Method/Component Forward solves Adjoint solves Precomputation
Naive Table 1 Summary of computational costs measured in terms of PDE solves. Here K = min{nsnt, n} and k ≤ K is the target rank. Furthermore, = k + p, where p ≥ 0 is an oversampling parameter and q is the number of subspace iterations (see Algorithm 1). Table 1: 1. The first two methods, "Naive" and Eig-K involve no approximations. In the "Naive" method, the matrix H H H m is constructed explicitly, using which the objective function and the gradient are computed. By contrast, the Eig-K approach only computes the nonzero eigenpairs of the matrix. The remaining three methods have some approximation built into them. 2. The three approximate methods require some form of precomputation. The Frozen approach requires O(k) PDE solves, whereas the Eig-k and Randomized methods require K PDE solves (to precompute {z j }); these costs can be lowered to O(k) by using a low-rank approximation F F F instead of F F F. 3. It was mentioned earlier that the derivative computations require an additional O(k) PDE solves. However, if the approximate truncated spectral representation is computed using either a Krylov subspace method, or the randomized approach, then the information content in the intermediate steps can be reutilized to avoid the additional PDE solves. 4. In the Eig-k approach, the k eigenpairs are computed using an iterative matrix-free method such as a Krylov subspace solver. The asymptotic computational cost of the randomized approach is comparable with that of Eig-k approach. However, in the randomized subspace iteration, the matvecs involving H H H m are readily parallelizable. In contrast, Krylov subspace methods are inherently sequential. 5. Based on the overall number of PDE solves, the Frozen approach is very attractive provided the low-rank approximation F F F remains sufficiently accurate throughout the optimization iterations. Table 1 assume that the target rank is the same for all the three approaches, but does not say anything about the accuracy of each approach. The Frozen approach may overestimate the rank k needed to accurately approximate F F F. To illustrate this, suppose during the iteration history (or at the optimal point), only one sensor is active. This means that the formal rank of F F F is exactly K, but the rank of H H H m is n t . The Eig-k and randomized approaches, on the other hand, target the rank of H H H m rather than F F F and therefore, can be much more accurate with a lower computational cost. 7. In the spectral and randomized methods, the target rank k remains fixed throughout the optimization iterations. Development of variations of these approaches, in which the rank is adaptively adjusted during the optimization iterations, is subject of our future work. Each of the proposed methods have their own advantages and disadvantages. The choice of the method used would depend on the context and the specific application under consideration. We advocate the randomized approach since it provides an excellent balance of accuracy, flexibility, efficiency, and ease of implementation.

The computational costs in
4.3. The optimization problem for finding D-optimal designs. The previous subsections were concerned with efficient methods for computation of the Doptimal criterion and its derivatives with respect to the design weights. Here we formulate the optimization problem to be solved for finding D-optimal designs. Let J(w) be as in (18). We formulate the optimization problem for finding a D-optimal experimental design as follows: where W = [0, 1] ns , γ > 0, and P (w) is a sparsifying penalty function. A typical approach is to to use an 1 -norm penalty; in this case, since components of w are non-negative, an 1 -penalty simplifies to P(w) = ns j=1 w j . To explicitly enforce binary designs we adopt two strategies. The first strategy is to use an 1 -norm penalty function, and then threshold the computed optimal weights based on a heuristic; that is, we place sensors only in locations whose corresponding weight is above a given threshold. The second strategy is to solve a sequence of optimization problems with penalty functions {P ε } ncont i=1 that successively approximate the 0 -"norm" as done in [3]. In the continuation approach, each successive optimization problem uses the result of the preceding problem as initial guess-a process commonly referred to as "warm starting." In our computations, we solve the OED optimization problem using Matlab's interior-point solver provided by the fmincon function; the objective function and its gradient, computed using our algorithms, are supplied to the optimization solver, and BFGS approximation to the Hessian is used.

Model problem.
In this section, we describe the model problem used to illustrate our proposed OED methods. The forward problem is a time-dependent advection-diffusion equation, and the inverse problem is the inference of the initial state from sensor measurements of the state variable at discrete points in time; see also, [1,3,19,35]. The setup of the model problem below is mainly based on [3].
The governing PDE and the parameter-to-observable map. Let D ⊂ R 2 , be a twodimensional domain depicted in Figure 1(left). The domain boundaries ∂D include the outer edges as well as the internal boundaries of the rectangles that model buildings/obstacles. Given an initial state θ, we solve a time-dependent advection-diffusion where, κ > 0 denotes the diffusion coefficient, and T > 0 is the final time. We use κ = 0.001 and T = 5 in the numerical experiments below. The velocity field v, shown in Figure 1(left), is computed by solving a steady Navier-Stokes equation as in [3,35]. In what follows, we also need the action of the adjoint F * of the parameterto-observable map F to vectors. For a given observation vectorȳ ∈ R ny , F * ȳ is computed by solving the adjoint equation (see [1,19,35]) for the adjoint variable p = p(x, t), and setting F * ȳ = −p(·, 0). Note that this equation is a final value problem that is solved backwards in time.
Bayesian inverse problem formulation. The Bayesian inverse problem seeks to infer the initial state θ from point measurements of u. Following the setup outlined in section 2, we utilize a Gaussian prior measure µ pr = N (θ pr , C pr ), with C pr as described before, and use an additive Gaussian noise model. Specifically, we use α = 2×10 −3 and β = 10 −1 in (6). The "true" initial state shown in Figure 1(left) is used to synthesize data, and the noise level is set as follows. We solve the forward model using the true initial state and record measurements of the state at the candidate sensor locations and at the observation times (see numerical results section for specifics), and define the noise level as two percent of the maximum of the recorded measurements.
Discretization. We discretize the forward and adjoint problems via linear triangular continuous Galerkin finite elements in the two-dimensional spatial domain, and use the implicit Euler method for time integration. As in [3,35], we follow a discretizethen-optimize approach, where the discrete adjoint equation is obtained as the adjoint of the discretized forward equation.
6. Numerical results. In this section, we test various aspects of our proposed numerical methods for D-optimal sensor placement. Section 6.1 is devoted to testing accuracy of our randomized estimators, and a simple illustration of D-optimal sensor placements using the randomized approach. In section 6.2 we consider a more realistic example in which we allow a denser grid of candidate sensor locations, and the unknown parameters are discretized on a grid of higher resolution. The D-optimal sensor placement methods are illustrated on this application.

Test of accuracy and basic illustrations.
In this subsection, we use the following parameters that specify the test problem: we pick the grid of n s = 35 candidate sensor locations, depicted in Figure 1 (middle), with n t = 3 observation times given by t = 1, t = 2, t = 3.5, for a total of 105 observations. The discretized parameter dimension for this example is n = 1018.
Accuracy of KL-divergence estimator. Here we demonstrate the accuracy of the estimators for the computation of the Kullback-Liebler divergence. We take w to be a 35 × 1 vector of all ones; i.e., all the sensors are active. We estimate the KL divergence using two different methods proposed in section 4.1. In the first method the KL divergence is computed using "exact" eigenvalues (21); the approximation is because only k ≤ K eigenvalues are used for estimating the KL divergence. The "exact" eigenvalues were computed using eigs function in MATLAB. In the second method we compute the KL divergence using the randomized estimator (24). The oversampling parameter was fixed to be p = 5, we increase the total number of random samples = k + p between 20 and 105; since the oversampling parameter is fixed, this amounts to increasing the target rank k. The error in the KL divergence is plotted in Figure 2. As we see, the error decays rapidly as the number of computed eigenvalue increases, and the randomized estimators is nearly as accurate as using the "exact" eigenvalues, but is, in general, considerably more efficient to compute.
Accuracy of D-optimal criterion and its derivative. The setup is the same as the previous experiment; now, we consider the accuracy of the D-optimal objective function and the gradient. We use Eig-k and randomized estimators and the results are plotted in Figure 3 (left, middle). Similar conclusions are drawn here as well: the error decreases with increasing target rank k and the accuracy of the randomized estimators is comparable with that of Eig-k. Note that since we used W = I, the frozen and Eig-k approaches are essentially equivalent, for the present test. In our next experiment, we plot the accuracy of the objective function as a function of mesh refinement in Figure 3 (right). The accuracy of the estimators does not degrade with increasing mesh refinement. This desirable property is also a numerical illustration of the fact that our formulation remains valid in the infinite-dimensional limit. Computing D-optimal designs. We have demonstrated that our estimators are accurate; we now show how the corresponding designs look like. As a first illustration, we use the randomized approach in section 4.2.2, where OED objective and gradient are computed according Algorithm 3, to compute a D-optimal sensor placement. The 1 -norm is chosen as the penalty function to ensure sparsification (see section 4.3). The resulting design has weights between 0 and 1. To enforce binary weights, the following thresholding criterion was used: if the sensor weight w i satisfies w i / j w j ≥ 3 × 10 −2 then it is set to be 1 else 0. The resulting designs are plotted in Figure 4 (left). We also plot the relative errors in objective and gradient computations during the course of the optimization iteration history in Figure 4 (right). The latter shows that the randomized estimators remain accurate over the course of the iterations.
To compare the different methods proposed in the present work, we compute optimal design weights using spectral, randomized, and frozen low-rank approaches; see Figure 5. To remain consistent across methods, for each method we use rank k = 70 approximations with an oversampling parameter of p = 5. We use a log-scale on the vertical axis, and plot the computed weights in descending order; reported also are the respective optimal objective values. As expected, the randomized and spectral approaches lead to nearly identical optimal weights; on the other hand, the solution obtained via the frozen approach, agree with that of the randomized and spectral approach only for the larger weights.
6.2. A higher resolution example. Here we illustrate computing a D-optimal sensor placement on a problem with parameter dimension n = 2605, and a grid of n s = 109 candidate sensor locations; see Figure 1(right). To enforce binary weights, . Left: D-optimal design with 15 active sensors; the sensor placement was obtained using an 1 -norm penalty approach, with penalty parameter γ = 9. Sensors are placed at locations whose corresponding sensor weight exceeded w i / j w j ≥ 3 × 10 −2 . Right: relative error in OED objective and gradient over optimization iterations.  we used the continuation approach described in section 4.3. We picked the penalty parameter γ = 1.8 and the penalty functions P ε is as in [3,Section 4.5]. We use six continuation steps with {P εi } ncont i=1 , ε i = 1/2 i . The resulting D-optimal sensor placement is shown in Figure 6 (left) and the convergence of the successive design vectors to a binary weight vector is shown in Figure 6 (right).
Testing the effectiveness of the computed D-optimal design. We illustrate the effectiveness of the computed D-optimal sensor placement by comparing it with an ensemble of randomly generated sensor placements. Since the continuation approach yielded 31 active sensors, the random sensor placements had the same number of sensors. For each design (random ones and the optimal one), we compute the objective value J(w), as well as the KL divergence from the posterior to prior. To compute the KL divergence, the inverse problem was solved for each of the sensor placements and the randomized estimator (24) was used.
In the first place, it is expected that the computed optimal design outperform random designs in terms of smaller OED objective value. Moreover, Theorem 1 indicates that maximizing J(w) results in maximizing the expected information gain. Both of these issues are illustrated in Figure 7. The red dots in that figure correspond to 500 randomly generated sensor configurations with 31 sensors and the black dot to the computed optimal design. That the black dot falls below all the red dots is to be expected, because the quantity in the vertical axis, −J(w) = − log det(I + H H H m (w)), is what we sought to minimize. On the other hand, the computed D-optimal sensor placement maximizes the information gain in the average sense specified in Theorem 1; this explains the location of the black dot in the horizontal axis. ) and the information gain, computed at the optimal design (black dot) with 500 randomly generated designs (red dots).
Solving the inverse problem using the computed optimal design. Finally, we report the results of solving the Bayesian inverse problem via the computed D-optimal sensor placement. Figure 8 shows the "true" parameter (initial state) and the computed MAP estimator. Figure 9 (top) compares the prior and posterior standard deviation fields, showing reduction in uncertainty and the information gain by using an optimal sensor placement. Figure 9 (bottom) shows three samples drawn from the posterior distribution.

Conclusion.
We have developed a computational framework for D-optimal experimental design for PDE-based Bayesian linear inverse problems with infinitedimensional parameters. Our methods exploit low-rank structure of the prior-preconditioned data misfit Hessian, and the prior-preconditioned forward operator for fast estimation of OED objective and its gradient, as well as the KL divergence from posterior to  prior. Our developments resulted in three different approaches: the low-rank spectral decomposition approach, the randomized approach, and the approach based on fixed low-rank SVD of the prior-preconditioned forward operator. Our numerical results, aiming at D-optimal sensor placement for initial state inversion in a time-dependent advection-diffusion equation, demonstrate the effectiveness of our methods.
The randomized and spectral approaches provide accurate and efficient estimators of the OED objective and gradient; the accuracy is controlled by setting the target rank k, which is specified a priori. Some open questions remain: can we adaptively determine the target rank k during the optimization iterations so as to ensure the objective function and the gradient are accurate to a given tolerance, with as few PDE solves as necessary? Can suitable error estimators be incorporated into an optimization routine with theoretical guarantees of convergence? Can the presented strategy be extended to OED for nonlinear inverse problems? Addressing these questions is subject of our future work. recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation.
Appendix A. Proofs.
A.1. Proof of Theorem 1. We define the operator S = (H H H m + I) −1 . Using the definition of θ y post we have, where we have used that Γ post = Γ This, along with (42) and (41) The stated result follows by averaging (15) over the prior and noise, as in the statement of the theorem, and using (43).
A.3. Proof of Theorems 3 and 5. First, we record the following basic lemma: Lemma 8. Let A be an n × n matrix and suppose B is a symmetric positive semidefinite matrix. Then, we have |tr (AB)| ≤ A 2 tr (B).
Proof (Theorem 3). The first statement is immediate; the second one follows from where we have abbreviated Z j ≡ M 1/2 F F F * E σ j F F FM −1/2 . The inequality is justified in Lemma (8) in the Appendix. Taking expectations and use (48) to obtain the desired bound. 26 A.4. Proof of Theorem 7. Before we turn to the proof of Theorem 7, we present the following Lemma. Since E is hpsd, it has a well-defined square root. Note that X has singular values at most 1, and is a contraction matrix. Furthermore, the multiplicative singular value inequalities [10, page 188] imply λ j (XEX) = σ 2 j (E 1/2 X) ≤ σ 2 j (E 1/2 ) = λ j (E) j = 1, . . . , n.
Writing the determinant as the product of the eigenvalues (1 + λ i (E)) = det(I + E).
Taking logarithms delivers the desired upper bound. The lower bound follows from the fact that since E is positive semidefinite, and X is Hermitian, I + XEX ≥ I. Therefore, log det(I + XEX) ≥ log det I = 0.
Proof (Theorem 7). Partition the singular value decomposition as F F F = F F F k +F F F k,⊥ , and observe F F FF F F * = F F F k F F F * k + F F F k,⊥ F F F * k,⊥ . The determinant identity [33, Corollary 2.1] gives (49) log det(I + F F F * WF F F) = log det(I + W 1/2 F F FF F F * W 1/2 ).
It can be readily verified that F F FF F F * ≥ F F F k F F F * k ; therefore, by the lower bound in Lemma 9 |J(w) − J froz (w)| = log det(I + F F F * WF F F) − log det(I + F F F * k WF F F k ).
Use (49) along with the upper bound of Lemma 9 log det(I + F F F * WF F F) − log det(I + F F F * k WF F F k ) ≤ log det(I + W 1/2 F F F k,⊥ F F F * k,⊥ W 1/2 ) = log det(I + F F F * k,⊥ WF F F k,⊥ ).
Since the weights are at most 1, then W ≤ I and log det(I + F F F * k,⊥ WF F F k,⊥ ) ≤ log det(I + F F F * k,⊥ F F F k,⊥ ) = log det(I + Σ 2 2 ).
The inequality is due to the properties of Loewner partial ordering of Hermitian matrices [25,Corollary 7.7.4]. The final equality employs the determinant identity.