Predictive Approaches For Gaussian Process Classifier Model Selection \\ September 10, 2008
\addr Yahoo! Labs
Bangalore
India
\addr Yahoo! Research
Santa Clara
USA
% <- trailing '%' for backward compatibility of .sty file In this paper we consider the problem of Gaussian process classifier (GPC) model selection with different Leave-One-Out (LOO) Cross Validation (CV) based optimization criteria and provide a practical algorithm using LOO predictive distributions with such criteria to select hyperparameters. Apart from the standard average negative logarithm of predictive probability (NLP), we also consider smoothed versions of criteria such as F-measure and Weighted Error Rate (WER), which are useful for handling imbalanced data. Unlike the regression case, LOO predictive distributions for the classifier case are intractable. We %Keerthi: use approximate LOO predictive distributions arrived from Expectation Propagation (EP) approximation. We conduct experiments on several real world benchmark datasets. When the NLP criterion is used for optimizing the hyperparameters, the predictive approaches show better or comparable NLP generalization performance with existing GPC approaches. On the other hand, when the F-measure criterion is used, the F-measure generalization performance improves significantly on several datasets. Overall, the EP-based predictive algorithm comes out as an excellent choice for GP classifier model selection with different optimization criteria.
cvgpcappendixv5.tex
\section*{Appendix: Optimization criteria, approximate predictive distributions and their derivatives}
From the definitions of various optimization criteria like NLP measure (\ref{nlp}) and (\ref{FMmooth}) etc, we see that the chosen measure and its derivatives can be obtained from the LOO predictive distributions  $p(y_i|x_i,{\it S}_{\setminus{i}},\bftheta)$ and their derivatives ${{\partial p(y_i|x_i,{\it S}_{\setminus{i}},\bftheta)} \over {\partial \theta_j}}$. Here, $\theta_j$ is $j^{th}$ component of the hyperparameter vector $\bftheta$. Note that the derivative of the NLP measure (\ref{nlp}) is given by:
\begin{equation}
	{{\partial G(\bftheta)} \over {\partial \theta_j}}\:=\:- {1\over n} \sum_{i=1}^n {1\over \Phi(y_i z_i)} {{\partial p(y_i|x_i,{\it S}_{\setminus{i}},\bftheta)} \over {\partial \theta_j}}
\label{nlpderv}
\end{equation}
and the derivative of the smoothed F-measure (\ref{FMmooth}) is given by:
\begin{equation}
	 	{{\partial F_{\zeta}(\bftheta)} \over {\partial \theta_j}}\:=\:{{\eta(\bftheta) {{\partial A(\bftheta)} \over {\partial \theta_j}} - A(\bftheta) (1-\zeta) {{\partial m_+(\bftheta)} \over {\partial \theta_j}}} \over {\eta^2(\bftheta)}}
\label{FMoothderv}
\end{equation}
where $\eta(\bftheta)\:=\:\zeta n_{+} + (1-\zeta)m_{+}(\bftheta)$. Note that the derivatives ${{\partial A(\bftheta)} \over {\partial \theta_j}}$ and ${{\partial m_+(\bftheta)} \over {\partial \theta_j}}$ are directly dependent on ${{\partial p(y_i=+1|x_i,{\it S}_{\setminus{i}},\bftheta)} \over {\partial \theta_j}}$. Now, from (\ref{Pyi}) we see that to define the LOO predictive distributions, we need the LOO mean $\mu_{\setminus{i}}$ and variance $\sigma^2_{\setminus{i}}$. In the case of EP approximation, analytical expressions to compute these quantities are already available (see eqn. (\ref{cavmu})). Next, we give details on how the derivatives of predictive distributions can be obtained with these approximations.  

%\subsection*{LOO Mean and Variance: Laplace Approximation} 
%The predictive distributions $p(f_i|x_i,{\it S},\bftheta)$ are approximated as $\mathcal{N}(\mu_{\setminus{i}},\sigma^2_{\setminus{i}})$. The expressions for the LOO mean $\mu_{\setminus{i}}$ and variance $\sigma^2_{\setminus{i}}$ are derived by relating the posterior of latent function $q({\bf f}|{\it S},\bftheta)$ from Laplace and EP approximations; this approach is simple as we have cavity distributions in EP directly related to the posterior. Specifically, the posterior from Laplace approximation is given by: $q({\bf f}|{\it S},\bftheta)\:=\:\mathcal{N}({\hat {\bf f}},({\bf K}^{-1}+{\bf W})^{-1})$ where ${\hat {\bf f}}\:=\:{\bf K} \bfalp$ and $\bfalp\:=\:\nabla \log p({\bf y}|{\hat {\bf f}})$. Similarly, the posterior from EP approximation is given by:$q({\bf f}|{\it S},\bftheta)\:=\:\mathcal{N}({\bf m},{\bf C})$ where ${\bf m}\:=\:{\bf C}\bfSigma^{-1}{\bfmu}$ and ${\bf C}\:=\:({\bf K}^{-1}\:+\:\bfSigma^{-1})^{-1}$. Thus, we see that ${\bf W}\:\leftrightarrow\:\bfSigma^{-1}$ and $w_{ii}$ with Laplace approximation is equivalent to $1 \over {\sigma^2_i}$ in EP approximation. From (\ref{cavmu}), we have $\sigma^2_{\setminus{i}}\:=\:(({\bf C}_{ii})^{-1}\:-\:\sigma^{-2}_i)^{-1}$. Let us denote variances of posterior marginals with Laplace approximation as ${\tilde \sigma^2_i}$; then, ${\tilde \sigma^2_i}\:=\:[({\bf K}^{-1}+{\bf W})^{-1}]_{ii}$. Note that ${\bf C} \leftrightarrow ({\bf K}^{-1}+{\bf W})^{-1}$; thus, ${\tilde \sigma^2_i}\:\leftrightarrow\:{\bf C}_{ii}$. Therefore, given the equivalences between $w_{ii}$ and $1 \over {\sigma^2_i}$ and, between ${\tilde \sigma^2_i}$ and ${\bf C}_{ii}$, we can write LOO-CV variance for the Laplace approximation case as: 
%\begin{equation}
%\sigma^2_{\setminus{i}}\:=\:{{\tilde \sigma^2_i} \over {1-w_{ii}{\tilde \sigma^2_i}}}
%\label{lapvarexp}
%\end{equation}
%To find $\mu_{\setminus{i}}$ we note that ${\hat {\bf f}}\:=\:{\bf K}\bfalp \leftrightarrow {\bf m}$. Further, $\bfmu\:=\:\bfSigma {\bf C}^{-1} {\bf m}\:=\:\bfSigma ({\bf K}^{-1}+\bfSigma^{-1}){\bf m}\:=\:\bfSigma {\bf K}^{-1}{\bf m}\:+\:{\bf m}.$  Thus, we find $\mu_i \leftrightarrow {{\alpha_i} \over {w_{ii}}}\:+\:{\hat f}_i$. This implies ${{\mu_i} \over {\sigma^2_i}} \leftrightarrow \alpha_i\:+\:w_{ii} {\hat f}_i$. From (\ref{cavmu}), we have $\mu_{\setminus{i}}\:=\:\sigma^2_{\setminus{i}}\bigl({m_i \over {\bf C}_{ii}}-{\mu_i \over \sigma^2_i}\bigr)$. With the equivalence of ${\tilde \sigma^2_i}\:\leftrightarrow\:{\bf C}_{ii}$ noted earlier and the above expression for $\sigma^2_{\setminus{i}}$, we have
%\begin{equation}
%\mu_{\setminus{i}}\leftrightarrow {\hat f}_i\:-\:\alpha_i \sigma^2_{\setminus{i}}
%\label{lapmuexp}
%\end{equation}
%Note that the LOO-CV mean is obtained by subtracting a correction term from ${\hat f}$ that was obtained from the complete set. The LOO-CV variance is obtained by scaling the variance of the posterior marginal ${\tilde \sigma^2_i}$. Also, we find that LOO-CV mean and variance can be computed from the mode ${\hat {\bf f}}$, diagonal elements of the matrix ${\bf W}$ and variances of the posterior marginals obtained from the diagonal elements of the matrix $({\bf K}^{-1}\:+\:{\bf W})^{-1}$. This completes the specification of predictive distributions $p(f_i|x_i,{\it S},\bftheta)$. 

\subsection*{Derivatives of Predictive Distributions}
For ease of reference we recall (\ref{Pyi}) here. $$p(y_i|x_i,{\it S}_{\setminus{i}},\bftheta)\:=\:\Phi \Bigl({{y_i(\mu_{\setminus{i}}+\gamma)} \over {\sqrt{1+\sigma^2_{\setminus{i}}}}}\Bigr).$$ Then, with $z_i\:=\:{{\mu_{\setminus{i}}+\gamma} \over {\sqrt{1+\sigma^2_{\setminus{i}}}}}$ and $N(z_i)\:=\:{1 \over \sqrt{2\pi}} \exp(-{{z^2_i} \over 2})$ we have $${{\partial p(y_i|{\bf x}_i,{\it S}_{\setminus{i}},\bftheta)} \over {\partial \theta_j}}\:=\:{{{\it N}(z_i) y_i} \over {\sqrt{1+\sigma^2_{\setminus{i}}}}}\Bigl({{\partial \mu_{\setminus{i}}}\over {\partial \theta_j}}\:-\:{1\over 2}{{z_i} \over  {\sqrt{1+\sigma^2_{\setminus{i}}}}} {{\partial \sigma^2_{\setminus{i}}} \over {\partial \theta_j}}\Bigr).$$ Here, $\theta_j$ represents any element of $\bftheta$ other than $\gamma$. Similarly, we have  
%$${{\partial p(y_i|{\bf x}_i,{\it S}_{\setminus{i}},\bftheta)} \over {\partial \gamma_1}}\:=\:{{{\it N}(z_i) y_i} \over {\sqrt{1+\gamma^2_1\sigma^2_{\setminus{i}}}}}\Bigl(\mu_{\setminus{i}}\:-\:{{{\gamma_1 z_i \sigma^2_{\setminus{i}}}} \over  {\sqrt{1+\gamma^2_1\sigma^2_{\setminus{i}}}}}\Bigr).$$ 
$${{\partial p(y_i|{\bf x}_i,{\it S}_{\setminus{i}},\bftheta)} \over {\partial \gamma}}\:=\:{{{\it N}(z_i) y_i} \over {\sqrt{1+\sigma^2_{\setminus{i}}}}}.$$  
Thus, we need ${\partial \mu_{\setminus{i}}}\over {\partial \theta_j}$ and ${\partial \sigma^2_{\setminus{i}}} \over {\partial \theta_j}$. Below, we give details on how they can be obtained with the EP approximation.    

%\subsubsection*{\bf Derivatives of LOO mean and variance: Laplace Approximation}
 %From (\ref{lapvarexp}) and (\ref{lapmuexp}) we see that their derivatives can be found, if we have ${\partial {\hat {\bf f}}} \over {\partial \theta_j}$, ${\partial {\bf W}} \over {\partial \theta_j}$ and ${\partial {\tilde \sigma^2_{i}}} \over {\partial \theta_j}$; note that ${\bf W}$ is a diagonal matrix and ${\tilde \sigma^2_{i}}\:=\:[({\bf K}^{-1}+{\bf W})^{-1}]_{ii}$. From stable computations perspective (Rasmussen and Williams, 2006), it is useful to rewrite $({\bf K}^{-1}+{\bf W})^{-1}\:=\:{\bf K}\:-\:{\bf K} {\bf W}^{1 \over 2}{\bf B}^{-1}{\bf W}^{1\over 2}{\bf K}$ where ${\bf B}\:=\:({\bf I}+{\bf W}^{1\over2}{\bf K}{\bf W}^{1\over2})$ using Sherman-Morrison-Woodbury formula and this representation avoids computation of ${\bf K}^{-1}$. As a result, we have ${\tilde \sigma^2_{i}}\:=\:{\bf K}(i,i)-{\bf k}({\bf x}_i)^T{\bf W}^{1\over 2}{\bf B}^{-1}{\bf W}^{1\over 2}{\bf k}({\bf x}_i)$. Note that ${{\partial {\bf B}^{-1}} \over {\partial \theta_j}}\:=\:-{\bf B}^{-1}{{\partial {\bf B}} \over {\partial \theta_j}}{\bf B}^{-1}$. Thus, to find ${\partial {\tilde \sigma^2_{i}}} \over {\partial \theta_j}$ we only need ${\partial {\bf K}} \over {\partial \theta_j}$ and ${\partial {\bf W}} \over {\partial \theta_j}$ and expressions for these derivatives are available in (Rasmussen and Williams, 2006).

\subsection*{Derivatives of LOO mean and variance with fixed site parameters: EP Approximation}
In the case of EP approximation, since we resort to EM type optimization, we derive expressions for ${{\partial \mu_{\setminus{i}}} \over {\partial \theta_j}}$ and ${{\partial \sigma^2_{\setminus{i}}} \over {\partial \theta_j}}$ assuming that the site parameters are fixed. From $\mu_{\setminus{i}}\:=\:\sigma^2_{\setminus{i}}\bigl({m_i \over {\bf C}_{ii}}-{\mu_i \over \sigma^2_i}\bigr)$, we have
$${{\partial \mu_{\setminus{i}}} \over {\partial \theta_j}}\:=\:{{\mu_{\setminus{i}}} \over {{\sigma^2}_{\setminus{i}}}} {{\partial \sigma^2_{\setminus{i}}} \over {\partial \theta_j}}\:+\:{{\sigma^2_{\setminus{i}}} \over {({\bf C}_{ii})^2}} \Bigl({\bf C}_{ii}{{\partial {\bf m}_i} \over {\partial \theta_j}}\:-\:{\bf m}_i {{\partial {\bf C}_{ii}} \over {\partial \theta_j}}\Bigr).$$ From $\sigma^2_{\setminus{i}}\:=\:(({\bf C}_{ii})^{-1}\:-\:\sigma^{-2}_i)^{-1}$, we have
$${{\partial \sigma^2_{\setminus{i}}} \over {\partial \theta_j}}\:=\: {{\sigma^4_{\setminus{i}}} \over {({\bf C}_{ii})^2}} {{\partial {\bf C}_{ii}} \over {\partial \theta_j}}.$$ Since ${\bf m}\:=\:{\bf C}{\bfSigma^{-1}}{\bfmu}$, we have ${{\partial {\bf m}} \over {\partial \theta_j}}\:=\:{{\partial {\bf C}} \over {\partial \theta_j}} {\bfSigma^{-1}}{\bfmu}$. Note that ${\bf C}$ can be re-written using Sherman-Morrison-Woodbury formula as: ${\bf C}\:=\:{\bf K}\:-\:{\bf K}({\bf K}\:+\:\bfSigma)^{-1}{\bf K}$ and it is useful to work with this expression to achieve improved numerical stability \citep{Raswil:06} as it avoids inversion of ${\bf K}$. Then we have 
$${{\partial {\bf C}} \over {\partial \theta_j}}=\:({\bf I}\:-\:({\bf K}+\bfSigma)^{-1}{\bf K})^T{{\partial {\bf K}} \over {\partial \theta_j}}({\bf I}\:-\:({\bf K}+\bfSigma)^{-1}{\bf K}).$$ Note that ${{\partial {\bf C}_{ii}} \over {\partial \theta_j}}$, $i \in {\tilde I}$ (where ${\tilde I}=\{1,2,\ldots,n\}$) are nothing but the diagonal entries of the above expression. Note also that ${{\partial {\bf m}} \over {\partial \theta_j}}$ can be efficiently computed by taking advantage of the presence of the vector ${\bfSigma^{-1}}{\bfmu}$. But, to compute ${{\partial {\bf C}_{ii}} \over {\partial \theta_j}}$, $i \in {\tilde I}$ we cannot avoid the matrix multiplication with ${{\partial {\bf K}} \over {\partial \theta_j}}$; this results in $O(n^3)$ for each $\theta_j$. Finally, it is useful to re-write $({\bf K}+\bfSigma)^{-1}\:=\:\bfSigma^{-{1\over2}}({\bf I}\:+\:\bfSigma^{-{1\over2}}{\bf K}\bfSigma^{-{1\over2}})^{-1}\bfSigma^{-{1\over2}}$ \citep{Raswil:06}.     
%$${{\partial {\bf C}} \over {\partial \theta_j}}=\:{{\partial {\bf K}} \over {\partial \theta_j}}({\bf I}\:-\:({\bf K}+\bfSigma)^{-1}{\bf K})\:-\:({\bf K}+\bfSigma)^{-1}{\bf K}{{\partial {\bf K}} \over {\partial \theta_j}}({\bf I}\:-\:({\bf K}+\bfSigma)^{-1}{\bf K}).$$   

          




cvgpcv16.tex
\documentclass[twoside,11pt,graphicx,subfigure]{article}

% Any additional packages needed should be included after jmlr2e.
% Note that jmlr2e.sty includes epsfig, amssymb, natbib and graphicx,
% and defines many common macros, such as 'proof' and 'example'.
%
% It also sets the bibliographystyle to plainnat; for more information on
% natbib citation styles, see the natbib documentation, a copy of which
% is archived at http://www.jmlr.org/format/natbib.pdf

\usepackage{jmlr2e}

% Definitions of handy macros can go here

\newcommand{\dataset}{{\cal D}}
\newcommand{\fracpartial}[2]{\frac{\partial #1}{\partial  #2}}

% Heading arguments are {volume}{year}{pages}{submitted}{published}{author-full-names}

%\jmlrheading{}{}{}{}{}{S. Sundararajan and S. Sathiya Keerthi}
%\jmlrheading{S. Sundararajan and S. Sathiya Keerthi}
%\jmlrheading{1}{2000}{1-48}{4/00}{10/00}{Marina Meil\u{a} and Michael I. Jordan}

% Short headings should be running head and authors last names

%\ShortHeadings{LOO CV based optimization criteria for Gaussian process classifier model selection}{S. Sundararajan and S. Sathiya Keerthi}
\ShortHeadings{Predictive approaches for Gaussian process classifier model selection}{S. Sundararajan and S. Sathiya Keerthi}
%\firstpageno{1}

\def\bfalp{{\mbox{\boldmath $\alpha$}}}
\def\bfgamma{{\mbox{\boldmath $\gamma$}}}
\def\bfSigma{{\mbox{\boldmath $\Sigma$}}}
\def\bfmu{{\mbox{\boldmath $\mu$}}}
\def\bfnu{{\mbox{\boldmath $\nu$}}}
\def\bfbeta{{\mbox{\boldmath $\beta$}}}
\def\bftheta{{\mbox{\boldmath $\theta$}}}
\def\bfeta{{\mbox{\boldmath $\eta$}}}
\def\bfzeta{{\mbox{\boldmath $\zeta$}}}
\def\kuu{{{\bf K}_{{\bf u},{\bf u}}}}
\def\kuf{{{\bf K}_{{\bf u},{\bf f}}}}
\def\kfu{{{\bf K}_{{\bf f},{\bf u}}}}
\def\kui{{\bf K}_{{\bf u},i}}
\def\kiu{{\bf K}_{i,{\bf u}}}
\def\kuj{{\bf K}_{{\bf u},j}}
\def\kju{{\bf K}_{j,{\bf u}}}
\def\bfAu{{\bf A}_{\bf u}}
\def\muii{mu_{-i}(i)}
\def\bfu{{\bf u}}
\def\bfl{\bf L}
\def\bfui{{\bar \bf u}_{i}}
\def\bfuj{{\bar {\bf u}}_{j}}
\def\bful{{\bar {\bf u}}_{l}}
\def\bfAu{{\bf A}_{\bfu}}
\def\bfinvGu{{\bf G}^{-1}_{\bfu}}
\def\bfGu{{\bf G}_{\bfu}}
\def\bfinvGut{{\bf G}^{-T}_{\bfu}}
\def\bfGut{{\bf G}^T_{\bfu}}
\def\bfinvGuj{{\bf G}^{-1}_{\bfuj}}
\def\bfGujt{{\bf G}^T_{\bfuj}}
\def\bfGuj{{\bf G}_{\bfuj}}
\def\bfinvGujt{{\bf L}^{-T}_{\bfuj}}
\def\bfLujt{{\bf L}^T_{\bfuj}}
\def\bfinvLu{{\bf L}^{-1}_{\bfu}}
\def\bfLu{{\bf L}_{\bfu}}
\def\bfinvLut{{\bf L}^{-T}_{\bfu}}
\def\bfLut{{\bf L}^T_{\bfu}}
\def\bfinvLuj{{\bf L}^{-1}_{\bfuj}}
\def\bfLuj{{\bf L}_{\bfuj}}
\def\bfinvLujt{{\bf L}^{-T}_{\bfuj}}
\def\bfLujt{{\bf L}^T_{\bfuj}}
\def\bfAuj{{\bf A}_{\bfuj}}
\def\kuuj{{{\bf K}_{{\bfuj},{\bfuj}}}}
\def\bfzjt{{\bf z}^T_j}
\def\bfzlt{{\bf z}^T_l}
\def\bfhjt{{\bf h}^T_j}
\def\bfhj{{\bf h}_j}
\def\bfhi{{\bf h}_i}
\def\bfzj{{\bf z}_j}
\def\bfzl{{\bf z}_l}
\def\bfzi{{\bf z}_i}
\def\dj{d_j}
\def\di{d_i}
\def\ej{e_j}
\def\ei{e_i}
\def\kfl{{\bf K}_{{\bf f},l}}
\def\kfi{{\bf K}_{{\bf f},i}}
\def\kfj{{\bf K}_{{\bf f},j}}
\def\kfl{{\bf K}_{{\bf f},l}}
\def\kif{{\bf K}_{i,{\bf f}}}
\def\kjf{{\bf K}_{j,{\bf f}}}
\def\klf{{\bf K}_{l,{\bf f}}}
\def\bfbjt{{\bf b}^T_j}
\def\bfblt{{\bf b}^T_l}
\def\cj{c_j}
\def\cl{c_l}
\def\bfbj{{\bf b}_j}
\def\bfbl{{\bf b}_l}
\def\bfwjt{{\bf w}^T_j}
\def\bfwj{{\bf w}_j}
\def\bfwi{{\bf w}_i}
\def\kii{{\bf K}_{i,i}}
\def\kil{{\bf K}_{i,l}}
\def\kij{{\bf K}_{i,j}}
\def\kji{{\bf K}_{j,i}}
\def\kli{{\bf K}_{l,i}}
\def\kjj{{\bf K}_{j,j}}
\def\kll{{\bf K}_{l,l}}
\def\bfy{{\bf y}}
\def\uj{u_j}
\def\ui{u_i}
\def\etaj{\eta_j}
\def\etai{\eta_i}
\def\bfs{\bf S}
\def\bfsi{{\bf S}_{-i}}


\begin{document}

\title{Predictive Approaches For Gaussian Process Classifier Model Selection \\ September 10, 2008}
%\title{Leave One Out Cross Validation Based Optimization Criteria For Gaussian Process Classifier Model Selection}

\author{\name S. Sundararajan \email ssrajan@yahoo-inc.com \\
       \addr Yahoo! Labs\\
       Bangalore, India\\
       \AND
       \name S. Sathiya Keerthi \email selvarak@yahoo-inc.com \\
       \addr Yahoo! Research \\
       Santa Clara, USA \\
       }

\editor{}

\maketitle

\begin{abstract}%   <- trailing '%' for backward compatibility of .sty file
In this paper we consider the problem of Gaussian process classifier (GPC) model selection with different Leave-One-Out (LOO) Cross Validation (CV) based optimization criteria and provide a practical algorithm using LOO predictive distributions with such criteria to select hyperparameters. Apart from the standard average negative logarithm of predictive probability (NLP), we also consider smoothed versions of criteria such as F-measure and Weighted Error Rate (WER), which are useful for handling imbalanced data. Unlike the regression case, LOO predictive distributions for the classifier case are intractable. We 
%Keerthi:
use 
approximate LOO predictive distributions arrived from Expectation Propagation (EP) approximation. We conduct experiments on several real world benchmark datasets. When the NLP criterion is used for optimizing the hyperparameters, the predictive approaches show better or comparable NLP generalization performance with existing GPC approaches. On the other hand, when the F-measure criterion is used, the F-measure generalization performance improves significantly on several datasets. Overall, the EP-based predictive algorithm comes out as an excellent choice for GP classifier model selection with different optimization criteria.

\end{abstract}

\begin{keywords}
  Gaussian process classification, Model Selection, LOO, Cross Validation, Predictive distributions, Smoothed F-measure, Weighted Error Rate, Precision, Recall, Imbalanced data
\end{keywords}

\section{Introduction}
Gaussian process (GP) models are flexible and powerful probabilistic models that 
%Keerthi:
are
used to solve classification problems in many areas of application \citep{Raswil:06}. In 
%Keerthi:
the
Bayesian GP setup, latent function values and hyperparameters that are involved in modeling are integrated with chosen priors. However, the required integrals are often not analytically tractable (due to various choices of likelihoods and priors) and closed form analytic expressions are not available. Two important problems in this context are finding good approximations for integrating over the latent function variables and the hyperparameters. There have been two approaches reported in the literature. In the first approach, both the latent function variables and the hyperparameters are integrated out within some approximations. \citet{Williams:98} used Laplace approximation to integrate over the latent function variables and Hybrid Monte Carlo (HMC) to integrate over the hyperparameters. \citet{Neal:98} used Gibbs sampling to integrate over the latent function variables and HMC to integrate over hyperparameters; this method is more accurate, but it is computationally expensive. In the second approach, only the latent function variables are integrated out and the hyperparameters are optimized on some well-defined objective function. This latter problem of choosing hyperparameters that define the model is essentially the Gaussian process model selection problem and in this paper we focus on this problem.  

%Gaussian process model selection is an important problem and 
%Keerthi:
%usually
%occurs in the form of choosing hyperparameters (instead of integrating over) that define the model. In practice, the choice is made by optimizing a well-defined objective function over the hyperparameters. 
There are two commonly used approaches to address this model selection problem. They are marginal likelihood or evidence maximization and minimization of LOO-CV based average negative logarithmic predictive probability (NLP). Both these approaches are available for GP regression model selection 
%Keerthi1:
%Include Sundararajan & Keerthi paper?
~\citep{Raswil:06,Sundar:01}. For GP classifier model selection, \citet{Gibbs:00} used a variational approximation method to integrate over the latent function values and estimated the hyperparameters by maximizing marginal likelihood (ML). 
Laplace approximation and Expectation Propagation (EP) approximation ~\citep{Raswil:06,Seeger:03} are other methods that have been used to integrate over the latent function variables. Then, the marginal likelihood is optimized using gradient information obtained with any one of these approximations \citep{Raswil:06,Seeger:03}. \citet{Kim:06} presented an approximate Expectation-Maximization (EM) algorithm to learn the hyperparameters. In the E-step, they used EP to estimate the joint density of latent function values and in the M-step, the hyperparameters were optimized by maximizing 
%Keerthi:
a variational lower bound on the marginal likelihood. 
%Also, the LOO-CV model selection method for the regression problem was extended to the GP classifier model selection problem (in the least squares framework) to develop the LOO-CV based probabilistic least squares (LS) classifiers (Rasmussen and Williams, 2006).     

In this paper, we consider the approach of using LOO-CV based predictive distributions to address the GP classifier model selection problem. In a related work, \citet{Opper:00} used LOO error estimate to choose rough hyperparameter values by scanning a range of values. In the EP framework \citep{Minka:01}, cavity distributions directly provide LOO error estimates during training and were used to assess predictive performance and select automatic relevance determination (ARD) type hyperparameters \citep{Qi:04}.  The LOO-CV based predictive distributions obtained from probabilistic least squares classifier were used in the minimization of NLP for GP classifier model selection \citep{Raswil:06}. 

In practice, while measures like marginal likelihood and average negative 
%Keerthi1: logarthmic 
logarithmic
predictive probability measures are very useful, other measures like F-measure \citep{van:74} and Weighted Error Rate (WER) are also important and useful, for instance in applications like medical diagnostics, image understanding, etc, where the number of positive examples 
%Keerthi:
is
much 
%Keerthi:
smaller
than the number of negative examples. 
%Keerthi1:
Several works that use such measures for hyperparameters optimization exist in the
non-GP literature.
\citet{Hong:07} proposed 
%Keerthi:
a
kernel classifier construction algorithm based on regularized orthogonal weighted least squares (ROWLS) estimation with LOO-Area Under the ROC Curve (AUC) as model selection criterion for handling imbalanced datasets. \citet{Jansche:05} proposed 
%Keerthi:
the
training of 
%Keerthi:
a
probabilistic classifier based on a logistic regression model by optimizing expected F-measure. Here an approximation to F-measure was made so that the F-measure is smoothed and becomes a smooth function of model weights.  
%Keerthi:
%Seeger's work on KLR also needs to be mentioned here.
\citet{Seeger:07} proposed a general framework for learning in probabilistic kernel classification models with a large or structured set of classes. He optimized the kernel parameters by minimizing the NLP over k-folds.  \citet{Keerthi:06} considered the task of tuning hyperparameters in SVM models based on minimizing smooth performance validation functions like smoothed k-fold CV error. 
%Keerthi:
The last three works did not use LOO-CV in their work.
%Keerthi: I am commenting the material below given the above modification.
%In the latter two cases, non-linear optimization techniques can be used to optimize model weights and hyperparameters respectively. However, they did not take LOO-CV approach in their work. 

%Besides these works, we are not aware of any well-defined LOO-CV based GP classifier model selection algorithm that {\it optimizes} the hyperparameters with classification problem specific likelihood functions directly for different criteria.

%Keerthi1:
This paper is aimed at addressing two issues. First, the proposed method is different from the work of \citet{Opper:00} and \citet{Qi:04} in that we optimize the hyperparameters directly in the continuous space and any standard non-linear optimization method can be used. It is also different from the LOO-CV based probabilistic LS method \citep{Raswil:06} in that we use the more accurate EP approximation than the LS approximation. Second, criteria such as F-measure and WER, which are needed for tackling imbalanced problems, have not been considered in GP classifier designs.
 
%The main contributions of this paper are as follows:  while cross-validation approach for model selection is very useful (Geisser, 1975; Stone, 1974; Wahba, 1998) and not new, we provide practical algorithms using LOO predictive distributions for GP classifier model selection with different LOO-CV based optimization criteria. We extend the above mentioned work addressed from various aspects as follows: 
%Keerthi1:
We define smoothed LOO-CV based measures using predictive distributions as a function of GP classifier model hyperparameters. Thus, the objective functions can be optimized using standard non-linear optimization techniques. We investigate usage of LOO-CV predictive distributions obtained from expectation propagation approximation. Actually, the proposed algorithm can also be used with Laplace approximation. However, \citet{Kuss:05} showed for binary classification problems that the EP approximation is better than the Laplace approximation. Therefore, we restrict our attention to using the EP approximation here.  
%Keerthi1:
%With EP approximation, we resort to EM type approach of alternating between optimization of site parameters and hyperparameters while fixing the other. This results in a very useful algorithm which we call EP-CV algorithm. Relevant expressions are derived and we propose algorithms to select model based on these criteria, by optimizing the hyperparameters. Specifically, apart from using standard average NLP optimization, we show how LOO-CV based predictive distributions can be used to optimize smoothed versions of criteria like F-measure, WER and AUC etc. Though we demonstrate using smoothed F-measure optimization, the approach can be readily used with other measures. 
%Thus, for demonstration we consider two criteria; the first one is standard average negative logarithm of predictive probability and the second one is smoothed version of F-measure that is widely used in the context of imbalanced data. 
%Keerthi1:

We conduct experiments on two criteria: the standard average negative logarithm of predictive probability (NLP) and a smoothed version of F-measure.
On the NLP criterion we compare our method (EP-CV(NLP)) against the LOO-CV based probabilistic least squares (LS) classifier \citep{Raswil:06} and standard GP classifier doing ML maximization using EP approximation. We refer the latter two methods as LS-CV(NLP) and EP(ML) respectively; the abbreviations in parentheses refer to the type of objective functions used. The experimental results on several real world benchmark datasets show that the proposed method is better than the LS-CV(NLP) method and is quite competitive to the EP-ML method. On the F-measure criterion we compare the EP-CV(NLP) method with the EP-CV(FM) method. In the latter method we optimize over the F-measure instead of the NLP measure. We also compare with a two-step method\footnote{This method was suggested by an anonymous reviewer.} where, in the first step we optimize over the hyperparameters using the NLP measure and, in the second step we optimize {\it only} the {\it bias} hyperparameter using the F-measure. Experimental results demonstrate that this method is also inferior to the EP-CV(FM) method.   
%Keerthi1:
%As the experimental results show, the EP-CV algorithm is an excellent choice for GP classifier model selection with different LOO-CV based optimization criteria.      

%Keerthi1:
The 
paper is organized as follows. We give a brief introduction to Gaussian process classification and ML optimization criteria in Section 2. In 
%Keerthi1: comma removed
Section 3 we give 
%Keerthi1:
a general set of smooth
LOO-CV based optimization criteria and illustrate their use with LOO-CV predictive distributions. Optimization aspects and 
%Keerthi1:
specific algorithm are given in Section 4. In 
%Keerthi1: comma removed
Section 5 we discuss related work on LOO-CV based GPC model selection. In 
%Keerthi1: comma removed
Section 6 we present experimental results and 
%Keerthi1:
then conclude 
%Keerthi1:
the paper in Section 7. 
     
\section{Gaussian Process Classification}
In binary classification problems, we are given a training data set ${\it S}$ composed of $n$ input-target pairs $({\bf x}_i, y_i)$ where ${\bf x}_i$ $\in$ $R^D$, $y_i$ $\in$ $\{+1,-1\}$, 
%Keerthi: Why not use I instead of Itilde? ---- To avoid confusion with {\bf I} that appears for identity matrix some where!
$i \in {\tilde I}$ and ${\tilde I}=\{1,2,\ldots,n\}$. The true function value at ${\bf x}_i$ is represented as a latent variable $f({\bf x}_i)$. The goal is to compute the predictive distribution of the class label $y_*$ at test location ${\bf x}_*$. In standard GPs for classification \citep{Raswil:06}, the latent variables $f({\bf x}_i)$ are modelled as random variables in a zero mean GP indexed by $\{{\bf x}_i\}$. The prior distribution of $\{{\bf f}({\bf X}_n)\}$ is a zero mean multivariate joint Gaussian, denoted as ${\it p}({\bf f})\:=\:{\mathcal N}({\bf 0},{\bf K})$, where ${\bf f}\:=\:[f({\bf x}_1),\ldots,f({\bf x}_n)]^T$, ${\bf X}_n\:=\:[{\bf x}_1,\ldots,{\bf x}_n]$ and ${\bf K}$ is the $n \times n$ covariance matrix whose $(i,j)^{th}$ element is $k({\bf x}_i,{\bf x}_j)$, denoted as ${\bf K}_{i,j}$. One of the most commonly used covariance 
%Keerthi1:
functions
is the squared exponential covariance function given by: $\mathrm{cov}(f({\bf x}_i),f({\bf x}_j))\:=\:k({\bf x}_i,{\bf x}_j)\:=\:\beta_0\:\exp(-{1 \over 2} \sum_{k=1}^D {{(x_{i,k}-x_{j,k})^2} \over {\beta_k}})$. Here, $\beta_0$ represents signal variance and 
the remaining $\beta_k$'s represent width parameters across different input dimensions. Let $\bfbeta\:=\:=[\beta_0, \beta_1, \ldots, \beta_D]$. These parameters are also known as {\it automatic relevance determination} (ARD) hyperparameters. We call this covariance function the ARD Gaussian kernel function. 
%Keerthi: The sentence below is too long and can be rewritten in a better way. --- modified ---
Next, it is assumed that the probability over class labels as a function of ${\bf x}$ depends only on the latent function value $f({\bf x})$. For the binary classification problem, given the value of $f({\bf x})$ the probability of class label is independent of all other quantities: $p(y=+1|f({\bf x}),{\it S})\:=\:p(y=+1|f({\bf x}))$ where ${\it S}$ is the dataset. 
%Keerthi:
The likelihood $p(y_i|f_i)$ can be modelled in several forms 
%Keerthi1: like 
such as a sigmoidal function or cumulative normal $\Phi(y_i f_i)$ where 
%Keerthi: Below, use \sqrt{2\pi} instead of \sqrt(2\pi)  --- done ---
$\Phi(z)\:=\:\int_{-\infty}^z {1 \over {\sqrt{2\pi}}} \exp(-{{w^2} \over 2}) dw$. Assuming that the examples are i.i.d, we have ${\it p}({\bf y}|{\bf f})\:=\:\prod_{i=1}^N p(y_i|f_i;\bfgamma)$. Here, $\bfgamma$ represents hyperparameters that characterize the likelihood.  The prior and likelihood along with the hyperparameters $\bftheta \:=\:[\bfbeta, \bfgamma]$ characterize the GP model. With these modelling assumptions, we can write the 
%Keerthi1:
inference probability 
%Keerthi: Leave out "of interest"
given $\bftheta$ as:
\begin{equation}
	p(y_*|{\bf x}_*,{\it S},\bftheta)\:=\: \int \:p(y_*|f_*,\bfgamma) \:p(f_*|{\it S},{\bf x}_*,\bftheta)\:\:df_*
\label{pred1}
\end{equation}
Here, the posterior predictive distribution of latent function $f_*$ is given by:
	%$$p(f_*|{\it S},{\bf x}_*,\bftheta)\:=\:\int \: p({\bf f},f_*|{\it S}, {\bf x}_*, \bftheta)\: d{\bf f}$$
%This can be rewritten as:
%\begin{equation}
	$$p(f_*|{\it S},{\bf x}_*,\bftheta)\:=\:\int\:p(f_*|{\bf x}_*,{\bf f},\bfbeta) \:p({\bf f}|{\it S},\bftheta)\:d{\bf f}.$$
%\label{pred2}
%\end{equation}
where $p({\bf f}|{\it S},\bftheta) \varpropto \prod_{i=1}^N p(y_i|f_i,\bfgamma)\: p({\bf f}|{\bf X},\bfbeta)$. In a Bayesian solution, the class probability at the test point $x_*$ would be obtained by integrating over the hyperparameters weighted by their posterior probability 
$$ p(y_*|x_*,{\it S})\:=\:\int \:p(y_*|x_*,{\it S},\bftheta)\:p(\bftheta|{\it S}) \:d\bftheta.$$
In general there is no closed form expression available for this integral and 
%Keerthi1:
it
is expensive to compute. Therefore, instead of integrating over the hyperparameters, 
%Keerthi1: they 
a single set of their values
is  
%Keerthi1:
usually
estimated from the dataset by optimizing various criteria as mentioned earlier and 
%Keerthi1:
then
used in (\ref{pred1}). 
%Keerthi1:
%Here, we consider the approach of selecting the model by optimizing the hyperparameters using certain well-defined objective function. 

\subsection{Marginal Likelihood Maximization}
Marginal likelihood or evidence maximization \citep{Raswil:06} is commonly used to estimate the hyperparameters during model selection. The marginal likelihood is given by: 
\begin{equation}
	p({\bf y}|{\bf X},\bftheta)\:=\:\int\:\prod_{i=1}^N p(y_i|f_i,\bfgamma)\: p({\bf f}|{\bf X},\bfbeta)\:d{\bf f}
\label{ML}
\end{equation}
This integral cannot be calculated analytically except for a special case like GP regression with Gaussian noise. Therefore, certain approximations are needed to compute these quantities. Laplace approximation and EP approximations are two popular methods used for this purpose. To gain more insight into the quality of the Laplace and EP approximations, \citet{Kuss:05} carried out comprehensive comparisons of these approximations with (the more exact) Markov Chain Monte Carlo (MCMC) sampling approach in terms of their predictive performance and marginal likelihood estimates. They found that EP is the method to be used for approximate inference in binary GPC models, when the computational cost of MCMC is prohibitive. Hence in our study we restrict ourselves to the more accurate EP approximation given in the next subsection. With the EP approximation the hyperparameters are learnt by optimizing marginal likelihood with gradient information \citep{Raswil:06, Seeger:03} using standard non-linear optimization techniques. 
%Variational bound on the marginal likelihood (Gibbs and MacKay, 2000; Kim and Ghahramani, 2006) 
%Keerthi1: is 
%has also been used 
%optimized to learn the hyperparameters.

%\subsection{Approximations}
%\subsubsection{Laplace Approximation}
%Keerthi:
%The Laplace method of Williams and Barber (1998) consists of integrating out
%latent function variables using a Taylor expansion of 
%$\psi({\bf f}|\bftheta)\:\triangleq\:\log\:p({\bf y}|{\bf f},\bftheta)\:+\:\log\:p({\bf f}|\bftheta)\:\simeq\:\psi(\hat {\bf f}|\bftheta)\:-\:{1\over 2}({\bf f}-\hat {\bf f})^T{\bf A}({\bf f}-\hat {\bf f})$. Here, $\hat {\bf f}\:=\:\arg \max_{\bf f}\:\psi({\bf f}|\bftheta)$ and ${\bf A}\:=\:-\nabla\nabla\psi({\bf f}|\bftheta)\:=\:{\bf W}+{\bf K}^{-1}$ where ${\bf W}\:=\:\nabla\nabla\log\:p({\bf y}|{\bf f})$ is diagonal since the likelihood factorizes over individual examples. If the likelihood $p({\bf y}|{\bf f})$ is log concave, the diagonal elements are non-negative and the Hessian $-\nabla\nabla\psi({\bf f})$ is negative definite, so that $\psi({\bf f}|\bftheta)$ is concave and has a unique maximum. Thus, the Laplace approximation to the posterior is specified as Gaussian with mean ${\hat {\bf f}}$ and covariance $({\bf K}^{-1}+{\bf W})^{-1}$. With this approximation, the integral (\ref{ML}) can be evaluated analytically to obtain an approximation to the log marginal likelihood:
%\begin{equation}
%	\log q({\bf y}|{\bf X},\bftheta)\:=\:-{1 \over 2}{\hat {\bf f}}^T{\bf K}^{-1}{\hat {\bf f}}\:+\:\log\:p({\bf y}|\hat {\bf f}) - {1\over 2}\:\log|{\bf B}|
%\label{ALML}
%\end{equation}
%where $|{\bf B}|\:=\:|{\bf I}_n\:+\:{\bf W}^{1\over2}{\bf K}{\bf W}^{1\over2}|$ (Rasmussen and Williams, 2006) and 
%Keerthi: Leave out "a notation of"
%Keerthi: Use \cdot
%$q(\cdot)$ is used as an approximation to 
%Keerthi: Use \cdot
%$p(\cdot)$. One problem with the Laplace approximation is that 
%Keerthi:
%the
%Hessian evaluated at $\hat {\bf f}$ may give a poor approximation to the true shape of the posterior (Kuss and Rasmussen, 2005). Hence in our study we restrict to the more accurate EP approximation given below.  

\subsection{Expectation Propagation}
The EP algorithm is an iterative algorithm which is used for approximate Bayesian inference \citep{Opper:00,Minka:01}. It has been applied to GP classification \citep{Raswil:06, Seeger:03} 
%Keerthi1: Why does the Seeger reference appear here? --- EP algorithm for GP classification is discussed at length in his thesis and some proof regarding marginal likelihood derivatives is available only in his thesis ----
. EP finds a Gaussian approximation $q({\bf f}|{\it S},\bftheta)\:=\:\mathcal{N}({\bf f}|{\bf m},{\bf C})$ to the posterior $p({\bf f}|{\it S},\bftheta)$ by moment matching of approximate marginal distribution and the posterior. Mean and covariance of the approximate Gaussian are given by: 
\begin{equation}
{\bf m}\:=\:{\bf C}\bfSigma^{-1}{\bfmu} \:\:\:\:\:\: {\bf C}\:=\:({\bf K}^{-1}\:+\:\bfSigma^{-1})^{-1}
\label{posmean}
\end{equation}
where 
%Keerthi: mu_2
${\bfmu}\:=\:(\mu_1,\mu_2,\ldots,\mu_n)^T$ and $\bfSigma\:=\:diag(\sigma^2_1,\sigma^2_2,\ldots,\sigma^2_n)$ are called site function parameters. As per the approximation, the posterior is written in terms of the site functions $t(f_i;\mu_i,\sigma^2_i,Z_i)\:=\:Z_i {\mathcal N}(f_i|\mu_i,\sigma^2_i)$ and prior $p({\bf f}|{\bf X},\bftheta)$ as 
$$ q({\bf f}|{\it S},\bftheta)\:=\:{{p({\bf f}|{\bf X},\bftheta)} \over q({\it S}|\bftheta)} \prod_{i=1}^N\:t(f_i;\mu_i,\sigma^2_i,Z_i).$$ The EP algorithm iteratively visits each site function in turn, and adjusts the site parameters to match moments of an approximation to the posterior marginals. This process requires replacement of intractable exact {\it cavity} distribution with a tractable approximation based on the site functions and is given by:
$$q_{\setminus{i}}(f_i)\:=\:\int \prod_{j\neq i} t(f_j;\mu_j,\sigma^2_j,Z_j) \:p({\bf f}|{\bf X},\bftheta)\:d{\bf f}^{\setminus{i}}$$ where $q_{\setminus{i}}(f_i) \varpropto \mathcal{N}(f_i|\mu_{\setminus{i}},\sigma^2_{\setminus{i}})$ is the approximate cavity function and is related to 
%Keerthi:
the
diagonal entries of the posterior $q({\bf f}|{\it S},\bftheta)$ as:
$ q_{\setminus{i}}(f_i)t(f_i;\mu_i,\sigma^2_i,Z_i) \:\varpropto \:\mathcal{N}(f_i|m_i,{\bf C}_{ii}).$ Here, 
%Keerthi:
the
${\bf C}_{ii}$ 
%Keerthi:
represent 
the diagonal entries of the matrix ${\bf C}$. Using Gaussian identities, the mean and variance of cavity distribution is related to the site parameters as:
\begin{equation}
	\mu_{\setminus{i}}\:=\:\sigma^2_{\setminus{i}}\bigl({m_i \over {\bf C}_{ii}}-{\mu_i \over \sigma^2_i}\bigr) \:\:\:\:\:\:\:\:\:\sigma^2_{\setminus{i}}\:=\:\bigl(({\bf C}_{ii})^{-1}-\sigma^{-2}_{i}\bigr)^{-1}
\label{cavmu}
\end{equation}              
Then, EP adjusts the site parameters $\mu_i$, $\sigma^2_i$ and $Z_i$ such that the approximate posterior marginal using the exact likelihood approximates the posterior marginal based on the site function well. That is, $q_{\setminus{i}}(f_i)\:p(y_i|f_i) \simeq q_{\setminus{i}}(f_i)\:t(f_i;\mu_i,\sigma^2_i,Z_i)$. This is done by matching the zeroth, first and second moments on both sides. Thus, the EP algorithm iteratively updates site parameters until convergence.  
%Keerthi: This LOO comment doesn't fit well here. --------- as discussed ---- 
%and LOO predictive distributions are available during training as well. 
Though there is no convergence proof, 
%Keerthi:
in practice
the EP algorithm converges 
%Keerthi: Instead of "in practice"
in most cases. See \citet{Raswil:06} for more details. 
%The moments of the left hand side can be computed analytically in the case of probit (cumulative Gaussian) likelihood function $p(y_i|f_i)\:=\:\Phi(y_if_i)$ and are given by $$m_0\:=\:\Phi(z),$$  $$m_1\:=\:\mu_{\setminus{i}}+{{\eta_i \phi(z)} \over y},$$ $$m_2\:=\:2\mu_{\setminus{i}}m_1-\mu^2_{\setminus{i}}+\sigma^2_{\setminus{i}}-{z\phi(z)\eta^2_i},$$ where $z={{y\mu_{\setminus{i}}} \over {\sqrt{1+\sigma^2_{\setminus{i}}}}}$, $\phi(z)\:=\:{{\mathcal{N}(z)} \over {\Phi(z)}}$ and $\eta_i\:=\:{{\sigma^2_{\setminus{i}}} \over {\sqrt{1+\sigma^2_{\setminus{i}}}}}$. Next, by equating the moments of the right hand side, the update equations for the site parameters turn out to be: 
%$$\sigma^2_i\:=\:((m_2-m_1)^{-1}-\sigma^2_{\setminus{i}})^{-1},$$ $$\mu_i\:=\:\sigma^2_i(m_1(\sigma^{-2}_{\setminus{i}}+\sigma^{-2}_{i})-{{\mu_{\setminus{i}}} \over {\sigma^2_{\i}}},$$ $$Z_i\:=\:m_0\sqrt{2\pi(\sigma^2_{\i}+\sigma^2_{i})}exp\Bigl({{(\mu_i-\mu_{\setminus{i})^2}} \over {2(\sigma^2_{\setminus{i}}+\sigma^2_{i})}}\Bigr).$$ 
%Thus, the EP algorithm iteratively updates site parameters using the above expressions until convergence. Though there is no convergence proof, the EP algorithm converges in practice.

Next, within some constant, the marginal likelihood with EP approximation \citep{Raswil:06} is given by: 
\begin{equation}
	%\log q({\bf y}|{\bf X},\bftheta)\:=\:\-{1\over2}\log|{\bf K}+\bfSigma|\:-\bfmu^T({\bf K}+\bfSigma)^{-1}\bfmu\:+G(\bftheta)\:+\:{1\over 2}\sum_{i=1}^n\log(\sigma^2_{\setminus{i}}+\sigma^2_i)\:+\:\sum_{i=1}^n{{(\mu_{\setminus{i}}-\mu_i)^2} \over {2(\sigma^2_{\setminus{i}}+\sigma^2_i)}}
	\log q({\bf y}|{\bf X},\bftheta)\:=\:-{1\over2}\log|{\bf K}+\bfSigma|\:-{1\over 2}\bfmu^T({\bf K}+\bfSigma)^{-1}\bfmu\:+\:\sum_{i=1}^n\log w_i
\label{EPML}
\end{equation}
where $w_i=\Phi(z_i)\exp\Bigl({{(\mu_{\setminus{i}}-\mu_i)^2} \over {2(\sigma^2_{\setminus{i}}+\sigma^2_i)}}\Bigr)\sqrt{\sigma^2_{\setminus{i}}+\sigma^2_i}$ and $z_i={{y_i \mu_{\setminus{i}}} \over \sqrt{{1+\sigma^2_{\setminus{i}}}}}$. 
%Keerthi: Leave out "Finally,"
The 
hyperparameters are optimized using gradient expressions with standard conjugate gradient or quasi-Newton type non-linear optimization techniques.  

\section{Leave-One-Out Cross Validation based Optimization Criteria}
In this section, we 
%Keerthi: Leave out "first"
give definitions of various LOO-CV based optimization criteria. 
%Keerthi: Instead of "Next"
In section 4
we give details on how these measures can be optimized using standard nonlinear optimization 
%Keerthi:
techniques. 
%Keerthi:
%Leave out "here, we assume that"
%Keerthi:
The LOO predictive distributions  $p(y_i|{\bf x}_i,{\it S}_{\setminus{i}},\bftheta)$, $i\in\tilde I$ 
%Keerthi: Leave out
%as they define the measures and their derivatives with respect to the hyperparameters are available within some approximations. Note that 
%Keerthi:
play a crucial role. 
%Keerthi:
Here
${\it S}_{\setminus{i}}$ represents the dataset without $i$th example.
%Keerthi: Leave out
% and we consider these distributions and their derivatives obtained from Laplace and EP approximations. Finally, we give expressions and some details for the LOO predictive distributions; we present more details on expressions for gradients and derivations in the appendix.   
%Keerthi:
Their exact computation is expensive. In section 4 we will also discuss how to approximate them efficiently.

%\subsection{EP-CV Algorithm}
%Given the predictive distributions with EP algorithm, their gradients with respect to the hyperparameters $\bftheta$ are required to use with gradient based non-linear optimization techniques. In the case of marginal likelihood, closed form expressions are fortunately available for the derivatives with only explicit derivatives (Rasmussen and Willam, 2006; Seeger, 2003). But, with only predictive distributions, closed form expressions are not available due to iterative nature of the EP algorithm associated with implicit derivatives arising with site parameters. Thus, we resort to expectation-maximization (EM) type approach. In the E-step, the hyperparameters are fixed and the site parameters are estimated using EP algorithm; in the M-step, the site parameters are fixed and the hyperparameters are estimated by optimizing LOO-CV based optimization criterion. This procedure is repeated until convergence. We call this algorithm EP-CV algorithm and can be used with any LOO-CV based optimization criterion. It works well and we have never faced any problem with convergence in all our experiments. Note that the EM approach was used by Kim and Ghahramani (2006) (which they called EM-EP algorithm) in the optimization of a lower bound on marginal likelihood.           

%\subsection{Optimization Criteria}
\subsection{NLP Measure}
The averaged negative logarithm of predictive probability (NLP) is defined as: 
\begin{equation}
	G(\bftheta)\:=\:-{1\over n}\sum_{i=1}^n \:\log\:p(y_i|x_i,{\it S}_{\setminus{i}},\bftheta)
\label{nlp}
\end{equation}
%Keerthi: Leave out
%Thus, NLP measure is LOO-CV based criterion and we can select the hyperparameters by optimizing this criterion. 
This 
%Keerthi:
LOO-CV based
measure is generic and has been used in the context of probabilistic LS classifiers 
%Keerthi:
(see section 5)
%Keerthi1:
%Include Sundararajan & Keerthi paper?
and GP regression \citep{Raswil:06, Sundar:01}. 

%In probabilistic LS classification, the posterior mean and variance are obtained by solving the classification problem as if it is a regression problem. Then, the LOO mean and variance are obtained efficiently from LOO GP regression formulation; they are substituted in (\ref{Pyi}) and finally, (\ref{nlp}) is used to optimize the hyperparameters. 
     
\subsection{Smoothed LOO Measures}
While 
%Keerthi: Leave out "the above"
%Keerthi:
%Leave out: measure is very useful like 
%Keerthi:
measures such as
marginal likelihood (\ref{EPML}) 
%Keerthi:
and NLP in (\ref{nlp}) are useful for normal situations, 
other measures like F-measure 
%Keerthi: etc 
and WER  
are important, for example, 
%Keerthi:
when dealing with
imbalanced datasets. 
%Keerthi: These 
Let us now define these measures.
%Keerthi: are defined as follows: c
%Keerthi:

Consider the binary classification problem with class labels +1 and -1. Assume that there are $n_+$ positive examples and $n_{-}$ negative examples. In general, the performance of the classifier may be evaluated using counts of data samples $\{a, b, c, d\}$ defined via the confusion matrix given in Table 1. 
%Keerthi: Also, 
Let $n_{+}\:=\:a+b$ and $n_{-}\:=\:c+d$. The true positive rate (TP) is the proportion of positive data samples that were correctly classified (that is, true positives) and the false positive rate (FP) is the proportion of negative data samples that were incorrectly classified (that is, false positives) \citep{Hong:07}. These rates are given by: TP$\:=\:{a \over {a+b}}\:=\:{a \over n_{+}}$ and FP$\:=\:{c \over{c+d}}\:=\:{c \over {n_{-}}}$. 
%Keerthi: Thus, t
The misclassification rate is given by: MCR$\:=\:{{b+c} \over n}$. Note that the true positive rate is also known as {\it Recall} (R). {\it Precision} is another important quantity defined as: P$\:=\:{a \over {a+c}}$. 

Now let us consider the imbalanced data case and assume that $n_{-}\gg n_{+}$. In this case if MCR is minimized then the classifier will be biased toward the negative class due to its effort in minimizing the false positives (that is c) more strongly than minimizing false negatives (that is b). In the worst case almost all the positive examples will be wrongly classified, that is $a\rightarrow0$ . This results in both $P\rightarrow 0$ and $R\rightarrow 0$. Thus MCR is not a good measure to use when the dataset is imbalanced. This problem can be addressed by optimizing other measures that we discuss next. 

The F-measure is one such measure and is defined \citep{van:74} as: 
%\begin{equation}
	$$F_{\zeta}(P,R)\:=\:\Bigl({\zeta \over R}+{{1-\zeta}\over P}\Bigr)^{-1}$$
%\label{Fmeasure}
%\end{equation}
\begin{table}
%Keerthi1:
\caption{Confusion Matrix for Binary classification}
\vskip 0.05in
\begin{small}
\begin{center}
\begin{tabular}{|l|c|c|} \hline
 & Positive (Predicted) & Negative (Predicted)  \\ \hline
Positive (Actual) & a & b \\ \hline 
Negative (Actual)  & c & d \\ \hline 
\end{tabular}
\end{center}
\end{small}
\end{table}
where $0\le \zeta \le 1$ and $0\le F_{\zeta}(P,R)\le 1$. It has been used in various applications like document retrieval \citep{van:74} and text classification \citep{Joachims:05}. It is particularly preferable over MCR when the dataset is highly imbalanced \citep{Joachims:05,Jansche:05}. 

Let us get into more details on the functioning of the F-measure. When $\zeta \rightarrow 0$, we get $F_{\zeta}(P,R) \rightarrow P$. Then optimizing the F-measure means we are interested {\it only} in maximizing the {\it Precision}. On the other hand, when $\zeta \rightarrow 1$, we get $F_{\zeta}(P,R) \rightarrow R$. In this case we are interested {\it only} in maximizing the {\it Recall}. The user can choose an appropriate value for $\zeta$ depending on how much of importance he/she wants to give to the precision and recall. Thus, the F-measure combines precision and recall into a single optimization criterion by taking their $\zeta$-weighted harmonic mean. In the imbalanced data case mentioned above MCR minimization can potentially result in $F_{\zeta}(P,R)\rightarrow 0$. By maximizing the F-measure we can prevent the classifier from being completely biased towards the negative class. Note that $F_{\zeta}(P,R)$ can be re-written in terms of $a$, $b$ and $c$ as:
\begin{equation}
	F_{\zeta}(a,b,c)\:=\:{a \over {a+\zeta b+(1-\zeta) c}}
\label{Fabcmeasure}
\end{equation} 
%Keerthi: From 
%Keerthi1:
In all our experiments we set $\zeta\:=\:0.5$. In this case, it becomes $F_{0.5}(P,R)\:=\:{{2PR} \over {P+R}}$ and can also be written as: $F_{0.5}(a,b,c)\:=\:{1 \over {1+{{b+c}\over {2a}}}}$.  Then, maximizing $F_{0.5}(a,b,c)$ is equivalent to maximizing ${a} \over {b+c}$. Thus, we can maximize $F_{0.5}(a,b,c)$ by both minimizing the error (that is, $b+c$) and maximizing the true positives. The trade-off kicks-in since maximizing the true positives tends to increase the false positives. 
%Now, if $n_{-}\gg n_{+}$ then with $a\:=\:n_{+}$ and $b\:=\:0$ (that is, $100\%$ recall), it is quite possible that ${c \over {2n_{+}}} \gg 1$. This happens when the number of false positives is relatively large (that is, the precision drops down). In such a case, we can approximate $F_{0.5}(a,b,c)\:\approx\:{{2n_{+}} \over c}\ll1$ and the F-score drops down significantly. 
Thus maximizing the F-measure controls both the true positives and the error appropriately. In general the F-measure summarizes a classifier's ability to identify the positive class and plays an important role in the evaluation of binary classifier. As a criterion for optimizing hyperparameters 
%Keerthi: point of view, we proceed as follows. Note that 
F-measure can be computed 
%Keerthi: for 
on
an evaluation or validation dataset. 
%Keerthi1:
%However, it may not be possible to have an evaluation dataset in many practical situations. 
%Keerthi1:
However, in practical situations involving small datasets\footnote{GP models are known to be particularly valuable for problems with small datasets.} it is wasteful to employ a separate evaluation set.
%Keerthi: Therefore, 
The
LOO-CV approach would be useful in such situations and we show next how the F-measure can be estimated with such approach. 

\citet{Hong:07} estimated $TP$ and $FP$ as:
%Keerthi1: I suggest that we leave out the "-" superscripts since they don't add much info. ---------- Removed ---------------
%Keerthi: Why not use \bar{TP}, \bar{FP} etc instead of TP^-, FP^-?
	$$\widehat{TP}\:=\:{1 \over {n_+}} \sum_{i=1}^n T({\hat y}_{\setminus{i}}y_i,y_i),$$
	$$\widehat{FP}\:=\:{1 \over {n_-}} \sum_{i=1}^n F({\hat y}_{\setminus{i}}y_i,y_i).$$
Here ${\hat y}_{\setminus{i}}$ represents predicted label for $i$th sample. Therefore ${\hat y}_{\setminus{i}}y_i$ takes value $+1$ when the prediction matches with the actual label and $-1$ otherwise. $T(u,v)$ is an indicator function which is 1 if $u=1$ and $v=1$. Similarly, $F(u,v)$ is one if $u=-1$ and $v=-1$. Otherwise, these functions take zero values. 
%Keerthi: They 
\citet{Hong:07} 
used these estimates to compute $\widehat{AUC}\:=\:{{1+\widehat{TP}-\widehat{FP}} \over 2}$ as an approximation of AUC and used this criterion to select 
%Keerthi:
a
subset of basis vectors in 
%Keerthi:
a
kernel classifier model construction procedure for imbalanced datasets. Note that this definition of AUC is applicable only for a hard classifier (fixed non-probabilistic classifier) with binary outputs. See \citet{Hong:07} for more details. In a strict sense such a definition of AUC is not suitable for a probabilistic classifier like GP classifier that provides continuous probabilistic output. 
%Keerthi: Note that t
However, we can make use of this approach of defining $TP$ and $FP$ as above to compute the quantities $a$, $b$ and $c$ that are needed to evaluate the F-measure 
%Keerthi:
in (\ref{Fabcmeasure}). There are two issues associated with these estimates. The first issue is 
%Keerthi:
that
these estimates are not 
%Keerthi1:
smooth (in fact, not even continuous) 
%Keerthi:
functions
of hyperparameters. Therefore they cannot be used directly in any  
%Keerthi1:context; 
approach that uses gradient-based nonlinear optimization methods to tune the hyperparameters.
%Keerthi1:
Secondly, these estimates do not use predictive probability values which is particularly important when we want to take variance also into account. 
%Keerthi1:
In non-GP contexts
\citet{Jansche:05} and \citet{Keerthi:06} addressed the first issue by defining smoothed F-measure or 
%Keerthi:
other
validation functions by replacing the indicator function with a sigmoid  
%Keerthi:
function,
%Keerthi: . This approach 
which makes the optimization criterion as a 
%Keerthi: continuous 
smooth
function of hyperparameters. However, they did not consider 
%Keerthi:
a
LOO approach
%Keerthi:
and used a validation set instead.
%Keerthi:; further, 
\citet{Jansche:05} considered maximum {\it a posterior} probabilities and \citet{Keerthi:06} used sigmoidal approximations for SVM models. Here, we propose to combine LOO based estimation and 
%Keerthi: smoothening 
smoothed version of the quantities $\{a, b, c, d\}$ 
%Keerthi: to define smoothed LOO measures. Let 
%Keerthi:
denoted as $A(\bftheta)$, $B(\bftheta)$, $C(\bftheta)$ and $D(\bftheta)$.
%Keerthi: be written as follows: 
We can set
\begin{equation}
	A(\bftheta)\:=\: \sum_{i:y_i=+1}\:p(y_i=+1|x_i,{\it S}_{\setminus{i}},\bftheta)
\label{Atheta}
\end{equation}
Since $n_{+}=a+b$, we can write $B(\bftheta)\:=\:n_{+}-A(\bftheta)$. With $m_{+}$ denoting the number of examples predicted as positive, we can parameterize it as $m_{+}(\bftheta)\:=\:A(\bftheta)+C(\bftheta)$. This can be rewritten as: 
\begin{equation}
m_{+}(\bftheta)\:=\:\sum_{i=1}^n p(y_i=+1|x_i,{\it S}_{\setminus{i}},\bftheta)
\label{Mtheta}
\end{equation}
Thus, the smoothed F-measure can be defined from (\ref{Fabcmeasure}) as: 
\begin{equation}
	F_{\zeta}(\bftheta)\:=\:{{A(\bftheta)} \over {\zeta n_{+} + (1-\zeta) m_{+}(\bftheta)}}
\label{FMmooth}
\end{equation} 
Note that $D(\bftheta)$ can be defined in a similar fashion as $m_{-}(\bftheta)\:=\:B(\bftheta)+D(\bftheta)$. Using these quantities, other derived quantities like $TP(\bftheta)$ and $FP(\bftheta)$ can be defined as LOO based estimates. Then, smoothed LOO estimates of WER can be obtained as shown below.

The WER measure is another useful measure for imbalanced datasets. Using the quantities defined above, its smoothed version can be written as:
\begin{equation} 
WER(\bftheta;\tau)\:=\:{{n_{+}(1-TP(\bftheta))+\tau n_{-} FP(\bftheta)} \over {n_{+}+\tau n_{-}}}.
\label{WER}
\end{equation}
where $\tau$ is the ratio of the cost of mis-classifications of the negative class to that of the positive class and $0\le \tau \le 1$. Thus by choosing a suitable $\tau$ value for a given problem and optimizing over the hyperparameters we can design classifiers without becoming biased toward one class. Note that for ease of notation we have omitted hat on $TP(\cdot)$ and $FP(\cdot)$. 
%The AUC of ROC curve is another important measure used in the binary classification problem. The ROC curve is drawn in a two-dimensional plane with the x-axis representing the false positive rate(FP) and the y-axis representing the true positive rate(TP). Now, the classifiers that produce binary outputs (hard classifiers) can be represented as individual points (corresponding to pairs (FP,TP)) in the ROC plane. On the other hand, classifiers that produce continuous numeric outputs can be thresholded at different values. Depending on the classifier used this can be a probabilistic output. Each threshold value produces a point corresponding to a pair of (FP,TP) value and the ROC curve is obtained from varying the threshold. The ROC curve shows the trade-off between TP and FP. The trade-off happens since the TP cannot be increased without increasing the FP. The AUC of ROC curve is thus a threshold independent measure of a continuous output classifier's performance. This is equivalent to the probability that for a pair of randomly drawn samples from the positive and negative classes respectively, the classifier ranks the positive example higher than the negative example. The goal is then to maximize the AUC of the ROC curve. Note that such a curve is not available for a hard classifier and only a point (FP, TP) is available. For such a classifier the ROC can be defined (Hong et al., 2007) as follows: a curve is constructed by drawing a line from the point (FP,TP) to the point (0,0) on one side and to the point (1,1) on the other side in the ROC plane. Now the AUC of this ROC can be readily obtained as (Hong et al., 2007): $\hat{AUC}\:=\:{{1+\hat{TP}-\hat{FP}} \over 2}$. While this measure is strictly not the AUC for a continuous classifier, it has the desirable property of trading-off between high TP and low FP. Within this restriction this measure can be used for selecting the hyperparameters from its smoothed version given by:   
%\begin{equation}
%AUC(\bftheta)\:=\:{{1+TP(\bftheta)-FP(\bftheta)} \over 2}.   
%\label{AUC}
%\end{equation}
%On comparing (\ref{WER}) and (\ref{AUC}) we see that they are equivalent when $\tau=1$ in the sense that maximizing AUC is just a special case of minimizing WER. 
Following the work of \citet{Hong:07}  one can also define 
\begin{equation}
AUC(\bftheta)\:=\:{{1+TP(\bftheta)-FP(\bftheta)} \over 2}
\label{AUC}
\end{equation}
 and optimize over the hyperparameters. As mentioned earlier, such a definition is not suitable for the GP classifier. Nevertheless it is interesting to note that it has the desirable property of trading-off between high TP and low FP. Also, on comparing this definition of AUC with (\ref{WER}) we see that they are related in the sense that maximizing AUC is same as minimizing WER when $\tau=1$ {\it and} $n_{+}\:=\:n_{-}$. 

Overall we see that the LOO-CV predictive distributions can be used to define various criteria that are 
%Keerthi: continous 
smooth
functions of hyperparameters resulting in smoothed LOO-CV measures. Now given that the LOO-CV predictive distributions are readily available from the EP algorithm, we can optimize the various smoothed LOO-CV measures directly using standard non-linear optimization techniques.  

\section{EP-CV Algorithm for Choosing Hyperparameters}
%Keerthi: See my mail on reorganization.
Various criteria 
%Keerthi1: like 
such as 
%Keerthi1:
(\ref{nlp}), (\ref{FMmooth}), (\ref{WER}) and (\ref{AUC})  
%Keerthi1:
depend 
on the hyperparameters $\bftheta$ via the predictive distributions 
%Keerthi2:
$p(y_i|x_i,{\it S}_{\setminus{i}},\bftheta)$.
%Keerthi2:, $\forall i$. 
With cumulative Gaussian likelihood, they can be written as:
\begin{equation}
	p(y_i|x_i,{\it S}_{\setminus{i}},\bftheta)\:=\:\Phi\Bigl({{y_i(\mu_{\setminus{i}}+\gamma)} \over {\sqrt{1+\sigma^2_{\setminus{i}}}}}\Bigr)
\label{Pyi}
\end{equation}
Note that (\ref{Pyi}) is obtained from (\ref{pred1}) with $p(f_{i}|x_i,S_{\setminus{i}},\bftheta)=\mathcal{N}(\mu_{\setminus{i}},\sigma^2_{\setminus{i}})$. The hyperparameter $\gamma$ is referred to as the bias parameter and it helps in shifting the decision  
%Keerthi2: Note that t
%The bias parameter $\gamma_2$ helps in shifting the decision 
%Keerthi2: function 
boundary with the probability value ${1\over2}$. In general, the bias hyperparameter $\gamma$ is very useful \citep{Raswil:06, Seeger:03} and can be optimized. 

\subsection*{EP Approximation}
To compute (\ref{Pyi}) we need the LOO mean $\mu_{\setminus{i}}$ and variance 
%Keerthi1: Removed comma
$\sigma^2_{\setminus{i}}$ $\forall i$. 
With the EP approximation, they can be computed using (\ref{cavmu}). 
%Keerthi1:
%Also, we employ slightly different 
%Keerthi1:
%hyperparameter 
%optimization procedures for Laplace and EP approximations. 
%Keerthi1: We 
%Therefore we discuss the details for them 
%Keerthi1: separately in subsequent sections.    
%in separate subsections below.
%Keerthi1:We give details of gradient calculations in the appendix. Given these expressions, we can optimize over the hyperparameters using standard nonlinear optimization methods. 
Full details of gradient calculations needed for implementing hyperparameter optimization are given in the appendix.
%These quantities can be computed using the expressions that we provide in a later section and appendix. An important aspect of these computations is: we first compute the posterior $q({\bf f}|{\bf S},\bftheta)$; then, we compute the LOO predictive distributions and their derivatives within certain approximations using the mean and covariance information of the posterior. This approach makes LOO computations efficient. Note that the posterior from Laplace and EP approximations can be obtained as detailed in (Rasmussen and Williams, 2006). 
%\subsection{Laplace Approximation}
%The LOO mean and variance are computed efficiently from the posterior solution $q({\bf f}|{\bf S},\bftheta)$ by removing contributions due to the $i^{th}$ example. As detailed in the appendix, this is done in a simple way by relating latent function posteriors obtained from Laplace and EP approximations, and approximate cavity distributions. Using this approach, the approximate posterior predictive distribution of latent function $f_i$ is given by:
%%Keerthi1: Shouldn't we use $f_{\setminus{i}}$ below? ---- $f_{\setminus{i}}$ is strictly speaking a vector $f$ with $i^{th}$ example removed --- where as, $f_i$ implies we are predicting $i^{th}$ latent value using the dataset $S_{\setminus{i}}$ and input $x_i$. Is it not? ----------
%\begin{equation}
%	p(f_{i}|x_i,S_{\setminus{i}},\bftheta)\:=\:\mathcal{N}(\mu_{\setminus{i}},\sigma^2_{\setminus{i}}) \:\:\:\:\:\:\mu_{\setminus{i}}\:=\:{\hat f}_i\:-\:\alpha_i\sigma^2_{\setminus{i}}  \:\:\:\:\:\:\:\:\:\sigma^2_{\setminus{i}}\:=\:{{\tilde \sigma^2_i} \over {1-w_{ii}{\tilde \sigma^2_i}}}
%\label{lapmu}
%\end{equation}
%where ${\hat {\bf f}}\:=\:{\bf K}\bfalp$, $\alpha_i$ is $i^{th}$ element of $\bfalp\:=\:\nabla \log\:p({\bf y}|\hat {\bf f},\bftheta)$, ${\tilde \sigma^2_i}\:=\:k(i,i)-{\bf k}({\bf x}_i)^T({\bf K}\:+\:{\bf W}^{-1})^{-1}{\bf k}({\bf x}_i)$ and $w_{ii}$ are 
%%Keerthi1:
%the
%diagonal elements of the matrix ${\bf W}$. These expressions have a good intuitive interpretation. Note that if 
%%Keerthi1:
%the
%$i^{th}$ example is well explained in the sense that its likelihood probability (at the posterior solution) is in the flat region, then $\alpha_i\:\approx\:0$ and $\mu_{\setminus{i}}\:\approx\:{\hat f}_i$. Similarly, $w_{ii}\:\approx\:0$ and $\sigma^2_{\setminus_{i}}\:\approx\:{\tilde \sigma^2}_{i}$. This implies that removing 
%%Keerthi1:
%the
%$i^{th}$ sample does not affect its prediction without that example. Thus, the correction or rescaling happens in a sense with reference to the uncertainty associated in labeling that example at the posterior solution. 

%%Keerthi1: Finally, w
%We optimize the hyperparameters using the following LAP-CV algorithm. 
%\subsubsection*{{\bf LAP-CV Algorithm}:}
%\begin{enumerate}
%\item Initialize the hyperparameters $\bftheta.$
%\item Find the hyperparameters $\bftheta$ by optimizing over any LOO-CV based measure like (\ref{nlp}) or (\ref{FMmooth}) using any standard gradient based optimization technique. We carry out 
%%Keerthi1: $ls_{max}$ 
%a
%number of line searches in this optimization process and terminate the algorithm when there is no significant change in the objective function value. 
%%Keerthi1: During line search as the hyperparameters change, 
%For each new set of hyperparameters that we come across,
%we perform the following sequence of operations
%%Keerthi1:. 
%to determine the function and gradient.
%\begin{enumerate}
%%Keerthi1:
%%\item Compute the posterior mean ${\hat {\bf f}}$ and covariance $({\bf K}^{-1}+{\bf W})^{-1}$ using Newton's method.
%\item Compute the posterior mean ${\hat {\bf f}}$ using Newton's method and 
%%Keerthi1:
%then the covariance $({\bf K}^{-1}+{\bf W})^{-1}$.
%\item Compute the LOO mean $\mu_{\setminus{i}}$ and variance $\sigma^2_{\setminus{i}}$ using (\ref{lapmu}).
%\item Compute the chosen objective function like (\ref{nlp}) or (\ref{FMmooth}) and its 
%%Keerthi1: derivatives.
%gradient with respect to hyperparameters.
%\end{enumerate}   
%\end{enumerate}   

%\subsection*{EP Approximation}
%With the EP approximation, the LOO predictive distributions $\mathcal{N}(\mu_{\setminus{i}},\sigma^2_{\setminus{i}})$ are readily available where $\mu_{\setminus{i}}$ and $\sigma^2_{\setminus{i}}$ are given in (\ref{cavmu}). 
We take the expectation-maximization (EM) type approach for hyperparameters optimization.  
%Keerthi1: I feel this might invite some criticism. --- after the discussion we had, I am just leaving it as it is --- I hope that is fine ---
 This is because gradient expressions involving implicit derivatives (with site parameters varying as a function of hyperparameters) are not available due to the iterative nature of 
%Keerthi2:
the EP algorithm. This approach results in the following algorithm.             
%Keerthi2: I suggest we remove "EM" and simply call it as "LOO-CV-EP" everywhere. In fact, even better, change "LAP-CV" to "LAP-CV" and "EP-CV" to "EP-CV" to make them well in tune with the notations of section 6. I haven't made these changes since you may have some good objection.
\subsubsection*{\bf {EP-CV Algorithm}:}
\begin{enumerate}
\item Initialize the hyperparameters $\bftheta.$
\item {\bf Perform E-Step}: Given the hyperparameters, we find the site parameters  ${\bfmu}$ and $\Sigma$ and the posterior $q({\bf f}|{\it S},\bftheta)\:=\:\mathcal{N}({\bf m},{\bf C})$ using EP algorithm.
\item {\bf Perform M-Step}: Find the hyperparameters $\bftheta$ by optimizing over any LOO-CV based measure like (\ref{nlp}), (\ref{FMmooth}), (\ref{WER}) or (\ref{AUC}) using any standard gradient based optimization technique. We carry out just one line search in this optimization process. During this line search as the hyperparameters change, we perform the following sequence of operations. 
\begin{enumerate}
\item Compute the posterior mean ${\bf m}$ and covariance ${\bf C}$ using (\ref{posmean}).
\item Compute the LOO mean $\mu_{\setminus{i}}$ and variance $\sigma^2_{\setminus{i}}$ using (\ref{cavmu}).
\item Compute the chosen objective function like (\ref{nlp}), (\ref{FMmooth}), (\ref{WER}) or (\ref{AUC})  and its derivatives.
\end{enumerate}
Note that through out this M-step, it is assumed that the site parameters are fixed and the values obtained from step (2) are used.   
\item repeat steps 2-3 until there is no significant change in the objective function value.
\end{enumerate}   
This algorithm 
%Keerthi1:
worked
well in our experiments. 
%Keerthi1: Note that the 
A similar 
EM approach was used by \citet{Kim:06} (which they called EM-EP algorithm) in the optimization of a lower bound on the marginal likelihood. 

\begin{figure}
%\begin{center}
%\setlength{\epsfxsize}{3in}
%\centerline{\epsfbox{car3.eps}}
%\hspace{1cm}
%\centerline{\epsfbox{car3.eps}}
%\end{center}
%\centering
{
	\includegraphics[width=7cm]{car3.eps}
}
\hspace{1cm}
{
	\includegraphics[width=7cm]{yeast7.eps}

}
\caption{F-Measure Optimization on Car3 and Yeast7 Datasets}
\label{fig1}
\end{figure}
Since the EP-CV algorithm optimizes the smoothed F-measure it is useful to study the behavior of the true F-measure as optimization proceeds. We do this study on two of the datasets described in Table 6 of section 6. The optimization algorithm was terminated when there was no significant change in the smoothed F-measure value. From Figure 1 we see that the smoothed F-measure monotonically increases in both cases. Also as expected, the true F-measure exhibits a non-smooth behavior expected of a discrete function, and also, the values of true and smoothed F-measures are not the same. The difference arises because the smoothed F-measure is based on probabilistic scores, which can take any value between 0.5 and 1 depending on the problem (even when correct classification occurs). The important point to observe is that, in general, there is an increasing trend in true F-measure value as the optimization progresses. In the case of Car3 dataset (left panel) we see that clearly. A similar trend is seen in the case of Yeast7 dataset (right panel) also, except for a small dip at the 10th iteration. Though such a behavior happens sometimes in early iterations, we observed that better true F-measure value is almost always obtained as the optimization progresses.         
%\subsection{Approximate LOO Predictive Distributions}
%In the case of Laplace approximation, we compute the LOO mean and variance from the posterior solution obtained using all the examples by removing contributions due to the $i^{th}$ example. This can be done in a simple way by relating latent function posteriors obtained from Laplace and EP approximations, and approximate cavity distributions. Using this approach, the approximate posterior predictive distribution of latent function $f_{\setminus{i}}$ with Laplace approximation are given as $\mathcal{N}(\mu_{\setminus{i}},\sigma^2_{\setminus{i}})$ where the mean and variance are given by:
%\begin{equation}
%	\mu_{\setminus{i}}\:=\:{\hat f}_i\:-\:\alpha_i\sigma^2_{\setminus{i}}  \:\:\:\:\:\:\:\:\:\sigma^2_{\setminus{i}}\:=\:{{\tilde \sigma^2_i} \over {1-w_{ii}{\tilde \sigma^2_i}}}
%\label{lapmu}
%\end{equation}
%where ${\hat {\bf f}}\:=\:{\bf K}\bfalp$, $\alpha_i$ is $i^{th}$ element of $\bfalp\:=\:\nabla \log\:p({\bf y}|\hat {\bf f},\bftheta)$, ${\tilde \sigma^2_i}\:=\:k(i,i)-{\bf k}({\bf x}_i)^T({\bf K}\:+\:{\bf W}^{-1})^{-1}{\bf k}({\bf x}_i)$ and $w_{ii}$ are diagonal elements of the matrix ${\bf W}$. These expressions have a good intutitive interpretation. Note that if $i^{th}$ example is well explained in the sense that its likelihood probability (at the posterior solution) is in the flat region, then $\alpha_i\:\approx\:0$ and $\mu_{\setminus{i}}\:\approx\:{\hat f}_i$. Similarly, $w_{ii}\:\approx\:0$ and $\sigma^2_{\setminus_{i}}\:\approx\:{\tilde \sigma^2}_{i}$. This implies that removing $i^{th}$ sample does not affect its prediction without that example. Thus, the correction or rescaling happens in a sense with reference to the uncertainty associated in labeling that example at the posterior solution.   
%the derivatives with respect to the hyperparameters can be found as shown in the appendix. These expressions can be used with standard gradient based non-linear optimization techniques to optimize any of the LOO-CV based optimization criteria given below. 

%With EP approximation, the predictive distributions of $f_{\setminus{i}}$ are readily available as the cavity distributions $\mathcal{N}(\mu_{\setminus{i}},\sigma^2_{\setminus{i}})$ where $\mu_{\setminus{i}}$ and $\sigma^2_{\setminus{i}}$ are given in (\ref{cavmu}). 
\subsection*{Computational and Storage Complexities}
The computational complexity of the EP-CV algorithm depends on the 
%Keerthi1:
number of ARD kernel parameters $D$. See appendix for more details. For a given problem with $D$ fixed, the complexity is $O(n^3)$. This complexity is same as that of the EP(ML) method (see equation (\ref{EPML})) and LS-CV(NLP) method given in the next section. 
%In the case of ARD Gaussian kernel function as given in section 2, the complexity is $O(n^3D)$ where $D$ is the number of input features. This complexity is same 
%Keerthi1: as 
%for
%the LOO-CV based probabilistic LS classifier (Rasmussen and Williams, 2006). However, it is higher compared to 
%Keerthi1:
%the
%marginal likelihood approximation technique which has complexity $O(n^3)$. This difference in complexity comes from a matrix multiplication requirement (in 
%Keerthi1:
%all
%the LOO-CV based methods) that involves 
%Keerthi1: derivative matrix 
%a matrix of derivatives
%with respect to the hyperparameters; see appendix for more details. 
%Keerthi1:
%Note that a global scaling hyperparameter for all the input features is sufficient in many practical problems and in such cases, the complexities of LOO-CV based and marginal likelihood algorithms are same. 
%In many practical problems a single global scaling hyperparameter for all the input features is sufficient; for such cases, the complexities of LOO-CV based and marginal likelihood algorithms are same. 
%Keerthi1: Why O(n^3) below? --- It should be $O(n^2)$ --- corrected ---
Also in many practical problems a single global scaling hyperparameter for all the input features is sufficient. Finally, the storage complexities of all the methods are $O(n^2)$.       

%\section{LOO-CV Based Probabilistic Least Squares Classifier}
%\section{Discussion}
%Keerthi1: \section{Related Work: LOO-CV based GPC Model Selection}
\section{Other LOO-CV based Methods}
%Keerthi1:
Having discussed our approach in detail it is useful to recall and discuss other LOO-CV based GPC model selection methods in relation to it.

\citet{Opper:00} derived a mean field algorithm for binary classification with GPs based on the TAP approach originally proposed in 
%Keerthi2:
the
statistical physics of disordered systems. They showed that this approach yields an approximate LOO estimator for the generalization error. 
%Keerthi1:Note that t
This estimate is equivalent to the 
%Keerthi1: LOO CV 
LOO-CV
error estimate obtained from EP \citep{Minka:01}. 
%Keerthi1: We also note that i
Instead 
of optimizing over the hyperparameters, 
%Keerthi1: they 
\citet{Opper:00} used the LOO-CV error count (using indicator functions) to choose rough hyperparameter values by scanning a range of values. 

\citet{Qi:04} used the LOO-CV error estimate obtained from EP to determine ARD hyperparameters. They worked with the Gaussian process classifier where each input feature is associated with a weight parameter $v_i$ with the prior $\mathcal{N}(0,{\zeta_i}^{-1})$. The hyperparameters $\zeta_i^{-1}$
%Keerthi1: , $\forall i$ 
were obtained by 
%Keerthi1:
{\it maximizing the evidence}
using a fast sequential update based on the work of \citet{Faul:02}. The outcome of this optimization is that many $\zeta_i$s would go to infinity such that only a few nonzero weights $v_i$ will be present. 
%Keerthi2: Note that e
Even though the ARD hyperparameters were optimized by maximizing the evidence, 
%Keerthi2: due 
to prevent
overfitting 
%Keerthi2: they 
\citet{Qi:04} proposed to select the final model as the one that gives the minimum LOO-CV error count or probability. As the LOO-CV error count is discrete, they chose the model with maximum evidence when there is a tie in the count. Compared to this approach, we work with 
%Keerthi2:
the
GP classifier model (without the weight parameters) detailed in Section 2 and optimize over the hyperparameters (including ARD) directly with various LOO-CV based measures (including F-measure, WER etc.) using gradient information. 

In this context, 
%Keerthi2:
the
LOO-CV based probabilistic LS classifier \citep{Raswil:06} is a 
%Keerthi2:
more
direct LOO-CV based GPC model selection approach. 
%Keerthi2: We 
For the sake of completion we give some details here and 
%Keerthi2: we 
later
compare our algorithms with this approach 
%Keerthi: also 
in our experiments. This approach treats classification as a regression problem. 
%Keerthi: Some rewriting below can help. --- done ---
Note that the probabilistic interpretation of LS criterion implies 
%Keerthi2:
a
Gaussian noise model. But the output ${\bf y}$ can take only $+1$ or $-1$ which is 
%Keerthi2:
slightly
odd. However, this approach is simple to implement and a probabilistic interpretation is given by passing the predictions through a sigmoid. 
%Note that though the output ${\bf y}$ can take only $+1$ or $-1$, the probabilistic interpretation of LS criterion implies Gaussian noise model, which is odd. However, this approach is simple to implement and a probabilistic interpretation is given by passing the predictions through a sigmoid. 

Specifically, the LOO mean $\mu_{\setminus{i}}$ and variance $\sigma^2_{\setminus{i}},\:i\in {\tilde I}$ are obtained from LOO-CV formulation of GP regression and the predictive distributions are obtained 
%Keerthi2: as 
via
(\ref{Pyi}). Here, $\mu_{\setminus{i}}$ and $\sigma^2_{\setminus{i}}$ are given by:
\begin{equation}
	\mu_{\setminus{i}}\:=\:y_i\:-\:{\tilde \alpha}_i \sigma^2_{\setminus{i}}  \:\:\:\:\:\:\sigma^2_{\setminus{i}}\:=\:{1 \over  {{\bar {\bf K}}_{ii}}}
\label{lsmu}
\end{equation}
where ${\tilde \bfalp}\:=\:{\bar {\bf K}}{\bf y}$ and ${\bar {\bf K}}\:=\:({\bf K}+\lambda{\bf I})^{-1}$. 
%Keerthi2: Note that $\lambda\ge0$ and 
Here $\lambda$ can 
%Keerthi2:
either
be set to a small 
%Keerthi2:
positive
value or 
%Keerthi2: may be 
treated as a regularization hyperparameter with 
%Keerthi2: this small value 
a small upper bound 
constraint. This is useful when ${\bf K}$ can become ill-conditioned during optimization. Finally, the hyperparameters are optimized using 
%Keerthi2: the NLP measure 
(\ref{nlp}). We call this method as LS-CV(NLP).
%Keerthi2:
%We can also consider optimizing (12) instead using the same approach.

%\begin{table}
%\caption{Confusion Matrix: Binary classification case}
%\vskip 0.05in
%\begin{small}
%\begin{center}
%\begin{tabular}{|l|c|c|} \hline
% & Positive (Predicted) & Negative (Predicted)  \\ \hline
%Positive (Actual) & a & b \\ \hline 
%Negative (Actual)  & c & d \\ \hline 
%\end{tabular}
%\end{center}
%\end{small}
%\end{table}
   

%\subsubsection*{Algorithm}
%\begin{enumerate}
%\item Initialize the hyperparameters $\bftheta.$
%\end{enumerate}
\begin{table}
\caption{Various methods and their descriptions. All the EP-CV methods are optimized using EP-CV algorithm.}
\vskip 0.05in
\begin{small}
\begin{sc}
\begin{center}
\begin{tabular}{|l|l|} \hline
Method & Description  \\ \hline
EP(ML)  & Marginal likelihood maximization within EP approximation. \\
        & That is, optimize (\ref{EPML}) over $\bftheta$. \\ \hline 
EP-CV(NLP)  & Negative Logarithmic Predictive loss minimization within EP approximation. \\
            & That is, optimize (\ref{nlp}) over $\bftheta$. For ease of notation, this \\
            & method is referred as {\bf NLP} in Table 7. \\ \hline 
LS-CV(NLP)  & Negative Logarithmic Predictive loss minimization within Least Squares \\
            & approximation as described in Section 5. Optimize (\ref{nlp}) over $\bftheta$. \\ \hline 
EP-CV(FM)  & F-Measure maximization within EP-CV approximation. \\ 
           & That is, optimize (\ref{FMmooth}) over $\bftheta$. For ease of notation, this \\
           & method is referred as {\bf FM} in Table 7. \\ \hline 
NLP-FM(BIAS)     & In the first step, hyperparameters are optimized using EP-CV(NLP) method. \\
           & In the second step, only the bias parameter ($\gamma$) is optimized using \\ 
	   & EP-CV(FM) method.We also refer this method as two-step method. \\ \hline 
\end{tabular}
\end{center}
\vskip -0.2in
\end{sc}
\end{small}
\end{table}

\begin{table}
\caption{Data sets description: NLP Experiment. Here, $n$, $D$, $p$ and $nr$ represent the numbers of training examples, input dimension, test examples and train/test partitions respectively.}
\vskip 0.05in
\begin{small}
\begin{sc}
\begin{center}
\begin{tabular}{|l|l|l|l|l|} \hline
Dataset & $n$ & $D$  & $p$ & $nr$ \\ \hline
Banana  & 400 & 2 & 4900 & 100 \\ \hline 
Breastcancer  & 200 & 9 & 77 & 100 \\ \hline 
Diabetes  & 468 & 8 & 300 & 100  \\ \hline 
German & 700 & 20 & 300 & 100  \\ \hline
Heart  & 170 & 13 & 100 & 100  \\ \hline 
Image & 1300 & 18 & 1010 & 20  \\ \hline
Ringnorm & 400 & 20 & 7000 & 100  \\ \hline
Splice & 1000 & 60 & 2175 & 20 \\ \hline
Thyroid & 140 & 5 & 75 & 100 \\ \hline
Titanic & 150 & 3 & 2051 & 100 \\ \hline
Twonorm & 400 & 20 & 7000 & 100 \\ \hline
Waveform & 400 & 21 & 4600 & 100 \\ \hline
\end{tabular}
\end{center}
\vskip -0.2in
\end{sc}
\end{small}
\end{table}

%\begin{table}
%\caption{Data sets description: F-Measure Experiment}
%\vskip 0.05in
%\begin{small}
%\begin{sc}
%\begin{center}
%\begin{tabular}{|l|c|c|c|c|c|} \hline
%Dataset & n & D  & p & nr & PPE\\ \hline
%Yeast7  & 297 & 8 & 1187 & 25 & 2 \\ \hline 
%Yeast5  & 320 & 8 & 1164 & 25 & 4 \\ \hline 
%Car3  & 350 & 6 & 1378 & 25 & 4 \\ \hline 
%Ecoli5 & 124 & 7 & 212 & 25 & 6 \\ \hline
%Yeast4  & 165 & 8 & 1319 & 25 & 11 \\ \hline 
%Breastcancer & 200 & 9 & 77 & 100 & 29 \\ \hline
%German & 700 & 20 & 300 & 100 & 30 \\ \hline
%Diabetes & 468 & 8 & 300 & 100 & 35 \\ \hline
%\end{tabular}
%\end{center}
%\vskip -0.2in
%\end{sc}
%\end{small}
%\end{table}

\begin{table}
\begin{center}
\caption{NLP Performance} 
\vskip 0.1in
\begin{small}
\begin{sc}
\begin{tabular}{|l|c|c|c|} \hline
 Dataset/Method & EP(ML) & EP-CV(NLP)& LS-CV(NLP)\\ \hline
Banana & 23.90 $\pm$ 0.81 & 24.26 $\pm$ 1.06 & 33.88 $\pm$ 1.89 \\ \hline
Breastcancer & 53.57 $\pm$ 4.75 & 54.12 $\pm$ 5.27 & 55.58 $\pm$ 5.03 \\ \hline 
Diabetes & 47.74 $\pm$ 1.96 & 47.97 $\pm$ 2.09 & 50.72 $\pm$ 2.13 \\ \hline
German & 48.67 $\pm$ 2.74 & 49.05 $\pm$ 2.76 & 50.51 $\pm$ 2.30 \\ \hline 
Heart &  40.16 $\pm$ 5.36 & 40.03 $\pm$ 5.00 & 45.11 $\pm$ 4.91 \\ \hline
Image & 8.26 $\pm$ 1.07 & 8.45 $\pm$ 0.97 & 22.70 $\pm$ 0.66 \\ \hline
Ringnorm & 16.88 $\pm$0.93 & 16.56 $\pm$ 1.01 & 28.48 $\pm$ 0.75 \\ \hline 
Solar & 57.25 $\pm$ 1.38 & 57.35 $\pm$ 1.42 & 59.61 $\pm$ 1.34 \\ \hline
Splice & 28.48 $\pm$ 0.88 & 29.60 $\pm$ 0.79 & 36.83 $\pm$ 0.42 \\ \hline
Thyroid & 10.21 $\pm$ 3.76 & 9.94 $\pm$ 3.69 & 25.33 $\pm$ 4.86 \\ \hline
Titanic & 66.86 $\pm$ 1.97 & 51.73 $\pm$ 1.73 & 53.78 $\pm$ 14.08 \\ \hline
Twonorm & 8.31 $\pm$ 0.88  & 9.08 $\pm$ 1.97 & 25.94 $\pm$ 0.53 \\ \hline
Waveform & 23.01 $\pm$ 0.89 & 22.97 $\pm$ 0.67 & 32.63 $\pm$ 0.59 \\ \hline
\end{tabular}
\end{sc}
\end{small}
\end{center}
\end{table}

\begin{table}
\begin{center}
\caption{Test Set Error Performance} 
\vskip 0.1in
\begin{small}
\begin{sc}
\begin{tabular}{|l|c|c|c|} \hline
 Dataset/Method & EP(ML) & EP-CV(NLP)& LS-CV(NLP)\\ \hline
Banana & 10.41 $\pm$ 0.65 & 10.51 $\pm$ 0.50 & 10.93 $\pm$ 0.67 \\ \hline
Breastcancer & 26.52 $\pm$ 4.89 & 26.61 $\pm$ 4.80 & 25.94 $\pm$ 4.59 \\ \hline 
Diabetes & 23.28 $\pm$ 1.82 & 23.41 $\pm$ 1.82 & 24.30 $\pm$ 2.51 \\ \hline
German & 23.36 $\pm$ 2.11 & 23.48 $\pm$ 2.00 & 23.94 $\pm$ 2.33 \\ \hline 
Heart & 16.65 $\pm$ 2.87 & 16.62 $\pm$ 3.08 & 17.91 $\pm$ 4.21 \\ \hline  
Image & 2.82 $\pm$ 0.54  & 2.77 $\pm$ 0.51 & 2.74 $\pm$ 0.65 \\ \hline
Ringnorm & 4.41 $\pm$ 0.64 & 4.29 $\pm$ 0.69 & 5.05 $\pm$ 0.99 \\ \hline
Solar & 34.20 $\pm$ 1.75 & 34.27 $\pm$ 1.80 & 35.03 $\pm$ 1.89 \\ \hline 
Splice & 11.61 $\pm$ 0.81 & 11.85 $\pm$ 0.83 & 11.83 $\pm$ 0.80 \\ \hline 
Thyroid & 4.37 $\pm$ 2.19 & 4.20 $\pm$ 2.17 & 6.97 $\pm$ 3.78 \\ \hline
Titanic & 22.64 $\pm$ 1.34 & 22.50 $\pm$ 0.98 & 22.99 $\pm$ 2.81 \\ \hline
Twonorm & 3.05 $\pm$ 0.34 & 3.19 $\pm$ 0.51 & 3.43 $\pm$ 0.43 \\ \hline 
Waveform & 10.10 $\pm$ 0.48 & 9.95 $\pm$ 0.48 &  11.70 $\pm$ 0.88 \\ \hline
\end{tabular}
\end{sc}
\end{small}
\end{center}
\end{table}

\begin{table}
\caption{Data sets description: F-Measure Experiment. Here, $n$, $D$, $p$, $nr$ and $PPE$ represent the numbers of training examples, input dimension, test examples, train/test partitions and approximate percentage of positive examples respectively.}
\vskip 0.05in
\begin{small}
\begin{sc}
\begin{center}
\begin{tabular}{|l|c|c|c|c|c|} \hline
Dataset & $n$ & $D$  & $p$ & $nr$ & $PPE$\\ \hline
Yeast7  & 297 & 8 & 1187 & 50 & 2 \\ \hline 
Yeast5  & 320 & 8 & 1164 & 50 & 4 \\ \hline 
Car3  & 350 & 6 & 1378 & 50 & 4 \\ \hline 
Ecoli5 & 124 & 7 & 212 & 50 & 6 \\ \hline
Yeast4  & 165 & 8 & 1319 & 50 & 11 \\ \hline 
Breastcancer & 200 & 9 & 77 & 100 & 29 \\ \hline
German & 700 & 20 & 300 & 100 & 30 \\ \hline
Diabetes & 468 & 8 & 300 & 100 & 35 \\ \hline
\end{tabular}
\end{center}
\vskip -0.2in
\end{sc}
\end{small}
\end{table}

\begin{table}
\begin{center}
\caption{F-Measure Performance} 
\vskip 0.1in
\begin{small}
\begin{sc}
\begin{tabular}{|l|c|c|c|c|} \hline
 Dataset/Method & NLP & FM  &  NLP-FM(bias) \\ \hline
yeast7 & 32.24 $\pm$ 15.54 & 42.58 $\pm$ 7.84 &  40.85 $\pm$ 8.59 \\ \hline
yeast5 & 19.79 $\pm$ 12.85 & 32.85 $\pm$ 8.56 &  28.45 $\pm$ 10.57 \\ \hline
car3 & 55.89 $\pm$ 9.30 & 64.35 $\pm$ 8.51 & 62.67 $\pm$ 8.49 \\ \hline
ecoli5 & 84.41 $\pm$ 6.25 & 83.79 $\pm$ 5.90 &  84.08 $\pm$ 5.86 \\ \hline
yeast4 & 71.35 $\pm$ 3.95 & 73.64 $\pm$ 2.97 &  73.23 $\pm$ 2.81 \\ \hline
BreastCancer & 38.42 $\pm$ 10.55 & 47.34 $\pm$ 6.44  & 46.98 $\pm$ 6.37 \\ \hline
German & 54.15 $\pm$ 4.17 & 57.53 $\pm$ 2.95 &  56.91 $\pm$ 2.87 \\ \hline
Diabetes & 62.69 $\pm$ 3.49 & 66.23 $\pm$ 2.87 & 66.12 $\pm$ 2.67 \\ \hline
\end{tabular}
\end{sc}
\end{small}
\end{center}
\end{table}

\section{Experiments}
We conducted two experiments with various methods. See Table 2 for a summary of the various methods. In the first experiment we compared the performance of EP-CV(NLP) method with that of EP(ML) and LS-CV(NLP) methods. In the second experiment we compared the performance of EP-CV(FM) method with that of EP-CV(NLP) and two step classifier methods. We used the {\it minimize} Matlab routine of the GPML Matlab code available at {\it http://www.gaussianprocess.org/gpml/code/matlab/doc/} for hyperparameters optimization. In all the experiments we used a single global scaling hyperparameter.

\subsection{NLP Experiment}
In this experiment we used the thirteen benchmark datasets available in the web{\footnote {\it http://ida.first.fraunhofer.de/projects/bench/benchmarks.htm}} summarized in Table 3. Let us first consider the results from the first experiment given in Table 4 and Table 5. For the EP(ML) method we used the GPML Matlab code available in the web{\footnote {\it http://www.gaussianprocess.org/gpml/code/matlab/doc/}}.   

We conducted Friedman test \citep{Demsar:06} with the corresponding post-hoc tests for comparison of classifiers over multiple datasets. The comparison over multiple datasets requires a performance score of each method on each dataset \citep{Demsar:06}. Here, we consider the mean over the partitions of a given dataset as the performance score. As pointed out in \citet{Demsar:06}, it is not clear how to make use of the standard deviation information when the datasets are not independent over the partitions. The Friedman test ranks the methods for each dataset separately based on the chosen performance score (mean performance in our case). The best performing method gets the rank of 1, the second best rank 2 and so on. In case of ties, average ranks are assigned. The Friedman test checks whether the measured average ranks (over the datasets) are significantly different from the mean rank under the null hypothesis. Under the null hypothesis all the methods are equivalent and so their ranks should be equal. 

In the case of NLP performance measure, the measured average ranks for the EP(ML), EP-CV(NLP) and LS-CV(NLP) methods were 
$1.46$, $1.62$ and $2.92$ respectively. With three methods and 13 datasets, the F-statistic comparison at a significance level of 0.05 rejected the null hypothesis. Since the null hypothesis was rejected we conducted the Nemenyi post-hoc test for pairwise comparisons. This test revealed that the results of the EP(ML) and EP-CV(NLP) methods are better than the LS-CV(NLP) method at the significance level of $0.05$. On the other hand, the post-hoc test did not detect any significant difference in the results of EP(ML) and EP-CV(NLP) methods. In the case of test set performance measure, the measure averaged ranks for the EP(ML), EP-CV(NLP) and LS-CV(NLP) methods were $1.62$, $1.77$ and $2.62$ respectively. Note that the average rank of LS-CV method has improved on the test set error performance. Here again, the null hypothesis is rejected at the same significance level and the post-hoc test did not detect any significant difference in the results of EP(ML) and EP-CV(NLP) methods. The results of the EP(ML) and EP-CV(NLP) methods are better than the LS-CV(NLP) method at the significance level of $0.05$ and $0.1$ respectively. Thus, we can conclude that EP(ML) and EP-CV(NLP) are competitive to each other. Further, both these methods perform better than the LS-CV(NLP) method in this experiment.  

We can also make other observations from the tables. We note that the NLP performance of the LS-CV(NLP) method is quite inferior on several datasets even though its test set error performances on most of these datasets (except {\it waveform and thyroid}) are relatively closer. Further, some kind of group behavior can be seen. For example, the NLP scores are high on {\it titanic}, {\it breast-cancer}, {\it diabetes}, {\it German} and {\it flare-solar}. Also, the test set errors are $>20\%$ on these datasets. Consequently, we may consider these datasets as difficult ones. From Table 4 we observe that the NLP performance of LS-CV(NLP) is {\it closer} to the other two methods on these datasets (compared to its performance on other datasets). Next we can order the remaining datasets {\it heart}, {\it splice}, {\it banana}, {\it waveform}, {\it ringnorm}, {\it thyroid}, {\it twonorm} and {\it image} in terms of descending difficulty. Note that the difference in the NLP performance seems to have an increasing trend as the dataset becomes easier. To understand this we looked at the predictive probabilities of these methods on both correctly and wrongly classified examples. In the case of thyroid dataset these average probability scores for the LS-CV(NLP) method were $0.84$ (for correct classification) and $0.35$ (wrong classification). Here, the averaging was done over all the correctly(wrongly) classified examples over all the partitions. On the other hand, the corresponding values for the EP-CV(NLP) were $(0.96, 0.39)$. In the case of banana dataset, these scores were $(0.79, 0.35)$ and $(0.92, 0.29)$ for the LS-CV(NLP) and EP-CV(NLP) methods respectively. In the case of German dataset, they were $(0.75, 0.33)$ and $(0.79, 0.32)$. The scores for the EP-ML method were very close to that of the EP-CV(NLP) method. In general, we observed that the predictive probability estimates from the LS-CV(NLP) method were relatively poor and resulted in poor NLP performance. We looked at the hyperparameter estimates of the different methods and observed that the LS-CV(NLP) method takes smaller width and signal variance hyperparameter values on most of the datasets except on the difficult datasets (mentioned above) compared to the other two methods. Apart from this we did not observe any specific pattern in the hyperparameter values chosen by these methods. Looking at the hyperparameter estimates of EP(ML) and EP-CV(NLP), it seems that several solutions in the space of hyperparameters that give close performances are possible.       

\subsection{F-measure Experiment}
 The datasets used in this experiment are described in Table 6.  The datasets {\it yeast}, {\it car} and {\it ecoli} are multi-class datasets and we converted them into binary classification datasets by considering examples belonging to the class label indicated by the number (for example, $7$ in yeast7) as positive class respectively and treating the rest of the examples as negative class. These datasets are available in the web  {\footnote {\it ftp://ftp.ics.uci.edu/pub/machine-learning-databases/}}. We created 50 partitions for these datasets in a stratified manner reflecting the class distributions.

Let us consider the results from the second experiment given in Table 7. In this experiment all the results were obtained within the EP-CV framework. The first and second columns represent results obtained using NLP and smoothed F-measure (i.e., eq.~(\ref{FMmooth})) as the optimization criterion respectively. The third column represents results obtained from the two step classifier described earlier. We looked at the hyperparameter estimates of the different methods. We observed that the bias estimates of the two step classifier were somewhat closer (within $10\%$) to those of the smoothed F-measure method on the {\it breast-cancer}, {\it diabetes} and {\it German} datasets. On the remaining datasets they were different by more than $40\%$. The width and signal variance hyperparameter estimates were also quite different. From Table 7, we observe that the two step method is also good and gives closer performance to the smoothed F-measure method on several datasets. Further analysis of the performance results revealed that even though the standard deviations are high, the smoothed F-measure method gave better performance than the two step method on majority of the partitions on several datasets. We believe that the larger standard deviations in the results arises from the sensitivity to the dataset with lesser number of positive examples. To carry out the statistical significance tests, again we used the mean (over the partitions) as the performance score for each of the methods. The measured average ranks for these three methods were $2.75$, $1.25$ and $2.00$ respectively. With three methods and 8 datasets, the F-statistic comparison at a significance level of 0.05 rejected the null hypothesis. The Nemenyi post-hoc test for pairwise comparisons revealed that {\it only} the results of smoothed F-measure is better than the NLP based method at the significance level of 0.05. In this experiment the post-hoc test did not detect any significant differences in the comparisons of smoothed F-measure method with the two step method and the two step method with the NLP method. However in these two comparisons the rank differences were closer to the required critical differences at the significance level of 0.1.  We also observed that if we were to conduct Wilcoxon signed-rank test on these methods (as if we were comparing only two classifiers) then the results were statistically significant at the significance level of 0.05 for all the three pairs. In summary, the results demonstrate the usefulness of direct optimization of smoothed F-measure.
  
\section{Conclusion}
In this paper, we considered the problem of Gaussian process classifier model selection with different LOO-CV based optimization criteria and provided a practical algorithm using LOO predictive distributions with criteria like standard NLP, smoothed F-measure and WER to select hyperparameters. More specifically, apart from optimization of standard NLP, we demonstrated its usefulness in direct optimization of smoothed F-measure, which is useful to handle imbalanced data. We considered predictive distribution arrived from the Expectation Propagation (EP) approximation. We derived relevant expressions and proposed a very useful EP-CV algorithm. The experimental results on several real world benchmark datasets showed comparable NLP generalization performance (with NLP optimization) with existing approaches. We demonstrated that the smoothed F-measure optimization method is a very useful method that improves the F-measure performance significantly. Overall, the EP-CV algorithm is an excellent choice for GP classifier model selection with different LOO-CV based optimization criteria.         

%\section*{\noindent References}
%\begingroup
%\parindent=0pt
%\normalsize
%Faul, A. C., and Tipping, M. E. (2002). Analysis of sparse Bayesian learning. {\em Advances in Neural Information Processing Systems 14}.

%Geisser, S. (1975). The predictive sample reuse methods with applications.  
%{\em Journal of American Statistical Association, 70} (35), 320-328.

%Gibbs, M., and MacKay, D. J. C. (2000). Variational Gaussian process classifiers. {\em IEEE Transactions on Neural Networks}, Vol. 11, no. 6, p. 1458. 

%Hong, X., Chen, S., and Harris, C.J (2007) A kernel based two-class classifier for imbalanced data sets. {\em IEEE Transactions on Neural Networks}, Vol. 18, No.1, 28-41. 

%Keerthi, S. S., Sindhwani, V., and Chapelle, O.  (2006). An efficient method for gradient based adaptation of hyperparameters in SVM models. {\em Advances in Neural Information Processing Systems 18}. 
 
%Kim, H. C., and Ghahramani, Z. (2006). Bayesian Gaussian process classification with the EM-EP algorithm. {\em IEEE Transactions on Pattern analysis and machine intelligence}, Vol. 28, No. 12, 1948-1959.  

%Kuss, M., and Rasmussen, C. E. (2005). Assessing approximate inference for Bayesian Gaussian process classification. {\em Journal of Machine Learning Research}, Vol. 6, 1679-1704. 

%Lim, T.-S., Loh, W.-Y., and Shih, Y.-S. (2000). A comparison of prediction accuracy, complexity, and training time of thirty-three old and new classification algorithms. {\it Machine Learning, 40}, 203-229.

%Jansche, M. (2005). Maximum expected F-measure training of logistic regression models. Proceedings of the conference on Human Language Technology and Empirical Methods in Natural Language Processing, pp. 692-699, Association for Computational Linguistics, NJ, USA. 

%Minka, T. P. (2001). Expectation propagation for approximate Bayesian inference. {\em Uncertainty in AI}.

%Neal, R. (1997). Regression and classification using Gaussian process priors. {\em Bayesian Statistics 6}, 475-501.

%Opper, M., and Winther, O. (2000). Gaussian processes for classification: Mean field algorithms. {\em Neural Computation}. 


%Qi, Y., Minka, T. P., Picard, R. W., and Ghahramani, Z. (2004). Predictive automatic relevance determination by expectation propagation. {\em Proceedings of 21st International Conference on Machine Learning.}

%Rasmussen, C. E., and Williams, C. K. I. (2006). {\em Gaussian processes for machine learning}. The MIT press. 

%Seeger, M. (2003). Bayesian Gaussian process models: PAC-Bayesian generalisation error bounds and sparse approximations. {\em PhD thesis, University of Edinburgh.} 

%Seeger, M. (2006). Cross-validation optimization for large scale hierarchical classification kernel methods. {\em Advances in Neural Information Processing Systems 19}.

%Stone, M. (1974). Cross-validatory choice and assessment of statistical predictions (with discussion). {\em Journal of Royal Statistical Society (series-B), 36}, 111-147.

%Sundararajan, S., and Keerthi, S. S. (2001). Predictive approaches for choosing hyperparameters in Gaussian Processes. {\em Neural Computation, 13}(5), 1103-1118. 

%van Rijsbergen, C. J. (1974). Foundation of evaluation. {\em Journal of Documentation}, {\bf 30(4)}: 365-373.

%Wahba, G. (1998). Support Vector machines, reproducing kernel Hilbert spaces and the randomized GACV. {\em Technical report no: 984, Dept. of statistics, University of Wisconsin}.

%Williams, C.K.I. and Barber, D. (1998). Bayesian classification with Gaussian processes. {\em IEEE Transactions on Pattern analysis and machine intelligence}, {\bf 20}, no. 12, 1342-1351. 

% Acknowledgements should go at the end, before appendices and references

%\acks{We would like to acknowledge support for this project
%from the National Science Foundation (NSF grant IIS-9988642)
%and the Multidisciplinary Research Program of the Department
%of Defense (MURI N00014-00-1-0637). }

% Manual newpage inserted to improve layout of sample file - not
% needed in general before appendices/bibliography.

%\newpage

%\appendix
\vskip 0.2in
\bibliography{cvgpc}
\input{cvgpcappendixv5.tex}
\end{document}


cvgpcv16.bbl
\begin{thebibliography}{19}
\providecommand{\natexlab}[1]{#1}
\providecommand{\url}[1]{\texttt{#1}}
\expandafter\ifx\csname urlstyle\endcsname\relax
  \providecommand{\doi}[1]{doi: #1}\else
  \providecommand{\doi}{doi: \begingroup \urlstyle{rm}\Url}\fi

\bibitem[Demsar(2006)]{Demsar:06}
J.~Demsar.
\newblock Statistical comparisons of classifiers over multiple data sets.
\newblock \emph{Journal of Machine Learning Research}, 7:\penalty0 1--30, 2006.

\bibitem[Faul and Tipping(2002)]{Faul:02}
A.~C. Faul and M.~E. Tipping.
\newblock Analysis of sparse {B}ayesian learning.
\newblock \emph{Advances In Neural Information Processing Systems}, 14, 2002.

\bibitem[Gibbs and MacKay(2000)]{Gibbs:00}
M.~Gibbs and D.~J.~C. MacKay.
\newblock Variational {G}aussian process classifiers.
\newblock \emph{IEEE Transactions on Neural Networks}, 11\penalty0
  (6):\penalty0 1458--1464, 2000.

\bibitem[Hong et~al.(2007)Hong, Chen, and Harris]{Hong:07}
X.~Hong, S.~Chen, and C.~J. Harris.
\newblock A kernel based two-class classifier for imbalanced data sets.
\newblock \emph{IEEE Transactions on Neural Networks}, 18\penalty0
  (1):\penalty0 28--41, 2007.

\bibitem[Jansche(2005)]{Jansche:05}
M.~Jansche.
\newblock Maximum expected f-measure training of logistic regression models.
\newblock In \emph{Proceedings of the Conference on Human Language Technology
  and Empirical Methods in Natural Language Processing}, pages 692--699.
  Association for Computational Linguistics, NJ, USA, 2005.

\bibitem[Joachims(2005)]{Joachims:05}
T.~Joachims.
\newblock A support vector method for multivariate performance measures.
\newblock In \emph{ICML}, pages 377--384, 2005.

\bibitem[Keerthi et~al.(2006)Keerthi, Sindhwani, and Chapelle]{Keerthi:06}
S.~S. Keerthi, V.~Sindhwani, and O.~Chapelle.
\newblock An efficient method for gradient based adaptation of hyperparameters
  in {SVM} models.
\newblock \emph{Advances In Neural Information Processing Systems}, 18, 2006.

\bibitem[Kim and Ghahramani(2006)]{Kim:06}
H.~C. Kim and Z.~Ghahramani.
\newblock {B}ayesian {G}aussian process classification with the {EM-EP}
  algorithm.
\newblock \emph{IEEE Transactions on Pattern Analysis and Maching
  Intelligence}, 28\penalty0 (12):\penalty0 1948--1959, 2006.

\bibitem[Kuss and Rasmussen(2005)]{Kuss:05}
M.~Kuss and C.~E. Rasmussen.
\newblock Assessing approximate inference for {B}ayesian {G}aussian process
  classification.
\newblock \emph{Journal of Machine Learning Research}, 6:\penalty0 1679--1704,
  2005.

\bibitem[Minka(2001)]{Minka:01}
T.~P. Minka.
\newblock Expectation propagation for approximate {B}ayesian inference.
\newblock In \emph{UAI}, pages 362--369. Morgan Kaufmann, 2001.

\bibitem[Neal(1998)]{Neal:98}
R.~Neal.
\newblock Regression and classification using {G}aussian process priors.
\newblock \emph{{B}ayesian Statistics}, 6:\penalty0 475--501, 1998.

\bibitem[Opper and Winther(2000)]{Opper:00}
M.~Opper and O.~Winther.
\newblock {G}aussian processes for classification: Mean-field algorithms.
\newblock \emph{Neural Computation}, 12:\penalty0 2655--2684, 2000.

\bibitem[Qi et~al.(2004)Qi, Minka, Picard, and Ghahramani]{Qi:04}
Y.~Qi, T.~P. Minka, R.~W. Picard, and Z.~Ghahramani.
\newblock Predictive automatic relevance determination by expectation
  propagation.
\newblock In \emph{ICML}, pages 362--369, 2004.

\bibitem[Rasmussen and Williams(2006)]{Raswil:06}
C.~E. Rasmussen and C.~Williams.
\newblock \emph{{G}aussian {P}rocesses for {M}achine {L}earning}.
\newblock The MIT Press, 2006.

\bibitem[Seeger(2003)]{Seeger:03}
M.~Seeger.
\newblock {B}ayesian {G}aussian process models: Pac-{B}ayesian generalization
  error bounds and sparse approximation.
\newblock In \emph{Phd Thesis}. University of Edinburgh, 2003.

\bibitem[Seeger(2007)]{Seeger:07}
M.~Seeger.
\newblock Cross-validation optimization for large scale hierarchical
  classification kernel methods.
\newblock \emph{Advances In Neural Information Processing Systems}, 19, 2007.

\bibitem[Sundararajan and Keerthi(2001)]{Sundar:01}
S.~Sundararajan and S.~S. Keerthi.
\newblock Predictive approaches for choosing hyperparameters in {G}aussian
  processes.
\newblock \emph{Neural Computation}, 13\penalty0 (5):\penalty0 1103--1118,
  2001.

\bibitem[van Rijsbergen(1974)]{van:74}
van Rijsbergen.
\newblock Foundation of evaluation.
\newblock \emph{Journal of Documentation}, 30\penalty0 (4):\penalty0 365--373,
  1974.

\bibitem[Williams and Barber(1998)]{Williams:98}
C.~Williams and D.~Barber.
\newblock {B}ayesian classification with {G}aussian processes.
\newblock \emph{IEEE Transactions on Pattern Analysis and Machine
  Intelligence}, 20\penalty0 (12):\penalty0 1342--1351, 1998.

\end{thebibliography}
