\documentclass{article}[12pt]

% useful packages
\usepackage{titlesec}
\usepackage{fullpage}
\usepackage{amsmath,amssymb,amsthm,amsfonts}
\usepackage{graphicx}
\usepackage{enumerate}
\usepackage{algorithm,algorithmic}
\usepackage{xcolor}
\usepackage{bbm}
\usepackage{url}
\usepackage{hyperref}
\usepackage{caption,subcaption}

% theorem type environments
\newtheorem{thm}{Theorem}
\newtheorem{prop}{Proposition}
\newtheorem{lemma}{Lemma}
\newtheorem{cor}{Corollary}
\newtheorem{defn}{Definition}
\newtheorem{assump}{Assumption}
\newtheorem{example}{Example}
\newtheorem{conjecture}{Conjecture}

% frequently used symbols
\newcommand{\bE}{\mathbb{E}}
\newcommand{\bP}{\mathbb{P}}
\newcommand{\bQ}{\mathbb{Q}}
\newcommand{\bR}{\mathbb{R}}
\newcommand{\bS}{\mathbb{S}}
\newcommand{\bN}{\mathbb{N}}
\newcommand{\bZ}{\mathbb{Z}}
\newcommand{\sC}{{\mathcal C}} 
\newcommand{\sD}{{\mathcal D}} 
\newcommand{\sE}{{\mathcal E}} 
\newcommand{\sF}{{\mathcal F}} 
\newcommand{\sL}{{\mathcal L}} 
\newcommand{\sH}{{\mathcal H}} 
\newcommand{\sN}{{\mathcal N}} 
\newcommand{\sO}{{\mathcal O}} 
\newcommand{\sP}{{\mathcal P}} 
\newcommand{\sR}{{\mathcal R}} 
\newcommand{\sS}{{\mathcal S}}
\newcommand{\sU}{{\mathcal U}} 
\newcommand{\sX}{{\mathcal X}} 
\newcommand{\sY}{{\mathcal Y}} 
\newcommand{\sZ}{{\mathcal Z}}

% operators
\newcommand{\sign}{\mathop{\mathrm{sign}}}
\newcommand{\supp}{\mathop{\mathrm{supp}}} % support
\newcommand{\argmin}{\operatornamewithlimits{arg\ min}}
\newcommand{\argmax}{\operatornamewithlimits{arg\ max}}
\newcommand{\dist}{\operatorname{dist}}
\newcommand{\tr}{\text{tr}}
\newcommand{\vecop}{\text{vec}}
\newcommand{\st}{\operatorname{s.t.}}
\newcommand{\cut}{\setminus}
\newcommand{\ind}[1]{\mathbbm{1}\left\{#1\right\}} 
\newcommand{\given}{\ | \ }

% grouping operators
\newcommand{\brac}[1]{\left[#1\right]}
\newcommand{\set}[1]{\left\{#1\right\}}
\newcommand{\abs}[1]{\left\lvert #1 \right\rvert}
\newcommand{\paren}[1]{\left(#1\right)}
\newcommand{\norm}[1]{\left\|#1\right\|}
\newcommand{\ip}[2]{\left\langle #1,#2 \right\rangle}

% code commands
\newcommand{\matlab}{\textsc{Matlab }}
\newcommand{\python}{\textsc{Python }}

% header command
\newcommand{\homework}[3]{
    \pagestyle{myheadings}
    \thispagestyle{plain}
    \newpage
    \setcounter{page}{1}
    \setlength{\headsep}{10mm}
    \noindent
    \begin{center}
    \framebox{
        \vbox{\vspace{2mm}
            \hbox to 6.28in { {\bf EE 516: Mathematical Foundations of Machine Learning
            \hfill Winter \the\year} }
        \vspace{4mm}
        \hbox to 6.28in { {\Large \hfill Homework #1 \hfill} }
        \vspace{2mm}
        \hbox to 6.28in { \Large \hfill Due: #2, 11:59PM PT \hfill }
        \vspace{2mm}
        \hbox to 6.28in { {\it Student Name: #3} \hfill {\it Instructor Name: John Lipor}}
        \vspace{2mm}}
   }
   \end{center}
   \markboth{Homework #1}{Homework #1}
   \vspace*{4mm}
}

% For problem titles
\titleformat{\section}{\normalfont\bf}{Problem \thesection}{1em}{}

\begin{document}

% PUT YOUR NAME BELOW WHERE I HAVE \X
\homework{4}{February 12, \the\year}{\X}
% PUT YOUR NAME BELOW WHERE I HAVE \X

\noindent Since this is an exam week, I'm giving you two extra days to complete the homework. However, Homework 5 will still be due on Friday, February 17, so it is best to complete this before the weekend if possible.

%\section{\normalfont{(4 pts)}}

%For $A \in \bR^{m \times n}$, $y \in \bR^{m}$, and $x \in \bR^{n}$, check directly that $(I - A^{+}A)x$ and $A^{+}y$ are orthogonal vectors, where $A^{+}$ denotes the Moore-Penrose pseudoinverse of the matrix $A$.

%\section{\normalfont{(5 pts each)}}

%Consider the set of measurements $\left( x_{i}, y_{i} \right)$: \quad $(1,2)$, \quad $(2,1)$, \quad $(3,3)$.

%\begin{enumerate}[(a)]
    %\item Find the best (in the 2-norm sense) line of the form $y = \alpha x + \beta$ that fits this data.
    %\item Find the best (in the 2-norm sense) line of the form $x = \gamma y + \delta$ that fits this data.
        
        %\textbf{Hint:} You may want to re-use your answer from (a).
%\end{enumerate}

\section{\normalfont{(4 pts each)}}

Let $f(t) = 0.5 e^{0.8t}$, $t \in [0,2]$.

\begin{enumerate}[(a)]
    \item Suppose 16 exact measurements of $f(t)$ are available to you, taken at the times $t$ in the array $T$ generated using the command \texttt{T = linspace(0,2,16)}. Use \matlab or \python to generate 16 exact measurements $y_{i} = f(t_{i})$, $i = 1,\dots,16$. Write a script to determine the coefficients of the least-squares polynomial approximation of the data for
        \begin{itemize}
            \item a polynomial of degree 15: $p_{15}(t)$;
            \item a polynomial of degree 2: $p_{2}(t)$.
        \end{itemize}
        Compare the resulting functions with the true function by plotting all three functions over a fine grid on the same figure. Turn in your plot and a sentence or two summarizing the results qualatatively.

    \item Now generate noisy measurements $y_{i} = f(t_{i}) + n_{i}$, where $n_{i}$ can be generated as \texttt{n = randn(length(T))}. Write another script to determine the coefficients of the least-squares polynomial approximation of the data for
        \begin{itemize}
            \item a polynomial of degree 15: $p_{15}(t)$;
            \item a polynomial of degree 2: $p_{2}(t)$.
        \end{itemize}
        Compare the resulting functions with the true function by plotting all three functions over a fine grid on the same figure. Turn in your plot and a sentence or two summarizing the results qualatatively.
\end{enumerate}

\section{\normalfont{(5 pts each)}}

Let
\begin{equation*}
    A = \begin{bmatrix} 1 & 0 \\ 0 & 0 \end{bmatrix}, \qquad b = \begin{bmatrix} 1 \\ 1 \end{bmatrix}.
\end{equation*}
The resulting least-squares solution is
\begin{equation*}
    \hat{x} = A^{+}b = \begin{bmatrix} 1 & 0 \\ 0 & 0 \end{bmatrix} \begin{bmatrix} 1 \\ 1 \end{bmatrix} = \begin{bmatrix} 1 \\ 0 \end{bmatrix}.
\end{equation*}
\begin{enumerate}[(a)]
    \item Consider a perturbation $E_{1} = \begin{bmatrix} 0 & \delta \\ 0 & 0 \end{bmatrix}$ of $A$, where $\delta$ is a small positive number. Solve the perturbed version of the problem: $\hat{z} = \argmin_{z} \norm{A_{1}z - b}_{2}^{2}$, where $A_{1} = A + E_{1}$. What happens to $\norm{\hat{x} - \hat{z}}_{2}$ as $\delta$ approaches 0?
    \item Now consider a perturbation $E_{2} = \begin{bmatrix} 0 & 0 \\ 0 & \delta \end{bmatrix}$, where again $\delta$ is a small positive number. Solve the perturbed version of the problem: $\hat{z} = \argmin_{z} \norm{A_{2}z - b}_{2}^{2}$, where $A_{2} = A + E_{2}$. What happens to $\norm{\hat{x} - \hat{z}}_{2}$ as $\delta$ approaches 0?
\end{enumerate}

\section{\normalfont{(5 pts each)}}

Recall the least-squares problem
\begin{equation*}
    \hat{x} = \argmin_{x} \norm{Ax - b}_{2}^{2}.
\end{equation*}
When $A$ is large, computing the pseudoinverse can be computationally prohibitive due to the cost of the SVD operation. In such settings, it is common to apply an optimization tool called \href{https://en.wikipedia.org/wiki/Gradient_descent}{\emph{gradient descent}}, which follows the iteration
\begin{equation*}
    x_{k+1} = x_{k} - \mu \nabla f(x_{k}),
\end{equation*}
where $\mu$ is the step size and $\nabla f(x_{k})$ is the direction of descent. In the case where $f(x) = \frac{1}{2} \norm{Ax - b}_{2}^{2}$, we have
\begin{equation*}
    \nabla f(x) = A^{T}Ax - A^{T}b.
\end{equation*}
Hence, the iteration given by
\begin{equation}
    x_{k+1} = x_{k} - \mu A^{T}(Ax_{k} - b)
    \label{eq:lsgd}
\end{equation}
will minimize $\norm{Ax - y}_{2}^{2}$ as the iteration $k \rightarrow \infty$ when $0 < \mu < 2/\sigma_{1}^{2}(A)$, where $\sigma_{1}(A)$ is the largest singular value of $A$. Note that the iteration has a \emph{fixed point}, i.e., $x_{k+1} = x_{k}$ when
\begin{equation*}
    A^{T}(Ax_{k} - b) = 0,
\end{equation*}
which are exactly the normal equations. So any fixed point minimizes the least-squares cost function!

\begin{enumerate}[(a)]
    \item Complete the function called \texttt{lsgd} (in files for download) that implements the above least-squares gradient descent algorithm described by \eqref{eq:lsgd}. Test your code using the included \texttt{prob5a} file and report the output. What do these values measure?
    \item After your implementation is complete, use it to generate a plot of $\norm{x_{k} - \hat{x}}_{2}$ as a function of $k$ for $A$ and $b$ as generated in \texttt{prob5b}. Test your algorithm for $\mu = c/\sigma_{1}^{2}(A)$ for $c \in \set{0.1, 1, 1.9, 2}$ and turn in one plot with all four curves on it. What do you observe about the convergence as a function of the step size?
\end{enumerate}

\section{\normalfont{(5 pts)}}

In this final installment of the photometric stereo problem, you will use your code from the previous two homework assignments to reconstruct a 3D surface from 2D images. In Homework 3, you computed the normal vectors to a particular surface from images by formulating and solving a least-squares problem. From vector calculus, we can express these surface normal vectors as
\begin{equation}
    n(x,y) = \frac{1}{\sqrt{1 + \left( \frac{\partial}{\partial x} f(x,y) \right)^{2} + \left( \frac{\partial}{\partial y} f(x,y) \right)^{2}}} \begin{bmatrix} - \frac{\partial}{\partial x} f(x,y) \\ - \frac{\partial}{\partial y} f(x,y) \\ 1 \end{bmatrix},
    \label{eq:normal}
\end{equation}
where $\frac{\partial}{\partial x} f$ and $\frac{\partial}{\partial y} f$ denote the partial derivatives of $f(x,y)$ with respect to $x$ and $y$, respectively. From \eqref{eq:normal}, we can compute the partial derivatives as
\begin{equation}
    \frac{\partial}{\partial x} f(x,y) = -\frac{n_{1}(x,y)}{n_{3}(x,y)} \qquad \frac{\partial}{\partial y} f(x,y) = - \frac{n_{2}(x,y)}{n_{3}(x,y)},
    \label{eq:partials}
\end{equation}
where
\begin{equation}
    n(x,y) = \begin{bmatrix} n_{1}(x,y) & n_{2}(x,y) & n_{3}(x,y) \end{bmatrix}^{T}.
    \label{eq:nxy}
\end{equation}

In Homework 2, you constructed the $A$ matrix satisfying
\begin{equation*}
    \begin{bmatrix} \texttt{dfdx} \\ \texttt{dfdy} \end{bmatrix} = A \; \texttt{fxy},
\end{equation*}
where \texttt{dfdx} and \texttt{dfdy} denote the vectorized approximations of the partial derivatives and \texttt{fxy} was the vectorized approximation of $f(x,y)$. Using \eqref{eq:partials}, we can compute \texttt{dfdx} and \texttt{dfdy} from our normal vectors, and using our $A$ matrix, we can solve the following least-squares problem
\begin{equation}
    \texttt{fxy} = \argmin_{f \in \bR^{mn}} \norm{\begin{bmatrix} \texttt{dfdx} \\ \texttt{dfdy} \end{bmatrix} - Af}_{2}^{2}
    \label{eq:fxy}
\end{equation}
to obtain the surface corresponding to our normal vectors.

\textbf{Your task} is to incorporate the functions \texttt{compute\_normals} and \texttt{first\_diffs\_2d\_matrix} into the \texttt{catdemo} script and formulate \eqref{eq:fxy} using the outputs. Once \texttt{A} and \texttt{b} are appropriately defined in the script, this will generate a surface view of the object, completing the project. Turn in this plot, which should match the below if using \matlab (the surface in \python is less sharp).

After $f(x,y)$ is estimated, one can separately generate a stereolithography file that can be rendered on a 3D display or printed by a 3D printer. The \matlab code to generate such a file is also included.

%\begin{figure}[h]
    %\centering
    %\includegraphics[width=6in]{catFigs.eps}
    %\caption{\matlab figure generated from Problem 8.}
    %\label{fig:catFigs}
%\end{figure}

\end{document}
