\documentclass[11pt]{article}
\usepackage[latin1]{inputenc}
\usepackage{amsmath}
\usepackage{epsfig}
\usepackage{pgf}
\usepackage{graphicx}
 \usepackage{amssymb}
\usepackage{url}
\def\urltilda{\kern -.15em\lower .7ex\hbox{\~{}}\kern .04em}

    \oddsidemargin  0.0in
    \evensidemargin 0.0in
    \textwidth      6.5in
    \headheight     0.0in
        \headsep        0.0in
    \topmargin      0.0in
    \textheight=9.0in

\newcommand{\Xv}{\mathbf{X}}
\newcommand{\Yv}{\mathbf{Y}}
\newcommand{\xv}{\mathrm{x}}
\newcommand{\bv}{\mathrm{\beta}}
\newcommand{\tv}{\mathrm{T}}
\newcommand{\ji}{{ji}}
\newcommand{\xvit}{{\mathrm{x_i}^{\mathrm{T}}}}
\newcommand{\xvt}{{\mathrm{x}^{\mathrm{T}}}}
\newcommand{\Xvt}{{\mathbf{X}^{\mathrm{T}}}}
\newcommand{\yv}{\mathrm{y}}
\newcommand{\zv}{\mathrm{z}}
\newcommand{\ev}{\mathrm{e}}
\newcommand{\betav}{\mathrm{\beta}}
\newcommand{\tauv}{\mathrm{\tau}}
\newcommand{\thetav}{\mathrm{\theta}}
\newcommand{\gammav}{\mathrm{\gamma}}
\newcommand{\mbs}{\textbf}
\newcommand{\DD}{\scriptscriptstyle}

\title{10-701 Machine Learning, Fall 2011: Homework 2}
\author{Due 10/17 at the beginning of class.}
\date{~}
\begin{document}
\maketitle


\section{Linear regression, model selction [25pt, Nan Li]}

In linear regression, we are given training data of the form,
$\mathcal{D}= (\Xv,\yv) = \{(\xv_{i},y_{i})\}, i=1,2,...,n,$ where
$\xv_{i} \in \mathcal{R}^{1 \times p}$, i.e. $\xv_{i}=
(x_{i,1},\cdots,x_{i,p})^{\tv}$, $y_{i} \in \mathcal{R}$, $\Xv \in
\mathcal{R}^{n \times p}$, where $n$ is number of samples, $p$ is
number of features. Each row $i$ of $\Xv$ is $\xvit$, and $\yv=
(y_1,\cdots,y_n)^{\tv}$.

Least square regression seeks to find $\bv$ that minimizes the
square-error, i.e.:

\begin{displaymath}
\mathrm{J}_1(\bv) = \sum_{i}{(y_i - \xvit \bv)^2} =
(\Xv\bv-\yv)^{\tv}(\Xv\bv-\yv).
\end{displaymath}
It has an unique solution
\begin{equation}\label{eq:1}
\hat{\bv} = (\Xv^T\Xv)^{-1}\Xv^T\yv.
\end{equation}
where $\hat{\beta}$ is called an estimator of $\beta$. An
``estimator'' is a statistic of your data (i.e. a function of your
data) which is intended to approximate a parameter of the underlying
distribution. There is a research field called ``estimation
theory'', which deals with constructing estimators that have nice
properties, like converging to the correct parameter given enough
data, and giving confidence intervals. In this problem we will
explore how regularization affects the bias and variance of the
least square regression model.

Let's assume that $p < n$, and $\Xv^T\Xv$ is invertible. Also assume
that our data is generated from a true model of the form: $y_i =
\xv_i \bv + \epsilon_{i}$ ( or in matrix form $\yv=\Xv \beta +
\epsilon$), where $\epsilon_1,...,\epsilon_n$ are IID and sampled
from a Gaussian with 0 mean and constant standard deviation, that is
$\epsilon_i \sim \mathcal{N}(0,\sigma^2)$ (or $\epsilon \sim
\mathcal{N}(0,\sigma^2 I )$).

%\subsection{The hat matrix[5 pt]}
%
%If we apply our model back to the training data, we have the
%estimated $y$ values as
%$\hat{\yv}=\Xv\hat{\beta}=\Xv(\Xv^T\Xv)^{-1}\Xv^T\yv$.
%$H=\Xv(\Xv^T\Xv)^{-1}\Xv^T$ is often called the hat matrix, because
%it can put a hat on $y$. Now please show that $H$ is the projection
%matrix onto $col(\Xv)$, the column space of $\Xv$. In other words,
%show that $H$ satisfies:
%
%(a) $H$ is symmetric.
%
%(b) $H$ is idempotent: $H^2 = H$.
%
%(c) For any vector $v$, $Hv=v$ iif  $v\in col(\Xv)$.
%
%\textbf{Hint:} proof is easier, if you do a Singular Value
%Decomposition (SVD) to $\Xv$ first.

% \subsection{Bias-Variance Trade-off}
Let's first check how regularization could change the bias and variance of the estimator $\hat{\beta}$.

%\subsection{Least square regression [6 pt]}

%Show that $\hat{\beta}$ has Gaussian distribution and write down its
%mean $\mu$ and covariance matrix $\Sigma$. You  will see that least
%square regression is unbiased $E[\hat{\beta}]=\beta$

\subsection{Ridge regression [10 pt]}

We know that the original linear regression without any regularization is unbiased, i.e.  $E[\hat{\beta}]=\beta$. What about ridge regression? The solution to ridge regression is
\begin{equation}\label{eq:1}
\hat{\bv} = (\Xv^T\Xv +\lambda I)^{-1}\Xv^T\yv.
\end{equation}
Show that $\hat{\beta}$ has Gaussian distribution and write down its
mean $\mu$ and covariance matrix $\Sigma$. You will see that ridge
regression is biased $E[\hat{\beta}] \neq \beta$.

\textbf{Hint:} if $\epsilon \sim \mathcal{N}(0,I)$ then $A\epsilon
\sim \mathcal{N}(0,AA^T)$, where $A$ is any matrix.

\textbf{Hint:} notation is simpler, if you do a Singular Value
Decomposition (SVD) to $\Xv$ first.

\subsection{Extra features [10 pt]}
Now, let's move on to see the relation between model complexity and overfitting. Can we keep increasing our model quality by adding more features in it?

%\subsubsection{Empirical Risk}
The objective of least square regression is to minimize $\mathrm{J}_1(\bv)$ given $\mathcal{D}$.  Assume that besides the existing feature data, $\Xv$, we have, we now have access to one more new feature. The feature values for the given samples are $X_{new}$, where $X_{new} \in \mathcal{R}^{n \times 1}$. Taken together we have a new feature matrix $\Xv_{new} = [\Xv~X_{new}]$. For simplicity, let's assume $\Xv_{new}^T\Xv_{new}=I$. Show that with this new feature, the minimized objective value $\mathrm{J}_1(\bv)$ will be decreased.

\subsection{Model Complexity [3 pt]}
Does this the above result mean that as long as we keep increasing features, we can always get a better model? Briefly explain why or why not.

\subsection{Overfitting [2 pt]}
Please suggest one method to address overfitting.

%\subsubsection{True Risk}
%But what we really want is to predict well on \emph{future} data, not on the training set. Hence, instead of minimizing the empirical risk, we would like to minimize the true risk, $E [ (\Xv\bv-\yv)^{\tv}(\Xv\bv-\yv) ]$.
%
%Suppose that only the first $p^*$ features among all features are relevant, that is, the underlying model for $y$ is $ y = x_{p*} \beta^*+\epsilon$, where $x_{p^*} = (x_1, x_2, ... x_p^*)$ is the feature values of the first $p^*$ dimensions.
%
%Suppose also that the underlying model for $y$ is $ y = x \beta^*+\epsilon$, where $\epsilon \sim\mathcal{N}(0,\sigma^2)$, 
%$\beta^*$ is a $p$-dimensional vector that is non-zero only on the first $p^* < p$ coordinates.
%
%Let $\hat{\beta} \in \mathcal{R}^{p \times 1}$ be any estimator of $\beta^*$. Suppose that $\hat{\beta}$ is non-zero on some coordinate $j > p^*$. Prove then that the average error on future data $E_{xy} [ (y - x \hat{\bv})^2 ]$ is minimized at $\bv^*$.
%
%Hint: $\mathbf{I-XX^T}$ is positive semidefinite.

%\subsubsection{Explanation}
%Please briefly explain what do the results of the two above questions mean, and why.


\section{Learning Theory [25 points + 3 extra points, Bin]}

\subsection{PAC Learning}
PAC stands for ``Probably Approximately Correct'' and concerns a nice formalism for deciding how much data you need to collect in order for a given classifier to achieve a given probability of correct predictions on a given fraction of future test data.
\begin{itemize}
\item  \textbf{[5 points]} True or false: (if true, give a 1 sentence justification; if false, give a counter example.)
Within the setting of PAC learning, it is \textbf{impossible} to assure with probability $1$ that the concept will be learned perfectly (i.e., with true error = 0), regardless of how many training examples are provided.
\end{itemize}


\subsection{VC Dimension}
For hypothesis class $\mathcal{H}$, its VC dimension is at least $d$ if there exists some samples $|S|=d$ which is shattered by $\mathcal{H}$. Please note that this does not mean all samples of size $d$ need to be shattered by $\mathcal{H}$. To show that the VC dimension is at most $d$, one must show that no sample of size $d+1$ could be shattered by $\mathcal{H}$. In this problem, you will calculate the VC-dimension of some hypothesis classes. Remember you need to show the following 2 steps to prove the VC dimension of $\mathcal{H}$ is $d$:
\begin{itemize}
\item There exists a set of $d$ points which can be shattered by $\mathcal{H}$;
\item There exists no set of $d+1$ points which can be shattered by $\mathcal{H}$.
\end{itemize}
You need to show details of how you get the VC dimension in order to get full credit.

\begin{itemize}
\item  [(1)] \textbf{[4 points]} $\mathcal{X}=\mathbb{R}$, $\mathcal{H}$ is the union of 2 intervals; %4
\item  [(2)] \textbf{[4 points]} $\mathcal{X}=\mathbb{R}^2$, $\mathcal{H}$ is the set of axis-parallel squares (with equal height and width); %3
\item  [(3)] \textbf{[4 points]} $\mathcal{X}=\mathbb{R}$, $\mathcal{H}$ is the set of functions defined as following:
\begin{equation}
h(x) = \mathrm{sign}(\sin(ax+b))
\end{equation}
where $a$ and $b$ are parameters in $h$, $\mathrm{sign}(x)=+1$ if $x\ge 0$ and $\mathrm{sign}(x)=-1$ otherwise. %\infty
\item  [(4)] \textbf{[3 extra points]} $\mathcal{H}$ is a finite hypothesis class such that $|\mathcal{H}|<\infty$. Show that the VC dimension of $\mathcal{H}$ is upper bounded by $\log_{2}|\mathcal{H}|$. (Hint: how many different labelings are there for a set of size $\log_{2}|\mathcal{H}|+1$?) % 4
\end{itemize}

\subsection{The Approximation-Estimation-Optimization Tradeoff}
In this problem, we consider a space of input-output pairs $(x,y)\in\mathcal{X}\times \mathcal{Y}$ generated by a probability distribution $P(x,y)$. We define a loss function $l(\hat{y},y)$ (for example, $l(\hat{y},y)=|\hat{y}-y|^2$ as in regression), to measure the discrepancy between the predicted output $\hat{y}$ and the real output $y$. Our target is to find the function $f^*$ that minimizes the expected risk
\begin{equation}
R(f) = \int l(f(x),y)dP(x,y)
\end{equation}
Usually the distribution $P(x,y)$ is unknown, we are instead given a sample $\mathcal{S}$ of $n$ independently drawn training examples $\{(x_i, y_i)\}$, $i=1,\ldots,n$. We define the empirical risk
\begin{equation}
R_n(f)=\frac{1}{n}\sum_{i=1}^{n}l(f(x_i),y_i)
\end{equation}
The learning principle we learned in class consists of first choosing a family $\mathcal{F}$ of candidate prediction functions, then finding the function $f_n=\mathop{\arg\min}_{f\in\mathcal{F}}R_n(f)$. Since the optimal function $f^*$ is usually unlikely to belong to the family $\mathcal{F}$, we also define $f_{\mathcal{F}}^*=\mathop{\arg\min}_{f\in\mathcal{F}}R(f)$. For simplicity, we assume that $f^*$, $f_{\mathcal{F}}^{*}$ and $f_n$ are well defined and unique. We can then decompose the excess error as
\begin{equation}
\mathbb{E}[R(f_n)-R(f^*)]=\mathbb{E}[R(f_{\mathcal{F}}^*)-R(f^*)]+\mathbb{E}[R(f_n)-R(f_{\mathcal{F}}^*)]=\varepsilon_{app}+\varepsilon_{est}
\end{equation}
where the expectation is taken with respect to the random choice of training set. The approximation
error $\varepsilon_{app}$ measures how closely functions in $\mathcal{F}$ can approximate the optimal
solution $f^*$. The estimation error $\varepsilon_{est}$ measures the effect of minimizing the empirical
risk $R_n(f)$ instead of the expected risk $R(f)$.

One flaw of the above decomposition of excess error is that it assumes we find $f_n$ that minimizes the empirical risk $R_n(f)$. However, this procedure is often a computationally expensive operation. Let us assume that our minimization algorithm returns an approximate solution $\tilde{f}_n$ that minimizes the objective function up to a predefined tolerance $\rho\ge0$
\begin{equation}
R_n(\tilde{f}_n)<R_n(f_n)+\rho
\end{equation}
We can then decompose the excess error $\varepsilon=\mathbb{E}[R(\tilde{f}_n)-R(f^*)]$ as
\begin{equation}
\varepsilon=\mathbb{E}[R(f_{\mathcal{F}}^*)-R(f^*)]+\mathbb{E}[R(f_n)-R(f_{\mathcal{F}}^*)]+\mathbb{E}[R(\tilde{f}_n)-R(f_n)]=\varepsilon_{app}+\varepsilon_{est}+\varepsilon_{opt}
\end{equation}
We call the additional $\varepsilon_{opt}$ the optimization error. It reflects the impact of the approximate
optimization on the generalization performance.
\begin{itemize}
\item  [(1)] \textbf{[8 points]} In this question, you will study how approximation error $\varepsilon_{app}$, estimation error $\varepsilon_{est}$, optimization error $\varepsilon_{opt}$ and computation time $T$ change when one of $\{\mathcal{F},n,\rho\}$ increases. (In creasing $\mathcal{F}$ means that the new $\mathcal{F}$ is a superset of the old $\mathcal{F}$) Fill in table \ref{table_bin_1} with $\uparrow$ indicating increase, $\downarrow$ indicating decrease, and $\times$ indicating not affected. Please briefly explain your answer. 

\begin{table}[h]
\begin{center}
\caption{Typical variations when $\mathcal{F}$, $n$ and $\rho$ increase.}
\begin{tabular}{|c||c|c|c|}\hline
&$\mathcal{F}$&$n$&$\rho$\\\hline
$\varepsilon_{app}$&&&$\times$\\\hline
$\varepsilon_{est}$&&&$\times$\\\hline
$\varepsilon_{opt}$&$\times$&$\times$&\\\hline
$T$&&&\\\hline
\end{tabular}\label{table_bin_1}
\end{center}
\end{table}
\end{itemize}


\section{Expectation-Maximization algorithm [Suyash, 25 points]}
In this question, we will derive some details about the expectation-maximization (EM) algorithm and work out an example application.
\subsection{True-False Questions [6 points]}
Explain whether the following statements are true or false with a single line of explanation.
\begin{enumerate}
 \item \textbf{(2 points)} The EM algorithm maximizes the complete log-likelihood.
\item \textbf{(2 points)} Even if the complete log-likelihood is multi-modal, EM cannot get stuck in a local minima.
\item \textbf{(2 points)} The free-energy functional that EM optimizes is a lower bound on the complete log-likelihood.
\end{enumerate}
\subsection{EM and KL divergence [7 points]}
\subsubsection{Optimizing the free-energy functional [4 points]}
Consider the free-energy functional discussed in class that the EM algorithm optimizes. The free-energy functional is given by:
\begin{equation}
 F(q,\theta)=\sum_{q(z|x)} q(z|x) \log \frac{p(z,x|\theta)}{q(z|x)} 
\end{equation}
for any probability distribution $q(z)$. (Note: In this setting, the only random variables are the $z$ variables). In class, we claimed that 
\begin{equation}
 \log p(x|\theta) - F(q,\theta)= KL(q \left \| \right. p(z|x;\theta) )
\label{eq:qandkl}
\end{equation}
where $ KL(a \left \| \right. b) = \sum_x a(x) \log \frac{a(x)}{b(x)} $ is called the Kullback-Leibler divergence between two probability distributions $a$ and $b$. Show that this result is true.
\subsubsection{Optimal value of q(z) in EM [3 points]}
Using the result in Equation~\ref{eq:qandkl}, find the value for $q(z)$ that maximizes $F(q,\theta)$. (Hint: $ KL(a \left \| \right. b) \geq 0$, with equality if and only if $a = b$.)	

\subsection{Example of coin-tossing [12 points]}
Consider the scenario where a sequence of heads and tails are being generated by tossing two coins independently. For each toss, the first coin is chosen with probability $\pi$ and the second coin is chosen with probability $1-\pi$. The head probabilities of the two coins are given by $p_1$ and $p_2$. The probabilities $\pi,p_1$ and $p_2$ are unknown to us and we want to find their values.

Let us represent the result of the $N$ coin tosses by $x_1,x_2,\cdots,x_N$. $x_i=1 $ indicates that that $i^{th}$ toss came out heads and $0$ otherwise. Let  $z_1,z_2,\cdots,z_N$ indicate the coin which was tossed each time. For ease of mathematical analysis, suppose that $z_i=1$ means the first coin was used for the $i^{th}$ toss and $z_i=0$ means the second coin was used for the $i^{th}$ toss. 

\subsubsection{All variables observed [5 points]}
Suppose that we are allowed to observe which coin was used for each toss, i.e, the values of the variables $z_1,\cdots,z_N$ are known.
\begin{enumerate}
\item \textbf{(2 points)} Write the expression for the complete log-likelihood of the N coin tosses. (Hint: For the $i^{th}$ toss, you can write the likelihood as
\begin{eqnarray}
p(x_i,z_i; \pi,p_1,p_2) = {\left[\pi p_1^{x_i} (1-p_1)^{1-x_i} \right]}^{z_i} {\left[(1-\pi) p_2^{x_i} (1-p_2)^{1-x_i}  \right]}^{1-z_i}).
\end{eqnarray}
\item \textbf{(3 points)} The maximum likelihood equations for the parameters in this setting are:
\begin{eqnarray}
\label{eq:pi}
 \pi &=& \frac{\sum_{i=1}^N z_i}{N} \\
 p_1 &=& \frac{\sum_{i=1}^N z_i x_i}{\sum_{i=1}^N z_i} \\
 p_2 &=& \frac{\sum_{i=1}^N (1-z_i)x_i}{\sum_{i=1}^N (1-z_i)}
\end{eqnarray}
Show that evaluating the derivative of the log-likelihood with respect to $\pi$ at the value specified in Equation~\ref{eq:pi} gives a value of zero.
\end{enumerate}

\subsubsection{Only $x$ variables observed [7 points]}
Suppose we are not allowed to observe which coin was used for each toss, i.e, the value of the variables $z_1, \cdots, z_N$ are unknown. We will now use the EM algorithm for finding the values of the unknown parameters.
\begin{enumerate}
\item \textbf{(2 points)} Write the expression for the incomplete log-likelihood.
\item \textbf{(3 points)} 
The E-step for our EM algorithm requires us to compute $q^*(z)=p(z|x;\pi^{(t)},p_1^{(t)},p_2^{(t)})$ where the superscript $t$ indicates the value of the parameters at step $t$. What is the expression for the conditional probability $p(z_i=1|x_i;\pi^{(t)},p_1^{(t)},p_2^{(t)})$ for the random variable $z_i$ associated with the $i^{th}$ toss? Note that your expression should only involve $\pi^{(t)},p_1^{(t)},p_2^{(t)}$ and $x_i$. For notational convenience, you can drop the superscript $t$ from your expression.
\item \textbf{(2 point)} We can show that the updates for the parameters $\pi,p_1$ and $p_2$ in the M-step are given by
\begin{eqnarray}
 \pi^{(t+1)} &=& \frac{\sum_{i=1}^N E[z_i|x;\pi^{(t)},p_1^{(t)},p_2^{(t)}]}{N} \\
 p_1^{(t+1)} &=& \frac{\sum_{i=1}^N E[z_i|x;\pi^{(t)},p_1^{(t)},p_2^{(t)}] x_i}{\sum_{i=1}^N E[z_i|x;\pi^{(t)},p_1^{(t)},p_2^{(t)}]} \\
 p_2^{(t+1)} &=& \frac{\sum_{i=1}^N (1-E[z_i|x;\pi^{(t)},p_1^{(t)},p_2^{(t)}]) x_i}{\sum_{i=1}^N (1-E[z_i|x;\pi^{(t)},p_1^{(t)},p_2^{(t)}])}
\end{eqnarray}
How are these updates different from the maximum-likelihood expressions derived for the parameters in the setting where all the variables are observed?
\end{enumerate}

\section{K-means (programming) {[}Qirong Ho, 25 points{]}}

The k-means algorithm is one of the most popular unsupervised machine
learning tools, thanks to a combination of simplicity, speed, effectiveness,
and perhaps most importantly, visual intuitiveness. Many of us can
{}``see'' how k-means works, but what is it really doing mathematically?
It turns out that k-means is connected to the following problem: given
an integer $k$ and a set of $d$-dimensional data points $x_{1},\dots,x_{n}$,
pick $k$ $d$-dimensional vectors $c_{1},\dots,c_{k}$ (referred
to as {}``centers'') that minimize the potential function\begin{eqnarray*}
\psi & = & \sum_{i=1}^{n}\min_{j\in\left\{ 1,\dots,k\right\} }\left\Vert x_{i}-c_{j}\right\Vert ^{2}.\end{eqnarray*}
The potential function $\psi$ measures the Euclidean distance between
every point $x_{i}$ and its closest center $c_{j}$. Minimizing $\psi$
essentially involves picking centers $c$ such that every data point
$x$ is close to some center --- a fairly intuitive notion of clustering.
Unfortunately, exactly minimizing $\psi$ happens to be NP-hard (don't
worry if you don't know what that means). We will soon prove that
k-means finds a local minimum of $\psi$ --- in other words, k-means
doesn't necessarily find the best solution to $\psi$, but we hope
it finds a good one at least.

First, recall the k-means algorithm:
\begin{enumerate}
\item Choose initial values for the centers $c_{1},\dots,c_{k}$.
\item For each data point $x_{i}$, let $z_{i}$ be the index of its closest
center, e.g. $z_{i}=j$ means $c_{j}$ is the closest center to $x_{i}$.
\item For each center $c_{j}$, set $c_{j}$ to be the mean of the data
points closer to it than any other center. In other words, set $c_{j}=\frac{\sum_{i=1}^{n}\delta\left(z_{i}=j\right)x_{i}}{\sum_{i=1}^{n}\delta\left(z_{i}=j\right)}$,
where $\delta\left(z_{i}=j\right)=1$ if $z_{i}=j$, and 0 otherwise.
\item Repeat steps 2 and 3 until $\psi$ no longer changes.
\end{enumerate}

\subsection{K-means finds a local minimum of $\psi$}

Let's convince ourselves that the k-means algorithm finds a local
minimum of $\psi$. We shall focus on a single iteration of k-means.
Consider the alternative potential function\begin{eqnarray*}
\phi & = & \sum_{i=1}^{n}\sum_{j=1}^{k}\delta\left(z_{i}=j\right)\left\Vert x_{i}-c_{j}\right\Vert ^{2}.\end{eqnarray*}
Notice that $\phi$ is almost equivalent to $\psi$, except that the
$\min_{j\in\left\{ 1,\dots,k\right\} }$ in $\psi$ has been replaced
by the closest center assignments $z_{i}$. The key insight is that
k-means performs a coordinate-wise minimization of $\phi$: in step
2, we minimize $\phi$ with respect to the closest center assignments
$z_{i}$ (while holding $c_{j}$ fixed), and then in step 3, we minimize
$\phi$ with respect to the centers $c_{j}$ (while holding $z_{i}$
fixed).
\begin{enumerate}
\item {[}\textbf{3 points}{]} Assume the $c_{j}$ are fixed to their values
from the previous iteration. Show that step 2 minimizes $\phi$ with
respect to all closest center assignments $z_{i}$. \textbf{Hint:}
if you don't know where to begin, we're really asking you to prove
the following statement: if you change one or more of the $z_{1},\dots,z_{n}$
produced by step 2, then the value of $\phi$ will either increase
or stay the same. Notice that the double summation in $\phi$ contains
$nk$ terms; you should begin by separating them into $k$ groups,
one for each $z_{i}$.
\item {[}\textbf{5 points}{]} Now assume that the $z_{i}$ are fixed to
their values from step 2. Show that step 3 minimizes $\phi$ with
respect to all centers $c_{j}$. \textbf{Hint:} try differentiating
$\phi$ with respect to a single $c_{j}$.
\end{enumerate}
The alternating minimization on $z_{i}$ and $c_{j}$ guarantees that
k-means finds a local minimum of $\phi$. Moreover, the values of
$c_{j}$ at this local minimum of $\phi$ (which is a function of
$z$ and $c$) also correspond to a local minimum of $\psi$ (which
is a function of $c$ only). We thus conclude that k-means finds a
local minimum of $\psi$.


\subsection{Implementing k-means\label{sub:Implementing-k-means}}

Having convinced ourselves that k-means does something sensible, we
shall now evaluate its performance on the dataset \texttt{faces.csv}. Each of the 2429 lines in this
file corresponds to a $19 \times 19$ image of a human face. Every image is
represented as a 361-dimensional vector of grayscale values, in column-major format (i.e. the first 19 numbers correspond
to the first column, then the next 19 numbers correspond to the second
column, etc.). If you are using MATLAB, you can load the images using
'\texttt{data = csvread('faces.csv')}', and from there you can display
the $i$-th image using '\texttt{imshow(uint8(reshape(data(i,:),{[}19
19{]})))}'.
\begin{enumerate}
\item {[}\textbf{4 points}{]} Implement the k-means algorithm as described
earlier. Note that there are two additional issues we need to consider:
how to choose the initial centers $c_{1},\dots,c_{k}$, and how to
decide if convergence has been reached in step 4. For this part, we
shall pick the initial centers \emph{uniformly at random} (with replacement)
from the set of data points $x_{1},\dots,x_{n}$. As for convergence,
we are not going to use a stopping criterion --- rather, we shall
visualize the performance of k-means by plotting the alternative potential
$\phi$ at every iteration (thus, your code will need to compute $\phi$).
\item {[}\textbf{4 points}{]} Run your implementation 10 times on the \texttt{faces.csv}
dataset, using $k=5$ centers and taking 50 iterations for each run
(one iteration is a single repeat of steps 2-3).
Save the centers learned by each run, you'll need them for the next question! Plot $\phi$ versus iteration
number for all 10 runs on the same graph. What do you observe?
\item {[}\textbf{4 points}{]} From your results in the previous part, find the run
that gives the lowest value of $\phi$ at the last (50th) iteration.
This is your algorithm's best run.
Using your favorite plotting software, draw the $k=5$ centers from
this run as $19\times19$ grayscale images. On MATLAB, the command
'\texttt{imshow(uint8(reshape(c,{[}19 19{]})))}' will plot the 361-dimensional
vector $c$ as a $19\times19$ grayscale image. What do you observe?
Do the 'center-faces' look like faces at all?
\end{enumerate}
Submit all written parts and plots as hardcopies at the beginning of lecture,
and submit all code by email to the grading TA (qho AT andrew DOT cmu
DOT edu). In particular, you are to provide scripts that generate all plots.


\subsection{Choosing good initial centers: the k-means++ algorithm}

In the paper {}``k-means++: The Advantages of Careful Seeding''
(Arthur and Vassilvitskii 2007), the authors describe a probabilistic
initialization scheme for the centers $c_{j}$. This scheme guarantees
that, in expectation, the potential function $\phi$ will take a good
value right from the start. The k-means++ initialization proceeds
as follows:
\begin{enumerate}
\item Choose an initial center $c_{1}$ uniformly at random from the data
points $x_{1},\dots,x_{n}$.
\item Choose the next center $c_{j}$ probabilistically, such that the probability
of choosing $c_{j}=x_{i}$ is $\frac{D\left(x_{i}\right)^{2}}{\sum_{h=1}^{n}D\left(x_{h}\right)^{2}}$,
where $D\left(x\right)$ is the Euclidean distance from $x$ to its
nearest center in $\left\{ c_{1},\dots,c_{j-1}\right\} $.
\item Repeat step 2 until all $k$ centers have been chosen.
\end{enumerate}
We want to compare this initialization scheme with the uniform random
initialization from earlier.
\begin{enumerate}
\item {[}\textbf{5 points}{]} Implement the k-means++ initialization, and
repeat parts 2 and 3 from question \ref{sub:Implementing-k-means}.
Compare your new plots to your old ones. What do you observe?
\end{enumerate}
Again, please email scripts for generating your plots.

\end{document}