\documentclass{article}

\usepackage[margin=1in]{geometry}
\usepackage{amsmath, amsfonts}
\usepackage{enumitem}
\usepackage{graphicx}
\usepackage{titling}
\usepackage{amssymb}
\usepackage{url}
\usepackage{xcolor}
\usepackage{subcaption}
\usepackage[colorlinks=true,urlcolor=blue]{hyperref}
\usepackage{minted}
\usepackage{etoolbox, environ, verbatim, textcomp, times}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Commands for customizing the assignment %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\newcommand \duedate {5:30 p.m. Wednesday, November 2, 2016}

\title{
  10-601 Machine Learning: Homework 5\\
  \vspace{0.2cm}
  \large{
    Out: Monday, October 24, 2016\\
    Due \duedate{} \\
    TAs: Ben Cowley, Varshaa Naganathan, Simon Shaolei Du
  }
}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Useful commands for typesetting the questions %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\newcommand \expect {\mathbb{E}}
\newcommand \mle [1]{{\hat #1}^{\rm MLE}}
\newcommand \map [1]{{\hat #1}^{\rm MAP}}
\newcommand \argmax {\operatorname*{argmax}}
\newcommand \argmin {\operatorname*{argmin}}
\newcommand \code [1]{{\tt #1}}
\newcommand \datacount [1]{\#\{#1\}}
\newcommand \ind [1]{\mathbb{I}\{#1\}}
\newcommand \norm [1]{\Vert#1\Vert}
\newcommand \reals {\mathbb{R}}


\def\u{\mathbf u}
\def\x{\mathbf x}
\def\y{\mathbf y}
\def\z{\mathbf z}
\def\A{\emph{A}}
\def\B{\emph{B}}
\def\C{\emph{C}}
\def\D{\emph{D}}
\def\U{\mathbf U}


%%%%%%%%%%%%%%%%%%%%%%
% Toggling Solutions %
%%%%%%%%%%%%%%%%%%%%%%

\newtoggle{solutions}
\toggletrue{solutions}

\NewEnviron{solution}{%
  \iftoggle{solutions}{%
    \color{blue} {\bf Solution:} \BODY%
  }{}%
}

%%%%%%%%%%%%%%%%%%%%%%%%%%
% Document configuration %
%%%%%%%%%%%%%%%%%%%%%%%%%%

% Don't display a date in the title and remove the white space
\predate{}
\postdate{}
\date{}

% Don't display an author and remove the white space
\preauthor{}
\postauthor{}
\author{}

%%%%%%%%%%%%%%%%%%
% Begin Document %
%%%%%%%%%%%%%%%%%%

\begin{document}



\section{\emph{k}-means Clustering [60 pts]}

In this problem you will implement Lloyd's method for the $k$-means clustering
problem and answer several questions about the $k$-means objective, Lloyd's
method, and $k$-means++.

Recall that given a set $X = \{x_1, \dots, x_n\} \subset \reals^d$ of $n$ points
in $d$-dimensional space, the goal of $k$-means clustering is to find a set of
centers $c_1$, \dots, $c_k \in \reals^d$ that minimize the $k$-means objective:
%
\begin{equation}
\label{eq:kmeans_objective}
\sum_{j = 1}^n \min_{i \in \{1, \dots, k\}} \norm{x_j - c_i}^2,
\end{equation}
%
which measures the sum of squared distances from each point $x_j$ to its nearest
center.




\vspace{0.2in}

In class we discussed that finding the optimal centers for the $k$-means
objective is NP-hard, which means that there is likely no algorithm that can
efficiently compute the optimal centers. Instead, we often use Lloyd's method,
which is a heuristic algorithm for minimizing the $k$-means objective that is
efficient in practice and often outputs reasonably good clusterings. Lloyd's
method maintains a set of centers $c_1$,~\dots,~$c_k$ and a partitioning of the
data $X$ into $k$ clusters, $C_1$,~\dots,~$C_k$. The algorithm alternates
between two steps: (i) improving the partitioning $C_1$,~\dots,~$C_k$ by
reassigning each point to the cluster with the nearest center, and (ii)
improving the centers $c_1$, \dots, $c_k$ by setting $c_i$ to be the mean of
those points in the set $C_i$ for $i = 1, \dots, k$. Typically, these two steps
are repeated until the clustering converges (i.e, the partitioning
$C_1$,~\dots,~$C_k$ remains unchanged after an update). Pseudocode is given
below:
%
\begin{enumerate}[label=\arabic*.]
\item Initialize the centers $c_1$, \dots, $c_k$ and the partition $C_1$, \dots, $C_k$ arbitrarily.
\item Do the following until the partitioning $C_1$, \dots, $C_k$ does not change:
\begin{enumerate}[label=\roman*.]
\item For each cluster index $i$, let $C_i = \{ x \in X \,:\, \hbox{$x$ is closer to $c_i$ than any other center}\}$, breaking ties arbitrarily but consistently.
\item For each cluster index $i$, let $c_i = \frac{1}{|C_i|} \sum_{x \in C_i} x$.
\end{enumerate}
\end{enumerate}

\subsection{Implementing Lloyd's Method}

In the remainder of this problem you will implement and experiment with Lloyd's
method and the $k$-means++ algorithm on the two dimensional dataset shown in
Figure~\ref{fig:2ddata}.  \textbf{Note:} Use only native Octave functions---higher-level
functions, such as \code{pdist2}, are not supported by Autolab.  Autolab should take \texttildelow 1~min to grade your submission.  If it is taking longer, there is an error, most likely with convergence.

\begin{figure}[h]
  \begin{center}
  \includegraphics[width=0.6\textwidth]{figures/clusters_proj.pdf}
  \end{center}
  \caption{The 2D dataset $X$ used for problem 1.}
  \label{fig:2ddata}
\end{figure}

\begin{enumerate}[label=(\alph*)]
%\setcounter{enumi}{1}

\item \label{q2:istart} \relax{[8 pts]} Complete the function {\tt [a] =
update\_assignments(X, C, a)}. The input {\tt X} is the $n \times d$ data
matrix, {\tt C} is the $k \times d$ matrix of current centers, and {\tt a} is
the $n \times 1$ vector of current cluster assignments. That is, {\tt C(i,:)} is
the center for cluster $i$, and {\tt a(j)} denotes the cluster label $(1, \ldots, k)$ for the $j^{\rm th}$ data point. Your function should output a new $n \times 1$
vector of cluster assignments so that each point is assigned to the cluster with
the center that has the smallest Euclidean distance with the point.

\begin{solution}
  Submit your code to Autolab.
\end{solution}

\item \relax{[8 pts]} Complete the function {\tt [C] = update\_centers(X, C,
a)}. The input arguments are as in part \ref{q2:istart}. Your function should
output a $k \times d$ matrix {\tt C} whose $i^{\rm th}$ row is the optimal cluster
center for those points in cluster $i$.

\begin{solution}
  Submit your code to Autolab.
\end{solution}

\item \relax{[8 pts]} Complete the function {\tt [C, a] = lloyd\_iteration(X,
$\texttt{C}_0$)}. This function takes a data matrix {\tt X}, initial centers $\texttt{C}_0$ and runs
Lloyd's method until convergence. Specifically, alternate between updating the
assignments and updating the centers until the the assignments stop changing.
Your function should output the final $k \times d$ matrix {\tt C} of centers and
final the $n \times 1$ vector {\tt a} of assignments.

\begin{solution}
  Submit your code to Autolab.
\end{solution}

\item \label{q2:iend} \relax{[8 pts]} Complete the function {\tt [obj] =
kmeans\_obj(X, C, a)}. This function takes the $n \times d$ data matrix {\tt X},
a $k \times d$ matrix {\tt C} of centers, and a $n \times 1$ vector {\tt a} of
cluster assignments. Your function should output the value of the $k$-means
objective of the provided clustering.

\begin{solution}
  Submit your code to Autolab.
\end{solution}

\end{enumerate}

We provide you with a function {\tt [C, a, obj] = kmeans\_cluster(X, k, init,
num\_restarts)} that takes an $n \times d$ data matrix {\tt X}, the number {\tt
k} of clusters, a string {\tt init} which must be either {\tt `random'} or {\tt
`kmeans++'}, and a number of restarts {\tt num\_restarts}. This function runs
Lloyd's method {\tt num\_restarts} times and outputs the best clustering it
finds. You may use this function when answering the following questions. When
{\tt init} is set to {\tt `random'} then the centers are initialized uniformly
at random, and when {\tt init} is set to {\tt `kmeans++'} the centers are
initialized using the $D^2$ weighting of $k$-means++.


\clearpage

\subsection{Experiment 1: The effect of \emph{k} on Lloyd's method}

In parts \ref{q2:lloydstart} and \ref{q2:lloydend} of this problem you will
investigate the effect of the number of clusters $k$ on Lloyd's method and a
simple heuristic for choosing a good value for $k$. You can load the data for
this problem by running the command {\tt load kmeans\_data} in Matlab or Octave.
This will introduce a $900 \times 2$ matrix {\tt X}, where each row corresponds
to a single point in the dataset shown in Figure~\ref{fig:2ddata}.

\begin{enumerate}[resume*]

\item \label{q2:lloydstart} \relax{[6 pts]} Finish coding the script {\tt plot\_kmeans\_sweep\_k.m}, which plots the $k$-means objective value obtained by running {\tt kmeans\_cluster(X, k,
`random', 10)} for each value of {\tt k} in $\{1, \dots, 20\}$. The $x$-axis of
your plot should be the value of {\tt k} used, and the $y$-axis should be the
{\tt k}-means objective.  Include only the plot in your written solutions.


\begin{solution}

    Include one plot in your solution.
%   \begin{center}
%     \includegraphics[width=0.7\textwidth]{./figures/ques1f_kmeans_sweep_k.pdf}
%   \end{center}
  

\end{solution}

\item \label{q2:lloydend} \relax{[6 pts]} One heuristic for choosing the value
of $k$ is to pick the value at the ``elbow'' or bend of the plot produced in
part \ref{q2:lloydstart}, since increasing $k$ beyond this point results in
relatively little gain. Based on your plot, what value of $k$ would you choose?
Does this agree with your intuitions when looking at the data?

\begin{solution}
  Brief explanation here.
\end{solution}

\end{enumerate}


\clearpage

\subsection{Experiment 2: The effect of initialization on Lloyd's method}

In parts \ref{q2:initstart} and \ref{q2:initend} you will investigate the effect
of the center initialization on Lloyd's method. In particular, you will compare
random initialization with the $k$-means++ initialization. Note that the
provided {\tt kmeans\_cluster} function has both kinds of initialization already
implemented. In this problem, you just need to run the code and answer
questions.

\begin{enumerate}[resume*]

\item \label{q2:initplot} \relax{[4 pts]} For both types of intialization, use the script {\tt plot\_kmeans\_init\_projs.m} to
plot projections similar to Figure~\ref{fig:2ddata} except that the points are colored based on
their class label, as well as the optimized centers.  Include both plots in your written
solutions.

\begin{solution}
    
    Include two plots in your solution.
    
%   \begin{center}
%     \includegraphics[width=0.7\textwidth]{./figures/ques1h_kmeans_random_proj.pdf}
%   \end{center}
  
%   \begin{center}
%     \includegraphics[width=0.7\textwidth]{./figures/ques1h_kmeans_plusplus_proj.pdf}
%   \end{center}

\end{solution}



\item \label{q2:initstart} \relax{[6 pts]} For each value of {\tt init} in $\{
\hbox{\tt `random'}, \hbox{\tt `kmeans++'} \}$, report the mean and standard deviation of the $k$-means
objective value returned by {\tt kmeans\_cluster(X, 9, init, 1)} over 1000 runs.  Note that {\tt num\_restarts} = 1. Since the initialization is
random, your answers will vary each time you perform this experiment. However, there should be a clear gap in performance between the two initialization methods.  Report the values and the script code.

\begin{solution}

  The output of the above code gives average objective value of
  $XX.X \pm XX.X$ for random initialization and $XX.X \pm XX.X$ for $k$-means++ initialization.
  
  Place your code here:
%   \color{black}
%   \inputminted{octave}{./code_solutions/init_study.m}
  
\end{solution}

\item \label{q2:initend} \relax{[6 pts]} Intuitively, a good center
initialization for the 2D dataset shown in Figure~\ref{fig:2ddata} is when we
pick one center from each of the $9$ natural Gaussian clusters. When two centers
land in the same Gaussian cluster, they will split that cluster in half and
leave only $7$ centers to cover the remaining $8$ Gaussian clusters. Explain in
1--2 sentences why using $k$-means++ initialization is more likely to choose one
sample from each cluster than random initialization.

\begin{solution}
  Brief explanation goes here.
\end{solution}

\end{enumerate}



\clearpage


\section{Z-scoring before applying PCA? [20 pts]}

\noindent Consider the $3 \times N$ matrix $X$ be a dataset with 3 variables and $N$ samples.  PCA follows these steps when applied to $X$:

\begin{enumerate}[leftmargin=2\parindent]
\item recenter $X$ such that the sample mean $\frac{1}{N}\sum_{n=1}^N X_n = \mathbf{0}$
\item identify ordered orthonormal principal components $\u_1, \u_2, \u_3 \in \mathbb{R}^3$
\item project recentered $X$ onto the principal components: $\u_i^T X$
\item identify ordered variances $\sigma_1, \sigma_2, \sigma_3$ for the principal components, \\ where $\sigma_i = \frac{1}{N} \sum\limits_{n=1}^N (\u_i^T X_n)^2$ and $\sigma_1 \geq \sigma_2 \geq \sigma_3$

\end{enumerate}
We can also think about PCA as performing an eigendecomposition on the $3 \times 3$ covariance matrix $\Sigma$.  Given the eigenvectors $\u_1, \u_2, \u_3$ (i.e., the basis vectors), we can compute each eigenvalue (i.e., the variances) as $\sigma_i = \u_i^T \Sigma \u_i$.  PCA orders the eigenvalues such that $\sigma_1 \geq \sigma_2 \geq \sigma_3$.\\


\noindent Now consider drawing the columns of $X = [X_1, \ldots, X_N]$ from the given multivariate Gaussian distribution:\\
\begin{center}
$
    X_n \sim
    \mathcal{N}\left(
        \mu = \begin{bmatrix}
            5 \\
            5 \\
            10 \\
        \end{bmatrix}
        ,
        \Sigma = \begin{bmatrix}
            6 & -3 & 0 \\
            -3 & 6 & 0 \\
            0 & 0 & 10
        \end{bmatrix}
    \right)
    $
\end{center}

\vspace{0.2in}

\noindent We apply PCA to $X$ and ask the following questions.\\

\noindent \textbf{Note:} For the following problems, feel free to use Matlab's built-in \texttt{pca}, \texttt{mvnrnd}, \texttt{cov2corr}, etc., to help guide your answers.  Use $N \geq 1000$ samples when generating the data $X$. Do not include code in your answers.  Provide analytical, exact solutions; do \emph{not} give numerical solutions from Matlab.

\noindent \textbf{Hint:} The fastest way to solve these problems is with intuition and to check with Matlab.

\begin{enumerate}


\item \relax{[4 pts.]} What are the ordered variances $\sigma_1, \sigma_2, \sigma_3$ for the ordered principal components?

\begin{solution}
    \\
    $\sigma_1 = XXX$\\
    $\sigma_2 = XXX$\\
    $\sigma_3 = XXX$\\
\end{solution}

\item \relax{[4 pts.]} What are the ordered basis vectors $\u_1, \u_2, \u_3 \in \mathbb{R}^3$ for the ordered principal components?


\begin{solution}
    
    \begin{center}
    $[\u_1, \u_2, \u_3] = 
        \begin{bmatrix}
            XXX & XXX & XXX\\
            XXX & XXX & XXX\\
            XXX & XXX & XXX
        \end{bmatrix}
    $
    \end{center}
\end{solution}



\end{enumerate}


\noindent One limitation of PCA is that it is scale variant.  For example, consider two variables, where one represents age (in years) and one represents height (in feet).  Because the range of ages (0 to 100 years) is much larger than the range of height (4 to 7 feet), PCA will trivially identify the age variable for the first principal component.  However, the first principal component should be able to reconstruct both age and height, which are likely correlated.\\

\noindent To overcome this, a typical procedure before applying PCA is to z-score each variable (i.e., for each $x_i$, subtract its mean and divide by its standard deviation).  This transforms the covariance matrix $\Sigma$ into a correlation matrix $\tilde{\Sigma}$.  Note that $\hat{\Sigma}_{ij} = \frac{\Sigma_{ij}}{\sqrt{\Sigma_{ii}} \sqrt{\Sigma_{jj}}}$.

\begin{enumerate}[resume]
\item \relax{[3 pts.]} Let's say we z-score $X$ to obtain $\tilde{X}$.  Then, $\tilde{X} = [\tilde{x}_1, \tilde{x}_2, \tilde{x}_3]^T \sim \mathcal{N}(\tilde{\mu}, \tilde{\Sigma})$.  What are the elements of $\tilde{\mu}$ and $\tilde{\Sigma}$?

\begin{solution}
Note that \begin{center}
\begin{align*}
\tilde{X} = \begin{bmatrix}
(x_1 - XXX)/XXX\\
(x_2-XXX)/XXX\\
(x_3-XXX)/XXX
\end{bmatrix}.
\end{align*}
\end{center}
Now using the definition of mean and covariance matrix, we obtain
    \begin{center}
        $\tilde{X} = 
            \begin{bmatrix}
                \tilde{x}_1 \\
                \tilde{x}_2 \\
                \tilde{x}_3
            \end{bmatrix}
            =   \mathcal{N} \left (
            \tilde{\mu} = \begin{bmatrix}
                XXX \\
                XXX \\
                XXX \\
            \end{bmatrix},
            \tilde{\Sigma} = \begin{bmatrix}
                XXX & XXX & XXX \\
                XXX & XXX & XXX \\
                XXX & XXX & XXX
            \end{bmatrix}
            \right )
            $
            
    \end{center}
\end{solution}


\end{enumerate}

\noindent We now apply PCA to $\tilde{X}$.  

\begin{enumerate}[resume]
\item \relax{[3 pts.]} What are the ordered variances $\tilde{\sigma}_1, \tilde{\sigma}_2, \tilde{\sigma}_3$ for the ordered principal components?

\begin{solution}
    \\
    $\sigma_1 = XXX$\\
    $\sigma_2 = XXX$\\
    $\sigma_3 = XXX$
\end{solution}

\item \relax{[3 pts.]} What are the ordered basis vectors $\tilde{\u}_1, \tilde{\u}_2, \tilde{\u}_3 \in \mathbb{R}^3$ for the ordered principal components?

\begin{solution}
    $[\tilde{u}_1, \tilde{u}_2, \tilde{u}_3] = 
        \begin{bmatrix}
            XXX & XXX & XXX\\
            XXX & XXX & XXX\\
            XXX & XXX & XXX
        \end{bmatrix}
    $
\end{solution}


\item \relax{[3 pts.]} Do your answers for 4 and 5 differ from your answers to 1 and 2?  Why or why not?

\begin{solution}
    Brief explanation here.
\end{solution}


\end{enumerate}

\clearpage

\section{Lots of dots [20 pts]}

\begin{enumerate}
    \item \relax{[10 pts]} For each of the two data sets shown below, draw vectors that correspond to the first and second principal components that would be produced by Principal Components Analysis (PCA). For each diagram, label which vector is the first component, and which is the second.
    
  \begin{center}
    \includegraphics[width=0.7\textwidth]{./figures/pca_dots_question.pdf}
  \end{center}
  
  \begin{solution}
     Include two plots with PCs drawn with red arrows.
     
    %   \begin{center}
    %     \includegraphics[width=0.7\textwidth]{./figures/pca_dots_solution.pdf}
    %   \end{center}
      
  \end{solution}
  
  \item \relax{[5 pts]} Is it possible for you to draw a third principal component for the first of these datasets (in $A$)? Please do, or explain in one sentence why it is impossible.
  
  \begin{solution}
      Brief explanation here (1 sentence).
  \end{solution}
  
  \item \relax{[5 pts]} PCA is typically performed before clustering high-dimensional data.  For this procedure, what assumption do we make about the data's covariance structure in relation to the distances between clusters?  Explain why this assumption holds or does not hold for the second dataset (in $B$).  Be brief (2-3 sentences).
  
  \begin{solution}
      Brief explanation here (2-3 sentences).
  \end{solution}
\end{enumerate}
\clearpage

\section{PCA two ways [Extra Credit: 10~pts]}

Like the primal and dual problems for SVM, we can have two perspectives of PCA: maximizing the variance or minimizing the
reconstruction error.  We will find that both these views lead to the same optimization problem.
This question uses a lot of linear algebra---make sure the sizes of each matrix and vector correctly match up!

\subsection{\relax{[Extra credit: 5 pts]} PCA: Maximizing the variance}
\label{sec:max_var}
Consider $N$ data points $X_1, \ldots, X_N \in \mathbb{R}^p$, such that the sample mean of $X$ is zero: $\frac{1}{N}\sum\limits_{i=1}^N X_i = \mathbf{0}$.  Also consider projection vector $\u \in \mathbb{R}^p$, where $\|\u\|_2 = 1$.  We would like to maximize the sample variance $\tilde{V}[\u^T X]$, which is the sample variance of $X$ projected onto $\u$.  The sample variance of $N$ samples of a variable $z$ is $\tilde{V}[z] = \frac{1}{N} \sum\limits_{i=1}^N z_i^2$, where $\frac{1}{N}\sum\limits_{i=1}^N z_i = 0$.  We can write an optimization problem to maximize the sample variance:

\begin{equation}
    \label{eqn:max_var}
    \max\limits_{\substack{\u \\ \textrm{s.t. } \|\u\|_2 = 1}} \tilde{V}[\u^T X]
\end{equation}

\noindent Reformulate Eqn.~\ref{eqn:max_var} to the following optimization problem, and define the $p \times p$ matrix $\Sigma$:
\begin{equation}
    \label{eqn:max_Sigma1}
    \max\limits_{\substack{\u \\ \textrm{s.t. } \|\u\|_2 = 1}} \u^T \Sigma \u
\end{equation}
\textbf{Show all work:}

\begin{solution}
    Solution goes here.  Include all steps.
\end{solution}


\clearpage

\subsection{\relax{[Extra credit: 5 pts]} PCA: Minimizing the reconstruction error}
\label{sec:min_error}

Consider the same variables as those defined in Question~\ref{sec:max_var}.
Instead of maximizing the projected variance, we now seek to minimize the reconstruction error.  In other words,
we would like to minimize the difference between $X$ and the reconstructed $\u\u^T X$. Note that $\u\u^T X$ first projects $X$ onto $\u$, and then projects this scalar back into the $p$-dimensional space along the $\u$ axis.  Minimzing the reconstruction error can be written as the following optimization problem:

\begin{equation}
    \label{eqn:min_error}
    \min\limits_{\substack{\u \\ \textrm{s.t. } \|\u\|_2 = 1}} \frac{1}{N} \sum\limits_{i=1}^N \|X_i - \u\u^T X_i\|_2^2
\end{equation}

\noindent Reformulate Eqn.~\ref{eqn:min_error} to the following optimization problem, with the same $\Sigma$ as defined in Section~\ref{sec:max_var}.
\begin{equation}
    \label{eqn:max_Sigma2}
    \max\limits_{\substack{\u \\ \textrm{s.t. } \|\u\|_2 = 1}} \u^T \Sigma \u
\end{equation}
\textbf{Show all work:}

\begin{solution}
    Solution goes here.  Include all steps.
\end{solution}


\end{document}