Randomization and reweighted $\ell_1$-minimization for A-optimal design of linear inverse problems
Elizabeth Herman, Alen Alexanderian, and Arvind K. Saibaba
TL;DR
This paper introduces randomized methods and reweighted $ ext{l}_1$-minimization techniques to efficiently design optimal experiments for Bayesian linear inverse problems governed by PDEs, focusing on A-optimality and sparse solutions.
Contribution
It develops structure-exploiting randomized estimators for the A-optimal design criterion and proposes a novel reweighted $ ext{l}_1$-minimization algorithm for sparse, binary design vectors.
Findings
Efficient computation of A-optimal design using randomized estimators.
Successful application to a contaminant source identification problem.
Demonstrated sparsity and optimality of the resulting experimental designs.
Abstract
We consider optimal design of PDE-based Bayesian linear inverse problems with infinite-dimensional parameters. We focus on the A-optimal design criterion, defined as the average posterior variance and quantified by the trace of the posterior covariance operator. We propose using structure exploiting randomized methods to compute the A-optimal objective function and its gradient, and provide a detailed analysis of the error for the proposed estimators. To ensure sparse and binary design vectors, we develop a novel reweighted -minimization algorithm. We also introduce a modified A-optimal criterion and present randomized estimators for its efficient computation. We present numerical results illustrating the proposed methods on a model contaminant source identification problem, where the inverse problem seeks to recover the initial state of a contaminant plume, using discrete…
| Method | and | Precompute | Storage Cost |
|---|---|---|---|
| “Exact” | |||
| Frozen | |||
| Eig- | |||
| Randomized |
| subproblem solves | function count | function value | ||
|---|---|---|---|---|
| 57 | 10 | 395 | 39.0176 | 38 |
| 67 | 9 | 610 | 44.7219 | 30 |
| 77 | 9 | 790 | 44.5484 | 29 |
| 87 | 9 | 756 | 44.1142 | 30 |
| 127 | 9 | 1015 | 44.1139 | 30 |
| 207 | 9 | 978 | 44.1139 | 30 |
| 307 | 9 | 970 | 44.1139 | 30 |
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Randomization and reweighted -minimization for A-optimal
design of linear inverse problems††thanks: This work was funded, in part, by the National Science Foundation through the grant DMS-1745654.
Elizabeth Herman Department of Mathematics, North Carolina State University, Raleigh, NC 27695-8205, USA
Alen Alexanderian22footnotemark: 2
Arvind K. Saibaba22footnotemark: 2
Abstract
We consider optimal design of PDE-based Bayesian linear inverse problems with infinite-dimensional parameters. We focus on the A-optimal design criterion, defined as the average posterior variance and quantified by the trace of the posterior covariance operator. We propose using structure exploiting randomized methods to compute the A-optimal objective function and its gradient, and provide a detailed analysis of the error for the proposed estimators. To ensure sparse and binary design vectors, we develop a novel reweighted -minimization algorithm. We also introduce a modified A-optimal criterion and present randomized estimators for its efficient computation. We present numerical results illustrating the proposed methods on a model contaminant source identification problem, where the inverse problem seeks to recover the initial state of a contaminant plume using discrete measurements of the contaminant in space and time.
keywords:
Bayesian Inversion; A-Optimal experimental design; Large-scale ill-posed inverse problems; Randomized matrix methods; Reweighted minimization; Uncertainty quantification
{AMS}
35R30; 62K05; 68W20; 35Q62; 65C60; 62F15.
1 Introduction
A central problem in scientific computing involves estimating parameters that describe mathematical models, such as initial conditions, boundary conditions, or material parameters. This is often addressed by using experimental measurements and a mathematical model to compute estimates of unknown model parameters. In other words, one can estimate the parameters by solving an inverse problem. Experimental design involves specifying the experimental setup for collecting measurement data, with the goal of accurately recovering the parameters of interest. As such, optimal experimental design (OED) is an important aspect of effective and efficient parameter estimation. Namely, in applications where collecting experimental measurements is expensive (e.g., because of budget, labor, or physical constraints), deploying experimental resources has to be done efficiently and in a parsimonious manner. Even when collecting large amounts of data is feasible, OED is still important; the computational cost of processing all the data may be prohibitive or a poorly designed experiment with many measurements may miss important information about the parameters of interest.
To make matters concrete, we explain the inverse problem and the experimental design in the context of an application. Consider the transport of a contaminant in an urban environment or the subsurface. The forward problem involves forecasting the spread of the contaminant, whereas the inverse problem involves using the measurements of the contaminant concentration at discrete points in space and time to recover the source of the contaminant (i.e., the initial state). In this application, OED involves optimal placement of sensors, at which measurement data is collected, to reconstruct the initial state.
We focus on OED for Bayesian linear inverse problems governed by PDEs. In our formulation of the OED problem, the goal is to find an optimal subset of sensors from a fixed array of candidate sensor sites. The experimental design is parameterized by assigning non-negative weights to each candidate sensor location. Ideally, we seek a binary weight vector ; if , a sensor will be placed at the th candidate location, and if , no sensor will be placed at the location. However, formulating an optimization problem over binary weight vectors leads to a problem with combinatorial complexity that is computationally prohibitive. A common approach to address this issue is to relax the binary requirement on design weights by letting the weights take values in the interval . The sparsity of the design will then be controlled using a penalty method; see e.g., [13, 2, 23]. This results in an optimization problem of the following form:
[TABLE]
where denotes the design criterion, is a penalty parameter, and is a penalty function.
Adopting a Bayesian approach to the inverse problem, the design criterion will be a measure of the uncertainty in the estimated parameters. In this article, we focus on a popular choice known as the A-optimal criterion [22, 8]. That is, we seek a sensor configuration that results in a minimized average posterior variance. The design criterion, in this case, is given by the trace of the posterior covariance operator.
One major challenge in solving Eq. 1, specifically for PDE-based inverse problems, is the computational cost of objective function and gradient evaluations. Namely, the posterior covariance operator is dense, high-dimensional, and computationally challenging to explicitly form—computing applications of to vectors requires solving multiple PDEs. Furthermore, these computations must be performed at each iteration of an optimization algorithm used to solve Eq. 1. To address this computational challenge, efficient and accurate matrix-free approaches for computing the OED objective and its gradient are needed. Another challenge in solving the OED problem Eq. 1 is the need for a suitable penalty method that is computationally tractable and results in sparse and binary optimal weight vectors. This article is about methods for overcoming these computational challenges.
Related work
For an extensive review of the OED literature, we refer the reader to [3]. We focus here on works that are closely related to the present article. Algorithms for A-optimal designs for ill-posed linear inverse problems were proposed in [11, 13, 10] and more specifically for infinite-dimensional Bayesian linear inverse problems in [2]. In majority of these articles, Monte-Carlo trace estimators are used to approximate the A-optimal design criterion and its gradient. Also, [13, 2] advocate use of low-rank approximations using the Lanczos algorithm or the randomized SVD [14]. We refer to our previous work [20] for comparison of Monte Carlo trace estimators and those based on randomized subspace iteration; it was shown that the latter are significantly more accurate than Monte Carlo trace estimators. Regarding sparsity control, various techniques have been used to approximate the -“norm” to enforce sparse and binary designs. For example, [11, 12, 13] use the -penalty function with an appropriate threshold to post-process the solution. In [2], a continuation approach is proposed that involves solving a sequence of optimization problems with non-convex penalty functions that approximate the -“norm”. More recently, in [23], a sum-up rounding approach is proposed to obtain binary optimal designs.
Our approach and contributions
In this article, we make the following advances in methods for A-optimal sensor placements in infinite-dimensional Bayesian linear inverse problems:
We present efficient and accurate randomized estimators of the A-optimal criterion and its gradient, based on randomized subspace iteration. This is accompanied by a detailed algorithm that guides efficient implementations, discussion of computational cost, as well as theoretical error analysis; see Section 3. Our estimators are structure exploiting, in that they use the low-rank structure embedded in the posterior covariance operator. To quantify the accuracy of the estimators we present rigorous error analysis, significantly advancing the methods in [20]. A desirable feature of our analysis is that the bounds are independent of the dimension of the discretized inversion parameter. Furthermore, the computational cost (measured in the number of PDE solves) of the A-optimal objective and gradient using our proposed estimators is independent of the discretized parameter dimension. 2. 2.
We propose a new algorithm for optimal sensor placement that is based on solving a sequence of reweighted -optimization problems; see Section 4. An important benefit of this approach is that one works with convex penalty functions, and since the A-optimal criterion itself is a convex function of , in each step of the reweighted algorithm a convex optimization problem is solved. We derive this algorithm by using the Majorization-Minimization principle applied to a novel penalty function that promotes binary designs. The solution of the reweighted -optimization problems is accelerated by the efficient randomized estimators for the optimality criterion and its gradient. To our knowledge, the presented framework, based on reweighted -minimization, is the first of its kind in the context of OED. 3. 3.
Motivated by reducing computational cost, we propose a new criterion known as modified A-optimal criterion; see Section 5. This criterion is derived by considering a suitably weighted A-optimal criterion. We present randomized estimators with complete error analysis for computing the modified A-optimal criterion and its gradient.
We illustrate the benefits of the proposed algorithms on a model problem from contaminant source identification. A comprehensive set of numerical experiments is provided so to test various aspect of the presented approach; see Section 6.
Finally, we remark that the randomized estimators and the reweighted approach for promoting sparse and binary weights are of independent interest beyond the application to OED.
2 Preliminaries
In this section, we recall the background material needed in the remainder of the article.
2.1 Bayesian linear inverse problems
We consider a linear inverse problem of estimating , using the model
[TABLE]
Here is a linear parameter-to-observable map (also called the forward operator), represents the measurement noise, and is a vector of measurement data. The inversion parameter is an element of , where is a bounded domain.
The setup of the inverse problem
To fully specify the inverse problem, we need to describe the prior law of and our choice of data likelihood. For the prior, we choose a Gaussian measure . We assume the prior mean is a sufficiently regular element of and that the covariance operator is a strictly positive self-adjoint trace-class operator. Following the developments in [21, 6, 9], we use with taken to be a Laplacian-like operator [21]. This ensures that is trace-class in two and three space dimensions.
We consider the case where represents a time-dependent PDE and we assume observations are taken at sensor locations at points in time. Thus, the vector of experimental data is an element of . An application of the parameter-to-observable map, , involves a PDE solve followed by an application of a spatiotemporal observation operator. We assume a Gaussian distribution on the experimental noise, {\mathchoice{\mbox{\boldmath\displaystyle{\eta}}}{\mbox{\boldmath\textstyle{\eta}}}{\mbox{\boldmath\scriptstyle{\eta}}}{\mbox{\boldmath\scriptscriptstyle{\eta}}}}\sim\mathcal{N}({{\mathchoice{\mbox{\boldmath\displaystyle{0}}}{\mbox{\boldmath\textstyle{0}}}{\mbox{\boldmath\scriptstyle{0}}}{\mbox{\boldmath\scriptscriptstyle{0}}}}},{{\mathchoice{\mbox{\boldmath\displaystyle{\Gamma}}}{\mbox{\boldmath\textstyle{\Gamma}}}{\mbox{\boldmath\scriptstyle{\Gamma}}}{\mbox{\boldmath\scriptscriptstyle{\Gamma}}}}_{\!\mathup{noise}}}). Given this choice of the noise model—additive and Gaussian—the likelihood probability density function (pdf) is
[TABLE]
Furthermore, the solution of the Bayesian inverse problem—the posterior distribution law \mu_{\mathup{post}}^{{{\mathchoice{\mbox{\boldmath\displaystyle{y}}}{\mbox{\boldmath\textstyle{y}}}{\mbox{\boldmath\scriptstyle{y}}}{\mbox{\boldmath\scriptscriptstyle{y}}}}}}—is given by the Gaussian measure \mu_{\mathup{post}}^{{{\mathchoice{\mbox{\boldmath\displaystyle{y}}}{\mbox{\boldmath\textstyle{y}}}{\mbox{\boldmath\scriptstyle{y}}}{\mbox{\boldmath\scriptscriptstyle{y}}}}}}=\mathcal{N}({m_{\mathup{post}}^{\mathchoice{\mbox{\boldmath\displaystyle{y}}}{\mbox{\boldmath\textstyle{y}}}{\mbox{\boldmath\scriptstyle{y}}}{\mbox{\boldmath\scriptscriptstyle{y}}}}},{\Gamma_{\mathup{post}}}) with
[TABLE]
Note that here the posterior mean m_{\mathup{post}}^{\mathchoice{\mbox{\boldmath\displaystyle{y}}}{\mbox{\boldmath\textstyle{y}}}{\mbox{\boldmath\scriptstyle{y}}}{\mbox{\boldmath\scriptscriptstyle{y}}}} coincides with the maximum a posteriori probability (MAP) estimator. We refer to [21] for further details.
Discretization
We use a continuous Galerkin finite element discretization approach for the governing PDEs, as well as the inverse problem. Specifically, our discretization of the Bayesian inverse problem follows the developments in [6]. The discretized parameter space in the present case is equipped with the inner product and norm , where is the finite element mass matrix. Note that is the discretized inner product. The discretized parameter-to-observable map is a linear transformation with adjoint discussed below. The discretized prior measure \mathcal{N}({{\mathchoice{\mbox{\boldmath\displaystyle{m}}}{\mbox{\boldmath\textstyle{m}}}{\mbox{\boldmath\scriptstyle{m}}}{\mbox{\boldmath\scriptscriptstyle{m}}}}_{\mathup{pr}}},{\mathbf{\Gamma}_{\text{pr}}}) is obtained by discretizing the prior mean and covariance operator, and the discretized posterior measure is given by \mathcal{N}({{\mathchoice{\mbox{\boldmath\displaystyle{m}}}{\mbox{\boldmath\textstyle{m}}}{\mbox{\boldmath\scriptstyle{m}}}{\mbox{\boldmath\scriptscriptstyle{m}}}}_{\text{post}}},{\mathbf{\Gamma}_{\text{post}}}), with
[TABLE]
We point out that the operator \mathbf{F}^{*}{\mathchoice{\mbox{\boldmath\displaystyle{\Gamma}}}{\mbox{\boldmath\textstyle{\Gamma}}}{\mbox{\boldmath\scriptstyle{\Gamma}}}{\mbox{\boldmath\scriptscriptstyle{\Gamma}}}}^{-1}_{\text{noise}}\mathbf{F} is the Hessian of the data-misfit cost functional whose minimizer is the MAP point, and is thus referred to as the data-misfit Hessian; see e.g., [2].
It is important to keep track of the inner products in the domains and ranges of the linear mappings, appearing in the above expressions, when computing the respective adjoint operators. For the readers’ convenience, in Fig. 1, we summarize the different spaces that are important, the respective inner products, and the adjoints of the linear transformations defined between these spaces.
Using the fact that is a self-adjoint operator on and the form of the adjoint operator (see Fig. 1), we can rewrite the expression for as follows:
[TABLE]
with
[TABLE]
Note that the operator \mathbfcal{H}_{\text{m}}={\mathbfcal{F}}^{\top}{{\mathchoice{\mbox{\boldmath\displaystyle{\Gamma}}}{\mbox{\boldmath\textstyle{\Gamma}}}{\mbox{\boldmath\scriptstyle{\Gamma}}}{\mbox{\boldmath\scriptscriptstyle{\Gamma}}}}}_{\text{noise}}^{-1}{\mathbfcal{F}} is a symmetric positive semidefinite matrix, and is a similarity transform of the prior-preconditioned data-misfit Hessian \mathbf{\Gamma}^{1/2}_{\text{pr}}\mathbf{F}^{*}{\mathchoice{\mbox{\boldmath\displaystyle{\Gamma}}}{\mbox{\boldmath\textstyle{\Gamma}}}{\mbox{\boldmath\scriptstyle{\Gamma}}}{\mbox{\boldmath\scriptscriptstyle{\Gamma}}}}^{-1}_{\text{noise}}\mathbf{F}\mathbf{\Gamma}^{1/2}_{\text{pr}}. In many applications (including the application considered in Section 6), has rapidly decaying eigenvalues and therefore, it can be approximated by a low-rank matrix. This is a key insight that will be exploited in our estimators for the OED criterion and its gradient.
2.2 Randomized subspace iteration algorithm
In this article, we develop and use randomized estimators to efficiently compute the design criteria and their derivatives. We first explain how to use randomized algorithms for computing low-rank approximations of a symmetric positive semidefinite matrix . To draw connection with the previous subsection, in our application will stand for . We first draw a random Gaussian matrix (i.e., the entries are independent and identically distributed standard normal random variables). We then perform steps of subspace iteration on with the starting guess to obtain the matrix . If, for example, the matrix has rank , or the eigenvalues decay sufficiently, then the range of is a good approximation to the range of under these suitable conditions. This is the main insight behind randomized algorithms. We now show how to obtain a low-rank approximation of . A thin-QR factorization of is performed to obtain the matrix , which has orthonormal columns. We then form the “projected” matrix and obtain the low-rank approximation
[TABLE]
This low-rank approximation can be manipulated in many ways depending on the desired application. An alternative low-rank approximation can be computed using the Nyström approximation, see e.g., [14].
In addition, once the matrix is computed, it can be used in various ways. In [20], was used as an estimator for , whereas was used as estimator for . The main idea behind these estimators is that the eigenvalues of are good approximations to the eigenvalues of , when is sufficiently low-rank or has rapidly decaying eigenvalues. Our estimators for the A-optimal criterion and its gradient utilize the same idea but in a slightly different form.
2.3 A-optimal design of experiments
As mentioned in the introduction, an experimental design refers to a placement of sensors used to collect measurement data for the purposes of parameter inversion. Here we describe the basic setup of the optimization problem for finding an A-optimal design.
Experimental design and A-optimal criterion
We seek to find an optimal subset of a network of candidate sensor locations, which collect measurements at points in time. The experimental design is parameterized by a vector of design weights {\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}}\in[0,1]^{n_{s}}. In the present work, we use the A-optimal criterion to find the optimal design. That is, we seek designs that minimize the average posterior variance, as quantified by \mathsf{tr}\,\left(\mathbf{\Gamma}_{\text{post}}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}})\right). (The precise nature of the dependence of on will be explained below.) Note that \mathsf{tr}\,\left(\mathbf{\Gamma}_{\text{post}}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}})\right)=\mathsf{tr}\,\left(\mathbf{\Gamma}_{\text{post}}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}})-\mathbf{\Gamma}_{\text{pr}}\right)+\mathsf{tr}\,\left(\mathbf{\Gamma}_{\text{pr}}\right) and thus, minimizing the trace of the posterior covariance operator is equivalent to minimizing
[TABLE]
This is the objective function we seek to minimize for finding A-optimal designs. As seen below, this formulation of the A-optimal criterion is well suited for approximations via randomized matrix methods Section 2.2. Note also that minimizing \Phi_{\text{aopt}}{(}{\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}}) amounts to maximizing \mathsf{tr}\,\left(\mathbf{\Gamma}_{\text{pr}}\right)-\mathsf{tr}\,\left(\mathbf{\Gamma}_{\text{post}}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}})\right), which can be thought of as a measure of uncertainty reduction.
We can also understand Eq. 5 from a decision theoretic point of view. It is well known [8, 1, 4] that for Bayesian linear inverse problems with Gaussian prior and additive Gaussian noise models, coincides with Bayes risk (with respect to the loss function):
[TABLE]
Here denotes the discretized prior measure, \mu_{\mathup{pr}}=\mathcal{N}({\mathchoice{\mbox{\boldmath\displaystyle{{\mathchoice{\mbox{\boldmath}}{\mbox{\boldmath}}{\mbox{\boldmath}}{\mbox{\boldmath}}}{\mathup{pr}}}}}{\mbox{\boldmath\textstyle{{\mathchoice{\mbox{\boldmath}}{\mbox{\boldmath}}{\mbox{\boldmath}}{\mbox{\boldmath}}}{\mathup{pr}}}}}{\mbox{\boldmath\scriptstyle{{\mathchoice{\mbox{\boldmath}}{\mbox{\boldmath}}{\mbox{\boldmath}}{\mbox{\boldmath}}}{\mathup{pr}}}}}{\mbox{\boldmath\scriptscriptstyle{{\mathchoice{\mbox{\boldmath}}{\mbox{\boldmath}}{\mbox{\boldmath}}{\mbox{\boldmath}}}{\mathup{pr}}}}}},\mathbf{\Gamma}_{\text{pr}}). Using
[TABLE]
we see that
[TABLE]
this provides an alternate interpretation of as a Bayes risk with respect to a modified loss function.
Dependence of the A-optimal criterion to the design weights
We follow the same setup as [2]. Namely, the design weights enter the Bayesian inverse problem through the data likelihood, resulting in the operator being dependent on . To weight the spatiotemporal observations, the matrix is defined as:
[TABLE]
where is the Kronecker product. Therefore, \mathbfcal{H}_{\text{m}}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}}) is expressed as follows:
[TABLE]
We refer to [2] for details. This results in the dependent posterior covariance operator:
[TABLE]
In our formulation, we assume uncorrelated observations across sensor locations and time which implies {\mathchoice{\mbox{\boldmath\displaystyle{\Gamma}}}{\mbox{\boldmath\textstyle{\Gamma}}}{\mbox{\boldmath\scriptstyle{\Gamma}}}{\mbox{\boldmath\scriptscriptstyle{\Gamma}}}}_{\text{noise}} is a block diagonal matrix, with blocks of the form ; here , , denote the measurement noise at individual sensor locations. Using this structure for {\mathchoice{\mbox{\boldmath\displaystyle{\Gamma}}}{\mbox{\boldmath\textstyle{\Gamma}}}{\mbox{\boldmath\scriptstyle{\Gamma}}}{\mbox{\boldmath\scriptscriptstyle{\Gamma}}}}_{\text{noise}}, we define
[TABLE]
with . Thus, we have \mathbfcal{H}_{\text{m}}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}})={\mathbfcal{F}}^{\top}{\mathbf{W}^{\text{noise}}}{\mathbfcal{F}} and the A-optimal criterion can be written as
[TABLE]
with
[TABLE]
Anticipating that we will use a gradient-based solver for solving Eq. 10, we also need the gradient of \Phi_{\text{aopt}}{({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}})} which we now derive. Using Theorems B.17 and B.19 in [22], the partial derivatives of Eq. 7 with respect to , , are
[TABLE]
(We have used the notation to denote .) Note that using the definition of \mathbfcal{H}_{\text{m}}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}}), we have \partial_{j}\mathbfcal{H}_{\text{m}}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}})={\mathbfcal{F}}^{\top}\mathbf{E}_{j}^{\text{noise}}{\mathbfcal{F}}, .
The optimization problem for finding an A-optimal design
We now specialize the optimization problem Eq. 1 to the case of A-optimal sensor placement for linear inverse problems governed by time-dependent PDEs:
[TABLE]
As explained before, to enable efficient solution methods for the above optimization problem we need (i) a numerical method for fast computation of \Phi_{\text{aopt}}{(}{\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}}) and its gradient, and (ii) a choice of penalty function that promotes sparse and binary weights. The former is facilitated by the randomized subspace iteration approach outlined earlier (see Section 3), and for the latter we present an approach based on rewighted minimization (see Section 4).
3 Efficient computation of A-optimal criterion and its gradient
The computational cost of solving Eq. 10 is dominated by the PDE solves required in OED objective and gradient evaluations; these operations need to be performed repeatedly when using an optimization algorithm for solving Eq. 10. Therefore, to enable computing A-optimal designs for large-scale applications, efficient methods for objective and gradient computations are needed. In this section, we derive efficient and accurate randomized estimators for Eq. 7 and Eq. 9. The proposed estimators are matrix-free—they require only applications of the (prior-preconditioned) forward operator and its adjoint on vectors. Moreover, the computational cost of computing these estimators does not increase with the discretized parameter dimension. This is due to the fact that our estimators exploit the low-rank structure of \mathbfcal{H}_{\text{m}}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}}), a problem property that is independent of the choice of discretization.
We introduce our proposed randomized estimators for the A-optimal design criterion \Phi_{\text{aopt}}{({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}})} and its gradient \nabla\Phi_{\text{aopt}}{({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}})} in Section 3.1. Additionally, we present a detailed computational method for computing the proposed estimators. We analyze the errors associated with our proposed estimators in Section 3.2.
3.1 Randomized estimators for \Phi_{\text{aopt}}{(}{\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}}) and its gradient
Consider the low-rank approximation of \mathbfcal{H}_{\text{m}}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}}) given by \widehat{\mathbfcal{H}}_{\text{m}}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}})=\mathbf{Q}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}})\mathbf{T}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}})\mathbf{Q}^{\top}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}}), with \mathbf{Q}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}}) and \mathbf{T}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}}) computed using Algorithm 1. Replacing \mathbfcal{H}_{\text{m}}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}}) by its approximation and using the cyclic property of the trace, we obtain the estimator for the A-optimal criterion Eq. 7:
[TABLE]
where is as in Eq. 8.
To derive an estimator for the gradient, once again, we replace \mathbfcal{H}_{\text{m}}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}}) with its low-rank approximation \widehat{\mathbfcal{H}}_{\text{m}}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}}) in Eq. 9 to obtain
[TABLE]
for .
Computational procedure
First, we discuss computation of the A-optimal criterion estimator using Algorithm 1. Typically, the algorithm can be used with , due to rapid decay of eigenvalues of \mathbfcal{H}_{\text{m}}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}}). In this case, Algorithm 1 requires applications of \mathbfcal{H}_{\text{m}}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}}). Since each application of \mathbfcal{H}_{\text{m}}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}}) requires one apply (forward solve) and one apply (adjoint solve), computing \widehat{\mathbfcal{H}}_{\text{m}}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}}) requires PDE solves. Letting the spectral decomposition of the matrix \mathbf{T}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}}) be given by \mathbf{T}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}})=\mathbf{U}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}})\mathbf{\Lambda}_{\mathbf{T}}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}})\mathbf{U}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}})^{\top} and denoting \mathbf{V}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}})\equiv\mathbf{Q}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}})\mathbf{U}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}}), we have \widehat{\mathbfcal{H}}_{\text{m}}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}})=\mathbf{V}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}})\mathbf{\Lambda}_{\mathbf{T}}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}})\mathbf{V}^{\top}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}}). Applying the Sherman–Morrison–Woodbury formula [19] and the cyclic property of the trace to Eq. 11, we obtain
[TABLE]
where \mathbf{D}_{\mathbf{T}}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}})=\mathbf{\Lambda}_{\mathbf{T}}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}})(\mathbf{I}+\mathbf{\Lambda}_{\mathbf{T}}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}}))^{-1}. To simplify notation, the dependence of is suppressed for the operators used in computing the estimators for the remainder of the article; however, the notation is retained for \mathbfcal{H}_{\text{m}}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}}) and \widehat{\mathbfcal{H}}_{\text{m}}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}}).
Next, we describe computation of the gradient estimator Eq. 12. Here we assume ; the extension to the case is straightforward and is omitted. Again, using the Woodbury formula and cyclic property of the trace, we rewrite Eq. 12 as
[TABLE]
for . Expanding this expression, we obtain
[TABLE]
Note that the first term in Eq. 15 does not depend on the design , for . As a result, this term can be precomputed and used in subsequent function evaluations. We expand to be
[TABLE]
where \widehat{{\mathchoice{\mbox{\boldmath\displaystyle{e}}}{\mbox{\boldmath\textstyle{e}}}{\mbox{\boldmath\scriptstyle{e}}}{\mbox{\boldmath\scriptscriptstyle{e}}}}}_{k} is the column of the identity matrix of size . Because there are only columns of with nonzero entries, the total cost to precompute for is PDE solves.
To compute the remaining terms in Eq. 15, we exploit the fact that has columns. Notice all the other occurrences of and occur as a combination of and (or of their transposes). Both of these terms require PDE solves to compute. As a result, the total cost to evaluate \widehat{\Phi}_{\text{aopt}}{({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}})} and \widehat{\nabla\Phi}_{\text{aopt}}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}};\ell) is PDE solves to apply Algorithm 1 and PDE solves to compute and . We detail the steps for computing our estimators for A-optimal criterion and its gradient in Algorithm 2.
Alternative approaches and summary of computational cost
A closely related variation of Algorithm 2 is obtained by replacing step 1 of the algorithm (i.e., randomized subspace iteration) by the solution of an eigenvalue problem to compute the dominant eigenvalues of \mathbfcal{H}_{\text{m}}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}}) “exactly”. We refer to this method as Eig-, where is the target rank of \mathbfcal{H}_{\text{m}}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}}). This idea was explored for computing Bayesian D-optimal designs in [3]. The resulting cost is similar to that of the randomized method: it would cost PDE solves per iteration to compute the spectral decomposition of \mathbfcal{H}_{\text{m}}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}}), plus PDE solves to precompute for . While both the randomized and Eig- methods provide a viable scheme for computing the A-optimal criteria, our randomized method can exploit parallelism to lower computational costs. Each matrix-vector application with \mathbfcal{H}_{\text{m}}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}}) in Algorithm 1 can be computed in parallel. However, if accurate eigenpairs of \mathbfcal{H}_{\text{m}}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}}) are of importance to the problem, one can choose to use the Eig- approach at the cost of computing a more challenging problem.
Another possibility suitable for problems where the forward model does not depend on the design (as is the case in the present work), is to precompute a low-rank SVD of , which can then be applied as necessary to compute the A-optimal criterion and its gradient. This frozen forward operator approach has been explored in [2, 13] for the A-optimal criterion and in [3] for the D-optimal criterion. The resulting PDE cost of precomputing a low-rank approximation of is , with indicating the target rank. The Frozen method is beneficial as no additional PDE solves are required when applying Algorithm 2; however, this approach would not favor problems where depends on , nor can the modeling errors associated with the PDE be controlled in subsequent evaluations without the construction of another operator.
Finally, if the problem size is not too large (i.e., in small-scale applications), one could explicitly construct the forward operator . This enables exact (excluding floating point errors) computation of the objective and its gradient. The total cost for evaluating involves an upfront cost of PDE solves. We summarize the computational cost of Algorithm 2 along with the other alternatives mentioned above in Table 1.
To summarize, the randomized methods for computing OED objective and its gradient present several attractive features; our approach is well suited to large-scale applications; it is matrix free, simple to implement and parallelize, and exploits low-rank structure in the inverse problem. Moreover, as is the case with the Eig-k approach, the performance of the randomized subspace iteration does not degrade as the dimension of the parameter increases due to mesh refinement. This is the case because the randomized subspace iteration relies on spectral properties of the prior-preconditioned data misfit Hessian—a problem structure that is independent of discretization.
3.2 Error Analysis
Here we analyze the error associated with our estimators computed using Algorithm 2. For fixed {\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}}\in[0,1]^{n_{s}}, since \mathbfcal{H}_{\text{m}}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}})\in\mathbb{R}^{n\times n} is symmetric positive semidefinite, we can order its eigenvalues as . Suppose are the dominant eigenvalues of \mathbfcal{H}_{\text{m}}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}}), we define and and we assume that the eigenvalue ratio satisfies
[TABLE]
We now present the error bounds for the objective function and its gradient. To this end, we define the constant as
[TABLE]
with r=\text{rank}(\mathbfcal{H}_{\text{m}}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}})) and .
Theorem 3.1**.**
Let \widehat{\Phi}_{\text{aopt}}{({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}};\ell)} and \widehat{\nabla\Phi}_{\text{aopt}}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}};\ell) be approximations of the A-optimal objective function \Phi_{\text{aopt}}{({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}})} and its gradient \nabla\Phi_{\text{aopt}}{({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}})}, respectively, computed using Algorithm 2 for fixed {\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}}\in[0,1]^{n_{s}}. Recall that is the target rank of \mathbfcal{H}_{\text{m}}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}}), is the oversampling parameter such that , and is the number of subspace iterations. Assume that . Then, with
[TABLE]
Furthermore, with , for ,
[TABLE]
Proof 3.2**.**
*See Section A.2. *
In Theorem 3.1, the estimators are unbiased when the target rank equals the rank of \mathbfcal{H}_{\text{m}}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}}). If the eigenvalues decay rapidly, the bounds suggest that the estimators are accurate. Recall that \text{rank}(\mathbfcal{H}_{\text{m}}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}}))\leq\min\{n_{s}n_{t},n\} is the number of nonzero eigenvalues. Consequently, it is seen that the bounds are independent of the dimension of the discretization .
4 An optimization framework for finding binary designs
We seek A-optimal designs by solving an optimization problem of the form,
[TABLE]
where the design criterion is either the A-optimal criterion \Phi_{\text{aopt}}{(}{\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}}) or the modified A-optimal criterion \Phi_{\text{mod}}{(}{\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}}) (see Section 5). In the previous sections, we laid out an efficient framework for computing accurate approximations to the A-optimal criterion and its gradient. We now discuss the choice of the penalty term and the algorithm for solving the optimization problem.
The choice of the penalty term must satisfy two conditions: enforcing sparsity, measured by the number of nonzeros of the design vector, and binary designs, i.e., designs vectors whose entries are either or [math]. One possibility for the penalty function is the -“norm”, P_{\ell_{0}}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}})=\|{\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}}\|_{0}, which measures the number of nonzero entries in the design. However, the resulting optimization problem is challenging to solve due to its combinatorial complexity. A common practice is to replace the -“norm” penalty by the -norm, P_{\ell_{1}}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}})=\|{\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}}\|_{1}. The penalty function has desirable features: it is a convex penalty function that promotes sparsity of the optimal design vector . However, the resulting design is sparse but not necessarily binary and additional post-processing in the form of thresholding is necessary to enforce binary designs.
In what follows, we introduce a suitable penalty function that enforces both sparsity and binary designs and an algorithm for solving the OED optimization problem based on the MM approach. The resulting algorithm takes the form of a reweighted -minimization algorithm [7].
4.1 Penalty functions
We propose the following penalty function
[TABLE]
for a user-defined parameter . This penalty function approximates for small values of ; however, as becomes smaller the corresponding optimization problem becomes harder. To illustrate the choice of penalty functions, in Fig. 2, we plot along with and , with . Using in the OED problem leads to the optimization problem,
[TABLE]
In Eq. 21, the absolute values in definition of can be dropped since we limit the search for optimal solutions in . Since P_{\epsilon}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}}) is concave, Eq. 21 is a non-convex optimization problem. To tackle this, we adopt the majorization-minimization (MM) approach.
4.2 MM approach and reweighted algorithm
The idea behind the MM approach is to solve a sequence of optimization problems whose solutions converge to that of the original problem [16, 17]. This sequence is generated by a carefully constructed surrogate that satisfies two properties—the surrogate must majorize the objective function for all values, and the surrogate must match the objective function at the current iterate. More specifically, suppose
[TABLE]
Then the surrogate function g({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}}|{\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}}^{(m)}) at the current iterate {\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}}^{(m)} must satisfy
[TABLE]
Granted the existence of this surrogate function, to find the next iterate {\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}}^{(m+1)} we solve the optimization problem
[TABLE]
To show the objective function decreases at the next iterate, observe that the next iterate {\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}}^{(m+1)} stays within the feasible region and use the two properties of the surrogate function
[TABLE]
To construct this surrogate function, we use the fact that a concave function is below its tangent [17, Equation (4.7)]. Applying this to our concave penalty P_{\epsilon}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}}), we have
[TABLE]
With this majorization relation, we define the surrogate function to be
[TABLE]
By dropping the terms that do not depend on , it can be readily verified that Eq. 22 can be replaced by the equivalent problem
[TABLE]
where {\mathchoice{\mbox{\boldmath\displaystyle{R}}}{\mbox{\boldmath\textstyle{R}}}{\mbox{\boldmath\scriptstyle{R}}}{\mbox{\boldmath\scriptscriptstyle{R}}}}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}})=\text{diag}\left(\frac{\epsilon}{(w_{1}+\epsilon)^{2}},\dots,\frac{\epsilon}{(w_{n_{s}}+\epsilon)^{2}}\right). We see that Eq. 23 is of the form of a reweighted -optimization problem. The details of the optimization procedure are given in Algorithm 3. (We remark that, in Algorithm 3, other metrics measuring the difference between the successive weight vectors can be used in step 3.)
We conclude this section with a few remarks regarding this algorithm. In our application, \Phi({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}}) is convex; therefore, each subproblem to update the design weights is also convex. To initialize the reweighted algorithm, we start with the weights , . This ensures that, in the first step, we are computing the solution of the -penalized optimization problem. The subsequent reweighted iterations further promote binary designs. To solve the subproblems in each reweighted iteration we use an interior point algorithm; however, any solver for appropriate convex optimization may be used. In Section 6, we will provide a discussion of our choice of for our application. It is also worth mentioning that besides the penalty function used above, another possible choice is which yields the weights .
5 Modifed A-optimal criterion
Motivated by reducing the computational cost of computing A-optimal designs, in this section, we introduce a modified A-optimality criterion. As mentioned in [8], we can consider a weighted A-optimal criterion \mathsf{tr}\,\left(\boldsymbol{\Gamma}\mathbf{\Gamma}_{\text{post}}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}})\right), where is a positive semidefinite weighting matrix. Similar to Eq. 5, we work with , since the term is independent of the weights . By choosing , we obtain the modified A-optimal criterion
[TABLE]
Note that the expression for \Phi_{\text{mod}}{({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}})} remains meaningful in the infinite-dimensional limit. This can be seen by noting that
[TABLE]
and using the fact that in the infinite-dimensional limit, for every {\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}}\in[0,1]^{n_{s}}, \mathbfcal{H}_{\text{m}}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}}) is trace class and \big{(}\mathbf{I}+\mathbfcal{H}_{\text{m}}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}})\big{)}^{-1} is a bounded linear operator.
We show in our numerical results that the modified A-optimality criterion can be useful in practice, if a cheaper alternative to the A-optimal criterion is desired, and minimizing can provide designs that lead to small posterior uncertainty.
5.1 Derivation of estimators
Here we seek to improve the efficiency of the modified A-optimal criterion by computing a randomized estimator for the modified A-optimal criterion and its gradient. As in the previous derivation of our estimators, we replace \mathbfcal{H}_{\text{m}}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}}) by its low rank approximation to obtain the randomized estimator for the modified A-optimal criterion
[TABLE]
Similarly, we use the same low-rank approximation {\widehat{\mathbfcal{H}}_{\text{m}}}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}}) in the gradient of the modified A-optimal criterion to obtain the randomized estimator
[TABLE]
for .
5.2 Computational procedure and cost
Using similar techniques as those described in Section 3.1, we can write the estimator for the modified A-optimal criterion in terms of the eigenvalues of :
[TABLE]
Moreover,
[TABLE]
with and defined as in Section 3.1.
The procedure for computing the estimators for the modified A-optimal criterion follows the steps in Algorithm 2 closely. Instead of presenting an additional algorithm, we provide an overview of the computation of the estimators for the modified A-optimal criterion along with the associated computational cost in terms of the number of PDE solves.
To evaluate the estimators for the modified A-optimal criterion, the only precomputation we perform is to obtain
[TABLE]
This term appears in the estimator for the gradient and accumulates a total cost of PDE solves. The remaining terms in the estimators depend on a design and, in particular, the eigenvalues and eigenvectors of {\widehat{\mathbfcal{H}}_{\text{m}}}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}}). As with the estimators for the A-optimal criterion, the eigenvalues and eigenvectors of \widehat{\mathbfcal{H}}_{\text{m}}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}}) are obtained by Algorithm 1. Recall that the cost associated with Algorithm 1, with , is PDE solves. Once the eigenvalues and eigenvectors are computed, Eq. 25 can be evaluated without any additional PDE solves. The remaining computational effort occurs in the evaluation of the gradient. Because of our modification to the A-optimal criterion, the expression Eq. 28 is efficiently evaluated by computing . Therefore, the total cost of evaluating the estimators for the modified A-optimal criterion and its gradient is PDE solves. From this computational cost analysis, we see the modified A-optimal estimators require less PDE solves than the A-optimal estimators.
5.3 Error analysis
We now quantify the absolute error of our estimators with the following theorem.
Theorem 5.1**.**
Let \widehat{\Phi}_{\text{mod}}{({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}};\ell)} and \widehat{\nabla\Phi}_{\text{mod}}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}};\ell) be the randomized estimators approximating the modified A-optimal objective function \Phi_{\text{mod}}{({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}})} and its gradient \nabla\Phi_{\text{mod}}{({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}})}, respectively. Using the notation and assumptions of Theorem 3.1, for fixed {\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}}\in[0,1]^{n_{s}}
[TABLE]
and for ,
[TABLE]
*where is defined in Eq. 16.
Proof 5.2**.**
*See Appendix A.3. *
Notice the bounds presented in Theorem 3.1 and 5.1 differ by a factor of . Since the modified A-optimal criterion removes one application of from the computation of the A-optimal criterion, the bounds related to the modified A-optimal criterion no longer have the factor of .
6 Numerical results
In this section, we present numerical results that test various aspects of the proposed methods. We begin by a brief description of the inverse advection-diffusion problem used to illustrate the proposed OED methods, in Section 6.1. The setup of our model problem is adapted from that in [2], where further details about the forward and inverse problem can be found. In Section 6.1, we also describle the numerical methods used for solving the forward problem, as well as the optimization solver for the OED problem. In Section 6.2, we test the accuracy of our randomized estimators and illustrate our error bounds. Then in Section 6.3, we investigate the performance of our proposed reweighted -optimization approach. Next, we utilize the proposed optimization framework to compute A-optimal designs in Section 6.4. Finally, in Section 6.5, we compare A-optimal sensor placements with those computed by minimizing the modified A-optimal criterion.
6.1 Model problem and solvers
We consider a two-dimensional time-dependent advection-diffusion equation
[TABLE]
which models the transport of contaminants (e.g., in the atmosphere or the subsurface). Here is the diffusion coefficient and is taken to be . The velocity field is computed by solving a steady-state Navier-Stokes equations, as in [2]. The domain , depicted in Fig. 3, is the unit square in with the gray rectangles, modeling obstacles/buildings, removed. The boundary is the union of the outer boundary and the boundaries of the obstacles. The PDE is discretized using linear triangular continuous Galerkin finite elements in space and implicit Euler in time. We let the final simulation time be .
The inverse problem involves reconstructing the initial state from space-time point measurements of u({\mathchoice{\mbox{\boldmath\displaystyle{x}}}{\mbox{\boldmath\textstyle{x}}}{\mbox{\boldmath\scriptstyle{x}}}{\mbox{\boldmath\scriptscriptstyle{x}}}},t). We consider sensor candidate locations distributed throughout the domain which are indicated by hollow squares in Fig. 3. Measurement data is collected from a subset of these locations at three observation times— and . To simulate noisy observations, 2% noise is added to the simulated data.
Recall, in Section 2, we define our prior covariance operator to be a Laplacian-like operator. Following [2], an application of the square root of the prior on is , which satisfies the following weak form
[TABLE]
Here and control the variance and correlation length and are chosen to be and , respectively.
The optimization solver used for the OED problem is a quasi-Newton interior point method. Specifically, to solve each subproblem of Algorithm 3, we use MATLAB’s interior point solver provided by the fmincon function; BFGS approximation to the Hessian is used for line search. We use a vector of all ones, {\mathchoice{\mbox{\boldmath\displaystyle{1}}}{\mbox{\boldmath\textstyle{1}}}{\mbox{\boldmath\scriptstyle{1}}}{\mbox{\boldmath\scriptscriptstyle{1}}}}\in\mathbb{R}^{n_{s}}, as the initial guess for the optimization solver.
For the numerical experiments presented in this section, the random matrix in Algorithm 1 is fixed during the optimization process. However, because of the randomness, the accuracy of the estimators, and the optimal design thus obtained, may vary with different realizations of . By conducting additional numerical experiments (not reported here) the stochastic nature of the estimators resulted in minor variability (at most one or two sensor locations) in the optimal experimental designs for a modest value of . On the other hand, if is sufficiently large, we observed that the same design was obtained with different realizations of
6.2 Accuracy of estimators
Here we examine the accuracy of the randomized estimators with respect to , the number of columns in the sampling matrix in the randomized subspace iteration algorithm. Specifically, we compute
[TABLE]
with taken to be a vector of all ones; that is, with all sensors activated. We let to vary from to , because the rank of \mathbfcal{H}_{\text{m}}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}}) is no larger than the number of observations taken, . Fig. 4 illustrates the relative error in the estimators for the A-optimal criterion and its gradient (left) and the modified A-optimal criterion and its gradient (right), as is varied. We observe that the error decreases rapidly with increasing . This illustrates accuracy and efficiency of our estimators.
Next, we consider the absolute error in the estimators for the objective function and compare them with the theoretical bounds derived in Theorems 3.1 and 5.1. In Fig. 5, we compare the absolute error in the estimators with bounds from Theorem 3.1 and Theorem 5.1. As before, we take {\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}}=[1,1,\dots,1]^{\top}\in\mathbb{R}^{n_{s}}. We observe that our error bound captures the general trend in the error. Moreover, the error bound for the modified A-optimality is better since it does not have the additional factor of .
6.3 Performance of the reweighted algorithm
We now consider solving Eq. 10 with Algorithm 2 and 3. We first consider the choice of the value of in Eq. 23. The user-defined parameter controls the steepness of the penalty function at the origin; see Fig. 2. However, we observed that if the penalty function is too steep, the optimization solvers took more iterations without substantially altering the optimal designs. We found to be sufficiently small in our numerical experiments, and we keep this fixed for the remainder of the numerical experiments.
We perform two numerical experiments examining the impact of changing and the penalty parameter (in Eq. 10). In the first experiment, we fix the penalty coefficient at an experimentally determined value of , and vary ; the results are recorded in Table 2. Notice when , the objective function value evaluated with the optimal solution, the number of active sensors, and number of subproblem solves do not change. This suggests that the randomized estimators are sufficiently accurate with and this yields a substantial reduction in computational cost.
The second experiment involves varying , which indirectly controls the number of sensors in computed designs, with kept fixed. Here we fix , which corresponds to an accuracy on the order of for the A-optimal criterion (cf. Fig. 4). In Fig. 6, we report the design weights sorted in descending order, as varies. This shows that the reweighted algorithm indeed produces binary designs for a range of penalty parameters. We also notice in Fig. 7 that as increases the sparsity increases (i.e., the number of active sensors decreases) and the number of function evaluations increases. The right panel compares the cost of solving an -penalized problem for the corresponding penalty parameter . Since this problem is the first iterate of the reweighted algorithm we see that an additional cost is required to obtain binary designs and this cost increases with increasing (i.e., more sparse designs).
6.4 Computing optimal designs
In Fig. 8 (left), we report an A-optimal sensor placement obtained using our optimization framework, with and ; the resulting optimal sensor locations, with active sensors, are superimposed on the posterior standard deviation field. While the design is computed to yield a minimal average variance of the posterior distribution, it is also important to consider the mean of this distribution. For completeness, in Fig. 8 we show the “true“ initial condition (middle panel), used to generate synthetic data, and the mean of the resulting posterior distribution for the active sensor design (right panel).
With the sensor design, we also illustrate the resulting uncertainty reduction by looking at the prior and posterior standard deviation fields; see Fig. 9.
We now compare the designs obtained using the reweighted algorithm and our estimators against designs chosen at random, illustrating the effectiveness of the proposed A-optimal design strategy. Recall that varying allows us to obtain optimal designs with different numbers of active sensors. For each value of , we use Algorithm 2 and 3 to compute an optimal design. We then draw random designs with the same number of active sensors as the optimal design obtained using our algorithms. To enable a consistent comparison, we evaluate the exact A-optimal criterion \Phi_{\text{aopt}}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}}) at the computed optimal designs and the random ones; the results are reported in the left panel of Fig. 10. The values corresponding to the computed optimal designs are indicated as dots on the black solid line. The values obtained from the random designs are indicated by the squares. We note that the designs computed with the reweighted algorithm consistently beat the random designs, as expected. This observation is emphasized in the right panel of Fig. 10. Here, we compared the computed optimal design (when ) with 1500 randomly generated designs using the exact trace of the posterior covariance. Again, the computed design results in a lower true A-optimal value than the random designs.
6.5 Comparing A-optimal and modified A-optimal designs
Here we provide a quantitative comparison of sensor placements obtained by minimizing A-optimal and modified A-optimal criteria using our proposed algorithms. Specifically, for various values of , we solve Eq. 10 with both the A-optimal and modified A-optimal estimators to obtain two sets of designs. By varying , the resulting designs obtained with the A-optimal and modified A-optimal estimators have different number of active sensors. Using both sets of designs, we evaluate the exact A-optimal criterion \Phi_{\text{aopt}}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}}); these are displayed in Fig. 11. Observe that in all cases the computed A-optimal and modified A-optimal designs lead to similar levels of average posterior variance. This suggests that the modified A-optimal criterion could be used as a surrogate for the A-optimal criterion. Using the modified A-optimal criterion decreases the overall number of PDE solves and yields designs that result in values of the average posterior variance close to those produced by A-optimal designs.
7 Conclusion
We have established an efficient and flexible computational framework for A-optimal design of experiments in large-scale Bayesian linear inverse problems. The proposed randomized estimators for the OED objective and its gradient are accurate, efficient, simple to implement and parallelize. Specifically, the randomized estimators exploit the low-rank structure in the inverse problem; namely, the low-rank structure of the prior-preconditioned data misfit Hessian—a common feature of ill-posed inverse problems. Our reweighted -minimization strategy is tailored to sensor placement problems, where finding binary optimal design vectors is desirable. We also presented the modified A-optimal criterion, which is more computationally efficient to compute and can provide designs that, while sub-optimal if the goal is to compute A-optimal designs, provide a systematic means for obtaining sensor placements with small posterior uncertainty levels.
Open questions that we seek to explore in our future work include adaptive determination of the target rank within the optimization algorithm, to further reduce computational costs, while ensuring sufficiently accurate estimates of the OED objective and gradient. Another possible line of inquiry is to use different low-rank approximations, such as Nyström’s method, and extending the randomized estimators to approximate trace of matrix functions. We also seek to incorporate the randomized estimators in a suitable optimization framework for Bayesian nonlinear inverse problems, in our future work.
Acknowledgments
We would like to thank Eric Chi for useful discussions regarding the MM algorithm. This material is based upon work supported in part by the National Science Foundation (NSF) award DMS-1745654.
Appendix A Proofs of bounds
A.1 Trace of matrix function
In the proofs below, we use the Loewner partial ordering [15, Chapter 7.7]; we briefly recapitulate some main results that will be useful in our proof. Let be symmetric positive definite; then means that is positive semidefinite. For any , it also follows that . Let be the eigendecomposition of . Then, and . If is monotonically increasing then since implies for .
The following bound allows us to bound the trace of a matrix function in terms of its diagonal subblocks.
Lemma 1**.**
Let
[TABLE]
be a symmetric positive definite matrix. Let be a nonnegative concave function on . Then
[TABLE]
Proof A.1**.**
See Theorem 2.1 and Remark 2.4 in [18].
We are ready to state and prove our main result of this section, which is the key in proving Theorem 3.1 and 5.1.
Theorem 2**.**
Let be a symmetric positive definite matrix with eigendecomposition
[TABLE]
where and contain the eigenvalues arranged in descending order. Assume that the eigenvalue ratio . Let be the target rank, be the oversampling parameter such that , and let be the number of subspace iterations. Furthermore, assume that and are computed using Algorithm 1 and define . Then
[TABLE]
*where , and the constant is defined in Eq. 16. *
Proof A.2**.**
Suppose . Then, has at most nonzero eigenvalues, and thus, we can define so that
[TABLE]
We split this proof into several steps.
Step 0: Lower bound
Let be the eigenvalues of (and also ). By the Cauchy interlacing theorem (see **[20, Lemma 1]** for the specific version of the argument), for . Using properties of the trace operator
[TABLE]
Since each term in the summation is nonnegative, the lower bound follows.
Step 1. Trace of matrix function
We first write , where is an orthogonal projection matrix onto the range of . Since has the same eigenvalues as **[15, Theorem 1.3.22]**,
[TABLE]
Also, since , it follows that and . Therefore, from the proof of **[5, Lemma X.1.4]** and Eq. 31,
[TABLE]
where was defined in the statement of the theorem.
Step 2. Reducing the dimensionality
Let . In Algorithm 1, we compute and let be the thin QR factorization of . Let be defined as in the proof of **[20, Theorem 6]**. It was also shown that has orthonormal columns, so that
[TABLE]
has orthonormal columns. The following sequence of identities also hold: , , and
[TABLE]
Since is a monotonic increasing function, .
Step 3. Split into the diagonal blocks
We can rewrite as
[TABLE]
where represents blocks that can be ignored and
[TABLE]
We can invoke 1, since is concave and nonnegative on . Therefore, we have
[TABLE]
Note that the matrix disappears, because the trace is unitarily invariant.
Step 4. Completing the structural bound
Using an SVD based argument it can be shown that , so that
[TABLE]
Therefore, since is monotonically increasing
[TABLE]
Note that is has at most nonzero singular values. Looking more into the structure of , and using Eq. 30, we can write
[TABLE]
where such that . Using the multiplicative singular value inequalities **[15, Equation (7.3.14)]** and repeated use of the submultiplicative inequality gives
[TABLE]
The analysis splits into two cases:
Case 1: .
Since is monotonically increasing, using Eqs. 32 and 33
[TABLE]
Case 2: .
Since has at most nonzero eigenvalues, use the fact that , along with Eqs. 32 and 33 to obtain
[TABLE]
To summarize, in both cases Similarly, since , we can show
[TABLE]
so that . Combine with step 3 to obtain
[TABLE]
Combine this with the results of steps 1 and 2, to obtain the structural bound
[TABLE]
Step 5. The expectation bound
Note that and . From the proof of **[20, Theorem 1]**, we have , where was defined in Eq. 16. By Jensen’s inequality, using the fact that is concave on we have
[TABLE]
*Combining this with the lower bound (step 0) completes the proof. *
A.2 Proof of Theorem 3.1
For the remaining discussion, recall the notation from Theorem 3.1
[TABLE]
where and are defined in Eq. 3 and Eq. 6, respectively. We will also need
Lemma 3** (See [3]).**
*Let and let be a symmetric positive semidefinite matrix. Then, we have .
Proof A.3** ( Theorem 3.1).**
Recall our estimator \widehat{\Phi}_{\text{aopt}}{({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}};\ell)} from Eq. 11. For fixed , using 3 we have
[TABLE]
Applying 2 establishes Eq. 17.
Next, we consider Eq. 18. Recall the estimator \widehat{\partial_{j}\Phi}_{\text{aopt}}{({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}};\ell)} from Eq. 12. We can write the absolute error as
[TABLE]
where \widehat{\mathbfcal{H}}_{\text{m}}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}})=\mathbf{Q}\mathbf{T}\mathbf{Q}^{\top}. We use the decomposition
[TABLE]
where \mathbf{D}\equiv(\mathbf{I}+\widehat{\mathbfcal{H}}_{\text{m}}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}}))^{-1}-(\mathbf{I}+\mathbfcal{H}_{\text{m}}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}}))^{-1}. Repeated application of 3 gives
[TABLE]
*Since \mathbf{I}+\mathbfcal{H}_{\text{m}}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}}) and \mathbf{I}+\widehat{\mathbfcal{H}}_{\text{m}}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}}) have eigenvalues greater than or equal to one, \|(\mathbf{I}+\mathbfcal{H}_{\text{m}}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}}))^{-1}\|_{2}+\|(\mathbf{I}+\widehat{\mathbfcal{H}}_{\text{m}}({\mathchoice{\mbox{\boldmath\displaystyle{w}}}{\mbox{\boldmath\textstyle{w}}}{\mbox{\boldmath\scriptstyle{w}}}{\mbox{\boldmath\scriptscriptstyle{w}}}}))^{-1}\|_{2}\leq 2. Finally, taking the expectation and applying 2, we have the desired result. *
A.3 Proof of Theorem 5.1
Proof A.4**.**
*The proof follows in similar lines as the proof of Theorem 3.1 except the fact that in Theorem 5.1 we do not have in the expressions. *
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1[1] A. Alexanderian, P. J. Gloor, O. Ghattas, et al. , On Bayesian A-and D-optimal experimental designs in infinite dimensions , Bayesian Analysis, 11 (2016), pp. 671–695.
- 2[2] A. Alexanderian, N. Petra, G. Stadler, and O. Ghattas , A-optimal design of experiments for infinite-dimensional Bayesian linear inverse problems with regularized ℓ 0 subscript ℓ 0 \ell_{0} -sparsification , SIAM Journal on Scientific Computing, 36 (2014), pp. pp. A 2122 – A 2148.
- 3[3] A. Alexanderian and A. K. Saibaba , Efficient D-optimal design of experiments for infinite-dimensional Bayesian linear inverse problems , SIAM Journal on Scientific Computing, (2018).
- 4[4] 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.
- 5[5] R. Bhatia , Matrix analysis , vol. 169 of Graduate Texts in Mathematics, Springer-Verlag, New York, 1997.
- 6[6] T. Bui-Thanh, O. Ghattas, J. Martin, and G. Stadler , A computational framework for infinite-dimensional Bayesian inverse problems Part I: The linearized case, with application to global seismic inversion , SIAM Journal on Scientific Computing, 35 (2013), pp. A 2494–A 2523.
- 7[7] E. J. Candès, M. B. Wakin, and S. P. Boyd , Enhancing sparsity by reweighted ℓ 1 subscript ℓ 1 \ell_{1} minimization , Journal of Fourier Analysis and Applications, 14 (2008), pp. 877–905.
- 8[8] K. Chaloner and I. Verdinelli , Bayesian experimental design: A review , Statistical Science, (1995), pp. 273–304.