\chapter{Optimal design for regularization}
\label{Ch4}

In this chapter, we talk about finding an optimal $\ell^{2}$ regularization that plays important roles in imaging science.
We consider a discrete linear ill-posed problem of the form
\begin{eqnarray}
\label{eq:linprob}
A x + \epsilon = b,
\end{eqnarray}
where $A:R^{m} \rightarrow R^{n}$ is a discretization of some linear (typically integral)
operator and $b$ is the observed data. The random vector $\epsilon$ is the noise contained in the data, which is 
assumed to be iid Gaussian with standard deviation $\sigma$ and $0$ mean.

The general goal is to recover the model $x$ from the observed noisy data $b$.
However, since the problem is ill-posed it is not likely to recover $x$ just from the data itself directly.
Thus, we recover $x$ by using a Tikhonov-like regularization and solve the optimization problem
\begin{eqnarray}
\label{eq:opt}
\min\ \  \hf \|A x - b\|^{2} + \alpha R(x),
\end{eqnarray}
where $R(x)$ is a regularization functional that penalizes unwanted solutions and $\alpha$
is a regularization parameter that is chosen to control the balance between the mis-fitting term
and the regularization term so that the solution does not over-fit the data.

For linear ill-posed problems, a regularization $R(x) = \hf \| Lx \|^{2}$ is commonly 
used \cite{Phillips, Tikhonov1977}, where $L$ is a discretization of a differential
operator based on smoothness of the solution.
Obviously, the Tikhonov regularization functional is quadratic and therefore, the 
solution to the problem can be easily obtained by solving the linear system
\begin{eqnarray}
\label{eq:optsol}
(A^{\top}A + \alpha L^{\top} L ) x = A^{\top} b,
\end{eqnarray} 
using either direct, but more typically, iterative methods such as the Preconditioned Conjugate Gradient method
and Lanczos tridiagonalization method \cite{Greenbaum1997, Hnetynkova}.

\section{An optimal regularization operator}
\label{sec41}

In order to develop an optimal regularization operator we need to define an optimality criteria.
In the following we develop such criteria. As we see next, the criteria heavily depends on the assumptions or a-priori
information known. Different assumptions lead to a different regularization operators.

Consider a solution to the quadratic regularization problem
$$
\hat{x} = (A^{\top}A + \alpha L^{\top}L)^{-1}A^{\top}b.
$$
where $b$ is given by equation~\eqref{eq:linprob}.
If we are given the ``true'' solution $x^{t}$, we may ask how well $\hat{x}$ reconstructs $x^{t}$.
The answer to this question is the well known Mean Square Error:
\begin{eqnarray*}
\label{eq:EbiasVariance}
\MSE = {\bf E} \parallel \hat{x} -x^{t} \parallel^{2} & = & \alpha ^{2} \parallel (A^{\top}A + \alpha L^{\top}L)^{-1}L^{\top} L x^{t} \parallel^{2} \\
											 ~ & ~ & + \sigma^{2}{\rm trace}\left[A(A^{\top}A + \alpha L^{\top}L)^{-2}A^{\top} \right].
\end{eqnarray*}
where the expected value is on the noise. The first term which depends on $x^{t}$ is referred to as the square of the bias and the second term is
referred to as  the variance. The decomposition of the MSE into the bias and variance is a major point in our discussion.
The bias is required if stable solutions to the problem are desired. This implies that
even for problems with no noise, the recovered solution $\widehat x$ will differ from the true solution $x^{t}$. An important
question is therefore, how should we decrease the error in recovery?

\subsection{The first complication: MSE is dependent on the true solution}

One could reduce the bias by decreasing $\alpha$ but, as is well documented, this will increase the variance.
Thus, the goal is to introduce the ``right'' $L$ that leads to the ``right'' bias, that is, a bias that will lead us 
closer to the true solution $x^{t}$ and decrease the overall MSE. 

This implies that we need some information on $x^{t}$, which leads to the first complication in finding the optimal $L$. Assume that one chose an $L$ that minimizes the MSE.
Since the MSE depends on the true solution $x^{t}$, minimizing the MSE that depends on $x^{t}$ may lead to a large MSE for
different ``true'' models. To overcome this disadvantage, a number of approaches can be used to eliminate the 
dependency of the MSE with respect to $x^{t}$.
In the following we explore two different scenarios where different information on $x^{t}$ is given.

\begin{itemize}
\item[A1:] If the covariance of the true $x^{t}$ exists and is known, then it is possible to use the covariance
to define an average MSE. The average MSE in this case, is identical to the so called Bayesian risk.
\item[A2:] If some training models $X = [x^{t}_{1},\ldots,x^{t}_{s}]$ are given, then
it is possible to obtain an empirical estimation of the average MSE. This case is identical to the empirical Bayesian risk.
% \item[A3] In some cases only upper bound on the norm of $x^{t}$ is given. In these cases one may want to use a minimax 
% estimator and to look at the maximal MSE.
\end{itemize}
We now treat each of these cases and discuss the optimization problems that stem from each assumption.

\subsubsection{The covariance design}

If the mean $\mu$ and the covariance matrix $\Sigma$ of the true model are known, then we can replace the norm of the bias 
by its average.
% Since the bias is quadratic in $x^{t}$ we have that
% $$  {\bf E}\left( \parallel C(L)^{-1}L^{\top} L x^{t} \parallel^{2} \right) = {\sf trace}
% \left( L^{\top}LC(L)^{-1} \Sigma C(L)^{-1}L^{\top}L \right) $$
% where
% $$ C(L) = A^{\top}A + \alpha L^{\top}L $$ and $\Sigma$ is the covariance matrix of the model. % and $M = \frac{1}{s} \displaystyle \sum_{j=1}^{s}x^{t}{x^{t}}^{\top}$.
An optimal solution, in this case, is the one that minimizes the average MSE. This leads to the following optimization problem
\begin{equation}
\begin{array}{l}
\displaystyle \min_{L} \hspace{2mm} \alpha ^{2} \parallel (A^{\top}A + \alpha L^{\top}L)^{-1}L^{\top} L \mu \parallel^{2} + \\
{\rm trace}\left[(A^{\top}A + \alpha L^{\top}L)^{-1}(\alpha^{2}L^{\top}L\Sigma L^{\top}L + \sigma^{2}A^{\top}A)(A^{\top}A + \alpha L^{\top}L)^{-1} \right].
\end{array}
\end{equation} 
At this point, it may seem surprising but $L$ can be found analytically.
For simplicity, we set the mean $\mu$ to be zero.
Setting $B = L^{\top}L$, $C = A^{\top}A + \alpha L^{\top}L$ and combining the traces we have
\begin{equation*}
\begin{array}{l}
{\rm AMSE}   =  {\rm trace}\left(C^{-1}(\alpha^{2}B\Sigma B + \sigma^{2}A^{\top}A)C^{-1}\right) = \\
           \displaystyle \sum^{m}_{i=1} {e_{i}}^{\top} \left[C^{-1}(\alpha^{2}B\Sigma B + \sigma^{2}A^{\top}A)C^{-1}\right] e_{i},
\end{array}
\end{equation*}
where $e_{i}, i = 1,\cdots, m$ are the standard unit basis vectors. 
In order to minimize the AMSE, we set its derivative with respect to $B$ to be $0$ as follows:
\begin{equation*}
\begin{array}{l}
\frac{\partial {\rm AMSE}}{\partial B} = \\
\displaystyle \sum^{m}_{i=1} \left[ -\alpha C^{-1}e_{i}{e_{i}}^{\top}C^{-1}(\alpha^{2}B\Sigma B + \sigma^{2}A^{\top}A)C^{-1}
                                        + \alpha^{2}C^{-1}e_{i}{e_{i}}^{\top}C^{-1}B\Sigma \right] = 0.
\end{array}
\end{equation*}
It is straight forward to verify that, if we chose $B$ to be
\begin{equation*}
B = \frac{\sigma^{2}}{\alpha}\Sigma^{-1} = L^{\top}L,
\end{equation*}
then the gradient of the ${\rm AMSE}$ vanishes. That is, the ``optimal'' $L^{\top}L$ is a proportion of the inverse covariance of the models.
This should not come as a surprise. The average MSE can be interpreted as the Bayesian risk and thus, if the model
is Gaussian with mean zero and covariance matrix $\Sigma$, then the Maximum A Posterior Estimate (MAP) also yields the minimal 
Bayes risk. The surprising result (at least for us) is the fact that, for {\em any} distribution of models with zero mean
and covariance matrix $\Sigma$, even for distributions that are far from Gaussian, the inverse covariance is the best regularization matrix.

Assume that we are given the covariance matrix and we would like to use its inverse for the solution of the 
problem. For small scale problems, this is straightforward. However, for large scale problems, working
with a general covariance matrix and its inverse is highly nontrivial. For example, in medical and geophysical
imaging, $x^{t}$ can be a vector with millions of entries. For these applications, working with a dense covariance matrix with $10^{12}$
entries makes computations infeasible unless some special structure exists. For example, one can assume that the
covariance matrix is space invariant which leads to a simple $\Sigma$. However, this assumption is rather unrealistic
and even if it is true, estimating $\Sigma^{-1}$ can be highly nontrivial.

\subsubsection{The training design}

When we do not have the covariance matrix, we use a set of related models which are used as training references.
These models are chosen for applications in different situations. For example, if the goal is to deblur some MRI 
images then $X$ can be chosen as a set of clean MRI examples.
In this case, the MSE to minimize is
\begin{equation}
\label{MSEtrain}
{\rm AMSE} = {\rm trace}\left[(A^{\top}A + \alpha L^{\top}L)^{-1}(\alpha^{2}L^{\top}LML^{\top}L + \sigma^{2}A^{\top}A)(A^{\top}A + \alpha L^{\top}L)^{-1} \right],
\end{equation}
where $M = \frac{1}{s} \displaystyle \sum_{j=1}^{s}x_{j}^{t}{x_{j}^{t}}^{\top}$ and $s$ is the number of training models that we use.

\subsection{The second complication: The computational complexity}

As we have stated above, using a dense regularization matrix for large scale problems will be very expensive in
computation. The weakness of the analytic optimal $L$ strengthens the motivation for our work.

To overcome the second disadvantage, the problem of the computational complexity, we can choose $L$ to have some 
structure that is easy to compute with. We return to the average MSE and minimize it over all $B = L^{\top}L$ assuming that $L$ has 
some given sparsity structure. Let ${\cal S}$ be the set of matrices with a specific sparsity structure and solve
the following {\em constrained} optimization problem
\begin{eqnarray}
\label{eq:OptLC}
\min_{L} && {\rm AMSE} \\
\nonumber
{\rm s.t.} && L \in {\cal S}.
\end{eqnarray} 
Here are a few comments:
\begin{itemize}
\item It is rather clear that, in general, the solution in this case is not that $\frac{\sigma^{2}}{\alpha}\Sigma^{-1} = L^{\top}L$ 
since $\Sigma^{-1}$
may not posses the appropriate sparsity pattern. Also, simple ``sparsification'' of $\Sigma^{-1}$ (that is, projecting
$\Sigma^{-1}$ into the constraint set) is in general, not the optimal solution.
\item As before, one could set $B=L^{\top}L$ and solve for $B$ directly, under the constraint that $B$ is symmetric
positive and semidefinite (PSD). However, note that, for this case, the ${\rm AMSE}$ function is non-convex. Furthermore, 
even if the function was convex, working with the constraint that $B$ is PSD is difficult for large scale problems. In fact,
while there is a body of work that deals with estimating the inverse of the covariance matrix, we are not aware of any papers
that work well in estimating the inverse covariance matrix when the size of the problem is very large.
% The only method that seems to be relevant for the size of problems we tackle here is presented in  \cite{SamBurer} where one computes
% the Cholesky factors of the inverse covariance. 
\item While being optimal is a novel goal, for most practical applications a significant improvement over existing regularization
methods will suffice. Thus, although we do not have a convex problem, if it is possible to obtain useful solutions, this will be welcome
for many of the applications we aim for.
\end{itemize}

In the rest of this chapter, we focus on solving the first complication, while the complication of the computational complexity
will be discussed in details in the next chapter.

\section{Numerical optimization of the optimal regularization}

In this section, we discuss numerical optimization methods of the optimal regularization.
Different approaches are explored in order to solve individual problems that arise from different designs.

\subsection{Matrix-based derivative techniques}

We start with the derivatives for the proposed designs. Most literature approaches optimization problems by minimizing vectors, 
while it is not common to minimize with respect to matrices. Hence, we will take some effort to review some matrix-based
derivative techniques \cite{Matrixcook} first. 

Here we list a couple basic derivative rules that we have used in deriving our derivatives. Let $N$ and $P$ be matrix variables
and $c$ be a constant. $i, j, k, l $ refer to different matrix entry index.
\begin{equation}
\begin{array}{rcl}
\partial (cN) &  = & c \partial (N) \\
\partial (N + P) & = & \partial (N) + \partial (P) \\
\partial (NP) & = & N \partial (P) + P \partial (N) \\
\partial (N^{-1}) & = & - N^{-1} \partial (N) N^{-1} \\
\partial (N^{\top}) & = & (\partial N)^{\top} \\
\frac{\partial N_{k,l}}{\partial N_{i, j}} & = & \delta_{i, k}\delta_{l, j}
\end{array}
\end{equation}

In particular, the following matrix differentiation rules have been very handy. Let $q$ be a vector variable and
$f, g$ be two constant vectors. $J^{ij}$ is a matrix with $1$ on the $ij^{\rm th}$ entry and $0$ elsewhere.
\begin{equation}
\begin{array}{rcl}
\frac{\partial (N^{-1})}{\partial q} & = & - N^{-1} \frac{\partial N}{\partial q} N^{-1} \\
\frac{\partial (N^{-1})_{k, l}}{\partial N_{i, j}} & = & -(N^{-1})_{k, i}(N^{-1})_{j, l} \\
\frac{\partial f^{\top}N^{-1}g}{\partial N} & = & -N^{-\top}fg^{\top}N^{-\top} \\
\frac{\partial f^{\top}Ng}{\partial N} & = & fg^{\top} \\
\frac{\partial f^{\top}N^{\top}g}{\partial N} & = & gf^{\top} \\
\frac{\partial f^{\top}Nf}{\partial N} & = & ff^{\top} \\
\frac{\partial N}{\partial N_{i, j}} & = & J^{ij} \\
\frac{\partial f^{\top}N^{\top}Ng}{\partial N} & = & N(fg^{\top} + gf^{\top})
\end{array}
\end{equation}

\subsection{The covariance design approach}

By applying the Stochastic trace approximation, the original covariance design can be rewritten as
\begin{equation}
\begin{array}{rcl}
\label{traceMSEcov}
{\rm AMSE_{Cov}} & = & \alpha ^{2} \parallel (A^{\top}A + \alpha L^{\top}L)^{-1}L^{\top} L \mu \parallel^{2} \\
~ & + & v^{\top}(A^{\top}A + \alpha L^{\top}L)^{-1}(\alpha^{2}L^{\top}L\Sigma L^{\top}L + \sigma^{2}A^{\top}A)(A^{\top}A + \alpha L^{\top}L)^{-1}v
\end{array}
\end{equation}
where $v$ is a random vector of $1$ and $-1$ with equal distribution.

We discuss the derivatives of the two terms in the AMSE with respect to the regularization matrix $L$ separately.
We denote the first quadratic term to be $U$ and the second term $V$. Hence, the objective functional can be expressed as
$$
{\rm AMSE} = \alpha^{2}U + V.
$$
We rewrite the $U$ as
$$
U = \mu^{\top}L^{\top}L(A^{\top}A + \alpha L^{\top}L)^{-2}L^{\top}L\mu.
$$
The derivation of its derivative is as follows:
\begin{eqnarray}
\frac{\partial U}{\partial L} & = & \label{dUdL1} \frac{\partial [\mu^{\top}L_{nf}^{\top}L(A^{\top}A + \alpha L^{\top}L)^{-1}(A^{\top}A + \alpha L^{\top}L)^{-1}L^{\top}L\mu]}{\partial L} \\
                              & + & \label{dUdL2} \frac{\partial [\mu^{\top}L^{\top}L_{nf}(A^{\top}A + \alpha L^{\top}L)^{-1}(A^{\top}A + \alpha L^{\top}L)^{-1}L^{\top}L\mu]}{\partial L} \\
                              & + & \label{dUdL3} \frac{\partial [\mu^{\top}L^{\top}L(A^{\top}A + \alpha L_{nf}^{\top}L)^{-1}(A^{\top}A + \alpha L^{\top}L)^{-1}L^{\top}L\mu]}{\partial L} \\
                              & + & \label{dUdL4} \frac{\partial [\mu^{\top}L^{\top}L(A^{\top}A + \alpha L^{\top}L_{nf})^{-1}(A^{\top}A + \alpha L^{\top}L)^{-1}L^{\top}L\mu]}{\partial L} \\
                              & + & \label{dUdL5} \frac{\partial [\mu^{\top}L^{\top}L(A^{\top}A + \alpha L^{\top}L)^{-1}(A^{\top}A + \alpha L_{nf}^{\top}L)^{-1}L^{\top}L\mu]}{\partial L} \\
                              & + & \label{dUdL6} \frac{\partial [\mu^{\top}L^{\top}L(A^{\top}A + \alpha L^{\top}L)^{-1}(A^{\top}A + \alpha L^{\top}L_{nf})^{-1}L^{\top}L\mu]}{\partial L} \\
                              & + & \label{dUdL7} \frac{\partial [\mu^{\top}L^{\top}L(A^{\top}A + \alpha L^{\top}L)^{-1}(A^{\top}A + \alpha L^{\top}L)^{-1}L_{nf}^{\top}L\mu]}{\partial L} \\
                              & + & \label{dUdL8} \frac{\partial [\mu^{\top}L^{\top}L(A^{\top}A + \alpha L^{\top}L)^{-1}(A^{\top}A + \alpha L^{\top}L)^{-1}L^{\top}L_{nf}\mu]}{\partial L}.
\end{eqnarray}
where we treat $L$ as fixed constant and $L_{nf}$ as non-fixed variable. Let's discuss these equations one by one. \\
For (\ref{dUdL1}) and (\ref{dUdL8}), we have
$$
(\ref{dUdL1}) = (\ref{dUdL8}) = L(A^{\top}A + \alpha L^{\top}L)^{-2}L^{\top}L\mu \mu^{\top}.
$$
For (\ref{dUdL2}) and (\ref{dUdL7}),
$$
(\ref{dUdL2}) = (\ref{dUdL7}) = L\mu (A^{\top}A + \alpha L^{\top}L)^{-2}L^{\top}L\mu^{\top}.
$$
For (\ref{dUdL3}) and (\ref{dUdL6}), we have
\begin{eqnarray*}
(\ref{dUdL3}) & = & \mu^{\top}L^{\top}L \frac{\partial (A^{\top}A + \alpha L_{nf}^{\top}L)^{-1}}{\partial L}(A^{\top}A + \alpha L^{\top}L)^{-1}L^{\top}L\mu \\
			  & = & - \mu^{\top}L^{\top}L(A^{\top}A + \alpha L^{\top}L)^{-1}\frac{\partial (A^{\top}A + \alpha L_{nf}^{\top}L)}{\partial L} \\
& ~ & (A^{\top}A + \alpha L^{\top}L)^{-1}(A^{\top}A + \alpha L^{\top}L)^{-1}L^{\top}L\mu \\
			  & = & - \alpha \frac{\partial [\mu^{\top}L^{\top}L(A^{\top}A + \alpha L^{\top}L)^{-1}L_{nf}^{\top}L(A^{\top}A + \alpha L^{\top}L)^{-2}L^{\top}L\mu ]}{\partial L} \\
			  & = & - \alpha L(A^{\top}A + \alpha L^{\top}L)^{-2}L^{\top}L\mu (A^{\top}A + \alpha L^{\top}L)^{-1}L^{\top}L\mu^{\top} = (\ref{dUdL6}).
\end{eqnarray*}
For (\ref{dUdL4}) and (\ref{dUdL5}), we have
\begin{eqnarray*}
(\ref{dUdL4}) & = & \mu^{\top}L^{\top}L \frac{\partial (A^{\top}A + \alpha L^{\top}L_{nf})^{-1}}{\partial L}(A^{\top}A + \alpha L^{\top}L)^{-1}L^{\top}L\mu \\
			  & = & - \mu^{\top}L^{\top}L(A^{\top}A + \alpha L^{\top}L)^{-1}\frac{\partial (A^{\top}A + \alpha L^{\top}L_{nf})}{\partial L} \\
& ~ & (A^{\top}A + \alpha L^{\top}L)^{-1}(A^{\top}A + \alpha L^{\top}L)^{-1}L^{\top}L\mu \\
			  & = & - \alpha \frac{\partial [\mu^{\top}L^{\top}L(A^{\top}A + \alpha L^{\top}L)^{-1}L^{\top}L_{nf}(A^{\top}A + \alpha L^{\top}L)^{-2}L^{\top}L\mu ]}{\partial L} \\
			  & = & - \alpha L(A^{\top}A + \alpha L^{\top}L)^{-1}L^{\top}L\mu (A^{\top}A + \alpha L^{\top}L)^{-2}L^{\top}L\mu^{\top} = (\ref{dUdL5}).
\end{eqnarray*}
Hence, we obtain
\begin{eqnarray*}
\frac{\partial U}{\partial L} & = & 2L(A^{\top}A + \alpha L^{\top}L)^{-2}L^{\top}L\mu \mu^{\top} \\
							  & + & 2L\mu (A^{\top}A + \alpha L^{\top}L)^{-2}L^{\top}L\mu^{\top} \\
							  & - & 2 \alpha L(A^{\top}A + \alpha L^{\top}L)^{-2}L^{\top}L\mu (A^{\top}A + \alpha L^{\top}L)^{-1}L^{\top}L\mu^{\top} \\
							  & - & 2 \alpha L(A^{\top}A + \alpha L^{\top}L)^{-1}L^{\top}L\mu (A^{\top}A + \alpha L^{\top}L)^{-2}L^{\top}L\mu^{\top}.
\end{eqnarray*}
With respect to the second term $V$, we define
$D = \alpha^{2}L^{\top}L\Sigma L^{\top}L + \sigma^{2}A^{\top}A$ and yield
$$
\grad_{L} V = \grad_{L} (v^{\top}C^{-1}DC^{-1}v).
$$
Setting 
$$
z = C^{-1}v \quad {\rm and} \quad y = C^{-1}Dz
$$
yields
$$
V = z^{\top}Dz.
$$
Following the same derivative techniques as described above, we obtain
$$
\grad_{L} V = - 2 \alpha L (yz^{\top} + zy^{\top}) + 2 \alpha^2 L (\Sigma L^{\top}Lzz^{\top} + zz^{\top}L^{\top}L\Sigma ). 
$$
Finally, the derivative of the objective functional is expressed as
$$
\frac{\partial {\rm AMSE}_{\rm Cov}}{\partial L} = \alpha^2 \frac{\partial U}{\partial L} + \frac{\partial V}{\partial L}.
$$

\subsection{The training design approach}

The training design can be rewritten as follows again using the Stochastic trace approximation
\begin{equation}
\label{MSEtrain}
{\rm AMSE} = v^{\top}(A^{\top}A + \alpha L^{\top}L)^{-1}(\alpha^{2}L^{\top}LML^{\top}L + \sigma^{2}A^{\top}A)(A^{\top}A + \alpha L^{\top}L)^{-1}v.
\end{equation}
This formulation looks a lot like the second term $V$ in the covariance design. Therefore, we will give the derivative directly.
$$
\grad_{L} {\rm AMSE} = - 2 \alpha L (y'z^{\top} + zy'^{\top}) + 2 \alpha^2 L (ML^{\top}Lzz^{\top} + zz^{\top}L^{\top}LM ),
$$
where $y' = C^{-1}D'z$ and $D' = \alpha^{2}L^{\top}LML^{\top}L + \sigma^{2}A^{\top}A$.

\subsection{Numerical Optimization}

The numerical optimization is carried out by the steepest descent method. One important thing to notice here is that, 
when checking if the derivatives are properly derived, we need to check if the relation
$$
{\rm AMSE} (L + hL_{s}) = {\rm AMSE} (x) + (\frac{\partial {\rm AMSE}}{\partial L}, hL_{s}) + o(h^{2})
$$
holds. In this case since both $\frac{\partial {\rm AMSE}}{\partial L}$ and $hL_{s}$ are matrices with the same size as $L$,
the normal inner product will give us another matrix which does not add up to the AMSE as it is a scalar.
Thus, we introduce the so called standard inner product for the two matrices \cite{Meyer}
$$
(\frac{\partial {\rm AMSE}}{\partial L}, hL_{s}) = {\rm trace}[(\frac{\partial {\rm AMSE}}{\partial L})^{\top}(hL_{s})].
$$