\chapter{Numerical methods for $E$-optimal design}
\label{Ch3}

In this chapter, we explore formulations of the Bayesian $E$-optimal design \cite{Bardow, ChalonerVerdinelli1995}, 
in which the the largest eigenvalue of the covariance matrix is minimized. For well-posed problems, this results 
in solving the optimization problem with sparsity control
\begin{eqnarray*}
\min_{w} \ \  \phi_{E}\left[\,\Sigma_\mh(W)\,\right] +\beta \|w\|_0 \quad
{\rm s.t.} \   \  0 \le w \le w_{\rm max},
\end{eqnarray*}
where
$$
\phi_{E}\left[\Sigma_{\mh}(W)\right] = \lambda_{\rm max} \left[\Sigma_{\mh}(W)\right],
$$
given the covariance matrix
$$
\Sigma_{\mh}(W) = \frac{\sigma^2}{N}(A_{n}^\top W A_{n})^{-1}.
$$

\section{The $E_{B}$ design}

For ill-posed problems, following the ideas of the $A_{B}$ and $A_{\pi}$ designs, again we propose two different designs
based on different prior information.

If the covariance matrix of the model is available, then the Bayes estimate is chosen as the model estimate, 
\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 $E_{B}$ design minimizes its Bayes risk, which is given by
\begin{equation}% \label{eq:phiab}
\phi_{E_B}(W)  = \lambda_{\rm max} \left[\,(\,A^\top WA + \Sigma_m^{-1}\,)^{-1}\,\right].
\end{equation}
The optimal sparsity is decided based on the Pareto curve in the same way as is described in Chapter~\ref{Ch2}.

\section{The $E_{\rm Tik}$ design}

When the covariance matrix is not available or its inverse is difficult to obtain, we use the Tikhonov solution
as the estimated model
\begin{equation}
\mh(w) = (A^\top W A + \alpha L^\top L)^{-1}A^\top Wd_w, % \label{eq:Tikhest},
\end{equation}
where $\alpha$ is a regularization parameter and $L$ is a chosen matrix (e.g., a discrete 
derivative operator).

In this case, the $E_{\rm Tik}$ design minimizes
\begin{equation}% \label{eq:phiab}
\phi_{E_{\rm Tik}}(W)  =  \lambda_{\rm max} \left[\,(\,A^\top WA + \alpha L^{\top}L\,)^{-1}\,\right].
\end{equation}

\section{Numerical optimization  of the $E_B$ and $E_{\rm Tik}$ designs}

The above approaches involve eigenvalue optimization techniques \cite{Lewis03themathematics}, which becomes very difficult
when the problem size is large. In general there are two ways to handle this situation: Optimize evaluate or Evaluate optimize.

In the Optimize evaluate approach, we must somehow compute the largest eigenvalue of the covariance matrix. After that,
the derivatives that are used to solve the optimization need to be found analytically. Both steps of this approach may be difficult.

Thus, we prefer the second option, the Evaluate optimize, in which we consider first approximating the largest eigenvalue of the 
covariance matrix and then evaluating the derivatives by differentiation. 

\subsection{Eigenvalue approximation}

The first step in solving the proposed designs is to approximate the largest eigenvalue of the covariance matrix.
By some simple linear algebra, it is obvious that minimizing the largest eigenvalue of the covariance matrix is equivalent
to maximizing the smallest eigenvalue of the information matrix, which is the inverse covariance matrix. Hence, the objective
functions to maximize are:
\begin{equation}
\label{EB1}
\phi_{E_B}(w) = \lambda_{\rm min} (\,A^\top WA + \Sigma_m^{-1}\,)
\end{equation}
and
\begin{equation}
\label{ETik1}
 \phi_{E_{\rm Tik}}(w) = \lambda_{\rm min} (\,A^\top WA + \alpha L^{\top}L\,).
\end{equation}

There are many eigenvalue approximation techniques in the literature, for example, methods based on power iteration
\cite{Trefethen} and SVD decomposition \cite{Nash}, etc. In this thesis we apply the inverse iteration to approximate
the smallest eigenvalue of the information matrix \cite{Haber2004}.
%------------------------------------------------------------------------------
\begin{algorithm}
\caption{Inverse iteration}
\label{alg:Inverse iteration}
\begin{algorithmic}
\STATE(1) Choose $u_{0}$ such that $\parallel u_{0} \parallel = 1$ and an integer $k$.
\STATE(2) For $j = 1, \cdots, k$, solve $Hu_{j} = \frac{1}{\sqrt{u_{j-1}^{\top}u_{j-1}}}u_{j-1}$.
\STATE(3) Set $\lambda_{\rm min} \approx \frac{1}{\sqrt{u_{k}^{\top}u_{k}}}$.
\end{algorithmic}
\end{algorithm}
%------------------------------------------------------------------------------

We start with a normalized vector $u_{0}$ and choose an appropriate integer $k$ to be the number of iterations we would 
like to perform in each optimization loop. In each inverse iteration, we solve a linear system, where the matrix $H$ is 
$A^\top WA + \Sigma_m^{-1}$ for the $E_{B}$ design and $A^\top WA + \alpha L^{\top}L$ for the $E_{\rm Tik}$ design.
Finally, we set the smallest eigenvalue $\lambda_{\rm min}$ that we are looking for to be the inverse of the norm of the 
converged eigenvector. 

Below is the matrix form of the inverse iteration process.  The original difficult eigenvalue optimization problem has 
been recast as a constrained optimization problem, which can be easily solved by standard optimization methods such as 
Steepest Descent.
$$
		\begin{array}{l}
		\displaystyle \min_{u,w} \frac{1}{2} u_{k}^{\top}u_{k} \\
		{\rm s.t.} \hspace{3mm} I(u)Bu - Eu - g = 0
		\end{array}
		$$
		where
		$$
		\begin{array}{l}
		B = {\rm diag}(H, \cdots, H) \\
		u = [u_{1}^{\top}, \cdots, u_{k}^{\top}]^{\top}, \hspace{5mm} g = [u_{0}^{\top}, 0, \cdots, 0]^{\top} \\
		\hspace{-8mm} E = \left( \begin{array}{cccc}
					0 &       &   & \\
					I & 0     &   & \\
					  & \ddots & \ddots  & \\
					  &       & I & 0
				  \end{array}
				  \right), \hspace{5mm} I(u) = \left( \begin{array}{cccc}
													    I & & & \\
													      & \sqrt{u_{1}^{\top}u_{1}}I & & \\
													      & & \ddots & \\
													      & & & \sqrt{u_{k-1}^{\top}u_{k-1}}I
													  \end{array}
														\right).												  
		\end{array}
$$

One important thing to notice here is the choice of the integer $k$, in other words, how many inverse iterations we 
want to perform in each optimization loop. According to our observation, usually setting $k$ to be between $5$ and $10$ 
should be enough. This will be demonstrated later in the experiments. 

\subsection{Evaluating the derivatives}

We use the $E_{\rm Tik}$ design as an example to show the derivation of the derivatives of the objective functions.
We rewrite the constraint in the optimization problem as
$$
\left[\begin{array}{cccc}
H & ~ & ~ & ~\\
~ & H & ~ & ~\\
~ & ~ & \ddots & ~ \\
~ & ~ & ~ & H\\
\end{array}\right] \left[\begin{array}{c}
u_{1}\\
u_{2}\\
\vdots\\
u_{k}\\
\end{array}\right] = \left[\begin{array}{cccc}
0 & ~ & ~ & ~\\
\frac{1}{\sqrt{u_{1}^{\top} u_{1}}} & 0 & ~ & ~\\
~ & \ddots & \ddots & ~ \\
~ & ~ & \frac{1}{\sqrt{u_{k - 1}^{\top} u_{k - 1}}} & 0\\
\end{array}\right] \left[\begin{array}{c}
u_{1}\\
u_{2}\\
\vdots\\
u_{k}\\
\end{array}\right] + \left[\begin{array}{c}
u_{0}\\
0 \\
\vdots\\
0\\
\end{array}\right]
,
$$
which gives us the linear system
$$
\left\{ \begin{array}{rcl}
A^{\top} \diag(Au_{1}) + H \frac{\partial u_{1}}{\partial w} & = & 0\\
A^{\top} \diag(Au_{2}) + H \frac{\partial u_{2}}{\partial w} & = & [(u_{1}^{\top}u_{1})^{-1/2} - (u_{1}^{\top}u_{1})^{-3/2} u_{1}u_{1}^{\top}] \frac{\partial u_{1}}{\partial w}\\
\vdots & ~ & ~\\
A^{\top} \diag(Au_{k}) + H \frac{\partial u_{k}}{\partial w} & = & [(u_{k-1}^{\top}u_{k-1})^{-1/2} - (u_{k-1}^{\top}u_{k-1})^{-3/2} u_{k-1}u_{k-1}^{\top}] \frac{\partial u_{1}}{\partial w}\\
\end{array}\right. .
$$
Defining
$$
B_{i} = (u_{i}^{\top}u_{i})^{-1/2} - (u_{i}^{\top}u_{i})^{-3/2} u_{i}u_{i}^{\top},
$$
we obtain
$$
\left[\begin{array}{cccc}
H & ~ & ~ & ~\\
-B_{1} & H & ~ & ~\\
~ & \ddots & \ddots & ~\\
~ & ~ & -B_{k-1} & H\\
\end{array} \right]
\left[ \begin{array}{c}
\frac{\partial u_{1}}{\partial w}\\
\frac{\partial u_{2}}{\partial w}\\
\vdots\\
\frac{\partial u_{k}}{\partial w}\\
\end{array} \right] = -\left[ \begin{array}{c}
A^{\top}\diag(Au_{1})\\
A^{\top}\diag(Au_{2})\\
\vdots\\
A^{\top}\diag(Au_{k})\\
\end{array} \right].
$$
Thus, the derivative of the objective function is
$$
\begin{array}{ll}
\displaystyle \frac{\partial}{\partial w}\left( \frac{1}{2} u_{k}^{\top}u_{k} \right) \\
= -\left[ \begin{array}{c}
A^{\top}\diag(Au_{1})\\
A^{\top}\diag(Au_{2})\\
\vdots\\
A^{\top}\diag(Au_{k})\\
\end{array} \right]^{\top} \left[ \begin{array}{ccccc}
H & -B_{1}^{\top} & ~ & ~ & ~\\
~ & H & -B_{2}^{\top} & ~ & ~\\
~ & ~ & \ddots & \ddots & ~\\
~ & ~ & ~ & H & -B_{k - 1}\\
~ & ~ & ~ & ~ & H\\
\end{array} \right]^{-1} \left[ \begin{array}{c}
0 \\
\vdots \\
0 \\
u_{k}\\
\end{array} \right].
\end{array}
$$

\section{Numerical experiments}

In this section, we will show the performance of the two $E$-optimal designs through first, the 1D magnetotelluric example
that we have used in the $A$-optimal designs and then a borehole ray tomography example.

\subsection{An ill-posed 1D magnetotelluric example}

We apply both $E$-optimal designs to the same 1D magnetotelluric example that we have used for the $A$-optimal designs.
For the purpose of comparison, all experiment and parameter setting as exactly the same as is in Chapter~\ref{Ch2}.

\subsubsection{Quality of eigenvalue approximation}

Before we start the experiments, we first test how well the eigenvalue approximation works for both designs.
Tables \ref{1DeigEB} and \ref{1DeigETik} show the difference between the true eigenvalues and the estimated ones 
for different values of $k$ with their corresponding CPU times that were needed to compute the eigenvalue approximation. 

\begin{table}[ht!]
\caption{\label{1DeigEB} Eigenvalue approximation for the $E_{B}$ design}
\centering
\begin{tabular}{c|c|c}
\hline\hline
$k$ & $\lambda_{\rm app} - \lambda_{\rm true}$ & CPU time \\
\hline
3   & $4.0500e-2$ & $3.120020e-2$ \\
5   & $1.3190e-4$ & $6.240040e-2$ \\
8   & $0$         & $6.240040e-2$ \\
\hline\hline
\end{tabular}
\end{table}

\begin{table}[ht!]
\caption{\label{1DeigETik} Eigenvalue approximation for the $E_{\rm Tik}$ design}
\centering
\begin{tabular}{c|c|c}
\hline\hline
$k$ & $\lambda_{\rm app} - \lambda_{\rm true}$ & CPU time \\
\hline
3   & $1.0800e-2$ & $0$ \\
6   & $1.0000e-5$ & $3.120020e-2$ \\
9   & $0$         & $6.240040e-2$ \\
\hline\hline
\end{tabular}
\end{table}

Based on the tables, we choose $k = 8$ for the $E_{B}$ design and $k = 9$ for the $E_{\rm Tik}$ design.
The approximation of eigenvalue is very accurate while the time cost is still reasonable.

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

First we show results from the $E_{B}$ design. Figure \ref{1DLcurveEB} shows the L-curve, in which the corner is located
at $\|w\|_0=20$.
\begin{figure}[ht!]
 	\centering
 	\includegraphics[scale=0.35]{LcurveEB.jpg}
 	\caption{\label{1DLcurveEB} The panel shows
the risk as a function of sparsity ${\rm nnz}(w)=\parallel
w\parallel_{0}$  of the optimal $w$ obtained with the $E_{B}$ design.}
\end{figure}

Figure \ref{1DEBL1} 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 20 selected experiments.

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

For comparison, we show the
model estimates obtained using all $100$ experiments and also the
estimates based on $20$ equally-spaced naively chosen
experiments in Figure \ref{1DEBnaive}. It is obvious that the $E_B$ design provides the best results.

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

Table \ref{tab:EBL1} shows the relative errors ($\|m_{\rm true}-\widehat{m}\|/\|m_{\rm true}\|$)
of the reconstructions. Again the estimates obtained using all $100$ experiments are worse 
than those using only the $20$ optimal experiments.

\begin{table}[ht!]
\caption{\label{tab:EBL1} Relative errors of the model reconstructions obtained with the $E_{B}$ design}
\centering
\begin{tabular}{c|c|c|c}
\hline\hline
design & error (noise 1) & error (noise 2) & error (noise 3) \\
\hline
$E_{B}$   & $7.32\%$  & $13.35\%$ & $12.44\%$ \\
all       & $10.29\%$ & $20.12\%$ & $19.26\%$ \\
naive     & $11.38\%$ & $24.40\%$ & $22.38\%$ \\
\hline\hline
\end{tabular}
\end{table}

\subsubsection{The $E_{\rm Tik}$ design}

\begin{figure}[ht!]
 \centering
 \includegraphics[scale=0.35]{LcurveEpi.jpg}
 \caption{\label{1DLcurveETik}The panel shows the risk as a function of sparsity 
${\rm nnz}(w)=\parallel w \parallel_{0}$ of the optimal $w$ obtained with the $E_{\rm Tik}$ design.}
\end{figure}

In Figure \ref{1DLcurveETik}, we choose the corner at $\|w\|_0=21$. This defines the optimal
solution used in the example.

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

\begin{figure}[ht!]
 	\centering
 	\includegraphics[scale=0.35]{mEpiall.jpg}
 	\includegraphics[scale=0.35]{mEpinaive.jpg}
 	\caption{\label{1DETiknaive}Model estimates obtained with the
          $E_{\rm Tik}$ design using all $100$ experiments (left) and the equally-spaced naive design (right).}
\end{figure}

Figure \ref{1DETikL1} shows the model estimates for three different
noise realizations obtained using the $\ell^1$ approach. Figure \ref{1DETiknaive} shows the estimates
obtained using all $100$ experiments and using $21$ equally-spaced naively chosen experiments.
The relative errors of the estimates are shown in Table
\ref{tab:ETikL1}. 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 $E_B$ design. 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:ETikL1} Relative errors of the model estimates
  obtained with the $E_{\rm Tik}$ ($\ell^{1}$) design}
\centering
\begin{tabular}{c|c|c|c}
\hline\hline
design & error (noise 1) & error (noise 2) & error (noise 3) \\
\hline
$E_{\rm Tik}$ $\ell^{1}$ & $40.05\%$ & $35.41\%$ & $46.70\%$ \\
all                  & $47.00\%$ & $41.07\%$ & $53.60\%$ \\
naive                & $46.81\%$ & $39.47\%$ & $53.71\%$ \\
\hline\hline
\end{tabular}
\end{table}

So far, we have worked on both the $A$-optimal designs and $E$-optimal designs. Although their results are comparable based on 
our results, in general the choice is application dependent. We know that the $A$-optimal design minimizes the variance in 
an average sense while the $E$-optimal design minimizes the situation that has the largest variance, in other words, 
it minimizes the worst case. The performance depends on whether your example is closer to the average case or to the 
worse case. For example, if a doctor wants to minimize the largest risk that a tumor could cause, then $E$-optimal design 
would be a good choice. On the other hand, if a car company wants to decide how powerful the air bag should be. 
Then, in order to fit people's needs on average, $A$-optimal design is needed.

\subsection{A borehole ray tomography example}

Next we apply our approaches to a borehole ray tomography example \cite{HaberHoreshTenorio08} that is often used to illustrate purposes in geophysical 
inverse problems. The purpose of borehole ray tomography is to determine the slowness, which is the inverse of the velocity 
of a medium. Sources and receivers are placed along the boreholes and travel times from sources to receivers are recorded.

In our experiment, the medium is discretized by the square region $[0, 1] \times [0, 1]$ and boreholes are covering 
both sides of the region. Totally we have $1600$ rays and our goal of an experimental design is to choose the 
optimal placement of sources and receivers. 

Since this is a real problem, we do not have the true covariance matrix, we will apply the example only to the $E_{\rm Tik}$ design,
for which the discretization of the gradient operator is used as a smoothing matrix.

\subsubsection{Quality of eigenvalue approximation}

Table \ref{RayeigETik} show the difference between the true eigenvalues and the estimated ones 
for different values of $k$ with their corresponding CPU times that were needed to compute the eigenvalue approximation. 

\begin{table}[ht!]
\caption{\label{RayeigETik} Eigenvalue approximation for the $E_{\rm Tik}$ design}
\centering
\begin{tabular}{c|c|c}
\hline\hline
$k$  & $\lambda_{\rm app} - \lambda_{\rm true}$ & CPU time \\
\hline
4    & $1.5406e-5$ & $9.874863e$ \\
7    & $2.5856e-6$ & $1.647371e+1$ \\
10   & $8.3150e-7$ & $2.266695e+1$ \\
\hline\hline
\end{tabular}
\end{table}

For solving the problem accurately, we choose $k = 10$. 

\subsubsection{The $E_{\rm Tik}$ design}

\begin{figure}[ht!]
 \centering
 \includegraphics[scale=0.35]{LcurveL0.jpg}
 \caption{\label{RayLcurveETik}The panel shows the risk as a function of sparsity 
${\rm nnz}(w)=\parallel w \parallel_{0}$ of the optimal $w$ obtained with the $E_{\rm Tik}$ design.}
\end{figure}

Figure \ref{RayLcurveETik} shows the Pareto curves for the $E_{\rm Tik}$ optimal design. 
The corner happens where $\|w\|_0 = 746$. Figures \ref{RayETik} show the optimal raypaths
reconstructions for the $E_{\rm Tik}$ design as well as the one obtained using all $1600$ raypaths.

\begin{figure}[ht!]
 	\centering
	\includegraphics[scale=0.35]{Raymtrue.jpg} \\
 	\includegraphics[scale=0.35]{Raym0.jpg}
 	\includegraphics[scale=0.35]{Raymall.jpg}
 	\caption{\label{RayETik}The true image (top), reconstructions using $1600$ raypaths (bottom left) and
	 those selected by the $E_{\rm Tik}$ design (bottom right).}
\end{figure}

It is evident that our optimal design yields images that are very close in quality to the ones
obtained using all $1600$ raypaths.

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

\begin{table}[ht!]
\caption{\label{table:RayETik}Relative errors of the reconstructed images: $E_{\rm Tik}$ design}
\centering
\begin{tabular}{cc}
\hline\hline
design & relative error \\
\hline
$E_{\rm Tik}$    & $69.63\%$  \\
full          & $65.06\%$  \\
\hline\hline
\end{tabular}
\end{table}

Based on \ref{RayETik} and \ref{table:RayETik}, we observe that the reconstruction from our optimal design
provides results almost as good as using all $1600$ raypaths. However, we have cut out more than half of the cost.