\documentclass[11pt]{article}
\usepackage{epsfig}
\usepackage{graphicx}
\usepackage{color}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{xypic}
\usepackage{nicematrix}
\usepackage{subcaption}
\usepackage{listings}
\lstset{language=Matlab}

\setlength{\textwidth}{6.5in}
\setlength{\topmargin}{0in}
\setlength{\headsep}{0in}
\setlength{\headheight}{0in}
\setlength{\oddsidemargin}{0in}
\setlength{\textheight}{9in}

\newcommand{\todo}[1]{\textcolor{red}{[#1]}}

\input xy
\xyoption{all}

\providecommand{\SolutionColumnOrientedBackwardSolve}{}
\providecommand{\SolutionBandedLUSolve}{}
\providecommand{\SolutionColumnDiagonallyDominant}{}
\providecommand{\SolutionRowDiagonallyDominantPartA}{}
\providecommand{\SolutionRowDiagonallyDominantPartB}{}
\providecommand{\SolutionRowDiagonallyDominantPartC}{}
\providecommand{\SolutionRowDiagonallyDominantPartD}{}
\providecommand{\SolutionMatrixReOrderingPartA}{}
\providecommand{\SolutionMatrixReOrderingPartB}{}
\providecommand{\SolutionImageDenoisingPartA}{}
\providecommand{\SolutionImageDenoisingPartB}{}
\providecommand{\SolutionImageDenoisingPartC}{}

%Used to display the marks available in the left margin.
\newcommand{\rmp}{\reversemarginpar\marginpar}

\pagestyle{empty}

\begin{document}

\begin{center}
{\Large {\bf CS 475/675 Spring 2025: Crowdmark Assignment 3}}\\
\end{center}

\noindent
{\bf Due June 16 at 11:59 pm Eastern.} \\
Submit all components of your solutions (written/analytical work, code/scripts, figures,
plots, output, etc.) to CrowdMark in PDF form in the section for each question.\\

\noindent
You must also separately submit a single zip file containing any and all code/scripts you
write to the A03 DropBox on LEARN, in runnable format (that is, .m).\\

\noindent
You might find these variants of the \textbf{triangle inequality} handy for Q3 and/or Q4:
\begin{eqnarray*}
|a|-|b| & \leq & |a+b| \leq |a|+|b|,\\
|a-b| & \geq & |a|-|b|.
\end{eqnarray*}

\noindent
\textbf{For full marks, be sure to show all your work!}

\begin{enumerate}

\item
\textbf{Column-oriented backward solve.} \rmp{[10]}
Consider solving a system $Ux=b$ where $U$ is upper triangular.
One approach to its solution is the following process.
First, resolve $x_{n}$.
It can then be removed from equations $1$ through $n-1$, and we can proceed with the reduced system in which the right-hand-side has been updated.
We next compute $x_{n-1}$ and remove it from equations $1$ through $n-2$, etc.
For example, if this approach is applied to
\begin{equation}
\begin{bmatrix} 3 & 1 & -1 \\ 0 & 7 & -2 \\ 0 & 0 & 2 \end{bmatrix}  \left[ \begin{array}{c} x_1 \\ x_2 \\ x_3 \end{array} \right] = \left[ \begin{array}{c} -2 \\ 17 \\ 4 \end{array} \right]
\end{equation}
we first compute $x_3 = 2$, and then solve the reduced $2 \times 2$ system:
\begin{equation}
\begin{bmatrix} 3 & 1 \\ 0 & 7 \end{bmatrix} \left[ \begin{array}{c} x_1 \\ x_2 \end{array} \right] = \left[ \begin{array}{c} -2 \\ 17 \end{array} \right] - 2\left[ \begin{array}{c} -1 \\ -2 \end{array} \right]
\end{equation}
Derive an efficient algorithm for the general case where the matrix $U$ is an $n \times n$ upper triangular matrix (not necessarily unit diagonal).
Present your algorithm in pseudocode.
The inputs to your algorithm should be the matrix $U$ and the right-hand side $b$, and the output should be the solution $x$.
What is the complexity of your algorithm?
Show your work, specifically the coefficients for all $O(n^{p})$ terms, $p \in \mathbb{N}$.
\SolutionColumnOrientedBackwardSolve

\newpage

\item
Consider \rmp{[10]}an $n \times n$ matrix, $A$, with upper and lower bandwidths equal to $q$.
Suppose you are given the $L$ and $U$ matrices comprising the LU decomposition of $A$, i.e., $A=LU$.
Derive an \textbf{efficient} algorithm (i.e., forward and backward solves) to solve $Ax=b$.
Give precise pseudocode for your algorithm.
What is the complexity of your algorithm, in terms of $n$ and $q$?
What is the coefficient of the leading order term (assuming $q \gg 1$)?
\SolutionBandedLUSolve

\newpage

\item
Let \rmp{[10]}$a_{ij}^{(k-1)}$ be the entries of $A$ after ($k-1$) steps of
Gaussian elimination.
Suppose $A^{(k-1)}$ is {\bf strictly column} diagonally dominant; i.e.
\begin{displaymath}
\left | a_{ii}^{(k-1)} \right | > \sum_{j \geq k, j \neq i} |a_{ji}^{(k-1)}|
\hspace{1cm} i=k,\ldots, n.
\end{displaymath}
If Gaussian elimination without pivoting is used; i.e.
\begin{displaymath}
a_{ji}^{(k)} = a_{ji}^{(k-1)} - \frac{a_{jk}^{(k-1)}a_{ki}^{(k-1)}}
{a_{kk}^{(k-1)}} \hspace{1cm} j=k+1, \ldots, n, \;\; i=k+1, \ldots, n,
\end{displaymath}
prove that
\begin{displaymath}
\left | a_{ii}^{(k)} \right | - \sum_{j \geq k+1, j \neq i} \left | a_{ji}^{(k)} \right | > 0
\hspace{1cm} i=k+1, \ldots, n.
\end{displaymath}
i.e. the submatrix $A^{(k)}$ is also column diagonally dominant.
(Hint: pay careful attention to the precise range of the index $j$.
Use the fact that:
$\sum_{j \geq k+1, j\neq i} \left | a_{ji}^{(k)} \right | = \sum_{j \geq k+1, j \neq i}
\left | a_{ji}^{(k-1)}-\frac{a_{jk}^{(k-1)}a_{ki}^{(k-1)}}{a_{kk}^{(k-1)}} \right |$
and then apply the triangle inequality.)
\SolutionColumnDiagonallyDominant

\newpage

\item
Let $A$ be a strictly {\bf row} diagonally dominant matrix; i.e. $A$'s entries $a_{ij}$ satisfy
\begin{displaymath}
\left | a_{ii} \right | > \sum_{j \neq i} \left | a_{ij} \right |.
\end{displaymath}
\begin{enumerate}
\item
As \rmp{[3]}we did in the lectures, let $A = D-L-U$, where $D$ is the diagonal, $-L$ contains only the entries below the diagonal of $A$, and $-U$ is the entries above the diagonal.
Let $G$ be the iteration matrix of the Gauss-Seidel method.

Suppose $x$ and $y$ are vectors such that $y = Gx.$
Show that $(D-L)y = Ux.$
\SolutionRowDiagonallyDominantPartA

\newpage

\item
Writing \rmp{[4]}out the expression in componentwise form yields
\begin{equation} 
       a_{jj} y_j + \sum_{i<j} a_{ji}y_i = -\sum_{i>j} a_{ji} x_i.
       \label{eqn:2b}
\end{equation} 
Recall that $\|y\|_{\infty} = \max_{j}|y_j|$, and consider the row $j$ corresponding to the largest magnitude entry.
Show that the absolute value of the left side of the equality~\eqref{eqn:2b} above satisfies
\begin{displaymath}
\left| a_{jj} y_j + \sum_{i<j} a_{ji}y_i \right| \geq \left( \left | a_{jj} \right | -\sum_{i<j} \left | a_{ji} \right | \right)\|y\|_{\infty}.
\end{displaymath}
\SolutionRowDiagonallyDominantPartB

\newpage

\item
Suppose \rmp{[4]}now that $x$ is a vector with infinity norm of 1.
Considering the right side of the
earlier equality~\eqref{eqn:2b}, show
\begin{displaymath}
\left|\sum_{i>j} a_{ji} x_i \right| \leq \left(\sum_{i>j} |a_{ji}|\right).
\end{displaymath}
\SolutionRowDiagonallyDominantPartC

\newpage

\item
Parts \rmp{[4]}(b)-(c) together imply that
\begin{displaymath}
\left( \left | a_{jj} \right | -\sum_{i<j} \left | a_{ji} \right | \right)\|y\|_{\infty} \leq \left ( \sum_{i>j} \left | a_{ji} \right | \right ).
\end{displaymath}
Recalling that an equivalent definition of the matrix infinity norm is
\begin{displaymath}
\|G\|_{\infty} \equiv \max_{\|x\|_{\infty}= 1}\|Gx\|_{\infty},
\end{displaymath}
show that Gauss-Seidel converges if the matrix $A$ is row-diagonally dominant, using the fact that $\rho(C) \leq \|C\|_{\infty}$ for any matrix $C$.
\SolutionRowDiagonallyDominantPartD

\end{enumerate}

\newpage

\item Consider a \textbf{symmetric} matrix whose graph is given by:
\begin{displaymath}
\entrymodifiers={++[o][F-]}
\SelectTips{cm}{}
\xymatrix @-1pc {
*\txt{} & 3 \ar@{-}[dr] & *\txt{} & *\txt{} & 6 \ar@{-}[r] \ar@{-}[ddr] & 7 \ar@{-}[dr] & *\txt{} \\
2 \ar@{-}[ur] \ar@{-}[dr] & *\txt{} & 1 \ar@{-}[r] & 5 \ar@{-}[ur] \ar@{-}[dr] \ar@{-}[urr] \ar@{-}[rrr] & *\txt{} & *\txt{} & 10 \\
*\txt{} & 4 \ar@{-}[ur] & *\txt{} & *\txt{} & 8 \ar@{-}[r] \ar@{-}[uur] & 9 \ar@{-}[ur] \\
}
\end{displaymath}
Ties will be broken by selecting the node with the numerically earlier original label (as given on the graph) first.
%In each case, redraw the final graph using the new labels.

\begin{enumerate}
\item
Perform \rmp{[5]}Cuthill-McKee (starting with node 1), and reverse Cuthill-McKee (starting with node 1) re-orderings on this graph.
\SolutionMatrixReOrderingPartA

\newpage

\item
Perform \rmp{[5]}minimum degree re-ordering on this graph.
\SolutionMatrixReOrderingPartB

\end{enumerate}

\newpage

\item
In PDE-based \textbf{image denoising}, using a simple Laplacian/diffusion approach produces overly-smoothed images that fail to preserve important edges in the data.
Approaches based on ``total variation regularization" tend to fare much better.
This can be expressed by solving the following equation,
repeatedly (for $k \geq 0$):
\begin{eqnarray}
\begin{array}{rcll}
-\alpha \nabla \cdot \frac{1}{|\nabla u^{k}|} \nabla
u^{k+1}(x,y) + u^{k+1}(x,y) & = & u^{0}(x,y) & \hspace{0.8cm}
\mbox{in $\; \Omega=(0,1) \times (0,1)$},
\end{array}
\label{eqn:image}
\end{eqnarray}
where we will use
\begin{eqnarray*}
\nabla u & = & \left ( \frac{\partial u}{\partial x}, \frac{\partial u}{\partial y} \right ) \textrm{, and}\\
|\nabla u| & = & \sqrt{ \left ( \frac{\partial u}{\partial x} \right )^2 + \left ( \frac{\partial u}{\partial y} \right )^2}.
\end{eqnarray*}
The parameter $\alpha > 0$ controls the degree of smoothing to be applied.
Figure 1 shows an example of an input image, the same image corrupted with noise, and a denoised image using the total variation approach.

\begin{figure}[htbp]
\centering
\includegraphics[height=3cm]{pic1.pdf}
\hspace{0.5cm}
\includegraphics[height=3cm]{pic2.pdf}
\hspace{0.5cm}
\includegraphics[height=3cm]{pic3.pdf}
\caption{(Left) original, (middle) noisy, and (right) denoised images.}
\end{figure}

One possible finite difference discretization of the PDE described above gives the following equation at each grid point
$(x_{i}, y_{j})$:
\begin{eqnarray}
AC u_{i,j}^{k+1} + AW u_{i-1,j}^{k+1} + AE u_{i+1,j}^{k+1} + AS u_{i,j-1}^{k+1} + AN u_{i,j+1}^{k+1} = u_{i,j}^0, \label{eqn:A}
\end{eqnarray}
where
\begin{eqnarray*}
AW &=& - \frac{\alpha}{h^2} \left (
       \frac{1}{2 \sqrt{ \left ( \frac{u_{i,j}^k-u_{i-1,j}^k}{h} \right )^2
       + \left ( \frac{u_{i,j}^k-u_{i,j-1}^k}{h} \right )^2 + \beta}} +
       \frac{1}{2\sqrt{ \left ( \frac{u_{i,j}^k-u_{i-1,j}^k}{h} \right )^2
       + \left ( \frac{u_{i-1,j+1}^k-u_{i-1,j}^k}{h} \right )^2 + \beta}} \right ) \\
AE &=& - \frac{\alpha}{h^2} \left(
       \frac{1}{2 \sqrt{ \left ( \frac{u_{i+1,j}^k-u_{i,j}^k}{h} \right )^2
       + \left ( \frac{u_{i+1,j}^k-u_{i+1,j-1}^k}{h} \right )^2 + \beta}} +
       \frac{1}{2 \sqrt{ \left ( \frac{u_{i+1,j}^k-u_{i,j}^k}{h} \right )^2
       + \left ( \frac{u_{i,j+1}^k-u_{i,j}^k}{h} \right )^2 + \beta}} \right ) \\
AS &=& - \frac{\alpha}{h^2} \left(
       \frac{1}{2 \sqrt{ \left ( \frac{u_{i,j}^k-u_{i-1,j}^k}{h} \right )^2
       + \left ( \frac{u_{i,j}^k-u_{i,j-1}^k}{h} \right )^2 + \beta}} +
       \frac{1}{2 \sqrt{ \left ( \frac{u_{i+1,j-1}^k-u_{i,j-1}^k}{h} \right )^2
       + \left ( \frac{u_{i,j}^k-u_{i,j-1}^k}{h} \right )^2 + \beta}} \right ) \\
AN &=& -\frac{\alpha}{h^2} \left(
       \frac{1}{2 \sqrt{ \left ( \frac{u_{i+1,j}^k-u_{i,j}^k}{h} \right )^2
       + \left ( \frac{u_{i,j+1}^k-u_{i,j}^k}{h} \right )^2 + \beta}} +
       \frac{1}{2 \sqrt{ \left ( \frac{u_{i,j+1}^k-u_{i-1,j+1}^k}{h} \right )^2
       + \left ( \frac{u_{i,j+1}^k-u_{i,j}^k}{h} \right )^2 + \beta}} \right) \\
AC &=& -(AW + AE + AS + AN) + 1.
\end{eqnarray*}
Note that the coefficients above depend on the data at step $k$ rather than $k+1$.
Here $h$ is grid cell size (width), and $\beta=10^{-6}$ is a constant parameter to circumvent numerical issues (division by zero).
The finite difference equations can be formulated into a linear system of the form
\begin{eqnarray}
A(u^{k}) u^{k+1} = u^{0},
\label{eqn:matrix}
\end{eqnarray}
where the coefficient matrix $A(u^{k})$ depends on the current values of $u^{k}$.
The matrix's sparsity pattern of nonzeros is the same as the standard 5-point stencil 2D Laplacian matrix.
You can assume zero boundary conditions.

\vspace{0.4cm}
\noindent
Our basic algorithm for image denoising process will be:

\begin{tabbing}
\hspace{3.5cm} \= Given a noisy image $u^{0}$. \\
\> for $k = 0, 1, \ldots, K$ \\
\> \hspace{0.5cm} Solve \hspace{0.1cm}
   $A(u^{k}) \; u^{k+1} = u^{0}$, \hspace{0.1cm} for $u^{k+1}$ \\
\> end
\end{tabbing}
where $K$ is some specified number of iterations of denoising to apply.

\noindent
To actually solve equation (\ref{eqn:matrix}) inside the loop, we will apply several of the iterative methods we have seen in class.
Let $u^{k+1,\ell}$ be the approximate solution of $u^{k+1}$ given by $\ell$ iterations of an iterative method.
Then the denoising algorithm can be written as:

\begin{tabbing}
\hspace{3.5cm} \= Given noisy image $u^{0}$. \\
\> for $k=0, 1, \ldots, K$ \\
\> \hspace{0.5cm}  \= $u^{k+1,0} = u^{k}$ \\
\> \> for $\ell=0, 1, \ldots$, until convergence \\
\> \> \hspace{0.5cm}
$
u^{k+1,\ell+1} = \mbox{IterativeMethodStep}(u^{k+1,\ell},A(u^{k}),u^{0})
$ \\
\> \> end \\
\> end
\end{tabbing}

\noindent
For this question, we will use $K=8$.
For $\alpha$, we will use $\alpha=6.4 \times 10^{-2}, 3.2 \times 10^{-2}$, $1.6 \times 10^{-2}$, and $8\times 10^{-3}$ for image sizes of $16 \times 16, 32 \times 32$, $64 \times 64$, and $128 \times 128$, respectively.

\begin{enumerate}
\item
Create \rmp{[10]}two MATLAB functions:
\begin{lstlisting}
A = FormMatrix(u, alpha)
u0 = FormRHS(X)
\end{lstlisting}
The input for \lstinline!FormMatrix! is an approximate solution vector $u$ (corresponding to the $u^{k}$ terms in the coefficients for equation (\ref{eqn:image})) and the
output is the matrix corresponding to equation (\ref{eqn:image}).
Be sure to make use of Matlab's sparse matrix functionality!
The input for \lstinline!FormRHS! is a noisy image \lstinline!X! and the output is the vector representation \lstinline!u0! of \lstinline!X!.
Note that if the image size of \lstinline!X! is $m \times m$, then the size of \lstinline!A! should be $n \times n$ where $n=m^{2}$ and the size of \lstinline!u0! is $n \times 1$. 
Noisy images of different sizes can be generated by the provided Matlab file \lstinline!set_image!, which takes an integer specifying the resolution along one dimension.
(You may assume the image is always square.)
\SolutionImageDenoisingPartA

\newpage

\item
Implement \rmp{[12]}the following iterative methods: Jacobi, Gauss-Seidel, SOR, and Conjugate
Gradient, by creating the following MATLAB functions:
\begin{lstlisting}
       [x,iter] = Jacobi(A,b,x_initial,maxiter,tol)
       [x,iter] = GS(A,b,x_initial,maxiter,tol)
       [x,iter] = SOR(omega,A,b,x_initial,maxiter,tol)
       [x,iter] = CG(A,b,x_initial,maxiter,tol)
\end{lstlisting}
These functions take as inputs the sparse matrix \lstinline!A!, the right-hand side
\lstinline!b!, the initial guess \lstinline!x_initial!, the maximum number of iterations
\lstinline!maxiter!, and the tolerance \lstinline!tol!, and computes the approximate solution
\lstinline!x! using \lstinline!iter! steps of the corresponding iterative method.
(SOR has one more input parameter, \lstinline!omega!.)
You may use MATLAB's built-in backslash operator ('\textbackslash') to perform the necessary triangular solves in the inner loops of Gauss-Seidel and SOR.
For Jacobi, you may explicitly construct and use $D^{-1}$ (i.e., breaking our usual rule of not inverting matrices). 

The solvers should be stopped either when they reach the maximum number of solver iterations or when the residual vector satisfies:
\begin{displaymath}
\|r\|_2 \leq tol \, \|b\|_2.
\end{displaymath}
\SolutionImageDenoisingPartB

\newpage

\item
Create \rmp{[18]}a MATLAB script, \lstinline!Denoise!, that solves the system (\ref{eqn:A}) to perform denoising, using different iterative methods.
In particular, for each iteration $k$, set up the matrix $A$ and right-hand side $b=u0$ by calling \lstinline!FormMatrix! and \lstinline!FormRHS!.
Then solve the linear system by calling one of the iterative methods in part (b).
(Tip: for easier debugging, you may want to temporarily use MATLAB's backslash operator, \lstinline!A\b!, for the solve first to make sure your outer loop and the two functions, \lstinline!FormMatrix! and \lstinline!FormRHS!, are correct, i.e., get the overall denoising working before implementing your own linear solvers.
Then, replace \lstinline!A\b! with the various iterative methods.)

In this assignment, use a relatively large solver tolerance, $tol=10^{-2}$.
Impose a maximum of 20000 solver iterations (note: this is separate from the outer/denoising iterations).
For SOR, determine the optimal $omega$ (requiring the fewest \emph{total} iterations for the whole denoising process) for each grid size up to
2 decimal places by trial-and-error (by hand or by code) for the three image sizes specified below.
Report the optimal values you found.

Record the CPU times and number of iterations.
Construct a table of execution
times and a table of \emph{total} number of iterations for the denoising process when using Jacobi,
Gauss-Seidel, SOR (with your optimal omega values), and Conjugate Gradient iterations.
Use image sizes $16 \times
16$, $32 \times 32$, and $64 \times 64$ (and \emph{optionally}, $128 \times 128$,
if you would like and time permits).
Comment on the results.
Convert the output vector back to a matrix (2D image), and use \lstinline!imagesc! with a grayscale color map to
display the denoised image. 
\SolutionImageDenoisingPartC

\end{enumerate}

Submit: 
\begin{itemize}
\item Your code for \lstinline!FormMatrix, FormRHS, Jacobi, GS, SOR, CG! and \lstinline!Denoise!;
\item Your list of optimal \lstinline!omega! values.
\item The plots of the denoised images for the grid sizes specified above;
\item Your table of CPU times and iteration counts;
\item Comments on your observations.
\end{itemize}

\emph{Once again, submit a copy of everything to CrowdMark in PDF/image format, and separately submit your [runnable] source code to the A03 LEARN DropBox in a single ZIP file.}

\end{enumerate}

\end{document}