\chapter{Optimal sparse regularization}
\label{Ch5}

In the previous chapter, we have discussed the formulation and derivative techniques for finding the optimal regularization
matrix. As you may see, computing the derivatives is very time consuming since PCG is used everywhere in the computation. 
Also, considering the structure of the image, in some cases, it is not really necessary to find $L$ with the dense pattern. 
Thus, we need to add certain sparsity constraints to $L$ in order to save the computational cost.

In this chapter, we propose three different sparsity patterns for $L$: The local diagonal pattern, the $\ell^{1}$ norm
pattern and the Kronecker product pattern.

\section{The local diagonal pattern}

The first one we call the local diagonal pattern, which means we only use information from the $k$-diagonals that are close 
to the main diagonal. In 1D problems, this indicates that there are only nonzero entries on the main, sub and super diagonals. 
In 2D problems, we choose the pattern of the Laplacian matrix as our structure for $L$, in other words, the nonzero pattern 
of $L$ is just $5$ diagonals. The reason for doing this will be easy to explain on an image. Images consist of pixels and 
each pixel is often more related to the pixels that surround it. Thus we choose the $5$ diagonals which exactly represent 
the $4$ pixels in the neighborhood, up, down, left and right as is shown as the black grids in \ref{pic:grid}.

\begin{figure}[ht!]
 	\centering
 	\includegraphics[scale=0.35]{grid.jpg}
 	\caption{\label{pic:grid}The picture shows that each pixel is often more related to the four pixels that surround it.}
\end{figure}

The formulation of the local diagonal pattern is the same as the dense pattern that we have discussed in Chapter \ref{Ch4},
except that this time the pattern of $L$ is restrained to only have nonzero entries in several diagonals.
The problem to solve for the local diagonal pattern is
\begin{eqnarray*}
	\min_{L} && {\rm AMSE} \\
	{\rm s.t.} && L \in {\rm local} \hspace{2mm} {\rm diagonal}.
\end{eqnarray*}

\section{The $\ell^{1}$ norm pattern}

A second approach is to have a sparse $L$ with only few nonzero entries so that the computational cost is much less. 
We consider adding an $\ell^{1}$ norm to the average MSE to control the sparsity of $L$. The number of nonzeros in 
$L$ is reduced by increasing the value of $\beta$. 

Thus, the optimization problem to solve for the $\ell^{1}$ norm pattern:
$$
\displaystyle \min_{L} {\rm AMSE} + \beta \parallel L \parallel_{1}.
$$
In order to find the best sparsity, we plot the AMSE as a function 
of the number of nonzeros in $L$, as is described in the $A$-optimal design work. The optimal sparsity is decided by 
the corner of the L-curve, which suggests an optimal compromise between the reconstruction quality and the sparsity of $L$.

\section{The Kronecker product pattern}

The last sparse pattern we consider is the Kronecker product pattern, for which we define
$$ 
L = L_{1} \otimes L_{2}, 
$$
where $L_{1}$ and $L_{2}$ are two smaller matrices. Hence, rather than minimizing with respect to a large dense matrix $L$ 
we now minimize the objective functional with respect to smaller matrices $L_{1}$ and $L_{2}$. Since all the derivatives 
and optimization process are now done on smaller matrices, the computational cost is greatly reduced.

In this case, the original covariance design turns into an optimization problem with respect to two small
matrices $L_{1}$ and $L_{2}$ as
\begin{equation}
\begin{array}{l}
\displaystyle \min_{L1, L2} \hspace{2mm} \alpha ^{2} \parallel (A^{\top}A + \alpha (L_{1} \otimes L_{2})^{\top}(L_{1} \otimes L_{2}))^{-1}(L_{1} \otimes L_{2})^{\top}(L_{1} \otimes L_{2})\mu \parallel^{2} + \\
{\rm trace} [ (A^{\top}A + \alpha (L_{1} \otimes L_{2})^{\top}(L_{1} \otimes L_{2}))^{-1} \\
(\alpha^{2}(L_{1} \otimes L_{2})^{\top}(L_{1} \otimes L_{2})\Sigma (L_{1} \otimes L_{2})^{\top}(L_{1} \otimes L_{2}) + \sigma^{2}A^{\top}A) \\
(A^{\top}A + \alpha (L_{1} \otimes L_{2})^{\top}(L_{1} \otimes L_{2}))^{-1} ].
\end{array}
\end{equation} 
Similarly, the training design is 
\begin{equation}
\begin{array}{l}
\displaystyle \min_{L1, L2} \hspace{2mm} {\rm trace} [ (A^{\top}A + \alpha (L_{1} \otimes L_{2})^{\top}(L_{1} \otimes L_{2}))^{-1} \\
(\alpha^{2}(L_{1} \otimes L_{2})^{\top}(L_{1} \otimes L_{2})M (L_{1} \otimes L_{2})^{\top}(L_{1} \otimes L_{2}) + \sigma^{2}A^{\top}A) \\
(A^{\top}A + \alpha (L_{1} \otimes L_{2})^{\top}(L_{1} \otimes L_{2}))^{-1} ].
\end{array}
\end{equation} 

\section{Numerical optimization of different sparse patterns}

As we have mentioned, in this section, totally we need to deal with three different patterns for the structure of $L$ and
the matrix derivatives techniques we use are different in each pattern. We will discuss them separately in this section.
Below is the formulation for the covariance design. The training design is very similar to the second term in the covariance
design, therefore we skip the details.

\subsection{The local diagonal pattern}

The computation of derivatives for the local diagonal pattern is very different from what we have seen in the dense pattern
in the previous chapter because the derivatives are taken entry-wise. Again we separate the AMSE into two terms $U$ and $V$.

We derive the derivative of the first term in details as an example. In order to distinguish the local diagonal pattern
from the dense pattern, we denote the regularization matrix as $L_{l}$.
\begin{eqnarray}
\frac{\partial U}{\partial L_{l}(i,j)} & = & \label{dUdL51} \frac{\partial [\mu^{\top}L_{lnf}^{\top}L_{l}(A^{\top}A + \alpha L_{l}^{\top}L_{l})^{-2}L_{l}^{\top}L_{l}\mu]}{\partial L_{l}(i,j)} \\
                              & + & \label{dUdL52} \frac{\partial [\mu^{\top}L_{l}^{\top}L_{lnf}(A^{\top}A + \alpha L_{l}^{\top}L_{l})^{-2}L_{l}^{\top}L_{l}\mu]}{\partial L_{l}(i,j)} \\
                              & + & \label{dUdL53} \frac{\partial [\mu^{\top}L_{l}^{\top}L_{l}(A^{\top}A + \alpha L_{lnf}^{\top}L_{l})^{-1}(A^{\top}A + \alpha L_{l}^{\top}L_{l})^{-1}L_{l}^{\top}L_{l}\mu]}{\partial L_{l}(i,j)} \\
                              & + & \label{dUdL54} \frac{\partial [\mu^{\top}L_{l}^{\top}L_{l}(A^{\top}A + \alpha L_{l}^{\top}L_{lnf})^{-1}(A^{\top}A + \alpha L_{l}^{\top}L_{l})^{-1}L_{l}^{\top}L_{l}\mu]}{\partial L_{l}(i,j)} \\
                              & + & \label{dUdL55} \frac{\partial [\mu^{\top}L_{l}^{\top}L_{l}(A^{\top}A + \alpha L_{l}^{\top}L_{l})^{-1}(A^{\top}A + \alpha L_{lnf}^{\top}L_{l})^{-1}L_{l}^{\top}L_{l}\mu]}{\partial L_{l}(i,j)} \\
                              & + & \label{dUdL56} \frac{\partial [\mu^{\top}L_{l}^{\top}L_{l}(A^{\top}A + \alpha L_{l}^{\top}L_{l})^{-1}(A^{\top}A + \alpha L_{l}^{\top}L_{lnf})^{-1}L_{l}^{\top}L_{l}\mu]}{\partial L_{l}(i,j)} \\
                              & + & \label{dUdL57} \frac{\partial [\mu^{\top}L_{l}^{\top}L_{l}(A^{\top}A + \alpha L_{l}^{\top}L_{l})^{-2}L_{lnf}^{\top}L_{l}\mu]}{\partial L_{l}(i,j)} \\
                              & + & \label{dUdL58} \frac{\partial [\mu^{\top}L_{l}^{\top}L_{l}(A^{\top}A + \alpha L_{l}^{\top}L_{l})^{-2}L_{l}^{\top}L_{lnf}\mu]}{\partial L_{l}(i,j)},
\end{eqnarray}
where we treat $L_{l}$ as fixed constant and $L_{lnf}$ as non-fixed variable. Let's discuss these equations one by one.
Denote $J_{i,j}$ to be the single-entry matrix with 1 at $(i,j)$ and 0 elsewhere. \\
For (\ref{dUdL51}) and (\ref{dUdL58}),
\begin{eqnarray*}
(\ref{dUdL51}) & = & \mu^{\top}\frac{\partial L_{lnf}^{\top}}{\partial L_{l}(i,j)}L_{l}(A^{\top}A + \alpha L_{l}^{\top}L_{l})^{-2}L_{l}^{\top}L_{l}\mu \\
			   & = & \mu^{\top}(\frac{\partial L_{lnf}}{\partial L_{l}(i,j)})^{\top}L_{l}(A^{\top}A + \alpha L_{l}^{\top}L_{l})^{-2}L_{l}^{\top}L_{l}\mu \\
			   & = & \mu^{\top}J_{i,j}^{\top}L_{l}(A^{\top}A + \alpha L_{l}^{\top}L_{l})^{-2}L_{l}^{\top}L_{l}\mu = (\ref{dUdL58}).
\end{eqnarray*}
For (\ref{dUdL52}) and (\ref{dUdL57}),
\begin{eqnarray*}
(\ref{dUdL52}) & = & \mu^{\top}L_{l}^{\top}\frac{\partial L_{lnf}}{\partial L_{l}(i,j)}(A^{\top}A + \alpha L_{l}^{\top}L_{l})^{-2}L_{l}^{\top}L_{l}\mu \\		   
			   & = & \mu^{\top}L_{l}^{\top}J_{i,j}(A^{\top}A + \alpha L_{l}^{\top}L_{l})^{-2}L_{l}^{\top}L_{l}\mu = (\ref{dUdL57}).
\end{eqnarray*}
For (\ref{dUdL53}) and (\ref{dUdL56}),
\begin{eqnarray*}
(\ref{dUdL53}) & = & \mu^{\top}L_{l}^{\top}L_{l}\frac{\partial (A^{\top}A + \alpha L_{lnf}^{\top}L_{l})^{-1}}{\partial L_{l}(i,j)}(A^{\top}A + \alpha L_{l}^{\top}L_{l})^{-1}L_{l}^{\top}L_{l}\mu \\
			   & = & - \mu^{\top}L_{l}^{\top}L_{l}(A^{\top}A + \alpha L_{l}^{\top}L_{l})^{-1}\frac{\partial (A^{\top}A + \alpha L_{lnf}^{\top}L_{l})}{\partial L_{l}(i,j)}(A^{\top}A + \alpha L_{l}^{\top}L_{l})^{-2}L_{l}^{\top}L_{l}\mu \\
			   & = & - \alpha \mu^{\top}L_{l}^{\top}L_{l}(A^{\top}A + \alpha L_{l}^{\top}L_{l})^{-1}J_{i,j}^{\top}L_{l}(A^{\top}A + \alpha L_{l}^{\top}L_{l})^{-2}L_{l}^{\top}L_{l}\mu = (\ref{dUdL56}).
\end{eqnarray*}
For (\ref{dUdL54}) and (\ref{dUdL55}),
\begin{eqnarray*}
(\ref{dUdL54}) & = & \mu^{\top}L_{l}^{\top}L_{l}\frac{\partial (A^{\top}A + \alpha L_{l}^{\top}L_{lnf})^{-1}}{\partial L_{l}(i,j)}(A^{\top}A + \alpha L_{l}^{\top}L_{l})^{-1}L_{l}^{\top}L_{l}\mu \\
			   & = & - \mu^{\top}L_{l}^{\top}L_{l}(A^{\top}A + \alpha L_{l}^{\top}L_{l})^{-1}\frac{\partial (A^{\top}A + \alpha L_{l}^{\top}L_{lnf})}{\partial L_{l}(i,j)}(A^{\top}A + \alpha L_{l}^{\top}L_{l})^{-2}L_{l}^{\top}L_{l}\mu \\
			   & = & - \alpha \mu^{\top}L_{l}^{\top}L_{l}(A^{\top}A + \alpha L_{l}^{\top}L_{l})^{-1}L_{l}^{\top}J_{i,j}(A^{\top}A + \alpha L_{l}^{\top}L_{l})^{-2}L_{l}^{\top}L_{l}\mu = (\ref{dUdL55}).
\end{eqnarray*}
Hence,
\begin{eqnarray*}
\frac{\partial U}{\partial L_{l}(i,j)} & = & 2 \mu^{\top}J_{i,j}^{\top}L_{l}(A^{\top}A + \alpha L_{l}^{\top}L_{l})^{-2}L_{l}^{\top}L_{l}\mu \\
									   & + & 2 \mu^{\top}L_{l}^{\top}J_{i,j}(A^{\top}A + \alpha L_{l}^{\top}L_{l})^{-2}L_{l}^{\top}L_{l}\mu \\
									   & - & 2 \alpha \mu^{\top}L_{l}^{\top}L_{l}(A^{\top}A + \alpha L_{l}^{\top}L_{l})^{-2}L_{l}^{\top}J_{i,j}(A^{\top}A + \alpha L_{l}^{\top}L_{l})^{-1}L_{l}^{\top}L_{l}\mu \\
									   & - & 2 \alpha \mu^{\top}L_{l}^{\top}L_{l}(A^{\top}A + \alpha L_{l}^{\top}L_{l})^{-1}L_{l}^{\top}J_{i,j}(A^{\top}A + \alpha L_{l}^{\top}L_{l})^{-2}L_{l}^{\top}L_{l}\mu .
\end{eqnarray*}
For $V$: \\
For each nonzero entry, the derivative is given by
$$
\grad_{L_{l}(i, j)} V = - 2 \alpha L_{l}^{i} (yz^{j} + zy^{j}) + 2 \alpha^2 L_{l}^{i}(\Sigma L_{l}^{\top}L_{l}zz^{j} + zz^{j}L_{l}^{\top}L_{l}\Sigma). 
$$
where $L_{l}^{i}$ denotes the $i^{th}$ low of $L$ and $z^{j}$ and $y^{j}$ denote the $j^{th}$ entry of the vectors $z$ and $y$.

Finally, we obtain
$$
\frac{\partial {\rm AMSE}}{\partial L_{l}(i,j)} = \alpha^2 \frac{\partial U}{\partial L_{l}(i,j)} + \frac{\partial V}{\partial L_{l}(i,j)}.
$$
Notice, in the derivative of ${\rm AMSE}$ with respect to $L_{l}$, we need to compute $L_{l}^{\top}J_{i,j}$, which is basically
just a single-column matrix generated by putting the transpose of the $i^{th}$ row of the matrix $L_{l}$ on the $j^{th}$ column. Also, we only need to do
the above derivatives on non-zero elements of $L_{l}$. Hence, the computational cost is greatly reduced.

In particular, for 1D problems, we only have nonzero entries on the tri-diagonals 
of the regularization matrix. Thus, the above computation is only done of those three diagonals. Since the cost of this 
derivative is less than $3n$ while it is $n^{2}$ in the case of the dense $L$, the local diagonal is considered to be 
one of the most efficient patterns. 

\subsection{The $\ell^{1}$ norm pattern}

There are two parts in the $\ell^{1}$ pattern: the AMSE and the $\ell^{1}$ norm. The only part that is new to us is the
$\ell^{1}$ norm of a matrix, which is simply the maximum absolute column sum of the matrix based on the theory of the 
induced matrix norm:
$$
\parallel L \parallel_{1} = \max_{1 \leq j \leq n} \sum_{i=1}^{m} |L_{i, j}|.
$$ 
However, it is very difficult to define a proper derivative for the above norm. Therefore, we treat the $\ell^{1}$ norm
as the sum of the absolute value of all entries in $L$.
$$
\parallel L \parallel_{1} = \sum_{i,j=1}^{m,n} |L_{i, j}|.
$$
The derivative of this norm is simply a matrix of the signs of the entries of $L$. In other words,
$$
\frac{\partial \parallel L \parallel_{1}}{\partial L_{i, j}} = \left\{ \begin{array}{rcl}
1 & \text{ if } & L_{i, j} > 0\\
-1 & \text{ if } & L_{i, j} < 0\\
0 & \text{ if } & L_{i, j} = 0 .\\
\end{array} 
\right.
$$

\subsection{The Kronecker product pattern}

The difficulty in deriving derivatives for this pattern is that the regularization matrix is a function of two other matrices.
Hence, the Chain rule for matrix needs to be applied \cite{Matrixcook}.

\subsubsection{The Chain rule}

Let $S = f(N)$, that is, the matrix $S$ is a function of another matrix $N$. The goal is to find the derivative of the
function $g(S)$ with respect to $N$:
$$
\frac{\partial g(S)}{\partial N} = \frac{\partial g(f(N))}{\partial N}.
$$
By applying the Chain rule, we obtain the entry-wise derivative
$$
\frac{\partial g(S)}{\partial N_{i, j}} = \sum_{k=1}^{n} \sum_{l=1}^{n} \frac{\partial g(S)}{\partial S_{k, l}} \frac{\partial S_{k, l}}{\partial N_{i, j}}.
$$
Using matrix notation, this can be written as:
$$
\frac{\partial g(S)}{\partial N_{i, j}} = {\rm trace} \left[ (\frac{\partial g(S)}{\partial S})^{\top} \frac{\partial S}{\partial N_{i, j}} \right].
$$

\subsubsection{Derivative of ${\rm AMSE}$ with respect to $L_{1}$}

Using the Chain Rule for matrix derivatives, we need to take derivative of ${\rm AMSE}$ with respect to each element of $L_{1}$, say $L_{1}(i,j)$
and then construct $\frac{\partial {\rm AMSE}}{\partial L_{1}} = \{\frac{\partial {\rm AMSE}}{\partial L_{1}(i,j)}\}_{i,j}^{n}$, where
$$
\frac{\partial {\rm AMSE}}{\partial L_{1}(i,j)} = {\rm trace} [ (\frac{\partial {\rm AMSE}}{\partial L})^{\top}\frac{\partial L}{\partial L_{1}(i,j)} ].
$$
From the formula above, we observe that both $\frac{\partial {\rm AMSE}}{\partial L}$ and $\frac{\partial L}{\partial L_{1}(i,j)}$ have the size of $n^2 \times n^2$.
Thus, the computational cost is pretty expensive when $n$ is large. In order to reduce the computational cost, we need to do some further observation.

As an example, for the derivative with respect to the element $L_{1}(i,j)$,
$$
\frac{\partial {\rm AMSE}}{\partial L_{1}(i,j)} = {\rm trace}(A_{j,i}L_{2})
$$
where $A_{j,i}$ is the $(j,i)^{th}$ $n \times n$ block of the matrix $A = \frac{\partial {\rm AMSE}}{\partial L}$.
Thus, the computational cost is reduced from $n^2 \times n^2$ to $n \times n$.

\subsubsection{Derivative of ${\rm AMSE}$ with respect to $L_{2}$}                      

Similarly, we take derivative of ${\rm AMSE}$ with respect to each element of $L_{2}$, say $L_{2}(i,j)$
and then construct $\frac{\partial {\rm AMSE}}{\partial L_{2}} = \{\frac{\partial {\rm AMSE}}{\partial L_{2}(i,j)}\}_{i,j}^{n}$, where
$$
\frac{\partial {\rm AMSE}}{\partial L_{2}(i,j)} = {\rm trace} [ (\frac{\partial {\rm AMSE}}{\partial L})^{\top}\frac{\partial L}{\partial L_{2}(i,j)} ].
$$
Again, we can reduce the computational cost after some observation.
As an example, for the derivative with respect to the element $L_{2}(i,j)$,
$$
\frac{\partial {\rm AMSE}}{\partial L_{2}(i,j)} = {\rm trace}(B_{j,i}L_{1})
$$
where $B_{j,i}$ is the matrix constructed by all the $(j,i)^{th}$ elements in each $n \times n$ block in $A$.
Again, the computational cost is reduced from $n^2 \times n^2$ to $n \times n$.

\section{Numerical experiments}

In this section, we experiment with different regularization matrices based on several sparsity patterns. 
Their performances are discussed and compared with other commonly used regularizations. The experiments will be
carried out on first a 1D toy example and then a real-data MRI example.

\subsection{The 1D magnetotelluric problem}

For the 1D example, we use the magnetotelluric problem as is described in chapter \ref{Ch2} again.  
Since the problem is discretized into $256$ points, the size of the regularization matrix is $256 \times 256$, which
will give us totally $65,536$ entries. Because this number is quite large, the goal of this problem is to find a sparse $L$ that 
still leads to a satisfactory recovery of the model $x$. 

For our numerical experiments, we provide both the covariance 
matrix and a set of training models so that the experiments can be done on both the covariance and training designs. 
In this experiment, we apply only the $\ell^{1}$ norm pattern to show the power of sparsity control.

\subsubsection{The covariance design}

We start by showing results from the covariance design. The mean of the model is set to be $0$. In order to check
how well our approaches perform, we design the covariance matrix such that the model will have some peaks.

\begin{figure}[ht!]
 	\centering
 	\includegraphics[scale=0.4]{CLcurve.jpg}
 	\caption{\label{covL1} The risk as a function of sparsity ${\rm nnz}(L)=\parallel L \parallel_{0}$ 
of the optimal $L$ obtained with the covariance design.}
\end{figure}  

Figure \ref{covL1} shows the L-curve, in which we plot the AMSE as a function of the sparsity of $L$.
As we can see, the AMSE decreases rapidly before the corner, which happens at $\parallel L \parallel_{0} = 16,467$, 
as opposed to the total $65,536$ entries in the dense $L$ pattern, and more slowly thereafter. 
Hence, this corner suggests an optimal compromise between the image reconstruction quality and the sparsity of $L$.  

\begin{figure}[ht!]
 	\centering
 	\includegraphics[scale=0.35]{goodtrue.jpg} \\
 	\includegraphics[scale=0.35]{goodmE.jpg}
 	\includegraphics[scale=0.35]{goodmD.jpg}
 	\caption{\label{ImagecovL1} The true image (top), reconstructed image using the optimal $16,467$ entries (bottom left) and
	 using the differential operator (bottom right).}
\end{figure} 

Figure \ref{ImagecovL1} shows the true model and we want to compare the reconstructed model from our approach to 
the result from the commonly used differential operator. Table \ref{table:covL1} displays the relative errors between 
the true model and the reconstructed ones. The advantage of our result is obvious.

\begin{table}[ht!]
\caption{\label{table:covL1}Relative errors of the reconstructed images: covariance design}
\centering
\begin{tabular}{cc}
\hline\hline
design & relative error \\
\hline
$\ell^{1}$    & $8.19\%$  \\
diff          & $17.02\%$  \\
\hline\hline
\end{tabular}
\end{table}

Also, we display the pattern of the sparse $L$ together with the pattern of the analytic $L$ in Figure \ref{fig:Lcov}.
We observe that the sparse $L$ has kept the parts that mainly contribute to the result and saved half of the effort.

\begin{figure}[ht!]
 	\centering
 	\includegraphics[scale=0.35]{LC.jpg}
 	\includegraphics[scale=0.35]{Ls.jpg}
 	\caption{\label{fig:Lcov} The left panel shows the pattern of the analytic $L$ and the right panel shows
the pattern of the sparse $L$ .}
\end{figure} 

\subsubsection{The training design}

In this section, we show results from the training design. In this case, the corner is located where the number 
of nonzeros of $L$ is $23,802$ as is shown in Figure \ref{trainL1}. 

\begin{figure}[ht!]
 	\centering
 	\includegraphics[scale=0.35]{ECLcurve.jpg}
 	\caption{\label{trainL1} The risk as a function of sparsity ${\rm nnz}(L)=\parallel L \parallel_{0}$ 
of the optimal $L$ obtained with the training design.}
\end{figure} 

\begin{figure}[ht!]
 	\centering
 	\includegraphics[scale=0.35]{goodtrue.jpg} \\
 	\includegraphics[scale=0.35]{ECmE.jpg}
 	\includegraphics[scale=0.35]{goodmD.jpg}
 	\caption{\label{ImagetrainL1} The true image (top), reconstructed image using the optimal $23,802$ entries (bottom left) and
	 using the differential operator (bottom right).}
\end{figure} 

\begin{table}[ht!]
\caption{\label{table:trainL1}Relative errors of the reconstructed images: training design}
\centering
\begin{tabular}{cc}
\hline\hline
design & relative error \\
\hline
$\ell^{1}$    & $9.60\%$  \\
diff          & $17.02\%$  \\
\hline\hline
\end{tabular}
\end{table}

From the resulting images in Figure \ref{ImagetrainL1} and Table \ref{table:trainL1}, we can see that again our 
result is much better than the one from the differential operator.

\subsection{An MRI example}

In the second experiment, we will apply our algorithms on some MRI examples. Since this is a real problem and we do not 
have the true covariance matrix, we will use the training design.

For the choice of the training models, we use the MRI data set in MATLAB which contains $27$ different image slices of 
a human head. The slices go from the nose gradually to the top of the head.

In each experiment, we choose four different sparse patterns of $L$: the dense pattern, the commonly used differential 
operator, the Kronecker product pattern and the 5-diagonal pattern since this is a 2D problem. Figure \ref{fig:LTL} 
show the matrix $L^{\top}L$ of the 5-diagonal pattern and the differential operator.

\begin{figure}[ht!]
 	\centering
 	\includegraphics[scale=0.35]{Ldiff.jpg}
 	\includegraphics[scale=0.35]{L5diag.jpg}
 	\caption{\label{fig:LTL} Images of $L^{\top}L$: the differential operator (left) and the 5-diagonal pattern (right).}
\end{figure} 

In the differential operator case, both the structure and the entry values are fixed. In the 5-diagonal case, 
we only keep the structure of the matrix while all these entry values are decided through optimization. For the choice 
of $\alpha$, we use the discrepancy principle to decide it before we start the iterations and set it as a constant 
during the optimization process.

\subsubsection{Experiment 1}

In the first experiment, we choose the $10^{th}$, $11^{th}$, $13^{th}$ and $14^{th}$ MRI examples to be our training models. 
Our goal is to use them to reconstruct the $12^{th}$ MRI, which is similar to the training set. Figure \ref{fig:true} show
the four training models we have used.

\begin{figure}[ht!]
 	\centering
 	\includegraphics[scale=0.35]{true10.jpg}
 	\includegraphics[scale=0.35]{true11.jpg} \\
 	\includegraphics[scale=0.35]{true13.jpg}
 	\includegraphics[scale=0.35]{true14.jpg}
 	\caption{\label{fig:true} Four training MRI models: $10^{th}$ (top-left), $11^{th}$ (top-right), $13^{th}$ (bottom-left)
and $14^{th}$ (bottom-right).}
\end{figure}

Here we display the true image and all the four reconstructed images from different patterns of $L$ in Figure \ref{fig:rec}. The top panel is the 
true image, the middle-left panel is the reconstructed image from the dense pattern, the middle-right panel is from the differential 
operator, the bottom-left panel is from the Kronecker product pattern and the bottom-right panel is from the 5-diagonal pattern. 

\begin{figure}[ht!]
 	\centering
 	\includegraphics[scale=0.35]{true12.jpg} \\
 	\includegraphics[scale=0.35]{dense12alpha0.jpg}
 	\includegraphics[scale=0.35]{diff12.jpg} \\
 	\includegraphics[scale=0.45]{kron12alpha0.jpg}
 	\includegraphics[scale=0.35]{diag12alpha0.jpg}
 	\caption{\label{fig:rec} The top panel is the 
true image, the middle-left panel is the reconstructed image from the dense pattern, the middle-right panel is from the differential 
operator, the bottom-left panel is from the Kronecker product pattern and the bottom-right panel is from the 5-diagonal pattern.}
\end{figure}

From the images, we observe that the dense pattern and the 5-diagonal pattern provide very nice reconstructions. 
The brain is very clear in these two images while it is just a big blur in the other two images.
Also, we observe that the Kronecker product pattern works a little bit better than the simple differential operator.

The relative errors between all the reconstructed images and the true image are listed in Table \ref{table:L}.

\begin{table}[ht!]
\caption{\label{table:L}Relative errors of the reconstructed images from four designs}
\centering
\begin{tabular}{cc}
\hline\hline
design & relative error \\
\hline
diff    & $25.27\%$  \\
Kron    & $21.65\%$  \\
5-diag  & $12.29\%$  \\
dense   & $9.96\%$  \\
\hline\hline
\end{tabular}
\end{table}

The explanation for the order of the reconstruction quality is straightforward: In the case of dense pattern, since no 
constraint is imposed on the regularization matrix $L$, the reconstruction should have the best possible quality for 
images that are similar to the training set. When it comes to the 5-diagonal pattern, based on the fact that
each pixel has stronger connection with the four pixels that surround it than other pixels, we only consider information 
in the neighborhood, therefore accuracy may be decreased a little bit because we completely ignore the influence from the rest 
of pixels in the image. The reconstruction will get better and better as we increase the number of diagonals.
For example, if we use nine diagonals, we should be able to get better results since now each pixel uses information from 
the eight pixels that surround it. However, although the reconstruction quality of the dense pattern is a bit better, 
the computational cost becomes tremendous. 

In the case of the Kronecker product pattern, each column in the image is treated
in the same way by the small matrix $L_{2}$ and each row is treated in the same way by $L_{1}$. For most image examples, 
this certainly will not work well. That is why the reconstruction from the Kronecker product pattern is much worse than 
the 5-diagonal pattern. The same thing happens in the case of differential operator. Moreover, it only considers the
one-sided difference while the Kronecker product pattern considers the centered difference. Thus, it is even worse than 
the Kronecker product pattern. Therefore, we have demonstrated that the order of reconstruction quality is as expected.

Therefore, we say our 5-diagonal pattern seems to be the optimal choice of the regularization matrix. 
It is very cheap in computation and the results are almost as good as the ones from the dense pattern.

One thing interesting to observe is the diagonal of the matrix $L^{\top}L$ from the optimal $L$ we
have obtained from the 5-diagonal pattern, as is plotted in Figure \ref{fig:diagLTL}.

\begin{figure}[ht!]
 	\centering
 	\includegraphics[scale=0.35]{5diagLTL.jpg}
 	\caption{\label{fig:diagLTL} The image of the diagonal of the matrix $L^{\top}L$.}
\end{figure}

We analyze the meaning of the figure as follows. Write the matrix $L$ as
$$
L = [ l_{1}, l_{2}, \cdots, l_{n}],
$$ 
where $l_{i}$ denotes a column of the matrix $L$. Hence,
$$
\diag(L^{\top}L) = [l_{1}^{\top}l_{1}, l_{2}^{\top}l_{2}, \cdots, l_{n}^{\top}l_{n}].
$$
On the other hand, we observe that
$$
\begin{array}{ccc}
\parallel Lx \parallel^{2} & = & L_{11}^{2}x_{1}^{2} + L_{12}^{2}x_{2}^{2} + \cdots + L_{1n}^{2}x_{n}^{2}\\
~ & ~ & \vdots\\
~ & ~ & + L_{n1}^{2}x_{1}^{2} + L_{n2}^{2}x_{2}^{2} + \cdots + L_{nn}^{2}x_{n}^{2} ,\\
\end{array}
$$
which gives us
$$
\begin{array}{ccl}
\parallel Lx \parallel^{2} & = & (L_{11}^{2} + \cdots + L_{n1}^{2})x_{1}^{2} + \cdots + (L_{1n}^{2} + \cdots + L_{nn}^{2})x_{n}^{2} \\
& = & l_{1}^{\top}l_{1}x_{1}^{2} + \cdots + l_{n}^{\top}l_{n}x_{n}^{2}.
\end{array}
$$
Hence, we have shown that image of $\diag (L^{\top}L)$ shows the amount of penalty that is added to each pixel of the image.

\subsubsection{Experiment 2}

We notice that, in the previous experiment, the training models we have chosen are very similar to the model we 
want to reconstruct. Unfortunately, for an unknown image, we do not know if it is similar to or different from 
the training models we have. What we could do is to gather lots of training models with relatively different 
information and see what will happen. `Relatively different' means they cannot be too different. For example, if we are trying 
to reconstruct the image of a patient's brain, then images containing different slices of the brain or views from 
different angles will be a good set of training models, while images of a tree certainly will not be helpful.

In this experiment, we randomly chose twenty images out of the twenty-seven MRI examples in Matlab to be our training 
models and the goal is to recover the remaining six images. In Table \ref{table:random} we list the relative errors 
between each reconstructed image from the simple differential operator, the 5-diagonal pattern and the true one.

\begin{table}[ht!]
\caption{\label{table:random}Relative errors of the reconstructed images from two designs}
\centering
\begin{tabular}{ccccccc}
\hline\hline
       & $1$       & $8$       & $9$       & $13$      & $19$      & $24$      \\
\hline
diff   & $34.57\%$ & $33.95\%$ & $29.43\%$ & $24.12\%$ & $25.28\%$ & $34.03\%$ \\
diag   & $32.42\%$ & $24.65\%$ & $23.38\%$ & $19.19\%$ & $20.33\%$ & $36.25\%$ \\
\hline\hline
\end{tabular}
\end{table}

By comparing the relative errors, we see that our approach works quite reasonably. The improvement differs in each 
image due to the different levels of similarity of that image to the training set. In real applications, we expect even 
better results because now we are taking different slices of the whole head while, in real cases, for example, a brain 
doctor will only take images of the brain part. Hence, the choice of training models will be more appropriate, which
will result in better reconstruction quality.