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

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

\title{Logistic Regression Using Stochastic Gradient Descent} 

\author{William Cohen}

\begin{document}

\maketitle

\section{Background: SGD for Logistic Regression}

One fairly simple way (and extremely scalable way) to implement logistic regression is stochastic gradient descent.

In the lecture I plan to follow Charles Elkan's
notes\footnote{http://cseweb.ucsd.edu/~elkan/250B/logreg.pdf} which
are for a binary classification task.  We estimate the probability $p$
that an example $\vek{x}=\langle x_1,\ldots,x_d\rangle$ is positive in
the log-odds form:
\begin{equation} \label{eq:logodds}
\log \frac{p}{1-p} = \alpha + \sum_{j=1}{d} \beta_j x_j
\end{equation}
If we assume there is a ``bias feature'' $x_0$ that is true for every
example, then you can simplify and drop the $\alpha$, leaving just the
$\beta_j$'s to estimate.

It's convenient to consider examples of the form $\vek{x},y$ where
$y=0$ or $y=1$.  The log of the conditional likelihood for example is
example will be $LCL(\vek{x},y)=\log p$ if $y=1$ and
$LCL(\vek{x},y)=\log (1-p)$ if $y=0$, where $p$ is computed as in
Eq.~\ref{eq:logodds}.  With a little calculus you can show that for a
positive example,
\[
\frac{\partial}{\partial \beta_j} LCL(\vek{x},y) = \frac{1}{p} \frac{\partial}{\partial \beta_j}p
\]
and for a negative example,
\[
\frac{\partial}{\partial \beta_j} LCL(\vek{x},y) = \frac{1}{1-p} ( - \frac{\partial}{\partial \beta_j}p )
\]
and that 
\[
\frac{\partial}{\partial \beta_j}p = p (1-p) x_j
\]
and putting this together we get that if $y=1$
\[
\frac{\partial}{\partial \beta_j} LCL(\vek{x},y) = (1-p) x_j
\]
and if $y=0$ then
\[
\frac{\partial}{\partial \beta_j} LCL(\vek{x},y) = -p x_j
\]
so in either case
\begin{equation} \label{eq:lcl}
\frac{\partial}{\partial \beta_j} LCL(\vek{x},y) = (y-p) x_j
\end{equation}
So an update to the $\beta$'s that would improve most LCL would be
along the gradient---i.e., for some small step size $\lambda$,
let 
\[ \beta_j = \beta_j + \lambda (y-p) x_i
\]
Notice that if $x_i=0$ then $\beta_j$ is unchanged.

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 $B$
\item For $t=1,\ldots,T$
  \begin{itemize}
  \item For each example $\vek{x}_i,y_i$:
     \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 $B$, set $B[j]=0$.
	\item Set $B[j] = B[j] + \lambda (y-p) x_i$
        \end{itemize}
     \end{itemize}
  \end{itemize}
\item Output the parameters $\beta_1,\ldots,\beta_d$.
\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{Efficient regularized SGD}

Logistic regression tends to overfit when there are many rare
features.  One fix is to penalize large values of $\beta$, by
optimizing, instead of $LCL$, some function such as
\( LCL - \mu \sum_{j=1}^d \beta_j^2
\)
Here $\mu$ controls how much weight to give to the penalty term.
The update for $\beta_j$ becomes
\[ \beta_j = \beta_j + \lambda ( (y-p) x_i - 2\mu\beta_j  )
\]
or equivalently
\[ \beta_j = \beta_j + \lambda (y-p) x_i - \lambda 2\mu\beta_j
\]
Experimentally this greatly improves overfitting - but unfortunately,
this makes the computation much more expensive, because now {\em every}
$\beta_j$ needs to be updated, not only the ones that are non-zero.

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_i$.  Let's
call this the ``LCL'' part of the update.  The second is the
``regularization part'' of the update, which is to replace $\beta$ by
\[ \beta_j = \beta_j - \lambda 2\mu\beta_j = \beta_j \cdot (1-2\lambda\mu)
\]
So we could perform our update of $\beta_j$ as follows:
\begin{itemize}
\item Set $\beta_j = \beta_j \cdot (1-2\lambda\mu)$
\item If $x_j\not=0$, set $\beta_j = \beta_j + \lambda (y-p) x_i$
\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 $\beta_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 $B$ be empty hashtables.  $A$ will record the
value of $k$ last time $B[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 $B$, set $B[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 
          \[ B[j] = B[j] \cdot (1-2\lambda\mu)^{k - A[j]}
          \]
	\item Set $B[j] = B[j] + \lambda (y-p) x_i$
        \item Set $A[j] = k$
        \end{itemize}
     \end{itemize}
  \end{itemize}
\item For each parameter $\beta_1,\ldots,\beta_d$, set
\[ B[j] = B[j] \cdot (1-2\lambda\mu)^{k - A[j]}
\]
\item Output the parameters $\beta_1,\ldots,\beta_d$.
\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}
