\documentclass[twoside,11pt]{article}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\usepackage{theorem,jmlr2e,amsmath,latexsym,float,alltt,epsfig,curves,epic,subfigure,url,multirow,afterpage,enumerate}
\usepackage{pstcol,pst-node,pstricks,color}
%\usepackage{color}
\definecolor{darkblue}{rgb}{0.1,0.1,0.5}
\definecolor{lightblue}{rgb}{.68,.85,.9}
\definecolor{problem}{rgb}{0.9,0.8,0.8}
\newcommand{\Proof}{\noindent{\bf Proof.}~}
% \newcommand{\qed}{$ \blacksquare $ \medskip}
\newcommand{\qed}{$ \Box $ \medskip}
\newcommand{\libsvm}{$\mbox{{\sf LIBSVM}}$}
\newcommand{\bsvm}{$\mbox{{\sf BSVM}}$}
\newcommand{\svmlight}{$SVM^{light}$}
%\numberwithin{equation}{section}
%\renewcommand{\thetheorem}{\thesection.\arabic{theorem}}
%\renewcommand{\thefigure}{\thesection.\arabic{figure}}
%\renewcommand{\thetable}{\thesection.\arabic{table}}
%\renewcommand{\thefootnote}{\fnsymbol{footnote}}
\newcommand{\tv}{\mbox{$\tilde{v}$}}
\newcommand{\tf}{\mbox{$\tilde{\nabla}$}}
\newcommand{\cK} {\mbox{$\cal K$}}
%\newenvironment{eqnarray*}{eqnarray}

\newcommand{\sech}{\mbox{sech}}

\theoremstyle{break}
%\newtheorem{theorembreak}[theorem]{Theorem}
\newtheorem{algbreak}{Algorithm}

\def\rij{r_{ij}}
\def\rji{r_{ji}}
\def\nij{n_{ij}}
\def\mij{\mu_{ij}}
\def\mji{\mu_{ji}}
\def\hmij{\hat\mu_{ij}}
\def\hmji{\hat\mu_{ji}}
\def\bx{{\bf x}}
\def\bp{{\bf p}}
\def\be{{\bf e}}
\def\bz{{\bf z}}
\def\bv{{\bf v}}
\def\f{{\bf f}}

\def\rf{random forest}
\def\dna{{\sf dna}}
\def\satimage{{\sf satimage}}
\def\waveform{{\sf waveform}}
\def\segment{{\sf segment}}
\def\letter{{\sf letter}}
\def\mnist{{\sf MNIST}}
\def\usps{{\sf USPS}}

\def\rf{random forest}
\def\dnak{{{\sf dna} ($k=3$) }}
\def\satimagek{{{\sf satimage} ($k=6$) }}
\def\waveformk{{{\sf waveform} ($k=3$) }}
\def\segmentk{{{\sf segment} ($k=7$) }}
\def\letterk{{{\sf letter} ($k=26$) }}
\def\mnistk{{{\sf MNIST} ($k=10$) }}
\def\uspsk{{{\sf USPS} ($k=10$) }}

%volumn year pages submitted published author-fullnames
\jmlrheading{5}{2004}{975--1005}{11/03; Revised 05/04}{8/04}{Ting-Fan Wu, Chih-Jen Lin, and Ruby C. Weng}

% Short headings should be running head and authors last names
\ShortHeadings{Probability Estimates for Multi-class Classification by Pairwise Coupling}{Wu, Lin, and Weng}
\firstpageno{975}
\begin{document}
%\setlength{\baselineskip}{32pt}

\title{Probability Estimates for Multi-class Classification by Pairwise Coupling}
\author{\name Ting-Fan Wu \email b89098@csie.ntu.edu.tw \\
{\bf Chih-Jen Lin} \email cjlin@csie.ntu.edu.tw \\
\addr Department of Computer Science \\
National Taiwan University\\
Taipei 106, Taiwan
\AND 
Ruby C. Weng \email chweng@nccu.edu.tw\\
\addr Department of Statistics \\
National Chengchi University \\
Taipei 116, Taiwan}

\editor{Yoram Singer}
\maketitle

\begin{abstract}%   <- trailing '%' for backward compatibility of .sty file
Pairwise coupling is a popular multi-class 
classification method that combines all 
comparisons for each pair of classes. 
This paper presents two approaches for obtaining
class probabilities. 
Both methods can be reduced to linear systems and
are easy to implement. 
We show conceptually and experimentally that
the proposed approaches are
more stable than the two existing popular methods: voting
and the method by \cite{TH98a}. 
\end{abstract}

\begin{keywords}
pairwise coupling, probability estimates, random forest, support vector machines
\end{keywords}

\section{Introduction}
The multi-class classification problem refers to assigning
each of the observations into one of $k$ classes.
As two-class problems are much easier to solve, 
many authors propose to use two-class classifiers for 
multi-class classification.
In this paper we focus on techniques 
that provide a multi-class probability estimate by
combining all pairwise comparisons. 

A common way to combine pairwise comparisons
is by voting \citep{SK90a,JF96a}.
It constructs 
a rule for discriminating 
between every pair of classes and then selecting the class 
with the most winning two-class decisions. 
Though the voting procedure requires just
pairwise decisions, it only predicts a class label. 
In many scenarios, however, probability estimates are desired.
As numerous (pairwise) classifiers do provide class probabilities,
several authors \citep{PR91a,DP95a,TH98a}
have proposed probability estimates by
combining the pairwise class probabilities.


Given the observation $\bx$ and the class label $y$,
we assume that the estimated pairwise class probabilities
$\rij$ of 
$\mij = P(y=i \mid y=i \mbox{ or } j, \bx) $
are available. 
From the $i$th and $j$th classes of a 
training  set, we obtain a model which,
for any new $\bx$,
calculates $\rij$ as an approximation of $\mij$.
Then, using all 
$\rij$, the goal is to estimate 
$p_i=P(y=i\mid \bx), i = 1, \ldots, k$.
In this paper, we first propose a method for
obtaining
probability estimates via an approximation
solution to an identity.
The existence of the solution is guaranteed by
theory in finite Markov Chains.
Motivated by the optimization formulation of
this method, we propose a second
approach.
Interestingly, it can also be regarded as an
improved
version of the coupling approach given by 
\cite{PR91a}.
Both of the proposed methods can be reduced 
to solving linear systems
and are simple in practical implementation.
Furthermore, from conceptual and experimental
points of view, we show that the two proposed methods
are more stable than voting and the method by \cite{TH98a}.

We organize the paper as follows.
In Section \ref{sec:surv-exist-meth}, we review
several existing methods.
Sections \ref{sec:method-linear-system}
and \ref{sec:meth-solv-line} detail the two proposed
approaches.
Section \ref{sec:relations-among-four} presents the relationship between 
different methods through their corresponding optimization
formulas.
In Section \ref{sec:experiments}, we compare these methods
using simulated data.
In Section \ref{real}, we conduct experiments using real data. 
The classifiers considered are support vector machines and random forest. 
A preliminary version of this paper was 
presented previously \citep{TFW03a}.


\section{Survey of Existing Methods}
\label{sec:surv-exist-meth}

For methods surveyed in this section and those proposed later,
each provides a vector of multi-class probability
estimates.
We denote it as $\bp^*$ according to method
$*$.
Similarly, there is an associated rule
$\arg \max_i [p^*_i]$ for prediction
and we denote the rule as 
$\delta_*$.

\subsection{Voting}
\label{sec:voting}
Let $\rij$ be the estimates of $\mij 
\equiv P(y=i \mid y=i \mbox{ or } j, \bx) $
and assume $\rij+\rji=1$.
The voting rule \citep{SK90a,JF96a} is
\begin{equation}
\delta_V=\arg\max_i [\sum_{j: j\neq
i} I_{\{r_{ij} > r_{ji}\}}],
\label{eq:deltaV}
\end{equation}
where $I$ is the indicator function:
$I\{x\}=1$ if $x$ is true, and 0 otherwise.
A simple estimate of probabilities can be derived as 
$$p^v_i=2 \sum_{j:j\neq i} I_{\{r_{ij} > r_{ji}\}}/(k(k-1)).$$

\subsection{Method by Refregier and  Vallet}
\label{sec:using-part-estimated}
With
$\mij = p_i/(p_i+p_j)$,
\cite{PR91a} consider that
\begin{equation}
  \label{eq:choosek}
  \frac{\rij}{\rji}
\approx
\frac{\mij}{\mji}
= 
\frac{p_i}{p_j}.
\end{equation}
Thus, making \eqref{eq:choosek}
an equality may be a way to solve $p_i$. However,
the number of equations, $k(k-1)/2$, is more than 
the number of unknowns $k$, so 
\cite{PR91a} propose to choose any
$k-1$ $\rij$.
Then, with 
the condition $\sum_{i=1}^k p_i = 1$,
$\bp^{RV}$ can be obtained by solving a linear system.
However, as pointed out previously by \cite{DP95a},
the results depend strongly on the selection of 
$k-1$ $\rij$.

In Section \ref{sec:meth-solv-line}, by considering 
\eqref{eq:choosek} as well, we propose a method 
which remedies this problem.

\subsection{Method by Price, Knerr, Personnaz, and Dreyfus}
\label{sec:PKPD}

\cite{DP95a} consider that
\begin{equation*}
\left(  \sum_{j:j\neq i}
P(y=i\mbox{ or } j \mid \bx)
\right)
- (k-2)
P(y=i\mid \bx)
= \sum_{j=1}^k P(y = j
\mid \bx) = 1.
\end{equation*}
Using 
\begin{equation*}
  \rij
\approx 
\mij = 
\frac{P(y=i\mid \bx)}
{P(y=i \mbox{ or } j\mid \bx)},
\end{equation*}
one obtains
\begin{equation}
\label{method2.3}
  p_i^{PKPD} =
\frac{1}{\sum_{j: j \neq i}\frac{1}{\rij}
-(k-2)}.
\end{equation}
As $\sum_{i=1}^k p_i =1$ does not hold, we must
normalize $\bp^{PKPD}$.
This approach is very simple and easy to implement.
In the rest of this paper, we refer to this method as PKPD.

\subsection{Method by Hastie and Tibshirani}
\cite{TH98a} propose to minimize
the Kullback-Leibler (KL) distance between 
$\rij$ and $\mij$:
\label{sec:analys-meth-hast}
\begin{eqnarray}
\label{eq:lp}
  l(\bp)& =& 
\sum_{i \neq j} \nij
\rij
\log\frac{\rij}{\mij},
\\
& = &
\nonumber
\sum_{i < j} \nij
\left(\rij
\log\frac{\rij}{\mij}
+
(1-\rij)
\log\frac{1-\rij}{1-\mij}\right),
\end{eqnarray}
where $\mij = p_i/(p_i+p_j)$,
$\rji=1-\rij$,
and $\nij$ is the number of training data in the
$i$th and $j$th classes. 

To minimize \eqref{eq:lp}, they first calculate
\begin{equation}
\nonumber
    \label{eq:lp'}
    \frac{\partial l(\bp)}{\partial p_i}
 =  
\sum_{j: j \neq i}
n_{ij}
\bigl(-\frac{\rij}{p_i} + 
\frac{1}{p_i + p_j}\bigr).
\end{equation}
Thus, letting
$   \partial l(\bp)/\partial p_i= 0,
i = 1, \ldots, k$ and
multiplying $p_i$ on each term,
\cite{TH98a} propose finding a point that satisfies
\begin{equation}
\label{eq:lp'=0}
  \sum_{j:j \neq i} \nij \mij
= \sum_{j:j \neq i} \nij \rij,
\qquad
\sum_{i=1}^k p_i = 1, \mbox{and}
\qquad
p_i > 0, i = 1, \ldots, k.
\end{equation}
Such a point is obtained by the following algorithm:
\begin{algbreak}
\label{alg:1}
  \begin{enumerate}
  \item \label{item:1}
Start with some initial $p_j > 0, \forall j$ and
corresponding $\mu_{ij}=p_i/(p_i+p_j)$.
\item \label{item:2}
Repeat $(i = 1, \ldots, k, 1, \ldots)$
\begin{eqnarray}
&& \alpha = 
\frac{\sum_{j:j \neq i} \nij \rij}
{\sum_{j:j \neq i} \nij \mij} 
  \label{eq:oldalg0}\\
&& \mij \leftarrow
\frac{\alpha \mij}{\alpha \mij + \mji},
\; \; \mji \leftarrow
1- \mij, \mbox{ for all } j \neq i
\label{eq:oldalg}
\\
&& p_i \leftarrow \alpha p_i
\label{eq:updatepi}
\\
&& normalize \; \bp \; (optional)
\label{eq:optional}
\end{eqnarray}
until $k$ consecutive $\alpha$ are all close to ones.
\item $\bp \leftarrow \bp/\sum_{i=1}^k p_i$
  \end{enumerate}
\end{algbreak}
\eqref{eq:updatepi}
implies that in each iteration, only the $i$th component
is updated and all others remain the same.
There are several remarks about this algorithm.
First, the initial $\bp$ must be positive so that all later
$\bp$ are positive and $\alpha$ is well defined
(i.e., no zero denominator in \eqref{eq:oldalg0}).
Second, \eqref{eq:optional} 
is an optional operation because
whether we normalize $\bp$ or not does not affect the values 
of $\mij$ and $\alpha$ in \eqref{eq:oldalg0} and 
\eqref{eq:oldalg}.

\cite{TH98a} prove that Algorithm \ref{alg:1}
generates a sequence of points at which the KL distance
is strictly decreasing.
However, 
\cite{DRH03a} indicates that
the strict decrease in $l(\bp)$
does not guarantee that any limit point
satisfies \eqref{eq:lp'=0}.
\cite{DRH03a} 
discusses the
convergence of algorithms for generalized Bradley-Terry
models where
Algorithm \ref{alg:1} is a special case.
It points out that \cite{EZ29a} has proved that,
if $\rij >0, \forall i \neq j$,
for any initial point,
the whole sequence
generated by Algorithm \ref{alg:1}
converges to a point satisfying \eqref{eq:lp'=0}.
Furthermore,
this point is the unique global minimum of $l(\bp)$
under the constraints
$\sum_{i=1}^k p_i = 1 \mbox{ and }
 0 \leq p_i \leq 1, i=1, \ldots, k.$

Let $\bp^{HT}$ denote the global minimum of
$l(\bp)$.
It is shown in \cite{EZ29a} and
Theorem 1 of \citep{TH98a} that
if weights $\nij$ in \eqref{eq:lp} are considered equal,
then
$\bp^{HT}$ satisfies
\begin{equation}
  \label{eq:pt}
  p^{HT}_i > p^{HT}_j
\mbox{ if and only if }
\tilde{p}^{HT}_i \equiv 
\frac{2\sum_{s:i\neq s} r_{is}}{k(k-1)}
>
\tilde{p}^{HT}_j \equiv
\frac{2\sum_{s:j\neq s} r_{js}}{k(k-1)}.
\end{equation}
Therefore, $\tilde{\bp}^{HT}$ is sufficient 
if one only requires the classification rule.
In fact, $\tilde{\bp}^{HT}$
can be derived as an approximation to the identity
\begin{equation}
  \label{eq:identity}
p_i = \sum_{j:j\neq i}
\left(\frac{p_i +p_j}{k-1}\right) 
\left(\frac{p_i}{p_i+p_j}\right)
= \sum_{j:j\neq i}
\left(\frac{p_i +p_j}{k-1}\right) \mij
\end{equation}
by replacing 
$p_i+p_j$ with $2/k$,
and $\mij$ with $\rij$ in \eqref{eq:identity}.
We refer to the decision rule as
$\delta_{HT}$, which is essentially
\begin{equation}
\arg\max_i [\tilde{p}^{HT}_i].\label{eq:HT}
\end{equation}
In the next two sections,
we propose two methods which are simpler in both
practical implementation and algorithmic analysis.

If the multi-class data are balanced, it is
reasonable to assume equal weighting (i.e., $\nij=1$) 
as the above.
In the rest of this paper, we restrict our
discussion under such an assumption. 


%A property of this algorithm is that $\bp$
%does not affect the convergence of the procedure.
%This is  why we only have to normalize 
%$\bp$
%in the end.
%However, we can derive the relation
%between $\bp$ and $\mu$ so an equivalent
%algorithm is as follows


%The equivalence comes from the following explanation.
%In Algorithm \ref{alg:1}, after the $t$th loop of 
%$i=1, \ldots, k$
%is done, we use $p^t_i,
%\mij^t$
%to represent values of 
%$p_i$ and $\mij$ at that time.
%If 
%\begin{equation}
%\label{eq:mijt}
%  \mij^t
%= \frac{p^t_i}{p^t_i+p^t_j},
%\mbox{ for all } (i, j) \mbox{ with } i \neq j,
%\end{equation}
%we will show that after $k$ more iterations
%$(i=1, \ldots, k$),
%for all $(i,j)$ with $i \neq j$,
% the new $\mij^{t+1}$
%and $p_i^{t+1}$
%satisfy
%\begin{equation}
%  \mij^{t+1}
%= \frac{p^{t+1}_i}{p^{t+1}_i+p^{t+1}_j}.
%\label{eq:t+1}
%\end{equation}
%In these $k$ iterations, each 
%$\mij$ is updated twice: one at the $i$th and 
%the other at the $j$th iteration. Without loss
%of generality, we assume $i <j$.
%Thus,
%with \eqref{eq:oldalg}, \eqref{eq:oldalg1},
%and
%\eqref{eq:mijt}, $\mij$ and $\mji$ are first updated as
%\begin{eqnarray*}
% &&\mij = 
%\frac{\alpha_i \mij^t}{\alpha_i \mij^t + \mji^t}
%= \frac{p^{t+1}_i}{p^{t+1}_i+p^{t}_j},\\
%&& \mji = 
% \frac{p^{t}_j}{p^{t+1}_i+p^{t}_j}.
%\end{eqnarray*}
%For the second time, we have
%\begin{eqnarray*}
% &&\mji^{t+1} = 
%\frac{\alpha_j \mji}{\alpha_j \mji + \mij}
%= \frac{ \alpha_j p^{t}_j}{\alpha_j p^{t}_j+p^{t+1}_i}
%= \frac{ p^{t+1}_j}{p^{t+1}_i+p^{t+1}_j},\\
%&& \mji^{t+1} = 
%= \frac{ p^{t+1}_i}{p^{t+1}_i+p^{t+1}_j}.
%\end{eqnarray*}
%Therefore, if initially 
%$\mij^1 = \frac{p^1_i}{p^1_i+p^1_j}$,
%by induction, \eqref{eq:t+1}
%always holds. Then Algorithms
%\ref{alg:1} and \ref{alg:2}
%are equivalent.


\section{Our First Approach}
\label{sec:method-linear-system}
As $\delta_{HT}$ relies on $p_i+p_j \approx 2/k$,
in Section \ref{sec:experiments} we use two examples to 
illustrate possible
problems with this rule.
In this section,
instead of replacing $p_i+p_j$ by $2/k$ in 
\eqref{eq:identity},
we propose to solve the
system:
\begin{equation}
p_i = \sum_{j: j \neq i}
(\frac{p_i + p_j}{k-1}) \rij, \forall i,
\label{m0}
\qquad \mbox{subject to } \sum^k_{i=1} p_i=1, p_i \geq 0, \forall i.
\end{equation}
Let $\bp^1$ denote the solution to (\ref{m0}). 
Then the resulting decision rule is
$$\delta_1=\arg\max_i [p^1_i].$$


\subsection{Solving (\ref{m0})}
\label{sec:solving-refm0}
To solve (\ref{m0}), we rewrite it as
\begin{equation}
\label{m1}
Q \bp=\bp, ~~\sum^k_{i=1} p_i=1,~~p_i \geq 0,
\forall i,
\mbox{ where }
  Q_{ij}\equiv
  \begin{cases}
     \rij/(k-1) & \mbox{ if } i \neq j,\\
     \sum_{s: s \neq i} r_{is}/(k-1) & \mbox{ if } i = j.
  \end{cases}
\end{equation}
Observe that $\sum^k_{i=1} Q_{ij}=1$ for $j=1,\ldots,k$
and $0 \leq Q_{ij} \leq 1$ for $i,j=1,\ldots,k$, so there exists
a finite Markov Chain whose transition matrix is $Q$.
Moreover, if $\rij>0$ for all $i \neq j$, then $Q_{ij}>0$,
which implies that this Markov Chain is irreducible and aperiodic.
From Theorem 4.3.3 of \cite{SR96a},
these conditions guarantee the existence of a unique stationary
probability and all states being positive recurrent.
Hence, 
we have the following theorem:
\begin{theorem}
\label{rij-thm}
If $\rij>0$, $i \neq j$, then (\ref{m1}) has a unique 
solution $\bp$ with $0<p_i<1$ $\forall i$. 
\end{theorem}

Assume the solution from Theorem \ref{rij-thm}
is $\bp^*$. We claim that
without the constraints $p_i \geq 0, \forall i$,
the linear system 
\begin{equation}
Q \bp=\bp, \qquad \sum^k_{i=1} p_i=1 \label{eq:k+1}
\end{equation}
still has the same unique solution
$\bp^*$. Otherwise, there is another
solution 
$\bar{\bp}^* (\neq \bp^*)$.
Then for any $0 \leq \lambda \leq 1$, 
$\lambda \bp^* + (1-\lambda)\bar{\bp}^*$
satisfies \eqref{eq:k+1} as well. 
As $p^*_i > 0, \forall i$, when 
$\lambda$ is sufficiently close to 1, 
$\lambda p^*_i + (1-\lambda)\bar{p}^*_i > 0,
i = 1, \ldots, k.$
This violates the uniqueness property in Theorem
\ref{rij-thm}.

Therefore, unlike the method in Section \ref{sec:analys-meth-hast}
where a special iterative procedure has to be implemented,
here we only  solve a simple linear system.
As \eqref{eq:k+1} has $k+1$ equalities but
only $k$ variables, 
practically we remove any one equality from
$Q \bp=\bp$ 
and obtain a square system.
Since the column sum of $Q$ is the vector
of all ones, the removed equality is a linear
combination of all remaining equalities.
Thus, any solution of the square system satisfies
\eqref{eq:k+1} and vice versa.
Therefore,
this square system
has the same unique solution as \eqref{eq:k+1}
and hence can be solved
by standard Gaussian elimination.

Instead of Gaussian elimination,
as the stationary solution of a Markov Chain 
can be derived by the limit of the $n$-step transition
probability matrix $Q^n$, we can solve \eqref{m0}
by repeatedly
multiplying $Q$ with any initial probability vector.


\subsection{Another Look at (\ref{m0})}
\label{another}
The following arguments show that the solution 
to (\ref{m0}) is a global minimum of
a meaningful optimization problem.
To begin, using the property that $\rij+\rji=1$, 
$\forall i \neq j$, we re-express 
$p_i = \sum_{j: j \neq i}
(\frac{p_i + p_j}{k-1}) \rij$
of \eqref{m1}
(i.e., $Q \bp=\bp$ of
\eqref{eq:k+1}) as
\begin{equation}
\nonumber
\label{re-m1}
\sum_{j:j \neq i} \rji p_i -\sum_{j:j \neq i} \rij p_j=0,
i=1,\ldots,k.
\end{equation}
Therefore, a solution of (\ref{m1}) is in fact 
the unique global minimum of the following convex problem:
\begin{eqnarray}
\min_{\bp} & & \sum^k_{i=1} 
(\sum_{j:j \neq i} \rji p_i -\sum_{j:j \neq i} \rij p_j)^2
\nonumber\\
\label{re-m1-op}
\mbox{subject to} & & \sum^k_{i=1} p_i=1, 
\qquad p_i \geq 0, i=1,\ldots,k.
\end{eqnarray}
The reason is that the object function is always nonnegative, 
and it attains zero under (\ref{m1}).
Note that the constraints $p_i \geq 0$ $\forall i$
are not redundant
following the discussion around Equation \eqref{eq:k+1}.

\section{Our Second Approach}
\label{sec:meth-solv-line}


Note that both approaches in 
Sections \ref{sec:analys-meth-hast}
and \ref{sec:method-linear-system} involve 
solving 
optimization problems using relations like
$p_i/(p_i+p_j) \approx \rij$ or 
$\sum_{j:j \neq i} \rji p_i \approx \sum_{j:j \neq i} \rij p_j$.
%This equality implies that for a good estimate $\bp$,
%\begin{equation}
%  \label{eq:rjipi}
%  p_i \approx (p_i + p_j)\rij
%\mbox{ or } \rji p_i \approx \rij p_j.
%\end{equation}
Motivated by (\ref{re-m1-op}), 
we suggest another
optimization formulation as follows:
\begin{equation}
\min_{\bp} 
\sum^k_{i=1} \sum_{j:j \neq i} (\rji p_i - \rij p_j)^2
\qquad \mbox{subject to } 
\sum_{i=1}^k p_i = 1,
\label{eq:method2qp}
p_i \ge 0, \forall i.
\end{equation}
Note that the method \citep{PR91a} described in 
Section \ref{sec:using-part-estimated}
considers a 
random selection of 
$k-1$ equations of the form $\rji p_i = \rij p_j$.
As (\ref{eq:method2qp}) considers all
$\rij p_j -\rji p_i$, not just $k-1$ of them, it
can be viewed as an improved version of 
the coupling approach by \cite{PR91a}. 

Let $\bp^2$ denote the corresponding solution.
We then define the classification rule as
$$\delta_2=\arg\max_i [p^2_i].$$

\subsection{A Linear System from  \eqref{eq:method2qp}}
\label{sec:linear-system-method2}

Since (\ref{re-m1-op}) has a unique solution, which
can be obtained by solving a simple linear system,
it is desirable to see whether the minimization problem 
(\ref{eq:method2qp}) has these nice properties.
In this subsection, we show that
(\ref{eq:method2qp})
has a unique solution and
 can be solved by a simple 
linear system.

First, the following theorem shows that the 
nonnegative 
constraints in \eqref{eq:method2qp} are redundant. 
\begin{theorem}
  Problem \eqref{eq:method2qp} is equivalent to 
  \begin{equation}
\min_{\bp}  \sum_{i=1}^k
\sum_{j: j \neq i}
(\rji p_i - \rij p_j)^2 \qquad
\mbox{\rm subject to} 
\sum_{i=1}^k p_i = 1.
  \label{eq:method2qp2} 
\end{equation}
\label{th:pnonneg}
\end{theorem}
The proof is in Appendix \ref{sec:proof-theor-refth:pn}.
Note that we can rewrite the objective function of
\eqref{eq:method2qp2} as 
\begin{equation}
\min_{\bp} \;\; 2 \bp^T Q \bp
\equiv
\min_{\bp} \;\; \frac{1}{2} \bp^T Q \bp, 
\label{eq:method2qp3}
\end{equation}
where
\begin{equation}
  Q_{ij}=
  \begin{cases}
    \mathop{\sum_{s: s \neq i} r_{si}^2} & \mbox{ if } i = j,\\
-\rji \rij & \mbox{ if } i \neq j.
  \end{cases}
\label{eq:Qij}
\end{equation}
We divide the objective function by a positive factor of four
so its derivative is a simple form 
$Q\bp$.
From \eqref{eq:Qij}, $Q$ is positive 
semi-definite
as for any  $\bv \neq {\bf 0},
\bv^T Q \bv
=1/2\sum_{i=1}^k \sum_{j=1}^k
(\rji v_i - \rij v_j)^2
\geq 0$.
Therefore, without constraints $p_i \geq 0, \forall i$,
(\ref{eq:method2qp3}) is a 
linear-equality-constrained convex quadratic programming
problem. Consequently, a point $\bp$
is a global minimum if and only if it satisfies
the optimality condition:
There is a scalar $b$ such that
\begin{equation}
  \label{eq:method2linearsys}
  \begin{bmatrix}
    Q & \be \\
\be^T & 0
  \end{bmatrix}
  \begin{bmatrix}
    \bp \\ b
  \end{bmatrix}
=
\begin{bmatrix}
  {\bf 0} \\ 1
\end{bmatrix}.
\end{equation}
Here 
$Q\bp$ is the derivative of \eqref{eq:method2qp3},
$b$ is the Lagrangian multiplier of the
equality constraint $\sum_{i = 1}^k p_i = 1$,
$\be$ is the $k \times 1$ vector of all ones, and
${\bf 0}$ is the $k \times 1$ vector of all zeros.
Thus, the solution to \eqref{eq:method2qp} can
be obtained by solving the simple linear system
(\ref{eq:method2linearsys}).


\subsection{Solving \eqref{eq:method2linearsys}}
\label{sec:solvemethod2}

Equation \eqref{eq:method2linearsys}
can be solved by some direct methods in numerical
linear algebra. 
Theorem \ref{th:Qpd}(\ref{item:Q+eetPD}) below shows that 
the matrix in \eqref{eq:method2linearsys}
is invertible; 
therefore, Gaussian elimination can be easily applied.

For symmetric positive definite systems, Cholesky
factorization reduces the time
for Gaussian elimination by half. Though
\eqref{eq:method2linearsys} is symmetric
but not positive definite,
if $Q$ is positive definite,
Cholesky factorization can be used to obtain
$b = -1/(\be^T Q^{-1} \be)$ first
and then
$\bp = -b Q^{-1} \be$.
Theorem \ref{th:Qpd}(\ref{item:Qpdwithassumption}) shows that
$Q$ is positive definite under quite general conditions.
Moreover, even if $Q$ is only positive semi-definite,
Theorem \ref{th:Qpd}(\ref{item:Q+eetPD}) proves that 
$Q+\Delta \be \be^T$ is positive definite for any 
constant $\Delta >0$. 
Along with the fact that \eqref{eq:method2linearsys} is 
equivalent to
\begin{equation}
\label{eq:method2linearsys1}
\nonumber
  \begin{bmatrix}
    Q+\Delta \be \be^T & \be \\
\be^T & 0
  \end{bmatrix}
  \begin{bmatrix}
    \bp \\ b
  \end{bmatrix}
=
\begin{bmatrix}
  \Delta \be \\ 1
\end{bmatrix},
\end{equation}
we can do Cholesky factorization on $Q+\Delta \be \be^T$
and solve $b$ and $\bp$ similarly, regardless
whether $Q$ is positive definite or not.
\begin{theorem}
\label{th:Qpd}
If $r_{tu}>0$ $\forall t \neq u$, we have
\begin{enumerate}[(i)]
\item \label{item:Q+eetPD}
For any $\Delta > 0$, $Q+\Delta \be \be^T$ is 
positive definite.
In addition,
$ \left[ \begin{smallmatrix}
 Q & \be\\ \be^T & 0   
  \end{smallmatrix}
\right]
$ is invertible, and hence
(\ref{eq:method2qp}) has a unique global minimum.
\item \label{item:Qpdwithassumption}
If for any $i = 1,\ldots, k$, there are $s \neq j$ for which $s \neq
  i$, $j \neq i$, and
\begin{equation}
  \frac{r_{si} r_{sj}}{r_{is}}\label{eq:hijneq}
\neq \frac{\rji r_{js}}{\rij},
\end{equation}
then $Q$ is positive definite.
\end{enumerate}
\end{theorem}
We leave the proof in Appendix \ref{sec:proof-theor-refth:qp}.

In addition to direct methods, next 
we propose a simple iterative method for solving
\eqref{eq:method2linearsys}:

\begin{algbreak}
\label{alg:2}
\begin{enumerate}
\item \label{item:2:1}
  Start with some initial $p_i \geq 0, \forall i$ and
  $\sum_{i=1}^k p_i =1$.
\item \label{item:2:2}
  Repeat $(t = 1, \ldots, k, 1, \ldots)$
  \begin{eqnarray}
  && p_t \leftarrow \frac{1}{Q_{tt}}[-\sum_{j:j \neq t}
  Q_{tj} p_j + \bp^T Q \bp] \label{eq:updatept}
  \\
  && normalize \; \bp \; \label{eq:normalizep}
  \end{eqnarray}
until (\ref{eq:method2linearsys}) is satisfied.
\end{enumerate}
\end{algbreak}
Equation \eqref{eq:Qij} and the assumption 
$\rij > 0, \forall i \neq j$, ensure that the
right-hand side of \eqref{eq:updatept} is always
nonnegative.
For
\eqref{eq:normalizep} to be well defined, we must ensure
that
$\sum_{i=1}^k p_i >0$ after the operation in 
\eqref{eq:updatept}. This property holds (see 
\eqref{eq:normalizepn} for more explanation).
With  $b = -\bp^T Q \bp$  obtained from
\eqref{eq:method2linearsys},
\eqref{eq:updatept} is motivated from the 
$t$th equality in \eqref{eq:method2linearsys}
with $b$ replaced by $-\bp^T Q \bp$. 
The convergence of Algorithm \ref{alg:2}
is established in the following theorem:
\begin{theorem}
\label{new}
If
$r_{sj} > 0, \forall s \neq j$,
then
$\{\bp^i\}^{\infty}_{i=1}$, the sequence generated by
Algorithm \ref{alg:2}, converges globally to 
the unique 
minimum of (\ref{eq:method2qp}).
\end{theorem}

The proof is in Appendix \ref{sec:proof-theorem-refnew}.
Algorithm \ref{alg:2}
is implemented in the software \libsvm\ developed
by \cite{CC01a}
for multi-class probability estimates. We discuss some
implementation issues of Algorithm \ref{alg:2}
in Appendix \ref{sec:impl-notes-algor}. 

\section{Relations Between Different Methods}
\label{sec:relations-among-four}
Among the methods discussed in this paper, the four decision rules $\delta_{HT}$, $\delta_1$,
$\delta_2$, and $\delta_V$ can be written as
$\arg\max_i [p_i]$, where $\bp$ is derived by the 
following four optimization formulations
under the constraints $\sum_{i=1}^k p_i = 1$
and  $p_i \ge 0, \forall i$:
\begin{alignat}{2}
\label{del-ht}
& \delta_{HT}: && \min_{\bp} \sum^k_{i=1} 
[\sum_{j:j \neq i} (\rij \frac{1}{k}-\frac{1}{2}p_i)]^2,\\
\label{del-1}
& \delta_1: && \min_{\bp} \sum^k_{i=1} 
[\sum_{j:j \neq i} (\rij p_j- \rji p_i)]^2,\\
\label{del-2}
& \delta_2: && \min_{\bp} \sum^k_{i=1} 
\sum_{j:j \neq i} (\rij p_j- \rji p_i)^2,\\
\label{del-v}
& \delta_V: && \min_{\bp} \sum^k_{i=1} 
\sum_{j:j \neq i} (I_{\{\rij > \rji\}} p_j- 
I_{\{\rji > \rij\}} p_i)^2.
\end{alignat}
Note that (\ref{del-ht}) can be easily verified
from \eqref{eq:pt},
and that (\ref{del-1}) and (\ref{del-2}) have been 
explained
in Sections \ref{sec:method-linear-system} and 
\ref{sec:meth-solv-line}.
For (\ref{del-v}), its solution is
\begin{equation}
\label{eq:pvopt}
p_i=\frac{c}{\sum_{j:j \neq i} I_{\{\rji > \rij\}}},
\end{equation}
where $c$ is the normalizing constant;
and therefore, $\arg\max_i [p_i]$ is the
same as (\ref{eq:deltaV}).\footnote{For
$I_{\{\rij > \rji\}}$ to be well defined, 
we consider $\rij \neq \rji$,
which is generally true. In addition,
if there is an $i$ for which
$\sum_{j:j \neq i} I_{\{\rji > \rij\}}=0$,
an optimal solution of (\ref{del-v}) is $p_i=1$,
and $p_j=0$, $\forall j\neq i$. The resulting
decision is the same as that of (\ref{eq:deltaV}).}
Detailed derivation of \eqref{eq:pvopt}
is in Appendix \ref{sec:deriv-pvopt}.

Clearly, (\ref{del-ht}) can be obtained from
(\ref{del-1}) by letting $p_j = 1/k$ and $\rji = 1/2$.
Such approximations ignore the differences between $p_i$. 
Next, (\ref{del-v}) is from (\ref{del-2}) 
with $\rij$ replaced by $I_{\{\rij > \rji\}}$,
and hence, (\ref{del-v}) may enlarge the differences
between $p_i$.
Moreover, compared with (\ref{del-2}), 
(\ref{del-1}) allows the difference between
$\rij p_j$ and $\rji p_i$ to be canceled first,
so (\ref{del-1}) may tend to underestimate
the differences between $p_i$. 
In conclusion, conceptually, (\ref{del-ht}) and
(\ref{del-v}) are more extreme -- the former tends to
underestimate the differences between
$p_i$, while the latter overestimates them.
These arguments will be supported by simulated and real 
data in the next two sections.

For PKPD approach (\ref{method2.3}),
the decision rule can be written as:
\begin{equation*}
  \delta_{PKPD}
=\arg\min_i [\sum_{j:j\neq i}
\frac{1}{r_{ij}}].
\end{equation*}
This form looks similar to 
$\delta_{HT} =\arg\max_i [\sum_{j:j\neq i} r_{ij}]$,
which can be obtained from (\ref{eq:pt}) and \eqref{eq:HT}. 
Notice that the differences among 
$\sum_{j:j\neq i} r_{ij}$
tend to be larger than those among $\sum_{j:j\neq i} \frac{1}{r_{ij}}$,
because $1/\rij  > 1 > \rij$.
More discussion on these two rules will be given in Section 
\ref{sec:experiments}.


\section{Experiments on Synthetic Data}
\label{sec:experiments}
In this section, we use synthetic data to
compare the performance of existing methods 
described  in Section 2 as well as two new approaches proposed in 
Sections 3 and 4. Here we do not include the method in 
Section \ref{sec:using-part-estimated} because its
results depend strongly on the choice of $k-1$ $r_{ij}$
and our second method is an improved version of it.

\cite{TH98a} design a simple experiment in which 
all $p_i$ are fairly close and their method 
$\delta_{HT}$ outperforms the voting strategy $\delta_V$.
We conduct this experiment first to assess
the performance of our proposed methods.
Following their settings, we define class probabilities
\begin{itemize}
\item [(a)]
$p_1= 1.5/k$, $p_j=(1-p_1)/(k-1)$, $j=2,\dots, k,$
\end{itemize}
and then set
\begin{alignat}{2}
\label{simple-rij-1}
  & r_{ij} = \frac{p_i}{p_i + p_j} + 0.1
z_{ij} && \mbox{ if } i>j,\\
  & r_{ji} = \frac{p_j}{p_i + p_j} + 0.1
z_{ji} = 1- r_{ij} && \mbox{ if } j>i,
\label{simple-rij}
\end{alignat}
where $z_{ij}$ are standard normal variates
and $z_{ji}=-z_{ij}$.
Since $\rij$ are required to be within (0,1),
we truncate $r_{ij}$ at 
$\epsilon$ below and $1-\epsilon$ above,
with $\epsilon = 10^{-7}$.
In this example, class 1 has the highest
probability and hence is the correct
class.

Figure \ref{fig:1a} shows accuracy
rates for each of the five methods 
when $k=2^2,\lceil 2^{2.5}\rceil,$ $ 2^3, \ldots, 2^7$,
where $\lceil x \rceil$ denotes the largest integer not
exceeding $x$.
The accuracy rates are averaged over 1,000 replicates.
Note that in this experiment all classes are 
quite competitive,
so, when using $\delta_V$, sometimes the 
highest vote occurs at two 
or more different classes.
We handle this problem by randomly selecting
one class from the ties. This partly explains the
poor performance of
$\delta_V$. 
Another explanation is that the $\rij$ here are
all close to 1/2, but
(\ref{del-v}) uses 1 or 0 instead, 
as stated in the previous section;
therefore, the solution may be severely biased. 
Besides $\delta_{V}$, the other four rules
have good performance in this example.


Since $\delta_{HT}$ relies on the approximation
$p_i+p_j \approx k/2$, this rule may suffer
some losses if the class probabilities
are not highly balanced. To examine this point,
we consider the following two sets of class
probabilities:
\begin{itemize}
\item[(b)] We let $k_1=k/2$ if $k$ is
even, and $(k+1)/2$ if $k$ is odd; then we define
$p_1=0.95 \times 1.5/k_1$, $p_i=(0.95-p_1)/(k_1-1)$
for $i=2,\ldots,k_1$, and
$p_i=0.05/(k-k_1)$ for $i=k_1+1,\ldots,k$.
\item[(c)] We define 
$p_1=0.95 \times 1.5/2$, $p_2=0.95-p_1$, 
and $p_i=0.05/(k-2), i = 3, \ldots, k$.
\end{itemize}
An illustration of these  three sets of class probabilities
is in Figure \ref{fig:prob}.

\begin{figure}[htbp]
\bigskip
        \begin{center}
          \psset{linewidth=0.3pt} \resizebox{0.8\linewidth}{!}{
                \begin{tabular}{ccccc}
                \input{piea.tex} &
                \mbox{\hspace{1cm}}&
                \input{pieb.tex} &
                \mbox{\hspace{1cm}} &
                \input{piec.tex}\\
                &&&&\\
                (a) & & (b) & & (c)     
                \end{tabular}
              }
            \end{center}
            \caption{Three sets of class probabilities}
        \label{fig:prob}
      \end{figure}

After setting $p_i$, we define the pairwise 
comparisons $\rij$ as in \eqref{simple-rij-1}-(\ref{simple-rij}).
Both experiments are repeated for 1,000 times.
The accuracy rates are shown in Figures
\ref{fig:1b} and \ref{fig:1c}.
In both scenarios, $p_i$ are not balanced.
As expected, $\delta_{HT}$ is quite sensitive
to the imbalance of $p_i$. The situation is
much worse in Figure \ref{fig:sim}(c) because the
approximation $p_i+p_j \approx k/2$ is more
seriously violated, especially when $k$ is large.

A further analysis of Figure \ref{fig:1c} shows that 
when $k$ is large, 
\begin{eqnarray*}
  && r_{12} = \frac{3}{4}+0.1z_{12}, r_{1j} \approx 1+0.1z_{1j},
j\geq 3,\\
  && r_{21} = \frac{1}{4}+0.1z_{21}, r_{2j} \approx 1+0.1z_{2j},
j \geq 3,\\
  && r_{ij} \approx 0 + 0.1z_{ij}, i \neq j, i \geq 3,
\end{eqnarray*}
where $z_{ji}=-z_{ij}$ are standard normal variates.
From (\ref{eq:pt}), the decision rule $\delta_{HT}$ in
this case is mainly on 
comparing $\sum_{j:j\neq 1} r_{1j}$
and $\sum_{j:j\neq 2} r_{2j}$.
The difference between these two sums is
$\frac{1}{2} + 0.1  (\sum_{j:j \neq 1} z_{1j}
-\sum_{j:j \neq 2} z_{2j})$, where
the second term has zero mean and, when
$k$ is large, high variance.
Therefore, for large $k$,
the decision depends strongly on these normal variates,
and the probability of choosing the first class is 
approaching half.
On the other hand, $\delta_{PKPD}$ 
relies on comparing
$\sum_{j:j\neq 1} 1/r_{1j}$
and $\sum_{j:j\neq 2} 1/r_{2j}$.
As the difference between 
$1/r_{12}$ and $1/r_{21}$ is larger than 
that between $r_{12}$ and $r_{21}$, 
though the accuracy rates decline when $k$ increases,
the situation is less serious.

We also analyze the mean square error (MSE) in 
Figure \ref{fig:mse}:
\begin{equation}
\label{mse}
\mbox{MSE} =  \frac{1}{1000}
\sum_{j=1}^{1000}
\frac{1}{k}\sum_{i=1}^k (\hat{p}_i^j - p_i)^2,
\end{equation}
where $\hat{\bp}^j$ is the probability
estimate obtained in the $j$th of the 1,000 replicates.
Overall, $\delta_{HT}$
and $\delta_{V}$ have higher MSE, confirming again that 
they are less stable.
Note that Algorithm \ref{alg:1}
and \eqref{eq:pt} give the same prediction for
$\delta_{HT}$, but their MSE are different.
Here we consider \eqref{eq:pt} 
as it is the one analyzed and compared in 
Section \ref{sec:relations-among-four}.

In summary,  $\delta_1$ and $\delta_2$ are less 
sensitive to $p_i$, and
their overall performance are fairly stable.
All observations about $\delta_{HT}$, $\delta_1$,
$\delta_2$, and $\delta_V$ here agree with our analysis
in Section \ref{sec:relations-among-four}.
Despite some similarity to $\delta_{HT}$,
$\delta_{PKPD}$ outperforms $\delta_{HT}$ in general.
Experiments in this study are conducted using MATLAB.

%        A possible explanation why {\bf markov} gives better results than
%min\_pair is by checking their corresponding optimization formulations.
%For min\_pair, if $\rij$ happens to be larger than $\rij$,
%a property that should not happen when $p_1$ is the largest, then for
%minimizing the term $(\rij p_i - \rji p_j)$, $p_1$ is forced to be 
%small. Therefore, such pairs of $(\rij, \rij)$ easily causes $p_1$
%to be not the largest. On the other hand, for the Markov approach,
%$r_{j1} p_1 - r_{1j} p_j, j\neq 1$, are summed first.
%
%The former minimize the distance between {\bf summations} 
%of $r_{ij}p_i$ and $r_{ji}p_j$ over all $j$, while the later minimize 
%{\bf each} $(r_{ij}p_i - r_{ji}p_j)$ respectively.
%As a result, $p_i$ in the {\bf markov} method is less sensitive to 
%outliers or inconsistence $\rij$, because the other $k-1$ terms in the
%same summation will smooth out the effects from outliers. 
%On the contrary, {\bf min\_pair} method 
%may give wrong estimation of $p_i$,
%due to variation of the corresponding $\rij$.

\afterpage{\clearpage
\begin{figure}[H]
\centering
\begin{tabular}{ccc}
\subfigure[]{\label{fig:1a}\epsfig{figure=fig-sim/figa.eps, width=.305\textwidth,height=.25\textheight}}&
\subfigure[]{\label{fig:1b}\epsfig{figure=fig-sim/figb.eps, width=.305\textwidth,height=.25\textheight}}&                       
\subfigure[]{\label{fig:1c}\epsfig{figure=fig-sim/figc.eps, width=.305\textwidth,height=.25\textheight}}
\end{tabular}
\caption{Accuracy of predicting the true class by
the methods: 
        $\delta_{HT}$ (solid line, cross marked),
        $\delta_V$ (dashed line, square marked),
        $\delta_1$ (dotted line, circle marked),
        $\delta_2$ (dashed line, asterisk marked), 
        and $\delta_{PKPD}$ (dashdot line, diamond marked).}
\label{fig:sim} 
\end{figure}
\begin{figure}[H]
\centering
\begin{tabular}{ccc}
\subfigure[]{\label{fig:mse1a}\epsfig{figure=fig-sim/figamse.eps, width=.305\textwidth,height=.25\textheight}}&
\subfigure[]{\label{fig:mse1b}\epsfig{figure=fig-sim/figbmse.eps, width=.305\textwidth,height=.25\textheight}}&                       
\subfigure[]{\label{fig:mse1c}\epsfig{figure=fig-sim/figcmse.eps, width=.305\textwidth,height=.25\textheight}}
\end{tabular}
\caption{MSE by the methods:
        $\delta_{HT}$ via \eqref{eq:pt} (solid line, cross marked),
        $\delta_V$ (dashed line, square marked),
        $\delta_1$ (dotted line, circle marked),
        $\delta_2$ (dashed line, asterisk marked), 
        and $\delta_{PKPD}$ (dashdot line, diamond marked).}
\label{fig:mse} 
\end{figure}
}
\section{Experiments on Real Data}
\label{real}
In this section we present experimental results on several 
multi-class problems:
\dna,  \satimage, \segment, and \letter\
 from the
Statlog collection \citep{DM94a},
\waveform\ from UCI Machine Learning Repository \citep{CLB98a},
\usps\ \citep{JJH94a}, and \mnist\ \citep{YL98a}.
The numbers of classes and features are reported 
in Table \ref{tbl:datasets}.
Except \dna, which takes two possible values 0 and 1,
each attribute of all other data is linearly 
scaled to $[-1,1]$.
In each scaled data, we randomly select 300 training and 
500 testing instances from thousands of data points. 
20 such selections are generated and the testing error rates
are averaged.
Similarly, we do experiments on larger sets 
(800 training and 1,000 testing).
All training and testing sets used are available at 
\url{http://www.csie.ntu.edu.tw/~cjlin/papers/svmprob/data}
and the code is available at
\url{http://www.csie.ntu.edu.tw/~cjlin/libsvmtools/svmprob/svmprob-1.0.tgz}.

For the implementation of the four probability estimates, 
        $\delta_1$ and $\delta_2$ are via solving linear systems.
        For $\delta_{HT}$, we implement Algorithm \ref{alg:1} with the following
        stopping condition
        $$
        \sum_{i=1}^k \left |\frac{\sum_{j:j\neq i} \rij}{\sum_{j:j \neq i}\mij}-1\right| \leq 10^{-3}.
        $$
We observe that the performance of $\delta_{HT}$ may downgrade
if the stopping condition is too loose.

\begin{table}[H]
\begin{center}
\label{tbl:datasets}
        \begin{tabular}{l|rrrrrrr}
        dataset & \dna & \waveform & \satimage & \segment &
        \usps & \mnist & \letter\\
        \hline 
        \#class & 3 & 3 & 6& 7& 10& 10& 26\\
        \#attribute & 180& 21& 36& 19& 256 & 784& 16
        \end{tabular}
\end{center}
\caption{Data set Statistics}
\end{table}

\subsection{SVM as the Binary Classifier}
\label{svm-binary}
We first consider support vector machines (SVM) 
\citep{BB92a,CC95a}
with the RBF kernel
$e^{-\gamma \|\bx_i - \bx_j\|^2}$ as the 
binary classifier. 
The regularization parameter $C$ 
and the kernel parameter 
$\gamma$ are selected by 
cross-validation (CV).
To begin, for each training set,
a five-fold cross-validation
is conducted on the following points of 
$(C,\gamma)$:
$[2^{-5}, 2^{-3}, \ldots,2^{15}] \times
[2^{-5}, 2^{-3}, \ldots,2^{15}]$.
This is done by modifying \libsvm\ \citep{CC01a},
a library for SVM.  At each $(C,\gamma)$,
sequentially four folds are used as the training set
while one fold as the validation set.
The training of the four folds consists of 
$k(k-1)/2$ binary SVMs.
For the binary SVM of the $i$th and $j$th classes,
we employ an improved implementation 
\citep{HTL03b} of Platt's posterior
probabilities \citep{JP00b} to estimate $\rij$:
\begin{equation}
\label{platt}
\rij=P(i \mid i \mbox{ or } j,\bx)=
\frac{1}{1+e^{A \hat{f} +B}},
\end{equation}
where $A$ and $B$ are estimated by minimizing
the negative log-likelihood function, and
$\hat{f}$ are the decision values of training 
data. \cite{JP00b} and \cite{ZT04a} observe that
SVM decision values are easily clustered at 
$\pm 1$, so the probability 
estimate \eqref{platt}
may be inaccurate.
Thus, it is better to use CV decision values
as we less overfit the model and values are not so
close to $\pm 1$.
In our experiments here, this requires a 
further CV
on the four-fold data (i.e.,
a second level CV).

Next, for each instance in the validation set,
we apply the pairwise coupling methods to obtain
classification decisions. 
The error of the five validation sets is 
thus the cross-validation error 
at $(C,\gamma)$.
From this, each rule obtains its
best $(C,\gamma)$.\footnote{If more than one parameter
sets return the smallest cross-validation error, we simply
choose the one with the smallest $C$.}
Then, the decision values from the five-fold 
cross-validation
at the best $(C,\gamma)$ are employed in (\ref{platt})
to find the final $A$ and $B$ for future use.
These two values and 
the model via applying the best parameters
on the whole training
set
are then used to predict  
testing data.
Figure \ref{fig:cvprocedure}
summarizes the procedure of getting validation accuracy
at each given $(C,\gamma)$.

\begin{figure}[htbp]
\centering
        \begin{pspicture}(-.5,-0.1)(3.3,4.0)
          \rput(1,3.8){\rnode{T}{Given $(C,\gamma)$}}
          \rput(0,3.0){\rnode{L1}{1 fold}} \rput(2,3.0){\rnode{R1}{4 folds}} 

                  \rput(2,2.1){\rnode{R15}{5-fold CV}}
          \rput(2,1.7){\rnode{R2}{decision values}}

          \rput(2,.8){\rnode{R3}{$r_{ij}$}}
          \rput(1,.1){\rnode{B}{validation accuracy}}
                  \psset{arrows=->} 
                  \ncline[nodesep=3pt]{T}{L1} \ncline[nodesep=3pt]{L1}{B}
          \ncline[nodesep=3pt]{T}{R1} \ncline[nodesep=3pt]{R1}{R15} \ncline[nodesep=3pt]{R2}{R3}
          \ncline[nodesep=3pt]{R3}{B}
        \end{pspicture}
      \caption{Parameter selection when using SVM as the binary classifier}
      \label{fig:cvprocedure}
    \end{figure}
The average of 20 MSEs are presented
on the left panel of Figure \ref{fig:svm-rf-mse},
where the solid line represents results of
small sets (300 training/500 testing),
and the dashed line of large sets (800 training/1,000 testing).
The definition of MSE here is similar to (\ref{mse}),
but as there is no correct $p_i$ for these problems, we let $p_i=1$
if the data is in the $i$th class, and 0 otherwise.
This measurement is called Brier Score \citep{GWB50a}, 
which is popular in meteorology.
The figures show that for smaller $k$,
$\delta_{HT}$, $\delta_1$, $\delta_2$ and $\delta_{PKPD}$ 
have similar MSEs, but for larger $k$, $\delta_{HT}$
has the largest MSE.
The MSEs of $\delta_V$ are much larger than those by all other methods, 
so they are not included in the figures.
In summary, the two proposed approaches,
$\delta_1$ and $\delta_2$, are fairly
insensitive to the values of $k$,
and all above observations agree well with previous findings
in Sections \ref{sec:relations-among-four} and 
\ref{sec:experiments}.

Next, left panels of Figures \ref{fig:300-500-accuracy} and 
\ref{fig:800-1000-accuracy}
present the average of 20 test errors
for problems with small size (300 training/500 testing) 
and large size (800 training/1,000 testing), respectively.
The caption of each sub-figure also shows the
average of 20 test errors of the multi-class implementation in
\libsvm.
This rule is voting using merely pairwise SVM decision values,
and is denoted as $\delta_{DV}$ for later discussion.
The figures show that the errors of
the five methods are fairly close for smaller $k$,
but quite different for larger $k$.
Notice that for smaller $k$
(Figures \ref{fig:300-500-accuracy}
and \ref{fig:800-1000-accuracy}
(a), (c), (e), and (g))
the differences of the averaged errors 
among the five methods are small,
and there is no particular trend in these figures. 
However, for problems with larger $k$
(Figures \ref{fig:300-500-accuracy}
and \ref{fig:800-1000-accuracy}
(i), (k), and (m)), the differences
are bigger and  $\delta_{HT}$ is less competitive.
In particular, for \letter\ problem
(Figure \ref{fig:300-500-accuracy} (m), $k=$26),
$\delta_2$ and $\delta_V$ outperform $\delta_{HT}$
by more than 4\%.
The test errors along with MSE seems to indicate 
that, for problems with larger $k$,
the posterior probabilities $p_i$
are closer to the setting of
Figure \ref{fig:1c}, rather than that of 
Figure \ref{fig:1a}.
Another feature consistent with earlier findings
is
that
when $k$ is larger
the results of $\delta_2$ are closer to those of 
$\delta_{V}$, and $\delta_1$ closer to $\delta_{HT}$,
for both small and large
training/testing sets.
As for $\delta_{PKPD}$, its overall performance 
is competitive, but we are not clear about
its relationships to the other methods.

Finally, we consider another criterion on evaluating
the probability estimates: the likelihood.
\begin{equation*}
 \prod_{j=1}^l p^j_{y_j}
\end{equation*}
In practice, we use its log likelihood
and divide the value by
a scaling factor $l$:
\begin{equation}
\frac{1}{l}  \sum_{j=1}^l \log p^j_{y_j},\label{eq:log}
\end{equation}
where $l$ is the number of test data,
$\bp^j$ is the probability estimates for the
$j$th data, and $y_j$ is its actual class label.

A larger value implies a possibly better estimate.
The left panel of Figure \ref{fig:svm-rf-lgl} presents the results
of using SVM as the binary classifier.
Clearly the trend is the same as MSE and accuracy. 
When $k$ is larger, $\delta_2$ and $\delta_V$ 
have larger values and hence better performance.
Similar to MSE, values of $\delta_V$ are 
not presented as they are too small.

%\vspace{1cm}

\subsection{Random Forest as the Binary Classifier}
\label{rf-binary}
In this subsection we 
consider random forest \citep{LB01a} as the binary classifier and
conduct experiments on the same data sets.
As random forest itself can provide multi-class 
probability estimates, we denote the corresponding rule as
$\delta_{RF}$ and also compare it
with the coupling methods.

For each two classes of data, we construct 500 trees as
the random forest classifiers.
Using $m_{\mbox{try}}$ randomly selected features,
a bootstrap sample (around two thirds) of training data
are employed to generate a full tree without pruning.
For each test instance, $r_{ij}$ is simply the proportion
out of the 500 trees that class $i$ wins over class $j$.
As we set the number of trees to be fixed at 500, 
the only parameter left for tuning is $m_{\mbox{try}}$.
Similar to \citep{VS03a}, we select $m_{\mbox{try}}$ 
from $\{1,\sqrt{m},m/3,m/2,m\}$ by five-fold cross
validation, where $m$ is the number of attributes.
The cross validation procedure 
first sequentially uses four folds as the training set
to construct $k(k-1)/2$ pairwise random forests,
next obtains the decision for each instance in the validation set
by the pairwise coupling methods, 
and then calculates the cross validation error at the given 
$m_{\mbox{try}}$ by the error of five validation sets.
This is similar to the procedure
in Section \ref{svm-binary},
but we do not need a second-level CV for
obtaining accurate two-class probabilistic estimates
(i.e., $\rij$).
Instead of CV, a more efficient ``out of bag" validation can 
be used for random forest, but here we
keep using CV
for consistency.
Experiments are conducted using an R-interface
\citep{AL02a} to the code from \citep{LB01a}.

The MSE presented in the right panel of
Figure \ref{fig:svm-rf-mse} shows that
$\delta_1$ and $\delta_2$ yield more stable 
results than $\delta_{HT}$ and $\delta_V$
for both small and large sets.
The right panels of Figures \ref{fig:300-500-accuracy} and 
\ref{fig:800-1000-accuracy} give
the average of 20 test errors. 
The caption of each sub-figure also shows the averaged error 
when using random forest as a multi-class classifier ($\delta_{RF}$).
Notice that random forest bears a resemblance to SVM: 
the errors are only slightly different among the
five methods for smaller $k$, 
but $\delta_V$ and $\delta_2$
tend to outperform $\delta_{HT}$ and $\delta_1$
for larger $k$.
The right panel of Figure \ref{fig:svm-rf-lgl}
presents the log likelihood value \eqref{eq:log}.
The trend is again the same.
In summary, the results by using random forest
as the binary classifier strongly support previous findings
regarding the four methods.

\subsection{Miscellaneous Observations and Discussion}
\label{sec:misc-observ}
Recall that in Section \ref{svm-binary} we
consider $\delta_{DV}$,
which does not use Platt's posterior probabilities. 
Experimental results in Figure \ref{fig:300-500-accuracy}
show that $\delta_{DV}$ 
is quite competitive
(in particular, 3\% better for \letter),
but is about 2\% worse than all
probability-based methods for \waveform. 
Similar observations on \waveform\
are also reported in \citep{KD03a}, where
the comparison is between $\delta_{DV}$ and $\delta_{HT}$.
We explain why the results by probability-based and
decision-value-based methods can be so distinct.
For some problems, the parameters selected by
$\delta_{DV}$ are quite different from those
by the other five rules.
In \waveform, at some 
parameters all
probability-based methods 
gives much higher cross validation accuracy
than $\delta_{DV}$.
We observe, for example, the decision values of 
validation sets are in [0.73, 0.97] and [0.93, 1.02]
for data in two classes;
hence, all data in the validation sets are classified
as in one class and the error is high.
On the contrary, the probability-based methods
fit the decision values by a sigmoid function,
which can better separate the two classes by
cutting at a decision value around 0.95.
This observation shed some light on the difference
between probability-based and decision-value based
methods.

Though the main purpose of this section is to compare
different probability estimates, here we check the
accuracy of another multi-class classification method:
exponential loss-based decoding by
\cite{ELA01a}.
In the pairwise setting, if 
$\hat{f}_{ij} \in R$ is the two-class hypothesis so that
$\hat{f}_{ij} > 0$ ($<0$) predicts the data to be
in the $i$th ($j$th) class, then
\begin{equation}
  \label{eq:ASS}
  \mbox{predicted label}
=
\arg \min_i
\left(
\sum_{j:j<i}
e^{\hat{f}_{ji}}
+
\sum_{j:j>i}
e^{-\hat{f}_{ij}}
\right).
\end{equation}
For SVM, we can simply use decision values
as $\hat{f}_{ij}$.
On the other hand, 
$\rij - 1/2$ is another choice.
Table \ref{tbl:ass-accuracy} presents the error of the seven problem using these two options. Results indicate that 
using decision values is worse than $\rij-1/2$
when $k$ is large
(\usps,\mnist, and \letter).
This observation seems to indicate that
large numerical ranges of $\hat{f}_{ij}$
may cause \eqref{eq:ASS} to have more erroneous
results ($\rij -1/2$ is always in $[-1/2, 1/2]$).
The results of using $\rij -1/2$ is competitive
with those in Figures \ref{fig:300-500-accuracy}
and \ref{fig:800-1000-accuracy} when 
$k$ is small. However, for larger $k$
(e.g., \letter), it is slightly worse than
$\delta_2$ and $\delta_V$.
We think this result is due to the similarity between
\eqref{eq:ASS} and $\delta_{HT}$.
When $\hat{f}_{ij}$ is close to zero,
$e^{\hat{f}_{ij}} \approx 
1 + \hat{f}_{ij}$, so \eqref{eq:ASS}
reduces to a ``linear loss-based encoding.''
When $\rij - 1/2$ is used, 
$\hat{f}_{ji}
= \rji-1/2 = 1/2 - \rij$.
Thus, the linear encoding is 
$
\arg \min_i [\sum_{j: j \neq i}
-\rij]
\equiv
\arg \max_i [\sum_{j: j \neq i}
\rij]$, exactly the same as \eqref{eq:pt}
of $\delta_{HT}$.

\begin{table}[H]
\begin{center}
        \begin{tabular}{@{}l|rrrrrrr@{}}
training/testing ($\hat{f}_{ij}$) & \dna & \waveform & \satimage & \segment & \usps & \mnist & \letter\\
        \hline 
300/500 (dec. values)    &  10.47&  16.23&  14.12&  6.21&  11.57&  14.99&  38.59\\
300/500 ($\rij-1/2$)     &  10.47&  15.11&  14.45&  6.03&  11.08&  13.58&  38.27\\
800/1000 (dec. values)   &  6.36&  14.20&  11.55&  3.35&  8.47&  8.97&  22.54\\
800/1000 ($\rij-1/2$)    &  6.22&  13.45&  11.6&  3.19&  7.71&  7.95&  20.29
        \end{tabular}
\caption{Average of
20 test errors using exponential loss-based decoding (in percentage)}
\label{tbl:ass-accuracy}
\end{center}
\end{table}

%\begin{table}[H]
%\begin{center}
%\caption{Test error using exponential loss-based decoding (in percentage)}
%\label{tbl:ass-accuracy}
%        \begin{tabular}{@{}l|rrrrrrr@{}}
%training/testing ($\hat{f}_{ij}$) & \dna & \waveform & \satimage & \segment & \usps & \mnist & \letter\\
%        \hline 
%300/500 (dec. values)    &12.23 & 15.25 &15.23 &  9.68&  14.65&  16.97 &   42.10\\
%300/500 ($\rij-1/2$)     &10.4  & 14.91 &14.56 &  6.18&  11.84&  14.08 &   40.29\\
%800/1000 (dec. values)   &7.39  & 14.24 &12.51 &  6.78&  10.83&  9.87  &   23.97\\
%800/1000 ($\rij-1/2$)    &6.23  & 13.48 &11.54 &  3.41&  8.07 &  8.13  &   21.70\\
%        \end{tabular}
%\end{center}
%\end{table}

Regarding the accuracy 
of pairwise  (i.e., $\delta_{DV}$) and non-pairwise 
(e.g., ``one-against-the-rest'')
multi-class classification methods,
there are already excellent comparisons.
As $\delta_V$ and $\delta_2$
have similar accuracy to $\delta_{DV}$,
roughly how non-pairwise methods compared
to $\delta_{DV}$ is the same as compared to
$\delta_V$ and $\delta_2$.


The results of random forest
as a multi-class classifier (i.e., 
$\delta_{RF}$) are reported in 
the caption of each sub-figure
in Figures \ref{fig:300-500-accuracy}
and \ref{fig:800-1000-accuracy}.
We observe from the figures that,
when the number of classes is larger, 
using random forest as a multi-class classifier is better
than coupling binary random forests. 
However, for \dna\ ($k=3$) the result is the other way around.
This observation for random forest is left as a future
research issue.

%One note  is that some papers (e.g. \citep{BZ01a})
%consider a different algorithm from Algorithm 
%\ref{alg:1} where
%the main loop is changed as follows:
%\begin{algbreak}
%\label{alg:2}

%\begin{enumerate}
%\item[\ref{item:2}'] Repeat

% \qquad \qquad for $i = 1, \ldots, k$
%\begin{eqnarray*}
%&& \alpha_i = 
%\frac{\sum_{j:j \neq i} \nij \rij}
%{\sum_{j:j \neq i} \nij \mij} \\  
%&& p_i \leftarrow \alpha_i p_i.
%\end{eqnarray*}
%\qquad \qquad normalize  $\bp$

%\qquad \qquad for all $(i,j)$ with $i \neq j$
%\begin{equation*}
%  \mij \leftarrow \frac{p_i}{p_i+ p_j}
%\end{equation*}
%until $\alpha_i$ is close to zero for all $i = 1, \ldots, k$.
%\end{enumerate}
%\end{algbreak}
%Therefore, $\mij$
%are updated only after all $p_i, i = 1,\ldots, k$
%are changed.
%However,  the convergence of
%this algorithm is not clear yet.

\acks{The authors thank S. Sathiya Keerthi for helpful comments.
They also thank the associate editor and anonymous referees for useful comments.
This work was supported in part by
the National Science Council of Taiwan via the grants
NSC 90-2213-E-002-111
and
NSC 91-2118-M-004-003.}

%---------------------bibliography


\appendix
%\section*{Appendix}
%Throughout the appendix, we denote 
%$T=\{1,2,\ldots\}$ and $T_i=\{i+nk:n=0,1,\ldots\}$ for $i=1,\ldots,k$.

%\section{Proof of Theorem \ref{th:convhastie}}
%\label{sec:proof-theor-refth:c}

%We first need to prove three lemmas:
%\begin{lemma}
%If $p^1_i > 0, i = 1, \ldots, k$,
%then
%$\{p^t_i\}^{\infty}_{t=1},
%i = 1, \ldots, k$ are bounded below from zero.
%\label{lemma:notzero}
%\end{lemma}
%\Proof
%If the result does not hold, there is an index $\bar{i}$
%and an infinite index set $T' \subset T$ such that
%\begin{equation*}
%\lim_{t\in T', t\rightarrow \infty}
%p^t_{\bar{i}}=0.
%\end{equation*}
%Since
%$\sum_{i=1}^k p^t_i = 1,
%\forall t$ and $k$ is finite, we can find an infinite
%subset $T^{''}$ of $T'$ and an index  $\bar{j}$ such that
%\begin{equation*}
%p^t_{\bar{j}}
%\geq \epsilon > 0,
%\forall t \in T^{''},
%\end{equation*}
%where $\epsilon$ is a fixed number.
%Then, with \eqref{eq:lpineqjform}
%and $\rij > 0, \forall i \neq j$,
%\begin{eqnarray*}
%  l(\bp^t) & = & 
%\mathop{\sum_{(i,j): i \neq j}}
%\nij \rij \log \frac{\rij}{p^t_i/(p^t_i+p^t_j)} \\
%& \geq & 
%\bigl(\mathop{\sum_{\substack{ 
%    (i,j): i \neq j\\ i \neq \bar{i}, j \neq \bar{j}
%}}
%}
%\nij \rij \log \rij
%\bigr)
%+ n_{_{\bar{i}\bar{j}}} r_{\bar{i}\bar{j}}
%\log \frac{r_{\bar{i}\bar{j}}}{p^t_{\bar{i}}/(p^t_{\bar{i}}+
%p^t_{\bar{j}})} \\
%& \rightarrow & \infty
%\mbox{ when } t \in T^{''}, t \rightarrow \infty.
%\end{eqnarray*}
%This contradicts \eqref{eq:decrease},
%which  implies $l(\bp^t)$
%is bounded above by $l(\bp^1)$.
%\qed

%\begin{lemma}
%Let $\alpha_t$ be as in (\ref{alpha_t}).
%Then, $\{\alpha^t\}_{t=1}^{\infty}$
%globally converges to 1.
%\label{lemma:alpha=1}
%\end{lemma}
%\Proof
%If this is not true, there is an $\epsilon > 0$
%and an infinite index set $T' \subset T$
%such that
%\begin{equation}
%  \label{eq:>1}
%  |\alpha^t-1| > \epsilon, \qquad \forall t \in T'.
%\end{equation}
%As $T'$ is an infinite set, there is an index $i$
%so that $T'\cap T_i$ is infinite.
%Since
%$\{\bp^t\}_{t=1}^{\infty}$
%are in a compact set,
%we can find an infinite subset $T^{''}$ of $T'\cap T_i$
%and $\bp^*$
%such that
%\begin{equation}
%\label{eq:bptconv}
%\lim_{t \in T^{''}, t \rightarrow \infty}
%\bp^t = \bp^*,
%\end{equation}
%where by the decreasing property of $l$,
%\begin{equation}
%l(\bp^*) <  l(\bp^t), ~\forall t.
%\label{add}
%\end{equation}
%Therefore, for all $t \in T^{''}$, the working index is 
%always $i$.
%Lemma \ref{lemma:notzero}
%implies $p^*_j > 0,
%j=1, \ldots, k$, so
%the following limit exists and is defined as 
%$\alpha^*$:
%\begin{equation}
%  \alpha^*
%\equiv \lim_{t\in T^{''}, t\rightarrow \infty}
%\alpha^t
%=
%\lim_{t\in T^{''}_i, t\rightarrow \infty}
%\frac{\sum_{j:j \neq i} \nij \rij}
%{\sum_{j:j \neq i} \nij \frac{p^t_i}{p^t_i + p^t_j}}
%=
%\frac{\sum_{j:j \neq i} \nij \rij}
%{\sum_{j:j \neq i} \nij \frac{p^*_i}{p^*_i + p^*_j}}.
%\label{eq:alpha*i}
%\end{equation}

%Next we define
%\begin{equation*}
%\bp^{*,i}
%\equiv[ p^{*}_1 , \ldots,
%p^*_{i-1},
%\alpha^* p^{*}_{i},
%p^{*}_{i+1}, \ldots,  
%p^{*}_k]^T. 
%\end{equation*}
%Using the explanation in \eqref{eq:pt+1beforenormalization},
%if we use 
%$\bp^*$ in the beginning of one iteration
%in Algorithm \ref{alg:1} (with the working index $i$),
%then
%$\bp^{*,i}$
%is the solution of the next iteration
%(without normalization).
%From \eqref{eq:>1} and \eqref{eq:alpha*i}
%$\alpha^* \neq 1$.
%Together with \eqref{eq:decrease}, we obtain
%\begin{equation}
%  \label{eq:p*idecrease}
%l(\bp^{*,i})
%<l(\bp^{*}).
%\end{equation}

%Similar to the definition of $\bp^{*,i}$,
%for each $t \in T'_i$, we define
%\begin{equation}
%\hat{\bp}^{t+1}
%\equiv 
%[p^{t}_1 , \ldots,
%p^{t}_{i-1},
%\alpha^t p^{t}_{i},
%p^{t}_{i+1}, \ldots,  
%p^{t}_k]^T,
%\label{eq:hatpt+1}
%\end{equation}
%from this we can calculate the value after normalization
%$\bp^{t+1} = \hat{\bp}^{t+1}/ \|\hat{\bp}^{t+1}\|$.
%From \eqref{eq:bptconv} and \eqref{eq:alpha*i},
%\begin{equation}
%  \label{eq:pt+1}
%\lim_{t \in T^{''}, t 
%\rightarrow \infty} \hat{\bp}^{t+1}
%= \bp^{*,i}.
%\end{equation}
%From \eqref{eq:p*idecrease}, \eqref{eq:pt+1}, 
%and the continuity of $l(\bp)$, 
%\begin{equation}
%\lim_{t\in T^{''}, t\rightarrow \infty}
%l(\hat{\bp}^{t+1})
%=l(\bp^{*,i}) < l(\bp^*).
%\label{eq:pt+1limit}
%\end{equation}
%Recall that $l(\hat{\bp}^t)$ is strictly decreasing in $t$,
%so \eqref{eq:pt+1limit} implies
%$l(\hat{\bp}^{t+1}) < l(\bp^*)$ for $t$ sufficiently large.
%However, by (\ref{add}) and (\ref{lp-alpha}),
%\begin{equation*}
%l(\bp^*) <  l(\bp^{t+1}) = l(\hat{\bp}^{t+1}),
%\end{equation*}
%contradicting \eqref{eq:pt+1limit}.
%\qed
%\begin{lemma}
%  For any 
%infinite set $T'_i \subset T_i$, $i \in \{1,\ldots,k\}$,
%if 
%\begin{equation*}
%  \lim_{t \in T'_i, t\rightarrow \infty}
%\bp^t = \bp^*,
%\end{equation*}
%then
%\begin{equation}
%  \label{eq:iandi+1}
%  \lim_{t \in T'_{i}, t\rightarrow \infty}
%\bp^{t+1} = \bp^*
%\mbox{ and }
%\frac{\sum_{j:j \neq i+1} n_{i+1,j} r_{i+1,j}}
%{\sum_{j:j \neq i+1} n_{i+1,j} \frac{p^*_{i+1}}{p^*_{i+1} + p^*_j}}
%=1.
%\end{equation}
%If $i=k$, then the same result holds with
%$i+1$ in the second equality replaced by 1.
%\label{lemma:iandi+1}
%\end{lemma}
%\Proof
%Since 
%$i$ is the working index for all $t \in T'_i$,
%the next iteration 
%$\bp^{t+1}, t \in T'_i$
%and its values before normalization
%$\hat{\bp}^{t+1}$
%can be defined 
%similar to \eqref{eq:hatpt+1}.
%With $\lim_{t\in T'_i, t\rightarrow \infty}
%\alpha^t=1$ by 
%Lemma \ref{lemma:alpha=1},
%the first equation of \eqref{eq:iandi+1} follows from
%a similar derivation in 
%\eqref{eq:pt+1}.

%For the second result, using Lemma \ref{lemma:alpha=1},
%\begin{eqnarray}
%1 & = & 
%\lim_{t \in T'_i, t\rightarrow \infty}
%\alpha^{t+1} \nonumber \\
%& = & 
%\lim_{t \in T'_i, t\rightarrow \infty}
%\frac{\sum_{j:j \neq i+1} n_{i+1,j} r_{i+1,j}}
%{\sum_{j:j \neq i+1} n_{i+1,j} \frac{p^{t+1}_{i+1}}
%{p^{t+1}_{i+1} + p^{t+1}_j}}
%\nonumber\\
%& = & 
%\frac{\sum_{j:j \neq i+1} n_{i+1,j} r_{i+1,j}}
%{
%\sum_{j:j \neq i+1} n_{i+1,j} \frac{p^*_{i+1}}{p^*_{i+1} + p^*_j}
%},
%  \label{eq:i+1limit}
%\end{eqnarray}
%where the last line follows from the first equation
%of \eqref{eq:iandi+1}. 
%\qed


%To prove the theorem, we consider any convergent subsequence
%$\{\bp^t\}_{t\in T'}$, where
%$T'$ is an infinite subset of $T$
%and $\bp^t \rightarrow \bp^*$
%as $t\in T', t\rightarrow \infty$.
%Since $k$ is finite, there is an 
%index $i$ such that 
%$T'_i \equiv
%T'\cap T_i$ is also an infinite set
%and
%\begin{equation*}
%  \lim_{t \in T'_i,  t\rightarrow \infty}
%\bp^t = \bp^*.
%\end{equation*}
%Then, by Lemma \ref{lemma:iandi+1},
%\begin{equation*}
%  \lim_{t \in T'_{i+1},  t\rightarrow \infty}
%\bp^t = \bp^*,
%\end{equation*}
%where $T'_{i+1} \equiv \{ t \mid t-1 \in T'_i\}$.
%If $i=k$, we simply replace $T'_{i+1}$ by $T'_1$. 
%Repeatedly apply this Lemma, we obtain
%\begin{equation*}
%\frac{\sum_{j:j \neq i} n_{ij} r_{ij}}
%{\sum_{j:j \neq i} n_{ij} \frac{p^*_{i}}{p^*_{i} + p^*_j}}
%=1,\forall i = 1, \ldots, k.
%\end{equation*}
%With $\sum_{i=1}^k p^*_i = 1$
%and $p^*_i > 0$ from Lemma \ref{lemma:notzero},
%$\bp^*$ satisfies \eqref{eq:lp'=0}, so
%the proof is complete.

%\section{Proof of Theorem \ref{th:localmin}}
%\label{sec:proof-localmin}

%For a point $\bp$ satisfying 
%\eqref{eq:lp'=0}, 
%$\nabla l(\bp)=0$.
%Thus, $(1/p_i)(\partial l(p)/\partial p_i)=0$
%for $i=1,\ldots,k$, which by \eqref{eq:lp'} can be rewritten as
%\begin{equation*}
%\sum_{j: j \neq i}
%\nij (\frac{\rij}{p_i^2})
%= 
%\sum_{j: j \neq i}
%\frac{\nij}{p_i(p_i+p_j)}.  
%\end{equation*}
%Using this equality and \eqref{eq:hessian},
%\begin{equation*}
%  \frac{\partial^2 l(\bp)}{\partial p_i
%\partial p_i}
%= \sum_{j: j \neq i}
%\nij \frac{p_j}{p_i(p_i+p_j)^2},~i=1,\ldots,k.
%\end{equation*}
%Then, at this $\bp$,
%\eqref{eq:hessian} is 
%positive semi-definite as for any 
%vector $\bv=[v_1,\ldots,v_k]^T \neq 0$,
%\begin{eqnarray}
%\nonumber
% \sum_{i=1}^k \sum_{j=1}^k v_i 
%  \frac{\partial^2 l(\bp)}{\partial p_i
%\partial p_j}
%v_j 
%& = & 
%\sum_{(i,j): i \neq j}
%v_i \frac{-\nij}{(p_i+p_j)^2}v_j
%+ \sum_{i=1}^k \sum_{j:j\neq i}
%\frac{\nij p_j}{p_i (p_i+p_j)^2}v_i^2
%\\
%&= &
%\sum_{(i,j):i< j}
%\nij 
%(\frac{\sqrt{p_j}}{\sqrt{p_i}(p_i+p_j)}
%v_i
%- \frac{\sqrt{p_i}}{\sqrt{p_j}(p_i+p_j)}
%v_j)^2 \geq 0,
%  \label{eq:vhessianv}
%\end{eqnarray}
%where equality holds if and only if 
%\begin{equation}
%  \label{eq:vandp}
%  \frac{v_1}{p_1}
%= \cdots
%= \frac{v_k}{p_k}.
%\end{equation}
%Suppose $\bp$
%is not a local minimum of \eqref{eq:minlp}. Then,
%there is a sequence of vectors $\{\bv^i\}_{i=1}^{\infty}$
%such that
%\begin{equation}
%\bv^i \neq 0, \forall i, 
%\qquad \sum_{j=1}^k v^i_j = 0,\forall i,
%\qquad \lim_{i \rightarrow \infty} \bv^i = 0, 
%\qquad
%l(\bp+\bv^i)
%< l(\bp), \forall i.\label{eq:limvi}
%\end{equation}
%For any $\bv$,
%\begin{eqnarray*}
%  l(\bp+\bv) & = &
%l(\bp)
%+ \nabla l(\bp)^T \bv
%+ \frac{1}{2}\bv^T \nabla^2 l(\bp)
%\bv
%+ 
%\frac{1}{6}
%\sum_{j,s,r}
% v_j v_s v_r
%\frac{\partial^3 l(\bz)}{
%\partial p_j \partial p_s \partial p_r}\\
%& = & 
%l(\bp)
%+ \frac{1}{2}\bv^T \nabla^2 l(\bp)
%\bv
%+ 
%\frac{1}{6}
%\sum_{j,s,r}
% v_j v_s v_r
%\frac{\partial^3 l(\bz)}{
%\partial p_j \partial p_s \partial p_r},
%\end{eqnarray*}
%where $\bz$ is on the line segment between $\bp$ and $\bp+\bv$.
%Thus, \eqref{eq:limvi}
%implies
%\begin{equation}
%  \label{eq:v/v}
% \frac{1}{2}(\frac{\bv^i}{\|\bv^i\|})^T
% \nabla^2 l(\bp)
%(\frac{\bv^i}{\|\bv^i\|})
%+ 
%\frac{1}{6}
%\sum_{j,s,r}
%\frac{v_j^i v_s^i v_r^i}{||\bv^i||^2}
%\frac{\partial^3 l(\bz^i)}{
%\partial p_j \partial p_s \partial p_r}
% < 0,
%\end{equation}
%where $\bz^i$ is on the segment between $\bp$ and $\bp+\bv^i$.

%Since $\{\frac{\bv^i}{\|\bv^i\|}
%\}_{i=1}^{\infty}$ are in a compact region, there is a subsequence
%which converges to
%$\bv^*$ with
%\begin{equation}
%  \label{eq:v*}
%  \|\bv^*\|=1 \mbox{ and }
%\sum_{j=1}^k \bv^*_j = 0,
%\end{equation}
%where the second equality comes from \eqref{eq:limvi}.
%Taking the limit of \eqref{eq:v/v},
%we have
%\begin{equation}
%  (\bv^*)^T \nabla^2 l(\bp) \bv^*
%\leq 0.
%\label{eq:>0}
%\end{equation}
%However, \eqref{eq:v*} implies that 
%\eqref{eq:vandp} does not hold, so by \eqref{eq:vhessianv},
%$ (\bv^*)^T \nabla^2 l(\bp) \bv^* > 0$.
%This contradicts \eqref{eq:>0}.
%Therefore, $\bp$ is a local
%minimum of \eqref{eq:lp}, so the proof is complete.

%\section{Proof of Theorem \ref{th:pilargest}}
%\label{sec:pilargest}

%Define
%\begin{equation*}
%  j \equiv \arg\max_l \{p_l\}.
%\end{equation*}
%If this theorem is not true, $j \neq i$ and
%$p_j > p_i$.
%Then,
%\begin{eqnarray}
%  0 & =&
%(  \sum_{l: l \neq i} r_{il}-1)p_i
%+ \sum_{l: l \neq i} r_{il} p_l \nonumber \\
%& \geq & 
%(  \sum_{l: l \neq i} r_{il}-1)p_i
%+  r_{ij} p_j \nonumber \\
%& > & 
%(  \sum_{l: l \neq i} r_{il}-1)p_j
%+  r_{ij} p_j \label{eq:sumril}\\
%& \geq & 
%(2 \sum_{l: l \neq j} r_{jl} -1)p_j \nonumber \\
%& = & 
%\sum_{l: l \neq j} r_{jl} p_j+
%( \sum_{l: l \neq j} r_{jl} -1)p_j \nonumber \\
%& \geq & \sum_{l: l \neq j} r_{jl} p_l+
%( \sum_{l: l \neq j} r_{jl} -1)p_j = 0, \nonumber
%\end{eqnarray}
%a contradiction.
%Here, the strict inequality \eqref{eq:sumril}
%uses the property
%$ \sum_{l: l \neq i} r_{il}-1 < 0$,
%which comes from the assumption that $\rij > 0, \forall i\neq j$.

\section{Proof of Theorem \ref{th:pnonneg}}
\label{sec:proof-theor-refth:pn}

It suffices to prove that any optimal
solution $\bp$ of \eqref{eq:method2qp2} 
satisfies 
$p_i \ge 0, 
i = 1, \ldots, k$.
If this is not true, without loss of 
generality, we assume
\begin{equation*}
  p_1 < 0, \ldots, p_r < 0, 
p_{r+1} \ge 0, \ldots, p_k \ge 0,
\end{equation*}
where $r <k$ because $\sum^k_{i=1}p_i=1$.
We can then define a new feasible solution of
\eqref{eq:method2qp2}:
\begin{equation*}
  p'_1 = 0, \ldots, p'_r = 0, 
p'_{r+1} = p_{r+1}/\alpha , \ldots, p'_k
= p_k/\alpha,
\end{equation*}
where $\alpha = 1 - \sum^r_{i=1} p_i > 1$.

With $\rij > 0$ and $\rji > 0$, we obtain
\begin{eqnarray}
\nonumber
  \label{eq:3cond1}
  && (\rji p_i - \rij p_j)^2 \geq 0
= (\rji p'_i - \rij p'_j)^2,
\mbox{ if } 1 \le i, j \le r,\\
\nonumber
&& (\rji p_i - \rij p_j)^2 > 
\frac{(\rij p_j)^2}{\alpha^2}
= (\rji p'_i - \rij p'_j)^2,
\mbox{ if } 1 \le i \le r, r+1\le j \le k,
\label{eq:3cond2}\\
\nonumber
&& (\rji p_i - \rij p_j)^2 \geq 
\frac{(\rji p_i -\rij p_j)^2}{\alpha^2}
= (\rji p'_i - \rij p'_j)^2,
\mbox{ if }  r+1\le i, j \le k.
\label{eq:3cond3}
\end{eqnarray}
Therefore,
\begin{equation*}
\sum^k_{i=1} \sum_{j: j \neq i}  (\rij p_i - \rji p_j)^2 >
\sum^k_{i=1} \sum_{j: j \neq i}  (\rij p'_i - \rji p'_j)^2. 
\end{equation*}
This contradicts the assumption that
$\bp$ is an optimal solution of \eqref{eq:method2qp2}.

\section{Proof of Theorem \ref{th:Qpd}}
\label{sec:proof-theor-refth:qp}
\begin{enumerate}[(i)]
\item
  If $Q+\Delta \be \be^T$ 
is not positive definite, there is a vector $\bv$ with $v_i
  \neq 0$ such that
\begin{equation}
\label{eq:sumvt}
  \bv^T(Q+\Delta \be \be^T) \bv
= \frac{1}{2}\sum^k_{t=1} \sum_{u: u \neq t} (r_{ut} v_t - r_{tu} v_u)^2 
+ \Delta (\sum_{t=1}^k v_t)^2 = 0.
\end{equation}
For all $t\neq i$, $r_{it}v_t - r_{ti}v_i =0$, so
\begin{equation*}
v_t = \frac{r_{ti}}{r_{it}} v_i
\neq 0.
\end{equation*}
Thus,
\begin{equation*}
  \sum_{t=1}^k v_t
= (1+\sum_{t: t \neq i}
\frac{r_{ti}}{r_{it}})v_i \neq 0,
\end{equation*}
which contradicts \eqref{eq:sumvt}.

The positive definiteness of $Q+\Delta \be \be^T$
implies that 
$
\left[
  \begin{smallmatrix}
    Q + \Delta \be \be^T & \be \\
\be^T & 0
  \end{smallmatrix}
\right]
$
is invertible.
As $
\left[
  \begin{smallmatrix}
    Q + \Delta \be \be^T & \be \\
\be^T & 0
  \end{smallmatrix}
\right]
$ is from adding 
the last row of 
$ \left[ \begin{smallmatrix}
    Q & \be \\
\be^T & 0
  \end{smallmatrix}
\right]
$ to its first $k$ rows
(with a scaling factor $\Delta$),
the two matrices have the same rank.
Thus, 
$ \left[ \begin{smallmatrix}
    Q & \be \\
\be^T & 0
  \end{smallmatrix}
\right]
$
is invertible as well.
Then
\eqref{eq:method2linearsys} 
has a unique solution, and so does 
\eqref{eq:method2qp}.

\item
  If $Q$ is not positive definite, there is a vector $\bv$ with $v_i
  \neq 0$ such that
\begin{equation*}
  \bv^TQ \bv
= 
\frac{1}{2}
\sum^k_{t=1} \sum_{u: u \neq t} (r_{ut} v_t - r_{tu} v_u)^2 = 0.
\end{equation*}
Therefore,
\begin{equation*}
(r_{ut} v_t - r_{tu} v_u)^2 = 0,
\forall t \neq u.
\end{equation*}
As $r_{tu} > 0, \forall t \neq u$, for any $s \neq j$ for which $s
\neq i$ and $j \neq i$, we have
\begin{equation}
\label{vsj}
  v_s = \frac{r_{si}}{r_{is}} v_i,\qquad 
  v_j = \frac{r_{ji}}{r_{ij}} v_i,\qquad 
  v_s = \frac{r_{sj}}{r_{js}} v_j.
\end{equation}
As $v_i \neq 0$, (\ref{vsj}) implies
\begin{equation*}
  \frac{r_{si} r_{sj}}{r_{is}}
= \frac{\rji r_{js}}{\rij},
  \end{equation*}
  which contradicts \eqref{eq:hijneq}.
\end{enumerate}
\section{Proof of Theorem \ref{new}}
\label{sec:proof-theorem-refnew}
First we need a lemma to show the strict decrease 
of the objective function:

\begin{lemma}
\label{new:lemma1}
If $\rij>0, \forall i \neq j$,
$\bp$ and $\bp^n$ are from two consecutive iterations of 
Algorithm \ref{alg:2}, and $\bp^n \neq \bp$,
then
\begin{equation}
\label{pnt}
\frac{1}{2}(\bp^{n})^T Q \bp^{n} < \frac{1}{2}\bp^T Q \bp.
\end{equation}
\end{lemma}
\Proof
Assume that $p_t$ is the component to be updated. 
Then, $\bp^n$ is obtained through the following calculation:
\begin{equation}
  \bar{p}_i=
  \begin{cases}
   p_i & \mbox{ if } i \neq t,\\
   \frac{1}{Q_{tt}}(-\sum_{j:j \neq t} Q_{tj} p_j + \bp^T Q \bp) & \mbox{ if }
   i =t,
  \end{cases}
\label{eq:pbari}
\end{equation}
and 
\begin{equation}
\bp^n =\frac{\bar{\bp}}{\sum_{i=1}^k \bar{p}_i}.\label{eq:normalizepn}
\end{equation}
For \eqref{eq:normalizepn} to be a valid operation,
$\sum_{i=1}^k \bar{p}_i$ must be strictly positive.
To show this, we first suppose that
the current solution $\bp$ satisfies 
$p_i \geq 0, i = 1, \ldots, l$,
but the next solution 
$\bar{\bp}$ has
$\sum_{i=1}^k \bar{p}_i=0$.
In Section \ref{sec:solvemethod2},
we have shown that 
$\bar{p}_t \geq 0$,
so with $\bar{p}_i = p_i\geq 0,
\forall i \neq t$, 
$\bar{p}_i=0$ for all $i$.
Next, from (\ref{eq:pbari}),
$p_i=\bar{p}_i=0$ for $i \neq t$,
which, together with the equality $\sum^k_{i=1} p_i=1$ implies
that $p_t=1$.
However, if $p_t=1$ and $p_i =0$ for $i \neq t$, then
$\bar{p}_t = 1$ from (\ref{eq:pbari}).
This contradicts the situation that 
$\bar{p}_i=0$ for all $i$.
Therefore, by induction, the only requirement is
to have nonnegative initial $\bp$.

To prove (\ref{pnt}),
first we rewrite
the update rule 
\eqref{eq:pbari}
as 
\begin{eqnarray}
\label{eq:pbari1}
\bar{p}_t & = & 
p_t + 
   \frac{1}{Q_{tt}}(-(Q\bp)_t + \bp^T Q \bp) \\
&=& p_t + \Delta. \nonumber
\end{eqnarray}
Since we keep 
$\sum_{i=1}^k p_i = 1$,
$\sum_{i=1}^k \bar{p}_i
= 1 + \Delta$. Then
\begin{eqnarray}
&& \bar{\bp}^T Q \bar{\bp} - (\sum_{i=1}^k \bar{p}_i)^2 \bp^T Q \bp 
\nonumber\\
&=& \bp^T Q \bp + 2 \Delta (Q\bp)_t 
+ Q_{tt} \Delta^2 -
(1+\Delta)^2 
\bp^T Q \bp
\nonumber\\
&=&  2 \Delta (Q\bp)_t 
+ Q_{tt} \Delta^2 -
(2\Delta + \Delta^2) 
\bp^T Q \bp \nonumber \\
&=&  \Delta \left( 2 (Q\bp)_t 
-2 \bp^T Q \bp
+ Q_{tt} \Delta -
\Delta  
\bp^T Q \bp\right) \nonumber \\
&=&  \Delta( -Q_{tt} \Delta -
\Delta  
\bp^T Q \bp) \label{eq:pbarT1}\\
& = & -\Delta^2(Q_{tt} +
\bp^T Q \bp) < 0.
\label{pbarT}
\end{eqnarray}
\eqref{eq:pbarT1} follows from the definition of 
$\Delta$ in \eqref{eq:pbari1}.
For \eqref{pbarT},
it uses 
$Q_{tt} = \sum_{j:j \neq t} r^2_{jt} >0$
and 
$\Delta \neq 0$,
which comes from the assumption
$\bp^n \neq \bp$.
\qed

Now we are ready to prove the theorem.
If this result does not hold, there is a
convergent sub-sequence $\{\bp^i\}_{i \in K}$
such that $\bp^*= \mbox{lim}_{i \in K, i \rightarrow \infty}
\bp^i$ is not optimal for (\ref{eq:method2qp}).
Note that there is at least one
index of $\{1,\ldots, k\}$
which is considered in infinitely many iterations.
Without loss of generality, we assume that
for all $\bp^i, i \in K$, 
$p^i_t$ is updated to generate  the next iteration
$\bp^{i+1}$.
As $\bp^*$ is not optimal for \eqref{eq:method2qp},
starting from $t, t+1, \ldots, k, 1, \ldots, t-1$,
there is a first component $\bar{t}$
for which
\begin{equation*}
  \sum_{j=1}^k Q_{\bar{t}j}p^*_j - (\bp^*)^T Q \bp^*\neq 0.
\end{equation*}
By applying one iteration of Algorithm \ref{alg:2} on 
$p^*_{\bar{t}}$, 
from an explanation similar to the proof of
Lemma \ref{new:lemma1},
we obtain $\bp^{*,n}$ satisfying
$p^{*,n}_{\bar{t}} \neq p^*_{\bar{t}}$.
Then by Lemma \ref{new:lemma1}, 
$$\frac{1}{2} (\bp^{*,n})^T Q \bp^{*,n} < \frac{1}{2} (\bp^{*})^T Q
\bp^{*}.$$
Assume it takes $\bar{i}$ steps from $t$ to $\bar{t}$
and $\bar{i} >1$,
\begin{eqnarray*}
\lim_{i \in K, i \rightarrow \infty}
p^{i+1}_t
&=& \lim_{i \in K, i \rightarrow \infty}
\frac{\frac{1}{Q_{tt}}(-\sum_{j:j \neq t} Q_{tj} p^i_t 
+(\bp^i)^T Q \bp^i)} 
{\frac{1}{Q_{tt}}(-\sum_{j:j \neq t} Q_{tj} p^i_t 
+(\bp^i)^T Q \bp^i) + \sum_{j: j \neq t} p^i_j} 
\\
&=& 
\frac{\frac{1}{Q_{tt}}(-\sum_{j:j \neq t} Q_{tj} p^*_t 
+(\bp^*)^T Q \bp^*)} 
{\frac{1}{Q_{tt}}(-\sum_{j:j \neq t} Q_{tj} p^*_t 
+(\bp^*)^T Q \bp^*) + \sum_{j: j \neq t} p^*_j} 
\\
&=& \frac{p^*_t}{\sum^k_{j=1} p^*_j}
=p^*_t,
\end{eqnarray*}
we have
\begin{equation*}
\lim_{i \in K, i \rightarrow \infty}
\bp^{i}
=
\lim_{i \in K, i \rightarrow \infty}
\bp^{i+1}
=
\cdots
=
\lim_{i \in K, i \rightarrow \infty}
\bp^{i+\bar{i}-1}
=\bp^*.
\end{equation*}
Moreover,
\begin{equation*}
\lim_{i \in K, i \rightarrow \infty}
\bp^{i+\bar{i}}
=\bp^{*,n}
\end{equation*}
and
\begin{eqnarray*}
\lim_{i \in K, i \rightarrow \infty} 
\frac{1}{2} (\bp^{i+\bar{i}})^T Q \bp^{i+\bar{i}}
&=& \frac{1}{2} (\bp^{*,n})^T Q \bp^{*,n}\\
&<& \frac{1}{2} (\bp^{*})^T Q \bp^{*}\\
&=& \lim_{i \in K, i \rightarrow \infty} \frac{1}{2} (\bp^{i})^T Q \bp^{i}.
\end{eqnarray*}
This contradicts the fact from Lemma \ref{new:lemma1}:
$$\frac{1}{2}(\bp^{1})^T Q \bp^{1} \geq \frac{1}{2}(\bp^{2})^T Q \bp^{2}
\geq \cdots
\geq \frac{1}{2} (\bp^{*})^T Q \bp^{*}.$$
Therefore, $\bp^*$ must be optimal for (\ref{eq:method2qp}).

\section{Implementation Notes of Algorithm \ref{alg:2}}
\label{sec:impl-notes-algor}

From Algorithm \ref{alg:2},
the main operation of each
iteration is on calculating 
$-\sum_{j:j\neq t} Q_{tj}p_j$
and 
$\bp^T Q\bp$,
both $O(k^2)$ procedures.
In the following, we show how to easily
reduce the cost per iteration to $O(k)$.

Following the notation in Lemma \ref{new:lemma1}
of Appendix \ref{sec:impl-notes-algor},
we consider $\bp$ the current solution.
Assume $p_t$ is the component to be updated, we 
generate $\bar{\bp}$ according to \eqref{eq:pbari}
and normalize $\bar{\bp}$ to the next iterate
$\bp^n$. Note that 
$\bar{\bp}$ is the same as $\bp$ except the 
$t$th component and we consider the 
form \eqref{eq:pbari1}.
Since 
$\sum_{i=1}^k p_i = 1$, 
\eqref{eq:normalizepn} is
$\bp^n = \bar{\bp}/(1+\Delta)$.
Throughout iterations, we keep the current 
$Q\bp$ and $\bp^T Q \bp$,
so $\Delta$ can be easily calculated. 
To obtain 
$Q\bp^n$ and $(\bp^n)^T Q \bp^n$, we use
\begin{eqnarray}
\nonumber
(Q\bp^n)_j
& = & \frac{(Q \bar{\bp})_j}{1+\Delta}\\
&=&
\frac{(Q \bp)_j + Q_{jt}\Delta}{1+\Delta},
j=1, \ldots, k,
  \label{eq:Qpj}
\end{eqnarray}
and
\begin{eqnarray}
  \label{eq:pqp}
(\bp^n)^T Q (\bp^n)&=&
\frac{\bar{\bp}^T Q \bar{\bp}
}
{(1+\Delta)^2}\\
\nonumber
&=&
\frac{\bp^T Q \bp
+2\Delta
\sum_{j=1}^k (Q\bp)_j
+ Q_{tt}\Delta^2
}
{(1+\Delta)^2}.
\end{eqnarray}
Both (and hence the whole iteration) takes
$O(k)$ operations.

Numerical inaccuracy may accumulate through
iterations, so gradually \eqref{eq:Qpj}
and \eqref{eq:pqp}
may be away from values directly calculated using
$\bp$. An easy prevention of this problem
is to recalculate 
$Q\bp$ and $\bp^T Q \bp$ directly using $\bp$
after several iterations (e.g., $k$ iterations).
Then, the average cost per iteration is still
$O(k) + O(k^2)/k = O(k)$.

\section{Derivation of \eqref{eq:pvopt}}
\label{sec:deriv-pvopt}
\begin{eqnarray*}
&& \sum^k_{i=1} 
\sum_{j:j \neq i} (I_{\{\rij > \rji\}} p_j- 
I_{\{\rji > \rij\}} p_i)^2\\
&=& \sum^k_{i=1} 
\sum_{j:j \neq i} (I_{\{\rij > \rji\}} p_j^2 
+ I_{\{\rji > \rij\}} p_i^2)\\
&=& 2 \sum^k_{i=1} 
(\sum_{j:j \neq i}
 I_{\{\rji > \rij\}}) p_i^2.
\end{eqnarray*}
If $\sum_{j:j \neq i}
 I_{\{\rji > \rij\}}\neq 0,
\forall i$,
then,
under the constraint $\sum_{i=1}^k p_i = 1$,
the optimal solution  satisfies
\begin{equation*}
  \frac{p_1}{
\sum_{j:j \neq 1}
 I_{\{r_{j1} > r_{1j}\}}}
=
\cdots 
=
  \frac{p_k}{
\sum_{j:j \neq k}
 I_{\{r_{jk} > r_{kj}\}}}.
\end{equation*}
Thus, \eqref{eq:pvopt}
is the optimal solution of \eqref{del-v}.
%-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
%300 ACC/MSE  SVM
%\afterpage{\clearpage
%\begin{figure}[H]
%\centering
%\resizebox{\linewidth}{!}{ \begin{tabular}{cc}
%       \subfigure[\dna by binary SVMs]{\input{fig-svm/300-svm-dna-acc.tex}}& 
%       \subfigure[\dna by binary SVMs]{\input{fig-svm/300-svm-dna-mse.tex}} \\ 
%
%       \subfigure[\waveform by binary SVMs]{\input{fig-svm/300-svm-waveform-acc.tex}}&
%       \subfigure[\waveform by binary SVMs]{\input{fig-svm/300-svm-waveform-mse.tex}} \\
%
%       \subfigure[\satimage by binary SVMs]{\input{fig-svm/300-svm-satimage-acc.tex}}&
%       \subfigure[\satimage by binary SVMs]{\input{fig-svm/300-svm-satimage-mse.tex}} \\
%
%       \subfigure[\segment by binary SVMs]{\input{fig-svm/300-svm-segment-acc.tex}}&
%       \subfigure[\segment by binary SVMs]{\input{fig-svm/300-svm-segment-mse.tex}}
%\end{tabular}}
%\end{figure}
%
%\begin{figure}[H]
%\resizebox{\linewidth}{!}{ \begin{tabular}{cc}
%       \subfigure[\USPS by binary SVMs]{\input{fig-svm/300-svm-usps-acc.tex}}&
%       \subfigure[\USPS by binary SVMs]{\input{fig-svm/300-svm-usps-mse.tex}} \\
%
%       \subfigure[\MNIST by binary SVMs]{\input{fig-svm/300-svm-mnist-acc.tex}}&
%       \subfigure[\MNIST by binary SVMs]{\input{fig-svm/300-svm-mnist-mse.tex}} \\
%
%       \subfigure[\letter by binary SVMs]{\input{fig-svm/300-svm-letter-acc.tex}}&
%       \subfigure[\letter by binary SVMs]{\input{fig-svm/300-svm-letter-mse.tex}}
%\end{tabular}}
%\caption{300-500-accuracy/mse by SVM}
%\end{figure}
%}
%
%%-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
%%300 ACC/MSE RF
%\afterpage{\clearpage
%\begin{figure}[H]
%\centering
%\resizebox{\linewidth}{!}{ \begin{tabular}{cc}
%       \subfigure[\dna by binary random forests]{\input{fig-rf/300-rf-dna-acc.tex}}& 
%       \subfigure[\dna by binary random forests]{\input{fig-rf/300-rf-dna-mse.tex}} \\ 
%
%       \subfigure[\waveform by binary random forests]{\input{fig-rf/300-rf-waveform-acc.tex}}&
%       \subfigure[\waveform by binary random forests]{\input{fig-rf/300-rf-waveform-mse.tex}} \\
%
%       \subfigure[\satimage by binary random forests]{\input{fig-rf/300-rf-satimage-acc.tex}}&
%       \subfigure[\satimage by binary random forests]{\input{fig-rf/300-rf-satimage-mse.tex}} \\
%
%       \subfigure[\segment by binary random forests]{\input{fig-rf/300-rf-segment-acc.tex}}&
%       \subfigure[\segment by binary random forests]{\input{fig-rf/300-rf-segment-mse.tex}}
%\end{tabular}}
%\end{figure}
%
%\begin{figure}[H]
%\resizebox{\linewidth}{!}{ \begin{tabular}{cc}
%       \subfigure[\USPS by binary random forests]{\input{fig-rf/300-rf-usps-acc.tex}}&
%       \subfigure[\USPS by binary random forests]{\input{fig-rf/300-rf-usps-mse.tex}} \\
%
%       \subfigure[\MNIST by binary random forests]{\input{fig-rf/300-rf-mnist-acc.tex}}&
%       \subfigure[\MNIST by binary random forests]{\input{fig-rf/300-rf-mnist-mse.tex}} \\
%
%       \subfigure[\letter by binary random forests]{\input{fig-rf/300-rf-letter-acc.tex}}&
%       \subfigure[\letter by binary random forests]{\input{fig-rf/300-rf-letter-mse.tex}}
%\end{tabular}}
%\caption{300-500-accuracy/mse by random forest}
%\end{figure}
%}
%XXX
%%--------------------
% 300-500 with 'SVM' 'RF' MSE
\afterpage{\clearpage
\begin{figure}[H]
%\centering
\resizebox{\linewidth}{!}{ \begin{tabular}{cc}
        \subfigure[\dnak by binary SVMs]{\input{fig-svm/combo-svm-dna-mse.tex}}& 
        \subfigure[\dnak by binary random forests]{\input{fig-rf/combo-rf-dna-mse.tex}} \\ 
        \subfigure[\waveformk by binary SVMs]{\input{fig-svm/combo-svm-waveform-mse.tex}}&
        \subfigure[\waveformk by binary random forests]{\input{fig-rf/combo-rf-waveform-mse.tex}} \\
        \subfigure[\satimagek by binary SVMs]{\input{fig-svm/combo-svm-satimage-mse.tex}}&
        \subfigure[\satimagek by binary random forests]{\input{fig-rf/combo-rf-satimage-mse.tex}} \\
        \subfigure[\segmentk by binary SVMs]{\input{fig-svm/combo-svm-segment-mse.tex}}&
        \subfigure[\segmentk by binary random forests]{\input{fig-rf/combo-rf-segment-mse.tex}}
\end{tabular}}
\end{figure}

\begin{figure}[H]
\resizebox{\linewidth}{!}{ \begin{tabular}{cc}
        \subfigure[\uspsk by binary SVMs]{\input{fig-svm/combo-svm-usps-mse.tex}}&
        \subfigure[\uspsk by binary random forests]{\input{fig-rf/combo-rf-usps-mse.tex}} \\
        \subfigure[\mnistk by binary SVMs]{\input{fig-svm/combo-svm-mnist-mse.tex}}&
        \subfigure[\mnistk by binary random forests]{\input{fig-rf/combo-rf-mnist-mse.tex}} \\
        \subfigure[\letterk by binary SVMs]{\input{fig-svm/combo-svm-letter-mse.tex}}&
        \subfigure[\letterk by binary random forests]{\input{fig-rf/combo-rf-letter-mse.tex}}
\end{tabular}}
\caption{MSE by using four probability estimates methods based on
binary SVMs (left) and binary random forests (right).
MSE of $\delta_V$ is too large and is not presented.
solid line: 300 training/500 testing points;
dotted line: 800 training/1,000 testing points.}
\label{fig:svm-rf-mse}
\end{figure}
}
%--------------------
%-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
% 300-500 with 'SVM' 'RF' accuracy
\afterpage{\clearpage
\begin{figure}[H]
\centering
\resizebox{\linewidth}{!}{ \begin{tabular}{cc}
                \subfigure[\dnak by binary SVMs; $\delta_{DV}=10.82\%$]{\input{fig-svm/300-svm-dna-acc.tex}}&
                \subfigure[\dnak by binary random forests; $\delta_{RF}=8.74\%$]{\input{fig-rf/300-rf-dna-acc.tex}} \\
                \subfigure[\waveformk by binary SVMs; $\delta_{DV}=16.47\%$ ]{\input{fig-svm/300-svm-waveform-acc.tex}}&
        \subfigure[\waveformk by binary random forests; $\delta_{RF}=17.39\%$]{\input{fig-rf/300-rf-waveform-acc.tex}} \\
        \subfigure[\satimagek by binary SVMs; $\delta_{DV}=14.88\%$]{\input{fig-svm/300-svm-satimage-acc.tex}}&
        \subfigure[\satimagek by binary random forests;$\delta_{RF}=14.74\%$]{\input{fig-rf/300-rf-satimage-acc.tex}} \\
        \subfigure[\segmentk by binary SVMs; $\delta_{DV}=5.82\%$]{\input{fig-svm/300-svm-segment-acc.tex}}&
        \subfigure[\segmentk by binary random forests; $\delta_{RF}=5.23\%$]{\input{fig-rf/300-rf-segment-acc.tex}}
\end{tabular}}
\end{figure}

\begin{figure}[H]
\resizebox{\linewidth}{!}{ \begin{tabular}{cc}
        \subfigure[\uspsk by binary SVMs; $\delta_{DV}=10.71\%$]{\input{fig-svm/300-svm-usps-acc.tex}}&
        \subfigure[\uspsk by binary random forests; $\delta_{RF}=14.28\%$]{\input{fig-rf/300-rf-usps-acc.tex}} \\
        \subfigure[\mnistk by binary SVMs; $\delta_{DV}=13.02\%$ ]{\input{fig-svm/300-svm-mnist-acc.tex}}&
        \subfigure[\mnistk by binary random forests; $\delta_{RF}=16.18\%$]{\input{fig-rf/300-rf-mnist-acc.tex}} \\
        \subfigure[\letterk by binary SVMs; $\delta_{DV}=32.29\%$]{\input{fig-svm/300-svm-letter-acc.tex}}&
        \subfigure[\letterk by binary random forests; $\delta_{RF}=32.55\%$]{\input{fig-rf/300-rf-letter-acc.tex}}
\end{tabular}}
\caption{Average of 20 test errors by five probability estimates
methods based on binary SVMs (left) and binary random forests (right).
Each of the 20 test errors is by 300 training/500 testing points.
Caption of each sub-figure shows the averaged error by voting using
pairwise SVM decision values ($\delta_{DV}$) and the multi-class
random forest ($\delta_{RF}$).}
\label{fig:300-500-accuracy}
\end{figure}
}
%%--------------------
%% 800-1000 with 'SVM' 'RF' accuracy
\afterpage{\clearpage
\begin{figure}[H]
\centering
\resizebox{\linewidth}{!}{ \begin{tabular}{cc}
        \subfigure[\dnak by binary SVMs; $\delta_{DV}=6.31\%$]{\input{fig-svm/800-svm-dna-acc.tex}}& 
        \subfigure[\dnak by binary random forests; $\delta_{RF}=6.91\%$]{\input{fig-rf/800-rf-dna-acc.tex}} \\ 

        \subfigure[\waveformk by binary SVMs; $\delta_{DV}=14.33\%$ ]{\input{fig-svm/800-svm-waveform-acc.tex}}&
        \subfigure[\waveformk by binary random forests; $\delta_{RF}=15.66\%$]{\input{fig-rf/800-rf-waveform-acc.tex}} \\

        \subfigure[\satimagek by binary SVMs; $\delta_{DV}=11.54\%$]{\input{fig-svm/800-svm-satimage-acc.tex}}&
        \subfigure[\satimagek by binary random forests; $\delta_{RF}=11.92\%$]{\input{fig-rf/800-rf-satimage-acc.tex}} \\

        \subfigure[\segmentk by binary SVMs; $\delta_{DV}=3.30\%$]{\input{fig-svm/800-svm-segment-acc.tex}}&
        \subfigure[\segmentk by binary random forests; $\delta_{RF}=2.77\%$]{\input{fig-rf/800-rf-segment-acc.tex}}
\end{tabular}}
\end{figure}

\begin{figure}[H]
\resizebox{\linewidth}{!}{ \begin{tabular}{cc}
        \subfigure[\uspsk by binary SVMs; $\delta_{DV}=7.54\%$]{\input{fig-svm/800-svm-usps-acc.tex}}&
        \subfigure[\uspsk by binary random forests; $\delta_{RF}=10.23\%$]{\input{fig-rf/800-rf-usps-acc.tex}} \\

        \subfigure[\mnistk by binary SVMs; $\delta_{DV}=7.77\%$ ]{\input{fig-svm/800-svm-mnist-acc.tex}}&
        \subfigure[\mnistk by binary random forests;  $\delta_{RF}=10.10\%$]{\input{fig-rf/800-rf-mnist-acc.tex}} \\

        \subfigure[\letterk by binary SVMs; $\delta_{DV}=18.61\%$]{\input{fig-svm/800-svm-letter-acc.tex}}&
        \subfigure[\letterk by binary random forests; $\delta_{RF}=20.25\%$]{\input{fig-rf/800-rf-letter-acc.tex}}
\end{tabular}}
\caption{Average of 20 test errors by five probability estimates
methods based on binary SVMs (left) and binary random forests (right).
Each of the 20 test errors is by 800 training/1,000 testing points.
Caption of each sub-figure shows the averaged error by voting using
pairwise SVM decision values ($\delta_{DV}$) and the multi-class
random forest ($\delta_{RF}$).}
\label{fig:800-1000-accuracy}
\end{figure}
}

% 800-1000 with 'SVM' 'RF' MSE
%\afterpage{\clearpage
%\begin{figure}[H]
%\centering
%\resizebox{\linewidth}{!}{ \begin{tabular}{cc}
%       \subfigure[\dnak by binary SVMs]{\input{fig-svm/800-svm-dna-mse.tex}}& 
%       \subfigure[\dnak by binary random forests]{\input{fig-rf/800-rf-dna-mse.tex}} \\ 
%
%       \subfigure[\waveformk by binary SVMs; $\delta_{DV}=??\%$ ]{\input{fig-svm/800-svm-waveform-mse.tex}}&
%       \subfigure[\waveformk by binary random forests]{\input{fig-rf/800-rf-waveform-mse.tex}} \\
%
%       \subfigure[\satimagek by binary SVMs; $\delta_{DV}=??\%$]{\input{fig-svm/800-svm-satimage-mse.tex}}&
%       \subfigure[\satimagek by binary random forests]{\input{fig-rf/800-rf-satimage-mse.tex}} \\
%
%       \subfigure[\segmentk by binary SVMs; $\delta_{DV}=??\%$]{\input{fig-svm/800-svm-segment-mse.tex}}&
%       \subfigure[\segmentk by binary random forests]{\input{fig-rf/800-rf-segment-mse.tex}}
%\end{tabular}}
%\end{figure}
%
%\begin{figure}[H]
%\resizebox{\linewidth}{!}{ \begin{tabular}{cc}
%       \subfigure[\uspsk by binary SVMs; $\delta_{DV}=??\%$]{\input{fig-svm/800-svm-usps-mse.tex}}&
%       \subfigure[\uspsk by binary random forests]{\input{fig-rf/800-rf-usps-mse.tex}} \\
%
%       \subfigure[\mnistk by binary SVMs; $\delta_{DV}=??\%$ ]{\input{fig-svm/800-svm-mnist-mse.tex}}&
%       \subfigure[\mnistk by binary random forests]{\input{fig-rf/800-rf-mnist-mse.tex}} \\
%
%       \subfigure[\letterk by binary SVMs; $\delta_{DV}=??\%$]{\input{fig-svm/800-svm-letter-mse.tex}}&
%       \subfigure[\letterk by binary random forests;  10.66 of generic random forest]{\input{fig-rf/800-rf-letter-mse.tex}}
%\end{tabular}}
%\caption{800-1000-mse}
%\end{figure}
%}
%-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

\afterpage{\clearpage
\begin{figure}[H]
%\centering
\resizebox{\linewidth}{!}{ \begin{tabular}{cc}
        \subfigure[\dnak by binary SVMs]{\input{fig-svm/combo-svm-dna-lgl.tex}}& 
        \subfigure[\dnak by binary random forests]{\input{fig-rf/combo-rf-dna-lgl.tex}} \\ 
        \subfigure[\waveformk by binary SVMs]{\input{fig-svm/combo-svm-waveform-lgl.tex}}&
        \subfigure[\waveformk by binary random forests]{\input{fig-rf/combo-rf-waveform-lgl.tex}} \\
        \subfigure[\satimagek by binary SVMs]{\input{fig-svm/combo-svm-satimage-lgl.tex}}&
        \subfigure[\satimagek by binary random forests]{\input{fig-rf/combo-rf-satimage-lgl.tex}} \\
        \subfigure[\segmentk by binary SVMs]{\input{fig-svm/combo-svm-segment-lgl.tex}}&
        \subfigure[\segmentk by binary random forests]{\input{fig-rf/combo-rf-segment-lgl.tex}}
\end{tabular}}
\end{figure}

\begin{figure}[H]
\resizebox{\linewidth}{!}{ \begin{tabular}{cc}
        \subfigure[\uspsk by binary SVMs]{\input{fig-svm/combo-svm-usps-lgl.tex}}&
        \subfigure[\uspsk by binary random forests]{\input{fig-rf/combo-rf-usps-lgl.tex}} \\
        \subfigure[\mnistk by binary SVMs]{\input{fig-svm/combo-svm-mnist-lgl.tex}}&
        \subfigure[\mnistk by binary random forests]{\input{fig-rf/combo-rf-mnist-lgl.tex}} \\
        \subfigure[\letterk by binary SVMs]{\input{fig-svm/combo-svm-letter-lgl.tex}}&
        \subfigure[\letterk by binary random forests]{\input{fig-rf/combo-rf-letter-lgl.tex}}
\end{tabular}}
\caption{Log likelihood \eqref{eq:log} by using four probability estimates methods based on
binary SVMs (left) and binary random forests (right).
MSE of $\delta_V$ is too small and is not presented.
solid line: 300 training/500 testing points;
dotted line: 800 training/1,000 testing points.}
\label{fig:svm-rf-lgl}
\end{figure}
}
%--------------------

\bibliography{sdp}

\end{document}

\pagebreak
\pagestyle{empty}

\subsection*{Associate editor}
\begin{center}
\textsf{Probability Estimates for Multi-class Classification by Pairwise Coupling} \\
Ting-Fan Wu, Chih-Jen Lin
and  Ruby C. Weng
\end{center}

Thank you for your carefully reading of our paper. We have tried to
address your comments in the revised manuscript. 

\begin{enumerate}
\item ``the wording is rough and it does not read
fluently''

We have largely improved the presentation of this paper.
A native speaker helps to check our writing.

\item ``the first reviewer expressed concerns regarding
the SVM binary probability estimates''

We agree with the first reviewer
that for SVM, the decision values
are clustered at $\pm 1$. 
We practically observe this on decision values
of training data. However, the situation of 
decision values by cross validation is not serious.
In all our experiments, only CV decision values 
are used, so in practice they might be fine
for providing probability information.
We now briefly address this 
issue in Section \ref{svm-binary}
and cite related references.

Furthermore, if the SVM probability estimate is erroneous,
our multi-class probability
estimates in Section \ref{svm-binary}
will be inaccurate.
However, the performance (accuracy and  MSE)
is similar to (or even better than) that by using
random forest in Section \ref{rf-binary}. 
Thus, the binary SVM
probability output is somewhat comparable to that
by binary random forest.

Note that this paper is for general multi-class 
probability estimates, where SVM is used only as
an example. Most experiments (e.g. Sections 
\ref{sec:experiments}
and \ref{rf-binary}) are not related to SVM.


\item ``please compare your approach to ASS''

We compare with ASS in Section \ref{sec:misc-observ}
by using both decision values and 
$\rij-1/2$ as the hypothesis in the exponential 
loss-based decoding. An interesting link between
this approach and $\delta_{HT}$ is also discussed.
We did not put the comparison
in Sections \ref{svm-binary} and \ref{rf-binary},
where the main purpose is to compare
different probability estimates. 

\item ``need to provide the log-likelihood by 
different approaches''

We have presented this information in Figure \ref{fig:svm-rf-lgl}.
Results are consistent with those using MSE and accuracy.

\end{enumerate}

\clearpage
\pagebreak



\subsection*{Reviewer 1}
\begin{center}
\textsf{Probability Estimates for Multi-class Classification by Pairwise Coupling} \\
Ting-Fan Wu, Chih-Jen Lin
and  Ruby C. Weng
\end{center}

Thank you for your carefully reading of our paper. 
We have tried to
address your comments in the revised manuscript. 

We agree that for SVM, the decision values
are clustered at $\pm 1$. Thus, they may
not carry reliable probability information.
We practically observe this on decision values
of training data. However, the situation of 
decision values by cross validation is not serious.
In all our experiments, only CV decision values 
are used, so in practice they might be fine
for providing probability information.
We now briefly address this issue in Section \ref{svm-binary}
and cite related references.


Furthermore, if the SVM probability estimate is erroneous,
our multi-class probability
estimates in Section \ref{svm-binary}
will be inaccurate.
However, the performance (accuracy and  MSE)
is similar to (or even better than) that by using
random forest in Section \ref{rf-binary}. 
Thus, the binary SVM
probability output is somewhat comparable to that
by binary random forest.

Note that this paper is for general multi-class 
probability estimates, where SVM is used only as
an example. Most experiments (e.g. Sections 
\ref{sec:experiments}
and \ref{rf-binary}) are not related to SVM.

\clearpage
\pagebreak

\subsection*{Reviewer 2}
\begin{center}
\textsf{Probability Estimates for Multi-class Classification by Pairwise Coupling} \\
Ting-Fan Wu, Chih-Jen Lin
and  Ruby C. Weng
\end{center}

Thank you for your carefully reading of our paper. 
We are glad that you find this paper
instructive and insightful.
We have tried to
address your comments in the revised manuscript. 

Specific comments:

\begin{enumerate}
\item ``$\sum_{j=1}^k Q_{ij} = 1$.''


Yes, you are right. It should be 
$\sum_{i=1}^k Q_{ij} = 1$.

\item ``in Eq. \eqref{eq:method2qp3}, 
$1/2$ should be 1''

Yes. We now remove $1/2$ from the first formula of
method 2 for consistence with method 1, so
\eqref{eq:method2qp3} should be 
$2\bp^T Q \bp$. We then explain that it is equivalent
to $\bp^T Q \bp/2$.


\item ``KKT when Lagrange is sufficient''

They are the same thing here. Constraint qualification
automatically satisfies for linear constraints
(both equalities and inequalities).
Such things depend on how we think about the optimization
theorem behind. If we consider problems with only
linear constraints, when stating KKT, there is
no need to mention constraint qualification. If
we consider problems with general constraints,
when deriving the Lagrangian dual, constraint 
qualification is needed.

Anyway, now I remove KKT and simply say ``
if and only if it satisfies the optimality condition.''
Then readers can interpret it in their own way.


\item ``Section \ref{sec:solvemethod2} is
not central to the discussion''

Section \ref{sec:linear-system-method2}
implies convexity but not strict convexity.
Hence it is possible to have multiple solutions.
Moreover, due to that $Q$ is only positive semi-definite,
numerical linear algebra techniques may not directly
apply. So it is quite important to have 
some discussion.
In Theorem \ref{th:Qpd}, we prove that even though
$Q$ may be only positive semi-definite, 
\eqref{eq:method2qp} has a unique minimum. We think
this is an important property.


\item ``How is \eqref{del-ht} verified and how you get
rid of $\nij$ ?''

Yes, there is a minor problem here and thank you
for the pointer. In fact, when \cite{TH98a}
proved \eqref{eq:pt}
in their Theorem 1, equal weighting (i.e., $\nij=1$)
is assumed. We forgot to mention this earlier.
All our discussion is under 
this assumption, so for connecting
\eqref{eq:pt} to \eqref{del-ht},
there is no need to worry about $\nij$.
We now elaborate more on this in the end of
Section \ref{sec:analys-meth-hast}.

By the way, our proposed methods can be extended
to a weighted version as well.

\end{enumerate}

\clearpage
\pagebreak

\subsection*{Reviewer 3}
\begin{center}
\textsf{Probability Estimates for Multi-class Classification by Pairwise Coupling} \\
Ting-Fan Wu, Chih-Jen Lin
and  Ruby C. Weng
\end{center}

Thank you for your carefully reading of our paper. 
We have tried to
address your comments in the revised manuscript. 


General comments:

\begin{enumerate}
\item ``derivation of the theoretical part is not detailed
enough''

We have provided all technical details. Please also see
the reply to ``comments to authors.''
\item ``Comparison to other multi-class classifiers.''

This paper is mainly for probability estimates, but not
for classification. In Section \ref{sec:experiments},
we can check the difference to the real
probability, but in Section \ref{real}
for real data, the true probability
is not available. Hence, we present
the accuracy, Brier score (difference
between probabilities and true labels),
and the log likelihood. For 
multi-class method which obtains only class labels,
the use of Brier score (now difference between predicted
and true labels) and log likelihood may not be very suitable.
Thus, we did not compare with other methods in Sections
\ref{svm-binary}
and \ref{rf-binary}.

Regarding SVM accuracy, in the caption of each subfigure,
we show the accuracy of pairwise voting 
via decision values (called $\delta_{DV}$),
which is one of the standard ways for SVM 
multi-class classification. Its accuracy is similar to 
using $\delta_V $ and $\delta_2$.
There are already excellent papers comparing 
$\delta_{DV}$ with methods like 1vsall or that by
Crammer and Singer. Roughly how these methods compared to 
$\delta_{DV}$ is the same as compared to
$\delta_V$ and $\delta_2$. We have now mentioned this point.

By the way, in Section \ref{sec:misc-observ},
Miscellaneous Observations and Discussion,
we now provide a simple comparison to another multi-class method 
based on pairwise coupling: loss-based decoding in 
\cite{ELA01a}. As it does not directly provide probabilities,
we compare the accuracy only. 
Table \ref{tbl:ass-accuracy}
shows that when the number of classes is larger,
this approach is slightly worse than our prediction
via probability estimates

For random forest, we find that pairwise 
prediction is sometimes worse than the built-in predictor
in RF. RF 
involves with decision trees which by default can
handle multi-class data. To study this observation,
we need to investigate RF itself and hence we list
it as a future issue.

\end{enumerate}


Comments to authors:

\begin{enumerate}
\item ``relation between $\rij$ and $\mij$''

We have included more details.


\item ``$\delta$ as the decision rule''

This is an excellent suggestion. We clarify
this point in the beginning of Section 
\ref{sec:surv-exist-meth}.


\item ``$p^v_i$''

Similar to $\delta_*$, we have provided suitable definitions.


\item ``Section \ref{sec:voting}: $I$ not defined''

Defined.
\item ``Section \ref{sec:using-part-estimated}: $\rij$'s''

We have removed all 's.


\item ``Section \ref{sec:analys-meth-hast}: model
different from \cite{TH98a}''

They are the same. Note that
\begin{eqnarray*}
  l(\bp)& =& 
\sum_{i \neq j} \nij
\rij
\log\frac{\rij}{\mij},
\\
& = &
\nonumber
\sum_{i < j} \nij
\left(\rij
\log\frac{\rij}{\mij}
+
(1-\rij)
\log\frac{1-\rij}{1-\mij}\right),
\end{eqnarray*}
where the latter is the form used in \cite{TH98a}.
We now have both forms in the paper.


A problem in the earlier version is that we did not
indicate that Eq. \eqref{eq:pt}
(i.e., Theorem 1 of \cite{TH98a})
and all discussion after then are under the assumption
of equal weighting (i.e., $\nij=1$). We have clarified
this in Section \ref{sec:analys-meth-hast}.



\item ``Section \ref{sec:analys-meth-hast}: how you 
derive \eqref{eq:lp'=0} from \eqref{eq:lp'} ?''

When $\partial l(\bp)/\partial p_i =0$, multiplying
$p_i$ on all terms of \eqref{eq:lp'}
gives \eqref{eq:lp'=0}. We now explain all
the details.

\item ``Section \ref{sec:method-linear-system}:
$k/2$ instead of $2/k$''

Yes, this should be $2/k$. Thanks for the pointer.
\item ``Section \ref{sec:solving-refm0}:
$\sum_{j=1}^k Q_{ij}=1$''

Yes, it should be 
$\sum_{i=1}^k Q_{ij}=1$.


\item ``Section \ref{sec:solving-refm0}: need more explanation''

We have explained more on the uniqueness of the Markov
system and discuss practically how to solve it.

\item ``Derivation of \eqref{re-m1} from 
$Q\bp=\bp$''

Elaborated.
\item ``Section \ref{sec:linear-system-method2}:
Maybe Theorem \ref{th:pnonneg}
as lemma''

We are fine to do so, but another reviewer thinks
that Theorem \ref{th:pnonneg}
is important (even more than Theorems
\ref{th:Qpd} and \ref{new}).
So we keep it as a theorem.
\item ``Explain \eqref{eq:method2linearsys}''

Done.
\item ``After Theorem \ref{new},
need to say the proof is in 
Appendix \ref{sec:proof-theorem-refnew}.''

Done.
\item ``increase marker size of Figure \ref{fig:sim}''

Done.
\end{enumerate}


\end{document}


\section{Other Details}

Regarding the property of the solution to the first method,
we prove the following theorem:
\begin{theorem}
Assume $\rij > 0, \forall i \neq j$.
  If there is an index $i$
such that for all 
$j \neq i$,
\begin{equation}
  \sum_{l: l \neq i} r_{il}
+ \rij 
\geq 2 \sum_{l: l \neq j} r_{jl},\label{eq:ril+rij}
\end{equation}
then the vector $\bp$ obtained by solving \eqref{m1}
satisfies
\begin{equation*}
  p_i > p_j, \forall j \neq i.
\end{equation*}
\label{th:pilargest}
\end{theorem}


A way to interpret the condition
\eqref{eq:ril+rij}
is as follows: if the data is in the $i$th
class, in average
$2\sum_{l:l\neq j} \approx k-1$.
If $r_{il}$ are all close to one,
$  \sum_{l: l \neq i} r_{il}
+ \rij 
\geq k-1$. Therefore, if \eqref{eq:ril+rij}
holds, we are quite sure that the data should be in the
$i$th class and our approach leads to such a conclusion. 
In contrast, Friedman's rule does not have this property.

Facts:
\begin{enumerate}
\item \eqref{eq:lp} is a convex function of $\mu$.
We have
\begin{gather*}
  f(\mu)
= -r \log \frac{r}{\mu} - (1-r) \log \frac{1-r}{1-\mu}\\
f'(\mu) = -r/\mu + (1-r)/(1-\mu)\\
f^{''}(\mu)
= r/\mu^2 + (1-r)/(1-\mu)^2 > 0,
\end{gather*}
if $0 < \mu  < 1$. Moreover, $f(\mu)$ has minimum at
$\mu = r$.
\item The gradient of $l(p)$ is 
  \begin{eqnarray}
    \frac{\partial l(p)}{\partial p_i}
& = & 
\frac{\partial}{\partial p_i}
\biggl(\sum_{j=i+1}^n \nij
(-\rij \log p_i  - (1-\rij)\log p_j
+ \log(p_i + p_j))\\
&& \qquad
\sum_{k=1}^{i-1}
n_{ki}(-r_{ki} \log p_k
- (1-r_{ki})\log p_i
+ \log (p_k + p_i)\biggr)
\nonumber \\
& = & 
\sum_{j=i+1}^n \nij
\bigl(-\frac{\rij}{p_i}
+ \frac{1}{p_i + p_j}\bigr)
+
\sum_{k=1}^{i-1}
n_{ki}\bigl(-\frac{r_{ki}}{p_k}
+ \frac{1}{p_k + p_i}\bigr) \nonumber \\
& = & 
\sum_{j: i \neq j}
n_{ij}
\bigl(-\frac{\rij}{p_i} + 
\frac{1}{p_i + p_j}\bigr) \nonumber
  \end{eqnarray}

%somi
\begin{eqnarray}
        \label{eq:lpdec}
        l(p')-l(p) &=& \sum_{i \neq j} \nij (\rij \log\frac{\rij}{\mij'}
         + ( 1- \rij) \log \frac{1-\rij}{1-\mij'})\\  &&
        -\sum_{i \neq j} \nij (\rij \log\frac{\rij}{\mij}
         + ( 1- \rij) \log \frac{1-\rij}{1-\mij}) \\
        &=& \sum_{i \neq j} \nij (\rij \log \frac{\mij}{\mij'} 
                        + (1 - \rij) \log\frac{1 - \mij}{1 - \mij'}) \\
        &=& \sum_{i \neq j} \nij 
        (\rij \log\frac{\mij}{\frac{\alpha\mij}{\alpha\mij + \mji}} 
                        + (1 - \rij) \log\frac{\mji}{\frac{\mji}{\alpha\mij+\mji}}) \\
        &=& \sum_{i \neq j} \nij 
        (\rij \log\frac{\alpha\mij + \mji}{\alpha} 
                        + (1 - \rij) \log(\alpha\mij+\mji)) \\
        &=& \sum_{i \neq j} \nij 
        (\log(\alpha\mij+\mji) -\rij \log \alpha) \\
        &=& \sum_{i \neq j} \nij 
        (\log(\alpha\mij + 1- \mij)-\rij \log \alpha ) \\
        &=& \sum_{i \neq j} \nij 
        ( \log((\alpha-1)\mij + 1) - \rij \log \alpha ) \\
        &\leq&\sum_{i \neq j} \nij 
        ( (\alpha-1)\mij - \rij \log \alpha ) \\
        &=& (\alpha-1)\sum_{i \neq j}\nij\mij - \log\alpha \sum_{i \neq j}\nij\rij \\
        &=&(\alpha-1)\frac{1}{\alpha}\sum_{i \neq j}\nij\rij - \log\alpha \sum_{i \neq j}\nij\rij  \\
        &=& (1 -\frac{1}{\alpha} - \log\alpha )\sum_{i \neq j}\nij\rij  \\
        &=& -C\sum_{i \neq j}\nij\rij  (const. C \geq 0)\\
        &\leq& 0
\end{eqnarray}

\end{enumerate}


\end{document}

% Local Variables:
% TeX-master: "jmlr_svmprob.tex"
% End:


Dear Yoram,

Enclosed please find the revision of MS # 03-541 which has earlier
been accepted subject to a revision.  

We have completed a thorough revision in which all suggestions by you
and referees have been addressed.  Detailed description of how we have
modified the manuscript can be found in the attached list.

We hope that you find all revisions satisfactory.
Thank you again for handling the paper.

Best,

Chih-Jen
