\chapter{Introduction}

This thesis develops numerical methods for optimal experimental designs for large-scale ill-posed inverse problems. 
Optimal experimental design has broad applications in various scientific areas such as geophysics, medical imaging and biology
 \cite{ChalonerVerdinelli1995, Curtis1999, HaberHoreshTenorio08}.
In recent years, the study of optimal experimental designs has become far more popular due to these and related applications
across all the natural and social sciences.

The term `experiment' is defined to be a random process that is conducted to support or disprove a hypothesis. It can be conducted to study
almost any object such as people, animals, materials, etc. For example, rolling a die to test whether or not the die is fair
is a very simple experiment. Other examples include taking an MRI for medical purposes \cite{Mine},
shooting a ray and studying the effects as it passes through certain media \cite{AjoFranklin2008} and even clinical trials 
for testing the effects of medications \cite{ChalonerVerdinelli1995}. 
 
An optimal experimental design is the design of an experiment, which is optimal with respect to some statistical criterion. 
As a specific example, researchers at OPTEC have been studying the optimal characteristics and positions for dampers used to 
minimize vibrations in footbridges. Unfortunately, the dampers are expensive and the optimal placement is very difficult to determine
due to the random nature of the load and patterns in foot traffic. Hence, they would like to design an experiment to determine 
what types of dampers to use and where to place them. Such optimally placed dampers would provide a balance between the chance of 
damage to the bridge and the financial cost.

We introduce the theory of optimal experimental design via studying the solution of a linear ill-posed inverse problem of the form
\begin{equation}
Am + \epsilon = b,
\end{equation} 
where $A$ is a discretization of some linear operator, $b$ is the data we have observed and $\epsilon$ is the noise contained in 
the data. Since the problem is ill-posed, there is no hope to recover the model $m$ directly from the observed data. Hence, we 
include the following two elements and solve this optimization problem instead:
\begin{equation}
\hat{m} =  \underset{m}{\operatorname{arg\,min}}\,(A m - b)^{\top}W(Am - b) + R(L, m).
\end{equation}

The first element we add is a diagonal matrix $W = {\rm diag}\{\,w\,\}$, in which each entry on the diagonal stands for the relative frequency that 
each experiment is chosen. Suppose the rows of the matrix $A$ represent $n$ different experiments. In order to obtain the best 
possible solution, we need to use all of the $n$ experiments. However, in most cases, the computation for this kind of problem 
is very expensive and time consuming, especially when $n$ is large. Also, some experiments contribute more to the solution than 
others. So, putting more weight on those entries in $w$ that correspond to the important experiments would be a good idea. 
On the other hand, if the contributions of some experiments are so small that leaving them out would not harm the result much, 
then we might consider not conducting those at all. This can be done by setting the corresponding $w$ entries to be $0$.
 
The other element we consider is a regularization functional that adds additional information to the problem in order to 
introduce stability or to incorporate some a priori information about the desired solution. 

It is obvious that, in order to 
obtain a good estimated solution, we would like to find the optimal $W$ or $L$. Therefore, there are two questions we want to 
answer throughout this thesis: what is the best weight matrix $W$? What is the best regularization matrix $L$?

With this thought in mind, the two goals of this thesis are to develop numerical methods for solving large-scale optimal 
experimental design problems efficiently and to apply optimal experimental design ideas to applications in regularization 
techniques and geophysics.

\section{Bayesian optimal experimental design}

Design of experiments leads to specifying all aspects of an experiment, in other words, controlling the values of variables
that are used to describe the experiment. This kind of control may include: choosing the set of experiments to study, deciding 
the sample size and determining a certain time length to perform the experiment.
  
Decisions need to be made before gathering information when designing experiments. In general, some information
is usually available before the experiments happen. For example, one usually has prior experience with the type of 
experiments to be conducted. One may even have collections of training models that provide information about the objects 
one can expect to recover. Thus, a Bayesian approach \cite{ChalonerVerdinelli1995} seems ideal for experimental design, where
such prior information is provided in terms of probability distribution functions.

Bayesian experimental design is based on using Bayesian inference to interpret the observed data from the experiment.
In Bayesian inference, observations are used to calculate the probability that 
a hypothesis is true. It is called `Bayesian' because of the use of the Bayes' theorem in the calculation process. 
As a consequence, the goal of the Bayesian design of an experiment is to obtain a high
probability of reaching a correct conclusion before conducting any experiment. 

We use Bayesian inference as opposed to frequentist inference, which uses only the probability of the observed data, thereby taking no
prior probability of the hypothesis into account. The usage of prior information is the major difference between 
the two ways of inference. 

Through the Bayesian approach, both prior knowledge of the parameters to be determined and uncertainties in the observations 
are taken into account. Hence, optimal decisions can be made 
based on the study of uncertainty, which can be controlled by properly choosing the values of the random variables of interest.
In order to reduce this uncertainty, it is desirable to obtain estimates of the variables with small variance.

\subsection{A brief review of Bayesian experimental designs}

The goal of designing an experiment is to maximize the expected utility that is chosen to reflect purposes of the experiment. 
The utility is usually defined by measuring the accuracy of the information provided by the experiments.
Various experimental design optimality criteria have been developed based on different choices of the utility function.

It is well known that the least squares estimator
minimizes the variance of the unbiased estimators. For single variable models, the reciprocal of the variance of an estimator 
is called the `Fisher information' of the estimator. Hence, minimizing the variance corresponds to maximizing this
`Fisher information'. However, when the statistical model consists of multiple variables, the mean of the estimator is a vector 
and its variance becomes a matrix. The inverse of this variance matrix is called the `information matrix'. Thus, minimizing the 
variance is equivalent to `maximizing' the `information matrix'.

The difficulty of maximizing the information matrix is: how do we determine whether or not a matrix is maximized? 
Various optimality criteria have been developed to measure the ‘largeness’ of this information matrix using statistical theory.
In practice, the information matrix has been compressed into different real-valued functions so that it can be easily maximized. 
Some popular criteria are Bayesian $A$-, $C$-, $D$- and $E$-optimality. This class of criteria is called Bayesian alphabetical optimality. 

\subsubsection*{Bayesian $D$-optimality}
Stone, DeGroot and Bernardo \cite{Bernardo1979, DeGroot1986, Stone1959a, Stone1959b} chose a utility function based on the Shannon information. This leads to the Bayesian 
$D$-optimality, in which the expected gain in the Shannon information is maximized. In this case, the determinant of the
inverse of the information matrix is minimized, which is equivalent to maximizing the determinant of the information matrix,
\begin{equation}
	\displaystyle \min \hspace{2mm} \phi_{D} = |(A^{\top}A)^{-1}|,
\end{equation}  
where $A^{\top}A$ is the information matrix.

The Bayesian $D$-optimality is the most common criterion for computer-generated optimal designs. It aims to maximize the geometric mean 
of the eigenvalues of $A^{\top}A$.  

\subsubsection*{Bayesian $A$-optimality}
The Bayesian $A$-optimality \cite{Brooks1974, DuncanDeGroot1976, Owen1970} is based on the sum of the variances of the estimated 
parameters for describing the model. Consequently, it minimizes 
the sum of the diagonal elements, the trace of the inverse of the information matrix,
\begin{equation}
	\displaystyle \min \hspace{2mm} \phi_{A} = {\rm trace}(A^{\top}A)^{-1}.
\end{equation}
This criterion results in minimizing the average variance of the estimates.

\subsubsection*{Bayesian $C$-optimality}
The Bayesian $C$-optimality \cite{Chaloner1984} is a special case of the $A$-optimality. It is used for estimating some linear combination
of the parameters of interest,
\begin{equation}
	\displaystyle \min \hspace{2mm} \phi_{C} = c^{\top}(A^{\top}A)^{-1}c,
\end{equation} 
where ${\rm rank}(A) = 1$ and $A = cc^{\top}$.

\subsubsection*{Bayesian $E$-optimality}
Another design is the Bayesian $E$-optimality \cite{ChalonerVerdinelli1995}, in which we maximize the smallest eigenvalue of 
the information matrix. This is another natural approach because the eigenvalue spectrum is also one way to measure the 
largeness of a matrix,
\begin{equation}
	\displaystyle \max \hspace{2mm} \phi_{E} = \lambda_{\rm min} (A^{\top}A).
\end{equation} 

\subsection{Numerical Challenges}

The above designs are difficult to solve for large-scale problems. While developing numerical methods for experimental design 
of small-scale problems is relatively straightforward, large-scale problems, with a large number of experiments and parameters,  
present a difficult challenge. Non-trivial matrix-functions and their derivatives need to be evaluated and large-scale 
constrained optimization problems need to be solved. In addition, since many of the matrices in experimental design of 
inverse problems are large and dense, one is restricted to perform only matrix-vector products. Furthermore, some of the 
matrices may be ill-conditioned which generates further difficulties. 

As far as we know, the computation of Bayesian experimental 
design, especially for large-scale problems, has not been studied extensively. Therefore, it is imperative to develop sophisticated 
algorithms to solve large scale problems efficiently. Nonetheless, experimental design of large-scale, ill-posed problems
is an emerging application whose treatment requires the development of new algorithms.

The goal of optimal experimental design (OED) is to control, in some sense, the quality of experimental
results. For well-posed problems this usually means controlling the variance of unbiased least-squares estimates
of the parameters of interest. Although the study of OED for well-posed problems is well established 
\cite{bauer, ChalonerVerdinelli1995, fedorov, Pukelsheim93}, its application to ill-posed problems has not yet received much attention. 
A difficulty that arises in ill-posed problems is that the estimates are biased and this
bias may dominate the overall error. Since the bias depends on the parameters to be recovered, its control hinges
on prior information about the plausible parameters. For example, properties of the bias can be learned from the training data. 
It then makes sense to choose a design that minimizes the average mean squared error (under the prior). That is, we look 
for a design that minimizes the Bayes risk. This is a generalization of the $A$-optimal experimental design for well-posed problems.

As is well known, the posterior mean minimizes the Bayes risk but, in some cases, the evaluation of its risk
may be computationally demanding. One may then choose to minimize the Bayes risk in a class of estimators
that are computationally tractable. We develop efficient numerical methods to minimize the Bayes risk over the
class of affine estimators. Note that this includes the case Gaussian-prior-Gaussian-likelihood that is often used.

Classical $A$-optimal designs have been treated in the literature (e.g.,
\cite{fedorov,Pukelsheim93}) but we do not know of algorithms capable of
dealing with large-scale problems. Furthermore, when applying the traditional designs to problems that arise
in geoscience and medical imaging, we have found them to be rather limited; new design criteria are clearly
needed. 

In this thesis, we start by presenting new design criteria that combine sparsity constraints with Bayes
risk minimization. Based on different assumptions of the prior information, two different designs were developed to 
handle ill-posed problems. If the covariance matrix of the models is available, we define the $A_{B}$ design. 
Otherwise, the $A_{\pi}$ design is suggested. We then develop the optimization and linear algebra techniques 
required to apply the designs to realistic large-scale problems.

\section{Optimal design for regularization}

In various areas such as medical imaging, geophysics and tomography, inverse problems frequently arise in which parameters 
of a model can be obtained from some noisy observed data. Many of these problems are ill-posed, that is, either there 
is no unique solution to the problem or a small perturbation of the data can lead to a large change in the recovered
solution. A variety of techniques for obtaining stable solutions have been developed; these techniques are so called 
regularization techniques, (see for example \cite{Bertero1998, Engl1996, Morozov1993, Vogel2002} and references within). 
Initially, regularization techniques aim to stabilize the problem and incorporate 
a priori information or assumptions about the desired solution.

\subsection{Review of regularization}

From a Bayesian point of view, many regularization techniques correspond to imposing certain prior information on the model 
to be recovered. Based on different prior knowledge, people have developed various regularization functions, among
which the most commonly used one is the Tikhonov regularization.

For linear ill-posed problems, the Tikhonov regularization takes the regularization functional as 
$R(x) = \frac{1}{2} \parallel Lx \parallel_{2}$. The regularization matrix $L$ usually plays a role as a penalty function, 
such as restrictions for smoothness or bounds on the vector space norm. If the goal is to obtain a solution with smaller norms, 
then $L$ can be chosen as an identity matrix. If prior information indicates that the solution needs to be smooth, then 
a discretization of a differential operator is commonly used \cite{Tikhonov1977} to force smoothness on the desired solution.
Other regularization techniques such as TSVD \cite{Hansen1998, Nagy}, for which the original ill-conditioned problem is replaced
by a well-conditioned rank-deficient problem, have also been considered.

In some cases, especially when the problem size is large and only matrix-vector product is available, iterative 
regularization methods based on Golub Kahan bidiagonalization \cite{Golub1996} and LSQR \cite{HaberThesis, Hansen1998} are very useful. 
In order to avoid drawbacks of the iterative methods, such as the convergence of small singular values in the Krylov
subspace, hybrid methods \cite{HaberThesis, OLeary} have become popular. Hybrid methods first use a subspace to reduce the
size of the problem while capturing all the large singular vectors. After that, a Tikhonov regularization is applied based
on the derived subspace.

Other regularization techniques that use $\ell^{1}$ norms are now commonly used (see for example \cite{Donoho2003, Elad2002}) 
but they require considerably more sophisticated algorithms. Although it is often claimed that the $\ell^{1}$ approach is much 
better than the quadratic $\ell^{2}$ approach (such as the Tikhonov regularization), in our experience, this may not be the case 
for many problems, especially if $L$ is chosen appropriately. The question that we believe has not been fully answered is, how well 
can algorithms that use $\ell^{2}$ regularization perform when the regularization operator is chosen judicially?

\subsection{Motivation of finding the optimal regularization}

Most of the literature focuses on improving the solution by choosing the optimal regularization parameter 
\cite{Hanke1993, Vogel2002, Wahba1977}. However, not much attention has been paid to the choice of the regularization functional. 
In many cases, the commonly used regularization functionals have limitations. For instance, the differential operator 
treats all pixels of the image equally. Hence, the spatial features at different locations of the image are 
ignored. In order to allow different treatments for different spatial features within the image, the goal of this thesis is 
to develop an optimal $\ell^{2}$ regularization technique in order to obtain inexpensive and good solutions to 
the inverse problem. 

The idea of finding the optimal regularization functional is not new. Lauter and Liero \cite{Lauter_ill-posedinverse} have analyzed the
optimal regularization for ill-posed problems by corrections of the data or the operator.
Sugiyama and Ogawa \cite{Sugiyama2001, Sugiyama2002a, Sugiyama2002b} developed a method of choosing the optimal regularization functional and regularization parameter
from given candidates based on the subspace information criterion. Their method gives an unbiased estimate of the
generalization error with finite samples under certain conditions.
Given a training set of feasible solutions, Haber and Tenorio \cite{HaberTenorio2003} proposed a supervised learning approach to estimate the regularization functional,
which belongs to a family of functions parameterized by a vector parameter.

However, none of the above approaches consider the sparsity structure, which is more practical for real applications,
when choosing the optimal regularization functional. In this work, we develop an
optimal regularization functional methodology for ill-posed problems that solves an optimization problem
based on the MSE of the estimate. In order to accomplish this goal, we assume the availability of a set of examples as 
our training models.

This methodology can be applied to many practical inverse problems by choosing different types of
training models in various applications. Although we only consider the linear case in this thesis,
it can be generalized to nonlinear problems. Moreover, the optimization problem we propose here
is formulated such that only standard unconstrained optimization techniques are needed.

\section{Optimal design in $CO_{2}$ injection monitoring}

With the rapid development of modern technology, pollutants have been collecting in the atmosphere. As a result of this accumulation, 
the radiation from the sun enters the atmosphere but cannot leave. This is called the
Greenhouse Effect. These pollutants are called Greenhouse Gases, which mainly contain carbon dioxide ($CO_{2}$), methane,
ozone, etc. 

The main complication of the Greenhouse Effect is the increase in average temperature of the Earth, which results in
climate change and the melting of the polar ice caps, among other terrifying consequences. Therefore, it is imperative to
obtain control over the amount of $CO_{2}$ near the Earth's surface.

\subsection{Background of $CO_{2}$ injection}

$CO_{2}$ is mainly produced by burning fossil fuels. Even if a hydrogen-based economy was technically feasible today, we would
still have a long way to go before the world would be able to live without carbon-based fossil fuels. Until that time comes,
carbon combustion will continue to produce Greenhouse Gases. Thus scientists need to find a way to reroute them.

Many techniques have been proposed for capturing or removing $CO_{2}$ from the atmosphere, among which the idea of piping liquefied
$CO_{2}$ deep under the ocean or underground has received the most interest lately. Depleted oil and gas reservoirs and saline aquifers 
are generally considered prime candidates for large scale storage of $CO_{2}$. Preliminary studies suggest that the 
underground storage capacity is sufficient for the storage of hundreds of years’ worth of $CO_{2}$ injection, and that potential 
storage sites exist worldwide.

Basically, $CO_{2}$ injection uses compressors to force compressed $CO_{2}$ down a long pipe drilled into the underground
reservoir. However, as time passes, the liquefied $CO_{2}$ will begin to evaporate and potentially escape to the surface,
which could result in acidifying the water, affecting the soil chemistry and suffocating animals or people.
Therefore, geologic studies need to be carried out to model and monitor the storage condition of the compressed 
underground $CO_{2}$.

\subsection{Motivation of design for $CO_{2}$ injection monitoring}

The injection of $CO_{2}$ causes a fluid substitution within the pore space. Studies have shown that time-lapse borehole 
and surface seismic surveys can be used to estimate the location of the injected $CO_{2}$ and monitor changes in reservoir 
properties \cite{Arts2002, Gritto2004, Hoversten2003}.

The goal of seismic surveys is to record sound waves that travel through the media in the underground reservoir.
It generates quantitative maps of variations in fluid saturation or pressure over spatial domains.
In particular, crosswell seismic methods have been successfully applied to $CO_{2}$ injection monitoring \cite{Harris1995, Marion1997}.
These methods detect changes in seismic velocity caused by $CO_{2}$ injection into reservoirs.
Based on the crosswell survey, the $CO_{2}$ saturation between the wells is spatially mapped using tomographic imaging. 

In order to obtain continuous crosswell seismic data, scientists need to install seismic sources and receivers using production 
tubing with a geochemical fluid sampling system \cite{Daley2007}. The idea is, if the $CO_{2}$ saturation and/or plume thickness 
increases along a given raypath, the travel-time would decrease. This would allow detection with some spatial resolution, 
especially in the vertical direction. The difference in the travel-times recorded at different depths in an observation well in
the close vicinity allows continuous monitoring of the growing $CO_{2}$ plume.

Using time-lapse borehole surveys, in stead of inverting for each velocity field, data were inverted based on the change in 
velocity. In other words, the data used in the tomographic inversion is the travel time difference between the post-injection
time and the pre-injection time for each source and receiver pair. By inverting the difference data, some potential errors, 
such as miscalculation of the source and receiver locations, are minimized or possibly even eliminated \cite{AjoFranklin2006, Spetzler2006}.
       
The major problem of the existing $CO_{2}$ injection monitoring is, geophysical experiments or surveys are expensive
to perform. For example, in order to enable repeatability and control the sampling distance of the monitoring surveys,
it is better to deploy permanent sensor arrays in boreholes and on the ocean floor. However, the cost of permanent
built-in sensors is usually extremely high. Furthermore, it is very difficult to reconfigure sensors after they are installed.

For this reason, it is imperative to design an optimal experiment based on the relation between the experimental cost
and the geophysical purpose. A good design should be an optimal trade-off between the expected information about the
parameters or models of interest and the cost of acquiring such information. In other words, an optimal experimental design
should provide maximum information about the target geophysical structure at minimum financial cost.

Since 1995, tremendous work has been done to develop optimal design theories in various applications. Famous examples
include determining optimal seismometer locations for locating earthquakes \cite{Steinberg1990}, searching optimal
placement for sensors in time-lapse traveltime tomography \cite{AjoFranklin2008, Curtis1999}, designing optimal
sensor geometries for acoustic tomography in detecting underwater conductivity structure \cite{Barth1990}, etc. 

Besides the geophysical purpose and parameterized model, the third component of an experimental design framework is
the experimental constraints, which are usually unique for each experiment. The biggest concern and difficulty of $CO_{2}$ injection
monitoring is spatial limitations of placing sensors because, due to the geophysical complexity, it is usually not possible 
to place sources and receivers everywhere.

Therefore, the main purpose of this work is to design an optimal placement of sources and receivers in order to obtain 
the best possible estimated model by minimizing the averaged MSE in the target zone through optimization approaches. 
The optimal reconstructed model is obtained based on the resulting optimal design thereafter.

\section {Overview of this thesis}
        
The goal of this thesis is to develop numerical methods to efficiently solve large-scale optimal experimental design problems
and to apply our design ideas to various applications such as regularization and $CO_{2}$ injection monitoring. The thesis
is mainly divided into three parts. 

The first part consists of Chapters 2-3, which cover the derivation and numerical methods of two new $A$-optimal designs.         
\begin{itemize}
	\item Chapter 2: \\ 
	In Chapter 2, we give a quick review of the classical $A$-optimal design for the well-posed case, which serves to review the
basic definitions and introduce the sparsity-constrained design. We generalize the design to ill-posed
cases where the estimates are biased. We discuss the algorithms to treat the sparsity constrained designs starting with 
the computation of the design functions and their derivatives. We then proceed to describe approximate solvers for the 
linear system subproblems that are the bottleneck of the computations. This chapter concludes with a discussion of the 
numerical optimization methods. We provide two examples: A one-dimensional example inspired by an actual magnetotelluric 
application and a two-dimensional, super-resolution example where the goal is to determine an optimal number of lower 
resolution MRI images required to recover a higher resolution one. \\
	\item Chapter 3: \\
	In this Chapter, we present numerical methods for solving the Bayesian $E$-optimal design. Derivative techniques based on inverse
iteration are studied. The performance of this design is examined on both a 1D problem and a borehole ray tomography experiment.
\end{itemize}

In the second part of this thesis (Chapters 4-5), an approach to find the optimal $\ell^{2}$ regularization is proposed
and some numerical techniques to solve the corresponding optimization problems are discussed.
\begin{itemize}
	\item Chapter 4: \\
	In Chapter 4, an optimality criteria for finding the optimal $\ell^{2}$ regularization matrix is developed. 
We show that this matrix can be obtained by solving an optimization problem that depends on our a priori information. Special
derivative techniques are explored to solve large-scale matrix-based optimization problems.\\
	\item Chapter 5: \\
	In Chapter 5 we impose sparsity constraints on the structure of the optimal regularization matrix, such as the local
diagonal pattern, the $\ell^{1}$ norm pattern and the Kronecker product pattern. We discuss the numerical solution of the 
problem and propose an effective algorithm for the recovery of that matrix. We experiment with different regularization 
matrices based on several sparsity patterns. Their performances are discussed and compared with other commonly used regularizations. 
\end{itemize}

The third part (Chapter 6) is dedicated to optimal design for $CO_{2}$ injection monitoring. A framework based on MSE is developed
and the crosswell constraints for possible locations of sources and receivers are described in detail. The constrained optimization
is solved by DIRECT which is commonly used for finding the global minimum. Our algorithms were tested on  
a synthetic geographic tomography example provided by Dr. Jonathan Ajo-Franklin.

Chapter 7 is a summary of this thesis. We conclude the main work in this chapter and give opinions about the future work to 
extend the approaches proposed. 