<html><head><meta name="color-scheme" content="light dark"></head><body><pre style="word-wrap: break-word; white-space: pre-wrap;">\documentclass[11pt,onecolumn]{article}
\usepackage{amssymb, amsmath, amsthm,graphicx, paralist,algpseudocode,algorithm,cancel,url,color}
\usepackage{sectsty}
\usepackage{fancyvrb}
\usepackage{mathrsfs}
\usepackage{multirow}
\usepackage{hhline}
\usepackage{booktabs}
\usepackage[table]{xcolor}
\usepackage{tikz}
% \usepackage[framed,numbered,autolinebreaks,useliterate]{mcode}
\usepackage{listings}
\usepackage{enumitem}
\usepackage{cleveref}
\usepackage{tcolorbox}


\newcommand{\bvec}[1]{\mathbf{#1}}
\newcommand{\R}{\mathbb{R}}
\newcommand{\C}{\mathcal{C}}
\newcommand{\Rn}{\R^{n\times n}}
\newcommand{\Rmn}{\R^{m\times n}}
\newcommand{\Cn}{\C^{n\times n}}
\newcommand{\Cmn}{\C^{m\times n}}
\newcommand{\cO}{\mathcal{O}}
\newcommand{\ls}{\textsc{ls}}
\DeclareMathOperator{\Tr}{Tr}
\DeclareMathOperator{\trace}{trace}
\DeclareMathOperator{\diag}{diag}
\DeclareMathOperator*{\argmin}{arg\,min}
\DeclareMathOperator*{\argmax}{arg\,max}
\sectionfont{\Large\sc}
\subsectionfont{\sc}
\usepackage[margin=1 in]{geometry}

\newcommand{\bluebox}[1]{
  \begin{tcolorbox}[colback=blue!5!white,colframe=blue!75!black,boxrule=0.5pt,boxsep=0pt,left=6pt,right=16pt,top=4pt,bottom=4pt]
  #1
  \end{tcolorbox}   
}

\newcommand{\redbox}[1]{
  \begin{tcolorbox}[colback=red!5!white,colframe=red!75!black,boxrule=0.5pt,boxsep=0pt,left=6pt,right=16pt,top=4pt,bottom=4pt]
  #1
  \end{tcolorbox}
} 

\newcommand{\greenbox}[1]{
  \begin{tcolorbox}[colback=green!5!white,colframe=red!75!black,boxrule=0.5pt,boxsep=0pt,left=6pt,right=16pt,top=4pt,bottom=4pt]
  #1
  \end{tcolorbox}
} 

\begin{document}
\noindent
\textsc{\Large Numerical analysis: Project 1}\\
Instructor: Anil Damle\\
Due: March 10, 2025\\

\subsection*{Policies}
You are permitted to work in groups of up to three for this project and each group only needs to submit a single report. While it is okay to discuss the project with other groups (similar to the HW collaboration policy), every group must produce their own code and report. The report should address all of the items referred to as \textbf{goals} or $\textbf{questions}$ in the project description. You may refer to inanimate sources (e.g., textbooks) to help you answer the questions and tackle the goals. However, your sources must be explicitly cited. For example, one of the goals requires writing pseudocode for an algorithm that exists, if you go to an external source to read more about it and help generate the pseudocoade you must cite where it came from.\\ 


\noindent
Part of the purpose of this project is to provide you with an opportunity to practice writing more free form reports and carefully choosing what plots, results, etc. are needed to convince us your solution is correct. Therefore, some of the goals are leading, but do not give a concrete list of exactly what has to be included. The questions tend to be more concrete.\\



\noindent
The report, including plots and requested output from your code,
should be typeset and submitted via the Gradescope as a pdf file. This file must be self contained for grading. Additionally, please submit any code written for the assignment as zip file to the separate Gradescope assignment for code.\\
\hrule
\vspace{1cm}
\noindent


\subsection*{Preamble}
For this project, we are going to work with images, and a simple form of matrix completion. We have picked a specific application and algorithm that you have all of the technical tools to complete. The state of the art in this field involves more complicated algorithms and, for images in particular, there are modifications that can be made to significantly improve the results. Furthermore, I have made some slight modifications to the underlying images to ensure things work out reasonablly here. If you are curious about this area I would be happy to chat more about it.

For the purposes of this project, all of the images we will work with are gray scale and represented by matrices with real valued entries between 0 and 1 (0 for black and 1 for white). So, an image that is $n_1 \times n_2$ pixels will be represented by an $n_1 \times n_2$ matrix. These techniques can be expanded to deal with color images, but here we will focus on the linear algebra underpinnings. When producing plots/images you should set the color space to gray so that it looks as you would expect for a gray-scale image. 

\subsection*{Compressing images}
Given a matrix $A \in \R^{n_1 \times n_2}$ we may want to store a compressed representation of it. One way to do this is to store a low-rank representation of $A,$ and this can be accomplished using the SVD. Say we compute the SVD
\[
A = USV^T,
\]
where $U\in \R^{n_1 \times p},$ $V\in \R^{n_2 \times p},$ $S\in \R^{p \times p},$ and $p = \min\left\{n_1,n_2\right\}.$ Let $U_k$ denote the first $k$ columns of $U,$ $V_k$ denote the first $k$ columns of $V,$ and $S_k$ denote the upper left $k \times k$ block of $S.$ Now, we can write $A\approx U_kS_kV_k^T$ and, in fact, this is the optimal approximation of $A$ in either the 2 or Frobenius norm. This means that of all rank-$k$ matrices (those with $k$ non-zero singular values), the one closest to $A$ in 2 or Frobenius norm is $U_kS_kV_k^T.$ The error of approximation is determined by the singular values of $A,$ specifically
\[
\left\|A-U_kS_kV_k^T\right\|_2 = \sigma_{k+1} \qquad \text{and} \qquad \left\|A-U_kS_kV_k^T\right\|_F = \left(\sum_{i=k+1}^p\sigma_{i}^2\right)^{1/2}.
\] 
So, if the singular values of $A$ decay rapidly there may be a good low-rank approximation for $A,$ and if $k$ is small enough, storing $U_k, V_k,$ and $S_k$ may be cheaper than storing $A.$

Now, while we have not talked about them, there are good algorithms to compute the SVD. However, in certain settings it may not be possible to use these algorithms directly. So, let's consider an alternative. We can write any rank $k$ approximation of $A$ as $WZ^T$ where $W\in \R^{n_1\times k}$ and $Z\in \R^{n_2\times k}.$ Our low-rank approximation problem  can then be cast as the optimization problem
\[
\min_{W,Z}\left\|A - WZ^T\right\|^2_F.
\] 
Unfortunately, robustly solving this problem can be a bit tricky if we do not resort to the SVD. However, one possible approach is known as alternating least squares. It is based on the observation that if we fix $W$ or $Z$ then we can solve for the other factor via a (matrix) least-squares problem. This gives Algorithm~\ref{alg:als}, though for our purposes we will run it for a small, fixed number of iterations rather than checking for convergence (which may be preferable in practice).\\
\begin{algorithm}
\caption{Alternating least squares, without regularization}
\label{alg:als}
\begin{algorithmic}[1]
\State {\bf initialize} $W \in \R^{n_1 \times k}$, $Z \in \R^{n_2 \times k}$
\While{not converged}
  \State $$Z \leftarrow \argmin_{Z}\left\|A - WZ^T\right\|^2_F$$ 
  \State $$W \leftarrow \argmin_{W}\left\|A - WZ^T\right\|^2_F$$ 
\EndWhile
\State {\bf return} $W$, $Z$
\end{algorithmic}
\end{algorithm}

\bluebox{
\paragraph{Goal 1:}
Implement the Householder scheme for computing QR factorizations. Your QR factorization routine should take in a $n \times k$ matrix $B$ with $k\leq n$ and returns a $n\times k$ matrix $Q$ with orthonormal columns (or sufficient information to be able to apply it to a vector, e.g., the Householder vectors), and a $k\times k$ upper triangular matrix $R.$ Please demonstrate that your algorithm behaves as expected (this means both that you are getting a valid QR factorization as output and that it achieves the desired computational scaling). \textbf{Demonstrate that your implementation achieves the expected scaling and explain how you tested your implementation.}}

\bluebox{
\paragraph{Goal 2:} 
Using your $QR$ factorization, implement Algorithm~\ref{alg:als}. This is written as a matrix least squares problem, but think about the Frobenius norm and how you might be able to split it up into problems we know how to solve. \textbf{Please explain how you do this in the project report.} Practically, it is possible to then block some of these operations together for efficiency; you may take advantage of this if you wish.}

\subsection*{A note on file formats.}
All matrices $A$ provided for this project are saved as .csv files where each line (delineated by a newline character) of the file denotes a row of $A$ and entries within a line are separated by commas. This should be easy enough to read in whatever language you are using.\\ 

\bluebox{
\paragraph{Goal 3:}
Using your implementation, load the file Cornell.csv which has an image stored in the matrix $C$ and try to compute the best rank 75 approximation of the image. Compare this with the result of the best rank 75 image from the SVD (you can use a built in routine to compute this), what do you observe qualitatively and quantitatively?}


\subsection*{What about if we do not know all the entries of our image?}
In the preceding section, alternating least squares seems like a roundabout way to get at something we know how to compute. But now we will explore a setting where we cannot compute the SVD as we might like. Specifically, we will consider being given an image $A$ where we only know the true values of some subset of the pixels, the rest are unobserved. In this case we can still use a variant of alternating least squares to try and recover an underlying low-rank approximation that should roughly resemble our image. 

The two key differences are that we will add regularization and only measure error on the observed pixels. Let $\Omega$ be a set of $(i,j)$ pairs that represents the set of pixels we observe in the image. So, for example, if $\Omega = \left\{(1,2),(10,20)\right\}$ then we are only observing the $(1,2)$ and $(10,20)$ pixels (matrix entries) of the image. We now define a restricted variant of the Frobenius norm as follows
\[
\left\|A-B\right\|_{F,\Omega}^2 = \sum_{(i,j)\in \Omega} \left(A_{i,j}-B_{i,j}\right)^2.
\]
This norm simply measures the difference between $A$ and $B$ on a subset of their entries defined by $\Omega.$ The second adjustment we will make is to penalize $W$ or $Z$ from becoming large using a parameter $\beta.$ This yields the following optimization problem
\[
\min_{W,Z}\left\|A - WZ^T\right\|^2_{F,\Omega} + \beta^2\left\|W\right\|^2_F + \beta^2\left\|Z\right\|^2_F. 
\]
The hope is that $WZ^T$ is then a good approximation of the underlying image.

Algorithm~\ref{alg:als2} is a version of alternating least squares with regularization and for incomplete information. Once again, you may simply run it for a fixed number of iterations (that you can tune based on the output) rather than checking for convergence as one would in practice.
\begin{algorithm}
\caption{Alternating least squares, with regularization and unknown entries}
\label{alg:als2}
\begin{algorithmic}[1]
\State {\bf initialize} $W \in \R^{n_1 \times k}$, $Z \in \R^{n_2 \times k}$
\While{not converged}
  \State $$Z \leftarrow \argmin_{Z}\left\|A - WZ^T\right\|^2_{F,\Omega} + \beta^2\left\|Z\right\|_F^2$$ 
  \State $$W \leftarrow \argmin_{W}\left\|A - WZ^T\right\|^2_{F,\Omega} + \beta^2\left\|W\right\|_F^2$$
\EndWhile
\State {\bf return} $W$, $Z$
\end{algorithmic}
\end{algorithm}

Turns out, you have all of the knowledge and machinery to implement this algorithm (save for one small item we develop in the first point below). 

\bluebox{
\paragraph{Question 1:}
Show that if we want to solve $$\min_x \|Ax-b\|_2^2 + \beta^2 \|x\|_2^2$$ we can instead solve 
\[
\min_x \left\|\begin{bmatrix}A \\ \beta I\end{bmatrix}x-\begin{bmatrix}b \\ 0\end{bmatrix}\right\|_2^2.
\]}

\bluebox{
\paragraph{Goal 4:}
Using your $QR$ factorization, implement Algorithm~\ref{alg:als2}. From the first set of goals, you should have worked out how to split this up into solving many least squares problems. Now, you have to think about what size those problems are, and which entries of the matrices they involve. Given that, you can then leverage your $QR$ factorization and the preceding item to implement the algorithm. \textbf{Please explain how you do this in the project report.}}

\bluebox{
\paragraph{Goal 5:}
Using your implementation, load each of the image files in Image*.csv and associated masks Mask*.csv, where * is 1, 2, and 3. Each file contains an image and the associated mask file contains an image of the same size whose entries are 1 if the corresponding entry of $C$ is observed and 0 if it is unobserved (so it defines the set $\Omega$). The unknown entries of $C$ have been set arbitrarily and you should not access or use them (as the will adversely affect your results). Using Algorithm~\ref{alg:als2}, try and recover each of the underlying images. You now have some parameters to consider and we encourage exploration of their values. You should not have to go to $k$ larger than 75 for any of the examples, and good $\beta$ to try are between $10^{-2}$ and 1. Report the images you recover, and discuss what you observe.}


\end{document}
</pre></body></html>