\documentclass[fleqn]{article}
\usepackage{graphicx}
\usepackage{latexsym}
\usepackage{html}
%\begin{latexonly}
\usepackage{bm}
\usepackage{mathrsfs}
%\end{latexonly}
\usepackage{natbib}
\bibpunct{(}{)}{;}{a}{,}{,}         % (Cheeseman, 1992)

\newcommand{\Old}[1]{#1^{\scriptstyle\rm old}}
\newcommand{\New}[1]{#1^{\scriptstyle\rm new}}

\newcommand{\maximize}{\mathop{\rm maximize}}
\newcommand{\minimize}{\mathop{\rm minimize}}
\newcommand{\bw}{{\mathbf w}}
\newcommand{\bx}{{\mathbf x}}
\newcommand{\ba}{{\mathbf a}}
\newcommand{\bzero}{{\mathbf 0}}
\begin{htmlonly}
\newcommand{\bm}[1]{#1}
\end{htmlonly}
\newcommand{\balpha}{\bm{\alpha}}
\newcommand{\bxi}{\bm{\xi}}
\newcommand{\bmu}{\bm{\beta}}
\newcommand{\Lag}{\mathscr{L}}
\newcommand{\Part}[2]{\frac{\partial {#1}}{\partial {#2}}}

%\begin{latexonly}
\newcommand{\CPLUSPLUS}{C\nolinebreak\hspace{-.05em}\raisebox{.4ex}{\tiny\bf +}\nolinebreak\hspace{-.10em}\raisebox{.4ex}{\tiny\bf +}}
%\end{latexonly}
\begin{htmlonly}
\newcommand{\CPLUSPLUS}{C++}
\end{htmlonly}
%\def\CPLUSPLUS{{C\nolinebreak[4]\hspace{-.05em}\raisebox{.4ex}{\tiny\bf ++}}}

%\font\ninerm=cmr9
%\font\sevenrm=cmr7
%\let\mc=\ninerm % medium caps
%\def\CPLUSPLUS/{{\mc C\PP\spacefactor1000}}
%\newbox\PPbox % symbol for ++
%\setbox\PPbox=\hbox{\kern.5pt\raise1pt\hbox{\sevenrm+\kern-1pt+}\kern.5pt}
%\def\PP{\copy\PPbox}
%\def\CEE/{{\mc C\spacefactor1000}}
%\def\SP{{\tt\char`\ }} % (visible) space in a string

\begin{document}
\title{Sequential Minimal Optimization for SVM}
\date{}
\maketitle

\tableofcontents
\newpage
\begin{abstract}
This is a \CPLUSPLUS\ implementation of John C. Platt's sequential minimal
optimization (SMO) for training a support vector machine (SVM).  This
program is based on the pseudocode in \citet{Platt:1998:SMO-BOOK}.

This is both the documentation and the \CPLUSPLUS\ code. It is a
\verb|NUWEB| document from which both the \LaTeX\ file and the \CPLUSPLUS\ file
can be generated. The documentation is essentially my notes when reading
the papers (most of them being {\it cut-and-paste} from the papers).
\end{abstract}


\section{Introduction to Support Vector Machine (SVM)}

This introductio to Support Vector Machine for
binary classification is based on \citet{Burges:1998}.

\subsection{Linear SVM}
First let us look at the linear support vector machine. It is based on
the idea of hyperplane classifier, or linearly separability.

\begin{figure}
\begin{center}
%\leavevmode
%\input{pic/fig1}
%    \centerline{\raise 1em\box\graph}
\includegraphics{pic/fig1}
\end{center}
\end{figure}

Suppose we have $N$ training data points $\{(\bx_1,y_1),
(\bx_2,y_2), \ldots, (\bx_N,y_N)\}$
where $ \bx_i \in {\cal R}^d$ and $ y_i \in \{\pm 1\}$.
We would like to learn a linear separating hyperplane classifier:
\[
   f(\bx) = sgn(\bw \cdot \bx - b).
\]

Furthermore, we want this hyperplane to have the maximum separating margin
with respect to the two classes.
Specifically, we want to find this hyperplane 
$H:\ y = \bw \cdot \bx - b = 0$
and two hyperplanes parallel to it and with equal distances to it,
\[
H_1:\ y = \bw \cdot \bx - b = +1,
\]
\[
H_2:\ y = \bw \cdot \bx - b = -1,
\]
with the condition that there are no data points between $H_1$ and $H_2$,
and the distance between $H_1$ and $H_2$ is maximized.

For any separating plane $H$ and the corresponding $H_1$ and $H_2$, we
can always ``normalize" the coefficients vector $\bw$ so that
$H_1$ will be $y = \bw \cdot \bx - b = +1$, and
$H_2$ will be $y = \bw \cdot \bx - b = -1$. See Appendix \ref{appendix:1}
for details.

We want to maximize the distance between $H_1$ and $H_2$. So there will be some
positive examples on $H_1$ and some negative examples on $H_2$. These examples
are called {\it support vectors} because only they participate in the definition
of the separating hyperplane, and other examples can be removed and/or moved
around as long as they do not cross the planes $H_1$ and $H_2$.

Recall that in 2-D, the distance from a point $(x_0,y_0)$ to a line
$Ax+By+C=0$ is $\frac{\mid Ax_0 + By_0 + C \mid }{\sqrt{A^2+B^2}}$.
Similarly, the distance of a point on $H_1$ to 
$H:\ \bw \cdot \bx - b = 0$ is 
$\frac{\mid \bw \cdot \bx -b
\mid}{\| \bw \|}
= \frac{1}{\| \bw \|}$, and the distance between $H_1$ and $H_2$
is $\frac{2}{\| \bw \|}$.  So, in order to maximize
the distance, we should minimize $\| \bw \| = \bw^T \bw$ with the condition
that there are no data points between $H_1$ and $H_2$:
\[
\bw \cdot \bx - b \geq +1,\quad
\mbox{for positive examples $y_i = +1$},
\]
\[
\bw \cdot \bx - b \leq -1,\quad
\mbox{for negative examples $y_i = -1$}.
\]
These two conditions can be combined into
\[
y_i ( \bw \cdot \bx_i - b) \geq 1.
\]

So our problem can be formulated as
\[
\min_{\bw,b} \frac{1}{2} \bw^T \bw \quad 
\mbox{ subject to $ y_i ( \bw \cdot \bx_i - b) \geq 1$}.
\]
This is a convex, quadratic programming problem (in $\bw,b$),
in a convex set.

Introducing Lagrange multipliers $\alpha_1,\alpha_2,\ldots,\alpha_N \geq 0$,
we have the following Lagrangian:
\[
\Lag(\bw,b,{\balpha}) \equiv \frac{1}{2} \bw^T\bw 
- \sum_{i=1}^{N} \alpha_i y_i ( \bw \cdot \bx_i - b) + \sum_{i=1}^N
  \alpha_i.
\]

\subsection{The dual problem}

We can solve the Wolfe dual instead:  {\it maximize}
$\Lag(\bw,b,\balpha)$ with respect to $\balpha$, 
subject to the constraints that the
gradient of $\Lag(\bw,b,\balpha)$ with respect to the primal variables 
$\bw$ and $b$ vanish:
\begin{equation}
\Part{\Lag}{\bw} = \bzero,
\label{eq:gradient-w}
\end{equation}
\begin{equation}
\Part{\Lag}{b} = 0,
\label{eq:gradient-b}
\end{equation}
and that
\[
\balpha \geq \bzero.
\]
From Equations \ref{eq:gradient-w} and \ref{eq:gradient-b}, we have
\[
\bw = \sum_{i=1}^N \alpha_i y_i \bx_i,
\]
\[
 \sum_{i=1}^N \alpha_i y_i = 0.
\]
Substitute them into $\Lag(\bw,b,\balpha)$, we have
\[
\Lag_D \equiv \sum_{i=1}^N \alpha_i - \frac{1}{2} \sum_{i,j}\alpha_i \alpha_j y_i y_j
\bx_i \cdot \bx_j,
\]
in which the primal variables are eliminated.

When we solve $\alpha_i$, we can get $\bw = \sum_{i=1}^N \alpha_i
y_i \bx_i$, (we will later show how to compute the threshold $b$), and
we can classify a new object $\bx$ with
\begin{eqnarray*}
f(\bx) & = & sgn(\bw \cdot \bx + b) \\
 & = & sgn((\sum_{i=1}^N \alpha_i y_i \bx_i) \cdot \bx + b) \\ 
 & = & sgn(\sum_{i=1}^N \alpha_i y_i (\bx_i \cdot  \bx) + b).
\end{eqnarray*}

Please note that in the objective function and the solution, the
training vectors $\bx_i$ occur only in the form of dot product.

Before going into the details to how to solve this quadratic
programming problem, let's extend it in two directions.

\subsection{Non-linear SVM}
What if the surface separating the two classes are not linear?
Well, we
can transform the data points to another high dimensional space such
that the data points will be linearly separable. Let the transformation
be $\Phi(\cdot)$. In the high dimensional space, we solve

\[
\Lag_D \equiv \sum_{i=1}^N \alpha_i - \frac{1}{2} \sum_{i,j}\alpha_i \alpha_j y_i y_j
\Phi(\bx_i) \cdot \Phi(\bx_j)
\]

Suppose, in addition, $\Phi(\bx_i) \cdot \Phi(\bx_j) =
k(\bx_i, \bx_j)$.  That is, the dot product in that high
dimensional space is equivalent to a {\it kernel} function of the
input space. So we need not be explicit about the transformation
$\Phi(\cdot)$ as long as we know that the kernel function 
$k( \bx_i, \bx_j)$ is equivalent to the dot product of some
other high dimensional space. There are many kernel functions that can
be used this way, for example, the radial basis function (Gaussian kernel)
\[
K(\bx_i, \bx_j) = e ^ {- \|  \bx_i - \bx_j \|^2 / 2 \sigma^2}.
\]

The Mercer's condition can be used to determine if a function can be used as a
kernel function:
\begin{quote}
There exists a mapping $\Phi$ and an expansion
\[
K(\bx,\mathbf{y}) = \sum_i \Phi(\bx)_i \Phi(\mathbf{y})_i,
\]
if and only if, for any $g(\bx)$ such that
\[
\int g(\bx)^2 d \bx \ \mbox{is finite},
\]
then
\[
\int K(\bx,\mathbf{y})g(\bx)g(\mathbf{y})d\bx d\mathbf{y} \geq 0.
\]
\end{quote}

\subsection{Imperfect separation}
The other direction to extend SVM is to allow for noise,
or imperfect separation. That is, we do not strictly enforce
that there be no data points between $H_1$ and $H_2$,
but we definitely want to penalize the data points that cross
the boundaries. The penalty $C$ will be finite. (If $C=\infty$,
we come back to the original perfect separating case.)

We introduce non-negative slack variables $\xi_i \geq 0$, so that
\[
\bw \cdot \bx_i - b \geq +1 - \xi_i \quad \mbox{for $y_i=+1$},
\]
\[
\bw \cdot \bx_i - b \leq -1 + \xi_i \quad \mbox{for $y_i=-1$},
\]
\[
\xi_i \geq 0, \quad \forall i.
\]
and we add to the objective function a penalizing term:
\[
\minimize_{\bw,b,\bxi} \frac{1}{2} \bw^T\bw + C(\sum_i \xi_i)^m
\]
where $m$ is usually set to 1, which gives us
\begin{eqnarray*}
    \minimize_{\bw,b,\xi_i} & \frac{1}{2}\bw^T \bw + C \sum_{i=1}^N \xi_i & \\
    \mbox{subject to } &  y_i (\bw^T \bx_i - b) + \xi_i - 1 \geq 0, & 
         1\leq i \leq N \\
                       & \xi_i \geq 0, & 1\leq i \leq N
\end{eqnarray*}
Introducing Lagrange multipliers $\balpha$, $\bmu$, the Lagrangian is
\begin{eqnarray*}
\Lag(\bw,b,\xi_i; \balpha, \bmu)  & = &
 \frac{1}{2}\bw^T \bw + C \sum_{i=1}^N \xi_i \\
& &
- \sum_{i=1}^N \alpha_i \bigl[ y_i (\bw^T \bx_i - b) + \xi_i - 1 \bigr]
- \sum_{i=1}^N \mu_i \xi_i \\
& = &  \frac{1}{2}\bw^T \bw 
+ \sum_{i=1}^N \left( C - \alpha_i - \mu_i \right) \xi_i \\
& &  
- \left( \sum_{i=1}^N \alpha_i y_i \bx_i^T \right) \bw
- \left( \sum_{i=1}^N \alpha_i y_i \right) b + \sum_{i=1}^N \alpha_i
\end{eqnarray*}

Neither the $\xi_i$'s, nor their Lagrange multipliers, appear in the
Wolfe dual problem:
\[
\maximize_{\balpha} \Lag_D \equiv \sum_i \alpha_i - \frac{1}{2} \sum_{i,j} \alpha_i \alpha_j y_i
y_j \bx_i \cdot \bx_j
\]
subject to:
\[
0 \leq \alpha_i \leq C,
\]
\[
\sum_i \alpha_i y_i = 0.
\]
The only difference from the perfectly separating case is that $\alpha_i$ is
now bounded above by $C$ instead of $\infty$.
For details, see Appendiex \ref{appendix:2}.

The solution is again given by
\[
\bw = \sum_{i=1}^N \alpha_i y_i \bx_i
\]

To train the SVM, we search through the feasible region of the dual problem
and maximize the objective function. The optimal solution can be checked
using the KKT conditions.

\subsection{The KKT conditions}
The KKT optimality conditions of the primal problem are, the
gradient of $\Lag(\bw,b,\balpha,\bmu)$ with respect to the primal variables 
$\bw$, $b$, $\bxi$ vanishes (this must always be satisfied by the dual
problem), and that
for $1 \leq i \leq N$,
\begin{equation}
\alpha_i (y_i (\bw^T \bx_i - b) + \xi_i - 1) = 0,
\label{eq:kkt:1}
\end{equation}
\begin{equation}
\mu_i \xi_i = 0.
\label{eq:kkt:2}
\end{equation}
Depending on the value of $\alpha_i$, we have three cases to consier:
\begin{enumerate}
\item
    If $\alpha_i = 0$, then $\mu_i = C - \alpha_i  = C > 0$. 
    From Equation \ref{eq:kkt:2}, $\xi_i = 0$. so we have
    \[
    y_i (\bw^T \bx_i - b) - 1 \geq 0.
    \]
\item
    If $0 < \alpha_i < C$, from Equation \ref{eq:kkt:1}, we have 
    \begin{equation}
    y_i (\bw^T \bx_i - b) + \xi_i - 1=0
    \label{eq:kkt:boundary}
    \end{equation}
    Note that $\mu_i = C - \alpha_i  > 0$, so
    $\xi_i = 0$ (Equation \ref{eq:kkt:2}). Substituting
    into Equation \ref{eq:kkt:boundary}, we have
    \[
    y_i(\bw^T \bx_i - b) - 1 = 0.
    \]
\item If $\alpha_i = C$, then from Equation \ref{eq:kkt:1}, we have
    \begin{equation}
    y_i (\bw^T \bx_i - b) + \xi_i - 1=0
    \label{eq:kkt:crossing}
    \end{equation}
    Note that $\mu_i = C - \alpha_i  = 0$, we have $\xi_i \geq 0$. So
    \[
    y_i (\bw^T \bx_i - b) - 1 \leq 0.
    \]
\end{enumerate}

The quantity $y_i (\bw^T \bx_i - b) - 1$ can be computed as
\[
R_i = y_i (\bw^T \bx_i - b) - y_i^2
= y_i(\bw^T \bx_i - b - y_i) = y_i E_i
\]
where $E_i=\bw^T\bx_i - b -y_i = u_i - y_i$ is the prediction error.

To summarize, the KKT condition implies:
\begin{eqnarray*}
  \alpha_i =0  & \Rightarrow & R_i \geq 0, \\
  0 < \alpha_i < C  & \Rightarrow & R_i \approx 0, \\
  \alpha_i = C  & \Rightarrow & R_i \leq 0.
\end{eqnarray*}

In the following two cases, the KKT condition is violated:
\begin{itemize}
 \item $\alpha_i < C$ and $R_i < 0$,
 \item $\alpha_i > 0$ and $R_i > 0$.
\end{itemize}

\subsection{Checking KKT condition without using threshold $b$}
As the dual problem does not solve for the threshold $b$ directly, it
would be beneficial to check the KKT condition without using threshold $b$.
This technique is due to \citet{Keerthi+Shevade+Bhattacharyya+Murthy:2001}.

The quantity $ y_i ( \bw^T \bx_i - b) -1$ (which must $\geq 0$
for all $i$ if the KKT condition is satisfied)
can also be written as
\begin{eqnarray*}
& & y_i( \bw^T \bx_i - b) - 1 \\
& =  & y_i( \bw^T \bx_i - b) - y_i^2 \\
& - & y_i(\bw^T \bx_i - y_i - b) \\
& = & y_i(F_i - b),
\end{eqnarray*}
where $F_i \equiv \bw^T \bx_i - y_i$.

Note for $E_i = F_i - b$, we have $E_i - E_j = F_i - F_j$.
(This equality is useful, as Platt's SMO algorithm uses $E_i-E_j$
when optimization the two Lagrange multipliers $\alpha_i$, $\alpha_j$.)

This notation is useful because the KKT conditions
\begin{eqnarray*}
  \alpha_i =0  & \Rightarrow & y_i(F_i - b) \geq 0 \\
  0 < \alpha_i < C  & \Rightarrow & y_i(F_i-b) \approx 0 \\
  \alpha_i = C  & \Rightarrow & y_i(F_i-b) \leq 0
\end{eqnarray*}
can be written as
\begin{eqnarray*}
  i\in I_0\cup I_1 \cup I_2  & \Rightarrow & F_i \geq b \\
  i\in I_0\cup I_3 \cup I_4  & \Rightarrow & F_i \leq b,
\end{eqnarray*}
where
\begin{eqnarray*}
  I_0 & \equiv & \{i: 0 < \alpha_i < C\} \\
  I_1 & \equiv & \{i: y_i = +1, \alpha_i = 0\} \\
  I_2 & \equiv & \{i: y_i = -1, \alpha_i = C\} \\
  I_3 & \equiv & \{i: y_i = +1, \alpha_i = C\} \\
  I_4 & \equiv & \{i: y_i = -1, \alpha_i = 0\}.
\end{eqnarray*}
So that $\forall i \in I_0\cup I_1 \cup I_2$, 
and $\forall j\in I_0\cup I_3 \cup I_4$, we should have $F_i \geq F_j$,
if KKT condition is satisfied.

\begin{center}
\includegraphics[width=3in,height=1.5in]{pic/I0-I4}
\end{center}

To check if this condition holds, we define
\begin{eqnarray*}
b_{\mbox{up}} &=& \min \{F_i: i \in I_0 \cup I_1 \cup I_2\},\\
b_{\mbox{low}} &=& \max \{F_i: i \in I_0 \cup I_3 \cup I_4\}.
\end{eqnarray*}
The KKT condition implies $b_{\mbox{up}} \geq b_{\mbox{low}}$, and
similarly, $\forall i \in I_0\cup I_1 \cup I_2$, $F_i \geq b_{\mbox{low}}$,
and $\forall i \in I_0\cup I_3 \cup I_4$, $F_i \leq b_{\mbox{up}}$.

These comparisons do not use the threshold $b$.

As an added benefit, given the first $\alpha_i$, 
these comparisons automatically finds the second
$\alpha_i$ for joint optimization in SMO.

\section{SMO Algorithm}

\subsection{Optimize two $\alpha_i$'s}

The SMO algorithm searches through the feasible region
of the dual problem and maxmizes the objective function
\[
\Lag_D \equiv \sum_{i=1}^N \alpha_i - \frac{1}{2} 
       \sum_{i,j}\alpha_i \alpha_j y_i y_j
\bx_i \cdot \bx_j,
\]
\[
0 \leq \alpha_i \leq C, \quad \forall i.
\]
It works by optimzing two $\alpha_i$'s at a time (with the other
$\alpha_i$'s fixed). 
It uses heuristics to choose the two $\alpha_i$' for optimization. This is essentially a hill-climbing.

Without loss of generality, suppose we are optimizing $\alpha_1$, $\alpha_2$,
from an old set of feasible solution: $\Old{\alpha}_1$, $\Old{\alpha}_2$, $\alpha_3$, $\ldots$, $\alpha_N$. (For initialization, we can set
$\Old{\balpha}=\bzero$.)

Because $\sum_{i=1}^N y_i \alpha_i = 0$, we have 
\[
y_1 \alpha_1+y_2 \alpha_2
= y_1 \Old{\alpha}_1+y_2 \Old{\alpha}_2.
\]
This confines the optimization to be on a line,
as shown in the following figure:
\begin{center}
%\unsetbox\graph
%\leavevmode
%\input{pic/fig2}
    %\centerline{\raise 1em\box\graph}
\includegraphics{pic/fig2}
\end{center}
\label{page:L-H-graph}

Let $s=y_1 y_2$. Multiply 
\[
y_1 \alpha_1+y_2 \alpha_2 = \mbox{Const.}
\]
by $y_1$, and we have
\[
\alpha_1 = \gamma - s \alpha_2.
\]
where $\gamma \equiv \alpha_1 + s\alpha_2 = \Old{\alpha_1} + s \Old{\alpha_2}$.

Fixing the other $\alpha_i$'s, the objective function can be written as
\begin{eqnarray*}
\Lag_D & = &  \alpha_1 + \alpha_2 + \mbox{Const.} \\
& & - \frac{1}{2} \Bigl( y_1 y_1 \bx_1^T \bx_1 \alpha_1^2 
        + y_2 y_2 \bx_2^T \bx_2  \alpha_2^2  
        + 2 y_1 y_2 \bx_1^T \bx_2 \alpha_1 \alpha_2 \\
& & \quad + 2 \left(\sum_{i=3}^N \alpha_i y_i \bx_i^T \right) (y_1 \bx_1 \alpha_1 + y_2 \bx_2 \alpha_2) + \mbox{Const.} \Bigr)
\end{eqnarray*}
Let $K_{11} = \bx_1^T \bx_1$, $K_{22}=\bx_2^T \bx_2$,
$K_{12}=\bx_1^T \bx_2$, and
\begin{eqnarray*}
v_j & \equiv & \sum_{i=3}^N \alpha_i y_i \bx_i^T \bx_j \\
 & = &  \bx_j^T \Old{\bw} - \Old{\alpha_1} y_1 \bx_1^T \bx_j 
   - \Old{\alpha_2} y_2 \bx_2^T \bx_j \\
 & = &  (\bx_j^T \Old{\bw} - \Old{b}) + \Old{b} - \Old{\alpha_1} y_1 \bx_1^T \bx_j 
    - \Old{\alpha_2} y_2 \bx_2^T \bx_j \\
 & = &  \Old{u_j} + \Old{b} - \Old{\alpha_1} y_1 \bx_1^T \bx_j 
     - \Old{\alpha_2} y_2 \bx_2^T \bx_j,
\end{eqnarray*}
where $\Old{u_j}=\bx_j^T \Old{\bw} - \Old{b}$ is the output of $\bx_j$
under old parameters.
\begin{eqnarray*}
\Lag_D & = & \alpha_1 + \alpha_2 - \frac{1}{2} \Bigl(
           K_{11}  \alpha_1^2 + K_{22} \alpha_2^2 
            + 2 s K_{12}  \alpha_1 \alpha_2 \\
& & \quad + 2 y_1 v_1 \alpha_1 + 2 y_2 v_2 \alpha_2 \Bigr) + \mbox{Const.} \\
 & = &  \gamma - s\alpha_2 + \alpha_2 - \frac{1}{2} \Bigl(
       K_{11} (\gamma - s\alpha_2)^2 + K_{22} \alpha_2^2  \\
  & & \quad
        + 2 s K_{12}  (\gamma - s\alpha_2) \alpha_2 \\
      & & \quad + 2 y_1 v_1 (\gamma - s\alpha_2) + 2 y_2 v_2 \alpha_2 \Bigr) +
        \mbox{Const.} \\
 &= & (1-s)\alpha_2 - \frac{1}{2} K_{11} (\gamma-s\alpha_2)^2
     - \frac{1}{2} K_{22} \alpha_2^2 - s K_{12} (\gamma - s\alpha_2) \alpha_2 \\
     & & \quad - y_1 v_1 (\gamma - s \alpha_2) 
        - y_2 v_2 \alpha_2 + \mbox{Const.} \\
 &= & (1-s)\alpha_2  - \frac{1}{2} K_{11}\gamma^2 + s K_{11} \gamma \alpha_2 
         - \frac{1}{2} K_{11}s^2 \alpha_2^2 - \frac{1}{2} K_{22} \alpha_2^2 \\ 
     & & \quad
         - s K_{12}\gamma \alpha_2 + s^2 K_{12} \alpha_2^2 - y_1 v_1 \gamma 
          + s y_1 v_1 \alpha_2 - y_2 v_2 \alpha_2 \\
     & & \quad + \mbox{Const.} \\
 &= & (1-s)\alpha_2  + s K_{11} \gamma \alpha_2 
         - \frac{1}{2} K_{11} \alpha_2^2 - \frac{1}{2} K_{22} \alpha_2^2 \\ 
     & & \quad - s K_{12}\gamma \alpha_2 +  K_{12} \alpha_2^2
          + y_2 v_1 \alpha_2 - y_2 v_2 \alpha_2 \\
     & & \quad + \mbox{Const.} \\
 & = & \left(- \frac{1}{2} K_{11}  - \frac{1}{2} K_{22} +  K_{12}\right)
     \alpha_2^2  \\
     & & \quad + \left(1-s + s K_{11} \gamma - s K_{12}\gamma + y_2 v_1 - y_2
     v_2 \right) \alpha_2 \\
     & & \quad + \mbox{Const.}\\
 & = & \frac{1}{2}\left(2K_{12}-K_{11} - K_{22}\right)
     \alpha_2^2  \\
     & & \quad + \left(1-s + s K_{11} \gamma - s K_{12}\gamma + y_2 v_1 - y_2
     v_2 \right) \alpha_2 \\
     & & \quad + \mbox{Const.}
\end{eqnarray*}
Let $\eta\equiv 2 K_{12}-K_{11}-K_{12}$. The coefficient of $\alpha_2$ is
\begin{eqnarray*}
& & 1-s+ s K_{11}\gamma - s K_{12}\gamma + y_2 v_1 - y_2 v_2 \\
& = & 1-s + s K_{11} (\Old{\alpha_1}+s\Old{\alpha_2})
 - s K_{12} (\Old{\alpha_1}+s\Old{\alpha_2}) \\
 & & 
 + y_2 (\Old{u_1} + \Old{b} - \Old{\alpha_1} y_1 K_{11} 
 - \Old{\alpha_2} y_2 K_{12} ) \\
 & & 
 - y_2 (\Old{u_2} + \Old{b} - \Old{\alpha_1} y_1 K_{12} 
 - \Old{\alpha_2} y_2 K_{22} ) \\
 & = & 1-s + s K_{11} \Old{\alpha_1} + K_{11} \Old{\alpha_2}
  - s K_{12} \Old{\alpha_1} - K_{12} \Old{\alpha_2} \\
  & & 
  + y_2 \Old{u_1} + y_2 \Old{b} - s K_{11} \Old{\alpha_1}
  - K_{12} \Old{\alpha_2} \\
  & & 
  - y_2 \Old{u_2} - y_2 \Old{b} + s K_{12} \Old{\alpha_1} 
  + K_{22} \Old{\alpha_2} \\
 & = & 1-s + (s K_{11} - s K_{12} - s K_{11} + s K_{12}) \Old{\alpha_1} \\
  & & 
    + (K_{11}  - 2K_{12} + K_{22}) \Old{\alpha_2}  \\
    & & 
    + y_2( \Old{u_1} - \Old{u_2}) \\
 & = & y_2^2-y_1 y_2 
 + (K_{11}  - 2K_{12} + K_{22}) \Old{\alpha_2}
 + y_2( \Old{u_1} - \Old{u_2})  \\
 & = & y_2( y_2-y_1+\Old{u_1} - \Old{u_2})
 - \eta \Old{\alpha_2} \\
 & = & y_2( (\Old{u_1}-y_1) - (\Old{u_2}-y_2))
 - \eta \Old{\alpha_2} \\
 & = & y_2(\Old{E_1} - \Old{E_2})
 - \eta \Old{\alpha_2}.
\end{eqnarray*}
So the objective function is
\[
\Lag_D = \frac{1}{2}\eta \alpha_2^2 + (y_2(\Old{E_1}-\Old{E_2})-\eta \Old{\alpha_2})\alpha_2 + \mbox{Const.}
\]
The first and second derivatives are
\[
\frac{\mbox{d} \Lag_D}{\mbox{d} \alpha_2}=
\eta \alpha_2 + (y_2(\Old{E_1}-\Old{E_2})-\eta \Old{\alpha_2}),
\]
\[
\frac{\mbox{d}^2 \Lag_D}{\mbox{d} \alpha_2^2}= \eta.
\]

Note that $\eta=2K_{12}-K_{11}-K_{22} \leq 0$. Proof: Let $K_{11}=\bx_1^T \bx_1
$, $K_{12}=\bx_1^T\bx_2$, $K_{22}=\bx_2^T \bx_2$. Then
$\eta = -(\bx_2 - \bx_1)^T (\bx_2 - \bx_1) = - \|\bx_2 - \bx_1\|^2 \leq 0$.

Let $\frac{\mbox{d} \Lag_D}{\mbox{d} \alpha_2}=0$, and we have
\begin{eqnarray*}
\New{\alpha_2} & = & - \frac{y_2(\Old{E_1}-\Old{E_2})
 - \eta \Old{\alpha_2}}{\eta}\\
  & = & \Old{\alpha_2} + \frac{y_2 (\Old{E_2} - \Old{E_1})}
  {\eta}
\end{eqnarray*}
If $\eta < 0$, the above equation gives us the unconstrained maximum point
$\New{\alpha_2}$. It must be checked against the feasible range.
Let $s=y_1 y_2$, and $\gamma = \Old{\alpha_1}+s \Old{\alpha_2}$.
The range of $\alpha_2$ is determined as follows:
\begin{itemize}
  \item If $s=1$, then $\alpha_1+\alpha_2 = \gamma$.
    \begin{itemize}
    \item If $\gamma > C$, then $\max \alpha_2 = C$,
          and $\min \alpha_2 = \gamma - C$.
            \begin{figure}[h]
            \begin{center}
            %\unsetbox\graph
            %\leavevmode
            %\input{pic/fig3}
            % \centerline{\raise 1em\box\graph}
            \includegraphics{pic/fig3}
            \end{center}
            \caption{$\alpha_1+\alpha_2=\gamma$, and $\gamma > C$.}
            \end{figure}
    \item If $\gamma < C$, then $\min \alpha_2 = 0$,
          and $\max \alpha_2 = \gamma$.
            \begin{figure}[h]
            \begin{center}
            %\unsetbox\graph
            %\leavevmode
            %\input{pic/fig4}
            %    \centerline{\raise 1em\box\graph}
            \includegraphics{pic/fig4}
            \end{center}
            \caption{$\alpha_1+\alpha_2=\gamma$, and $\gamma < C$.}
            \end{figure}
    \end{itemize}
  \item If $s=-1$, then $\alpha_1-\alpha_2 = \gamma$.
    \begin{itemize}
    \item If $\gamma > 0$, then $\min \alpha_2 = 0$,
          and $\min \alpha_2 = C - \gamma$.
            \begin{figure}[h]
            \begin{center}
            %\unsetbox\graph
            %\leavevmode
            %\input{pic/fig5}
            %    \centerline{\raise 1em\box\graph}
            \includegraphics{pic/fig5}
            \end{center}
            \caption{$\alpha_1-\alpha_2=\gamma$, and $\gamma > 0$.}
            \end{figure}
    \item If $\gamma < 0$, then $\min \alpha_2 = -\gamma$,
          and $\max \alpha_2 = C$.
            \begin{figure}[h]
            \begin{center}
            %\unsetbox\graph
            %\leavevmode
            %\input{pic/fig6}
            %    \centerline{\raise 1em\box\graph}
            \includegraphics{pic/fig6}
            \end{center}
            \caption{$\alpha_1-\alpha_2=\gamma$, and $\gamma < 0$.}
            \end{figure}
    \end{itemize}
\end{itemize}
Let the minimum feasible value of $\alpha_2$ be $L$,
maximum be $H$. Then
\[
\alpha_2^{\mbox{new,clipped}} = \left\{
\begin{array}{ll}
H, & \mbox{if } H < \New{\alpha_2},\\
\New{\alpha_2}, & \mbox{if } L \leq \New{\alpha_2} \leq H \\
L, & \mbox{if } \New{\alpha_2} < L.
\end{array}
\right.
\]

To summarize, given $\alpha_1$, $\alpha_2$ (and the corresponding
$y_1$, $y_2$, $K_{11}$, $K_{12}$, $K_{22}$, $\Old{E_2}-\Old{E_1}$),
we can optimize the two $\alpha$'s by the following procedure:
\begin{enumerate}
 \item $\eta = 2K_{12}-K_{11}-K_{22}$.
 \item If $\eta < 0$,
       \[
          \Delta\alpha_2 = \frac{y_2(\Old{E_2}-\Old{E_1})}{\eta},
       \]
       and clip the solution within the feasible region.
       Then
       \[
          \Delta\alpha_1 = -s \Delta\alpha_2.
       \]
 \item If $\eta=0$, we need to evaluate the objective function at the
       two endpoints, and set $\New{\alpha_2}$ to be the one with larger
       objective function value. The objective function is
       \begin{equation}
         \Lag_D = \frac{1}{2}\eta \alpha_2^2 + 
         (y_2(\Old{E_1}-\Old{E_2})-\eta\Old{\alpha_2})\alpha_2 + \mbox{Const.}
         \label{eq:ld-alpha2}
       \end{equation}
\end{enumerate}
       
\subsection{SMO Algorithm: Updating after a successful optimization step}

When $\alpha_1$, $\alpha_2$ are changed by $\Delta\alpha_1$, $\Delta\alpha_2$,
we can update $E_i$'s, $F_i$'s, $\bw$ (for linear kernel), and $b$.
Let $E(\bx,y)$ be the prediction error on $(\bx,y)$:
\[
E(\bx,y) = \sum_{i=1}^N \alpha_i y_i \bx_i^T \bx - b - y,
\]
The change in $E$ is
\begin{equation}
\Delta E(\bx,y) = \Delta\alpha_1 y_1 \bx_1^T \bx
+ \Delta\alpha_2 y_2 \bx_2^T \bx - \Delta b.
\label{eq:update:e}
\end{equation}

The change in the threshold can be computed by forcing
$\New{E_1}=0$ if $0<\New{\alpha_1}<C$ (
or $\New{E_2}=0$ if $0<\New{\alpha_2}<C$).
From 
\begin{eqnarray*}
 0 & = & \New{E(\bx,y)}  \\
 & = & \Old{E(\bx,y)} + \Delta E(\bx,y) \\
 & = & \Old{E(\bx,y)} +
 \Delta\alpha_1 y_1 \bx_1^T \bx
  + \Delta\alpha_2 y_2 \bx_2^T \bx - \Delta b
\end{eqnarray*}
we have
\begin{equation}
\Delta b =  \Old{E(\bx,y)} +
 \Delta\alpha_1 y_1 \bx_1^T \bx
   + \Delta\alpha_2 y_2 \bx_2^T \bx.
\label{eq:update:b}
\end{equation}
If $\alpha_1 = 0$, we can only say $y_1 \New{E_1} \geq 0$;
similarly, if $\alpha_1 = C$, we have $y_1\New{E_2} \leq 0$.
If both $\alpha_1$ and $\alpha_2$ take values 0 or $C$,
the original SMO algorithm computes two values of
the new $b$ for $\alpha_1$
and $\alpha_2$ using Equation \ref{eq:update:b},
and takes the average. This is regarded as problematic
by \citet{Keerthi+Shevade+Bhattacharyya+Murthy:2001}.

Similarly, from
\[
F(\bx,y) = \sum_{i=1}^N \alpha_i y_i \bx_i^T \bx - y
\]
we have
\begin{equation}
\Delta F(\bx,y) =  \Delta\alpha_1 y_1 \bx_1^T \bx
 + \Delta\alpha_2 y_2 \bx_2^T \bx.
\label{eq:update:f}
\end{equation}

For the weight vector  of linear kernels,
\[
\bw = \sum_{i=1}^N \alpha_i y_i \bx_i,
\]
\begin{equation}
\Delta\bw = \Delta\alpha_1 y_1 \bx_1 + \Delta\alpha_2 y_2 \bx_2.
\label{eq:update:w}
\end{equation}

\subsection{SMO Algorithm: Pick two $\alpha_i$'s for optimization}

The heuristics for picking two $\alpha_i$'s for optimization
in the original SMO paper are as follows:
\begin{itemize}
 \item The outer loop selects the first $\alpha_i$, the inner loop selects
       the second $\alpha_i$ that maximizes $|E_2-E_1|$.
 \item The outer loop alternates between one sweep through all examples
       and as many sweeps as possible through the non-boundary examples 
       (those with $0 < \alpha_i < C$), 
       selecting the example that violates the KKT condition.
 \item Given the first $\alpha_i$, the inner loop looks for a non-boundary
       that maximizes $|E_2-E_1|$. If this does not make progress, it
       starts a sequential scan through the non-boundary examples, starting
       at a random  position; if this fails too, it startis a sequential
       scan through all the examples, also starting at a random postion.
\end{itemize}
Because the algorithm spends most of the time adjusting the non-boundary
examples, the $E_i$'s of these examples are cached.

The improvement proposed in
\citet{Keerthi+Shevade+Bhattacharyya+Murthy:2001}
avoids the use of the
threshold $b$ in checking KKT condition, and compares two $F_i$'s, which 
also automatically selects the second $\alpha_i$ for joint optimization.
There are two variantions when the outer loop deals only with the non-boundary
examples:
\begin{itemize}
  \item The first $\alpha$ is selected sequentially
        from all the non-boundary examples.
        If the first $\alpha_i$ violates the KKT condition
        when compared with $\alpha_j$ with $F_j = b_{\mbox{low}}$,
        or $F_j = b_{\mbox{up}}$, then select $\alpha_j$ as the second $\alpha$.
  \item The two $\alpha$'s are also those
        with $F_i = b_{\mbox{low}}$
          or $F_i = b_{\mbox{up}}$.
\end{itemize}
After a successful step using a pair of indices, $(i_2, i_1)$, let
$\tilde{I} = I_0 \cup \{i_1,i_2\}$. We claim that each of the two sets,
$\tilde{I}\cap (I_0 \cup I_1 \cup I_2)$ and 
$\tilde{I}\cap (I_0 \cup I_3 \cup I_4)$, is non-empty, hence we can compute
partial $b_{\mbox{low}}$ and $b_{\mbox{up}}$ from the two sets.
\begin{quote}
Proof: The two sets are non-empty if $I_0 \neq \emptyset$.
       If $I_0 = \emptyset$, then $\alpha_1$, $\alpha_2$ can only take values
       from 0 or $C$. They cannot take the same value and $y_1=y_2$ at the
       same time, otherwise we have $0+0=\gamma$ or $C+C=\gamma$ for $\alpha_1
       + s \alpha_2 = \gamma$, in which $\alpha_1$ and $\alpha_2$ cannot be
       changed, which contradicts the fact that we just had a successful step.
       So if they take the same value, with different $y_i$'s, they will
       belong to two different sets. If they take different values, with the
       same $y_i$'s, they will also belong to two different sets.
\end{quote}

\section{\protect\CPLUSPLUS\ Implementation}

Now let's write the {\CPLUSPLUS} code, based on the pseudocode
in \citet{Platt:1998:SMO-BOOK}.

@o c/smo.cc -d
@{@<Header files to include@>
using namespace std;
@<Global variables@>
@<Functions@>
@<Main routine@>
@}

\subsection{The {\tt main} routine}

The \verb|main| routine implements the outer loop that selects
the first $\alpha_i$ for optimization. It alternates between
a sweep through all the examples (\verb|examineAll==1|)
and as many sweeps as possible 
through the non-boundary examples (\verb|examineAll==0|).
If an example \verb|k| violates the KKT condition more than 
\verb|eps|=$\epsilon$, it is selected
as the first $\alpha_i$, and \verb|examineExample(k)| is
called, which returns 1 if positive progress is made to
improve the objective function (with two changed $\alpha_i$'s).

@D Main routine
@{int main(int argc, char *argv[]) {
  @<Variables local to main@>
  int numChanged;
  int examineAll;

  @<Get in parameters@>
  @<Read in data@>

  if (!is_test_only) {
      alph.resize(end_support_i, 0.);

      /* initialize threshold to zero */
      b = 0.;

      /* E_i = u_i - y_i = 0 - y_i = -y_i */
      error_cache.resize(N);

      if (is_linear_kernel)
         w.resize(d,0.);
  }

  @<Initialization@>

  if (!is_test_only) {
      numChanged = 0;
      examineAll = 1;
      while (numChanged > 0 || examineAll) {
          numChanged = 0;
          if (examineAll) { 
              for (int k = 0; k < N; k++)
                numChanged += examineExample (k);
            }
          else { 
              for (int k = 0; k < N; k++)
                if (alph[k] != 0 && alph[k] != C)
                  numChanged += examineExample (k);
            }
          if (examineAll == 1)
            examineAll = 0;
          else if (numChanged == 0)
            examineAll = 1;

          //cerr << error_rate() << endl;
          @<Diagnostic info@>
      }
      @<Write model parameters@>
      cerr << "threshold=" << b << endl;
  }
  cout << error_rate() << endl;
  @<Write classification output@>
}
@}

Let's define the global variables.
@D Global variables
@{
  int N = 0;                    /* N points(rows) */
  int d = -1;                   /* d variables */
  float C=0.05;
  float tolerance=0.001;
  float eps=0.001;
  float two_sigma_squared=2;

  vector<float> alph;           /* Lagrange multipliers */
  float b;                      /* threshold */
  vector<float> w;              /* weight vector: only for linear kernel */

  vector<float> error_cache;

  struct sparse_binary_vector {
      vector<int> id;
    };
  struct sparse_vector {
      vector<int> id;
      vector<float> val;
    };
  typedef vector<float> dense_vector;

  bool is_sparse_data = false;
  bool is_binary = false;
  /* use only one of these */
  vector<sparse_binary_vector> sparse_binary_points;
  vector<sparse_vector> sparse_points;
  vector<dense_vector> dense_points;

  vector<int> target;           /* class labels of training data points */
  bool is_test_only = false;
  bool is_linear_kernel = false;

  /* data points with index in [first_test_i .. N)
   * will be tested to compute error rate
   */
  int first_test_i = 0;

  /*
   * support vectors are within [0..end_support_i)
   */
  int end_support_i = -1;
@}

@D Header files to include
@{
#include <vector>
#include <algorithm>
#include <functional>
@}

\subsection{The {\tt examineExample} routine}

Given the first $\alpha_i$ (with index \verb|i1|),
\verb|examineExample(i1)| first checks if it violates
the KKT condition by more than \verb|tolerance|,
if it does, then looks for the second $\alpha_i$
(with index \verb|i2|) and jointly optimize the two $\alpha_i$'s
by calling \verb|takeStep(i1,i2)|.

@D Functions
@{int examineExample(int i1)
{
  float y1, alph1, E1, r1;

  y1 = target[i1];
  alph1 = alph[i1];

  if (alph1 > 0 && alph1 < C)
     E1 = error_cache[i1];
  else 
     E1 = learned_func(i1) - y1;

  r1 = y1 * E1;
  if ((r1 < -tolerance && alph1 < C)
     || (r1 > tolerance && alph1 > 0))
    { 
      /* Try i2 by three ways; if successful, then immediately return 1; */
      @<Try argmax \verb|E1 - E2|@>
      @<Try iterating through the non-bound examples@>
      @<Try iterating through the entire training set@>
    }

  return 0;
}@}

Use the heuristic to choose the second example from non-bound examples,
so that \verb|E1-E2| is maximized.
@D Try argmax \verb|E1 - E2|
@{{
  int k, i2;
  float tmax;

  for (i2 = (-1), tmax = 0, k = 0; k < end_support_i; k++)
    if (alph[k] > 0 && alph[k] < C) {
        float E2, temp;

        E2 = error_cache[k];
        temp = fabs(E1 - E2);
        if (temp > tmax)
          {
            tmax = temp;
            i2 = k;
          }
      }

  if (i2 >= 0) {
      if (takeStep (i1, i2))
        return 1;
    }
}@}

@D Global variables
@{int takeStep(int i1, int i2);
@}

If we cannot make progress with the best non-bound example,
then try any non-bound examples.
@D Try iterating through the non-bound examples
@{{
  int k, k0;
  int i2;

  for (k0 = (int) (drand48 () * end_support_i), k = k0; k < end_support_i + k0; k++) {
      i2 = k % end_support_i;
      if (alph[i2] > 0 && alph[i2] < C) {
          if (takeStep(i1, i2))
              return 1;
        }
    }
}@}

If we cannot make progress with the non-bound examples,
then try any example.

@D Try iterating through the entire training set
@{{
  int k0, k, i2;

  for (k0 = (int)(drand48 () * end_support_i), k = k0; k < end_support_i + k0; k++) {
      i2 = k % end_support_i;
      if (takeStep(i1, i2))
          return 1;
    }
}@}

\subsection{The {\tt takeStep} routine}

Now let's write \verb|takeStep| which optimizes two Lagrange multipliers.
If successful, return 1, else return 0.

@D Functions
@{int takeStep(int i1, int i2) { 
  int y1, y2, s;
  float alph1, alph2; /* old_values of alpha_1, alpha_2 */
  float a1, a2;       /* new values of alpha_1, alpha_2 */
  float E1, E2, L, H, k11, k22, k12, eta, Lobj, Hobj;

  if (i1 == i2) return 0;

  @<Look up \verb|alph1, y1, E1, alph2, y2, E2|@>
  s = y1 * y2;

  @<Compute L, H@>
  if (L == H)
    return 0;

  @<Compute \verb|eta|@>
  if (eta < 0) {
      a2 = alph2 + y2 * (E2 - E1) / eta;
      if (a2 < L)
        a2 = L;
      else if (a2 > H)
        a2 = H;
    }
  else {
    @<Compute Lobj, Hobj: objective function at a2=L, a2=H@>
    if (Lobj > Hobj+eps)
      a2 = L;
    else if (Lobj < Hobj-eps)
      a2 = H;
    else
      a2 = alph2;
  }

  if (fabs(a2-alph2) < eps*(a2+alph2+eps))
    return 0;

  a1 = alph1 - s * (a2 - alph2);
  if (a1 < 0) {
        a2 += s * a1;
        a1 = 0;
  }
  else if (a1 > C) {
        float t = a1-C;
        a2 += s * t;
        a1 = C;
  }

  @<Update threshold to reflect change in Lagrange multipliers@>
  @<Update weight vector to reflect change in a1 and a2, if linear SVM@>
  @<Update error cache using new Lagrange multipliers@>

  alph[i1] = a1;  /* Store a1 in the alpha array.*/
  alph[i2] = a2;  /* Store a2 in the alpha array.*/

  return 1;
}@}

As the SMO algorithm spends most of its time on adjusting
the $\alpha_i$'s of the non-boundary examples,
an error cache is maintained for them.
Each time after a successful optimization step, for the two $\alpha_i$'s,
if $0<\alpha_i<C$, the corresponding $E_i$ is set zero.
The $E_i$'s for other $\alpha_i$'s (that have been kept fixed during
the optimization step) is updated using Equation \ref{eq:update:e}.

@D Look up \verb|alph1, y1, E1, alph2, y2, E2|
@{  alph1 = alph[i1];
  y1 = target[i1];
  if (alph1 > 0 && alph1 < C)
     E1 = error_cache[i1];
  else 
     E1 = learned_func(i1) - y1;

  alph2 = alph[i2];
  y2 = target[i2];
  if (alph2 > 0 && alph2 < C)
     E2 = error_cache[i2];
  else 
     E2 = learned_func(i2) - y2;
@}
@D Global variables
@{ float (*learned_func)(int) = NULL;
@}

Compute the feasible range of $\New{\alpha_2}$. See the graphs
on Page \pageref{page:L-H-graph}.

@D Compute L, H
@{
   if (y1 == y2) {
     float gamma = alph1 + alph2;
     if (gamma > C) {
         L = gamma-C;
         H = C;
     }
     else {
         L = 0;
         H = gamma;
     }
   }
   else {
       float gamma = alph1 - alph2;
       if (gamma > 0) {
           L = 0;
           H = C - gamma;
       }
       else {
           L = -gamma;
           H = C;
       }
   }
@}

@D Compute \verb|eta|
@{  k11 = kernel_func(i1, i1);
  k12 = kernel_func(i1, i2);
  k22 = kernel_func(i2, i2);
  eta = 2 * k12 - k11 - k22;
@}

@D Global variables
@{
 float (*kernel_func)(int,int)=NULL;
@}

See Equation \ref{eq:ld-alpha2} on Page \pageref{eq:ld-alpha2} for evaluating $\Lag_D$ at $\alpha_2$.
@D Compute Lobj, Hobj: objective function at a2=L, a2=H
@{{
    float c1 = eta/2;
    float c2 = y2 * (E1-E2)- eta * alph2;
    Lobj = c1 * L * L + c2 * L;
    Hobj = c1 * H * H + c2 * H;
}@}

See Equation \ref{eq:update:b} on Page \pageref{eq:update:b} 
for updating the threshold $b$
when either of $\alpha_1$ and $\alpha_2$ are non-boundary.

@D Update threshold to reflect change in Lagrange multipliers
@{{
  float b1, b2, bnew;

  if (a1 > 0 && a1 < C)
    bnew = b + E1 + y1 * (a1 - alph1) * k11 + y2 * (a2 - alph2) * k12;
  else {
      if (a2 > 0 && a2 < C)
          bnew = b + E2 + y1 * (a1 - alph1) * k12 + y2 * (a2 - alph2) * k22;
      else {
          b1 = b + E1 + y1 * (a1 - alph1) * k11 + y2 * (a2 - alph2) * k12;
          b2 = b + E2 + y1 * (a1 - alph1) * k12 + y2 * (a2 - alph2) * k22;
          bnew = (b1 + b2) / 2;
      }
  }

  delta_b = bnew - b;
  b = bnew;
}@}

@D Global variables
@{ float delta_b;
@}

A linear SVM can be sped up by only using the weight vector
(rather than all of the training examples that correspond to non-zero Lagrange
multipliers)
when evaluating the learned classification function.

If the joint optimization succeeds, this stored weight vector
must be updated to reflect the new Lagrange multiplier values.
See Equation \ref{eq:update:w} on Page \pageref{eq:update:w}.

@D Update weight vector to reflect change in a1 and a2, if linear SVM
@{
  if (is_linear_kernel) {
    float t1 = y1 * (a1 - alph1);
    float t2 = y2 * (a2 - alph2);

    if (is_sparse_data && is_binary) {
        int p1,num1,p2,num2;

        num1 = sparse_binary_points[i1].id.size();
        for (p1=0; p1<num1; p1++)
            w[sparse_binary_points[i1].id[p1]] += t1;

        num2 = sparse_binary_points[i2].id.size();
        for (p2=0; p2<num2; p2++)
            w[sparse_binary_points[i2].id[p2]] += t2;
    }
    else if (is_sparse_data && !is_binary) {
        int p1,num1,p2,num2;

        num1 = sparse_points[i1].id.size();
        for (p1=0; p1<num1; p1++)
            w[sparse_points[i1].id[p1]] += 
                t1 * sparse_points[i1].val[p1];

        num2 = sparse_points[i2].id.size();
        for (p2=0; p2<num2; p2++)
            w[sparse_points[i2].id[p2]] +=
                t2 * sparse_points[i2].val[p2];
    }
    else
        for (int i=0; i<d; i++)
            w[i] += dense_points[i1][i] * t1 + dense_points[i2][i] * t2;
}@}

See Equation \ref{eq:update:e} on Page \pageref{eq:update:e}.
@D Update error cache using new Lagrange multipliers
@{{
    float t1 = y1 * (a1-alph1);
    float t2 = y2 * (a2-alph2);

    for (int i=0; i<end_support_i; i++)
    if (0 < alph[i] && alph[i] < C)
        error_cache[i] +=  t1 * kernel_func(i1,i) + t2 * kernel_func(i2,i)
                          - delta_b;
    error_cache[i1] = 0.;
    error_cache[i2] = 0.;
}@}

\subsection{Evaluating classification function}

We use a function pointer \verb|learned_func| to represent
the learned function $f(\bx)=\bw \cdot \bx - b$.
It takes the index of the data point $k$,
and computes $f(\bx_k)$.

According to the kernel type and input data type,
we define the following functions for evaluating the learned classification
function.

@D Functions
@{float learned_func_linear_sparse_binary(int k) {
  float s = 0.;

  for (int i=0; i<sparse_binary_points[k].id.size(); i++)
    s += w[sparse_binary_points[k].id[i]];

  s -= b;
  return s;
}@}

@D Functions
@{float learned_func_linear_sparse_nonbinary(int k) {
  float s = 0.;

  for (int i=0; i<sparse_points[k].id.size(); i++)
  {
      int j = sparse_points[k].id[i];
      float v = sparse_points[k].val[i];
      s += w[j] * v;
  }
  s -= b;
  return s;
}@}

@D Functions
@{float learned_func_linear_dense(int k) {
  float s = 0.;

  for (int i=0; i<d; i++)
    s += w[i] * dense_points[k][i];

  s -= b;
  return s;
}@}

@D Functions
@{float learned_func_nonlinear(int k) {
  float s = 0.;
  for (int i=0; i<end_support_i; i++)
      if (alph[i] > 0)
          s += alph[i]*target[i]*kernel_func(i,k);
  s -= b;
  return s;
}@}

During initialization, we point \verb|learned_func| one of
functions below.
@D Initialization
@{
    if (is_linear_kernel && is_sparse_data && is_binary)
        learned_func = learned_func_linear_sparse_binary;
    if (is_linear_kernel && is_sparse_data && !is_binary)
        learned_func = learned_func_linear_sparse_nonbinary;
    if (is_linear_kernel && !is_sparse_data)
        learned_func = learned_func_linear_dense;
    if (!is_linear_kernel)
        learned_func = learned_func_nonlinear;
@}

\subsection{Functions to compute dot product}

Accroding to the input data type, we have different functions
to compute the dot product of two data points.

@D Global variables
@{
 float (*dot_product_func)(int,int)=NULL;
@}

@D Initialization
@{
    if (is_sparse_data && is_binary)
        dot_product_func = dot_product_sparse_binary;
    if (is_sparse_data && !is_binary)
        dot_product_func = dot_product_sparse_nonbinary;
    if (!is_sparse_data)
        dot_product_func = dot_product_dense;
@}

@D Functions
@{
float dot_product_sparse_binary(int i1, int i2)
{
  int p1=0, p2=0, dot=0;
  int num1 = sparse_binary_points[i1].id.size();
  int num2 = sparse_binary_points[i2].id.size();

  while (p1 < num1 && p2 < num2) {
      int a1 = sparse_binary_points[i1].id[p1];
      int a2 = sparse_binary_points[i2].id[p2];
      if (a1 == a2) {
          dot++;
          p1++;
          p2++;
        }
      else if (a1 > a2)
        p2++;
      else
        p1++;
    }
  return (float)dot;
}@}

@D Functions
@{
float dot_product_sparse_nonbinary(int i1, int i2)
{
  int p1=0, p2=0;
  float dot = 0.;
  int num1 = sparse_points[i1].id.size();
  int num2 = sparse_points[i2].id.size();

  while (p1 < num1 && p2 < num2) {
      int a1 = sparse_points[i1].id[p1];
      int a2 = sparse_points[i2].id[p2];
      if (a1 == a2) {
          dot += sparse_points[i1].val[p1] * sparse_points[i2].val[p2];
          p1++;
          p2++;
        }
      else if (a1 > a2)
        p2++;
      else
        p1++;
    }
  return (float)dot;
}@}

@D Functions
@{
float dot_product_dense(int i1, int i2)
{
  float dot = 0.;
  for (int i=0; i<d; i++)
    dot += dense_points[i1][i] * dense_points[i2][i];

  return dot;
}
@}

\subsection{Kernel functions}

The linear kernel is simply the dot product.
Currently, we have only one nonlinear kernel: 
radial basis function kernel.

@D Initialization
@{
    if (is_linear_kernel)
        kernel_func = dot_product_func;
    if (!is_linear_kernel)
        kernel_func = rbf_kernel;
@}

The calculation of $\|\bx_1 - \bx_2 \|^2$ in a Gaussian kernel
can be sped up using the following equation:
\begin{eqnarray*}
\|\bx_1 - \bx_2 \|^2 & = & (\bx_1 - \bx_2)^T (\bx_1 - \bx_2) \\
& = & \bx_1^T \bx_1 + \bx_2^T \bx_2 - 2 \bx_1^T \bx_2,
\end{eqnarray*}
where $\bx_i^T \bx_i$ can be pre-computed. For each of the $d$ dimensions,
directly computing $\|\bx_1 - \bx_2 \|^2$ needs 3 operations:
\begin{enumerate}
 \item
     $a = x_{1j} - x_{2j}$
 \item
     $b = a \times a$
 \item
     $s = s + b$.
\end{enumerate}
In comparison, the new method needs only 2 operations:
\begin{enumerate}
 \item
     $a = x_{1j} \times x_{2j}$
 \item
     $s = s + a$.
\end{enumerate}

@D Functions
@{
float rbf_kernel(int i1, int i2)
{
    float s = dot_product_func(i1,i2);
    s *= -2;
    s += precomputed_self_dot_product[i1] + precomputed_self_dot_product[i2];
    return exp(-s/two_sigma_squared);
}
@}
@D Header files to include
@{#include <cmath>
@}

@D Initialization
@{
    if (!is_linear_kernel) {
        precomputed_self_dot_product.resize(N);
        for (int i=0; i<N; i++)
            precomputed_self_dot_product[i] = dot_product_func(i,i);
    }
@}

@D Global variables
@{
    vector<float> precomputed_self_dot_product;
@}

\subsection{Input and output}

\subsubsection{Get parameters by command line}

Finally, we have to get the parameters \verb|C, eps, tolerance|, etc. (And
the input/output file names.)
We use the \verb|getopt()| routine to handle this.

@D Get in parameters
@{{
  extern char *optarg;
  extern int optind;
  int c;
  int errflg = 0;

  while ((c = getopt (argc, argv, "n:d:c:t:e:p:f:m:o:r:lsba")) != EOF)
    switch (c)
      {
      case 'n':
        N = atoi(optarg);
        break;
      case 'd':
        d = atoi(optarg);
        break;
      case 'c':
        C = atof (optarg);
        break;
      case 't':
        tolerance = atof(optarg);
        break;
      case 'e':
        eps = atof (optarg);
        break;
      case 'p':
        two_sigma_squared = atof (optarg);
        break;
      case 'f':
        data_file_name = optarg;
        break;
      case 'm':
        svm_file_name = optarg;
        break;
      case 'o':
        output_file_name = optarg;
        break;
      case 'r':
        srand48 (atoi (optarg));
        break;
      case 'l':
        is_linear_kernel = true;
        break;
      case 's':
        is_sparse_data = true;
        break;
      case 'b':
        is_binary = true;
        break;
      case 'a':
        is_test_only = true;
        break;
      case '?':
        errflg++;
      }

  if (errflg || optind < argc)
    {
      cerr << "usage: " << argv[0] << " " <<
                "-f  data_file_name\n"
                "-m  svm_file_name\n"
                "-o  output_file_name\n"
                "-n  N\n"
                "-d  d\n"
                "-c  C\n"
                "-t  tolerance\n"
                "-e  epsilon\n"
                "-p  two_sigma_squared\n"
                "-r  random_seed\n"
                "-l  (is_linear_kernel)\n"
                "-s  (is_sparse_data)\n"
                "-b  (is_binary)\n"
                "-a  (is_test_only)\n"
                ;
      exit (2);
    }
}@}

@D Variables local to main
@{
char *data_file_name = "svm.data";
char *svm_file_name = "svm.model";
char *output_file_name = "svm.output";
@}

@D Header files to include
@{#include <iostream>
#include <cstdlib>
#include <unistd.h>
@}

\subsubsection{Read in data}

The data file is a flat text file, each data point occupies
one line in which the class label ($+1$ or $-1$) 
follows the attribute values.
Ordinarily, a line will be
\begin{verbatim}
attribute_1_value attribute_2_value ... attribute_d_value target_value
\end{verbatim}
For sparse format, a line will be
\begin{verbatim}
id_1 val_1 id_2 val_2 ... id_m val_m target_value
\end{verbatim}
where \verb|id_j| should be between 1 and $d$.
For sparse binary format, a line will be
\begin{verbatim}
id_1 id_2 ... id_m target_value
\end{verbatim}
here, too, \verb|id_j| should be between 1 and $d$.

The data is read into \verb|dense_points|, or \verb|sparse_points|,
or \verb|sparse_binary_points|, according to the input format.

@D Read in data
@{{
  int n;
  if (is_test_only) {
      ifstream svm_file(svm_file_name);
      end_support_i = first_test_i = n = read_svm(svm_file);
      N += n;
  }
  if (N > 0) {
    target.reserve(N);
    if (is_sparse_data && is_binary)
        sparse_binary_points.reserve(N);
    else if (is_sparse_data && !is_binary)
        sparse_points.reserve(N);
    else
        dense_points.reserve(N);
  }
  ifstream data_file(data_file_name);
  n = read_data(data_file);
  if (is_test_only) {
      N = first_test_i + n;
  }
  else {
      N = n;
      first_test_i = 0;
      end_support_i = N;
  }
}@}

The actually reading of data is handled by \verb|read_data(istream&)|;
it appends data points from the input stream to \verb|dense_points|
(or \verb|sparse_points|, or \verb|sparse_binary_points|, depending
on input format).

The function \verb|read_data(istream&)| may be called by \verb|read_svm()|
to read in the support vectors of previous trained model (if non-linear
kernel is used). These support vectors are read into \verb|dense_points|
(or \verb|sparse_points|, or \verb|sparse_binary_points|, depending
on input format) before the data points in the data file. The starting index
of the data points in the data file is \verb|first_test_i|.

@D Functions
@{int read_data(istream& is)
{
  string s;
  int n_lines;

  for (n_lines = 0; getline(is, s, '\n'); n_lines++) {
    istrstream line(s.c_str());
    vector<float> v;
    float t;
    while (line >> t)
       v.push_back(t);
    target.push_back(v.back());
    v.pop_back();
    int n = v.size();
    if (is_sparse_data && is_binary) {
        sparse_binary_vector x;
        for (int i=0; i<n; i++) {
            if (v[i] < 1 || v[i] > d) {
                cerr << "error: line " << n_lines+1
                    << ": attribute index " << int(v[i]) << " out of range."<<endl;
                exit(1);
            }
            x.id.push_back(int(v[i])-1);
        }
        sparse_binary_points.push_back(x);
    }
    else if (is_sparse_data && !is_binary) {
        sparse_vector x;
        for (int i=0; i<n; i+=2) {
            if (v[i] < 1 || v[i] > d) {
                cerr << "data file error: line " << n_lines+1
                     << ": attribute index " << int(v[i]) << " out of range."
                     << endl;
                exit(1);
            }
            x.id.push_back(int(v[i])-1);
            x.val.push_back(v[i+1]);
        }
        sparse_points.push_back(x);
    }
    else {
        if (v.size() != d) {
                cerr << "data file error: line " << n_lines+1 
                     << " has " << v.size() << " attributes; should be d=" << d
                     <<endl;
                exit(1);
        }
        dense_points.push_back(v);
    }
  }
  return n_lines;
}@}

@D Header files to include
@{#include <string>
#include <vector>
#include <iostream>
#include <fstream>
#include <strstream>
@}

\subsubsection{Saving and loading model parameters}

The output order of the model paramters will be
\begin{enumerate}
  \item The number of attributes $d$.
  \item The flag \verb|is_sparse_data|
  \item The flag \verb|is_binary|
  \item The flag \verb|is_linear_kernel|
  \item The threshold $b$
  \item If the linear kernel is used:
    \begin{enumerate}
      \item The weight vector $\bw$
    \end{enumerate}
  \item If non-linear kernel is used
    \begin{enumerate}
      \item Kernel paramters (e.g., $2\sigma^2$ for radial basis function 
                             kernel)
      \item The number of support vectors
      \item The Lagrange multipliers of the support vectors
      \item The support vectors, one per line
    \end{enumerate}
\end{enumerate}

@D Functions
@{
void write_svm(ostream& os) {
    os << d << endl;
    os << is_sparse_data << endl;
    os << is_binary << endl;
    os << is_linear_kernel << endl;
    os << b << endl;
    if (is_linear_kernel) {
        for (int i=0; i<d; i++)
            os << w[i] << endl;
    }
    else {
        os << two_sigma_squared << endl;
        int n_support_vectors=0;
        for (int i=0; i<end_support_i; i++)
            if (alph[i] > 0)
                n_support_vectors++;
        os << n_support_vectors << endl;
        for (int i=0; i<end_support_i; i++)
            if (alph[i] > 0)
                os << alph[i] << endl;
        for (int i=0; i<end_support_i; i++)
            if (alph[i] > 0) {
                if (is_sparse_data && is_binary) {
                    for (int j=0; j<sparse_binary_points[i].id.size(); j++)
                        os << (sparse_binary_points[i].id[j]+1) << ' ';
                }
                else if (is_sparse_data && !is_binary) {
                    for (int j=0; j<sparse_points[i].id.size(); j++)
                        os << (sparse_points[i].id[j]+1) << ' ' 
                           << sparse_points[i].val[j] << ' ';
                }
                else {
                    for (int j=0; j<d; j++)
                        os << dense_points[i][j] << ' ';
                }
                os << target[i];
                os << endl;
            }
    }
}@}

@D Functions
@{
int read_svm(istream& is) {
    is >> d;
    is >> is_sparse_data;
    is >> is_binary;
    is >> is_linear_kernel;
    is >> b;
    if (is_linear_kernel) {
        w.resize(d);
        for (int i=0; i<d; i++)
            is >> w[i];
    }
    else {
        is >> two_sigma_squared;
        int n_support_vectors;
        is >> n_support_vectors;
        alph.resize(n_support_vectors, 0.);
        for (int i=0; i<n_support_vectors; i++)
                is >> alph[i];
        string dummy_line_to_skip_newline;
        getline(is, dummy_line_to_skip_newline, '\n');
        return read_data(is);
    }
    return 0;
}@}

@D Write model parameters
@{{
if (!is_test_only && svm_file_name != NULL) {
    ofstream svm_file(svm_file_name);
    write_svm(svm_file);
}
}@}

\subsection{Compute error rate}

@D Functions
@{
float
error_rate()
{
    int n_total = 0;
    int n_error = 0;
    for (int i=first_test_i; i<N; i++) {
       if (learned_func(i) > 0 != target[i] > 0)
           n_error++;
       n_total++;
    }
    return float(n_error)/float(n_total);
}
@}

The classification output is $\bw \bx_i - b$ for each data point
$\bx_i$, one per line.
@D Write classification output
@{{
    ofstream output_file(output_file_name);
    for (int i=first_test_i; i<N; i++)
       output_file << learned_func(i) << endl;
}@}

@D Diagnostic info
@{
/* L_D */
{
#if 0
  float s = 0.;
  for (int i=0; i<N; i++)
     s += alph[i];
  float t = 0.;
  for (int i=0; i<N; i++)
      for (int j=0; j<N; j++)
         t += alph[i]*alph[j]*target[i]*target[j]*kernel_func(i,j);
  cerr << "Objective function=" << (s - t/2.) << endl;
  for (int i=0; i<N; i++)
      if (alph[i] < 0)
         cerr << "alph[" << i << "]=" << alph[i] << " < 0" << endl;
  s = 0.;
  for (int i=0; i<N; i++)
    s += alph[i] * target[i];
  cerr << "s=" << s << endl;
  cerr << "error_rate=" << error_rate() << '\t';
#endif
  int non_bound_support =0;
  int bound_support =0;
  for (int i=0; i<N; i++)
    if (alph[i] > 0) {
       if (alph[i] < C)
          non_bound_support++;
       else
          bound_support++;
    }
  cerr << "non_bound=" << non_bound_support << '\t';
  cerr << "bound_support=" << bound_support << endl;
}
@}

@D Is the objective function increasing?
@{{
    float c1 = eta/2;
    float c2 = y2 * (E1-E2)- eta * alph2;
    float t1 = c1 * alph2 * alph2 + c2 * alph2;
    float t2 = c1 * a2 * a2 + c2 * a2;
    if (t2-t1 < 0)
        cerr <<  "change=" << t2 - t1 << endl;
}@}
\subsection{Multiclass}

The SMO code handles only binary classification. To handle
the multiclass case, we use the following script.
The input data format is similar to that of SMO, except
the class labels with be 0, 1, $\ldots$, $n-1$, if there
are $n$ classes.
The script \verb|smo_multi_class| builds $n$ binary classifiers,
$f_c(\bx)=sgn(\bw_c \bx - b_c)$,
one for each of the $c$ classes.
The classification rule of the multiclass classifier is
\[
\hat{c} = \arg\max_c \bw_c\bx - b_c.
\]

The models are saved in
\verb|${svm_file_name_prefix}|.$c$,
where $c$ = 0, 1, \ldots, $n-1$.

The $i$th line of the classification output file contains 
the $n$ values of $\bw_c\bx_i - b_c$ of the data point $\bx_i$.

@O scripts/smo_multi_class -x
@{#!/bin/sh
## smo_multi_class: multi-class wrapper for SMO
##   Usage: smo_multi_class options -- smo-options
##      options must include:
##           -c number-of-classes
##           -f data-file-name
##           -o output-file-name
##           -m svm-file-name-prefix
##      The `smo-options' after `--' are passed to smo.

if [ $# -lt 8 ]
then
    sed -n '/^##/s/^## //p' $0 >&2
    exit 1
fi

number_of_classes=0
data_file_name=NULL
output_file_name=NULL
svm_file_name_prefix=NULL

while getopts c:f:o:m: c
do
    case $c in
        c) number_of_classes=$OPTARG;;
        f) data_file_name=$OPTARG;;
        o) output_file_name=$OPTARG;;
        m) svm_file_name_prefix=$OPTARG;;
       \?) sed -n '/^##/s/^## //p' $0 >&2
           exit 1;;
    esac
done
shift `expr $OPTIND - 1`

if [ $output_file_name = NULL ] || [ $svm_file_name_prefix = NULL ]
then
       sed -n '/^##/s/^## //p' $0 >&2
       exit 1
fi

if [ $number_of_classes -ge 2 ]
then
    :
else
    echo "error: invalid number of classes ($number_of_classes); should be >= 2" >&2
    exit 1
fi

if [ ! -f $data_file_name ]
then
    echo "error: cannot open data file: $data_file_name" >&2
    exit 1
fi

tmp_data_file_name=../tmp/multiclasstmpsvm.data
tmp_output_file_name=../tmp/multiclasstmpsvm.output
all_target_file_name=../tmp/multiclasstmpsvm.all_target
cat $data_file_name | awk '{ print $NF }' > $all_target_file_name
printf "" > $output_file_name
i=0
while [ $i -lt $number_of_classes ]
do
    printf "class $i: "
    individual_svm_file_name=${svm_file_name_prefix}.$i
    cat $data_file_name |
    awk '{ for (i=1; i<NF; i++)
             printf("%s ",$i);
           if ($NF == '$i')
             printf("1\n");
           else
             printf("-1\n");
         }' > $tmp_data_file_name
    ../c/smo "$@@" \
         -f $tmp_data_file_name \
         -o $tmp_output_file_name \
         -m $individual_svm_file_name
    paste $output_file_name  $tmp_output_file_name > ${output_file_name}.tmp
    mv ${output_file_name}.tmp $output_file_name 
    rm $tmp_data_file_name $tmp_output_file_name
    i=`expr $i + 1`
done
printf "multi-class: "
paste $output_file_name $all_target_file_name |
awk 'BEGIN { n_total = 0.
             n_error = 0.
           }
     {
        best_val = $1
        best_i = 1
        for (i=2; i<NF; i++)
            if ($i > best_val) {
                best_val = $i
                best_i = i
            }
        best_i--
        if (best_i != $NF)
            n_error++
        n_total++
     }
     END { print n_error/n_total }'
rm $all_target_file_name
@}

\subsection{Makefiles}

@O Makefile -t
@{
all: smo.tex smo.dvi smo.ps ccode
ccode:
	cd c; make
smo.dvi:        smo.tex \
                pic/fig1.eps pic/fig2.eps pic/fig3.eps \
                pic/fig4.eps pic/fig5.eps pic/fig6.eps \
                pic/I0-I4.eps
pic/fig1.eps:   pic/fig1.pic
pic/fig2.eps:   pic/fig2.pic
pic/fig3.eps:   pic/fig3.pic
pic/fig4.eps:   pic/fig4.pic
pic/fig5.eps:   pic/fig5.pic
pic/fig6.eps:   pic/fig6.pic

smo.pdf:        smo.tex \
                pic/fig1.pdf pic/fig2.pdf pic/fig3.pdf \
                pic/fig4.pdf pic/fig5.pdf pic/fig6.pdf \
                pic/I0-I4.pdf
pic/I0-I4.pdf:  pic/I0-I4.eps
pic/fig1.pdf:   pic/fig1.eps
pic/fig2.pdf:   pic/fig2.eps
pic/fig3.pdf:   pic/fig3.eps
pic/fig4.pdf:   pic/fig4.eps
pic/fig5.pdf:   pic/fig5.eps
pic/fig6.pdf:   pic/fig6.eps

dist:
	sh scripts/create-dist-tar
include $(HOME)/doc/rules.mk
@}

@O c/Makefile -t
@{
all: smo

# CXXFLAGS=-g
CXXFLAGS=-O3
@}

@o tmp/dummy -n
@{@}

\bibliographystyle{plainnat}
%\bibliographystyle{abbrv}
\bibliography{xge}

\appendix
\section{The weight vectors of the parallel supporting planes}
\label{appendix:1}

Suppose $H_1:  \ba \cdot \bx - b_1 = 0$, and
$H_2:  \ba \cdot \bx - b_2 = 0$
are the two parallel planes. Because they are parallel, they
can have the same weight vector $\ba$.
Let $b'=\frac{b_1 + b_2}{2}$, and $\delta = b_1 - b'$.
So $b_1 = b' + \delta$ and $b_2 = b' - \delta$.
We can rewrite the equations as
\begin{eqnarray*}
&  H_1: & \ba \cdot \bx - (b' + \delta)  = 0 \\
& H_2: & \ba \cdot \bx - (b' - \delta)  = 0
\end{eqnarray*}
or
\begin{eqnarray*}
&  H_1: & \ba \cdot \bx - b' =  \delta \\
& H_2: & \ba \cdot \bx - b' = -\delta 
\end{eqnarray*}
Divide the equations by $\delta$, we have
\begin{eqnarray*}
&  H_1: & \frac{1}{\delta}\ba \cdot \bx -\frac{b'}{\delta} =  1
\\
& H_2: &\frac{1}{\delta} \ba \cdot \bx -\frac{b'}{\delta} = -1 
\end{eqnarray*}

Let $\bw'= \frac{1}{\delta} \ba$ and
$b= \frac{b'}{\delta}$, we have
\begin{eqnarray*}
&  H_1: & \bw \cdot \bx - b = + 1 \\
&  H_2: & \bw \cdot \bx - b = - 1
\end{eqnarray*}

\section{The objective function of the dual problem}
\label{appendix:2}

For the convex quadratic primal problem
\begin{eqnarray*}
    \minimize_{\bw,b,\xi_i} & \frac{1}{2}\bw^T \bw + C \sum_{i=1}^N \xi_i & \\
    \mbox{subject to } &  y_i (\bw^T \bx_i - b) + \xi_i - 1 \geq 0, & 
         1\leq i \leq N \\
                       & \xi_i \geq 0, & 1\leq i \leq N,
\end{eqnarray*}
the Lagrangian is
\begin{eqnarray*}
\Lag(\bw,b,\xi_i; \balpha, \bmu)  & = &
 \frac{1}{2}\bw^T \bw + C \sum_{i=1}^N \xi_i \\
& &
- \sum_{i=1}^N \alpha_i \bigl[ y_i (\bw^T \bx_i - b) + \xi_i - 1 \bigr]
- \sum_{i=1}^N \mu_i \xi_i \\
& = &  \frac{1}{2}\bw^T \bw 
+ \sum_{i=1}^N \left( C - \alpha_i - \mu_i \right) \xi_i \\
& &  
- \left( \sum_{i=1}^N \alpha_i y_i \bx_i^T \right) \bw
- \left( \sum_{i=1}^N \alpha_i y_i \right) b + \sum_{i=1}^N \alpha_i,
\end{eqnarray*}
where $\balpha$, $\bmu$ are the Lagrange multipliers,
the Wolfe dual problem is
\begin{eqnarray*}
    \maximize_{\balpha, \bmu} & \Lag(\bw,b,\xi_i; \balpha, \bmu) & \\
    \mbox{subject to } &  \Part{\Lag}{\bw} = \bzero & \\
                       &  \Part{\Lag}{b} = 0 & \\
                       &  \Part{\Lag}{\xi_i} = 0 & 1\leq i \leq N \\
                       &  \balpha \geq \bzero & \\
                       &  \bmu \geq \bzero.
\end{eqnarray*}
The constraint $ \Part{\Lag}{\bw} = \bzero$ implies
\[
\bw^T - \sum_{i=1}^N  \alpha_i y_i \bx_i^T = \bzero,
\]
or, equivalently,
\[
\bw =  \sum_{i=1}^N \alpha_i y_i \bx_i.
\]
The constraint $\Part{\Lag}{b} = 0$ implies
\[
\sum_{i=1}^N \alpha_i y_i = 0.
\]
The constraints $\Part{\Lag}{\xi_i} = 0$ imply
\[
C-\alpha_i - \mu_i = 0, \qquad 1\leq i \leq N.
\]
Note that $ \balpha \geq \bzero$, $ \bmu \geq \bzero$, and we have
\[
0 \leq \alpha_i \leq C.
\]
Substituting these results into $\Lag(\bw,b,\xi_i; \balpha, \bmu)$:
\begin{eqnarray*}
\Lag(\bw,b,\xi_i; \balpha, \bmu)  & = &
\frac{1}{2}\bw^T \bw + \sum_{i=1}^N 0 \times \xi_i - \bw^T \bw - 0 \times b 
 + \sum_{i=1}^N \alpha_i \\
& = & - \frac{1}{2}\bw^T \bw + \sum_{i=1}^N \alpha_i \\
& = & \sum_{i=1}^N \alpha_i - \frac{1}{2} \sum_{i=1}^{N} \sum_{j=1}^N
        y_i y_j \bx_i^T \bx_j \alpha_i \alpha_j
\end{eqnarray*}
To summarize, the dual problem is
\begin{eqnarray*}
    \maximize_{\balpha} & \Lag_D = \sum\limits_{i=1}^N \alpha_i 
        - \frac{1}{2} \sum\limits_{i=1}^{N} \sum\limits_{j=1}^N
        y_i y_j \bx_i^T \bx_j \alpha_i \alpha_j & \\
    \mbox{subject to } &  \sum\limits_{i=1}^N y_i \alpha_i = 0 & \\
                       &  0 \leq \alpha_i \leq C & 1 \leq i \leq N.
\end{eqnarray*}
\end{document}