\chapter{Numerical methods for $A$-optimal design}
\label{Ch2}

In the first part of this thesis, we aim at finding the optimal weight matrix $W$. We start by reviewing the classical 
framework of experimental design for discrete linear problems.

The data vector $d \in {\rm I\!R}^N$ and unknown model $m \in {\rm I\!R}^k$ are related via indirect noisy measurements
\begin{eqnarray}
\label{eq:for}
d = Am + \epsilon.
\end{eqnarray}
The rows of the matrix $A$ represent experiments; repeated rows correspond to repeated observations of the same experiment. 
Hence, if there are $n$ different experiments and each is observed $k_i$ times, then $\displaystyle \sum_{i=1}^n k_i = N$ is the total number 
of experiments. The random noise vector $\epsilon$ is assumed to have zero mean and covariance matrix $\sigma^2 I$.

Assuming that each experiment has its own variance, by defining $w$ to be the inverse variance, $w = \frac{1}{\sigma^2}$, 
we are able to plant the information of variances into the inverse problem as follows:
$$
\sqrt{W}Am + \sqrt{W}\epsilon = \sqrt{W}b.
$$
This implies that $\sqrt{W}\epsilon \sim N(0, 1)$. Now the problem of finding the best $W$ really becomes the problem of 
finding the optimal variances for conducting the experiments. We see that, if a certain entry on the diagonal of $W$ is 
very small, which corresponds to large variance, we say that particular experiment might be very unreliable. On the other 
hand, if the entry is very large, which corresponds to small variance, then we can trust the corresponding experiment more. 

Therefore, the goal of this work is to find the optimal $W$, thereby telling us which experiments are important for giving 
optimal estimates that provide reliable conclusions. The optimality is defined by choosing a function that, in some way, 
measures the performance of the estimate; an optimal experiment is one that minimizes such optimality function subject 
to appropriate constraints.

\section{The well-posed case}
\label{sec:wellposed}

In classical well-posed linear experimental design \cite{bauer,ChalonerVerdinelli1995,fedorov,Pukelsheim93} the matrix 
$A^{\top}A$ is nonsingular and $m$ is typically estimated using least-squares:
$$
\widehat{m} = \arg\min \frac{1}{2} \parallel Am-d \parallel^2. 
$$
This is equivalent to minimizing 
$$
S_k(m)= \frac{1}{2} (\bar{d} - A_{n} m)^\top W (\bar{d} - A_{n} m),
$$
where $\bar{d}_i$ is the average of the $k_i$ observations of the $i$th experiment,
the rows of $A_n$ are the $n$ rows of $A$ corresponding to different 
experiments and $W = \diag\{k_1/N,...,k_n/N\}$. 
Note that $(\sigma^2/N)W^{-1}$ is the covariance matrix of $\bar d$. The selection of $k_i$
controls the variances of the observations of different experiments. 

Since $\widehat{m}$ is an unbiased estimator, the optimality criteria used to choose $k$ are 
usually based only on characteristics of the covariance matrix of $\mh$ given by
$$
\Sigma_{\mh}(W) = \frac{\sigma^2}{N}(A_{n}^\top W A_{n})^{-1}.
$$
The proof of the above covariance matrix is not difficult.
\begin{eqnarray*}
\Sigma_{\mh}(W) & = & \Ex (\mh \mh^{\top}) \\
				& = & \Ex \left[ (A_{n}^\top W A_{n})^{-1}A_{n}^{\top}W\bar{d}\bar{d}^{\top}WA_{n}(A_{n}^\top W A_{n})^{-1} \right] \\
				& = & (A_{n}^\top W A_{n})^{-1}A_{n}^{\top}W\Ex (\bar{d}\bar{d}^{\top})WA_{n}(A_{n}^\top W A_{n})^{-1} \\
				& = & (A_{n}^\top W A_{n})^{-1}A_{n}^{\top}W{\rm Cov}(\bar{d})WA_{n}(A_{n}^\top W A_{n})^{-1} \\
				& = & \frac{\sigma^2}{N}(A_{n}^\top W A_{n})^{-1}.
\end{eqnarray*}

To this end, one defines an optimality function $\phi$ that scalarizes
covariance matrices; different choices of $\phi$ may define different designs. 
For example, an $A$-optimal design is defined by the function 
$$
\phi_{A}\left(\Sigma_{\mh}\right) = {\rm E} \parallel m-\mh\parallel^2 = {\rm MSE}(\mh) = \sigma^{2} \mbox{trace}(\Sigma_{\mh}),
$$
which is proved as follows
\begin{eqnarray*}
\phi_{A}\left(\Sigma_{\mh}\right) & = & \Ex (\parallel (A_{n}^\top W A_{n})^{-1}A_{n}^{\top}W\bar{d} - m \parallel_{2}^{2}) \\
								  & = & \Ex (\parallel (A_{n}^\top W A_{n})^{-1}A_{n}^{\top}W(A_{n}m+\bar{\epsilon}) - m \parallel_{2}^{2}) \\
								  & = & \Ex (\parallel (A_{n}^\top W A_{n})^{-1}A_{n}^{\top}W\bar{\epsilon} \parallel_{2}^{2}) \\
								  & = & \Ex (\bar{\epsilon}^{\top}WA_{n}(A_{n}^{\top} W A_{n})^{-2}A_{n}^{\top}W\bar{\epsilon}).
\end{eqnarray*}
By denoting $H = WA_{n}(A_{n}^{\top} W A_{n})^{-2}A_{n}^{\top}W$, we have
\begin{eqnarray*}
\phi_{A}\left(\Sigma_{\mh}\right) & = & \Ex (\bar{\epsilon}^{\top} H \bar{\epsilon} ) \\
                                  & = & \Ex (\sum^{n}_{j=1} \sum^{n}_{i=1} \bar{\epsilon}_{i}\bar{\epsilon}_{j}H_{ij}) \\
                                  & = & \Ex (\sum^{n}_{i=1} \bar{\epsilon}_{i}^{2} H_{ii} + \sum^{n}_{i,j=1;i \neq j} \bar{\epsilon}_{i}\bar{\epsilon}_{j}H_{ij}).
\end{eqnarray*}
In general, the expected value operator is not multiplicative, i.e. $\Ex (xy) \neq \Ex (x) \Ex (y)$.
The lack of multiplicativity gives rise to the study of covariance. \\
For $i = j$,
$$
\Ex (\sum^{n}_{i=1} \bar{\epsilon}_{i}^{2} H_{ii}) = \sum^{n}_{i=1} \Ex (\bar{\epsilon}_{i}^{2}) H_{ii} =  \sum^{n}_{i=1} {\rm Cov}(\bar{\epsilon}_{i}, \bar{\epsilon}_{i})H_{ii} = \frac{\sigma^{2}}{N} {\rm trace}(W^{-1}H).
$$
For $i \neq j$,
$$
\Ex (\sum^{n}_{i,j=1;i \neq j} \bar{\epsilon}_{i}\bar{\epsilon}_{j}H_{ij}) = \sum^{n}_{i,j=1;i \neq j} \Ex (\bar{\epsilon}_{i}\bar{\epsilon}_{j})H_{ij} = \sum^{n}_{i,j=1;i \neq j}{\rm Cov}(\bar{\epsilon}_{i}, \bar{\epsilon}_{j})H_{ij} = 0,
$$
since ${\rm Cov}(x,y) = 0$ if $x$ and $y$ are independent.
Therefore,
\begin{eqnarray*}
\phi_{A}\left(\Sigma_{\mh}\right) & = & {\rm trace}\left[ \frac{\epsilon^{2}}{N} W^{-1}WA_{n}(A_{n}^{\top} W A_{n})^{-2}A_{n}^{\top}W \right] \\
                                  & = & \frac{\epsilon^{2}}{N} {\rm trace}\left[ W^{-1}WA_{n}(A_{n}^{\top} W A_{n})^{-2}A_{n}^{\top}W \right] \\
                                  & = & \sigma^{2} \mbox{trace}(\Sigma_{\mh}).
\end{eqnarray*}
There are other popular optimalities  as well. For example, the $D$-optimal design is defined by $\phi_{D} \left(\Sigma_{\mh} \right) = \mbox{det}(\Sigma_{\mh})$,
and if the goal is to estimate the linear functional $c^\top m$ using the estimator
$c^\top\mh$, then a $C$-design is defined by 
$$
\phi_C\left(\Sigma_{\mh}\right) = {\rm MSE}(c^\top\mh)=c^\top\Sigma_{\mh} c.
$$
However, we will only consider the $A$-design because it has a natural generalization to ill-posed problems.

Given a chosen optimality function $\phi$, $k$ is selected by solving
the integer optimization problem:
$$
\widehat{k} = \arg\min \phi\left[\,\Sigma_{\mh}(W)\,\right]
\quad \mbox{s.t.}\quad \sum_i k_i = N, \,\, W = \diag\{\,k_i/N\,\}.
$$
The issue here is that $k$ is a vector of integers. So solving this problem requires integer optimization techniques. 
That is difficult so an approximation that is much easier to solve is defined by a
`relaxation' of the problem; one solves for real fractions that sum to
one instead of integers summing to $N$ (in the literature this design
goes by different names such as  `relaxed design' \cite{boydBook} and `continuous design' \cite{fedorov,Pukelsheim93}).
For example, the $A$-design is changed to:
\begin{equation}\label{eq:wwellposed}
\widehat{w} = \arg\min {\rm trace} \left[\, \frac{\sigma^2}{N}\,(\,A_n^\top W A_n\,)^{-1}\,\right]
\quad \mbox{s.t.}\quad \sum_i w_i = 1,\,\, w\geq 0,
\end{equation}
with $W = \diag\{\,w\,\}$.
That is, the non-negative weights $w_i$ replace $k_i/N$. The 
estimated weights are used as follows: the $i$th experiment is to be conducted 
$\widehat{k}_i=\lbrack N \widehat{w}_i\rbrack$ times and $m$ is estimated by minimizing
\begin{equation}\label{eq:Whls}
\widehat{S}_k(m)=\,  \frac{1}{2}\, (\,\widehat{d} - A_n m\,)^\top \widehat{W} (\,\widehat{d} - A_n m\,),
\end{equation}
where $\widehat{d}_i$ is the average of the $\widehat{k}_i$ observations of the $i$th experiment
and $\widehat{W} = \diag\{\,\widehat{k}_1/N,...,\widehat{k}_n/N\,\}$.
Note that the solution of the optimization problem (\ref{eq:wwellposed}) 
provides a proportional allocation of each experiment; the fractions
$\widehat{w}_i$ do not depend on the noise level $\sigma$ or the total number of experiments $N$.
For a given $\sigma$, the total $N$ can be chosen to match a target MSE.

In many experimental settings one may be able to control the variance of the experiments
by means other than replication; for example, by changing the exposure time of a particular experiment.
In this case $N$ is a continuous variable that can be used as a tuning parameter to obtain optimal
target variances $\sigma^2/\widehat{k}_i$. 

While it is clear that the MSE decreases to zero
as $N\to\infty$, in practice there is usually a smallest
achievable variance $\sigma^2_{\min}$ for the experiments. In this case
the MSE cannot be arbitrarily small without violating this lower bound.
We now define an experimental design that focuses on finding optimal
variances subject to a lower bound.

Let $\sigma_i$ be the target variance of the $i$th experiment.
Set $w_i=1/\sigma_i^2$ and define $d_w$ to be the $n\times 1$ data vector
where the $i$th experiment is implemented to have variance $1/w_i$ for each $i$. 
An estimate of $m$ is obtained by minimizing
\begin{equation}\label{eq:Vw}
V_w(m)=\,  \frac{1}{2}\, (\,{d}_w - A_n m\,)^\top W (\,{d}_w - A_n m\,), 
\end{equation}
where again $W=\diag\{w_i\}$. This time the covariance matrix of $\mh$ is
\[
\Sigma_\mh(W) = (\,A_n^\top W A_n\,)^{-1}.
\]
The corresponding optimization problem for the $A$-design is:
\begin{equation}\label{eq:wwellposedwmax}
\widehat{w} = \arg\min {\rm trace} \left[\,(\,A_n^\top W A_n\,)^{-1}\,\right]
\quad \mbox{s.t.}\quad 0\leq w \leq w_{\max}.
\end{equation}

\section{The sparsity control design}
\label{sec:SCD}

Clearly the solution of (\ref{eq:wwellposedwmax}) is $w=w_{\max}$. However, in order to save the experimental cost, 
one needs to control the total number of different experiments. Such control is particularly important when the 
cost of conducting the experiments a second time is negligible compared to the cost of doing them the first time. 
For example, drilling a bore-hole to measure data once has a substantial cost but measuring again
in the same bore-hole has a negligible cost. 
There are also situations where only a few experiments are actually needed or 
reasonable to realize, so performing all experiments is obviously wasteful.
In this case it is better to select a few different experimental configurations
to be conducted a number of times. 

This implies that we would like the vector $w$ to be sparse. 
This is the idea studied in Haber's previous work \cite{HaberHoreshTenorio08} where they considered controlling the sparsity
of $w$ by minimizing 
\begin{eqnarray*}
\min_{w} \ \  \phi\left[\,\Sigma_\mh(W)\,\right] +\beta \|w\|_0 \quad
{\rm s.t.} \   \  0 \le w \le w_{\rm max},
\end{eqnarray*}
with $ \|w\|_{0} = \#\{\,w_i\neq 0\,\}.$ They have used this design in the context of 
ill-posed problems and called it sparsity control design.
However, since the problem with $\|w\|_0$ is of combinatorial complexity, it is often 
approximated using the $\ell^{1}$-norm $\|w\|_1$ instead of  $\|w\|_{0}$ 
\cite{HaberHoreshTenorio08}.
A sparsity controlled modification of (\ref{eq:wwellposedwmax}) is
\begin{equation}\label{eq:wspwellposed}
\widehat{w} = \arg\min\ \ {\rm trace} \left[\, (\,A_n^\top W A_n\,)^{-1}\,\right]
+\beta \sum w_{i}
\quad \mbox{s.t.}\quad 0\leq w\leq w_{\max}.
\end{equation}
The sparsity of the design is controlled with $\beta$.
To determine a reasonable value of this parameter one can study the trade-off in MSE 
reduction as a function of $\|w\|_0$, which is, in turn, controlled by $\beta$. 
The plot of ${\rm trace} [\,(A_n^\top \widehat{W} A_n)^{-1}\,]$ 
vs $\|w\|_0$, is often referred to as a Pareto curve.
If this curve has an L-shape one can argue that the $\beta$ corresponding to
the corner is a reasonable choice. Note, however, that an optimal $w_i$ may give
a variance that is larger than that of the actual instrument. If this happens
one needs to make a decision of whether to set $w_i$ to match the instrument
variance or to leave out the experiment.

\medskip
An even more realistic approach to determine $w$ would be to also consider the actual financial
cost of the experiments. For example, bore-hole experiments
are typically much more expensive than surface experiments, and deeper bore-holes 
are more expensive than shallow ones. However, such
cost-efficient methods are application dependent and will not be discussed here.

\section{Formulation for ill-posed problems}
\label{sec:illposed}

We now consider the case where the recovery of $m$ given the indirect noisy data 
is an ill-posed inverse problem. 

\subsection{The $A_{B}$ design}

To design an experiment one needs to have some information about the class of models
to be recovered and about the type of noise to be expected. It is thus reasonable
to consider a Bayesian framework where such information is given in terms of
distribution functions. We assume that the distribution of $d_w$ given
$m$ is Gaussian $N(Am,W^{-1})$ and $m$ is a random vector with prior distribution $\pi$.
Since we use MSE as the risk function, the idea is then to use the posterior mean $\mh$, which minimizes the Bayes risk 
among all estimators, as the Bayes estimate and choose $w$ that minimizes
its Bayes risk $\Ex_\pi \MSE(\mh)$. However, given the large-scale problems
we want to address, we reduce the computational cost by choosing $\mh$
to be the linear (or affine) function of the data $d_w$ that minimizes
the Bayes risk. To compute this estimator and its Bayes risk, only the first two
moments of the prior $\pi$ are required:
\begin{equation}\label{eq:bestlinear}
\mh(w) =  (A^\top W A + \Sigma_m^{-1})^{-1}(A^\top Wd_w + \Sigma^{-1}\mu),
\end{equation}
where $\mu$ and $\Sigma_m$ are, respectively, the prior mean and covariance matrix.
The Bayes risk of $\mh$ is
\begin{eqnarray}\label{eq:phiab}
\phi_{A_B}(W)  &=& {\rm trace} \left[\,(\,A^\top WA + \Sigma_m^{-1}\,)^{-1}\,\right],
\end{eqnarray}
in which we assume that the prior information is given by the covariance matrix. 

A modification of the optimal design (\ref{eq:wspwellposed}) for the ill-posed
case is obtained by  minimizing the Bayes risk with an $\ell^{1}$ penalty:
\begin{eqnarray}\label{eq:L1illposed}
\widehat{w} &=& \arg\min \phi_{A_B}(W)  + \beta\|w\|_1\\
&\mbox{s.t.}& 0\leq w_i\leq w_{\rm max},\quad  W = \diag\{\,w_i\,\}.
\end{eqnarray}
To select appropriate values of $\beta$, we plot the Pareto curve and make a plot of 
$\phi_{A_B}(\widehat{W})$ as a function of $\|w\|_0$.

\subsection{The $A_{\pi}$ design}

In some cases, it is difficult to compute the $A_{B}$ design because it requires the inverse of the covariance matrix $\Sigma_m$, 
which will be nontrivial when the problem size is large. Thus some approximations are necessary. 
We may use the same ideas to choose the optimal linear estimator of a 
particular type, for example, the Tikhonov estimate $\mh$ in (\ref{eq:Tikhest}):
\begin{eqnarray}
\mh(w) &=& \arg\min\,\, ({d}_w-A_n m)^\top W ({d}_w-A_n m) + 
\alpha \|Lm\|^2 \nonumber\\
&=& (A^\top W A +  \alpha L^\top L)^{-1}A^\top Wd_w\label{eq:Tikhest},
\end{eqnarray}
where $\alpha$ is a regularization parameter, $L$ is a chosen matrix (e.g., a discrete 
derivative operator) and $W=\diag\{w_1,...,w_n\}$ with $\Var({d_w}_i)= 1/w_i$.

However, one of the problems we have this time is that the estimator $\mh$ is biased and thus its
covariance matrix does not provide all the information required to assess its
performance. Furthermore, its mean squared error (MSE),
depends on the unknown $m$: The MSE of $\mh$ for a fixed $\alpha$ is
\begin{equation}\label{eq:MSEtik}
\MSE(\mh) = \alpha^2\|\,C(W)^{-1} L^\top L m\,\|^2 + 
{\rm trace} \left[\,C(W)^{-2} A^\top W A\,\right],
\end{equation}
where $C(W) = A^\top W A + \alpha L^\top L$. Hence we have two problems: (i) without knowing 
$m$ one cannot control the MSE; (ii) the MSE requires an appropriate selection of
$\alpha$ which should be adapted to $m$ and the noise level.

The objective is then to choose $w$ to minimize its Bayes risk given by the expectation of
the MSE in (\ref{eq:MSEtik}) with respect to the prior:
\begin{equation}\label{eq:Apieq}
\phi_{A_\pi}(w) = \alpha^2\|\,C(W)^{-1} L^\top L \,\mu\,\|^2 + \\
{\rm trace} \left[\,\left(\,C(W)^{-2} \left[\,\alpha^2\,L^\top L \,\Sigma_m L^\top L  + 
A^\top W A\,\right]\,\right)\,\right]\hspace{-1mm}.
\end{equation}
Here the $\alpha$ and $L$ are chosen by the experimenter based on previous experience.
Note that the function $\phi_{A_B}$ is used in the usual Bayesian $A$-design and
$\phi_{A_\pi}$ is just the Bayes risk of a Tikhonov estimator with
prior moment conditions \cite{ChalonerVerdinelli1995,fedorov,Pukelsheim93}.

We provide the proofs of \ref{eq:MSEtik} and \ref{eq:Apieq} in the following.
\begin{eqnarray*}
\MSE(\mh) & = & \Ex (\parallel \mh - m \parallel^{2}_{2})\\
		  & = & \Ex (\parallel (A^{\top}WA + \alpha L^{\top}L)^{-1} A^{\top} Wd_{w} - m \parallel^{2}_{2}) \\
		  & = & \Ex (\parallel (A^{\top}WA + \alpha L^{\top}L)^{-1} A^{\top} W(Am + \epsilon_{w}) - m \parallel^{2}_{2}) \\
		  & = & \Ex (\parallel (A^{\top}WA + \alpha L^{\top}L)^{-1} A^{\top} WAm + \\
		  & ~ & (A^{\top}WA + \alpha L^{\top}L)^{-1} A^{\top} W\epsilon_{w} - m \parallel^{2}_{2}) \\
		  & = & \Ex (\parallel (A^{\top}WA + \alpha L^{\top}L)^{-1}( A^{\top} WA + \alpha L^{\top}L - \alpha L^{\top}L)m + \\
		  & ~ & (A^{\top}WA + \alpha L^{\top}L)^{-1} A^{\top} W\epsilon_{w} - m \parallel^{2}_{2}) \\
		  & = & \Ex (\parallel -\alpha (A^{\top}WA + \alpha L^{\top}L)^{-1} L^{\top}Lm + \\
		  & ~ & (A^{\top}WA + \alpha L^{\top}L)^{-1} A^{\top} W\epsilon_{w} \parallel^{2}_{2}) \\
		  & = & \alpha^{2} \parallel (A^{\top}WA + \alpha L^{\top}L)^{-1} L^{\top}Lm \parallel_{2}^{2} + \\
		  & ~ & \Ex (\parallel (A^{\top}WA + \alpha L^{\top}L)^{-1}A^{\top}W\epsilon_{w}\parallel_{2}^{2}) \\
		  & = & \alpha^2\|\,C(W)^{-1} L^\top L m\,\|^2 + {\rm trace} \left[\,C(W)^{-2} A^\top W A\,\right].
\end{eqnarray*}
\begin{eqnarray*}
\phi_{A_\pi}(w) & = & \Ex_{\pi} \MSE(\mh) \\
				& = & \Ex_{\pi} \left[ \alpha^2\|\,C(W)^{-1} L^\top L m\,\|^2 \right] + {\rm trace} \left[\,C(W)^{-2} A^\top W A\,\right] \\
				& = & \alpha^2 \Ex_{\pi}(m^{\top}L^{\top}LC(W)^{-2}L^{\top}Lm) + {\rm trace} \left[\,C(W)^{-2} A^\top W A\,\right] \\
				& = & \alpha^2 \Ex_{\pi} \left[ {\rm trace} (m^{\top}L^{\top}LC(W)^{-2}L^{\top}Lm) \right] + {\rm trace} \left[\,C(W)^{-2} A^\top W A\,\right] \\ 
				& = & \alpha^2 \Ex_{\pi} \left[ {\rm trace} (mm^{\top}L^{\top}LC(W)^{-2}L^{\top}L) \right] + {\rm trace} \left[\,C(W)^{-2} A^\top W A\,\right] \\
				& = & \alpha^2 {\rm trace} \left[ \Ex_{\pi}(mm^{\top})L^{\top}LC(W)^{-2}L^{\top}L \right] + {\rm trace} \left[\,C(W)^{-2} A^\top W A\,\right] \\
				& = & \alpha^2 {\rm trace} \left[ (\Sigma_{m} + \mu \mu^{\top})L^{\top}LC(W)^{-2}L^{\top}L \right] + {\rm trace} \left[\,C(W)^{-2} A^\top W A\,\right] \\
				& = & \alpha^2 {\rm trace} (C(W)^{-2}L^{\top}L \Sigma_{m} L^{\top}L) + \\
				& ~ & \alpha^{2}\parallel C(W)^{-1}L^{\top}L\mu \parallel_{2}^{2} + {\rm trace} \left[\,C(W)^{-2} A^\top W A\,\right] \\
				& = & \alpha^2\|\,C(W)^{-1} L^\top L \,\mu\,\|^2 + \\
				& ~ & {\rm trace} \left[\,\left(\,C(W)^{-2} \left[\,\alpha^2\,L^\top L \,\Sigma_m L^\top L + A^\top W A\,\right]\,\right)\,\right].
\end{eqnarray*}
 
\subsection{More about the above designs}

Notice that if $\pi$ is Gaussian with $\mu=0$ and $\Sigma_m = (L^\top L)^{-1}/\alpha$,
then (\ref{eq:bestlinear}) is equal to (\ref{eq:Tikhest}) and (\ref{eq:bestlinear})
is the posterior mean, which minimizes the Bayes risk among all functions of
the data.

To understand the different nature of the designs, 
it is important to say a few words about the different interpretation of the results depending
on the type of chosen prior distribution. In the well-posed case the weights 
$w$ control the
trace of the covariance matrix of the estimator that the experimenter will use.
Similarly, in the ill-posed case with a Gaussian prior, the weights control the trace
of the covariance matrix of the posterior distribution that the experimenter
will obtain. When the prior is not Gaussian and
the weights are based on the affine estimator with smallest risk, then this risk
only provides a bound for the expected value of the trace of the posterior covariance matrix.
More precisely, let $\mh$ be any function of
the data (estimator) and $\Sigma_{m|y}$ the posterior covariance matrix
of $m$ given $y$. It is easy to see that
\[
\Ex \|\mh - m\|^2 = \Ex \|\mh - \Ex(m|y)\|^2 + \Ex\,{\rm trace} \,\Sigma_{m|y}. 
\]
Consider now the case of a Gaussian prior. Let $\mh_a$ and $\mh_L$ be, respectively,
the affine estimator that minimizes the risk and the Tikhonov estimator (\ref{eq:Tikhest}). Then
$\mh_a$ is the posterior mean $\Ex(m|y)$, the trace of $\Sigma_{m|y}$
does not depend on $y$ (so $\Ex\,{\rm trace} \,\Sigma_{m|y} = {\rm trace} \,\Sigma_{m|y}$)
and
\[
{\rm trace} \,\Sigma_{m|y} = \Ex \|\mh_L - m\|^2 - \Ex \|\mh_L - \mh_a\|^2
\leq \Ex \|\mh_L - m\|^2.
\]
We see that the risk of $\mh_L$ provides a bound for the trace of the posterior
matrix but it could be a poor bound if $\mh_L$ and $\mh_a$ are not close enough.
In the general case of non-Gaussian priors, the trace of $\Sigma_{m|y}$
may depend on $y$ and all we may be able to say is that
\[
\Ex\,{\rm trace} \,\Sigma_{m|y}\leq  \Ex \|\mh_a - m\|^2\leq \Ex \|\mh_L - m\|^2.
\]
Thus, the risks of $\mh_a$ and $\mh_L$ only constrain an average value
of ${\rm trace} \,\Sigma_{m|y}$ but, if the sampling variability
of the posterior covariance matrix is large, then neither
$\mh_a$ nor $\mh_L$ may provide a good control of the posterior distribution of
the actual experiment.

\section{Numerical optimization  of the $A_B$ and $A_{\pi}$ designs}
\label{sec3}

The $A_B$ and $A_{\pi}$ designs are difficult to compute for large-scale problems as
they involve the inverse of large matrices. In particular, it is
difficult to evaluate the objective function $\phi$ and its derivatives.  
We now discuss methods to do this efficiently.

\subsection{Evaluating the traces in the objective function}

The definition of the objective function $\phi$ includes the traces of large dense matrices.
Efficient approximations of such traces can be obtained using stochastic trace estimators
\cite{BaiFaheyGolub1996,gvmat2,hutchnison}. If the vectors $U_1$,...,$U_r$ are independent 
and each with independent entries taking the values
$1$ and $-1$ with equal probability, then the trace estimator 
\[
\widehat{T}_r(H) = \frac{1}{r}\sum_{i=1}^r U_i^{\top}HU_i
\]
 is an unbiased estimator of ${\rm trace} (\,H\,)$ with variance
\begin{eqnarray}\label{eq:trbnd}
\Var(\,\widehat{T}_r(H) \,) = \frac{2[\,{\rm trace} (H^2) - \|\,h\,\|^2\,]}{r},
\end{eqnarray}
where $h=\diag (H) $.
To estimate the quality of this estimator we would like to evaluate the relative error, 
that is we would like to estimate the variance \eqref{eq:trbnd}.
Let $\lambda_{1},\ldots,\lambda_{n}$ be the eigenvalues of $H$.
Also let
$$\bar{\lambda}= {\frac 1n}\sum_i\lambda_i$$.

We have that
$$ {\rm trace}(H^2) = \sum \lambda_{i}^{2}\ \ \  \ \  {\rm trace}(H)^{2} = \left(\sum \lambda_{i}\right)^{2} = n^{2} \bar \lambda^{2}$$.
We can therefore write
\begin{eqnarray*}
 \frac{{\rm trace}(H^2) - \|\,h\,\|^2}{{\rm trace}(H)^2} = \frac{ \sum \lambda_{i}^{2} - \|\,h\,\|^2}{n^{2} \bar \lambda^{2}}.
\end{eqnarray*}
To have an upper bound, we need to bound $\|h\|^{2}$ from below. We have that
$ h = \diag(H)$ and  ${\rm trace}(H) = \sum \lambda_{i} = n \bar \lambda = \sum h_{i}$.
We now seek the smallest $\|h\|^{2}$ subject to the above equality constraint. This can be easily achieved by minimizing
$$\min\ \|h\|^{2}\ \ \ {\rm s.t}\ \ n \bar \lambda = \sum h_{i}. $$
The solution to this problem is obviously
$h_{i} = \bar \lambda$.
We can therefore write
\begin{eqnarray*}
\frac{{\rm trace}(H^2) - \|\,h\,\|^2}{{\rm trace}(H)^2} \le \frac{ \sum \lambda_{i}^{2} - n^{2} \bar \lambda^{2}}{n^{2} \bar \lambda^{2}}.
\end{eqnarray*}
This shows that the relative variance of the trace estimate will be small if the 
scatter of the eigenvalues is small compared to their mean. 
In our applications, $H$ stems from an ill-posed problem. If no regularization is used, then
the eigenvalues are bounded from above and tend to cluster at $0$ as $n \rightarrow \infty$.
For a regularized version
in standard form, the eigenvalues cluster around the regularization parameter $\alpha$.
This implies that the total spread of the eigenvalues is typically small, which explains why stochastic
trace estimators are so successful in ill-posed problems. It is also clear that, for well posed problems
where $\lambda_{\rm max}/\lambda_{\rm min}$ is close to $1$, such estimators can be very effective.

Using randomized trace estimators with $r=1$, the $\phi_{A_{B}}$ and $\phi_{A_{\pi}}$ objective functions are 
replaced by the following approximations
\begin{equation}
\label{aopt}
\phi_{A_B}(w) = v^{\top}(\,A^\top WA + \Sigma_m^{-1}\,)^{-1}\,v
\end{equation}
and
\begin{equation}
\label{fropt}
 \phi_{A_\pi}(w) = \alpha^2\|\,C(W)^{-1} L^\top L \,\mu\,\|^2 
  +v^\top\left(\,C(W)^{-2} \left[\,\alpha^2\,L^\top L \Sigma_m L^\top L  + A^\top W A\,\right]\,
\right)v,
\end{equation}
where $v$ is a random vector whose entries take the values $\pm 1$ with equal probability. 
Note that the computation of $C(W)^{-1}$ is not needed explicitly. Instead, we need to evaluate 
the action of $C(W)^{-1}$ on a vector, say $q$, and this can be done 
by solving the linear system $Cz=q$, which can be done by conjugate gradient. We discuss 
this step at length in Section \ref{sec4}.

\subsection{Evaluating the derivatives} 

We start with the derivatives for the $A_{B}$ design:
$$ \grad_{w}\phi_{{A}_{B}} = \grad_{w}\left(v^{\top} (A^{\top}WA  + \Sigma_m^{-1})^{-1} v \right) $$
and setting
$$ z = (A^{\top}WA  + \Sigma_m^{-1})^{-1} v  \quad \leftrightarrow \quad
 (A^{\top}WA  + \Sigma_m^{-1}) z - v = 0$$ 
yields
$$ A^{\top}\, {\rm diag}(Az) + (A^{\top}WA  + \Sigma_m^{-1}) \grad_{w}z = 0. $$
Therefore
$$ \grad_{w}z = -(\,A^\top WA + \Sigma_m^{-1}\,)^{-1}A^{\top} {\rm diag}(Az), $$
which implies 
\begin{eqnarray}
\label{eq:gradphiA}
\grad_{w}\phi_{A_{B}} = -{\rm diag}(Az)\, A(\,A^\top WA + \Sigma_m^{-1}\,)^{-1}v =  -{\rm diag}(Az) Az = -Az\odot Az.
\end{eqnarray}

The computation of the objective function and its derivatives for the  $A_{B}$ design is 
summarized in Algorithm~\ref{alg:AC}
%------------------------------------------------------------------------------
\begin{algorithm}
\caption{Objective function and gradient for $A_{B}$ design}%
\label{alg:AC}
\begin{algorithmic}
\STATE(1) Solve the system $ (A^{\top}WA  + \Sigma_m^{-1}) z = v $
\STATE(2) Set $\phi_{A_{B}} = v^{\top} z$ and $\grad_{w}\phi_{A_{B}} =  -Az\odot Az$
\end{algorithmic}
\end{algorithm}
%------------------------------------------------------------------------------

An interesting  property of the $A_{B}$ design is that the computation of its gradient is almost
`free' compared to the computation of the objective function. After the vector $z$ 
is obtained, which requires a substantial amount of work, the objective function 
is computed  by a simple inner product and its derivative requires only one more matrix 
vector product with the matrix $A$.

We now derive similar expressions for the $A_{\pi}$ design. 
We rewrite \eqref{fropt} in a symmetric form
\begin{eqnarray*}  
\phi_{A_\pi}(w) & = & \alpha^{2}\mu^{\top}L^{\top}LC(W)^{-2}L^{\top}L\mu \\
& + & v^\top C(W)^{-1} \left(\, A^\top W A + \alpha^2\,L^\top L \,\Sigma_m L^\top L \,\right)\, C(W)^{-1}v. 
\end{eqnarray*}
Defining as above $z =  C(W)^{-1} v$, we can write
$ \grad_{w} z = -C(w)^{-1}A^{\top} {\rm diag}(Az) $.
Similarly, we define $y = C(W)^{-1}L^{\top}L\mu$ and write $\grad_{w}y = -C(w)^{-1}A^{\top} {\rm diag}(Ay)$.
We thus obtain
\begin{align}
\grad_{w} &\phi_{A_\pi}(w) = -2\alpha^{2}{\rm diag}(Ay)AC(W)^{-1}y + \diag(Az)Az\nonumber \\
&-2\,{\rm diag}(Az)\, A C(W)^{-1} \left(\, A^\top W A + 
\alpha^2\,L^\top L\, \Sigma_m L^\top L \,\right)z.\label{gradAFR}
\end{align}

% The algorithm for the evaluation the function and its derivatives is summarized in algorithm
% \ref{alg:FR}.
%% ------------------------------------------------------------------------------
% \begin{algorithm}
% \caption{Objective function and gradient for A and C design}%
% \label{alg:FR}
% \begin{algorithmic}
% \STATE(1) Solve the systems  $(A^{\top}WA  + L^{\top}L) p_{j} = L^{\top}Lm_{j}  \quad {\rm and} \quad (A^{\top}WA  + L^{\top}L) q = Av$
% \STATE(2) Set $\phi_{FR} = \alpha^{2} \sum \|p_{j}\|^{2} + \sigma^{2} \|q\|^{2}$
% \STATE(3) Solve the systems $C y_{j} = p_{j}$ and $C z = q$
% \STATE(4) Set $\grad_{w} \phi_{FR} = - \alpha^{2 }\sum{\rm diag}(Ap_{j}) A y_{j} - \sigma^{2}{\rm diag}(Aq) A z$ 
% \end{algorithmic}
% \end{algorithm}
%% ------------------------------------------------------------------------------

\subsection{Solving the linear systems } 
\label{sec4}

Regardless of the chosen optimization criteria, a key component in the
solution of the problem is the solution of the system
\begin{eqnarray}
\label{eq:linsys} 
 (\,A^{\top}WA + \Sigma_{m}^{-1}\,)m = A^{\top}Wb,
\end{eqnarray}
% which is equivalent to the least square problem
% \begin{eqnarray}
% \label{eq:linsysls} 
% \begin{pmatrix} W^{\hf}A \\  L \end{pmatrix} x = \begin{pmatrix} W^{\hf }b \\ 0 \end{pmatrix}
% \end{eqnarray}
where, in the $A_\pi$ design, $\Sigma_{m}^{-1}=\alpha L^{\top}L$.
This system is large and dense and therefore, iterative methods are typically used
for its solution. This system needs to be solved every time $w$ is updated, that is,
whenever we evaluate the objective function or its
derivatives. Therefore, we need a method that enables us to quickly evaluate $m$ given a new $w$.
We now use the SVD of $A$ to better understand the properties of the solution,
and then, we propose an approximation of the SVD using Lanczos decomposition.

\subsubsection{SVD analysis}

To understand how to solve the large linear system quickly, we first turn
to a simple SVD analysis.
Assume first that $\Sigma_{m}^{-1}=I$ (this case is referred to as `standard form').
The more general case can be transformed to this form \cite{Hansen1998}.

Write the SVD of $A$  as $A = U \Lambda V^{\top}$,
where $U$ is an orthogonal $n \times n$ matrix, $V$ is a $k \times n$
matrix with orthonormal columns, and $\Lambda$ is diagonal $n \times n$ with singular values
$\lambda_{1} \ge \lambda_{2} \ge \ldots \ge \lambda_{n}$.
Using this decomposition the linear system is rewritten as
\begin{eqnarray}
\label{eq:linsyt} 
 (\Lambda U^{\top}WU\Lambda  + I)\xi = \Lambda U^{\top}Wb,\quad  \xi = V^{\top}m.
\end{eqnarray}

We now start by reviewing the  `usual' unweighted case, which is commonly used to analyze the
solution of the inverse problem (see, for example, \cite{Hansen1998}). When $W= \mu I$ we obtain 
$\Lambda U^{\top}WU\Lambda =  \mu \Lambda^{2}$ and the system decouples:
$$  (\mu \Lambda^{2}  + I) \xi = \mu \Lambda U^{\top}b, \quad \quad \xi_{j} = {\frac {\mu \lambda_{j} }{\mu \lambda_{j}^{2} + 1}}\, u_{j}^{\top}b. $$ 
The ratio $\mu \lambda_{j}/(\mu \lambda_{j}^{2}+1)$ decreases to $1/\lambda$ for large $\lambda$
and is close to zero for small $\lambda$. It stands to reason that an
approximate solution can be obtained by truncating the singular values:
$$ \xi_{j} = \left\{ \begin{matrix}{\frac {1}{\lambda_{j}}}\, u_{j}^{\top}b  & \lambda_{j} \gg 0  \\
0 & {\rm otherwise.}
\end{matrix} \right. $$ 
This is the truncated SVD (TSVD) solution. The problem is that it may
be difficult to determine where to truncate the singular values. 
One option is to use hybrid methods that were proposed first in \cite{OLeary} 
and have been studied substantially in recent years. In this case the
solution will be $0$ only when its corresponding singular value is much smaller than the largest singular value $\lambda_1$. 
That is, we keep the first $s$ singular values such that $\lambda_s$ is less than some tolerance times the largest singular 
value and truncate everything after that. Hence, 
$$ \xi_{j} = \left\{ \begin{matrix}{\frac {\mu \lambda_{j} }{\mu \lambda_{j}^{2} + 1}}\, u_{j}^{\top}b  & 
j = 1, \cdots, s \hspace{3mm} {\rm where} \hspace{3mm} \lambda_{s} \le {\rm tol} \lambda_{1}  \\ 0 & {\rm otherwise.}
\end{matrix} \right. $$
This simple idea reduces the dimensionality of the problem and makes 
the solution less sensitive to the selection of a truncation level.

We now show that a similar analysis can be used in our case. If $W$ is not a multiple of the identity,
then
\begin{eqnarray*}
\Lambda U^{\top} W U \Lambda = 
\begin{pmatrix}
\lambda_{1}^{2} (w \odot u_{1})^{\top}u_{1} & \lambda_{1} \lambda_{2} (w \odot u_{1})^{\top}u_{2} & \cdot   & \cdot  &   \lambda_{1} \lambda_{n} (w \odot u_{1})^{\top}u_{n} \\
\lambda_{1} \lambda_{2} (w \odot u_{1})^{\top}u_{2} &  \lambda_{2}^{2} (w \odot u_{2})^{\top}u_{2} & \cdot   & \cdot  &   \lambda_{2} \lambda_{n} (w \odot u_{2})^{\top}u_{n} \\
\cdot   &      &    &   &  \\
\cdot   &      &    &   &  \\
\lambda_{1} \lambda_{n} (w \odot u_{1})^{\top}u_{n} &  \lambda_{2} \lambda_{n} (w \odot u_{2})^{\top}u_{n} & \cdot   & \cdot  &  \lambda_{n}^{2} (w \odot u_{n})^{\top}u_{n} \\
\end{pmatrix},
\end{eqnarray*}
where $\odot$ represents the Hadamard product, for which the product is taken entry-wise.
Thus, the system does not decouple as before. Nevertheless, 
% $$ |(w \odot u_{i})^{\top}u_{j}| \le \|w\|\,   \|u_{i}\|\,  \|u_{j}\| \le \sqrt{n} w_{\rm max}. $$
typically  $\lambda_{\ell} \ll \lambda_{1}$ for $\ell$ greater than
some index $s$, and since in addition $w$ is bounded, it makes sense
to use a TSVD as before; we obtain a reduced $s \times s$ dense system
$\Lambda_{s} U_{s}^{\top} W U_{s} \Lambda_{s}$ where 
$ U_{s} = \begin{pmatrix} u_{1}, \ldots, u_{s} \end{pmatrix}$ and 
$\Lambda_{s} = \diag ( \lambda_{1},\ldots , \lambda_{s})$.
If $s$ is not too large (say, up to a few thousand) then it is possible to use direct methods to
quickly solve the system for different $w$'s. The algorithm is described in Algorithm~\ref{alg:solSVD}.
\begin{algorithm}
\caption{SVD approximate solution of the system $(A^{\top}WA + I) m = A^{\top}Wb$ }%
\label{alg:solSVD}
\begin{algorithmic}
\STATE(1) Compute the SVD: $A = U \Lambda V^{\top}$
\STATE(2) Choose an index $s$ to truncate the SVD, $\lambda_{s} \le {\rm tol}\, \lambda_{1}$
\STATE(3) Solve the systems $(\Lambda_{s} U_{s}^{\top} W U_{s} \Lambda_{s} + I) \xi = \Lambda_{s} U_{s}^{\top} W b $
\STATE(4) Set $m = V_{s} \xi$ 
\end{algorithmic}
\end{algorithm}

The important observation is that the SVD is computed only
once at the beginning of the optimization. We then use the decomposition to solve the optimization
problem at a negligible cost. Thus, if the computation of the SVD is
not prohibitively expensive, then it is possible to quickly solve the problem.

\subsubsection{Approximation using Lanczos decomposition}

The computation of the SVD is typically not practical for large-scale problems; it has a computational
complexity of ${\cal O}(n^{2}k)$ and large storage requirements. We therefore turn to approximating
the SVD using Lanczos methods. 

Lanczos bidiagonalization has been extensively studied in the context
of inverse problems \cite{HaberThesis,Hanke1993,Hansen1998,nagypalmer2003}. In its simplest formulation, $s$ iterations of the Lanczos 
process are computed yielding an approximate decomposition
\begin{equation}\label{eq:decBs}
A \approx U_{s} B_{s} V_{s}^{\top}, 
\end{equation}
where $U_{s} =  [u_{1} , \ldots, u_{s}]$ is  an $n \times s$ matrix and  $V_{s} =  [v_{1} ,\ldots, v_{s+1}]$
is a $k \times (s+1)$ matrix, and both have orthonormal columns.
 The matrix $B_{s}$ is an $s \times (s+1)$ bidiagonal matrix. The computation of the
Lanczos decomposition can be done with or without re-orthogonalization of the vectors.
Without re-orthogonalization the vectors tend to loose the orthogonality, especially for large $s$.
This motivated the study of a number of re-orthogonalization techniques \cite{Hansen1998}. 
An important property of the decomposition (\ref{eq:decBs}) is that the singular values of $B_{s}$ approximate
the singular values  of $A$.  
In numerical experiments, it has been widely observed that the
large singular values of $A$ are approximated first, yielding an approximation to the truncated SVD solution.
This observation has motivated research on iterative regularization as well as hybrid regularization
techniques \cite{Nagy}. 

As we have done with the SVD, using the truncated bidiagonalization we obtain
 \begin{eqnarray}
\label{eq:linsyslanc} 
 (V_{s}B_{s}^{\top} U_{s}^{\top}WU_{s}B_{s} V_{s}^{\top} + I)m = V_{s}B_{s}^{\top} U_{s}^{\top}Wb,
\end{eqnarray}
and multiplying both sides of \eqref{eq:linsyslanc} by $V_{s}^{\top}$  from the left and setting $\xi = V_{s}^{\top}m$ yields
\begin{eqnarray}
\label{eq:linsylanct} 
 (B_{s}^{\top} U_{s}^{\top}WU_{s}B_{s}  + I)\xi = B_{s}^{\top}U_{s}^{\top}Wb.
\end{eqnarray}
Again, the system \eqref{eq:linsylanct} is a small $s \times s$ dense system and its solution
can be obtained quickly using direct methods. The algorithm for solving the system using the Lanczos process is
described in Algorithm~\ref{alg:solLanc}.  
%------------------------------------------------------------------------------
\begin{algorithm}
\caption{Lanczos approximate solution of the system $(A^{\top}WA + I) m = A^{\top}Wb$}
\label{alg:solLanc}
\begin{algorithmic}
\STATE(1) Compute the Lanczos decomposition of $A$ with starting vector $b$
\STATE(2) For each step compute the SVD of $B_{j}$, stopping at step $s$ when $\lambda_{s}(B_{s}) \le  {\rm tol}\, \lambda_{1}(B_{s}) $
\STATE(3) Solve the system $(B_{s}^{\top} U_{s}^{\top} W U_{s} B_{s} + I) \xi = B_{s}^{\top} U_{s}^{\top} W b $
\STATE(4) Set $m = V_{s} \xi$ 
\end{algorithmic}
\end{algorithm}
%------------------------------------------------------------------------------
We have experimented using the Lanczos process with
and without re-orthogonalization and found that re-orthogonalization did not yield any benefit.
We believe this is because, unlike the case $W = I$, we do not require $U_{s}^{\top}U_{s} = I$ and,
anyway, the product $U_{s}^{\top} W U_{s} $ is computed when obtaining the solution.

\subsection{Numerical optimization}
\label{sec5}

Given the tools described in the last section, we can now use standard
box-constrained optimization to solve the problem. We start by solving the optimization problem
\begin{eqnarray}
\label{eq:scd1norm2}
\min && \phi_{\beta} = \phi(w) + \beta \|w\|_1 \\
\nonumber
{\rm s.t} && 0 \le w \le w_{\rm max}.
\end{eqnarray}
% As in our previous work, we have used the projected BFGS method based
% on the algorithm of More and Lin (see \cite{kelley2, LinMore1999}) 
% with minor modifications as described in \cite{givoba, vogelbook}.
To solve the optimization problem (\ref{eq:scd1norm2}) we use
projected steepest descent as this method can easily deal with the non-differentiability of the 
$\ell^1$ norm. For further discussion see \cite{friedlander,wrightL1}.

We also consider the $\ell^0$ penalty $\beta \|w\|_0$ as it is
expected to provide the sparsest solution. To approximate $\|w\|_0$ and the sparsest solution, we 
use an approximation described in \cite{boydBook}.

We divide all experiments into two sets: ${\cal I}_{0}$ and ${\cal I}_{A}$. The set ${\cal I}_{0}$ contains all the indices for 
the zero entries of the solution to \eqref{eq:scd1norm2}, namely $w_{\cal I} = 0$ and ${\cal I}_{A}$ contains the rest. 
Assuming that ${\cal I}_{0}$ is known a priori, 
which can be done by solving the $\ell^{1}$ approach, the $\ell^{0}$ solution could be obtained by solving this 
un-regularized optimization problem only on the set ${\cal I}_{A}$. We do not need any regularization term here because 
the zero set is already known. It has been shown that the approximated $\ell^{0}$ solution may improve upon the $\ell^{1}$ solution. 

The $\ell^{0}$ penalty
is approximated by the solution of the problem
\begin{eqnarray}
\label{eq:scd1norm0}
\min_{{\cal I}_{A}} && \phi_{\beta} = \phi(w)  \\
\nonumber
{\rm s.t} && w_{{\cal I}_{0}} = 0, \hspace{3mm} w_{{\cal I}_{A}} \geq 0 \\
               &&          0 \le w_{{\cal I}_{A}} \le w_{\rm max}.
\end{eqnarray}

As we have shown in \cite{HaberHoreshTenorio08}, an important aspect of the optimization is that to study the
design space we need to solve \eqref{eq:scd1norm2} for different
values of $\beta$ that lead to different sparsity structure and thus provide
information about the design cost as a function of the number
of experiments. Here we have used a simple continuation strategy 
to achieve this goal \cite{hao}.

\section{Numerical experiments}
\label{sec6}

We illustrate the performance of our algorithms using a small-scale 1D problem, where
exact quantities can be easily computed and a realistic, large-scale super-resolution 
inverse problem.

\subsection{An ill-posed 1D magnetotelluric example}\label{sec:1Dilltoy}

The data are modeled as
\begin{eqnarray}
\label{eq:1D}
d_{j} = \int_{0}^{L}\, {\rm exp}(-\alpha_{j} x) \cos(\gamma_{j}x)\, m(x) \, dx + \epsilon_{j} \hspace{5mm} j = 1, \cdots, n.
\end{eqnarray}
This kernel mimics the linearized 1D magnetotelluric experiment \cite{oldenburg, wo}. 
The parameters $\gamma_{j}$ correspond to the recording frequencies
and the $\alpha_{j}$ depend on the frequencies and the background
conductivity. The objective is to select optimal values of these
parameters to best evaluate a 1D conductivity structure by solving the following inverse problem
$$
Am + \epsilon = d,
$$
where $A$ is a matrix generated by the kernel, $m$ is the conductivity and $d$ is the data we have observed.

The function $m$ is discretized using 256 points. Figure \ref{1Dkernel} shows an example of the 
magnetotelluric kernel that is used in our simulations.  The values of $\alpha_{j}$ 
are equally spaced between $0$ and $3$, namely $\alpha_{j} = {3j}/{100}$
for $0 \leq j \leq 100$ and the values of $\gamma_{j}$ are equally spaced between 
$1$ and $10$: $\gamma_{j} = 1 + {9j}/{100}$. Each row in the kernel matrix represents 
one pair $(\alpha_{j},\gamma_{j})$, which gives a total of 
$100$ different choices as our bank of experiments. The objective is to choose a subset of these $100$ 
experiments to best estimate $m$.
\begin{figure}[ht!]
 	\centering
 	\includegraphics[scale=0.4]{1Dkernelgray.jpg}
 	\caption{\label{1Dkernel}Geometry of the magnetotelluric kernel. The model is 
discretized using $256$ points. Each row of the kernel represents one pair
$(\alpha_{j},\gamma_{j})$.}
\end{figure}

The noise variance $\sigma^2$ is chosen for a SNR $\approx 10$.
The true test model is a realization from a Gaussian $N(0, \Sigma)$ where the covariance matrix 
$\Sigma$ is defined as $\Sigma_{i, j} = \exp{({-(x_{i}-x_{j})^{2}}/{2\tau^{2}})}$
with $x_i,x_j\in(0,2\pi)$. For simplicity, in this example, $\mu$ is set to be $0$. The parameter $\tau^2$ is chosen so that the
correlation between $m(x_i)$ and $m(x_j)$ is less than 50\% for $|x_i-x_j|>1$.
We set $w_{\max}=10^2/\sigma^2$. The optimization is done using steepest descent.

Notice that it is possible for the optimization to yield
$w_i<1/\sigma^2$, indicating that the $i$th experiment is to be
conducted with a variance larger than that of the instrument's noise. To
avoid this, the $w_i$ are thresholded as follows: tolerance
values $w_{\rm t}<w_{\min}$ are chosen and $w_i$ is set to zero if $w_i\leq
w_{\rm t}$, basically, we will not conduct that experiment and to $w_{\min}$ if $w_{\rm t}<w_i\leq w_{\min}$.
For the example we have used $w_{{\rm t}} = 10^{-5}/\sigma^2$ and $w_{{\rm min}} = 1/\sigma^2$.
From now on it will be understood that $w$ has been thresholded in this fashion.

\subsubsection{The $A_{B}$ design}

We start by illustrating the controlled sparsity of the $A_B$ design.
Figure \ref{1DLcurveAB} shows the risk as a function of $\|w\|_0$.
As is shown here, the risk decreases rapidly before
the corner of the L-curve  at $\|w\|_0=25$ and more slowly thereafter. 
Hence, this corner of the L-curve suggests an optimal compromise between the image reconstruction 
quality and the experimental cost. This implies that, even if you do more experiments, you are not going to improve the result much. 

\begin{figure}[ht!]
 	\centering
 	\includegraphics[scale=0.35]{1DLcurveModABL1black.jpg}
	\includegraphics[scale=0.35]{1DwwModABL1Riskgray.jpg}
 	\caption{\label{1DLcurveAB} The left panel shows
the risk as a function of sparsity ${\rm nnz}(w)=\parallel
w\parallel_{0}$  of the optimal $w$ obtained with the $A_{B}$ design.
The columns and rows in the image on the right correspond, respectively, to different 
experiments and different values of the risk. The color refers to the
values of the $w_i$.}
\end{figure}

The right panel in Figure \ref{1DLcurveAB} depicts all the $w$ vectors obtained from the
experiments. Each column represents one of the 100 different
experiments and the rows correspond to the different values of the
risk for the chosen values of $\beta$ used
to draw the L-curve. The color shows the values of $w_i$ for the corresponding experiment and risk.
The figure shows which experiments should be conducted to achieve the
corresponding risk; one can use the plot to decide if the reduction in risk is worth
the increase in the number of experiments. This figure tells us not only the number of experiments but which ones 
they are as well. This information can be quite helpful to practitioners.

Figure \ref{1DABL1} shows examples of model estimates using the
optimal $w_i$ determined by the corner in the L-curve.
The figure shows the true model as well as three reconstructions from
three different noise realizations. We see good reconstruction results using only the 25 selected 
experiments but we also see large sampling variability for $x>3$.

\begin{figure}[ht!]
 	\centering
 	\includegraphics[scale=0.35]{1Dtrue.jpg}
  	\includegraphics[scale=0.35]{1DAB25L1S3black.jpg}
 	\caption{\label{1DABL1}The left panel shows the true test
          model. The right panel shows examples
          of model estimates obtained with the $A_{B}$ design using
          $25$ optimal experiments for three different noise realizations.}
\end{figure}

To show the advantage of our optimal design, Figure \ref{1DABnaive} displays the
model estimates obtained using all $100$ experiments and also the
estimates based on $25$ equally-spaced naively chosen
experiments. Note also that the sampling variability for $x>3$ seems even larger.
The figures show that the $A_B$ design provides the best
results. Compare, for example, the amplitudes of the last two peaks of the curve.

\begin{figure}[ht!]
 	\centering
 	\includegraphics[scale=0.35]{1DABallS3black.jpg}
 	\includegraphics[scale=0.35]{1DABnaiveS3black.jpg}
 	\caption{\label{1DABnaive}Model estimates obtained using all
          $100$ experiments (left) and the equally-spaced naive design (right).}
\end{figure}

Table \ref{tab:ABl1} shows the relative errors ($\|m_{\rm true}-\widehat{m}\|/\|m_{\rm true}\|$)
of the reconstructions. Note that the estimates obtained using all $100$ experiments are worse 
than those using only the $25$ optimal experiments. This is because
the full design does not use weights to compute the estimates. The
optimal design selects a subset of the experiments and assigns to them
optimal weights (variances).

\begin{table}[ht!]
\caption{\label{tab:ABl1} Relative errors of the model reconstructions obtained with the $A_{B}$ design}
\centering
\begin{tabular}{c|c|c|c}
\hline\hline
design & error (noise 1) & error (noise 2) & error (noise 3) \\
\hline
$A_{B}$   & $7.75\%$  & $11.91\%$ & $14.81\%$ \\
all       & $10.29\%$ & $20.12\%$ & $19.26\%$ \\
naive     & $11.84\%$ & $23.94\%$ & $22.88\%$ \\
\hline\hline
\end{tabular}
\end{table}

\subsubsection{The $A_{\pi}$ design}

For the $A_{\pi}$ design we choose the matrix $L$ to be the discrete 1D Laplacian
operator. The parameter $\alpha$ is chosen by trial and error
so that given all the data we obtain the best possible recovery error.

\begin{figure}[ht!]
 \centering
 \includegraphics[scale=0.35]{1DLcurveModApiL1black.jpg}
 \includegraphics[scale=0.35]{1DwwModApiL1Riskgray.jpg}
 \caption{\label{1DLcurveApi}The left panel shows the risk as a function of sparsity 
${\rm nnz}(w)=\parallel w \parallel_{0}$ of the optimal $w$ obtained with the $A_{\pi}$ design.
The right panel is similar to that in Figure \ref{1DLcurveAB} but
          for the $A_\pi$ design using the $\ell^{1}$-optimization.}
\end{figure}

In Figure \ref{1DLcurveApi}, we choose the corner at $\|w\|_0=24$. This defines the optimal
solution used in the example. Figure \ref{1DLcurveApi} shows the
equivalent of Figure \ref{1DLcurveAB} for the $A_\pi$ design.

\begin{figure}[ht!]
 	\centering
 	\includegraphics[scale=0.35]{1Dtrue.jpg}
  	\includegraphics[scale=0.35]{1DApi24L1S3black.jpg}
 	\caption{\label{1DApiL1}The true test model (left) and the
          estimated models obtained with $\ell^1$- $A_{\pi}$ design
          using the $24$ optimal experiments with three different
          noise realizations (right).}
\end{figure}

\begin{figure}[ht!]
 	\centering
 	\includegraphics[scale=0.35]{1DApi24allS3black.jpg}
 	\includegraphics[scale=0.35]{1DApi24naiveS3black.jpg}
 	\caption{\label{1DApinaiveL1}Model estimates obtained with the
          $A_\pi$ design using all $100$ experiments (left) and the equally-spaced naive design (right).}
\end{figure}

Figure \ref{1DApiL1} shows the model estimates for three different
noise realizations obtained using the $\ell^1$ approach. Figure \ref{1DApinaiveL1} shows the estimates
obtained using all $100$ experiments and using $24$ equally-spaced naively chosen experiments.
The relative errors of the estimates are shown in Table
\ref{tab:Apil1}. The optimal design still leads to smaller errors
than the full and naive designs but the results are not as good as
those obtained with the $A_B$ design. This is to be expected as the
the $A_{B}$ design uses the Bayes estimate while the regularization for $A_{\pi}$ design is restrained to a particular type. 
We omit the rest of the results for the $\ell^{0}$ approach as they are very similar to those obtained
with the $\ell^{1}$ penalty.

\begin{table}[ht!]
\caption{\label{tab:Apil1} Relative errors of the model estimates
  obtained with the $A_{\pi}$ ($\ell^{1}$) design}
\centering
\begin{tabular}{c|c|c|c}
\hline\hline
design & error (noise 1) & error (noise 2) & error (noise 3) \\
\hline
$A_{\pi}$ $\ell^{1}$ & $39.57\%$ & $33.64\%$ & $45.98\%$ \\
all                  & $47.00\%$ & $41.07\%$ & $53.60\%$ \\
naive                & $47.00\%$ & $39.60\%$ & $53.87\%$ \\
\hline\hline
\end{tabular}
\end{table}

\subsection{Super-resolution}

In this section we consider a large-scale super-resolution problem that has a wide range of applications 
\cite{CHN06,ElFe97}. Super-resolution methods are techniques to enhance 
the resolution of an imaging system, in particular, to construct
a high-resolution image by combining a set of lower resolution images.

The data are modeled as
\begin{eqnarray}
\label{eq:srdata}
d_j = K\, S(u_j) m+\epsilon_{j}, \hspace{3mm} j = 1, \cdots, k;
\end{eqnarray}
where $d_{j}$ represent the low resolution images 
that have been collected, $K$ is the sparse matrix approximating the averaging process,
$u_{j}$ are the relative displacements among the low resolution images and $S(u_{j})$ is a 
sparse matrix representing the bilinear interpolation operation that connects the point in the
displaced image to the four pixel values in the reference image 
that surround it (see \cite{CHN06} for details). 
In order to obtain high-quality, high-resolution images, it is desirable to have many low
resolution ones. Nonetheless, generating and/or using many low
resolution images can be very costly. 
Thus, we would like to find an optimal subset of low resolution
images to generate one of high resolution and satisfactory quality.

We apply our algorithms to one of the MRI examples provided by
J. Orchard \cite{JeffOrchard}.
In order to reconstruct this $64 \times 64$ image, $100$ slightly different $32\times32$
low resolution images are generated. For displacements, both shifting
and rotation are considered. Figure \ref{lowSR} shows one of the
$100$ low resolution images that we generated together with the original high resolution image. 
\begin{figure}[ht!]
 	\centering
 	\includegraphics[scale=0.35]{EDSRlow.jpg}
 	\includegraphics[scale=0.35]{EDSRtrue.jpg}
 	\caption{\label{lowSR}One of the $100$ low resolution images (left) and the true high resolution image (right).}
\end{figure}

We apply the discretization of the gradient operator as a
smoothing matrix to obtain the solution of the inverse problem.
Since this is a realistic problem, we do not have the true covariance matrix of the MRI images
but do have other MRI images that can be used to estimate it assuming that they
are iid random realizations of the same stochastic process. We use
the method described in \cite{FHT2008}. For the prior mean $\mu$ we
have used the sample mean of all MRI images.

\begin{figure}[ht!]
 	\centering
 	\includegraphics[scale=0.35]{EDSRLcurveABL1black.jpg}
 	\includegraphics[scale=0.35]{EDSRLcurveApiL1black.jpg}
 	\caption{\label{LcurveSR}Risk as a function of sparsity
          $\| w \|_{0}$ of the optimal $w$ obtained
          with the $A_{B}$ design (left) and with the $A_{\pi}$ design (right).}
\end{figure}
Figure \ref{LcurveSR} shows the Pareto curves for the $A_{B}$ and $A_{\pi}$ optimal designs.
This time the curves are not as much `$L$-like' as in the 1D
example. We have chosen the points corresponding, respectively, to
$\|w\|_0=31$ and $\|w\|_0=30$ for the the $A_B$ and $A_\pi$ designs.
Hence, we use $31$ low resolution images in the $A_{B}$ design and $30$ with the $A_{\pi}$ design
to reconstruct a high resolution image.
Figures \ref{EDSRreconAB} and \ref{EDSRreconApi} show the optimal image
reconstructions for the $A_B$ and $A_\pi$ designs as well as the one
obtained using 100 low resolution images.

\begin{figure}[ht!]
 	\centering	
	\includegraphics[scale=0.35]{EDSRABall.jpg} 
 	\includegraphics[scale=0.35]{EDSRABL1.jpg}
 	\caption{\label{EDSRreconAB}Reconstruction using 100 low resolution images (left) and
	 those selected by the $A_{B}$ design (right).}
\end{figure}
\begin{figure}[ht!]
 	\centering
	\includegraphics[scale=0.35]{EDSRApiall.jpg} 
 	\includegraphics[scale=0.35]{EDSRApiL1.jpg}
 	\caption{\label{EDSRreconApi}Reconstruction using 100 low resolution images (left) and
	 those selected by the $A_{\pi}$ design (right).}
\end{figure}

As expected, the $A_{B}$ design provides a better reconstruction;
it has fewer artifacts in the background. It is also evident that
both optimal designs yield images that are very close to the ones
obtained using all $100$ low resolution images.
This illustrates the advantage of our method, it enables us to obtain
comparable results with far less data.

To show the difference between all designs, we list the relative
2-norm reconstruction errors in Tables \ref{table:SRAB} and \ref{table:SRApi}.

\begin{table}[ht!]
\caption{\label{table:SRAB}Relative errors of the reconstructed images: $A_{B}$ design}
\centering
\begin{tabular}{cc}
\hline\hline
design & relative error \\
\hline
$A_{B}$      & $4.94\%$  \\
full          & $4.95\%$  \\
\hline\hline
\end{tabular}
\end{table}

\begin{table}[ht!]
\caption{\label{table:SRApi}Relative errors of the reconstructed images: $A_{\pi}$ design}
\centering
\begin{tabular}{cc}
\hline\hline
design & relative error \\
\hline
$A_{\pi}$    & $12.15\%$  \\
full          & $9.42\%$  \\
\hline\hline
\end{tabular}
\end{table}