\documentclass[12pt]{article}  
\usepackage{epsfig}

\newcommand{\vek}[1]{{\bf {#1}}}

\title{Stochastic Gradient Descent: Lecture Notes} 

\author{William W. Cohen}

\begin{document}

\maketitle

\section{A Warmup}

If you have a continuous mathematical expression for a probability,
you can often learn the parameter values by simply maximizing the
probability of the data with respect to the parameter values.

For example, suppose we have a binomial---a coin with some unknown
probability $\theta$ of heads (i.e., $x_i=1$)---and $n$ draws from it,
$x_1,\ldots,x_n$, of which $k$ came up heads.  We can write the
probability of a single $x_i$ as
\[ P(X=x_i|\theta) = \left\{ \begin{array}{ll} \theta & \mbox{if $x_i=1$} \\ 1-\theta & \mbox{if $x_i=0$} \end{array} \right.
\]
or better, since we'd like something we can differentiate (to find the max), 
\[
 P(X=x_i|\theta) = \theta^{x_i}(1-\theta)^{1-x_i}
\]
Let $D$ be the dataset, $D=\{ x_1,\ldots,x_n \}$.  The
\[ P(D|\theta)= \prod_i P(x_i|\theta) = \prod_i \theta^{x_i}(1-\theta)^{1-x_i} = \theta^{k}(1-\theta)^{n-k} 
\] 
So we can differentiate, using $(fg)'=f'g+fg'$ and the chain rule for $g$....
\[
\frac{\partial}{\partial\theta} P(D) = k\theta^{k-1}(1-\theta)^{n-k} + \theta^k(n-k)(1-\theta)^{n-k-1}(-1)
\]
Now simplify by pulling out the common factors of $\theta^{k-1}$ and $(1-\theta)^{n-k-1}$ to get
\[
\frac{\partial}{\partial\theta} P(D) = \theta^{k-1}(1-\theta)^{n-k-1}\left(k(1-\theta)-\theta(n-k)\right)
\]
Setting this to zero, we need one of the factors to be zero, so either
$\theta^{k-1}=0$ at $\theta=0$, or $(1-\theta)^{n-k-1}=0$ at
$\theta=1$, or $k(1-\theta)-\theta(n-k)=0$, which we can solve to give
$\theta=n/k$, the maximum.

\section{Logistic regression}

\newcommand{\ee}{\frac{1}{1+e^{-\vek{x}\cdot\vek{w}}}}
\newcommand{\eee}{\frac{1}{1+\exp(-\sum_j x^j w^j)}}
\newcommand{\pwj}{\frac{\partial}{\partial w^j}}

The model we'd like to optimize next is for an example $(\vek{x},y)$,
where $y=0$ or $y=1$, and $\vek{x}$ is a feature vector.  We will let
the predicted value of $y$ be
\[ P(Y=y|X=\vek{x},\vek{w}) = \ee
\]
So then the probability of a single example $(\vek{x},y)$ is
\[ P(Y=y|X=\vek{x},\vek{w}) = 
   \left\{ \begin{array}{ll} \ee & \mbox{if $y=1$} \\ 1-\ee & \mbox{if $y=0$} \end{array} \right.
\]
To motivate this formula, note that this similar in form to the other
linear classifiers we've looked at, like Naive Bayes, Rocchio, and the
perceptron---it's just that instead instead of using an argmax or a
sign function on the inner product $\vek{x}\cdot\vek{w}$, we're using
the easy-to-differentiate logistic function $f(z)=\frac{1}{1+e^{-z}}$,
which approximates a step function between 0 and 1.

What we'll do is look at the gradient of $\log P_\vek{w}(y|\vek{x})$
for a single example $(\vek{x},y)$.  The algorithm we will use will be
to pick a random example, compute this gradient, and adjust the
parameters to take a small step in that direction - so we're using a
randomized approximation to the true gradient over all the data.  This
turns out to be fast and surprisingly effective---and algorithmically
a lot like the perceptron method.

As notation, let $p$ be
\[ p \equiv \ee = \eee 
\]
and consider the log probability (which, since it's monotonic with the
probability, will have the same maximal values):
\[ \log P(Y=y|X=\vek{x},\vek{w}) = 
   \left\{ \begin{array}{ll} \log p & \mbox{if $y=1$} \\ \log(1-p) & \mbox{if $y=0$} \end{array} \right.
\]
Let's take the gradient with respect to some parameter value $w^j$.
We'll do the two cases separately for now: using $(\log{}f)'=\frac{1}{f}f'$, we have
\[ \pwj \log P(Y=y|X=\vek{x},\vek{w}) = 
   \left\{ 
       \begin{array}{ll} 
           \frac{1}{p} \pwj p  & \mbox{if $y=1$} \\ 
           \frac{1}{1-p} (-\pwj p) & \mbox{if $y=0$} 
       \end{array} \right.
\]
Now we can churn through and get the gradient for $p$, using $(e^{f})'=e^f f'$ and the chain rule:
\begin{eqnarray*}
\pwj p & = & \pwj (1 + \exp(-\sum_j x^j w^j))^{-1} \\
       & = & (-1) (1 + \exp(-\sum_j x^j w^j))^{-2} \pwj \exp(-\sum_j x^j w^j) \\
       & = & (-1) (1 + \exp(-\sum_j x^j w^j))^{-2} \exp(-\sum_j x^j w^j) (-x^j)\\
       & = & \frac{1}{1 + \exp(-\sum_j x^j w^j)} \frac{\exp(-\sum_j x^j w^j)}{1 + \exp(-\sum_j x^j w^j)} x^j
\end{eqnarray*}
Note that 
\[ 1-p = \frac{1 + \exp(-\sum_j x^j w^j)}{1 + \exp(-\sum_j x^j w^j)} - \frac{1}{1 + \exp(-\sum_j x^j w^j)} 
  = \frac{\exp(-\sum_j x^j w^j)}{1 + \exp(-\sum_j x^j w^j)}
\]
so the final line can be rewritten to give
\[
\pwj p  = p (1-p) x^j
\]
which we can plug into the cases above and get
\[ \pwj \log P(Y=y|X=\vek{x},\vek{w}) = 
   \left\{ 
       \begin{array}{ll} 
           \frac{1}{p} p (1-p) x^j  = (1-p)x^j & \mbox{if $y=1$} \\ 
           \frac{1}{1-p} (-1) p (1-p) x^j = -p x^j  & \mbox{if $y=0$} 
       \end{array} \right.
\]
Finally, these cases can be combined to give the very simple gradient
\[
\pwj \log P(Y=y|X=\vek{x},\vek{w}) = (y-p) x^j
\]
So, taking a small step in this direction would be to increment $\vek{w}$ as follows:
\[
 \vek{w}^{(t+1)} = \vek{w}^{(t)} + \lambda (y-p) \vek{x} 
\]
Compare this to the perceptron update rule: it's not very different.

So this leads to this algorithm, which is very fast (assuming you have enough memory to hash all the parameter values).
\begin{enumerate}
\item Initialize a hashtable $W$
\item For $t=1,\ldots,T$
  \begin{itemize}
  \item For each example $\vek{x}_i,y_i$:
     \begin{itemize}
     \item Compute the prediction for $\vek{x}_i$:
       \[
       p_i = \frac{1}{1+\exp(-\sum_{j:x_i^j>0} x_i^j w^j)}
       \]
     \item For each non-zero feature of $\vek{x_i}$ with index $j$ and value $x^j$:
        \begin{itemize}
	\item If $j$ is not in $W$, set $W[j]=0$.
	\item Set $W[j] = W[j] + \lambda (y_i-p_i) x^j$
        \end{itemize}
     \end{itemize}
  \end{itemize}
\item Output the hash table $W$.
\end{enumerate}

The time to run this is $O(nT)$, where $n$ is the total number of
non-zero features for each example and $T$ is the number of
iterations.

\section{An observation}

Consider averaging the gradient
\[
\pwj \log P(Y=y|X=\vek{x},\vek{w}) = (y-p) x^j
\]
over many examples $D=\{(\vek{x}_1,y_1),\ldots,(\vek{x}_n,y_n)\}$---which is the usual non-stochastic gradient procedure---and assume
that the $x^j$'s are binary.  If the examples are independent, then
\[
\pwj \log P(D|\vek{w}) = \frac{1}{n} \sum_i (y_i-p_i) x_i^j = \frac{1}{n} \sum_{i:x_i^j=1} y_i - \frac{1}{n} \sum_{i:x_i^j=1} p_i 
\]
This can be viewed as the difference between the actual fraction of
labels $y=1$ in the subset of the data for which $x^j=1$, and the
predicted fraction of labels $y=1$ in the subset of the data for which
$x^j=1$.

\section{Efficient regularized SGD}

Logistic regression tends to overfit when there are many rare
features.  One fix is to penalize large values of $w^j$, by
optimizing, instead of LCL (log-conditional likelihood), some function such as
\[ LCL - \mu \sum_{j=1}^d (w^j)^2
\]
Here $\mu$ controls how much weight to give to the penalty term.
Then instead of
\[
\pwj \log P(Y=y|X=\vek{x},\vek{w}) = (y-p) x^j
\]
we have
\[
\pwj \log P(Y=y|X=\vek{x},\vek{w}) - \mu \sum_{j=1}^d (w^j)^2 = (y-p) x^j - 2\mu w^j
\]
and the update for $w^j$ becomes
\[ w^j = w^j + \lambda ( (y-p) x^j - 2\mu w^j  )
\]
or equivalently
\[ w^j = w^j + \lambda (y-p) x^j - \lambda 2\mu w^j
\]
Experimentally this greatly improves overfitting - but unfortunately,
this makes the computation much more expensive, because now {\em every}
$w^j$ needs to be updated, not only the ones that are associated with non-zero $x^j$'s.

The trick to making this efficient is to break the update into two
parts.  One is the usual update of adding $\lambda (y-p) x^j$.  Let's
call this the ``LCL'' part of the update.  The second is the
``regularization part'' of the update, which is to replace $w^j$ by
\[ w^j = w^j - \lambda 2\mu w^j = w^j \cdot (1-2\lambda\mu)
\]
So we could perform our update of $w^j$ as follows:
\begin{itemize}
\item Set $w^j = w^j \cdot (1-2\lambda\mu)$
\item If $x_j\not=0$, set $w^j = w^j + \lambda (y-p) x^j$
\end{itemize}

Following this up, we note that we can perform $m$ successive
``regularization'' updates by letting 
\( B_j = B_j\cdot(1-2\lambda\mu)^m
\). The basic idea of the new algorithm is to not
perform regularization updates for zero-valued $x_j$'s, but instead to
simply keep track of how many such updates would need to be performed
to update $w^j$, and perform them only when we would normally
perform ``LCL'' updates (or when we output the parameters at the end
of the day).

Here's the final algorithm (for more detail, see ``Lazy sparse
stochastic gradient descent for regularized multinomial logistic
regression'', Bob Carpenter,
http://lingpipe.files.wordpress.com/2008/04/lazysgdregression.pdf.)


\begin{enumerate}
\item Let $k=0$, and let $A$ and $W$ be empty hashtables.  $A$ will record the
value of $k$ last time $W[j]$ was updated.
\item For $t=1,\ldots,T$
  \begin{itemize}
  \item For each example $\vek{x}_i,y_i$:
    \item Let $k=k+1$
     \begin{itemize}
     \item For each non-zero feature of $\vek{x_i}$ with index $j$ and value $x^j$:
        \begin{itemize}
	\item If $j$ is not in $W$, set $W[j]=0$.
	\item If $j$ is not in $A$, set $A[j]=0$.
        \item Simulate the ``regularization'' updates that would have been performed 
          for the $k-A[j]$ examples since the last time a non-zero $x_j$ was encountered
          by setting 
          \[ W[j] = W[j] \cdot (1-2\lambda\mu)^{k - A[j]}
          \]
	\item Set $W[j] = W[j] + \lambda (y-p) x^j$
        \item Set $A[j] = k$
        \end{itemize}
     \end{itemize}
  \end{itemize}
\item For each parameter $w^1,\ldots,w^j,\ldots,w^d$, set
\[ W[j] = W[j] \cdot (1-2\lambda\mu)^{k - A[j]}
\]
\item Output the hash table $W$.
\end{enumerate}

\section{Other hints}

The learning rate $\lambda$ is often decreased over time.  On the
$t$-th sweep through the data, set $\lambda=\frac{\eta}{t^2}$.  I used
$\eta=0.5$.  (Sometimes $\lambda$ is also scaled by $1/n_e$, where
$n_e$ is the number of examples.)  I also used a value of $2\mu=0.1$.

When running stochastic gradient descent, it usual to randomize the
order of examples, and scale the feature values so that they are
comparable (if they are not already binary).  However, randomization
is not trivial to do for a large dataset.  I recommend implementing
SGD as a process that streams once through a data stream, with the
number of examples $n_e$ being passed in separately as a command-line
argument so that the algorithm is aware of what the current value of
$t$ is.  Then write a separate module that will input a file of
examples and then stream the individual examples out in approximately
random order.


\end{document}
