\documentclass[unicode,11pt,a4paper,oneside,numbers=endperiod,openany]{scrartcl}

\usepackage{graphicx}
\usepackage{subcaption}
\usepackage{amsmath}
\input{assignment.sty}
\usepackage{pgfplots}
\pgfplotsset{compat=newest}
\usetikzlibrary{plotmarks}
\usetikzlibrary{arrows.meta}
\usepgfplotslibrary{patchplots}
\usepackage{grffile}
\usepackage{amsmath}
\usepackage{subcaption}
\usepgfplotslibrary{external}
\tikzexternalize

\begin{document}

\setassignment
\setduedate{Wednesday, December 02, 2020, 11:59 PM}

\serieheader{Numerical Computing}{2020}{Student: Claudio Maggioni}{Discussed with: Gianmarco De Vita (3.3)}{Solution for Project 5}{}
\newline

\assignmentpolicy

The purpose of this assignment is to gain insight on the theoretical and
numerical properties of the Conjugate Gradient method. Here we use this method
in an image processing application with the goal of deblur an image given the
exact (noise-free) blurred image and the original transformation matrix. Note
that the ``noise-free" simplification is essential for us to solve this problem
in the scope of this assignment.

\section{General Questions [10 points]}
\subsection{What is the size of the matrix $A$?}
$A$ is an $n^2$ by $n^2$ matrix, where $n$ is the width and height in pixels
of the image to transform.

\subsection{How many diagonals bands does $A$ have?}
$A$ as $d^2$ diagonal bands, where $d$ is strictly an order of magnitude below
$n$.

\subsection{What is the length of the vectorized blur image $b$?}
$b$ is the row-vectorized form of the image pixel matrix, and thus has
dimensions 1 by $n^2$.

\section{Properties of A [10 points]}
\subsection{ If $A$ is not symmetric, how would this affect $\tilde{A}$?}
If $A$ were not symmetric, then $\tilde{A}$ would not be positive definite since
by definition $\tilde{A} = A A^T$, thus not satifying the assumptions taken when
solving the system.

\subsection{ Explain why solving $Ax = b$ for $x$ is equivalent to
minimizing $x^T A x - b^T x$ over $x$.}

First, we can say that:

\[f(x) = \frac12 x^T A x - b^T x = \frac12\langle Ax,x \rangle -
\langle b,x \rangle\]

Then by taking the derivative of $f(x)$ w.r.t. $x$ we have (assuming $A$ is
\textit{spd}, which it is):

\[f'(x) = \frac12 A^T x + \frac12Ax - b = Ax - b\]

which for $f'(x) = 0$ will be equivalent to solving $Ax = b$. By taking the
second derivative we have:

\[f''(x) = A > 0\]

since $A$ is positive definite. Therefore, we can say that the absolute minima
of $f(x)$ is the solution for $Ax = b$.

\section{Conjugate Gradient [40 points]}
\subsection{ Write a function for the conjugate gradient solver
\texttt{[x,rvec]=myCG(A,b,x0,max\_itr,tol)}, where \texttt{x}
and \texttt{rvec} are, respectively, the solution value and a
vector containing the residual at every iteration.}

The implementation can be found in file \texttt{myCG.m} in the source directory.
The test code for the function \texttt{myCG} can be found in the \texttt{test.m} file.

\subsection{ In order to validate your implementation, solve the system
defined by \texttt{A\_test.mat} and \texttt{b\_test.mat}. Plot
the convergence (residual vs iteration).}

The plot of the squared residual 2-norms over all iterations can be found in Figure
\ref{fig:plot1}.

\begin{figure}[h]
	\centering
	\resizebox{0.6\textwidth}{!}{\input{test_semilogy}}
	\caption{Semilog plot of the plot of the squared residual 2-norms over all iterations}\label{fig:plot1}
\end{figure}

\subsection{ Plot the eigenvalues of \texttt{A\_test.mat} and comment on the
condition number and convergence rate.}

The eigenvalues of A can be found in figure
\ref{fig:plot2}.

\begin{figure}[h]
	\centering
	\resizebox{0.6\textwidth}{!}{\input{A_eig}}
	\caption{Semilog plot of the eigenvalues of A}\label{fig:plot2}
\end{figure}

The condition number for matrix $A$ according to \texttt{cond(...)} is $\approx 1.3700 \cdot 10^6$,
which is rather ill conditioned. The eigenvalue plot agrees with the condition number, by showing that there are significant differences in the magnitude of the eigenvalues of \texttt{A\_test}.

The two methods agree due to the fact that the condition number is computed by using singular values,
which in turn are derived by eigenvalues. This fact was demonstrated practically by Dr. Edoardo Vecchi, \textit{future PhD}
using this MATLAB snippet.

\begin{verbatim}
X = A_test'*A_test;
singular_values = sqrt(eig(X));
norm(sort(singular_values,'descend') - svd(A_test))/numel(A_test)
\end{verbatim}

The final norm is less than $10^{12}$, showning that the definition of
\texttt{singular\_values}, which is a generalization of the computation of
eigenvalues for rectangular matrices, is equivalent to MATLAB's \texttt{svd(...)} other than approximation errors.


\section{Debluring problem [40 points]}
\subsection{ Solve the debluring problem for the blurred image matrix
\texttt{B.mat} and transformation matrix \texttt{A.mat} using
your routine \texttt{myCG} and Matlab's preconditioned conjugate
gradient \texttt{pcg}. As a preconditioner, use \texttt{ichol}
to get the incomplete Cholesky factors and set routine type to
\texttt{nofill} with $\alpha=0.01$ for the diagonal shift (see
Matlab documentation). Solve the system with both solvers using
$max\_iter=200$ $tol= 10^{-6}$. Plot the convergence (residual
vs iteration) of each solver and display the original and final
deblurred image.}

My implementation is in the file \texttt{deblurring.m}. Plots for the original image, the deblurred image from
\texttt{myCG}, the deblurred image from \texttt{pcg}, and a semi-logarithmic plot on the y-axis of the residuals
from the two conjugate gradient functions over the iteration count can be found respectively in figure \ref{fig:orig},
\ref{fig:mycg}, \ref{fig:pcg}, and \ref{fig:rvec}.

As a terminating condition for the \texttt{myCG} implementation, I chose to check if the residual divided by the 2-norm of
$b$ is less than the given tolerance. This mimicks the behaviour of MATLAB's \texttt{pcg} to allow for a better comparison.

\begin{figure}[h]
	\begin{subfigure}{0.33\textwidth}
	\centering
		\resizebox{0.8\textwidth}{!}{\input{img_orig}}
	\caption{Original image grayscale matrix}\label{fig:orig}
	\end{subfigure}
	\begin{subfigure}{0.33\textwidth}
	\centering
		\resizebox{0.8\textwidth}{!}{\input{img_my}}
	\caption{Deblurred image using \texttt{myCG}}\label{fig:mycg}
	\end{subfigure}
	\begin{subfigure}{0.33\textwidth}
	\centering
		\resizebox{0.8\textwidth}{!}{\input{img_rcg}}
	\caption{Deblurred image using \texttt{rcg}}\label{fig:pcg}
	\end{subfigure}
	\caption{Blurred and deblurred images}
\end{figure}
\begin{figure}[h]
	\centering
	\resizebox{0.6\textwidth}{!}{\input{res_log}}
	\caption{Residuals of \texttt{myCG} (in blue) and \texttt{rcg} (in orange) over iteration count (y axis is a log scale)}\label{fig:rvec}
\end{figure}
\subsection{ When would \texttt{pcg} be worth the added computational cost?
What about if you are debluring lots of images with the same
blur operator?}

The \textit{pcg} algorithm provided by MATLAB would be worth for many deblurring operations useing the same blur operator, since the cost of computing the incomplete cholesky decomposition (i.e. \texttt{ichol}) can be payed only once and thus amortized. \texttt{myCG} is better for few iterations thanks to not needing any seed that is expensive to compute.


\end{document}