\documentclass[11pt]{article} \usepackage{packages/mlapa} \usepackage{packages/fullpage} \usepackage{graphics} \usepackage{graphicx} \usepackage{psfrag} \usepackage{amsmath} \usepackage{amsthm} \usepackage{amssymb} \usepackage{psfig} \usepackage{amsfonts} \usepackage{latexsym} \usepackage{amscd} \renewcommand{\baselinestretch}{1.1} \newcommand{\ignore}[1]{} \newcommand{\deff}{\mbox{$\;\stackrel{\mbox{\tiny\rm def}}{=}\;$}} \newcommand{\bb}{\boldsymbol} \newcommand{\dd}{\mbox} \def\norm#1{\left\|#1\right\|} \newtheorem{prop}{Proposition} \newtheorem{lemma}{Lemma} \newtheorem{thrm}{Theorem} \newtheorem{defn}{Definition} \newtheorem{corr}{Corollary} \def\cA{\mathcal{A}} % caligraphic math letters \def\cX{\mathcal{X}} \def\cY{\mathcal{Y}} \def\cM{\mathcal{M}} \def\cN{\mathcal{N}} \def\cP{\mathcal{P}} \def\cF{\mathcal{F}} \def\cE{\mathcal{E}} \def\cQ{\mathcal{Q}} \def\cD{\mathcal{D}} \def\cL{\mathcal{L}} \def\Rn{\mathbb{R}^n} % blackboard math letters \def\R{\mathbb{R}} \def\Pn{\mathbb{P}_n} \def\Pm{\mathbb{P}_m} \def\P{\mathbb{P}} \def\H{\mathbb{H}} \def\N{\mathbb{N}} \def\cJ{\mathcal{J}} \def\Rn1{\mathbb{R}^{n+1}} \def\Sn{\mathbb{S}^{n}} \def\Snp{\mathbb{S}^{n}_{+}} \def\ss#1{\mbox{\it\scriptsize #1}} % font sizes \def\st#1{\mbox{\it\tiny #1}} \def\argmin{\operatornamewithlimits{arg\, min}} % new operators \def\argmax{\operatornamewithlimits{arg\, max}} \def\dotp#1#2{\left\langle{#1}, {#2}\right\rangle} \def\dvol{\text{dvol}\,} \def\seq{\;=\;} \def\given{\,|\,} \def\det{\text{det}\,} \def\div{\text{div}\,} \def\grad{\text{grad}\,} %----------------------------------------------------------------- %----------------------------------------------------------------- %----------------------------------------------------------------- \title{Riemannian Geometry and Statistical Machine Learning \\ \vspace{0.5in} \normalsize{\textsc{Doctoral Thesis\\ \vspace{0.3in} }}} \vspace{1in} \author{Guy Lebanon\\ \small{Language Technologies Institute}\\ \small{School of Computer Science}\\ \small{Carnegie Mellon University}\\ \small{\texttt{lebanon@cs.cmu.edu}}\\} \date{\today} \begin{document} \maketitle \vspace{1in} \begin{abstract} Statistical machine learning algorithms deal with the problem of selecting an appropriate statistical model from a model space $\Theta$ based on a training set $\{x_i\}_{i=1}^N \subset \cX$ or $\{(x_i,y_i)\}_{i=1}^N \subset \cX\times \mathcal{Y}$. In doing so they either implicitly or explicitly make assumptions on the geometries of the model space $\Theta$ and the data space $\cX$. Such assumptions are crucial to the success of the algorithms as different geometries are appropriate for different models and data spaces. By studying these assumptions we are able to develop new theoretical results that enhance our understanding of several popular learning algorithms. Furthermore, using geometrical reasoning we are able to adapt existing algorithms such as radial basis kernels and linear margin classifiers to non-Euclidean geometries. Such adaptation is shown to be useful when the data space does not exhibit Euclidean geometry. In particular, we focus in our experiments on the space of text documents that is naturally associated with the Fisher information metric on corresponding multinomial models. \end{abstract} \vspace{1in} \begin{tabular}{ll} Thesis Committee:& John Lafferty (chair)\\ &Geoffrey J. Gordon, Michael I. Jordan, Larry Wasserman \end{tabular} \newpage \section*{Acknowledgements} This thesis contains work that I did during the years 2001-2004 at Carnegie Mellon University. During that period I received a lot of help from faculty, students and friends. However, I feel that I should start by thanking my family: Alex, Anat and Eran Lebanon, for their support during my time at Technion and Carnegie Mellon. Similarly, I thank Alfred Bruckstein, Ran El-Yaniv, Michael Lindenbaum and Hava Siegelmann from the Technion for helping me getting started with research in computer science. \paragraph{} At Carnegie Mellon University I received help from a number of people, most importantly my advisor John Lafferty. John helped me in many ways. He provided excellent technical hands-on assistance, as well as help on high-level and strategic issues. Working with John was a very pleasant and educational experience. It fundamentally changed the way I do research and turned me into a better researcher. I also thank John for providing an excellent environment for research without distractions and for making himself available whenever I needed. \paragraph{} I thank my thesis committee members Geoffrey J. Gordon, Michael I. Jordan and Larry Wasserman for their helpful comments. I benefited from interactions with several graduate students at Carnegie Mellon University. Risi Kondor, Leonid Kontorovich, Luo Si, Jian Zhang and Xioajin Zhu provided helpful comments on different parts of the thesis. Despite not helping directly on topics related to the thesis, I benefited from interactions with Chris Meek at Microsoft Research, Yoram Singer at the Hebrew University and Joe Verducci and Douglas Critchlow at Ohio State University. These interactions improved my understanding of machine learning and helped me write a better thesis. \paragraph{} Finally, I want to thank Katharina Probst, my girlfriend for the past two years. Her help and support made my stay at Carnegie Mellon very pleasant, despite the many stressful factors in the life of a PhD student. I also thank her for patiently putting up with me during these years. \newpage\tableofcontents\newpage \section*{Mathematical Notation} Following is a list of the most frequent mathematical notations in the thesis. \paragraph{} \begin{tabular}{ll} $\cX$ & Set/manifold of data points\\ $\Theta$ & Set/manifold of statistical models\\ $X\times Y$ & Cartesian product of sets/manifold\\ $X^k$ & The Cartesian product of $X$ with itself $k$ times\\ $\|\cdot\|$ & The Euclidean norm \\ $\dotp{\cdot}{\cdot}$ & Euclidean dot product between two vectors\\ $p(x\,;\theta)$ & A probability model for $x$ parameterized by $\theta$\\ $p(y|x\,;\theta)$ & A probability model for $y$, conditioned on $x$ and parameterized by $\theta$\\ $D(\cdot,\cdot),D_r(\cdot,\cdot)$ & $I$-divergence between two models\\ $T_x \cM$ & The tangent space to the manifold $\cM$ at $x\in\cM$\\ $g_x(u,v)$ & A Riemannian metric at $x$ associated with the tangent vectors $u,v\in T_x \cM$\\ $\{\partial_i\}_{i=1}^n,\{e_i\}_{i=1}^n$ & The standard basis associated with the vector space $T_x\cM\cong \R^n$\\ $G(x)$ & The Gram matrix of the metric $g$, $[G(x)]_{ij}=g_x(e_i,e_j)$\\ $\cJ_{\theta}(u,v)$ & The Fisher information metric at the model $\theta$ associated with the vectors $u,v$\\ $\delta_{x}(u,v)$ & The induced Euclidean local metric $\delta_x(u,v)=\dotp{u}{v}$\\ $\delta_{x,y}$ & Kronecker's delta $\delta_{x,y}=1$ if $x=y$ and 0 otherwise\\ $\iota:A\to X$ & The inclusion map $\iota(x)=x$ from $A\subset X$ to $X$.\\ $\mathbb{N},\mathbb{Q},\mathbb{R}$ & The natural, rational and real numbers respectively\\ $\R_{+}$ & The set of positive real numbers\\ $\R^{k\times m}$ & The set of real $k\times m$ matrices\\ $[A]_i$ & The $i$-row of the matrix $A$\\ $\{\partial_{ab}\}_{a,b=1}^{k,m}$ & The standard basis associated with $T_x \R^{k\times m}$\\ $\overline X$ & The topological closure of $X$\\ $\H^n$ & The upper half plane $\{x\in\R^n:x_n\geq 0\}$\\ $\Sn$ & The $n$-sphere $\Sn=\{x\in\R^{n+1}:\sum_i x_i^2=1\}$\\ $\Snp$ & The positive orthant of the $n$-sphere $\Snp=\{x\in\R^{n+1}_{+}:\sum_i x_i^2=1\}$\\ $\Pn$ & The $n$-simplex $\Pn=\{x\in\R_{+}^{n+1}:\sum_i x_i=1\}$\\ $f\circ g$ & Function composition $f\circ g(x)=f(g(x))$\\ $C^{\infty}(\cM,\cN)$ & The set of infinitely differentiable functions from $\cM$ to $\cN$\\ $f_*u$ & The push-forward map $f_*:T_x\cM\to T_{f(x)}\cN$ of $u$ associated with $f:\cM\to\cN$\\ $f^*g$ & The pull-back metric on $\cM$ associated with $(\cN,g)$ and $f:\cM\to \cN$\\ $d_g(\cdot,\cdot),d(\cdot,\cdot)$ & The geodesic distance associated with the metric $g$\\ $\dvol g(\theta)$ & The volume element of the metric $g_{\theta}$, $\dvol g(\theta)=\sqrt{\det g_{\theta}}=\sqrt{\det G(\theta)}$\\ $\nabla_X Y$ & The covariant derivative of the vector field $Y$ in the direction \\ & of the vector field $X$ associated with the connection $\nabla$ \end{tabular} \newpage \begin{tabular}{ll} $\ell(\theta)$ & The log-likelihood function\\ $\cE(\theta)$ & The AdaBoost.M2 exponential loss function\\ $\tilde p$ & Empirical distribution associated with a training set $D\subset \cX\times \cY$\\ $\hat u$ & The $L^2$ normalized version of the vector $u$\\ $d(x,S)$ & The distance from a point to a set $d(x,S)=\inf_{y\in S}d(x,y)$\\ \end{tabular} \newpage %----------------------------------------------------------------- %----------------------------------------------------------------- %----------------------------------------------------------------- \section{Introduction} \label{sec:introduction} There are two fundamental spaces in machine learning. The first space $\cX$ consists of data points and the second space $\Theta$ consists of possible learning models. In statistical learning, $\Theta$ is usually a space of statistical models, $\{p(x\,;\theta):\theta\in\Theta\}$ in the generative case or $\{p(y\given x\,;\theta):\theta\in\Theta\}$ in the discriminative case. The space $\Theta$ can be either a low dimensional parametric space as in parametric statistics, or the space of all possible models as in non-parametric statistics. \paragraph{} Learning algorithms select a model $\theta\in\Theta$ based on a training sample $\{x_i\}_{i=1}^n\subset \cX$ in the generative case or $\{(x_i,y_i)\}_{i=1}^n\subset \cX\times\mathcal{Y}$ in the discriminative case. In doing so they, either implicitly or explicitly, make certain assumptions about the geometries of $\cX$ and $\Theta$. In the supervised case, we focus on the classification setting, where $\mathcal{Y}=\{y_1,\ldots,y_k\}$ is a finite set of unordered classes. By this we mean that the space \[\cX\times\mathcal{Y}=\cX\times\cdots\times\cX =\cX^k\] has the product geometry over $\cX$. This is a common assumption that makes sense in many practical situations, where there is no clear relation or hierarchy between the classes. As a result, we will mostly ignore the role of $\mathcal{Y}$ and restrict our study of data spaces to $\cX$. \paragraph{} Data and model spaces are rarely Euclidean spaces. In data space, there is rarely any meaning to adding or subtracting two data points or multiplying a data point by a real scalar. For example, most representations of text documents as vectors are non-negative and multiplying them by a negative scalar does not yield another document. Similarly, images $I$ are usually represented as matrices whose entries are in some bounded interval of the real line $I\in [a,b]^{k\times m}$ and there is no meaning to matrices with values outside that range that are obtained by addition or scalar multiplication. The situation is similar in model spaces. For example, typical parametric families such as normal, exponential, Dirichlet or multinomial, as well as the set of non-negative, normalized distributions are not Euclidean spaces. \paragraph{} In addition to the fact that data and model spaces are rarely $\R^n$ as topological spaces, the geometric structure of Euclidean spaces, expressed through the Euclidean distance \[ \|x-y\| \deff \sqrt{\sum_i |x_i-y_i|^2}\] is artificial on most data and model spaces. This holds even in many cases when the data or models are real vectors. To study the geometry of $\cX$ and $\Theta$, it is essential to abandon the realm of Euclidean geometry in favor of a new, more flexible class of geometries. The immediate generalization of Euclidean spaces, Banach and Hilbert spaces, are still vector spaces, and by the arguments above, are not a good model for $\cX$ and $\Theta$. Furthermore, the geometries of Banach and Hilbert spaces are quite restricted, as is evident from the undesirable linear scaling of the distance \[ d_E(cx,cy)=|c|\, d_E(x,y) \qquad \forall c\in\R.\] \paragraph{} Despite the fact that most data and model spaces are not Euclidean, they share two important properties: they are smooth and they are locally Euclidean. Manifolds are the natural generalization of Euclidean spaces to locally Euclidean spaces and differentiable manifolds are their smooth counterparts. Riemannian geometry is a mathematical theory of geometries on such smooth, locally Euclidean spaces. In this framework, the geometry of a space is specified by a local inner product $g_x(\cdot,\cdot), x\in\cX$ between tangent vectors called the Riemannian metric. This inner product translates into familiar geometric quantities such as distance, curvature and angles. Using the Riemannian geometric approach to study the geometries of $\cX$ and $\Theta$ allows us to draw upon the vast mathematical literature in this topic. Furthermore, it is an adequate framework as it encompasses most commonly used geometries in machine learning. For example, Euclidean geometry on $\cX=\mathbb{R}^n$ is achieved by setting the local metric $g_x(u,v), x\in\cX$ to be the Euclidean inner product \[\delta_x(u,v) \deff \dotp{u}{v} \deff \sum_{i=1}^n u_iv_i.\] The information geometry on a space of statistical models $\Theta\subset\mathbb{R}^n$ is achieved by setting the metric $g_{\theta}(u,v), \theta\in\Theta$ to be the Fisher information $\cJ_{\theta}(u,v)$ \begin{align} \cJ_{\theta}(u,v) &\deff \sum_{i=1}^n\sum_{j=1}^n u_i v_j \,\int p(x\,;\theta) \,\frac{\partial}{\partial \theta_i}\log p(x\,;\theta) \,\frac{\partial}{\partial \theta_j}\log p(x\,;\theta)\,\mbox{d} x \end{align} where the above integral is replaced with a sum if $\cX$ is discrete. \paragraph{} This thesis is concerned with applying tools from Riemannian geometry to study the relationship between statistical learning algorithms and different geometries of $\cX$ and $\Theta$. As a result, we gain considerable insight into current learning algorithms and we are able to design powerful new techniques that often outperform the current state-of-the-art. \paragraph{} We start the thesis with Section~\ref{sec:RiemannianGeometry} that contains background from Riemannian geometry that is relevant to most subsequent sections. Additional background that is relevant to a specific section will appear in that section alone. The treatment in Section~\ref{sec:RiemannianGeometry} is short and often not rigorous. For a more complete description most textbooks in Riemannian geometry will suffice. Section~\ref{sec:previousWork} gives an overview of relevant research in the interface of machine learning and statistics, and Riemannian geometry and Section~\ref{sec:probSimplexGeom} applies the mathematical theory of Section~\ref{sec:RiemannianGeometry} to spaces of probability models. \paragraph{} In Section~\ref{sec:boostMl} we study the geometry of the space $\Theta$ of conditional models underlying the algorithms logistic regression and AdaBoost. We prove the surprising result that both algorithms solve the same primal optimization problem with the only difference being that AdaBoost lacks normalization constraints, hence resulting in non-normalized models. Furthermore, we show that both algorithms implicitly minimize the conditional $I$-divergence \begin{align*} D(p,q) = \sum_{i=1}^n \sum_y \left( p(y| x_i) \log \frac{p(y| x_i)}{q(y| x_i)} - p(y| x_i) + q(y| x_i) \right) \end{align*} to a uniform model $q$. Despite the fact that the $I$-divergence is not a metric distance function, it is related to a distance under a specific geometry on $\Theta$. This geometry is the product Fisher information geometry whose study is pursued in Section~\ref{sec:axiomaticGeometry}. \paragraph{} By generalizing the theorems of \v{C}encov and Campbell, Section~\ref{sec:axiomaticGeometry} shows that the only geometry on the space of conditional distributions consistent with a basic set of axioms is the product Fisher information. The results of Sections~\ref{sec:boostMl} and \ref{sec:axiomaticGeometry}, provide an axiomatic characterization of the geometries underlying logistic regression and AdaBoost. Apart from providing a substantial new understanding of logistic regression and AdaBoost, this analysis provides, for the first time a theory of information geometry for conditional models. \paragraph{} The axiomatic framework mentioned above provides a natural geometry on the space of distributions $\Theta$. It is less clear what should be an appropriate geometry for the data space $\cX$. A common pitfall shared by many classifiers is to assume that $\cX$ should be endowed with a Euclidean geometry. Many algorithms, such as radial basis machines and Euclidean $k$-nearest neighbor, make this assumption explicit. On the other hand, Euclidean geometry is implicitly assumed in linear classifiers such as logistic regression, linear support vector machines, boosting and the perceptron. Careful selection of an appropriate geometry for $\cX$ and designing classifiers that respect it should produce better results than the naive Euclidean approach. \paragraph{} In Section~\ref{sec:embedding} we propose the following principle to obtain a natural geometry on the data space. By assuming the data is generated by statistical models in the space $\cM$, we embed data points $x\in\cX$ into $\cM$ and thus obtain a natural geometry on the data space -- namely the Fisher geometry on $\cM$. For example, consider the data space $\cX$ of text documents in normalized term frequency (tf) representation embedded in the space of all multinomial models $\cM$. Assuming the documents are generated from multinomial distributions we obtain the maximum likelihood embedding $\hat{\theta}_{MLE}:\cX\to\cM$ which is equivalent to the inclusion map \[\iota:\cX\hookrightarrow\cM, \qquad \iota(x)=x.\] Turning to \v{C}encov's theorem and selecting the Fisher geometry on $\cM$ we obtain a natural geometry on the closure of the data space $\overline{\cX}=\cM$. \paragraph{} The embedding principle leads to a general framework for adapting existing algorithms to the Fisher geometry. In Section~\ref{sec:hyperplane} we generalize the notion of linear classifiers to non-Euclidean geometries and derive in detail the multinomial analogue of logistic regression. Similarly, in Section~\ref{sec:diffusion} we generalize radial basis kernel to non-Euclidean geometries by approximating the solution of the geometric diffusion equation. In both cases, the resulting non-Euclidean generalizations outperform its Euclidean counterpart, as measured by classification accuracy in several text classification tasks. \paragraph{} The Fisher geometry is a natural choice if the only known information is the statistical family that generates the data. In the presence of actual data $\{x_1,\ldots,x_n\}\subset \cX$ it may be possible to induce a geometry that is better suited for it. In Section~\ref{sec:metricLearning} we formulate a learning principle for the geometry of $\cX$ that is based on maximizing the inverse volume of the given training set. When applied to the space of text documents in tf representation, the learned geometry is similar to, but outperforms the popular tf-idf geometry. \paragraph{} We conclude with a discussion in Section~\ref{sec:discussion}. The first two appendices contain technical information relevant to Sections~\ref{sec:boostMl} and \ref{sec:metricLearning}. Appendix~\ref{sec:contributionSummary} contains a summary of the major contributions included in this thesis, along with relevant publications. %---------------------------------------------------------------- \section{Relevant Concepts from Riemannian Geometry} \label{sec:RiemannianGeometry} In this section we describe concepts from Riemannian geometry that are relevant to most of the thesis, with the possible exception of Section~\ref{sec:boostMl}. Other concepts from Riemannian geometry that are useful only for a specific section will be introduced later in the thesis as needed. For more details refer to any textbook discussing Riemannian geometry. \emcite{Milnor:63} and \emcite{Spivak1975} are particularly well known classics and \emcite{Lee2002} is a well-written contemporary textbook. \paragraph{} Riemannian manifolds are built out of three layers of structure. The topological layer is suitable for treating topological notions such as continuity and convergence. The differentiable layer allows extending the notion of differentiability to the manifold and the Riemannian layer defines rigid geometrical quantities such as distances, angles and curvature on the manifold. In accordance with this philosophy, we start below with the definition of topological manifold and quickly proceed to defining differentiable manifolds and Riemannian manifolds. \subsection{Topological and Differentiable Manifolds} A homeomorphism between two topological spaces $X$ and $Y$ is a bijection $\phi:X\to Y$ for which both $\phi$ and $\phi^{-1}$ are continuous. We then say that $X$ and $Y$ are homeomorphic and essentially equivalent from a topological perspective. An $n$-dimensional topological manifold $\cM$ is a topological subspace of $\mathbb{R}^m, m\geq n$ that is locally equivalent to $\mathbb{R}^n$ i.e. for every point $x\in\cM$ there exists an open neighborhood $U\subset \cM$ that is homeomorphic to $\mathbb{R}^n$. The local homeomorphisms in the above definition $\phi_U:U\subset\cM\to\mathbb{R}^n$ are usually called charts. Note that this definition of a topological manifold makes use of an ambient Euclidean space $\mathbb{R}^m$. While sufficient for our purposes, such a reference to $\mathbb{R}^m$ is not strictly necessary and may be discarded at the cost of certain topological assumptions\footnote{The general definition, that uses the Hausdorff and second countability properties, is equivalent to the ambient Euclidean space definition by Whitney's embedding theorem. Nevertheless, it is considerably more elegant to do away with the excess baggage of an ambient space.} \cite{Lee2000}. Unless otherwise noted, for the remainder of this section we assume that all manifolds are of dimension $n$. \paragraph{} An $n$-dimensional topological manifold with a boundary is defined similarly to an $n$-dimensional topological manifold, except that each point has a neighborhood that is homeomorphic to an open subset of the upper half plane \begin{align*} \label{eq:upperHalfPlane} \H^n \deff \left\{x\in\mathbb{R}^n: x_n\geq 0\right\}. \end{align*} It is possible to show that in this case some points $x\in\cM$ have neighborhoods homeomorphic to $U\subset \mathbb{H}^n$ such that $\forall y\in U, y_n>0$ while other points are homeomorphic to a subset $U\subset \mathbb{H}^n$ that intersects the line $y_n=0$. These two sets of points are disjoint and are called the interior and boundary of the manifold $\cM$ and are denoted by $\text{Int}\,\cM$ and $\partial \cM$ respectively \cite{Lee2000}. \paragraph{} Figure \ref{fig:manifoldWithBoundary} illustrates the concepts associated with a manifold with a boundary. Note that a manifold is a manifold with a boundary but the converse does not hold in general. However, if $\cM$ is an $n$-dimensional manifold with a boundary then $\text{Int}\cM$ is an $n$ dimensional manifold and $\partial \cM$ is an $n-1$ dimensional manifold. The above definition of boundary and interior of a manifold may differ from the topological notions of boundary and interior, associated with an ambient topological space. When in doubt, we will specify whether we refer to the manifold or topological interior and boundary. \begin{figure} \centering { \psfrag{m}{$\cM$} \psfrag{g}{$x$} \psfrag{h}{$y$} \psfrag{f}{$\mathbb{H}^2$} \includegraphics[scale=0.57]{RiemannianGeomFigures/manifoldWithBoundary.eps}} \caption{A 2-dimensional manifold with a boundary. The boundary $\partial\cM$ is marked by a black contour. For example, $x$ is a boundary point $x\in\partial \cM$ while $y\in\text{Int}\,\cM$ is an interior point.} \label{fig:manifoldWithBoundary} \end{figure} We return to manifolds with boundary at the end of this Section. \paragraph{} We are now in a position to introduce the differentiable structure. First recall that a mapping between two open sets of Euclidean spaces $f:U\subset \mathbb{R}^k\to V\subset\mathbb{R}^l$ is infinitely differentiable, denoted by $f\in C^{\infty}(\mathbb{R}^k,\mathbb{R}^l)$ if $f$ has continuous partial derivatives of all orders. If for every pair of charts $\phi_U,\phi_V$ the transition function defined by \[\psi:\phi_V(U\cap V)\subset \mathbb{R}^n\to\mathbb{R}^n, \quad \psi=\phi_U\circ\phi_V^{-1}\] (when $U\cap V\neq \emptyset$) is a $C^{\infty}(\mathbb{R}^n,\mathbb{R}^n)$ differentiable map then $\cM$ is called an $n$-dimensional differentiable manifold. The charts and transition function for a 2-dimensional manifold are illustrated in Figure~\ref{fig:transitionFunction}. \begin{figure} \centering { \psfrag{m}{$\cM$} \psfrag{h}{$\phi_V:V\to\mathbb{R}^2$} \psfrag{g}{$\phi_U:U\to\mathbb{R}^2$} \psfrag{f}{$\text{dom}\,\psi=\phi_V(U\cap V)$} \psfrag{k}{$\psi=\phi_U\circ\phi_V^{-1}$} \psfrag{u}{$U$} \psfrag{v}{$V$} \psfrag{l}{$\text{range}\,\psi=\phi_U(U\cap V)$} \includegraphics[scale=0.57]{RiemannianGeomFigures/transitionFunction.eps}} \caption{Two neighborhoods $U,V$ in a 2-dimensional manifold $\cM$, the coordinate charts $\phi_U,\phi_V$ and the transition function $\psi$ between them.} \label{fig:transitionFunction} \end{figure} \paragraph{} Differentiable manifolds of dimensions 1 and 2 may be visualized as smooth curves and surfaces in Euclidean space. Examples of $n$-dimensional differentiable manifolds are the Euclidean space $\mathbb{R}^n$, the $n$-sphere \begin{align}\label{eq:sphereDef} \Sn \deff \left\{x\in\mathbb{R}^{n+1}:\sum_{i=1}^n x_i^2=1\right\}, \end{align} its positive orthant \begin{align}\label{eq:posSphereDef} \Snp \deff \left\{x\in\mathbb{R}^{n+1}:\sum_{i=1}^n x_i^2=1, \,\, \forall i\,\, x_i>0\right\}, \end{align} and the $n$-simplex \begin{align}\label{eq:simplexDef} \mathbb{P}_n \deff \left\{x\in\mathbb{R}^{n+1}: \sum_{i=1}^n x_i=1, \,\, \forall i\,\, x_i>0\right\}. \end{align} We will keep these examples in mind as they will keep appearing throughout the thesis. \paragraph{} Using the charts, we can extend the definition of differentiable maps to real valued functions on manifolds $f:\cM\to\mathbb{R}$ and functions from one manifold to another $f:\cM\to\cN$. A continuous function $f:\cM\to\mathbb{R}$ is said to be $C^{\infty}(\cM,\mathbb{R})$ differentiable if for every chart $\phi_U$ the function $f\circ \phi_U^{-1}\in C^{\infty}(\mathbb{R}^n,\mathbb{R})$. A continuous mapping between two differentiable manifolds $f:\cM\to\cN$ is said to be $C^{\infty}(\mathcal{M,N})$ differentiable if \[\forall r\in C^{\infty}(\cN,\mathbb{R}), \qquad r\circ f\in C^{\infty}(\cM,\mathbb{R}).\] A diffeomorphism between two manifolds $\mathcal{M,N}$ is a bijection $f:\cM\to\cN$ such that $f\in C^{\infty}(\mathcal{M,N})$ and $f^{-1}\in C^{\infty}(\mathcal{N,M})$. \subsection{The Tangent Space} For every point $x\in\cM$, we define an $n$-dimensional real vector space $T_x\cM$, isomorphic to $\R^n$, called the tangent space. The elements of the tangent space, the tangent vectors $v\in T_x \cM$, are usually defined as directional derivatives at $x$ operating on $C^{\infty}(\cM,\mathbb{R})$ differentiable functions or as equivalence classes of curves having the same velocity vectors at $x$ \cite{Spivak1975,Lee2002}. Intuitively, tangent spaces and tangent vectors are a generalization of geometric tangent vectors and spaces for smooth curves and two dimensional surfaces in the ambient $\mathbb{R}^3$. For an $n$-dimensional manifold $\cM$ embedded in an ambient $\R^m$ the tangent space $T_x \cM$ is a copy of $\R^n$ translated so that its origin is positioned at $x$. See Figure~\ref{fig:tangentSpaces} for an illustration of this concept for two dimensional manifolds in $\R^3$. \begin{figure} \centering { \psfrag{p}{$\mathbb{P}_2$} \psfrag{s}{$\mathbb{S}_2$} \psfrag{x}{$T_x\mathbb{P}_2$} \psfrag{t}{$T_x\mathbb{S}_2$} \includegraphics[scale=0.7]{RiemannianGeomFigures/tangentSpaces.eps}} \caption{Tangent spaces of the 2-simplex $T_x\mathbb{P}_2$ and the 2-sphere $T_x\mathbb{S}_2$.} \label{fig:tangentSpaces} \end{figure} \paragraph{} In many cases the manifold $\cM$ is a submanifold of an $m$-dimensional manifold $\cN$, $m\geq n$. Considering $\cM$ and its ambient space $\mathbb{R}^m, m\geq n$ is one special case of this phenomenon. For example, both $\Pn$ and $\Sn$ defined in \eqref{eq:sphereDef},\eqref{eq:simplexDef} are submanifolds of $\mathbb{R}^{n+1}$. In these cases, the tangent space of the submanifold $T_x\cM$ is a vector subspace of $T_x\cN\cong \mathbb{R}^m$ and we may represent tangent vectors $v\in T_x\cM$ in the standard basis $\{\partial_i\}_{i=1}^m$ of the embedding tangent space $T_x\mathbb{R}^m$ as $v=\sum_{i=1}^m v_i \partial_i$. For example, for the simplex and the sphere we have (see Figure~\ref{fig:tangentSpaces}) \begin{align}\label{eq:simplexTangentSpace} T_x\mathbb{P}_n=\left\{ v\in\mathbb{R}^{n+1}:\sum_{i=1}^{n+1} v_i=0\right\} \qquad \qquad T_x\mathbb{S}_n=\left\{ v\in\mathbb{R}^{n+1}:\sum_{i=1}^{n+1} v_ix_i=0\right\}. \end{align} \paragraph{} A $C^{\infty}$ vector field\footnote{The precise definition of a $C^{\infty}$ vector field requires the definition of the tangent bundle. We do not give this definition since it is somewhat technical and does not contribute much to our discussion.} $X$ on $\cM$ is a smooth assignment of tangent vectors to each point of $\cM$. We denote the set of vector fields on $\cM$ as $\mathfrak{X}(\cM)$ and $X_p$ is the value of the vector field $X$ at $p\in\cM$. Given a function $f\in C^{\infty}(\cM,\R)$ we define the action of $X\in \mathfrak{X}(\cM)$ on $f$ as \[ Xf\in C^{\infty}(\cM,\R) \qquad (Xf)(p) = X_p(f)\] in accordance with our definition of tangent vectors as directional derivatives of functions. \subsection{Riemannian Manifolds} A Riemannian manifold $(\cM,g)$ is a differentiable manifold $\cM$ equipped with a Riemannian metric $g$. The metric $g$ is defined by a local inner product on tangent vectors \[g_x(\cdot,\cdot): T_x\cM\times T_x\cM\to\mathbb{R}, \quad x\in\cM\] that is symmetric, bi-linear and positive definite\\ \begin{align*} g_x(u,v)&=g_x(v,u)\\ g_x\left(\sum_{i=1}^n u_i,\sum_{i=1}^n v_i\right)&=\sum_{i=1}^n\sum_{j=1}^n g_x(u_i,v_j)\\ g_x(u,u) &\geq 0\\ g_x(u,u) &=0 \Leftrightarrow u=0 \end{align*} and is also $C^{\infty}$ differentiable in $x$. By the bi-linearity of the inner product $g$, for every $u,v\in T_x\cM$ \[ g_x(v,u) = \sum_{i=1}^n\sum_{j=1}^n v_i u_j g_x(\partial_i,\partial_j)\] and $g_x$ is completely described by $\{g_x(\partial_i,\partial_j):1\leq i,j\leq n\} $ -- the set of inner products between the basis elements $\{\partial_i\}_{i=1}^n$ of $T_x\cM$. The Gram matrix $[G(x)]_{ij}=g_x(\partial_i,\partial_j)$ is a symmetric and positive definite matrix that completely describes the metric $g_x$. \paragraph{} The metric enables us to define lengths of tangent vectors $v\in T_x \cM$ by $\sqrt{g_x(v,v)}$ and lengths of curves $\gamma:[a,b]\to\cM$ by \[L(\gamma)=\int_a^b \sqrt{g_x( \dot{\gamma}(t), \dot{\gamma}(t))} %\|\dotp{\gamma}(t)\| \text{d}t\] where $\dot{\gamma}(t)$ is the velocity vector of the curve $\gamma$ at time $t$. Using the above definition of lengths of curves, we can define the distance $d_{g}(x,y)$ between two points $x,y\in\cM$ as \[ d_{g}(x,y) = \inf_{\gamma\in \Gamma(x,y)} \int_a^b \sqrt{g_x( \dot{\gamma}(t), \dot{\gamma}(t))} %\|\dotp{\gamma}(t)\| \text{d}t\] where $\Gamma(x,y)$ is the set of piecewise differentiable curves connecting $x$ and $y$. The distance $d_{g}$ is called geodesic distance and the minimal curve achieving it is called a geodesic curve\footnote{It is also common to define geodesics as curves satisfying certain differential equations. The above definition, however, is more intuitive and appropriate for our needs.}. Geodesic distance satisfies the usual requirements of a distance and is compatible with the topological structure of $\cM$ as a topological manifold. If the manifold in question is clear from the context, we will remove the subscript and use $d$ for the geodesic distance. \paragraph{} A manifold is said to be geodesically complete if any geodesic curve $c(t)$, $t\in [a,b]$, can be extended to be defined for all $t\in \R$. It can be shown, that the following are equivalent \begin{itemize} \item $(\cM,g)$ is geodesically complete \item $d_g$ is a complete metric on $\cM$ \item closed and bounded subsets of $\cM$ are compact. \end{itemize} In particular, compact manifolds are geodesically complete. The Hopf-Rinow theorem asserts that if $\cM$ is geodesically complete, then any two points can be joined by a geodesic. \paragraph{} Given two Riemannian manifolds $(\cM,g)$, $(\cN,h)$ and a diffeomorphism between them $f:\cM\to\cN$ we define the push-forward and pull-back maps below\footnote{ The push-forward and pull-back maps may be defined more generally using category theory as covariant and contravariant functors \cite{Lee2000}.} \begin{defn} The push-forward map $f_*:T_x\cM\to T_{f(x)}\cN$, associated with the diffeomorphism $f:\cM\to\cN$ is the mapping that satisfies \[ v(r\circ f)=(f_*v)r, \qquad \forall r\in C^{\infty}(\cN,\mathbb{R}).\] \end{defn} The push-forward is none other than a coordinate free version of the Jacobian matrix $J$ or the total derivative operator associated with the local chart representation of $f$. In other words, if we define the coordinate version of $f:\cM\to\cN$ \[\tilde f=\phi\circ f\circ \psi^{-1}:\mathbb{R}^n\to\mathbb{R}^m\] where $\phi,\psi$ are local charts of $\cN,\cM$ then the push-forward map is \[ f_*u = Ju = \sum_i \left(\sum_j \frac{\partial \tilde f_i}{\partial x_j} u_j \right) e_i \] where $J$ is the Jacobian of $\tilde f$ and $\tilde f_i$ is the $i$-component function of $\tilde f:\R^m\to\R^n$. Intuitively, as illustrated in Figure~\ref{fig:push-forward}, the push-forward transforms velocity vectors of curves $\gamma$ to velocity vectors of transformed curves $f(\gamma)$. \begin{figure} \centering { \psfrag{m}{$\cM$} \psfrag{n}{$\cN$} \psfrag{f}{$f$} \psfrag{p}{$T_x\cM$} \psfrag{q}{$T_{f(x)}\cN$} \includegraphics[scale=0.6]{RiemannianGeomFigures/push-forward.eps} \caption{The map $f:\cM\to\cN$ defines a push forward map $f_*:T_x\cM\to T_{f(x)}\cN$ that transforms velocity vectors of curves to velocity vectors of the transformed curves.} \label{fig:push-forward}} \end{figure} \begin{defn} \label{def:pull-back} Given $(\cN,h)$ and a diffeomorphism $f:\cM\to\cN$ we define a metric $f^*h$ on $\cM$ called the pull-back metric by the relation \[(f^* h)_x(u,v)=h_{f(x)}(f_*u,f_*v).\] \end{defn} \begin{defn} An isometry is a diffeomorphism $f:\cM\to \cN$ between two Riemannian manifolds $(\cM,g),(\cN,h)$ for which \[g_x(u,v)=(f^*h)_x(u,v) \qquad \forall x\in\cM,\qquad \forall u,v\in T_x\cM.\] \end{defn} Isometries, as defined above, identify two Riemannian manifolds as identical in terms of their Riemannian structure. Accordingly, isometries preserve all the geometric properties including the geodesic distance function $d_{g}(x,y)=d_h(f(x),f(y))$. Note that the above definition of an isometry is defined through the local metric in contrast to the global definition of isometry in other branches of mathematical analysis. \paragraph{} A smooth and Riemannian structure may be defined over a topological manifold with a boundary as well. The definition is a relatively straightforward extension using the notion of differentiability of maps between non-open sets in $\mathbb{R}^n$ \cite{Lee2002}. Manifolds with a boundary are important for integration and boundary value problems; our use of them will be restricted to Section~\ref{sec:diffusion}. \paragraph{} In the following section we discuss some previous work and then proceed to examine manifolds of probability distributions and their Fisher geometry. %---------------------------------------------------------------- \section{Previous Work}\label{sec:previousWork} Connections between statistics and Riemannian geometry have been discovered and studied for over fifty years. We give below a roughly chronological overview of this line of research. The overview presented here is not exhaustive. It is intended to review some of the influential research that is closely related to this thesis. For more information on this research area consult the monographs \cite{Murray1993,Kass1997,Amari2000}. Additional previous work that is related to a specific section of this thesis will be discussed inside that section. \paragraph{} \emcite{Rao1945} was the first to point out that a statistical family can be considered as a Riemannian manifold with the metric determined by the Fisher information quantity. \emcite{Efron1975} found a relationship between the geometric curvature of a statistical model and Fisher and Rao's theory of second order efficiency. In this sense, a model with small curvature enjoys nice asymptotic properties and a model with high curvature implies a break-down of these properties. Since Efron's result marks a historic breakthrough of geometry in statistics we describe it below in some detail. \paragraph{} Recall that by Cram\'er-Rao lower bound, the variance of unbiased estimators is bounded from below by the inverse Fisher information. More precisely, in matrix form we have that $\forall \theta$, $V(\theta)-G^{-1}(\theta)$ is positive semi-definite, where $V$ is the covariance matrix of the estimator and $G$ is the Fisher information matrix\footnote{The matrix form of the Cram\'er-Rao lower bound may be written as $V(\theta)-G^{-1}(\theta)\succeq 0$ with equality $V(\theta)=G^{-1}(\theta)$ if the bound is attained.}. Estimators that achieve this bound asymptotically, for all $\theta$, are called asymptotically efficient, or first order efficient. The prime example for such estimators, assuming some regularity conditions, is the maximum likelihood estimator. The most common example of statistical families for which the regularity conditions hold is the exponential family. Furthermore, in this case, the MLE is a sufficient statistic. These nice properties of efficiency and sufficiency of the MLE do not hold in general for non-exponential families. \paragraph{} The term second order efficiency was coined by Rao, and refers to a subset of consistent and first order efficient estimators $\hat\theta(x_1,\ldots,x_n)$ that attain equality in the following general inequality \begin{align} \label{eq:secOrder} \lim_{n\to\infty} n(i(\theta^{\text{true}})-i(\hat\theta(x_1,\ldots,x_n))) \geq i(\theta^{\text{true}}) \gamma^2(\theta^{\text{true}}) \end{align} where $i$ is the one-dimensional Fisher information and $\gamma$ some function that depends on $\theta$. Here $x_1,\ldots,x_n$ are sampled iid from $p(x\,;\theta^{\text{true}})$ and $\hat\theta(x_1,\ldots,x_n)$ is an estimate of $\theta^{\text{true}}$ generated by the estimator $\hat\theta$. The left hand side of \eqref{eq:secOrder} may be interpreted as the asymptotic rate of the loss of information incurred by using the estimated parameter rather than the true parameter. It turns out that the only consistent asymptotically efficient estimator that achieves equality in \eqref{eq:secOrder} is the MLE, thereby giving it a preferred place among the class of first order efficient estimators. \paragraph{} The significance of Efron's result is that he identified the function $\gamma$ in \eqref{eq:secOrder} as the Riemannian curvature of the statistical manifold with respect to the exponential connection. Under this connection, exponential families are flat and their curvature is 0. For non-exponential families the curvature $\gamma$ may be interpreted as measuring the breakdown of asymptotic properties which is surprisingly similar to the interpretation of curvature as measuring the deviation from flatness, expressing in our case an exponential family. Furthermore, Efron showed that the variance of the MLE in non-exponential families exceeds the Cram\'er-Rao lower bound in approximate proportion to $\gamma^2(\theta)$. \paragraph{} \emcite{Dawid1975} points out that Efron's notion of curvature is based on a connection that is not the natural one with respect to the Fisher geometry. A similar notion of curvature may be defined for other connections, and in particular for the Riemannian connection that is compatible with the Fisher information metric. Dawid's comment about the possible statistical significance of curvature with respect to the metric connection remains a largely open question, although some results were obtained by \emcite{Brody1998}. \emcite{Chentsov1982} introduced a family of connections, later parameterized by $\alpha$, that include as special cases the exponential connection and the metric connection. Using Amari's parametrization of this family, $\alpha=1$ corresponds to the exponential connection, $\alpha=0$ corresponds to the metric connection and the $\alpha=-1$ corresponds to the mixture connection, under which mixture families enjoy 0 curvature \cite{Amari2000}. \paragraph{} \emcite{Chentsov1982} proved that the Fisher information metric is the only Riemannian metric that is preserved under basic probabilistic transformations. These transformations, called congruent embeddings by a Markov morphism, represent transformations of the event space that is equivalent to extracting a sufficient statistic. Later on, \emcite{Campbell86} extended \v{C}encov's result to non-normalized positive models, thus axiomatically characterizing a set of metrics on the positive cone $\mathbb{R}^n_{+}$. \paragraph{} In his short note \emcite{Dawid1977} extended these ideas to infinite dimensional manifolds representing non-parametric sets of densities. More rigorous studies include \cite{Lafferty1988} that models the manifold of densities on a Hilbert space and \cite{Pistone1995} that models the same manifold on non-reflexive Banach spaces called Orlicz spaces. This latter approach, that does not admit a Riemannian structure, was further extended by Pistone and his collaborators and by \emcite{Grasselli2001}. \paragraph{} Additional research by Barndorff-Nielsen and others considered the connection between geometry and statistics from a different angle. Below is a brief description of some these results. In \cite{Barndorff1986} the expected Fisher information is replaced with the observed Fisher information to provide an alternative geometry for a family of statistical models. The observed geometry metric is useful in obtaining approximate expansion of the distribution of the MLE, conditioned on an ancillary statistic. This result continues previous research by Efron and Amari that provides a geometric interpretation for various terms appearing in the Edgeworth expansion of the distribution of the MLE for curved exponential models. \emcite{Barndorff1983} studied the relation between certain partitions of an exponential family of models and geometric constructs. The partitions, termed affine dual foliations, refers to a geometric variant of the standard division of $\R^n$ into copies of $\R^{n-1}$. Based on differential geometry, \emcite{Barndorff1993} define a set of models called orthogeodesic models that enjoy nice higher-order asymptotic properties. Orthogeodesic models include exponential models with dual affine foliations as well as transformation models and provides an abstract unifying framework for such models. \paragraph{} The geodesic distance under the Fisher metric has been examined in various statistical studies. It is used in statistical testing and estimation as an alternative to Kullback Leibler or Jeffreys divergence \cite{Atkinson1981}. It is also essentially equivalent to the popular Hellinger distance \cite{Beran1977} that plays a strong role in the field of statistical robustness \cite{Tamura1986,Lindsay1994,Cutler1996}. \paragraph{} In a somewhat different line of research, \emcite{Csiszar1975} studied the geometry of probability distributions through the notion of $I$-divergence. In a later paper \emcite{Csiszar1991} showed that $I$-divergence estimation, along with least squares, enjoy nice axiomatic frameworks. While not a distance measure, the $I$-divergence bears close connection to the geodesic distance under the Fisher information metric \cite{Kullback1968}. This fact, together with the prevalence of the $I$-divergence and its special case the Kullback-Leibler divergence, brings these research directions together under the umbrella of information geometry. \emcite{Amari2000} contains some further details on the connection between $I$-divergence and Fisher geometry. \paragraph{} Information geometry arrived somewhat later in the machine learning literature. Most of the studies in this context were done by Amari's group. Amari examines several geometric learning algorithms for neural networks \cite{Amari1995} and shows how to adapt the gradient descent algorithm to information geometry in the context of neural networks \cite{Amari1998} and independent component analysis (ICA) \cite{Amari1999}. \emcite{Ikeda2004} interprets several learning algorithms in graphical models such as belief propagation using information geometry and introduces new variations on them. Further information on Amari's effort in different applications including information theory and quantum estimation theory may be obtained from \cite{Amari2000}. \emcite{Saul1997} interpolates between different models based on differential geometric principles and \emcite{Jaakkola1998} use the Fisher information to enhance a discriminative classifier with generative qualities. \emcite{Gous1998} and \emcite{Hall2000} use information geometry to represent text documents by affine subfamilies of multinomial models. \paragraph{} In the next section we apply the mathematical framework developed in Section~\ref{sec:RiemannianGeometry} to manifolds of probability models. %---------------------------------------------------------------- \section{Geometry of Spaces of Probability Models}\label{sec:probSimplexGeom} Parametric inference in statistics is concerned with a parametric family of distributions $\{p(x\,;\theta):\theta\in\Theta\subset\mathbb{R}^n\}$ over the event space $\cX$. If the parameter space $\Theta$ is a differentiable manifold and the mapping $\theta\mapsto p(x\,;\theta)$ is a diffeomorphism we can identify statistical models in the family as points on the manifold $\Theta$. The Fisher information matrix $E\{s s^{\top}\}$ where $s$ is the gradient of the score $[s]_i=\partial \log p(x\,;\theta)/\partial \theta_i$ may be used to endow $\Theta$ with the following Riemannian metric \begin{align} \nonumber \cJ_{\theta}(u,v) &\deff \sum_{i,j}u_i v_j \int p(x\,;\theta) \frac{\partial}{\partial \theta_i} \log p(x\,;\theta) \frac{\partial}{\partial \theta_j} \log p(x\,;\theta) \text{d}x\\ &= \sum_{i,j}u_i v_j \,\, E\left\{\frac{\partial \log p(x\,;\theta)}{\partial\theta_i} \frac{\partial \log p(x\,;\theta)}{\partial\theta_j}\right\}.\label{eq:FisherInformationMetric} \end{align} If $\cX$ is discrete the above integral is replaced with a sum. An equivalent form of \eqref{eq:FisherInformationMetric} for normalized distributions that is sometimes easier to compute is \begin{align} \nonumber \cJ_{\theta}(u,v) &= \cJ_{\theta}(u,v) - \sum_{ij} u_i v_j \frac{\partial^2}{\partial\theta_i\partial\theta_j} \int p(x\,; \theta) \dd d x \\ \nonumber &= \sum_{ij} u_i v_j\int p(x\,;\theta) \left( \left(\frac{1}{p(x\given \theta)} \frac{\partial p(x\,; \theta)}{\partial \theta_j}\right) \left(\frac{1}{p(x\given \theta)} \frac{\partial p(x\,; \theta)}{\partial \theta_i}\right) -\frac{1}{p(x\,;\theta)} \frac{\partial^2}{\partial\theta_i\partial\theta_j} p(x\,; \theta) \right) \, \dd d x \\ \nonumber &= -\sum_{ij} u_i v_j\int p(x\,;\theta) \frac{\partial }{\partial \theta_j}\frac{1}{p(x\given \theta)}\frac{\partial p(x\,; \theta)}{\partial \theta_i} \, \dd d x \\ &=-\sum_{ij} u_i v_j\int p(x\,;\theta) \frac{\partial^2 }{\partial \theta_j \partial \theta_i}\log p(x\,; \theta) \, \dd d x\nonumber \\ &=\sum_{ij} u_i v_j E\left\{-\frac{\partial^2 }{\partial \theta_j \partial \theta_i}\log p(x\,; \theta)\right\} \label{eq:secDerivFish} \end{align} assuming that we can change the order of the integration and differentiation operators. \paragraph{} In the remainder of this section we examine in detail a few important Fisher geometries. The Fisher geometries of finite dimensional non-parametric space, finite dimensional conditional non-parametric space and spherical normal space are studied next. \subsection{Geometry of Non-Parametric Spaces} In the finite non-parametric setting, the event space $\cX$ is a finite set with $|\cX|=n$ and $\Theta=\P_{n-1}$, defined in \eqref{eq:simplexDef}, which represents the manifold of all positive probability models over $\cX$. The positivity constraint is necessary for $\Theta=\P_{n-1}$ to be a manifold. If zero probabilities are admitted, the appropriate framework for the parameter space $\Theta=\overline{\P_{n-1}}$ is a manifold with corners \cite{Lee2002}. Note that the above space $\Theta$ is precisely the parametric space of the multinomial family. Hence, the results of this section may be interpreted with respect to the space of all positive distributions on a finite event space, or with respect to the parametric space of the multinomial distribution. \paragraph{} The finiteness of $\cX$ is necessary for $\Theta$ to be a finite dimensional manifold. Relaxing the finiteness assumption results in a manifold where each neighborhood is homeomorphic to an infinite dimensional vector space called the model space. Such manifolds are called Frechet, Banach or Hilbert manifolds (depending on the model space) and are the topic of a branch of geometry called global analysis \cite{Lang1999}. \emcite{Dawid1977} remarked that an infinite dimensional non-parametric space may be endowed with multinomial geometry leading to spherical geometry on a Hilbert manifold. More rigorous modeling attempts were made by \emcite{Lafferty1988} that models the manifold of densities on a Hilbert space and by \emcite{Pistone1995} that model it on a non-reflexive Banach space. See also the brief discussion on infinite dimensional manifolds representing densities by \emcite{Amari2000} pp. 44-45. \paragraph{} Considering $\mathbb{P}_{n-1}$ as a submanifold of $\mathbb{R}^n$, we represent tangent vectors $v\in T_{\theta}\mathbb{P}_{n-1}$ in the standard basis of $T_{\theta}\mathbb{R}^{n}$. As mentioned earlier \eqref{eq:simplexTangentSpace}, this results in the following representation of $v\in T_{\theta}\mathbb{P}_{n-1}$ \[v=\sum_{i=1}^n v_i\partial_i \quad \text{ subject to } \quad \sum_{i=1}^n v_i=0.\] Using this representation, the loglikelihood and its derivatives are \begin{align*} \log p(x\,; \theta) &= \sum_{i=1}^{n} x_i \log \theta_i \\ \frac{\partial \log p(x\,; \theta)}{\partial \theta_i} &= \frac{x_i}{\theta_i} \\ \frac{\partial^2 \log p(x\,; \theta)}{\partial \theta_i\partial\theta_j} &= -\frac{x_i}{\theta_i^2}\,\delta_{ij} \end{align*} and using equation \eqref{eq:secDerivFish} the Fisher information metric on $\P_{n-1}$ becomes \begin{align*} \cJ_{\theta}(u,v) &= - \sum_{i=1}^{n}\sum_{j=1}^{n} u_i v_j E\left[ \frac{\partial^2 \log p(x\given \theta)} {\partial \theta_i\partial \theta_j} \right] = - \sum_{i=1}^{n} u_i v_i E\left\{ -x_i/\theta_i^2\right\} = \sum_{i=1}^{n} \frac{u_i v_i}{\theta_i} \end{align*} since $E x_i=\theta_i$. Note that the Fisher metric emphasizes coordinates that correspond to low probabilities. The fact that the metric $\cJ_{\theta}(u,v)\to \infty$ when $\theta_i\to 0$ is not problematic since length of curves that involves integrals over $g$ converge. \paragraph{} While geodesic distances are difficult to compute in general, in the present case we can easily compute the geodesics by observing that the standard Euclidean metric on the surface of the positive $n$-sphere is the pull-back of the Fisher information metric on the simplex. More precisely, the transformation $F(\theta_1,\ldots,\theta_{n+1})=(2\sqrt{\theta_1},\ldots,2\sqrt{\theta_{n+1}})$ is a diffeomorphism of the $n$-simplex $\Pn$ onto the positive portion of the $n$-sphere of radius~2 \[\tilde{\mathbb{S}}_{+}^n = \left\{ \theta\in\R^{n+1} \;:\; \sum_{i=1}^{n+1}\theta_i^2=4,\; \theta_i > 0 \right\}.\] The inverse transformation is \[F^{-1}:\tilde{\mathbb{S}}_{+}^n\to\Pn, \qquad F^{-1}(\theta)=\left(\frac{\theta_1^2}{4},\ldots,\frac{\theta_{n+1}^2}{4}\right)\] and its push-forward is \[ F^{-1}_*(u) = \left(\frac{u_1}{2},\ldots,\frac{u_{n+1}}{2}\right).\] The metric on $\tilde{\mathbb{S}}_{+}^n$ obtained by pulling back the Fisher information on $\Pn$ through $F^{-1}$ is \begin{align*} h_{\theta}(u,v)&= \cJ_{\theta^2/4}\left(F^{-1}_*\sum_{k=1}^{n+1} u_k e_k,F^{-1}_*\sum_{l=1}^{n+1} v_l e_l\right) =\sum_{k=1}^{n+1}\sum_{l=1}^{n+1} u_k v_l\, g_{\theta^2/4}(F^{-1}_*e_k,F^{-1}_*e_l)\\ &= \sum_{k=1}^{n+1} \sum_{l=1}^{n+1} u_k v_l \sum_i \frac{4}{\theta_i^2}\, (F_*^{-1}e_k)_i\,(F_*^{-1}e_l)_i = \sum_{k=1}^{n+1}\sum_{l=1}^{n+1} u_k v_l \sum_i \frac{4}{\theta_i^{2}}\,\frac{\theta_k \delta_{ki}}{2}\, \frac{\theta_l\delta_{li}}{2} =\sum_{i=1}^{n+1} u_i v_i\\ &= \delta_{\theta}(u,v) \end{align*} the Euclidean metric on $\tilde{\mathbb{S}}_{+}^n$ inherited from the embedding Euclidean space $\R^{n+1}$. \paragraph{} Since the transformation $F:(\Pn,\cJ)\to(\tilde{\mathbb{S}}_{+}^n,\delta)$ is an isometry, the geodesic distance $d_{\cJ}(\theta,\theta')$ on $\Pn$ may be computed as the shortest curve on $\tilde{\mathbb{S}}_{+}^n$ connecting $F(\theta)$ and $F(\theta')$. These shortest curves are portions of great circles -- the intersection of a two dimensional subspace and $\tilde{\mathbb{S}}_{+}^n$ whose lengths are \begin{equation} \label{eq:arccos} d_{\cJ}(\theta,\theta') = d_{\delta}(F(\theta),F(\theta')) = 2 \arccos \left(\sum_{i=1}^{n+1} \sqrt{\theta_i\, \theta'_i}\right)\,. \end{equation} We illustrate these geodesic distances in Figures~\ref{fig:vis}-\ref{fig:equalDist2}. Figure~\ref{fig:vis} shows how to picture $\P_2$ as a triangle in $\R^2$ and Figure~\ref{fig:equalDist2} shows the equal distant contours for both Euclidean and Fisher geometries. Will often ignore the factor of 2 in \eqref{eq:arccos} to obtain a more compact notation for the geodesic distance. \begin{figure} \begin{center} {\psfrag{x}{$(1,0,0)$} \psfrag{y}{$(0,1,0)$} \psfrag{z}{$(0,0,1)$} \includegraphics[scale=0.6]{probSpacesFigures/simplex2D-3D.eps}} \end{center} \caption{The 2-simplex $\P_2$ may be visualized as a surface in $\R^3$ (left) or as a triangle in $\R^2$ (right).} \label{fig:vis} \end{figure} \begin{figure} \begin{center} \begin{tabular}{cc} \includegraphics[scale=0.4]{diffusionFigures/heat1}& \includegraphics[scale=0.4]{diffusionFigures/rbf1}\\[16pt] \includegraphics[scale=0.4]{diffusionFigures/heat2}& \includegraphics[scale=0.4]{diffusionFigures/rbf2}\\[16pt] \includegraphics[scale=0.4]{diffusionFigures/heat3} & \includegraphics[scale=0.4]{diffusionFigures/rbf3} \end{tabular} \end{center} \caption{Equal distance contours on $\mathbb{P}_2$ from the upper right edge (top row), the center (center row), and lower right corner (bottom row). The distances are computed using the Fisher information metric (left column) or the Euclidean metric (right column).} \label{fig:equalDist2} \end{figure} \paragraph{} The geodesic distances $d_{\cJ}(\theta,\theta')$ under the Fisher geometry and the Kullback-Leibler divergence $D(\theta,\theta')$ agree up to second order as $\theta\to\theta'$ \cite{Kullback1968}. Similarly, the Hellinger distance \cite{Beran1977} \begin{equation} d_H(\theta,\theta') = \sqrt{\sum_i \left( \sqrt{\theta_i} - \sqrt{\theta'_i}\right)^2} \end{equation} is related to $d_{\cJ}(\theta,\theta')$ by \begin{equation} d_H(\theta,\theta') = 2 \sin\left(d_{\cJ}(\theta,\theta')/4\right) \end{equation} and thus also agrees with the distance up to second order as $\theta'\rightarrow \theta$. \subsection{Geometry of Non-Parametric Conditional Spaces} Given two finite event sets $\cX,\cY$ of sizes $k$ and $m$ respectively, a conditional probability model $p(y|x)$ reduces to an element of $\mathbb{P}_{m-1}$ for each $x\in\cX$. We may thus identify the space of conditional probability models associated with $\cX$ and $\cY$ as the product space \[\P_{m-1}\times\cdots\times\P_{m-1}=\P_{m-1}^{k}.\] \paragraph{} For our purposes, it will be sometimes more convenient to work with the more general case of positive non-normalized conditional models. Dropping the normalization constraints $\sum_i p(y_i|x_j)=1$ we obtain conditional models in the cone of $k\times m$ matrices with positive entries, denoted by $\mathbb{R}_{+}^{k\times m}$. Since a normalized conditional model is also a non-normalized one, we can consider $\P_{m-1}^{k}$ to be a subset of $\mathbb{R}_{+}^{k \times m}$. Results obtained for non-normalized models apply then to normalized models as a special case. In addition, some of the notation and formulation is simplified by working with non-normalized models. \paragraph{} In the interest of simplicity, we will often use matrix notation instead of the standard probabilistic notation. A conditional model (either normalized or non-normalized) is described by a positive matrix $M$ such that $M_{ij}=p(y_j|x_i)$. Matrices that correspond to normalized models are (row) stochastic matrices. We denote tangent vectors to $\mathbb{R}_{+}^{k\times m}$ using the standard basis \[T_M\mathbb{R}_{+}^{k\times m}=\text{span} \{\partial_{ij}:i=1,\ldots,k, j=1,\ldots,m\}.\] Tangent vectors to $\P_{m-1}^{k}$, when expressed using the basis of the embedding tangent space $T_M\mathbb{R}_{+}^{k\times m}$ are linear combinations of $\{\partial_{ij}\}$ such that the sums of the combination coefficients over each row are 0, e.g. \begin{align*} \frac{1}{2}\partial_{11}+\frac{1}{2}\partial_{12}-\partial_{13}+\frac{1}{3}\partial_{21}-\frac{1}{3} \partial_{22}&\in T_M\P_{2}^{3}\\ \frac{1}{2}\partial_{11}+\frac{1}{2}\partial_{12}-\partial_{13}+\frac{1}{3}\partial_{21}-\frac{1}{2} \partial_{22}&\not\in T_M\P_{2}^{3}. \end{align*} \paragraph{} The identification of the space of conditional models as a product of simplexes demonstrates the topological and differentiable structure. In particular, we do not assume that the metric has a product form. However, it is instructive to consider as a special case the product Fisher information metric on $\P_{n-1}^k$ and $\mathbb{R}_{+}^{k\times m}$. Using the above representation of tangent vectors $u,v\in T_M\mathbb{R}_{+}^{k\times m}$ or $u,v\in T_M\mathbb{P}_{m-1}^{k}$ the product Fisher information \begin{align*} \cJ^k_M(u_1\oplus\cdots\oplus u_k,v_1\oplus\cdots\oplus v_k)&\deff (\cJ\otimes \cdots \otimes\cJ)_M(u_1\oplus\cdots\oplus u_k,v_1\oplus\cdots\oplus v_k)\\ &\deff \sum_{i=1}^k \cJ_{[M]_i}(u_i,v_i), \end{align*} where $[A]_i$ is the $i$-row of the matrix $A$, $\otimes$ is the tensor product and $\oplus$ is the direct sum decomposition of vectors, reduces to \begin{align}\label{eq:productFisherMetric1} \cJ^k_M(u,v) &= \sum_{i=1}^k\sum_{j=1}^m \frac{u_{ij} v_{ij}}{M_{ij}}. \end{align} A different way of expressing \eqref{eq:productFisherMetric1} is by specifying the values of the metric on pairs of basis elements \begin{align}\label{eq:productFisherMetric2} g_M(\partial_{ab},\partial_{cd}) &=\delta_{ac}\delta_{bd}\frac{1}{M_{ab}} \end{align} where $\delta_{ab}=1$ if $a=b$ and 0 otherwise. %---------------------- \subsection{Geometry of Spherical Normal Spaces} \label{sec:geomSphericalNormal} Given a restricted parametric family $\Theta\subset \Pn$ the Fisher information metric on $\Theta$ agrees with the induced metric from the Fisher information metric on $\Pn$. When $\cX$ is infinite and $\Theta$ is a finite dimensional parametric family, we can still define the Fisher information metric $\cJ$ on $\Theta$, however without a reference to an embedding non-parametric space. We use this approach to consider the Fisher geometry of the spherical normal distributions on $\cX=\R^{n-1}$ \[ \{\cN(\mu,\sigma I): \mu\in\R^{n-1},\sigma\in \R_{+}\}\] parameterized by the upper half plane $\Theta=\H^n\cong \R^{n-1}\times \R_{+}$. \paragraph{} To compute the Fisher information metric for this family, it is convenient to use the expression given by equation \eqref{eq:secDerivFish}. Then simple calculations yield, for $1 \leq i,j \leq n-1$ \begin{align*} [G(\theta)]_{ij} &= -\int_{\R^{n-1}} \frac{\partial^2}{\partial {\mu_i}\partial {\mu_j}} \left( - \sum_{k=1}^{n-1} \frac{(x_k-\mu_k)^2}{2\sigma^2}\right) p(x\given\theta) \, dx = \frac{1}{\sigma^2} \,\delta_{ij} \\ [G(\theta)]_{ni} &= -\int_{\R^{n-1}} \frac{\partial^2}{\partial \sigma\partial {\mu_i}} \left( - \sum_{k=1}^{n-1} \frac{(x_k-\mu_k)^2}{2\sigma^2}\right) p(x\given\theta) \, dx = \frac{2}{\sigma^3} \int_{\R^{n-1}} (x_i - \mu_i) \, p(x\given\theta)\, dx = 0 \\ [G(\theta)]_{nn} &= -\int_{\R^{n-1}} \frac{\partial^2}{\partial \sigma\partial \sigma} \left( - \sum_{k=1}^{n-1} \frac{(x_k-\mu_k)^2}{2\sigma^2} - (n-1)\log \sigma\right) p(x\given\theta) \, dx \\ &= \frac{3}{\sigma^4} \int_{\R^{n-1}} \sum_{k=1}^{n-1}(x_k - \mu_k)^2 \, p(x\given\theta)\, dx - \frac{n-1}{\sigma^2} = \frac{2(n-1)}{\sigma^2}. \end{align*} \paragraph{} Letting $\theta'$ be new coordinates defined by $\theta'_i = \mu_i$ for $1\leq i\leq n-1$ and $\theta'_n = \sqrt{2(n-1)}\, \sigma$, we see that the Gram matrix is given by \begin{equation} [G(\theta')]_{ij} = \frac{1}{\sigma^2} \,\delta_{ij} \end{equation} and the Fisher information metric gives the spherical normal manifold a hyperbolic geometry\footnote{The manifold $\H^n$ with with hyperbolic geometry is often referred to as Poincar\'{e}'s upper half plane and is a space of constant negative curvature.}. It is shown by \emcite{Kass1997} that any location-scale family of densities \begin{equation*} p(x\,; \mu,\sigma) = \frac{1}{\sigma} f\left(\frac{x-\mu}{\sigma}\right) \quad (\mu,\sigma)\in\R\times\R_+ \quad f:\R\rightarrow \R \end{equation*} have a similar hyperbolic geometry. The geodesic curves in the two dimensional hyperbolic space are circles whose centers lie on the line $x_2=0$ or vertical lines (considered as circles whose centers lie on the line $x_2=0$ with infinite radius) \cite{Lee1997}. An illustration of these curves appear in Figure~\ref{fig:hypGeodesics}. \begin{figure} \begin{center} \includegraphics[scale=0.5]{probSpacesFigures/hyperbolic-Geodesics.eps} \caption{ Geodesic curves in $\H^2$ with the hyperbolic metric are circles whose centers lie on the line $x_2=0$ or vertical lines.} \label{fig:hypGeodesics} \end{center} \end{figure} To compute the geodesic distance on $\H^n$ we transform points in $\H^n$ to an isometric manifold known as Poincar\'e's ball. We first define the sphere inversion of $x$ with respect to a sphere $S$ with center $a$ and radius $r$ as \[I_S(x)=\frac{r^2}{\|x-a\|^2}(x-a)+a.\] The Cayley transform is the sphere inversion with respect to a sphere with center $(0,\ldots,0,-1)$ and radius $\sqrt{2}$. We denote by $\eta$ the inverse of the Cayley's transform that maps the hyperbolic half plane to Poincar\'e's ball. \[ \eta(x) = -I_{S'}(-x)\qquad x\in\H^n \] where $S'$ is a sphere with center at $(0,\ldots,0,1)$ and radius $\sqrt{2}$. The geodesic distance in $\H^n$ is then given by \begin{align} \label{eq:hypDist} d(x,y) &= \text{acosh}\, \left(1+2\frac{\|\eta(x)-\eta(y)\|^2}{(1-\|\eta(x)\|^2)(1-\|\eta(y)\|^2)}\right) \qquad x,y\in\H^n. \end{align} For more details see \cite{Bridson1999} pages 86--90. \paragraph{} The following sections describe the main contributions of the thesis. The next two sections deal with the geometry of the model space $\Theta$ in the context of estimation of conditional models. The later sections study the geometry of the data space $\cX$. %---------------------------------------------------------------- %---------------------------------------------------------------- %---------------------------------------------------------------- %---------------------------------------------------------------- \section{Geometry of Conditional Exponential Models and AdaBoost} \label{sec:boostMl} Several recent papers in statistics and machine learning have been devoted to the relationship between boosting and more standard statistical procedures such as logistic regression. In spite of this activity, an easy-to-understand and clean connection between these different techniques has not emerged. \emcite{Friedman2000} note the similarity between boosting and stepwise logistic regression procedures, and suggest a least-squares alternative, but view the loss functions of the two problems as different, leaving the precise relationship between boosting and maximum likelihood unresolved. \emcite{Kivinen1999} note that boosting is a form of ``entropy projection,'' and \emcite{Lafferty1999} suggests the use of Bregman distances to approximate the exponential loss. \emcite{Mason:99} consider boosting algorithms as functional gradient descent and \emcite{Duffy+Helmbold:00} study various loss functions with respect to the PAC boosting property. More recently, \emcite{Collins2002} show how different Bregman distances precisely account for boosting and logistic regression, and use this framework to give the first convergence proof of AdaBoost. However, in this work the two methods are viewed as minimizing different loss functions. Moreover, the optimization problems are formulated in terms of a reference distribution consisting of the zero vector, rather than the empirical distribution of the data, making the interpretation of this use of Bregman distances problematic from a statistical point of view. \paragraph{} In this section we present a very basic connection between boosting and maximum likelihood for exponential models through a simple convex optimization problem. In this setting, it is seen that the only difference between AdaBoost and maximum likelihood for exponential models, in particular logistic regression, is that the latter requires the model to be normalized to form a probability distribution. The two methods minimize the same $I$-divergence objective function subject to the same feature constraints. Using information geometry, we show that projecting the exponential loss model onto the simplex of conditional probability distributions gives precisely the maximum likelihood exponential model with the specified sufficient statistics. In many cases of practical interest, the resulting models will be identical; in particular, as the number of features increases to fit the training data the two methods will give the same classifiers. We note that throughout the thesis we view boosting as a procedure for minimizing the exponential loss, using either parallel or sequential update algorithms as presented by \emcite{Collins2002}, rather than as a forward stepwise procedure as presented by \emcite{Friedman2000} and \emcite{Freund1996}. \paragraph{} Given the recent interest in these techniques, it is striking that this connection has gone unobserved until now. However in general, there may be many ways of writing the constraints for a convex optimization problem, and many different settings of the Lagrange multipliers that represent identical solutions. The key to the connection we present here lies in the use of a particular non-standard presentation of the constraints. When viewed in this way, there is no need for special-purpose Bregman divergence as in \cite{Collins2002} to give a unified account of boosting and maximum likelihood, and we only make use of the standard $I$-divergence. But our analysis gives more than a formal framework for understanding old algorithms; it also leads to new algorithms for regularizing AdaBoost, which is required when the training data is noisy. In particular, we derive a regularization procedure for AdaBoost that directly corresponds to penalized maximum likelihood using a Gaussian prior. Experiments on UCI data support our theoretical analysis, demonstrate the effectiveness of the new regularization method, and give further insight into the relationship between boosting and maximum likelihood exponential models. \paragraph{} The next section describes an axiomatic characterization of metrics over the space of conditional models and the relationship between the characterized metric and the $I$-divergence. In this sense this section and the next one should be viewed as one unit as they provide an axiomatic characterization of the geometry underlying conditional exponential models such as logistic regression and AdaBoost. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \subsection{Definitions} Let $\cX$ and $\cY$ be finite sets of sizes $k$ and $m$ and $\P_{m-1}^k,\mathbb{R}_{+}^{k\times m}$ be the sets of normalized and non-normalized conditional models as defined in Section~\ref{sec:probSimplexGeom}. Their closures $\overline{\P_{m-1}^k},\overline{\mathbb{R}_{+}^{k\times m}}$ represent the sets of non-negative conditional models. Let \[f=(f_1,\ldots,f_l), \qquad f_j:\cX\times \cY \to \R\] be a sequence functions, which we will refer to as a feature vector. These functions correspond to the weak learners in boosting, and to the sufficient statistics in an exponential model. The empirical distribution associated with a training set $\{(x_i,y_i)\}_{i=1}^n$ is ${\tilde{p}}(x,y) = \frac{1}{n} \sum_{i=1}^{n} \delta_{x_i, x} \, \delta_{y_i, y}.$ Based on ${\tilde{p}}(x,y)$ we define marginal and conditional distributions ${\tilde{p}}(x),{\tilde{p}}(y|x)$ as usual. We assume that ${\tilde{p}}(x)>0$ and that for all $x$ there is a unique $y\in\cY$, denoted by $\tilde{y}(x)$, for which ${\tilde{p}}(y\given x) > 0$. This assumption, referred to as the consistent data assumption, is made to obtain notation that corresponds to the conventions used to present boosting algorithms; it is not essential to the correspondence between AdaBoost and conditional exponential models presented here. We will use the notation $f(x,y)$ to represent the real vector $(f_1(x,y),\ldots,f_l(x,y))$ and $\langle \cdot,\cdot\rangle$ to be the usual scalar or dot product between real vectors. \paragraph{} The conditional exponential model $q(y\given x\,;\theta)$ associated with the feature vector $f$ is defined by \begin{equation} q(y\given x\,;\theta) = \frac{e^{\dotp{\theta}{f(x,y)}}} {\sum_{y} e^{\dotp{\theta}{f(x,y)}}} \qquad \theta\in\R^l \end{equation} where $\langle \cdot,\cdot\rangle$ denotes the standard scalar product between two vectors. The maximum likelihood estimation problem is to determine a parameter vector $\theta$ that maximize the conditional log-likelihood \[\ell(\theta) \deff \sum_{x,y} {\tilde{p}}(x,y) \log q(y\given x\,;\theta).\] \paragraph{} The objective function to be minimized in the multi-label boosting algorithm AdaBoost.M2 \cite{Collins2002} is the exponential loss given by \begin{equation} \cE(\theta) \deff \sum_{i=1}^n \sum_{y\neq y_i} e^{\,\dotp{\theta}{f(x_i, y) - f(x_i, y_i)}}. \end{equation} In the binary case $\cY=\{-1,+1\}$ and taking $f_j(x,y) = \frac{1}{2} y f_j(x)$ the exponential loss becomes \begin{equation} \cE(\theta) = \sum_{i=1}^n e^{-y_i \dotp{\theta}{f(x_i)}}, \label{eq:exploss} \end{equation} the conditional exponential model becomes the logistic model \begin{equation} q(y\given x\,;\theta) \seq \frac{1}{1+e^{-y\dotp{\theta}{f(x)}}}, \end{equation} for which the maximum likelihood problem becomes equivalent to minimizing the logistic loss function \begin{equation} -\ell(\theta)=\sum_{i=1}^n \log\left(1+e^{-y_i \dotp{\theta}{f(x_i)}} \right). \label{eq:logloss} \end{equation} As has been often noted, the log-loss (\ref{eq:logloss}) and the exponential loss (\ref{eq:exploss}) are qualitatively different. The exponential loss (\ref{eq:exploss}) grows exponentially with increasing negative margin $-y \dotp{\theta}{f(x)}$, while the log-loss grows linearly. %------------------------------------------------------------------------------- \subsection{Correspondence Between AdaBoost and Maximum Likelihood} We define the conditional $I$-divergence with respect to a distribution $r$ over $\mathcal{X}$ as \begin{equation} \label{eq:I-div} D_{r}(p, q) \deff \sum_{x} r(x) \sum_y \left( p(y\given x) \log \frac{p(y\given x)}{q(y\given x)} - p(y\given x) + q(y\given x) \right). \end{equation} It is a non-negative measure of discrepancy between two conditional models $p,q\in\overline{\R_{+}^{k\times m}}$. If $q\not\in \R_{+}^{k\times m}$, $D_r(p,q)$ may be $\infty$. The $I$ divergence is not a formal distance function as it does not satisfy symmetry and the triangle inequality. In this section we will always take $r(x)={\tilde{p}}(x)$ and hence we omit it from the notation and write $D(p,q)=D_{{\tilde{p}}}(p,q)$. For normalized conditional models, the $I$-divergence $D(p,q)$ is equal to the Kullback-Leibler divergence \cite{Kullback1968}. The formula presented here \eqref{eq:I-div} is a straightforward adaptation of the non-conditional form of the $I$-divergence studied by \emcite{Csiszar1991}. The $I$-divergence comes up in many applications of statistics and machine learning. See \cite{Kullback1968,OSullivan:98,Amari2000} for many examples of such connections. \paragraph{} We define the feasible set $\cF({\tilde{p}}, f)$ of conditional models associated with $f=(f_1,\ldots,f_l)$ and an empirical distribution ${\tilde{p}}(x,y)$ as \begin{equation} \label{eq:constraints} \cF({\tilde{p}},f) \deff \left\{ p\in\overline{\R_{+}^{k\times m}} : \sum_x {\tilde{p}}(x) \sum_{y} p(y\given x)\, \left(f_j(x,y) - E_{{\tilde{p}}}[f_j\given x]\right) = 0,\; j=1,\ldots,l\right\}. \end{equation} Note that this set is non-empty since ${\tilde{p}}\in\cF({\tilde{p}},f)$ and that under the consistent data assumption $E_{{\tilde{p}}}[f\given x] = f(x,\tilde{y}(x))$. The feasible set represents conditional models that agree with ${\tilde{p}}$ on the conditional expectation of the features. \paragraph{} Consider now the following two convex optimization problems, labeled $P_1$ and $P_2$. \begin{center} \begin{tabular}{cc} \begin{tabular}{rrl} $(P_1)$ & {\it minimize} & $D(p,q_{0})$ \\[2pt] & {\it subject to} & $p \in \cF({\tilde{p}}, f)$ \\ & & \end{tabular} &\hskip20pt \begin{tabular}{rrl} $(P_2)$ & {\it minimize} & $D(p,q_{0})$ \\[2pt] & {\it subject to} & $p \in \cF({\tilde{p}}, f)$ \\ & & $p \in \overline{\P_{m-1}^k}$. \end{tabular} \end{tabular} \end{center} Thus, problem $P_2$ differs from $P_1$ only in that the solution is required to be normalized. As we will show, the dual problem $P_1^*$ corresponds to AdaBoost, and the dual problem $P_2^*$ corresponds to maximum likelihood for exponential models. \paragraph{} This presentation of the constraints is the key to making the correspondence between AdaBoost and maximum likelihood. The constraint $\sum_{x}{\tilde{p}}(x)\sum_y p(y\given x)\, f(x,y) = E_{{\tilde{p}}}[f]$, which is the usual presentation of the constraints for maximum likelihood (as dual to maximum entropy), doesn't make sense for non-normalized models, since the two sides of the equation may not be on the same scale. Note further that attempting to re-scale by dividing by the mass of $p$ to get \[\sum_{x}{\tilde{p}}(x)\frac{\sum_y p(y\given x)\, f(x,y)}{\sum_y p(y\given x)} = E_{{\tilde{p}}}[f]\] would yield nonlinear constraints. \paragraph{} Before we continue, we recall the dual formulation from convex optimization. For mode details refer for example to Section 5 of \cite{Boyd2004}. Given a convex optimization problem \begin{align}\label{eq:primalProb} \min_{x\in\R^n} f_0(x) \quad \text{subject to} \quad h_i(x)=0 \quad i=1,\ldots,r \end{align} the Lagrangian is defined as \begin{align} \label{eq:Lagrangian} \cL(x,\theta)\deff f_0(x)-\sum_{i=1}^r \theta_i h_i(x). \end{align} The vector $\theta\in\mathbb{R}^r$ is called the dual variable or the Lagrange multiplier vector. The Lagrange dual function is defined as $h(\theta)=\inf_x \cL(x,\theta)$ and the dual problem of the original problem \eqref{eq:primalProb}, is $\max_{\theta} h(\theta)$. The dual problem and the original problem, called the primal problem, are equivalent to each other and typically, the easier of the two problems is solved. Both problems are useful, however, as they provide alternative views of the optimization problem. \subsubsection{The Dual Problem $(P_1^*)$} Applying the above definitions to $(P_1)$, and noting that the term $q(y\given x)$ in \eqref{eq:I-div} does not play a role in the minimization problem, the Lagrangian is \begin{align*} \cL_1(p,\theta) &= \sum_x {\tilde{p}}(x) \sum_y p(y\given x) \left(\log \frac{p(y\given x)}{q_{0}(y\given x)} - 1\right) - \sum_i \theta_i \sum_x {\tilde{p}}(x) \sum_y p(y\given x) (f_i(x,y) - E_{{\tilde{p}}}[f_i\given x])\\ &= \sum_x {\tilde{p}}(x) \sum_y p(y\given x) \left(\log \frac{p(y\given x)}{q_{0}(y\given x)} - 1 - \dotp{\theta}{f(x,y) - E_{{\tilde{p}}}[f\given x]}\right). \end{align*} The first step is to minimize the Lagrangian with respect to $p$, which will allow us to express the dual function. Equating the partial derivatives $\frac{\partial \cL_1(p,\theta)}{\partial p(y\given x)}$ to 0 gives \begin{align*} 0&= {\tilde{p}}(x)\left(\log \frac{p(y\given x)}{q_{0}(y\given x)} - 1 - \dotp{\theta}{f(x,y) - E_{{\tilde{p}}}[f\given x]}+p(y\given x)\frac{1}{p(y\given x)}\right)\\ &= {\tilde{p}}(x)\left(\log \frac{p(y\given x)}{q_{0}(y\given x)} - \dotp{\theta}{f(x,y) - E_{{\tilde{p}}}[f\given x]}\right) \end{align*} and we deduce that $z(y\given x\,;\theta) \deff \mbox{arg}\min_{p} \cL_1(p, \theta)$, is \begin{align*} z(y\given x\,;\theta) = q_{0}(y\given x)\, \exp\left(\sum_{j} \theta_j \left(f_j(x,y)-E_{{\tilde{p}}}[f_j\given x]\right)\right). \end{align*} Thus, the dual function $h_1(\theta) = \cL_1(z(y\given x\,;\theta), \theta)$ is given by \begin{align*} h_1(\theta) &= \sum_x {\tilde{p}}(x)\sum_y q_{0}(y\given x) e^{\dotp{\theta}{f(x,y)-E_{{\tilde{p}}}[f\given x]}} \left(\dotp{\theta}{f(x,y)-E_{{\tilde{p}}}[f\given x]}-1-\dotp{\theta}{f(x,y)-E_{{\tilde{p}}}[f\given x]}\right) \\ &=- \sum_x {\tilde{p}}(x) \sum_{y} q_{0}(y\given x) e^{\dotp{\theta}{f(x,y)-E_{{\tilde{p}}}[f\given x]}}. \end{align*} The dual problem $(P_1^*)$ is to determine \begin{align} \theta^\star &= \argmax_\theta h_1(\theta) \nonumber =\argmin_{\theta}\sum_x {\tilde{p}}(x) \sum_{y} q_{0}(y\given x) e^{\dotp{\theta}{f(x,y)-E_{{\tilde{p}}}[f\given x]}}\\ &= \argmin_{\theta}\sum_x {\tilde{p}}(x) \sum_{y} q_{0}(y\given x) e^{\dotp{\theta}{f(x,y)-f(x,\tilde{y}(x)}}\nonumber \\ &= \argmin_{\theta}\sum_x {\tilde{p}}(x) \sum_{y\neq \tilde{y}(x)} q_{0}(y\given x) e^{\dotp{\theta}{f(x,y)-f(x,\tilde{y}(x)}}. \label{eq:dualProb1} \end{align} \subsubsection{The Dual Problem $(P_2^*)$} To derive the dual for $P_2$, we add additional Lagrange multipliers $\mu_x$ for the constraints $\sum_{y} p(y\given x) = 1$ and note that if the normalization constraints are satisfied then the other constraints take the form \begin{align*} \sum_{x} {\tilde{p}}(x) \sum_y p(y\given x) f_j(x,y) = \sum_{x, y} {\tilde{p}}(x, y) f_j(x,y). \end{align*} The Lagrangian becomes \begin{align*} \cL_2(p,\theta,\mu) =& D(p,q_0)-\sum_j \theta_j \left(\sum_{x} {\tilde{p}}(x) \sum_y p(y\given x) f_j(x,y)-\sum_{x, y} {\tilde{p}}(x, y) f_j(x,y)\right)\\ &-\sum_x \mu_x\left(1-\sum_y p(y\given x)\right). \end{align*} Setting the partial derivatives $\frac{\partial \cL_2(p,\theta)}{\partial p(y|x)}$ to 0 and noting that in the normalized case we can ignore the last two terms in the $I$-divergence, we get \begin{align*} 0&= {\tilde{p}}(x)\left(\log \frac{p(y|x)}{q_0(y|x)}+1-\dotp{\theta}{f(x,y)}\right)+ \mu_x \end{align*} from which the minimizer $z(y\given x\,;\theta) \deff \mbox{arg}\min_{p} \cL_2(p, \theta)$ is seen to be \begin{align*} z(y\given x\,;\theta) = q_{0}(y\given x) e^{\dotp{\theta}{f(x,y)}-1-\mu_x/{\tilde{p}}(x)}. \end{align*} Substituting $z(y\given x\,;\theta)$ in $\cL_2$ we obtain the dual function $h_2(\theta,\mu)$. Maximizing the dual function with respect to $\mu$ results in a choice of $\mu_x$ that ensure the normalization of $z$ \[z(y\given x\,;\theta) = \frac{1}{Z_x} q_{0}(y\given x) e^{\dotp{\theta}{f(x,y)}}\] and maximizing $h_2$ with respect to $\theta$ we get the following dual problem \begin{align} \nonumber \theta^*&=\argmax_{\theta} \sum_x {\tilde{p}}(x)\sum_y \frac{1}{Z_x(\theta)} q_{0}(y\given x) e^{\dotp{\theta}{f(x,y)}} (\dotp{\theta}{f(x,y)}- \log Z_x(\theta))\\ \nonumber &-\sum_x {\tilde{p}}(x)\sum_y \dotp{\theta}{f(x,y)} \frac{1}{Z_x(\theta)} q_{0}(y\given x) e^{\dotp{\theta}{f(x,y)}} +\sum_{xy} {\tilde{p}}(x,y)\dotp{\theta}{f(x,y)}-0 \sum_x \mu_x\\ \nonumber &= \argmax_{\theta} -\sum_x {\tilde{p}}(x)\log Z_x(\theta) \sum_y \frac{1}{Z_x(\theta)} q_{0}(y\given x) e^{\dotp{\theta}{f(x,y)}} +\sum_{xy} {\tilde{p}}(x,y)\dotp{\theta}{f(x,y)}\\ \nonumber &=\argmax_{\theta} \sum_{x,y} {\tilde{p}}(x,y) \dotp{\theta}{f(x,y)} -\sum_{x,y} {\tilde{p}}(x,y) \log Z_x(\theta)\\ \nonumber &=\argmax_{\theta}\sum_{x,y} {\tilde{p}}(x,y) \log \frac{1}{Z_x(\theta)} q_{0}(y \given x) e^{\dotp{\theta}{f(x,y)}}\\ &=\argmax_{\theta}\sum_{x} {\tilde{p}}(x) \log \frac{1}{Z_x(\theta)} q_{0}(\tilde y(x) \given x) e^{\dotp{\theta}{f(x,\tilde y(x))}} \label{eq:dualProb2}. \end{align} \subsubsection{Special cases} It is now straightforward to derive various boosting and conditional exponential models problems as special cases of the dual problems \eqref{eq:dualProb1} and \eqref{eq:dualProb2}. \vspace{0.1in} \newline \emph{Case 1: AdaBoost.M2}.\enspace The dual problem $(P_1^*)$ with $q_{0}(y\given x) = 1$ is the optimization problem of AdaBoost.M2 \begin{align*} \theta^\star &= \argmin_{\theta} \sum_{x}{\tilde{p}}(x) \sum_{y \neq y_i} e^{\dotp{\theta}{f(x_i,y)-f(x_i, y_i)}}=\argmin_{\theta} \cE(\theta). \end{align*} \vspace{0.1in} \newline \emph{Case 2: Binary AdaBoost}.\enspace The dual problem $(P_1^*)$ with $q_{0}(y\given x) = 1$, $\mathcal{Y}=\{-1,1\}$ and $f_j(x,y) = \frac{1}{2} y\,f_j(x)$ is the optimization problem of binary AdaBoost \begin{align*} \theta^\star &= \argmin_{\theta} \sum_{x}{\tilde{p}}(x) e^{-y_i \dotp{\theta}{f(x_i)}}. \end{align*} \vspace{0.1in} \newline \emph{Case 3: Maximum Likelihood for Exponential Models}.\enspace The dual problem $(P_2^*)$ with $q_{0}(y\given x) = 1$ is is maximum (conditional) likelihood for a conditional exponential model with sufficient statistics $f_j(x,y)$ \begin{align*} \theta^\star &= \argmax_{\theta}\sum_{x} {\tilde{p}}(x) \log \frac{1}{Z_x}e^{\dotp{\theta}{f(x,\tilde y(x))}} = \argmax_{\theta} \ell(\theta). \end{align*} \vspace{0.1in} \newline \emph{Case 4: Logistic Regression}.\enspace The dual problem $(P_2^*)$ with $q_{0}(y\given x) = 1$, $\mathcal{Y}=\{-1,1\}$ and $f_j(x,y) = \frac{1}{2} y\,f_j(x)$ is maximum (conditional) likelihood for binary logistic regression. \begin{align*} \theta^\star = \argmax_{\theta}\sum_{x} {\tilde{p}}(x) \frac{1}{1 + e^{-{\tilde{y}(x)\dotp{\theta}{f(x)}}}}. \end{align*} We note that it is not necessary to scale the features by a constant factor here, as in \cite{Friedman2000}; the correspondence between logistic regression and boosting is direct. \vspace{0.1in} \newline \emph{Case 5: Exponential Models with Carrier Density}.\enspace Taking $q_0(y|x)\neq 1$ to be a non-parametric density estimator in $(P_2^*)$ results in maximum likelihood for exponential models with a carrier density $q_0$. Such semi-parametric models have been proposed by \emcite{Efron1996} for integrating between parametric and nonparametric statistical modeling and by \emcite{DellaPietra1992} and \emcite{Rosenfeld1996} for integrating exponential models and $n$-gram estimators in language modeling. \paragraph{} Making the Lagrangian duality argument rigorous requires care, because of the possibility that the solution may lie on the (topological) boundary of $\overline{\mathbb{R}_{+}^{k\times m}}$ or $\overline{\P_{m-1}^k}$. \paragraph{} Let $\cQ_1$ and $\cQ_2$ be \begin{align*} \cQ_1(q_{0}, f) &= \left\{ q\in\overline{\mathbb{R}_{+}^{k\times m}} \;\given\; q(y\given x) \,=\, q_{0}(y\given x)\, e^{\dotp{\theta}{f(x,y)-f(x,\tilde{y}(x))}} ,\,\theta\in\R^l\right\} \\ \cQ_2(q_{0}, f) &= \left\{ q\in\overline{\P_{m-1}^k} \;\given\; q(y\given x) \,\propto\, {q_{0}(y\given x)\, e^{\dotp{\theta}{f(x,y)}}} ,\,\theta\in\R^l\right\} \end{align*} and the boosting solution $q^\star_{\mbox{\small\it boost}}$ and maximum likelihood solution $q^\star_{\mbox{\small\it ml}}$ be \begin{align*} q^\star_{\mbox{\small\it boost}} &= \argmin_{q\in \overline{\cQ_1}} \sum_{x}{\tilde{p}}(x) \sum_y q(y\given x) \\ q^\star_{\mbox{\small\it ml}} &= \argmax_{q\in \overline{\cQ_2}} \sum_{x}{\tilde{p}}(x) \log q(\tilde{y}(x)\given x). \end{align*} \paragraph{} The following proposition corresponds to Proposition~4 of \cite{DellaPietra1997} for the usual Kullback-Leibler divergence; the proof for the $I$-divergence carries over with only minor changes. In \cite{DellaPietra2001} the duality theorem is proved for a class of Bregman divergences, including the extended Kullback-Leibler divergence as a special case. Note that we do not require divergences such as $D(0, q)$ as in \cite{Collins2002}, but rather $D({\tilde{p}}, q)$, which is more natural and interpretable from a statistical point-of-view. \begin{prop} Suppose that $D({\tilde{p}},q_{0}) < \infty$. Then $q^\star_{\mbox{\small\it boost}}$ and $q^\star_{\mbox{\small\it ml}}$ exist, are unique, and satisfy \begin{eqnarray} q^\star_{\mbox{\small\it boost}} &=& \argmin_{p\in\cF} D(p,q_{0}) \seq \argmin_{q\in\overline{\cQ_1}} D({\tilde{p}}, q) \\ q^\star_{\mbox{\small\it ml}} &=& \argmin_{p\in\cF\cap\P_{m-1}^k} D(p,q_{0}) \seq \argmin_{q\in\overline{\cQ_2}} D({\tilde{p}}, q) \end{eqnarray} Moreover, $q^\star_{\mbox{\small\it ml}}$ is computed in terms of $q^\star_{\mbox{\small\it boost}}$ as $ q^\star_{\mbox{\small\it ml}} = \displaystyle \argmin_{p\in\cF\cap\P_{m-1}^k} D(p,q^\star_{\mbox{\small\it boost}}) $. \end{prop} \paragraph{} This result has a simple geometric interpretation. The non-normalized exponential family $\cQ_1$ intersects the feasible set of measures $\cF$ satisfying the constraints (\ref{eq:constraints}) at a single point. The algorithms presented in \cite{Collins2002} determine this point, which is the exponential loss solution $q^\star_{\mbox{\small\it boost}} = \argmin_{q\in\overline{\cQ_1}} D({\tilde{p}}, q)$ (see Figure~\ref{fig:geom}, left). On the other hand, maximum conditional likelihood estimation for an exponential model with the same features is equivalent to the problem $q^\star_{\mbox{\small\it ml}} = \argmin_{q\in\overline{\cQ_1'}} D({\tilde{p}}, q)$ where $\cQ_1'$ is the exponential family with additional Lagrange multipliers, one for each normalization constraint. The feasible set for this problem is $\cF \cap \overline{\P_{m-1}^k}$. Since $\cF \cap\overline{\P_{m-1}^k} \subset \cF$, by the Pythagorean equality we have that $q^\star_{\mbox{\small\it ml}} = \argmin_{p\in\cF\cap\overline{\P_{m-1}^k}} D(p, q^\star_{\mbox{\small\it boost}})$ (see Figure \ref{fig:geom}, right). \begin{figure} \centering \psfrag{f}{$\cF$} \psfrag{qb}{$q^\star_{\mbox{\small\it boost}}$} \psfrag{qm}{$q^\star_{\mbox{\small\it ml}}$} \psfrag{g}{$\mathcal{G}$} {\psfrag{q}{$\mathcal{Q}_1$} \includegraphics[scale=0.45]{boostingLogRegFigures/geom1.eps}} \hskip10pt \psfrag{f}{$\cF$} \psfrag{qb}{$q^\star_{\mbox{\small\it boost}}$} \psfrag{qm}{$q^\star_{\mbox{\small\it ml}}$} \psfrag{g}{$\cF\cap \overline{\P_{m-1}^k}$} {\psfrag{q}{$\mathcal{Q}_1'$} \includegraphics[scale=0.45]{boostingLogRegFigures/geom2.eps}} \caption{\small Geometric view of duality. Minimizing the exponential loss finds the member of $\cQ_1$ that intersects the feasible set of measures satisfying the moment constraints (left). When we impose the additional constraint that each conditional distribution $q_\theta(y\given x)$ must be normalized, we introduce a Lagrange multiplier for each training example $x$, giving a higher-dimensional family $\cQ_1'$. By the duality theorem, projecting the exponential loss solution onto the intersection of the feasible set with the simplex of conditional probabilities, $\cF\cap\overline{\P_{m-1}^k}$, we obtain the maximum likelihood solution. In many practical cases this projection is obtained by simply normalizing by a constant, resulting in an identical model.} \label{fig:geom} \vskip.4in \end{figure} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \subsection{Regularization} Minimizing the exponential loss or the log-loss on real data often fails to produce finite parameters. Specifically, this happens when for some feature $f_j$ \begin{eqnarray} \label{eq:cond0} && f_j(x,y) - f_j(x,\tilde{y}(x)) \geq 0 \; \mbox{for all $y$ and $x$ with ${\tilde{p}}(x)>0$} \\ \mbox{or} && f_j(x,y) - f_j(x,\tilde{y}(x)) \leq 0 \; \mbox{for all $y$ and $x$ with ${\tilde{p}}(x)>0$} \nonumber \end{eqnarray} This is especially harmful since often the features for which (\ref{eq:cond0}) holds are the most important for the purpose of discrimination. The parallel update in \cite{Collins2002} breaks down in such cases, resulting in parameters going to $\infty$ or $-\infty$. On the other hand, iterative scaling algorithms work in principle for such features. In practice however, either the parameters $\theta$ need to be artificially capped or the features need to be thrown out altogether, resulting in a partial and less discriminating set of features. Of course, even when (\ref{eq:cond0}) does not hold, models trained by maximizing likelihood or minimizing exponential loss can overfit the training data. The standard regularization technique in the case of maximum likelihood employs parameter priors in a Bayesian framework. \paragraph{} In terms of convex duality, a parameter prior for the dual problem corresponds to a ``potential'' on the constraint values in the primal problem. The case of a Gaussian prior on $\theta$, for example, corresponds to a quadratic potential on the constraint values in the primal problem. Using this correspondence, the connection between boosting and maximum likelihood presented in the previous section indicates how to regularize AdaBoost using Bayesian MAP estimation for non-normalized models, as explained below. \paragraph{} We now consider primal problems over $(p,c)$ where $p\in\R_{+}^{k\times m}$ and $c\in\R^m$ is a parameter vector that relaxes the original constraints. The feasible set $\cF({\tilde{p}}, f, c)\subset \R_{+}^{k\times m}$ allows a violation of the expectation constraints, represented by the vector $c$ \begin{equation} \label{eq:constraintsreg} \cF({\tilde{p}},f,c) \seq \left\{ p\in\R_{+}^{k\times m} \;\given\; \sum_x {\tilde{p}}(x) \sum_{y} p(y\given x)\, \left(f_j(x,y) - E_{{\tilde{p}}}[f_j\given x]\right) = c_j\right\}. \end{equation} A regularized problem for non-normalized models is defined by \begin{center} \begin{tabular}{rrl} $(P_{1,\mbox{\small reg}})$ & {\it minimize} & $D(p,q_{0}) + U(c)$ \\[2pt] & {\it subject to} & $p \in \cF({\tilde{p}}, f, c)$ \\ \end{tabular} \end{center} where $U:\mathbb{R}^l\to \mathbb{R}$ is a convex function whose minimum is at $0\in\R^l$. Intuitively $(P_{1,\mbox{\small reg}})$ allows some trade-off between achieving low $I$ divergence to $q_0$ and some constraint violation, with the exact form of the trade-off represented by the function $U$. Note that it is possible to choose $U$ in a way that considers some feature constraints more important than others. This may be useful when the values of the features $(f_1,\ldots,f_l)$ are known to be corrupted by noise, where the noise intensity varies among the features. \paragraph{} The dual function of the regularized problem $(P_{1,\text{\small reg}})$, as derived in Appendix~\ref{sec:regulLossFunc}, is \[h_{1,\mbox{\small reg}}(\theta) = h_{1}(\theta) + U^*(\theta)\] where $U^*(\theta)$ is the convex conjugate of $U$. For a quadratic penalty $U(c)=\sum_j \frac{1}{2}\sigma_j^2\,c_j^2 $, we have ${U}^*(\theta)=-\sum_j\frac{1}{2}\sigma_j^{-2}\,\theta_j^{2}$ and the dual function becomes \begin{eqnarray} \label{eq:dual} h_{1,\mbox{\small reg}}(\theta) \seq - \sum_x {\tilde{p}}(x) \sum_{y} q_{0}(y\given x)\, e^{\dotp{\theta}{f(x,y)-f(x,\tilde{y}(x))}} - \sum_j\frac{\theta_j^2}{2\sigma_j^2}. \end{eqnarray} A sequential update rule for (\ref{eq:dual}) incurs the small additional cost of solving a nonlinear equation by Newton's method every iteration. See Appendix~\ref{sec:regulLossFunc} for more details. \emcite{Chen2000} contains a similar regularization for maximum likelihood for exponential models in the context of statistical language modeling. \subsection{Experiments} We performed experiments on some of the UCI datasets \cite{UCI} in order to investigate the relationship between boosting and maximum likelihood empirically. The weak learner was \texttt{FindAttrTest} as described in \cite{Freund1996}, and the training set consisted of a randomly chosen 90\% of the data. Table~\ref{table:reg} shows experiments with regularized boosting. Two boosting models are compared. The first model $q_1$ was trained for 10 features generated by \texttt{FindAttrTest}, excluding features satisfying condition (\ref{eq:cond0}). Training was carried out using the parallel update method described in \cite{Collins2002}. The second model, $q_2$, was trained using the exponential loss with quadratic regularization. The performance was measured using the conditional log-likelihood of the (normalized) models over the training and test set, denoted $\ell_{\ss{train}}$ and $\ell_{\ss{test}}$, as well as using the test error rate $\epsilon_{\ss{test}}$. The table entries were averaged by 10-fold cross validation. \paragraph{} For the weak learner \texttt{FindAttrTest}, only the Iris dataset produced features that satisfy (\ref{eq:cond0}). On average, 4 out of the 10 features were removed. As the flexibility of the weak learner is increased, (\ref{eq:cond0}) is expected to hold more often. On this dataset regularization improves both the test set log-likelihood and error rate considerably. In datasets where $q_1$ shows significant over-fitting, regularization improves both the log-likelihood measure and the error rate. In cases of little over-fitting (according to the log-likelihood measure), regularization only improves the test set log-likelihood at the expense of the training set log-likelihood, however without affecting much the test set error. \begin{table} \centering \begin{tabular}{|l|c|c|c||c|c|c|} \multicolumn{1}{c}{} & \multicolumn{3}{c}{\it Unregularized} & \multicolumn{3}{c}{\it Regularized} \\ \hline Data & $\ell_{\ss{train}}(q_1)$ & $\ell_{\ss{test}}(q_1)$&$\epsilon_{\ss{test}}(q_1)$ & $\ell_{\ss{train}}(q_2)$ & $\ell_{\ss{test}}(q_2)$ & $\epsilon_{\ss{test}}(q_2)$ \\ \hline\hline Promoters& -0.29 &-0.60 &0.28 &-0.32 &-0.50 &0.26\\ \hline Iris&-0.29 &-1.16 &0.21 &-0.10 &-0.20&0.09\\ \hline Sonar& -0.22 &-0.58 &0.25 &-0.26 &-0.48&0.19\\ \hline Glass&-0.82 &-0.90&0.36& -0.84 &-0.90 &0.36 \\ \hline Ionosphere&-0.18 &-0.36 &0.13 &-0.21 &-0.28 &0.10 \\ \hline Hepatitis&-0.28 &-0.42 &0.19 &-0.28 &-0.39&0.19 \\ \hline Breast Cancer Wisconsin&-0.12 &-0.14 &0.04 &-0.12 &-0.14 &0.04 \\ \hline Pima-Indians& -0.48&-0.53 &0.26 &-0.48 &-0.52 &0.25 \\ \hline \end{tabular} \caption{\small Comparison of unregularized to regularized boosting. For both the regularized and unregularized cases, the first column shows training log-likelihood, the second column shows test log-likelihood, and the third column shows test error rate. Regularization reduces error rate in some cases while it consistently improves the test set log-likelihood measure on all datasets. All entries were averaged using 10-fold cross validation.} \label{table:reg} \end{table} \begin{figure} \begin{center} \vbox to 4.0in{ \def\yleft{\lower.4ex\hbox{\tiny{$\ell_{\st{train}}(q^{\star}_{\st{boost}})$}}} \def\yright{\lower.3ex\hbox{\tiny{$D_{\st{train}}(q^{\star}_{\st{ml}},q^{\star}_{\st{boost}})$}}} \def\xleft{\raise.8ex\hbox{\tiny{$\ell_{\st{train}}(q^{\star}_{\st{ml}})$}}} \def\xright{\lower1ex\hbox{\tiny{$\ell_{\st{train}}(q^{\star}_{\st{ml}})$}}} \begin{tabular}{ccl} \hskip.2in {\psfrag{xl}{\xleft} \psfrag{yl}{\yleft} \includegraphics[scale=0.35]{boostingLogRegFigures/hepatitisD.eps} } &&\hskip0.3in\mbox{\ } {\psfrag{xl}{\xright} \psfrag{yl}{\yright} \includegraphics[scale=0.35]{boostingLogRegFigures/hepatitis2D.eps} } \\[10pt] \hskip.2in {\psfrag{xl}{\xleft} \psfrag{yl}{\yleft} \includegraphics[scale=0.35]{boostingLogRegFigures/promotersD.eps} } &&\hskip0.3in\mbox{\ } {\psfrag{xl}{\xright} \psfrag{yl}{\yright} \includegraphics[scale=0.35]{boostingLogRegFigures/promoters2D.eps} } \\[10pt] \hskip.2in {\psfrag{xl}{\xleft} \psfrag{yl}{\yleft} \includegraphics[scale=0.35]{boostingLogRegFigures/sonarD.eps} } &&\hskip0.3in\mbox{\ } {\psfrag{xl}{\xright} \psfrag{yl}{\yright} \includegraphics[scale=0.35]{boostingLogRegFigures/sonar2D.eps} } \\[10pt] \hskip.2in {\psfrag{xl}{\xleft} \psfrag{yl}{\yleft} \includegraphics[scale=0.35]{boostingLogRegFigures/glassD.eps} } && \hskip0.3in\mbox{\ } {\psfrag{xl}{\xright} \psfrag{yl}{\yright} \includegraphics[scale=0.35]{boostingLogRegFigures/glass2D.eps} } \end{tabular} } \end{center} \vskip5pt \caption{\small Comparison of AdaBoost and maximum likelihood on four UCI datasets: Hepatitis (top row), Promoters (second row), Sonar (third row) and Glass (bottom row). The left column compares $\ell_{\ss{train}}(q^{\star}_{\ss{ml}})$ to $\ell_{\ss{train}}(q^{\star}_{\ss{boost}})$, and the right column compares $\ell_{\ss{train}}(q^{\star}_{\ss{ml}})$ to $D_{\ss{train}}(q^{\star}_{\ss{ml}},q^{\star}_{\ss{boost}})$.} \label{fig:ll} \end{figure} \begin{figure} \begin{center} \vbox to 3.8in { \def\yleft{\raise1ex\hbox{\tiny{$\ell_{\st{test}}(q^{\star}_{\st{boost}})$}}} \def\xright{\lower.4ex\hbox{\tiny{$\epsilon_{\st{test}}(q^{\star}_{\st{ml}})$}}} \def\xleft{\raise.2ex\hbox{\tiny{$\ell_{\st{test}}(q^{\star}_{\st{ml}})$}}} \def\yright{\raise1ex\hbox{\tiny{$\epsilon_{\st{test}}(q^{\star}_{\st{boost}})$}}} \begin{tabular}{ccl} {\psfrag{xl}{\xleft} \psfrag{yl}{\yleft} \includegraphics[scale=0.35]{boostingLogRegFigures/hepatitis3D.eps} } && {\psfrag{xl}{\xright} \psfrag{yl}{\yright} \includegraphics[scale=0.35]{boostingLogRegFigures/hepatitis4D.eps} } \\[4pt] {\psfrag{xl}{\xleft} \psfrag{yl}{\yleft} \includegraphics[scale=0.35]{boostingLogRegFigures/promoters3D.eps} } && {\psfrag{xl}{\xright} \psfrag{yl}{\yright} \includegraphics[scale=0.35]{boostingLogRegFigures/promoters4D.eps} } \\[4pt] {\psfrag{xl}{\xleft} \psfrag{yl}{\yleft} \includegraphics[scale=0.35]{boostingLogRegFigures/sonar3D.eps} } && {\psfrag{xl}{\xright} \psfrag{yl}{\yright} \includegraphics[scale=0.35]{boostingLogRegFigures/sonar4D.eps} } \\[4pt] {\psfrag{xl}{\xleft} \psfrag{yl}{\yleft} \includegraphics[scale=0.35]{boostingLogRegFigures/glass3D.eps} } && {\psfrag{xl}{\xright} \psfrag{yl}{\yright} \includegraphics[scale=0.35]{boostingLogRegFigures/glass4D.eps} } \end{tabular} } \end{center} \vskip20pt \caption{\small Comparison of AdaBoost and maximum likelihood on the same UCI datasets as in the previous figure. The left column compares the test likelihoods, $\ell_{\ss{test}}(q^{\star}_{\ss{ml}})$ to $\ell_{\ss{test}}(q^{\star}_{\ss{boost}})$, and the right column compares test error rates, $\epsilon_{\ss{test}}(q^{\star}_{\ss{ml}})$ to $\epsilon_{\ss{test}}(q^{\star}_{\ss{boost}})$. In each plot, the color represents the {\it training\/} likelihood $\ell_{\ss{train}}(q^{\star}_{\ss{ml}})$; red corresponds to fitting the training data well.} \label{fig:color} \end{figure} \ignore{ \begin{figure} \begin{center} \begin{tabular}{ccl} {\psfrag{yl}{\tiny{$\ell_{\st{test}}(q^{\star}_{\st{ml}})$}} \psfrag{xl}{\tiny{$\ell_{\st{test}}(q^{\star}_{\st{boost}})$}} \includegraphics[scale=0.35]{boostingLogRegFigures/ionoshpere3D.eps} } && {\psfrag{xl}{\tiny{$\epsilon_{\st{test}}(q^{\star}_{\st{ml}})$}} \psfrag{yl}{\tiny{$\epsilon_{\st{test}}(q^{\star}_{\st{boost}})$}} \includegraphics[scale=0.35]{boostingLogRegFigures/ionoshpere4D.eps} } \\[3pt] {\psfrag{yl}{\tiny{$\ell_{\st{test}}(q^{\star}_{\st{ml}})$}} \psfrag{xl}{\tiny{$\ell_{\st{test}}(q^{\star}_{\st{boost}})$}} \includegraphics[scale=0.35]{boostingLogRegFigures/pima3D.eps} } && {\psfrag{xl}{\tiny{$\epsilon_{\st{test}}(q^{\star}_{\st{ml}})$}} \psfrag{yl}{\tiny{$\epsilon_{\st{test}}(q^{\star}_{\st{boost}})$}} \includegraphics[scale=0.35]{boostingLogRegFigures/pima4D.eps} } \end{tabular} \end{center} \caption{\small Comparison of AdaBoost and maximum likelihood on two UCI datasets: Ionosphere (top row) and Pima-Indians (second row). Each row has two plots: $\ell_{\ss{test}}(q^{\star}_{\ss{ml}})$, $\ell_{\ss{test}}(q^{\star}_{\ss{boost}})$ (left) and test error rates $\epsilon_{\ss{test}}(q^{\star}_{\ss{ml}})$, $\epsilon_{\ss{test}}(q^{\star}_{\ss{boost}})$ (right). In all figures the color represent $\ell_{\ss{train}}(q^{\star}_{\ss{ml}})$.} \label{fig:color2} \end{figure} } \paragraph{} Next we performed a set of experiments to test how much $q_{\ss{boost}}^{\star}$ differs from $q_{\ss{ml}}^{\star}\,$, where the boosting model is normalized (after training) to form a conditional probability distribution. For different experiments, \texttt{FindAttrTest} generated a different number of features (10--100), and the training set was selected randomly. The plots in Figure~\ref{fig:ll} show for different datasets the relationship between $\ell_{\ss{train}}(q_{\ss{ml}}^{\star})$ and $\ell_{\ss{train}}(q_{\ss{boost}}^{\star})$ as well as between $\ell_{\ss{train}}(q_{\ss{ml}}^{\star})$ and $D_{\ss{train}}(q^{\star}_{\ss{ml}},q^{\star}_{\ss{boost}})$. The trend is the same in each data set: as the number of features increases so that the training data is more closely fit ($\ell_{\ss{train}}(q_{\ss{ml}}) \rightarrow 0$), the boosting and maximum likelihood models become more similar, as measured by the $I$-divergence. \paragraph{} The plots in Figure~\ref{fig:color} show the relationship between the test set log-likelihoods, $\ell_{\ss{test}}(q_{\ss{ml}}^{\star})$ to $\ell_{\ss{test}}(q_{\ss{boost}}^{\star})$, together with the test set error rates $\epsilon_{\ss{test}}(q_{\ss{ml}}^{\star})$ and $\epsilon_{\ss{test}}(q_{\ss{boost}}^{\star})$. In these figures the testing set was chosen to be 50\% of the total data. The color represents the {\it training\/} data log-likelihood, $\ell_{\ss{train}}(q_{\ss{ml}}^{\star})$, with the color red corresponding to high likelihood. In order to indicate the number of points at each error rate, each circle was shifted by a small random value to avoid points falling on top of each other. \paragraph{} While the plots in Figure~\ref{fig:ll} indicate that $\ell_{\ss{train}}(q^{\star}_{\ss{ml}}) > \ell_{\ss{train}}(q^{\star}_{\ss{boost}})$, as expected, on the test data the linear trend is reversed, so that $\ell_{\ss{test}}(q^{\star}_{\ss{ml}}) < \ell_{\ss{test}}(q^{\star}_{\ss{boost}})$. This is a result of the fact that for the above data-sets and features, little over-fitting occurs and the more aggressive exponential loss criterion is superior to the more relaxed log-loss criterion. However, as $\ell(q^{\star}_{\ss{ml}})\longrightarrow 0$, the two models come to agree. Appendix~\ref{sec:divExpModels} shows that for any exponential model $q_{\theta}\in\mathcal{Q}_2,$ \begin{equation} \label{eq:KLeqll} D_{\ss{train}}(q^{\star}_{\ss{ml}},q_{\theta}) \;\;=\;\; \ell(q^{\star}_{\ss{ml}})-\ell(q_{\theta}). \end{equation} By taking $q_{\theta}=q^{\star}_{\ss{boost}}$ it is seen that as the difference between $\ell(q^{\star}_{\ss{ml}})$ and $\ell(q^{\star}_{\ss{boost}})$ gets smaller, the divergence between the two models also gets smaller. %Furthermore, since the correlation %coefficient $\rho$ is close to 1, we can use the approximation %$\ell(q^{\star}_{\ss{boost}})\approx \alpha\, \ell(q^{\star}_{\ss{ml}})$ %to obtain %\begin{equation} %D_{\ss{train}}(q^{\star}_{\ss{ml}},q^{\star}_{\ss{boost}})\;\;\approx\;\; %(1-\alpha)\,\ell(q^{\star}_{\ss{ml}}). %\end{equation} \paragraph{} The results are consistent with the theoretical analysis. As the number of features is increased so that the training data is fit more closely, the model matches the empirical distribution ${\tilde{p}}$ and the normalizing term $Z(x)$ becomes a constant. In this case, normalizing the boosting model $q^{\star}_{\ss{boost}}$ does not violate the constraints, and results in the maximum likelihood model. \paragraph{} In Appendix \ref{sec:parUpdate},\ref{sec:secUpdate} we derive update rules for exponential loss minimization. These update rules are derived by minimizing an auxiliary function that bounds from above the reduction in loss. See \cite{Collins2002} for the definition of an auxiliary function and proofs that these functions are indeed auxiliary functions. The derived update rules are similar to the ones derived by \emcite{Collins2002}, except that we do not assume that $M=\mbox{max}_{i,y} \sum_j|f_j(x,y)-f_j(x,\tilde{y})|<1$. In Appendix~\ref{sec:regulLossFunc} the regularized formulation is shown in detail and a sequential update rule is derived. Appendix~\ref{sec:divExpModels} contains a proof for (\ref{eq:KLeqll}). \paragraph{} The next section derives an axiomatic characterization of the geometry of conditional models. It then builds on the results of this section to given an axiomatic characterization of the geometry underlying conditional exponential models and AdaBoost. %---------------------------------------------------------------- %---------------------------------------------------------------- %---------------------------------------------------------------- %---------------------------------------------------------------- \section{Axiomatic Geometry for Conditional Models} \label{sec:axiomaticGeometry} A fundamental assumption in the information geometric framework, is the choice of the Fisher information as the metric that underlies the geometry of probability distributions. The choice of the Fisher information metric may be motivated in several ways the strongest of which is \v{C}encov's characterization theorem (Lemma 11.3, \cite{Chentsov1982}). In his theorem, \v{C}encov proves that the Fisher information metric is the only metric that is invariant under a family of probabilistically meaningful mappings termed congruent embeddings by a Markov morphism. Later on, Campbell extended \v{C}encov's result to include non-normalized positive models \cite{Campbell86}. \paragraph{} The theorems of \v{C}encov and Campbell are particularly interesting since the Fisher information is pervasive in statistics and machine learning. It is the asymptotic variance of the maximum likelihood estimators under some regularity conditions. Cramer and Rao used it to compute a lower bound on the variance of arbitrary unbiased estimators. In Bayesian statistics, it was used by Jeffreys to define non-informative prior. It is tightly connected to the Kullback-Leibler divergence which the cornerstone of maximum likelihood estimation for exponential models as well as various aspects of information theory. \paragraph{} While the geometric approach to statistical inference has attracted considerable attention, little research was conducted on the geometric approach to conditional inference. The characterization theorems of \v{C}encov and Campbell no longer apply in this setting and the different ways of choosing a geometry for the space of conditional distributions, in contrast to the non-conditional case, are not supported by theoretical considerations. \paragraph{} In this section we extend the results of \v{C}encov and Campbell to provide an axiomatic characterization of conditional information geometry. We derive the characterization theorem in the setting of non-normalized conditional models from which the geometry for normalized models is obtained as a special case. In addition, we demonstrate a close connection between the characterized geometry and the conditional $I$-divergence which leads to a new axiomatic interpretation of the geometry underlying the primal problems of logistic regression and AdaBoost. This interpretation builds on the connection between AdaBoost and constrained minimization of $I$-divergence described in Section~\ref{sec:boostMl}. \paragraph{} Throughout the section we consider spaces of strictly positive conditional models where the sample spaces of the explanatory and response variable are finite. Moving to the infinite case presents some serious difficulties. The positivity constraint on the other hand does not play a crucial role and may by discarded at some notational cost. \paragraph{} In the characterization theorem we will make use of the fact that $\P_{m-1}^k\cap \mathbb{Q}^{k\times m}$ and $\mathbb{R}_{+}^{k\times m}\cap \mathbb{Q}^{k\times m}= \mathbb{Q}_{+}^{k\times m}$ are dense in $\P_{m-1}^k$ and $\mathbb{R}_{+}^{k\times m}$ respectively. Since continuous functions are uniquely characterized by their values on dense sets, it is enough to compute the metric for positive rational models $\mathbb{Q}_{+}^{k\times m}$. The value of the metric on non-rational models follows from its continuous extension to $\mathbb{R}_{+}^{k\times m}$. \paragraph{} In Section~\ref{sec:congEmbed} we define a class of transformations called congruent embeddings by a Markov morphism. These transformations set the stage for the axioms in the characterization theorem of Section~\ref{sec:charTheorem}. %----------------------------------------------------------- \subsection{Congruent Embeddings by Markov Morphisms of Conditional Models} \label{sec:congEmbed} The characterization result of Section~\ref{sec:charTheorem} is based on axioms that require geometric invariance through a set of transformations between conditional models. These transformations are a generalization of the transformations underlying \v{C}encov's theorem. For consistency with the terminology of \emcite{Chentsov1982} and \emcite{Campbell86} we refer to these transformations as Congruent embeddings by Markov morphisms of conditional models. \begin{defn} Let $\mathcal{A}=\{A_1,\ldots,A_{m}\}$ be a set partition of $\{1,\ldots,n\}$ with $00 & j\in A_i\\ 0 & j\not\in A_i \end{cases}. \end{align*} \end{defn} In other words, $\mathcal{A}$-stochastic matrices are stochastic matrices whose rows are concentrated on the sets of the partition $\mathcal{A}$. For example, if $\mathcal{A}=\{\{1,3\},\{2,4\},\{5\}\}$ then the following matrix is $\mathcal{A}$-stochastic \begin{align} \label{eq:matrixExample} \begin{pmatrix} 1/3&0 &2/3&0 &0\\ 0 &1/2&0 &1/2&0\\ 0 & 0 &0 &0 &1 \end{pmatrix}. \end{align} Obviously, the columns of any $\mathcal{A}$-stochastic matrix have precisely one non-zero element. If $m=n$ then an $\mathcal{A}$-stochastic matrix is a permutation matrix. \paragraph{} Multiplying a row probability vector $p\in\mathbb{R}_{+}^{1\times m}$ with an $\mathcal{A}$-stochastic matrix $Q\in\mathbb{R}^{m\times n}$ results in a row probability vector $q\in\mathbb{R}_{+}^{1\times n}$. The mapping $p\mapsto pQ$ has the following statistical interpretation. The event $x_i$ is split into $|A_i|$ distinct events stochastically, with the splitting probabilities given by the $i$-row of $Q$. The new event space, denoted by $\mathcal{Z}=\{z_1,\ldots,z_n\}$ may be considered a refinement of $\mathcal{X}=\{x_1,\ldots,x_m\}$ (if $m 0$. \end{enumerate} If in addition $M$ is compact, then $K_t$ can be expressed in terms of the eigenvalues and eigenfunctions of the Laplacian as $K_t(x,y) = \sum_{i=0}^{\infty} e^{-\lambda_i t} \phi_i(x)\, \phi_i(y)$. \end{thrm} Properties 2 and 3 imply that $K_t(x,y)$ solves the heat equation in $x$, starting from $y$. It follows that $e^{t\Delta} f(x) = f(x,t) = \int_{M} K_t(x,y)\,f(y)\,\dd dy$ solves the heat equation with initial conditions $f(x,0) = f(x)$, since \begin{eqnarray} \frac{\partial f(x,t)}{\partial t} &=& \int_M \frac{\partial K_t(x,y)}{\partial t}\, f(y)\,\dd dy \\ &=& \int_M \Delta K_t(x,y)\, f(y)\,\dd dy \\ &=& \Delta \int_M K_t(x,y)\, f(y)\, \dd dy \\ &=& \Delta f(x) \end{eqnarray} and $\lim_{t\rightarrow 0} f(x,t) = \int_M \lim_{t\rightarrow 0} K_t(x,y)\, dy = f(x)$. Property~4 implies that $e^{t\Delta}\kern.1ex e^{s\Delta} = e^{(t+s)\Delta}$, which has the physically intuitive interpretation that heat diffusion for time $t$ is the composition of heat diffusion up to time~$s$ with heat diffusion for an additional time $t-s$. Since $e^{t\Delta}$ is a positive operator, \begin{eqnarray} \int_M\int_M K_t(x,y) f(x) f(y) \, \dd dx\, \dd dy &=& % \int_M f(x) \,e^{t\Delta} f(x) \, \dd dx \\ \dotp{f}{e^{t\Delta} f} \,\geq\, 0 \end{eqnarray} and $K_t(x,y)$ is positive-definite. In the compact case, positive-definiteness follows directly from the expansion $K_t(x,y) = \sum_{i=0}^{\infty} e^{-\lambda_i t} \phi_i(x)\, \phi_i(y)$, which shows that the eigenvalues of $K_t$ as an integral operator are $e^{-\lambda_i t}$. Together, these properties show that $K_t$ defines a Mercer kernel. \paragraph{} The heat kernel $K_t(x,y)$ is a natural candidate for measuring the similarity between points between $x,y\in\cM$, while respecting the geometry encoded in the metric $g$. Furthermore it is, unlike the geodesic distance, a Mercer kernel -- a fact that enables its use in statistical kernel machines. When this kernel is used for classification, as in our text classification experiments presented in Section~\ref{sec:experimentsDiff}, the discriminant function $y_t(x) = \sum_i \alpha_i \kern.1ex y_i \kern.1ex K_t(x, x_i)$ can be interpreted as the solution to the heat equation with initial temperature $y_0(x_i) = \alpha_i\, y_i$ on labeled data points $x_i$, and $y_0(x)=0$ elsewhere. \subsubsection{The parametrix expansion} \label{sec:expansion} For most geometries, there is no closed form solution for the heat kernel. However, the short time behavior of the solutions can be studied using an asymptotic expansion, called the {\it parametrix expansion\/}. In fact, the existence of the heat kernel, as asserted in the above theorem, is most directly proven by first showing the existence of the parametrix expansion. Although it is local, the parametrix expansion contains a wealth of geometric information, and indeed much of modern differential geometry, notably index theory, is based upon this expansion and its generalizations. In Section~\ref{sec:experimentsDiff} we will employ the first-order parametrix expansion for text classification. \paragraph{} Recall that the heat kernel on $n$-dimensional Euclidean space is given by \begin{equation} K_t^{\mbox{\tiny Euclid}}(x,y)\; =\; (4\pi t)^{-\frac{n}{2}} \exp\left( - \frac{\|x - y\|^2}{4t}\right) \end{equation} where $\|x-y\|^2 = \sum_{i=1}^n \left| x_i - y_i\right|^2$ is the squared Euclidean distance between $x$ and $y$. The parametrix expansion approximates the heat kernel locally as a correction to this Euclidean heat kernel. It is given by \begin{equation} P_t^{(m)}(x,y) = (4\pi t)^{-\frac{n}{2}} \exp\left( - \frac{d^2(x,y)}{4t}\right) \left( \psi_0(x,y) + \psi_1(x,y) t + \cdots + \psi_m(x,y) t^m\right) \label{eq:param} \end{equation} where $d$ is the geodesic distance and $\psi_k$ are recursively obtained by solving the heat equation approximately to order $t^m$, for small diffusion time $t$. Denoting $K_t^{(m)}(x,y)=P_t^{(m)}(x,y)$ we thus obtain an approximation for the heat kernel, that converges as $t\to 0$ and $x\to y$. For further details refer to \cite{Schoen1994,Rosenberg:97}. \ignore{ \paragraph{} Let $r = d(x,y)$ denote the length of the radial geodesic from $x$ to $y\in V_x$ in the normal coordinates defined by the exponential map. For any functions $f(r)$ and $h(r)$ of $r$, it can be shown that \begin{eqnarray} \Delta f &=& \frac{d^2 f}{d r^2} + \frac{d\left(\log \sqrt{\det g}\right)}{d r} \frac{d f}{dr} \\ \Delta(f h) &=& f\Delta h + h \Delta f + 2 \frac{df}{dr} \frac{dh}{dr} \end{eqnarray} Starting from these basic relations, some calculus shows that \begin{equation} \left(\Delta - \frac{\partial}{\partial t}\right) P_t^{(m)} = \left(t^m \, \Delta \psi_m \right) \left(4\pi t\right)^{-\frac{n}{2}} \exp\left(-\frac{r^2}{4t}\right) \label{eq:approxheat} \end{equation} when $\psi_k$ are defined recursively as \begin{eqnarray} \label{eq:psi} \psi_0 &=& \left(\frac{\sqrt{\det g}}{r^{n-1}}\right)^{-\frac{1}{2}} \\ \psi_k &=& r^{-k}\psi_0 \int_{0}^r \psi_0^{-1} \left(\Delta \phi_{k-1} \right)\, s^{k-1} ds \kern10pt \mbox{for $k > 0$} \end{eqnarray} With this recursive definition of the functions $\psi_k$, the expansion \eqref{eq:param}, which is defined only locally, is then extended to all of $M\times M$ by smoothing with a ``cut-off function'' $\eta$, with the specification that $\eta:\R_+\longrightarrow [0,1]$ is $C^\infty$ and \begin{equation} \eta(r) = \begin{cases} 0 & r \geq 1 \\ 1 & r \leq \frac{1}{2} \end{cases} \end{equation} Thus, the order-$m$ parametrix is defined as \begin{equation} K^{(m)}_t(x,y) = \eta(d(x,y)) \, P_t^{(m)}(x,y) \end{equation} As suggested by equation \eqref{eq:approxheat}, $K_t^{(m)}$ is an approximate solution to the heat equation, and satisfies $K_t(x,y) = K_t^{(m)}(x,y) + O(t^m)$ for $x$ and $y$ sufficiently close; in particular, the parametrix is not unique. For further details we refer to \emcite{Schoen1994,Rosenberg:97}. } \paragraph{} While the parametrix ${K}_t^{(m)}$ is not in general positive-definite, and therefore does not define a Mercer kernel, it is positive-definite for $t$ sufficiently small. In particular, define $f(t) = \min\text{spec} \left(K_t^{(m)}\right)$, where $\min\text{spec}$ denotes the smallest eigenvalue. Then $f$ is a continuous function with $f(0) = 1$ since $K_0^{(m)} = I$. Thus, there is some time interval $[0,\epsilon)$ for which $K_t^{(m)}$ is positive-definite in case $t\in [0,\epsilon)$. \paragraph{} The following two basic examples illustrate the geometry of the Fisher information metric and the associated diffusion kernel it induces on a statistical manifold. Under the Fisher information metric, the spherical normal family corresponds to a manifold of constant negative curvature, and the multinomial corresponds to a manifold of constant positive curvature. The multinomial will be the most important example that we develop, and we will report extensive experiments with the resulting kernels in Section~\ref{sec:experimentsDiff}. \paragraph{} \def\myrho{r} The heat kernel on the hyperbolic space $\H^n$ has the following closed form \cite{Grigoryan98}. For odd $n=2m+1$ it is given by \begin{equation} \label{eq:hyper_a} K_t(x,x') = \displaystyle \frac{(-1)^{m}}{2^m\pi^m}\frac{1}{\sqrt{4\pi t}} \left(\frac{1}{\sinh \myrho}\frac{\partial}{\partial \myrho}\right)^m \exp\left(-m^2 t - \frac{\myrho^2}{4t}\right) \end{equation} and for even $n=2m+2$ it is given by \begin{equation} K_t(x,x') = \label{eq:hyper_b} \displaystyle \frac{(-1)^{m}}{2^m\pi^m}\frac{\sqrt{2}}{\sqrt{4\pi t}^3} \left(\frac{1}{\sinh \myrho}\frac{\partial}{\partial \myrho}\right)^m \int_{\myrho}^\infty \displaystyle \frac{s \exp\left(-\frac{(2m+1)^2t}{4}-\frac{s^2}{4t}\right)}{\sqrt{\cosh s - \cosh \myrho}} \, \text{d}s \end{equation} where $\myrho = d(x,x')$ is the geodesic distance between the two points in $\H^n$ given by equation \eqref{eq:hypDist}. If only the mean $\theta=\mu$ is unspecified, then the associated kernel is the standard RBF or Gaussian kernel. The hyperbolic geometry is illustrated in Figure~\ref{fig:normalexamples} where decision boundaries of SVM with the diffusion kernel are plotted for both Euclidean and hyperbolic geometry. \begin{figure} \begin{center} \begin{tabular}{ccc} \includegraphics[width=2.7in]{diffusionFigures/normalExample2} && \includegraphics[width=2.7in]{diffusionFigures/normalExample1} \end{tabular} \caption{Example decision boundaries for a kernel-based classifier using information diffusion kernels for spherical normal geometry with $d=2$ (right), which has constant negative curvature, compared with the standard Gaussian kernel for flat Euclidean space (left). Two data points are used, simply to contrast the underlying geometries. The curved decision boundary for the diffusion kernel can be interpreted statistically by noting that as the variance decreases the mean is known with increasing certainty.} \label{fig:normalexamples} \end{center} \end{figure} \paragraph{} Unlike the explicit expression for the Gaussian geometry discussed above, there is no explicit form for the heat kernel on the sphere, nor on the positive orthant of the sphere. We will therefore resort to the parametrix expansion to derive an approximate heat kernel for the multinomial geometry. \paragraph{} \ignore{ Recall from Section~\ref{sec:expansion} that the parametrix is obtained according to the local expansion given in equation \eqref{eq:param}, and then extending this smoothly to zero outside a neighborhood of the diagonal, as defined by the exponential map. As we have just derived, this results in the following parametrix for the multinomial family: \begin{equation} P_t^{(m)}(\theta,\theta') = (4\pi t)^{-\frac{n}{2}} \exp\left( - \frac{\arccos^2( \sqrt{\theta}\cdot \sqrt{\theta'})}{t}\right) \left( \psi_0(\theta,\theta') + \cdots + \psi_m(\theta,\theta') t^m\right) \end{equation} The first-order expansion is thus obtained as \begin{equation} K_t^{(0)}(\theta,\theta') = \eta(d(\theta,\theta'))\, P_t^{(0)}(\theta,\theta') \end{equation} } For the $n$-sphere it can be shown \cite{Berger:71} that the function $\psi_0$ of in the parametrix expansion, which is the leading order correction of the Gaussian kernel under the Fisher information metric, is given by \begin{align*} \psi_0(r) &= \left(\frac{\sqrt{\det g}}{r^{n-1}}\right)^{-\frac{1}{2}} = \left(\frac{\sin r }{r}\right)^{-\frac{(n-1)}{2}} \\ &= 1 + \frac{(n-1)}{12}\,r^2 + \frac{(n-1)(5n-1)}{1440}\,r^4 + O(r^6). %\label{sinexp} \end{align*} \ignore{ Thus, the leading order parametrix for the multinomial diffusion kernel is \begin{equation} P_t^{(0)}(\theta,\theta') = (4\pi t)^{-\frac{n}{2}} \exp\left( - \frac{1}{4t} d^2(\theta,\theta')\right) \left(\frac{\sin d(\theta,\theta')}{d(\theta,\theta')}\right)^{-\frac{(n-1)}{2}} \end{equation} } In our experiments we approximate the diffusion kernel using $\psi_0\equiv 1$ and obtain \begin{align} K(\theta,\theta') = (4\pi t)^{-\frac{n}{2}} e^{-\frac{1}{t}\, \arccos^2\left(\sum_i \sqrt{\theta_i\theta_i'}\right)}. \label{nosin} \end{align} In Figure~\ref{fig:examples} the kernel \eqref{nosin} is compared with the standard Euclidean space Gaussian kernel for the case of the trinomial model, $d=2$, using an SVM classifier. \paragraph{} \begin{figure} \begin{center} \begin{tabular}{ccc} \includegraphics[width=2.7in]{diffusionFigures/multExample2} && \includegraphics[width=2.7in]{diffusionFigures/multExample1} \end{tabular} \end{center} \caption{Example decision boundaries using support vector machines with information diffusion kernels for trinomial geometry on the 2-simplex (top right) compared with the standard Gaussian kernel (left).} \label{fig:examples} \end{figure} \subsection{Rounding the Simplex} The case of multinomial geometry poses some technical complications for the analysis of diffusion kernels, due to the fact that the open simplex is not complete, and moreover, its closure is not a differentiable manifold with boundary. Thus, it is technically not possible to apply several results from differential geometry, such as bounds on the spectrum of the Laplacian, as adopted in Section~\ref{sec:bounds}. We now briefly describe a technical ``patch'' that allows us to derive all of the needed analytical results, without sacrificing in practice any of the methodology that has been derived so far. The idea is to ``round the corners'' of $\overline{\Pn}$ to obtain a compact manifold with boundary, and that closely approximates the original simplex $\Pn$. \paragraph{} For $\epsilon > 0$, let $B_\epsilon(x) = \{ y\given \|x-y\| < \epsilon\}$ be the open Euclidean ball of radius $\epsilon$ centered at~$x$ and $C_\epsilon(\P_n)$ be \begin{equation} C_\epsilon(\P_n) = \left\{ x\in\overline{\Pn} \;:\; B_\epsilon(x) \subset \Pn\right\} \end{equation} The $\epsilon$-rounded simplex is then defined as the closure of \begin{equation} \P_n^\epsilon = \bigcup_{x\,\in\, C_\epsilon(\P_n)} B_\epsilon(x). \end{equation} The above rounding procedure that yields $\overline{\P_2^\epsilon}$ is suggested by Figure~\ref{fig:rounding}. Note that in general the $\epsilon$-rounded simplex $\overline{\P_n^\epsilon}$ will contain points with a single, but not more than one component having zero probability and it forms a compact manifold with boundary whose image under the isometry $F$ described above is a compact submanifold with boundary of the $n$-sphere. Since, by choosing $\epsilon$ small enough we can approximate $\overline{\Pn}$ arbitrarily well (both in the Euclidean and geodesic distances), no harm is done by assuming that we are dealing with a rounded compact manifold with a boundary. \begin{figure} \begin{center} \begin{tabular}{ccc} \includegraphics[scale=0.22]{diffusionFigures/rounding1}& \includegraphics[scale=0.22]{diffusionFigures/rounding2}& \includegraphics[scale=0.22]{diffusionFigures/rounding3} \end{tabular} \end{center} \caption{Rounding the simplex. Since the closed simplex is not a manifold with boundary, we carry out a ``rounding'' procedure to remove edges and corners. The $\epsilon$-rounded simplex is the closure of the union of all $\epsilon$-balls lying within the open simplex.} \label{fig:rounding} \end{figure} %------------------------------------------------------------- \subsection{Spectral Bounds on Covering Numbers and Rademacher Averages} \label{sec:bounds} We now turn to establishing bounds on the generalization performance of kernel machines that use information diffusion kernels. We begin by adopting the approach of \emcite{Guo2002}, estimating covering numbers by making use of bounds on the spectrum of the Laplacian on a Riemannian manifold, rather than on VC dimension techniques; these bounds in turn yield bounds on the expected risk of the learning algorithms. Our calculations give an indication of how the underlying geometry influences the entropy numbers, which are inverse to the covering numbers. We then show how bounds on Rademacher averages may be obtained by plugging in the spectral bounds from differential geometry. The primary conclusion that is drawn from these analyses is that from the point of view of generalization error bounds, information diffusion kernels behave essentially the same as the standard Gaussian kernel. \subsubsection{Covering Numbers} We begin by recalling the main result of \emcite{Guo2002}, modifying their notation slightly to conform with ours. Let $\cX\subset\R^d$ be a compact subset of $d$-dimensional Euclidean space, and suppose that $K:\cX\times\cX\longrightarrow\R$ is a Mercer kernel. Denote by $\lambda_1\geq \lambda_2\geq \cdots \geq 0$ the eigenvalues of $K$, i.e., of the mapping $f\mapsto \int_{\cX} K(\cdot, y)\, f(y)\, dy$, and let $\psi_j(\cdot)$ denote the corresponding eigenfunctions. We assume that $C_K \deff \sup_j \|\psi_j\|_\infty < \infty$. \paragraph{} Given $m$ points $x_i\in\cX$, the kernel hypothesis class for $\bb x = \{x_i\}$ with weight vector bounded by $R$ is defined as the collection of functions on $\bb x$ given by \begin{equation} \cF_R(\bb x) = \left\{f : f(x_i) = \dotp{w}{\Phi(x_i)}\; \mbox{for some $\|w\|\leq R$}\right\} \end{equation} where $\Phi(\cdot)$ is the mapping from $M$ to feature space defined by the Mercer kernel, and $\dotp{\cdot}{\cdot}$ and $\|\cdot\|$ denote the corresponding Hilbert space inner product and norm. It is of interest to obtain uniform bounds on the covering numbers $\cN(\epsilon,\cF_R(\bb x))$, defined as the size of the smallest $\epsilon$-cover of $\cF_R(\bb x)$ in the metric induced by the norm $\|f\|_{\infty, \bb x} = \max_{i=1,\ldots,m} |f(x_i)|$. The following is the main result of \cite{Guo2002}. \begin{thrm} Given an integer $n\in\N$, let $j_n^*$ denote the smallest integer $j$ for which \begin{equation} \lambda_{j+1} < \left(\frac{\lambda_1 \cdots \lambda_j}{n^2}\right)^{\textstyle \frac{1}{j}} \end{equation} and define \begin{equation} \epsilon_n^* = 6 C_K R \sqrt{j_{n}^{*} \left({\frac{\lambda_1 \cdots \lambda_{j_n^*}}{n^2}}\right)^{\textstyle\frac{1}{j_n^*}} + \sum_{i={j_n^*}}^{\infty} \lambda_i}\,. \end{equation} Then $\sup_{\{x_i\}\in M^m} \cN(\epsilon_n^*, \cF_R(\bb x)) \leq n$. \end{thrm} \vskip5pt To apply this result, we will obtain bounds on the indices $j_n^*$ using spectral theory in Riemannian geometry. The following bounds on the eigenvalues of the Laplacian are due to \cite{LiYau:80}. \def\dsq{D^2} \def\dsq{V^{\frac{2}{d}}} \begin{thrm} \label{thm:eigenbounds} Let $M$ be a compact Riemannian manifold of dimension $d$ with non-negative Ricci curvature, %$\mbox{Ric}(M)\geq 0$. %and assume that the boundary of $M$ is convex. and let $\,0 < \mu_1 \leq \mu_2\leq \cdots$ denote the eigenvalues of the Laplacian with Dirichlet boundary conditions. Then \begin{equation} \label{libounds} c_1(d) \left(\frac{j}{V}\right)^{\frac{2}{d}} \leq\; \mu_j \;\leq \; c_2(d) \left(\frac{j+1}{V}\right)^{\frac{2}{d}} \end{equation} %\begin{equation} %c'_1(d) \left(\frac{j+1}{V}\right)^{\frac{2}{d}} \leq \; \eta_j \leq \; %c'_2(d) \left(\frac{j+2}{V}\right)^{\frac{2}{d}} %\end{equation} where $V$ is the volume of $M$ and $c_1$ and $c_2$ are constants depending only on the dimension. \end{thrm} Note that the manifold of the multinomial model satisfies the conditions of this theorem. Using these results we can establish the following bounds on covering numbers for information diffusion kernels. We assume Dirichlet boundary conditions; a similar result can be proven for Neumann boundary conditions. We include the constant $V=\mbox{vol}(M)$ and diffusion coefficient $t$ in order to indicate how the bounds depend on the geometry. \def\ts{\textstyle} \begin{thrm} Let $M$ be a compact Riemannian manifold, with volume $V$, satisfying the conditions of Theorem~\ref{thm:eigenbounds}. Then the covering numbers for the Dirichlet heat kernel $K_{\kern-1pt t}$ on $M$ satisfy \begin{equation} \label{eq:nbound} \log \cN(\epsilon, \cF_R(\bb x)) = O\left(\left(\frac{V}{\ts t^{\frac{d}{2}}}\right) \log^{\ts\frac{d+2}{2}}\left(\frac{1}{\epsilon}\right) \right) \end{equation} \end{thrm} \begin{proof} By the lower bound in Theorem~\ref{thm:eigenbounds}, the Dirichlet eigenvalues of the heat kernel $K_t(x,y)$, which are given by $\lambda_j = e^{-t\mu_j}$, satisfy $ \log \lambda_j \leq - {t c_1(d)} \left(\frac{j}{V}\right)^{\frac{2}{d}}$. Thus, \begin{equation} - \frac{1}{j} \log \left(\frac{\lambda_1\cdots\lambda_j}{n^2}\right) \;\geq\; \frac{t c_1}{j} \sum_{i=1}^{j} \left(\frac{i}{V}\right)^{\ts \frac{2}{d}} + \frac{2}{j} \log n \;\geq\; {t c_1} \frac{d}{d+2} \left(\frac{j}{V}\right)^{\ts \frac{2}{d}} + \frac{2}{j} \log n \end{equation} where the second inequality comes from $\sum_{i=1}^{j} i^p \geq \int_{0}^j x^p \, dx = \frac{j^{p+1}}{p+1} $. Now using the upper bound of Theorem~\ref{thm:eigenbounds}, the inequality $j_n^* \leq j$ will hold if \begin{equation} {t c_2} \left(\frac{j+2}{V}\right)^{\ts \frac{2}{d}} \;\geq\; -\log \lambda_{j+1} \;\geq\; {t c_1} \frac{d}{d+2} \left(\frac{j}{V}\right)^{\ts \frac{2}{d}} + \frac{2}{j} \log n \end{equation} or equivalently \begin{equation} \frac{t c_2}{\dsq} \left( j (j+2)^{\ts \frac{2}{d}} - \frac{c_1}{c_2}\frac{d}{d+2} j^{\ts \frac{d+2}{d}}\right) \;\geq\; 2 \log n \end{equation} The above inequality will hold in case \begin{equation} j \geq \left\lceil\left( \frac{2\dsq}{t (c_2 - c_1\frac{d}{d+2})} \log n\right)^{\ts\frac{d}{d+2}}\right\rceil \;\geq \left\lceil\left( \frac{\dsq (d+2)}{t c_1} \log n\right)^{\ts\frac{d}{d+2}}\right\rceil \end{equation} since we may assume that $c_2\geq c_1$; thus, $ j_n^* \leq \left\lceil \overline c_1\left( \frac{\dsq}{t} \log n \right)^{\frac{d}{d+2}}\right\rceil$ for a new constant $\overline c_1(d)$. Plugging this bound on $j_n^*$ into the expression for $\epsilon_n^*$ in Theorem~2 and using \begin{equation} \sum_{i={j_n^*}}^{\infty} e^{-i^\frac{2}{d}} = O\left(e^{-{j_n^*}^\frac{2}{d}}\right) \end{equation} we have after some algebra that \begin{equation} \log\left(\frac{1}{\epsilon_n}\right) = \Omega\left(\left(\frac{t}{V^\frac{2}{d}}\right)^{\frac{d}{d+2}} \log^{\frac{2}{d+2}} n\right) \end{equation} Inverting the above expression in $\log n$ gives equation (\ref{eq:nbound}). \end{proof} We note that Theorem~4 of \cite{Guo2002} can be used to show that this bound does not, in fact, depend on $m$ and $\bb x$. Thus, for fixed $t$ the covering numbers scale as $\log\cN(\epsilon, \cF) = O\left(\log^{\frac{d+2}{2}} \left(\frac{1}{\epsilon}\right)\right)$, and for fixed $\epsilon$ they scale as $\log\cN(\epsilon, \cF) = O\left(t^{-\frac{d}{2}}\right)$ in the diffusion time $t$. \subsubsection{Rademacher Averages} We now describe a different family of generalization error bounds that can be derived using the machinery of Rademacher averages \emcite{Bartlett:02,Bartlett:02b}. The bounds fall out directly from the work of \cite{Mendelson:03} on computing local averages for kernel-based function classes, after plugging in the eigenvalue bounds of Theorem~3. \paragraph{} As seen above, covering number bounds are related to a complexity term of the form \begin{equation} C(n) = \sqrt{j_{n}^{*} \left({\frac{\lambda_1 \cdots \lambda_{j_n^*}}{n^2}}\right)^{\textstyle\frac{1}{j_n^*}} + \sum_{i={j_n^*}}^{\infty} \lambda_i} \end{equation} In the case of Rademacher complexities, risk bounds are instead controlled by a similar, yet simpler expression of the form \begin{equation} C(r) = \sqrt{ j^*_r \,r + \sum_{i={j^*_r}}^{\infty} \lambda_i} \end{equation} where now $j^*_r$ is the smallest integer $j$ for which $\lambda_{j} < r$ \emcite{Mendelson:03}, with $r$ acting as a parameter bounding the error of the family of functions. To place this into some context, we quote the following results from \cite{Bartlett:02b} and \cite{Mendelson:03}, which apply to a family of loss functions that includes the quadratic loss; we refer to \cite{Bartlett:02b} for details on the technical conditions. \paragraph{} Let $(X_1,Y_1), (X_2,Y_2)\ldots,(X_n,Y_n)$ be an independent sample from an unknown distribution~$P$ on ${\mathcal X}\times {\mathcal Y}$, where ${\mathcal Y}\subset \R$. For a given loss function $\ell:{\mathcal Y}\times {\mathcal Y}\rightarrow \R$, and a family $\cF$ of measurable functions $f:{\mathcal X}\rightarrow {\mathcal Y}$, the objective is to minimize the expected loss $E[\ell(f(X),Y)]$. Let $E\ell_{f^*} = \inf_{f\in\cF} E\ell_f$, where $\ell_f(X,Y) = \ell(f(X),Y)$, and let $\hat f$ be any member of $\cF$ for which $E_n \ell_{\hat f} = \inf_{f\in\cF} E_n \ell_f$ where $E_n$ denotes the empirical expectation. The {\it Rademacher average\/} of a family of functions $\mathfrak{G} = \left\{g: {\mathcal X}\rightarrow \R\right\}$ is defined as the expectation $ER_n \mathfrak{G} = E\left[\sup_{g\in\mathfrak{G}} R_n g\right]$ with $R_n g = \frac{1}{n} \sum_{i=1}^n \sigma_i \,g(X_i)$, where $\sigma_1,\ldots,\sigma_n$ are independent {\it Rademacher\/} random variables; that is, $p(\sigma_i=1) = p(\sigma_i=-1) = \frac{1}{2}$. \begin{thrm} Let $\cF$ be a convex class of functions and define $\psi$ by \begin{equation} \psi(r) = a\, ER_n\left\{ f\in\cF : E(f-f^*)^2 \leq r\right\} + \frac{b\, x}{n} \end{equation} where $a$ and $b$ are constants that depend on the loss function $\ell$. Then when $r\geq \psi(r)$, \begin{equation} E\left(\ell_{\hat f} - \ell_{f^*}\right) \leq c\,r + \frac{d\, x}{n} \end{equation} with probability at least $1-e^{-x}$, where $c$ and $d$ are additional constants. Moreover, suppose that $K$ is a Mercer kernel and $\cF = \left\{ f\in {\mathcal H}_K : \|f\|_K \leq 1\right\}$ is the unit ball in the reproducing kernel Hilbert space associated with $K$. Then \begin{eqnarray} \psi(r) \leq a \sqrt{\frac{2}{n} \sum_{j=1}^\infty \min\{r,\lambda_j\}} + \frac{b x}{n} \end{eqnarray} \end{thrm} Thus, to bound the excess risk for kernel machines in this framework it suffices to bound the term \begin{equation} \widetilde\psi(r) = \sqrt{\sum_{j=1}^\infty \min\{r,\lambda_j\}} \end{equation} involving the spectrum. Given bounds on the eigenvalues, this is typically easy to do. \begin{thrm} Let $M$ be a compact Riemannian manifold, satisfying the conditions of Theorem~\ref{thm:eigenbounds}. Then the Rademacher term $\widetilde \psi$ for the Dirichlet heat kernel $K_{\kern-1pt t}$ on $M$ satisfies \begin{equation} \widetilde\psi(r) \leq C \sqrt{\left(\frac{r}{t^\frac{d}{2}}\right) \log^{\frac{d}{2}} \left(\frac{1}{r}\right)} \end{equation} for some constant $C$ depending on the geometry of $M$. \end{thrm} \begin{proof} We have that \begin{eqnarray} \widetilde\psi^2(r) &=& \sum_{j=1}^\infty \min\{r, \lambda_j\} \\ &= & j_r^* \, r + \sum_{j=j_r^*}^\infty e^{-t\mu_j} \\ &\leq & j_r^*\, r + \sum_{j=j_r^*}^\infty e^{-t c_1 j^{\frac{2}{d}}} \\ &\leq & j_r^*\, r + C e^{-t c_1 {j_r^*}^{\frac{2}{d}}} \end{eqnarray} for some constant $C$, where the first inequality follows from the lower bound in Theorem~\ref{thm:eigenbounds}. But $j_r^* \leq j$ in case $\log \lambda_{j+1} > r$, or, again from Theorem~\ref{thm:eigenbounds}, if \begin{eqnarray} t\,c_2(j+1)^{\frac{2}{d}} \leq -\log \lambda_j < \log\frac{1}{r} \end{eqnarray} or equivalently, \begin{equation} j_r^* \leq \frac{C'}{t^\frac{d}{2}}\log^{\frac{d}{2}}\left(\frac{1}{r}\right) \end{equation} It follows that \begin{eqnarray} \widetilde \psi^2(r) &\leq& C{''} \left(\frac{r}{t^{\frac{d}{2}}}\right)\, \log^{\frac{d}{2}}\left(\frac{1}{r}\right) \end{eqnarray} for some new constant $C{''}$. \end{proof} From this bound, it can be shown that, with high probability, \begin{eqnarray} E\left(\ell_{\hat f} - \ell_{f^*}\right) = O\left(\frac{\log^{\frac{d}{2}} n}{n}\right) \end{eqnarray} which is the behavior expected of the Gaussian kernel for Euclidean space. \paragraph{} Thus, for both covering numbers and Rademacher averages, the resulting bounds are essentially the same as those that would be obtained for the Gaussian kernel on the flat $d$-dimensional torus, which is the standard way of ``compactifying'' Euclidean space to get a Laplacian having only discrete spectrum; the results of \cite{Guo2002} are formulated for the case $d=1$, corresponding to the circle. While the bounds for information diffusion kernels were derived for the case of positive curvature, which apply to the special case of the multinomial, similar bounds for general manifolds with curvature bounded below by a negative constant should also be attainable. %--------------------------------------------------------------- \subsection{Experimental Results for Text Classification} \label{sec:experimentsDiff} In this section we present the application of multinomial diffusion kernels to the problem of text classification. Text processing can be subject to some of the ``dirty laundry'' referred to in the introduction---documents are cast as Euclidean space vectors with special weighting schemes that have been empirically honed through applications in information retrieval, rather than inspired from first principles. However for text, the use of multinomial geometry is natural and well motivated; our experimental results offer some insight into how useful this geometry may be for classification. \paragraph{} We consider several embeddings $\hat{\theta}:\cX\to\Theta$ of documents in bag of words representation into the probability simplex. The term frequency (tf) representation uses normalized counts; the corresponding embedding is the maximum likelihood estimator for the multinomial distribution \begin{equation} \widehat{\theta}_{\text{tf}}(x) = \left(\frac{x_1}{\sum_i x_i},\ldots,\frac{x_{n+1}}{\sum_i x_i} \right). \end{equation} Another common representation is based on term frequency, inverse document frequency (tf-idf). This representation uses the distribution of terms across documents to discount common terms; the document frequency $df_v$ of term $v$ is defined as the number of documents in which term $v$ appears. Although many variants have been proposed, one of the simplest and most commonly used embeddings is \begin{equation} \widehat{\theta}_{\text{tf-idf}}(x) = \left(\frac{x_1 \log (D/df_1)}{\sum_i x_i \log (D/df_i)} ,\ldots, \frac{x_{n+1} \log (D/df_{n+1})}{\sum_i x_i \log (D/df_i)}\right) \end{equation} where $D$ is the number of documents in the corpus. \paragraph{} In text classification applications the tf and tf-idf representations are typically normalized to unit length in the $L_2$ norm rather than the $L_1$ norm, as above \emcite{Joachims2000}. For example, the tf representation with $L_2$ normalization is given by \begin{equation} x \mapsto \left(\frac{x_1}{\sum_i x_i^2},\ldots,\frac{x_{n+1}}{\sum_i x_i^2} \right) \end{equation} and similarly for tf-idf. When used in support vector machines with linear or Gaussian kernels, $L_2$-normalized tf and tf-idf achieve higher accuracies than their $L_1$-normalized counterparts. However, for the diffusion kernels, $L_1$ normalization is necessary to obtain an embedding into the simplex. These different embeddings or feature representations are compared in the experimental results reported below. \paragraph{} The three kernels that we compare are the linear kernel \begin{equation} K^{\text{Lin}}(\theta,\theta') \sum_{v=1}^{n+1} \theta_v \,\theta'_v, \end{equation} the Gaussian kernel \begin{equation} K_\sigma^{\text{Gauss}}(\theta',\theta')\; =\; (2\pi \sigma)^{-\frac{n+1}{2}} \exp\left( - \frac{\sum_{i=1}^{n+1}|\theta_i-\theta_i'|^2}{2\sigma^2}\right) \end{equation} and the multinomial diffusion kernel approximation \begin{equation} K_t^{\mbox{\tiny Mult}}(\theta,\theta') = (4\pi t)^{-\frac{n}{2}} \exp\left( -\frac{1}{t}\, \arccos^2\left( \sum_{i=1}^{n+1} \sqrt{\theta_i\theta_i'}\right)\right). \end{equation} \paragraph{} In our experiments, the multinomial diffusion kernel using the tf embedding was compared to the linear or Gaussian kernel with tf and tf-idf embeddings using a support vector machine classifier on the WebKB and Reuters-21578 collections, which are standard data sets for text classification. \paragraph{} The WebKB dataset contains web pages found on the sites of four universities \emcite{Craven1998}. The pages were classified according to whether they were student, faculty, course, project or staff pages; these categories contain 1641, 1124, 929, 504 and 137 instances, respectively. Since only the \texttt{student}, \texttt{faculty}, \texttt{course} and \texttt{project} classes contain more than 500 documents each, we restricted our attention to these classes. The Reuters-21578 dataset is a collection of newswire articles classified according to news topic \emcite{Lewis94}. Although there are more than 135 topics, most of the topics have fewer than 100 documents; for this reason, we restricted our attention to the following five most frequent classes: \texttt{earn}, \texttt{acq}, \texttt{moneyFx}, \texttt{grain} and \texttt{crude}, of sizes 3964, 2369, 717, 582 and 578 documents, respectively. \paragraph{} For both the WebKB and Reuters collections we created two types of binary classification tasks. In the first task we designate a specific class, label each document in the class as a ``positive'' example, and label each document on any of the other topics as a ``negative'' example. In the second task we designate a class as the positive class, and choose the negative class to be the most frequent remaining class (\texttt{student} for WebKB and \texttt{earn} for Reuters). In both cases, the size of the training set is varied while keeping the proportion of positive and negative documents constant in both the training and test set. \begin{figure} \begin{center} \begin{tabular}{cc} \includegraphics[scale=0.35]{diffusionFigures/webkb-1-vs-all/webkb-course-TF.eps}& \includegraphics[scale=0.35]{diffusionFigures/webkb-1-vs-all/webkb-course-TFIDF.eps}\\ \includegraphics[scale=0.35]{diffusionFigures/webkb-1-vs-all/webkb-faculty-TF.eps}& \includegraphics[scale=0.35]{diffusionFigures/webkb-1-vs-all/webkb-faculty-TFIDF.eps}\\ \includegraphics[scale=0.35]{diffusionFigures/webkb-1-vs-all/webkb-project-TF.eps}& \includegraphics[scale=0.35]{diffusionFigures/webkb-1-vs-all/webkb-project-TFIDF.eps}\\ \includegraphics[scale=0.35]{diffusionFigures/webkb-1-vs-all/webkb-student-TF.eps}& \includegraphics[scale=0.35]{diffusionFigures/webkb-1-vs-all/webkb-student-TFIDF.eps}\\ \end{tabular} \end{center} \caption{Experimental results on the WebKB corpus, using SVMs for linear (dotted) and Gaussian (dash-dotted) kernels, compared with the diffusion kernel for the multinomial (solid). Classification error for the task of labeling \texttt{course} (top row), \texttt{faculty} (second row), \texttt{project} (third row), and \texttt{student}(bottom row) is shown in these plots, as a function of training set size. The left plot uses tf representation and the right plot uses tf-idf representation. The curves shown are the error rates averaged over 20-fold cross validation. The results for the other ``1 vs.~all'' labeling tasks are qualitatively similar, and therefore not shown.} \label{fig:experiments1} \end{figure} \begin{figure} \begin{center} \begin{tabular}{cc} \includegraphics[scale=0.35]{diffusionFigures/webkb-1-vs-student/course-vs-student-CV-TF.eps}& \includegraphics[scale=0.35]{diffusionFigures/webkb-1-vs-student/course-vs-student-CV-TFIDF.eps}\\ \includegraphics[scale=0.35]{diffusionFigures/webkb-1-vs-student/faculty-vs-student-CV-TF.eps}& \includegraphics[scale=0.35]{diffusionFigures/webkb-1-vs-student/faculty-vs-student-CV-TFIDF.eps}\\ \includegraphics[scale=0.35]{diffusionFigures/webkb-1-vs-student/project-vs-student-CV-TF.eps}& \includegraphics[scale=0.35]{diffusionFigures/webkb-1-vs-student/project-vs-student-CV-TFIDF.eps}\\ \end{tabular} \end{center} \caption{Results on the WebKB corpus, using SVMs for linear (dotted) and Gaussian (dash-dotted) kernels, compared with the diffusion kernel (solid). The tasks are \texttt{course} vs. \texttt{student} (top row), \texttt{faculty} vs. \texttt{student} (top row) and \texttt{project} vs. \texttt{student} (top row). The left plot uses tf representation and the right plot uses tf-idf representation. Results for other label pairs are qualitatively similar. The curves shown are the error rates averaged over 20-fold cross validation.} \label{fig:experiments4} \end{figure} \begin{figure} \begin{center} \begin{tabular}{cc} \includegraphics[scale=0.35]{diffusionFigures/reuters-1-vs-all/reuters-earn-vs-all-CV-TF.eps}& \includegraphics[scale=0.35]{diffusionFigures/reuters-1-vs-all/reuters-earn-vs-all-CV-TFIDF.eps}\\ \includegraphics[scale=0.35]{diffusionFigures/reuters-1-vs-all/reuters-acq-vs-all-CV-TF.eps}& \includegraphics[scale=0.35]{diffusionFigures/reuters-1-vs-all/reuters-acq-vs-all-CV-TFIDF.eps}\\ \includegraphics[scale=0.35]{diffusionFigures/reuters-1-vs-all/reuters-moneyFx-vs-all-CV-TF.eps}& \includegraphics[scale=0.35]{diffusionFigures/reuters-1-vs-all/reuters-moneyFx-vs-all-CV-TFIDF.eps}\\ \end{tabular} \end{center} \caption{Experimental results on the Reuters corpus, using SVMs for linear (dotted) and Gaussian (dash-dotted) kernels, compared with the diffusion kernel (solid). The tasks are classifying \texttt{earn} (top row), \texttt{acq} (second row), \texttt{moneyFx} (bottom row) vs. the rest. Plots for the other classes are qualitatively similar. The left column uses tf representation and the right column uses tf-idf. The curves shown are the error rates averaged over 20-fold cross validation.} \label{fig:experiments2} \end{figure} \begin{figure} \begin{center} \begin{tabular}{cc} \includegraphics[scale=0.35]{diffusionFigures/reuters-1-vs-earn/acq-vs-earn-TF.eps}& \includegraphics[scale=0.35]{diffusionFigures/reuters-1-vs-earn/acq-vs-earn-TFIDF.eps}\\ \includegraphics[scale=0.35]{diffusionFigures/reuters-1-vs-earn/moneyFx-vs-earn-TF.eps}& \includegraphics[scale=0.35]{diffusionFigures/reuters-1-vs-earn/moneyFx-vs-earn-TFIDF.eps}\\ \includegraphics[scale=0.35]{diffusionFigures/reuters-1-vs-earn/grain-vs-earn-TF.eps}& \includegraphics[scale=0.35]{diffusionFigures/reuters-1-vs-earn/grain-vs-earn-TFIDF.eps}\\ \includegraphics[scale=0.35]{diffusionFigures/reuters-1-vs-earn/crude-vs-earn-TF.eps}& \includegraphics[scale=0.35]{diffusionFigures/reuters-1-vs-earn/crude-vs-earn-TFIDF.eps}\\ \end{tabular} \end{center} \caption{Experimental results on the Reuters corpus, using SVMs for linear (dotted) and Gaussian (dash-dotted) kernels, compared with the diffusion (solid). The tasks are \texttt{acq} vs. \texttt{earn} (top row), \texttt{moneyFx} vs. \texttt{earn} (top row), \texttt{grain} vs. \texttt{earn} (top row), \texttt{crude} vs. \texttt{earn} (top row). The left column uses tf representation and the right column uses tf-idf representation. The left column uses tf representation and the right column uses tf-idf. The curves shown are the error rates averaged over 20-fold cross validation.} \label{fig:experiments5} \end{figure} \paragraph{} Figure~\ref{fig:experiments1} shows the test set error rate for the WebKB data, for a representative instance of the one-versus-all classification task; the designated class was \texttt{course}. The results for the other choices of positive class were qualitatively very similar; all of the results are summarized in Table~\ref{tab:webkb1}. Similarly, Figure~\ref{fig:experiments2} shows the test set error rates for two of the one-versus-all experiments on the Reuters data, where the designated classes were chosen to be \texttt{acq} and \texttt{moneyFx}. All of the results for Reuters one-versus-all tasks are shown in Table~\ref{tab:reuters1}. \def\nb#1{\hbox to 30pt{\hfill #1 \hfill}} \def\bb#1{\hbox to 30pt{\hfill \bf #1 \hfill}} \def\nb#1{#1} \def\bb#1{{\bf #1}} \def\str{${}^*$} \renewcommand{\arraystretch}{1.1} \begin{table} \renewcommand{\arraystretch}{1.1} \begin{center} \begin{small} \begin{tabular}{|c|r||ccc||ccc|} \multicolumn{2}{c}{} & \multicolumn{3}{c}{tf Representation} & \multicolumn{3}{c}{tf-idf Representation} \\[1pt] \hline Task & $\;L\;$ & Linear & Gaussian & Diffusion & Linear & Gaussian & Diffusion \\ % WEBKB 1 vs all \hline & 40 & \nb{0.1225} & \nb{0.1196} & \bb{0.0646} & \nb{0.0761} & \nb{0.0726} & \bb{0.0514} \\ & 80 & \nb{0.0809} & \nb{0.0805} & \bb{0.0469} & \nb{0.0569} & \nb{0.0564} & \bb{0.0357} \\ course vs.~all & 120 & \nb{0.0675} & \nb{0.0670} & \bb{0.0383} & \nb{0.0473} & \nb{0.0469} & \bb{0.0291} \\ & 200 & \nb{0.0539} & \nb{0.0532} & \bb{0.0315} & \nb{0.0385} & \nb{0.0380} & \bb{0.0238} \\ & 400 & \nb{0.0412} & \nb{0.0406} & \bb{0.0241} & \nb{0.0304} & \nb{0.0300} & \bb{0.0182} \\ & 600 & \nb{0.0362} & \nb{0.0355} & \bb{0.0213} & \nb{0.0267} & \nb{0.0265} & \bb{0.0162} \\ \hline \hline & 40 & \nb{0.2336} & \nb{0.2303} & \bb{0.1859} & \nb{0.2493} & \nb{0.2469} & \bb{0.1947} \\ & 80 & \nb{0.1947} & \nb{0.1928} & \bb{0.1558} & \nb{0.2048} & \nb{0.2043} & \bb{0.1562} \\ faculty vs.~all & 120 & \nb{0.1836} & \nb{0.1823} & \bb{0.1440} & \nb{0.1921} & \nb{0.1913} & \bb{0.1420} \\ & 200 & \nb{0.1641} & \nb{0.1634} & \bb{0.1258} & \nb{0.1748} & \nb{0.1742} & \bb{0.1269} \\ & 400 & \nb{0.1438} & \nb{0.1428} & \bb{0.1061} & \nb{0.1508} & \nb{0.1503} & \bb{0.1054} \\ & 600 & \nb{0.1308} & \nb{0.1297} & \bb{0.0931} & \nb{0.1372} & \nb{0.1364} & \bb{0.0933} \\ \hline \hline & 40 & \nb{0.1827} & \nb{0.1793} & \bb{0.1306} & \nb{0.1831} & \nb{0.1805} & \bb{0.1333} \\ & 80 & \nb{0.1426} & \nb{0.1416} & \bb{0.0978} & \nb{0.1378} & \nb{0.1367} & \bb{0.0982} \\ project vs.~all & 120 & \nb{0.1213} & \nb{0.1209} & \bb{0.0834} & \nb{0.1169} & \nb{0.1163} & \bb{0.0834} \\ & 200 & \nb{0.1053} & \nb{0.1043} & \bb{0.0709} & \nb{0.1007} & \nb{0.0999} & \bb{0.0706} \\ & 400 & \nb{0.0785} & \nb{0.0766} & \bb{0.0537} & \nb{0.0802} & \nb{0.0790} & \bb{0.0574} \\ & 600 & \nb{0.0702} & \nb{0.0680} & \bb{0.0449} & \nb{0.0719} & \nb{0.0708} & \bb{0.0504} \\ \hline \hline & 40 & \nb{0.2417} & \nb{0.2411} & \bb{0.1834} & \nb{0.2100} & \nb{0.2086} & \bb{0.1740} \\ & 80 & \nb{0.1900} & \nb{0.1899} & \bb{0.1454} & \nb{0.1681} & \nb{0.1672} & \bb{0.1358} \\ student vs.~all & 120 & \nb{0.1696} & \nb{0.1693} & \bb{0.1291} & \nb{0.1531} & \nb{0.1523} & \bb{0.1204} \\ & 200 & \nb{0.1539} & \nb{0.1539} & \bb{0.1134} & \nb{0.1349} & \nb{0.1344} & \bb{0.1043} \\ & 400 & \nb{0.1310} & \nb{0.1308} & \bb{0.0935} & \nb{0.1147} & \nb{0.1144} & \bb{0.0874} \\ & 600 & \nb{0.1173} & \nb{0.1169} & \bb{0.0818} & \nb{0.1063} & \nb{0.1059} & \bb{0.0802} \\ \hline \end{tabular} \end{small} \end{center} \caption{Experimental results on the WebKB corpus, using SVMs for linear, Gaussian, and multinomial diffusion kernels. The left columns use tf representation and the right columns use tf-idf representation. The error rates shown are averages obtained using 20-fold cross validation. The best performance for each training set size $L$ is shown in boldface. All differences are statistically significant according to the paired $t$ test at the 0.05 level.} \label{tab:webkb1} \end{table} \begin{table} \renewcommand{\arraystretch}{1.1} \begin{center} \begin{small} \begin{tabular}{|c|r||ccc||ccc|} \multicolumn{2}{c}{} & \multicolumn{3}{c}{tf Representation} & \multicolumn{3}{c}{tf-idf Representation} \\[1pt] \hline Task & $\;L\;$ & Linear & Gaussian & Diffusion & Linear & Gaussian & Diffusion \\ % WEBKB 1 vs student \hline & 40 & \nb{0.0808} & \nb{0.0802} & \bb{0.0391} & \nb{0.0580} & \nb{0.0572} & \bb{0.0363} \\ & 80 & \nb{0.0505} & \nb{0.0504} & \bb{0.0266} & \nb{0.0409} & \nb{0.0406} & \bb{0.0251} \\ course vs.~student & 120 & \nb{0.0419} & \nb{0.0409} & \bb{0.0231} & \nb{0.0361} & \nb{0.0359} & \bb{0.0225} \\ & 200 & \nb{0.0333} & \nb{0.0328} & \bb{0.0184} & \nb{0.0310} & \nb{0.0308} & \bb{0.0201} \\ & 400 & \nb{0.0263} & \nb{0.0259} & \bb{0.0135} & \nb{0.0234} & \nb{0.0232} & \bb{0.0159} \\ & 600 & \nb{0.0228} & \nb{0.0221} & \bb{0.0117} & \nb{0.0207} & \nb{0.0202} & \bb{0.0141} \\ \hline \hline & 40 & \nb{0.2106} & \nb{0.2102} & \bb{0.1624} & \nb{0.2053} & \nb{0.2026} & \bb{0.1663} \\ & 80 & \nb{0.1766} & \nb{0.1764} & \bb{0.1357} & \nb{0.1729} & \nb{0.1718} & \bb{0.1335} \\ faculty vs.~student & 120 & \nb{0.1624} & \nb{0.1618} & \bb{0.1198} & \nb{0.1578} & \nb{0.1573} & \bb{0.1187} \\ & 200 & \nb{0.1405} & \nb{0.1405} & \bb{0.0992} & \nb{0.1420} & \nb{0.1418} & \bb{0.1026} \\ & 400 & \nb{0.1160} & \nb{0.1158} & \bb{0.0759} & \nb{0.1166} & \nb{0.1165} & \bb{0.0781} \\ & 600 & \nb{0.1050} & \nb{0.1046} & \bb{0.0656} & \nb{0.1050} & \nb{0.1048} & \bb{0.0692} \\ \hline \hline & 40 & \nb{0.1434} & \nb{0.1430} & \bb{0.0908} & \nb{0.1304} & \nb{0.1279} & \bb{0.0863} \\ & 80 & \nb{0.1139} & \nb{0.1133} & \bb{0.0725} & \nb{0.0982} & \nb{0.0970} & \bb{0.0634} \\ project vs.~student & 120 & \nb{0.0958} & \nb{0.0957} & \bb{0.0613} & \nb{0.0870} & \nb{0.0866} & \bb{0.0559} \\ & 200 & \nb{0.0781} & \nb{0.0775} & \bb{0.0514} & \nb{0.0729} & \nb{0.0722} & \bb{0.0472} \\ & 400 & \nb{0.0590} & \nb{0.0579} & \bb{0.0405} & \nb{0.0629} & \nb{0.0622} & \bb{0.0397} \\ & 600 & \nb{0.0515} & \nb{0.0500} & \bb{0.0325} & \nb{0.0551} & \nb{0.0539} & \bb{0.0358} \\ \hline \end{tabular} \end{small} \end{center} \caption{Experimental results on the WebKB corpus, using SVMs for linear, Gaussian, and multinomial diffusion kernels. The left columns use tf representation and the right columns use tf-idf representation. The error rates shown are averages obtained using 20-fold cross validation. The best performance for each training set size $L$ is shown in boldface. All differences are statistically significant according to the paired $t$ test at the 0.05 level.} \label{tab:webkb2} \end{table} \begin{table} \renewcommand{\arraystretch}{1.1} \begin{center} \begin{small} \begin{tabular}{|c|r||ccc||ccc|} \multicolumn{2}{c}{} & \multicolumn{3}{c}{tf Representation} & \multicolumn{3}{c}{tf-idf Representation} \\[1pt] \hline Task & $\;L\;$ & Linear & Gaussian & Diffusion & Linear & Gaussian & Diffusion \\ % Reuters 1 vs all \hline & 80 & \nb{0.1107} & \nb{0.1106} & \bb{0.0971} & \nb{0.0823} & \nb{0.0827} & \bb{0.0762} \\ & 120 & \nb{0.0988} & \nb{0.0990} & \bb{0.0853} & \nb{0.0710} & \nb{0.0715} & \bb{0.0646} \\ earn vs.~all & 200 & \nb{0.0808} & \nb{0.0810} & \bb{0.0660} & \nb{0.0535} & \nb{0.0538} & \bb{0.0480} \\ & 400 & \nb{0.0578} & \nb{0.0578} & \bb{0.0456} & \nb{0.0404} & \nb{0.0408} & \bb{0.0358} \\ & 600 & \nb{0.0465} & \nb{0.0464} & \bb{0.0367} & \nb{0.0323} & \nb{0.0325} & \bb{0.0290} \\ \hline \hline & 80 & \nb{0.1126} & \nb{0.1125} & \bb{0.0846} & \nb{0.0788} & \nb{0.0785} & \bb{0.0667} \\ & 120 & \nb{0.0886} & \nb{0.0885} & \bb{0.0697} & \nb{0.0632} & \nb{0.0632} & \bb{0.0534} \\ acq vs.~all & 200 & \nb{0.0678} & \nb{0.0676} & \bb{0.0562} & \nb{0.0499} & \nb{0.0500} & \bb{0.0441} \\ & 400 & \nb{0.0506} & \nb{0.0503} & \bb{0.0419} & \nb{0.0370} & \nb{0.0369} & \bb{0.0335} \\ & 600 & \nb{0.0439} & \nb{0.0435} & \bb{0.0363} & \nb{0.0318} & \nb{0.0316} & \bb{0.0301} \\ \hline \hline & 80 & \nb{0.1201} & \nb{0.1198} & \bb{0.0758} & \nb{0.0676} & \nb{0.0669} & \bb{0.0647\str} \\ & 120 & \nb{0.0986} & \nb{0.0979} & \bb{0.0639} & \nb{0.0557} & \nb{0.0545} & \bb{0.0531\str} \\ moneyFx vs.~all & 200 & \nb{0.0814} & \nb{0.0811} & \bb{0.0544} & \nb{0.0485} & \nb{0.0472} & \bb{0.0438} \\ & 400 & \nb{0.0578} & \nb{0.0567} & \bb{0.0416} & \nb{0.0427} & \nb{0.0418} & \bb{0.0392} \\ & 600 & \nb{0.0478} & \nb{0.0467} & \bb{0.0375} & \nb{0.0391} & \nb{0.0385} & \bb{0.0369\str} \\ \hline \hline & 80 & \nb{0.1443} & \nb{0.1440} & \bb{0.0925} & \nb{0.0536} & \bb{0.0518\str} & \nb{0.0595} \\ & 120 & \nb{0.1101} & \nb{0.1097} & \bb{0.0717} & \nb{0.0476} & \bb{0.0467\str} & \nb{0.0494} \\ grain vs.~all & 200 & \nb{0.0793} & \nb{0.0786} & \bb{0.0576} & \nb{0.0430} & \bb{0.0420\str} & \nb{0.0440} \\ & 400 & \nb{0.0590} & \nb{0.0573} & \bb{0.0450} & \nb{0.0349} & \bb{0.0340\str} & \nb{0.0365} \\ & 600 & \nb{0.0517} & \nb{0.0497} & \bb{0.0401} & \nb{0.0290} & \bb{0.0284\str} & \nb{0.0306} \\ \hline \hline & 80 & \nb{0.1396} & \nb{0.1396} & \bb{0.0865} & \nb{0.0502} & \bb{0.0485\str} & \nb{0.0524} \\ & 120 & \nb{0.0961} & \nb{0.0953} & \bb{0.0542} & \nb{0.0446} & \bb{0.0425\str} & \nb{0.0428} \\ crude vs.~all & 200 & \nb{0.0624} & \nb{0.0613} & \bb{0.0414} & \nb{0.0388} & \nb{0.0373} & \bb{0.0345\str} \\ & 400 & \nb{0.0409} & \nb{0.0403} & \bb{0.0325} & \nb{0.0345} & \nb{0.0337} & \bb{0.0297} \\ & 600 & \nb{0.0379} & \nb{0.0362} & \bb{0.0299} & \nb{0.0292} & \nb{0.0284} & \bb{0.0264\str} \\ \hline \end{tabular} \end{small} \end{center} \caption{Experimental results on the Reuters corpus, using SVMs for linear, Gaussian, and multinomial diffusion kernels. The left columns use tf representation and the right columns use tf-idf representation. The error rates shown are averages obtained using 20-fold cross validation. The best performance for each training set size $L$ is shown in boldface. An asterisk (*) indicates that the difference is not statistically significant according to the paired $t$ test at the 0.05 level.} \label{tab:reuters1} \end{table} \begin{table} \renewcommand{\arraystretch}{1.1} \begin{center} \begin{small} \begin{tabular}{|c|r||ccc||ccc|} \multicolumn{2}{c}{} & \multicolumn{3}{c}{tf Representation} & \multicolumn{3}{c}{tf-idf Representation} \\[1pt] \hline Task & $\;L\;$ & Linear & Gaussian & Diffusion & Linear & Gaussian & Diffusion \\ % Reuters 1 vs earn \hline & 40 & \nb{0.1043} & \nb{0.1043} & \bb{0.1021\str} & \nb{0.0829} & \nb{0.0831} & \bb{0.0814\str} \\ & 80 & \nb{0.0902} & \nb{0.0902} & \bb{0.0856\str} & \nb{0.0764} & \nb{0.0767} & \bb{0.0730\str} \\ acq vs.~earn & 120 & \nb{0.0795} & \nb{0.0796} & \bb{0.0715} & \nb{0.0626} & \nb{0.0628} & \bb{0.0562} \\ & 200 & \nb{0.0599} & \nb{0.0599} & \bb{0.0497} & \nb{0.0509} & \nb{0.0511} & \bb{0.0431} \\ & 400 & \nb{0.0417} & \nb{0.0417} & \bb{0.0340} & \nb{0.0336} & \nb{0.0337} & \bb{0.0294} \\ \hline \hline & 40 & \nb{0.0759} & \nb{0.0758} & \bb{0.0474} & \nb{0.0451} & \nb{0.0451} & \bb{0.0372\str} \\ & 80 & \nb{0.0442} & \nb{0.0443} & \bb{0.0238} & \nb{0.0246} & \nb{0.0246} & \bb{0.0177} \\ moneyFx vs.~earn & 120 & \nb{0.0313} & \nb{0.0311} & \bb{0.0160} & \nb{0.0179} & \nb{0.0179} & \bb{0.0120} \\ & 200 & \nb{0.0244} & \nb{0.0237} & \bb{0.0118} & \nb{0.0113} & \nb{0.0113} & \bb{0.0080} \\ & 400 & \nb{0.0144} & \nb{0.0142} & \bb{0.0079} & \nb{0.0080} & \nb{0.0079} & \bb{0.0062} \\ \hline \hline & 40 & \nb{0.0969} & \nb{0.0970} & \bb{0.0543} & \nb{0.0365} & \nb{0.0366} & \bb{0.0336\str} \\ & 80 & \nb{0.0593} & \nb{0.0594} & \bb{0.0275} & \nb{0.0231} & \nb{0.0231} & \bb{0.0201\str} \\ grain vs.~earn & 120 & \nb{0.0379} & \nb{0.0377} & \bb{0.0158} & \nb{0.0147} & \nb{0.0147} & \bb{0.0114\str} \\ & 200 & \nb{0.0221} & \nb{0.0219} & \bb{0.0091} & \nb{0.0082} & \nb{0.0081} & \bb{0.0069\str} \\ & 400 & \nb{0.0107} & \nb{0.0105} & \bb{0.0060} & \nb{0.0037} & \nb{0.0037} & \bb{0.0037\str} \\ \hline \hline & 40 & \nb{0.1108} & \nb{0.1107} & \bb{0.0950} & \bb{0.0583\str} & \nb{0.0586} & \nb{0.0590} \\ & 80 & \nb{0.0759} & \nb{0.0757} & \bb{0.0552} & \nb{0.0376} & \nb{0.0377} & \bb{0.0366\str} \\ crude vs.~earn & 120 & \nb{0.0608} & \nb{0.0607} & \bb{0.0415} & \nb{0.0276} & \bb{0.0276\str} & \nb{0.0284} \\ & 200 & \nb{0.0410} & \nb{0.0411} & \bb{0.0267} & \bb{0.0218\str} & \nb{0.0218} & \nb{0.0225} \\ & 400 & \nb{0.0261} & \nb{0.0257} & \bb{0.0194} & \nb{0.0176} & \bb{0.0171\str} & \nb{0.0181} \\ \hline \end{tabular} \end{small} \end{center} \caption{Experimental results on the Reuters corpus, using support vector machines for linear, Gaussian, and multinomial diffusion kernels. The left columns use tf representation and the right columns use tf-idf representation. The error rates shown are averages obtained using 20-fold cross validation. The best performance for each training set size $L$ is shown in boldface. An asterisk (*) indicates that the difference is not statistically significant according to the paired $t$ test at the 0.05 level.} \label{tab:reuters2} \end{table} \begin{table} \renewcommand{\arraystretch}{1.1} \begin{small} \begin{center} \begin{tabular}{|l|c|c|c|}\hline Category&Linear&RBF&Diffusion \\ \hline \hline \texttt{earn} &0.01159& 0.01159& {\bf 0.01026}\\ \texttt{acq} &0.01854& 0.01854& {\bf 0.01788}\\ \texttt{money-fx}&0.02418& 0.02451& {\bf 0.02219}\\ \texttt{grain} &0.01391& 0.01391& {\bf 0.01060}\\ \texttt{crude} &0.01755& 0.01656& {\bf 0.01490}\\ \texttt{trade} &0.01722& {\bf 0.01656} & 0.01689\\ \texttt{interest}&0.01854& 0.01854& {\bf 0.01689}\\ \texttt{ship} &0.01324& 0.01324& {\bf 0.01225}\\ \texttt{wheat} &0.00894& 0.00794& {\bf 0.00629}\\ \texttt{corn} &0.00794& 0.00794& {\bf 0.00563}\\ \hline \end{tabular} \end{center} \caption{Test set error rates for the Reuters top 10 classes using tf features. The train and test sets were created using the Mod-Apt split.} \label{tab:modapte-split1} \end{small} \end{table} \begin{table}[t] \renewcommand{\arraystretch}{1.1} \begin{small} \begin{center} \begin{tabular}{|l|c|c|c|c|}\hline Category&Linear&RBF&Diffusion & $\pm$ \\ \hline \hline \texttt{earn} & 0.9781 & 0.9781 & {\bf 0.9808 } & $-$\\ \texttt{acq} & 0.9626 & 0.9626 & {\bf 0.9660 } & $+$ \\ \texttt{money-fx} & 0.8254 & 0.8245 & {\bf 0.8320 } & $+$ \\ \texttt{grain} & 0.8836 & 0.8844 & {\bf 0.9048 } & $-$\\ \texttt{crude} & 0.8615 & 0.8763 & {\bf 0.8889 } & $+$ \\ \texttt{trade} & 0.7706 & 0.7797 & {\bf 0.8050 } & $+$\\ \texttt{interest} & {\bf 0.8263 } & {\bf 0.8263 } & 0.8221 & $+$ \\ \texttt{ship} & 0.8306 & 0.8404 & {\bf 0.8827 } & $+$ \\ \texttt{wheat} & 0.8613 & 0.8613 & {\bf 0.8844 } & $-$\\ \texttt{corn} & 0.8727 & 0.8727 & {\bf 0.9310 } & $+$ \\ \hline \end{tabular} \end{center} \caption{F1 measure for the Reuters top 10 classes using tf features. The train and test sets were created using the Mod-Apte split. The last column compares the presented results with the published results of \cite{Zhang:01}, with a $+$ indicating the diffusion kernel F1 measure is greater than the result published in \cite{Zhang:01} for this task.} \label{tab:modapte-split2} \end{small} \end{table} \paragraph{} Figure~\ref{fig:experiments4} and Figure~\ref{fig:experiments5} show representative results for the second type of classification task, where the goal is to discriminate between two specific classes. In the case of the WebKB data the results are shown for \texttt{course} vs.~\texttt{student}. In the case of the Reuters data the results are shown for \texttt{moneyFx} vs.~\texttt{earn} and \texttt{grain} vs.~\texttt{earn}. Again, the results for the other classes are qualitatively similar; the numerical results are summarized in Tables~\ref{tab:webkb2} and~\ref{tab:reuters2}. \paragraph{} In these figures, the leftmost plots show the performance of tf features while the rightmost plots show the performance of tf-idf features. As mentioned above, in the case of the diffusion kernel we use $L_1$ normalization to give a valid embedding into the probability simplex, while for the linear and Gaussian kernels we use $L_2$ normalization, which works better empirically than $L_1$ for these kernels. The curves show the test set error rates averaged over 20 iterations of cross validation as a function of the training set size. The error bars represent one standard deviation. For both the Gaussian and diffusion kernels, we test scale parameters ($\sqrt{2}\sigma$ for the Gaussian kernel and $2t^{1/2}$ for the diffusion kernel) in the set $\{0.5,1,2,3,4,5,7,10\}$. The results reported are for the best parameter value in that range. \paragraph{} We also performed experiments with the popular Mod-Apte train and test split for the top 10 categories of the Reuters collection. For this split, the training set has about 7000 documents and is highly biased towards negative documents. We report in Table~\ref{tab:modapte-split1} the test set accuracies for the tf representation. For the tf-idf representation, the difference between the different kernels is not statistically significant for this amount of training and test data. The provided train set is more than enough to achieve outstanding performance with all kernels used, and the absence of cross validation data makes the results too noisy for interpretation. \paragraph{} In Table~\ref{tab:modapte-split2} we report the F1 measure rather than accuracy, since this measure is commonly used in text classification. The last column of the table compares the presented results with the published results of \emcite{Zhang:01}, with a $+$ indicating the diffusion kernel F1 measure is greater than the result published by \emcite{Zhang:01} for this task. \paragraph{} Our results are consistent with previous experiments in text classification using SVMs, which have observed that the linear and Gaussian kernels result in very similar performance \cite{Joachims:01}. However the multinomial diffusion kernel significantly outperforms the linear and Gaussian kernels for the tf representation, achieving significantly lower error rate than the other kernels. For the tf-idf representation, the diffusion kernel consistently outperforms the other kernels for the WebKb data and usually outperforms the linear and Gaussian kernels for the Reuters data. The Reuters data is a much larger collection than WebKB, and the document frequency statistics, which are the basis for the inverse document frequency weighting in the tf-idf representation, are evidently much more effective on this collection. It is notable, however, that the multinomial information diffusion kernel achieves at least as high an accuracy without the use of any heuristic term weighting scheme. These results offer evidence that the use of multinomial geometry is both theoretically motivated and practically effective for document classification. %--------------------------------------------------------------- \subsection{Experimental Results for Gaussian Embedding} \label{sec:experimentsGaussian} In this section we report experiments on synthetic data demonstrating the applicability of the heat kernel on the hyperbolic space $\H^n$ that corresponds to the manifold of spherical normal distributions. \paragraph{} Embedding data points in the hyperbolic space is more complicated than the multinomial case. Recall that in the multinomial case, we embedded points by computing the maximum likelihood estimator the the data. A similar method would fail for $\H^n$ embedding since all the data points will be mapped to normal distributions with variance 0, which will result in a degenerate geometry that is equivalent to the Euclidean one (see Section~\ref{sec:geomSphericalNormal}). \paragraph{} A more realistic embedding can be achieved by sampling from the posterior of a Dirichlet Process Mixture Model (DPMM) \cite{Ferguson73,Blackwell73,Antoniak74}. A Dirichlet Process Mixture model, based on the spherical Normal distribution, associates with the data $x_1,\ldots,x_m\in\R^n$ a posterior $p(\theta_1,\ldots,\theta_m|x_1,\ldots,x_m)$ where $\theta_i\in\H^{n+1}$. Instead of going into the description and properties of Dirichlet Process Mixture Model we refer the interested reader to the references above, and to \cite{Escobar95,Neal2000} for relevant Monte Carlo approximations. \paragraph{} Obtaining $T$ samples from the DPMM posterior $\{\theta_i^{(t)}\}_{i=1,t=1}^{m,T}$ we can measure the similarity between points $x_i,x_j$ as \begin{align} \tilde{K}_t(x_i,x_j) = \frac{1}{T}\,\,\sum_{t=1}^T K_t\left(\theta_i^{(t)},\theta_j^{(t)}\right)\end{align} where $K_t$ is the heat kernel on the hyperbolic space given by equations \eqref{eq:hyper_a}--\eqref{eq:hyper_b}. The above definition of $\tilde K$ has the interpretation of being approximately the mean of the heat kernel -- which is now a random variable under the DPMM posterior $p(\theta|x)$. The details of Gibbs sampling from the posterior of a spherical Normal based DPMM is given in Appendix~\ref{sec:appDPMM}. \paragraph{} Some intuition may be provided by Figure~\ref{fig:DPMMvis}\footnote{This Figure is better displayed in color.}. A sample from the posterior of a DPMM based on data from two Gaussians $N((-1,-1)^{\top},I)$ (solid blue dots) and $N((1,1)^{\top},I)$ (hollow red dots). The sample is illustrated by the circles that represent one standard deviation centered around the mean. Blue dotted circles represent a parameter associated with a point from $N((-1,-1)^{\top},I)$, red dashed circles represent a parameter associated with a point from $N((1,1)^{\top},I)$ and solid black circles represent parameters associated with points from both Gaussians. \begin{figure} \begin{center} \includegraphics[scale=0.6]{diffusionFigures/DPMMvis.eps} \end{center} \caption{A sample from the posterior of a DPMM based on data from two Gaussians $N((-1,-1)^{\top},I)$ (solid blue dots) and $N((1,1)^{\top},I)$ (hollow red dots). The sample is illustrated by the circles that represent one standard deviation centered around the mean. Blue dotted circles represent a parameter associated with a point from $N((-1,-1)^{\top},I)$, red dashed circles represent a parameter associated with a point from $N((1,1)^{\top},I)$ and solid black circles represent parameters associated with points from both Gaussians.} \label{fig:DPMMvis} \end{figure} \paragraph{} After generating data from two significantly overlapping Gaussians $N((-2,-2)^{\top},5I),N((2,2)^{\top},5I)$, and sampling 20 samples from the posterior we compared the performance of a standard RBF kernel and $\tilde K$. The test set accuracy for SVM, after 20-fold cross validation, as a function of the train set size is displayed in Figure~\ref{fig:DPMMacc}. The mean hyperbolic heat kernel outperforms the RBF kernel consistently. \begin{figure} \begin{center} \includegraphics[scale=0.6]{diffusionFigures/DPMMaccPlot.eps} \end{center} \caption{Test set error rate for SVM based on standard RBF kernel and the mean hyperbolic heat kernel $\tilde{K}$ as a function of the train set size, after 20-fold cross validation. The data was generated from two significantly overlapping Gaussians $N((-2,-2)^{\top},5I),N((2,2)^{\top},5I)$.} \label{fig:DPMMacc} \end{figure} %---------------------------------------------------------------- \subsection{Discussion} \label{sec:discuss} In this section we introduced a family of kernels that is intimately based on the geometry of the Riemannian manifold associated with a statistical family through the Fisher information metric. The metric is canonical in the sense that it is uniquely determined by requirements of invariance \cite{Chentsov1982}, and moreover, the choice of the heat kernel is natural because it effectively encodes a great deal of geometric information about the manifold. While the geometric perspective in statistics has most often led to reformulations of results that can be viewed more traditionally, the kernel methods developed here clearly depend crucially on the geometry of statistical families. \paragraph{} The main application of these ideas has been to develop the multinomial diffusion kernel. Our experimental results indicate that the resulting diffusion kernel is indeed effective for text classification using support vector machine classifiers, and can lead to significant improvements in accuracy compared with the use of linear or Gaussian kernels, which have been the standard for this application. The results of Section~\ref{sec:experimentsDiff} are notable since accuracies better or comparable to those obtained using heuristic weighting schemes such as tf-idf are achieved directly through the geometric approach. In part, this can be attributed to the role of the Fisher information metric; because of the square root in the embedding into the sphere, terms that are infrequent in a document are effectively up-weighted, and such terms are typically rare in the document collection overall. The primary degree of freedom in the use of information diffusion kernels lies in the specification of the mapping of data to model parameters. For the multinomial, we have used the maximum likelihood mapping. The use of other model families and mappings remains an interesting direction to explore. \paragraph{} While kernel methods generally are ``model free,'' and do not make distributional assumptions about the data that the learning algorithm is applied to, statistical models offer many advantages, and thus it is attractive to explore methods that combine data models and purely discriminative methods. Our approach combines parametric statistical modeling with non-parametric discriminative learning, guided by geometric considerations. In these aspects it is related to the methods proposed by~\emcite{Jaakkola1998}. However, the kernels proposed in the current section differ significantly from the Fisher kernel of \emcite{Jaakkola1998}. In particular, the latter is based on the score $\grad \theta \log p(X\given\hat\theta)$ at a single point $\hat\theta$ in parameter space. In the case of an exponential family model it is given by a covariance $K_F(x,x') = \sum_i \left(x_i - E_{\hat\theta}[X_i]\right) \left(x_i' - E_{\hat\theta}[X_i]\right)$; this covariance is then heuristically exponentiated. In contrast, information diffusion kernels are based on the full geometry of the statistical family, and yet are also invariant under re-parameterizations of the family. In other conceptually related work, \emcite{Belkin:2003} suggest measuring distances on the data graph to approximate the underlying manifold structure of the data. In this case the underlying geometry is inherited from the embedding Euclidean space rather than the Fisher geometry. \paragraph{} While information diffusion kernels are very general, they will be difficult to compute in many cases -- explicit formulas such as equations (\ref{eq:hyper_a}--\ref{eq:hyper_b}) for hyperbolic space are rare. To approximate an information diffusion kernel it may be attractive to use the parametrices and geodesic distance between points, as we have done for the multinomial. In cases where the distance itself is difficult to compute exactly, a compromise may be to approximate the geodesic distance between nearby points in terms of the Kullback-Leibler divergence. In effect, this approximation is already incorporated into the kernels recently proposed by \emcite{Moreno:03} for multimedia applications, which have the form $K(\theta,\theta') \propto \exp(-\alpha D(\theta,\theta')) \approx \exp(-2\alpha d^2(\theta,\theta'))$, and so can be viewed in terms of the leading order approximation to the heat kernel. The results of \emcite{Moreno:03} are suggestive that diffusion kernels may be attractive not only for multinomial geometry, but also for much more complex statistical families. %---------------------------------------------------------------- %---------------------------------------------------------------- %---------------------------------------------------------------- %---------------------------------------------------------------- \section{Hyperplane Margin Classifiers} \label{sec:hyperplane} Linear classifiers are a mainstay of machine learning algorithms, forming the basis for techniques such as the perceptron, logistic regression, boosting, and support vector machines. A linear classifier, parameterized by a vector $w\in\mathbb{R}^n$, classifies examples according to the decision rule $\hat y(x) = \text{sign} \left( \sum_i w_i \phi_i(x) \right) = \text{sign}(\langle w,x\rangle) \in \{-1,+1\}$, following the common practice of identifying $x$ with the feature vector $\phi(x)$. The differences between different linear classifiers lie in the criteria and algorithms used for selecting the parameter vector $w$ based on a training set. \paragraph{} Geometrically, the decision surface of a linear classifier is formed by a hyperplane or linear subspace in $n$-dimensional Euclidean space, $\{x\in\mathbb{R}^n: \langle x,w \rangle =0\}$ where $\langle\cdot,\cdot\rangle$ denotes the Euclidean inner product. (In both the algebraic and geometric formulations, a bias term is sometimes added; we prefer to absorb the bias into the notation given by the inner product, by setting $x_n=1$ for all $x$.) The linearity assumption made by such classifiers can be justified on purely computational grounds; linear classifiers are generally easy to train, and the linear form is simple to analyze and compute. \paragraph{} Modern learning theory emphasizes the tension between fitting the training data well and the more desirable goal of achieving good generalization. A common practice is to choose a model that fits the data closely, but from a restricted class of models. The model class needs to be sufficiently rich to allow the choice of a good hypothesis, yet not so expressive that the selected model is likely to overfit the data. Hyperplane classifiers are attractive for balancing these two goals. Indeed, linear hyperplanes are a rather restricted set of models, but they enjoy many unique properties. For example, given two points $x,y\in\mathbb{R}^n$, the set of points equidistant from $x$ and $y$ is a hyperplane; this lies behind the intuition that a hyperplane is the correct geometric shape for separating sets of points. Similarly, a hyperplane is the best decision boundary to separate two Gaussian distributions of equal covariance. Another distinguishing property is that a hyperplane in $\mathbb{R}^n$ is isometric to $\mathbb{R}^{n-1}$, and can therefore be thought of as a reduced dimension version of the original feature space. Finally, a linear hyperplane is the union of straight lines, which are distance minimizing curves, or geodesics, in Euclidean geometry. \paragraph{} However, a fundamental assumption is implicitly associated with linear classifiers, since they are based crucially on the use of the Euclidean geometry of $\mathbb{R}^n$. If the data or features at hand lack a Euclidean structure, the arguments above for linear classifiers break down; arguably, there is lack of Euclidean geometry for the feature vectors in most applications. This section studies analogues of linear hyperplanes as a means of obtaining simple, yet effective classifiers when the data can be represented in terms of a natural geometric structure that is only locally Euclidean. This is the case for categorical data that is represented in terms of multinomial models, for which the associated geometry is spherical. \paragraph{} Because of the complexity of the notion of linearity in general Riemannian spaces, we focus our attention on the multinomial manifold, which permits a relatively simple analysis. Hyperplanes in multinomial manifold is discussed in Section \ref{sec:snp}. The construction and training of margin based models is discussed in Section \ref{sec:logReg}, with an emphasis on spherical logistic regression. A brief examination of linear hyperplanes in general Riemannian manifolds appears in Section~\ref{sec:general} followed by experimental results for text classification given in Section~\ref{sec:exp}. Concluding remarks are made in Section~\ref{sec:sum}. %-------------------------------------------------------------------------- \subsection{Hyperplanes and Margins on $\Sn$} \label{sec:sn} This section generalizes the notion of linear hyperplanes and margins to the $n$-sphere $\Sn=\{x\in\mathbb{R}^{n+1}:\sum_i x_i^2=1\}$. A similar treatment on the positive $n$-sphere $\Snp$ is more complicated, and is postponed to the next section. In the remainder of this section we denote points on $\Pn,\Sn$ or $\Snp$ as vectors in $\mathbb{R}^{n+1}$ using the standard basis of the embedding space. The notation $\langle \cdot,\cdot\rangle$ and $\|\cdot\|$ will be used for the Euclidean inner product and norm. A hyperplane on $\Sn$ is defined as $H_u=\Sn\cap E_{u}$ where $E_u$ is an $n$-dimensional linear subspace of $\mathbb{R}^{n+1}$ associated with the normal vector $u$. We occasionally need to refer to the unit normal vector (according to the Euclidean norm) and denote it by $\hat{u}$. $H_u$ is an $n-1$ dimensional submanifold of $\Sn$ which is isometric to $\mathbb{S}^{n-1}$ \cite{Bridson1999}. Using the common notion of the distance of a point from a set $d(x,S) = \inf_{y\in S} d(x,y)$ we make the following definitions. \begin{defn} Let $X$ be a metric space. A {\sl decision boundary} is a subset of $X$ that separates $X$ into two connected components. The {\sl margin\/} of $x$ with respect to a decision boundary $H$ is $d(x,H)=\inf_{y\in H} d(x,y)$. \end{defn} Note that this definition reduces to the common definition of margin for Euclidean geometry and affine hyperplanes. \paragraph{} In contrast to \emcite{Gous1998}, our submanifolds are intersections of the sphere with linear subspaces, not affine sets. One motivation for the above definition of hyperplane as the correct generalization of a Euclidean hyperplane is that $H_u$ is the set of points equidistant from $x,y\in\Sn$ in the spherical metric. Further motivation is given in Section~\ref{sec:general}. \paragraph{} Before we can obtain a closed form expression for margins on $\Sn$ we need the following definitions. \begin{defn} Given a point $x\in\mathbb{R}^{n+1}$, we define its {\sl reflection} with respect to $E_u$ as \[r_u(x)=x-2\langle x,\hat u\rangle \hat u.\] \end{defn} Note that if $x\in\Sn$ then $r_u(x)\in\Sn$ as well, since $\norm{r_u(x)}^2 = \norm{x}^2-4\langle x,\hat u\rangle^2+4\langle x,\hat u\rangle^2 =1.$ \paragraph{} \begin{defn} The {\sl projection\/} of $x\in\Sn\setminus\{\hat u\}$ on $H_u$ is defined to be \[p_u(x)=\frac{x-\langle x,\hat u\rangle \hat u}{\sqrt{1-\langle x,\hat u\rangle^2}}.\] \end{defn} Note that $p_u(x)\in H_u$, since $\norm{p_u(x)}=1$ and $\langle x-\langle x,\hat u\rangle \hat u,\hat u\rangle =\langle x,\hat u\rangle - \langle x,\hat u\rangle \norm{\hat u}^2=0.$ The term projection is justified by the following proposition. \paragraph{} \begin{prop} \label{prop:sphereProj} Let $x\in\Sn\setminus (H_u\cup\{\hat u\})$. Then \begin{align*} &(a)\quad d(x,q)=d(r_u(x),q) \quad \forall q\in H_u\\ &(b)\quad d(x,p_u(x))= \arccos\left(\sqrt{1-\langle x,\hat u\rangle^2}\right)\\ %&\qquad\qquad =\frac{1}{2}\arccos(1-2\langle x,\hat u\rangle^2) \\ &(c)\quad d(x,H_u) = d(x,p_u(x)). \end{align*} \end{prop} \begin{proof} Since $q\in H_u$, \begin{align*} \cos d(r_u(x),q) &= \langle x-2\langle x,\hat u\rangle \hat u,q\rangle \\ & =\langle x,q\rangle-2\langle x,\hat u\rangle \langle\hat u,q\rangle\\ &=\langle x,q\rangle = \cos d(x,q) \end{align*} and $(a)$ follows. Assertion $(b)$ follows from \begin{align*} \cos d(x,p_u(x)) &= \left\langle x,\frac{x-\langle x,\hat u\rangle \hat u}{\sqrt{1-\langle x,\hat u\rangle^2}}\right\rangle =\frac{1-\langle x,\hat u\rangle^2}{\sqrt{1-\langle x,\hat u\rangle^2}}. \end{align*} Finally, to prove $(c)$ note that by the identity $\cos 2\theta=2\cos^2\theta-1$, \begin{align*} \cos(2d(x,p_u(x)))&=2\cos^2(d(x,p_u(x)))-1\\ &=1-2\langle x,\hat u\rangle^2 =\cos(d(x,r_u(x))) \end{align*} and hence $d(x,p_u(x))=\frac{1}{2} d(x,r_u(x))$. The distance $d(x,q), q\in H_u$ cannot be any smaller than $d(x,p_u(x))$ since this would result in a path from $x$ to $r_u(x)$ of length shorter than the geodesic $d(x,r_u(x))$. \end{proof} \paragraph{} Parts $(b)$ and $(c)$ of Proposition~\ref{prop:sphereProj} provide a closed form expression for the $\Sn$ margin analogous to the Euclidean unsigned margin $|\langle x,\hat u\rangle|$. Similarly, the $\Sn$ analogue of the Euclidean signed margin $y\langle \hat u,x\rangle$ is \[ y\frac{\langle x,\hat u\rangle}{|\langle x,\hat u\rangle|} \,\,\,\arccos\left(\sqrt{1-\langle x,\hat u\rangle^2}\right).\] A plot of the signed margin as a function of $\langle x,\hat u\rangle$ and a geometric interpretation of the spherical margin appear in Figure~\ref{fig:snMargin}. \begin{figure} \begin{center} \includegraphics[scale=0.45]{hyperplaneFigures/acosPlot.eps} \hspace{.45in} { \psfrag{x}{{\scriptsize $x$}} \psfrag{u}{{\scriptsize $\hat u$}} \psfrag{p(x)}{{\scriptsize $p_u(x)$}} \psfrag{x-2u}{{\scriptsize $r_u(x)$}} \psfrag{}{{\scriptsize $\langle x,\hat u\rangle$}} \psfrag{n}{{\scriptsize $d(x,p_u(x))=\frac{1}{2}d(x,r_u(x))$}} \includegraphics[scale=0.37]{hyperplaneFigures/sphere.eps} } \caption{ The signed margin $\text{sign}(\langle x,\hat y\rangle)d(x,H_u)$ as a function of $\langle x,\hat{u}\rangle$, which lies in the interval $[-\frac{\pi}{2},\frac{\pi}{2}]$ (left) and a geometric interpretation of the spherical margin (right).} \label{fig:snMargin} \end{center} \end{figure} %----------------------------------------------------------------------- \subsection{Hyperplanes and Margins on $\Snp$} \label{sec:snp} A hyperplane on the positive $n$-sphere $\Snp$ is defined as $H_{u+}=E_u\cap \Snp$, assuming it is non-empty. This definition leads to a margin concept $d(x,H_{u+})$ different from the $\Sn$ margin $d(x,H_u)$ since \begin{align*} d(x,H_{u+})&=\inf_{y\in E_u\cap \Snp} d(x,y) \\&\geq \inf_{y\in E_u\cap \Sn} d(x,y)=d(x,H_u). \end{align*} The above infimum is attained by the continuity of $d$ and compactness of $E_u\cap\Snp$ justifying the notation $q=\argmin_{y\in E_u\cap\Snp}d(x,y)$ as a point realizing the margin distance $d(x,H_{u+})$. \begin{figure} \begin{center} { \psfrag{x}{$x$} \psfrag{r}{$r$} \psfrag{q}{$q$} \psfrag{p(x)}{$p_u(x)$} \includegraphics[scale=0.33]{hyperplaneFigures/triangles.eps} }\end{center} \caption{The spherical law of cosines implies $d(r,x)\leq d(q,x)$.} \label{fig:triangles} \end{figure} \paragraph{} The following theorem will be useful in computing $d(x,H_{u+})$. For a proof see \emcite{Bridson1999} page 17. \begin{thrm} \emph{(\textbf{The Spherical Law of Cosines})}\\ Consider a spherical triangle with geodesic edges of lengths $a,b,c$, where $\gamma$ is the vertex angle opposite to edge $c$. Then \[ \cos c = \cos a\cos b + \sin a\sin b\cos \gamma.\] \end{thrm} \paragraph{} \paragraph{} We have the following corollaries of Proposition \ref{prop:sphereProj}. \begin{prop} If $x\in\Snp$ and $p_u(x)\in\Snp$ then \begin{align*} p_u(x)&=\argmin_{y\in\Snp\cap E_u} d(x,y)\\ d(x,H_u)&=d(x,H_{u+}) \end{align*} %$p_u(x)=\argmin_{y\in\Snp\cap E_u} d(x,y)$ and %$d(x,H_u)=d(x,H_{u+})$. \end{prop} \begin{proof} \label{prop:triangles} This follows immediately from the fact that $p_u(x)=\argmin_{y\in\Sn\cap E_u} d(x,y)$ and from $\Snp\cap E_u\subset \Sn\cap E_u$. \end{proof} \begin{prop} \label{prop:bdP1} For $x\in\Snp$ and $p_u(x)\not\in\Snp$ we have \[q=\argmin_{y\in\Snp\cap E_u} d(x,y) \in \partial \Snp\] where $\partial\Snp$ is the boundary of $\Snp$. \end{prop} \begin{proof} Assume that $q\not\in\partial\Snp$ and connect $q$ and $p_u(x)$ by a minimal geodesic $\alpha$. Since $p_u(x)\not\in\Snp$, the geodesic $\alpha$ intersects the boundary $\partial\Snp$ at a point $r$. Since $q,p_u(x)\in H_u$ and $H_u$ is geodesically convex, $\alpha \subset H_u$. Now, since $p_u(x)=\argmin_{y\in \alpha} d(y,x)$, the geodesic from $x$ to $p_u(x)$ and $\alpha$ intersect orthogonally (this is an elementary result in Riemannian geometry, e.g. \emcite{Lee1997} p. 113). Using the spherical law of cosines, applied to the spherical triangles $(q,x,p_u(x))$ and $(r,x,p_u(x))$ (see Figure \ref{fig:triangles}), we deduce that \begin{align*} \cos d(x,q) &= \cos d(q,p_u(x)) \cos d(x,p_u(x))\\ &\leq \cos d(r,p_u(x)) \cos d(x,p_u(x)) \\ &= \cos d(x,r) \end{align*} Hence $r$ is closer to $x$ than $q$. This contradicts the definition of $q$; thus $q$ can not lie in the interior of $\Snp$. \end{proof} \paragraph{} Before we proceed to compute $d(x,H_{u+})$ for $p_u(x)\not\in\Snp$ we define the following concepts. \begin{defn}\label{def:bd} The {\sl boundary} of $\Sn$ and $\Snp$ with respect to $A\subset\{1,\ldots,n+1\}$ is \begin{align*} \partial_A\Sn &= \Sn \cap \{x\in\mathbb{R}^{n+1}: \forall i\in A,\, x_i=0\}\cong \mathbb{S}^{n-|A|}\\ \partial_A\Snp &= \Snp \cap \{x\in\mathbb{R}^{n+1}: \forall i\in A,\, x_i=0\}\cong \mathbb{S}^{n-|A|}_{+} \end{align*} \end{defn} Note that if $A\subset A'$ then $\partial_{A'}\Snp\subset \partial_{A}\Snp$. We use the notation $\langle \cdot,\cdot \rangle_A$ and $\norm{\cdot}_A$ to refer to the Euclidean inner product and norm, where the summation is restricted to indices \emph{not} in $A$. \begin{defn} Given $x\in\Sn$ we define $x|_A\in\partial_A\Sn$ as \[(x|_A)_i = \begin{cases} 0 & i\in A \\ x_i/\norm{x}_{A} & i\not\in A\,. \end{cases}\] \end{defn} We abuse the notation by identifying $x|_A$ also with the corresponding point on $\mathbb{S}^{n-|A|}$ under the isometry $\partial_A \Sn\cong\mathbb{S}^{n-|A|}$ mentioned in Definition \ref{def:bd}. Note that if $x\in\Snp$ then $x|_{A}\in\partial_{A}\Snp$. The following proposition computes the $\Snp$ margin $d(x,H_{u+})$ given the boundary set of $q=\argmin_{y\in \Snp\cap E_u}d(y,x)$. \begin{prop} Let $\hat u\in\Rn1$ be a unit vector, $x\in\Snp$ and $q=\argmin_{y\in \Snp\cap E_u}d(y,x)\in \partial_A\Snp$ where $A$ is the (possibly empty) set $A = \{1\leq i\leq n+1 \,:\, q_i=0\}.$ Then \begin{align*} d(x,H_{u+}) = \arccos \left(\norm{x}_A \sqrt{1-\langle x|_A,\hat u|_A\rangle^2}\right) \end{align*} \end{prop} \begin{proof} If $p_u(x)\in\Snp$ then the proposition follows from earlier propositions and the fact that when $A=\emptyset$, $\norm{x}_A=\norm{x}=1$ and $v|_A=v$. We thus restrict our attention to the case of $A\neq\emptyset$. For all $I\subset \{1,\ldots,n+1\}$ we have \begin{eqnarray*} \argmin_{y\in \partial_I\Snp\cap E_u} d(x,y) &=&\argmax_{y\in \partial_I\Snp\cap E_u} \langle x,y\rangle \\ &=&\argmax_{y\in \partial_I\Snp\cap E_u} \langle x,y\rangle_I \\ &=& \argmax_{y\in \mathbb{S}^{n-|I|}_{+}\cap E_{u|I}} \norm{x}_I \langle x|_I,y\rangle \\ &=& \argmin_{y\in \mathbb{S}^{n-|I|}_{+}\cap E_{u|I}} d(x|_I,y)\,. \end{eqnarray*} It follows that \begin{align} \label{eq:compQ} q|_A= \argmin_{y\in \mathbb{S}^{n-|A|}_{+}\cap E_{u|A}} d(x|_A,y). \end{align} By Proposition~\ref{prop:bdP1} applied to $\mathbb{S}^{n-|A|}$ we have that since $q|_A$ lies in the interior of $\mathbb{S}^{n-|A|}$ then so does \begin{eqnarray*} p_{u|A}(x|_A) = \frac{x|_A-\langle x|_A,\hat u|_A\rangle \hat u|_A}{\sqrt{1-\langle x|_A,\hat u|_A\rangle^2}}\,, \quad x|_A, \hat u|_A \in\mathbb{S}^{n-|A|}_{+}. \end{eqnarray*} Using Proposition \ref{prop:sphereProj} applied to $\mathbb{S}^{n-|A|}$ we can compute $d(x,H_{u+})$ as \begin{align*} d(x, p_{u|A}(x|_A)) &= \arccos \left\langle x, \frac{x|_A-\langle x|_A,\hat u|_A\rangle \hat u|_A}{\sqrt{1-\langle x|_A,\hat u|_A\rangle^2}} \right\rangle \\ %= \arccos \frac{\langle x,x|_A-\langle x|_A,\hat u|_A\rangle \hat u|_A\rangle}{ % \sqrt{1-\langle x|_A,\hat u|_A\rangle^2}} \\ &= \arccos \frac{\norm{x}_A- \langle x|_A,\hat u|_A\rangle \langle x, \hat u|_A\rangle}{ \sqrt{1-\langle x|_A,\hat u|_A\rangle^2} }\\ %&= \arccos \frac{\norm{x}_A(1- % \langle x|_A,\hat u|_A\rangle^2)} % {\sqrt{1-\langle x|_A,\hat u|_A\rangle^2} % }\\ &= \arccos \left(\norm{x}_A \sqrt{1-\langle x|_A,\hat u|_A\rangle^2}\right). \end{align*} \end{proof} \paragraph{} In practice the boundary set $A$ of $q$ is not known. In our experiments we set $A= \{i:(p_u(x))_i \leq 0\}$; in numerical simulations in low dimensions, the true boundary never lies outside of this set. %----------------------------------------------------------------------- \subsection{Logistic Regression on the Multinomial Manifold} \label{sec:logReg} The logistic regression model $p(y\given x)=\frac{1}{Z}\exp(y\langle x,w\rangle)$, with $y\in\{-1,1\}$, assumes Euclidean geometry. It can be re-expressed as \begin{eqnarray*} p(y\given x\,; u) &\propto& \exp\left(y\norm{u}\langle x,\hat{u}\rangle\right)\\ &=& \exp\left(y\,\text{sign}(\langle x,\hat u\rangle)\,\theta d(x,H_u)\right) \end{eqnarray*} where $d$ is the Euclidean distance of $x$ from the hyperplane that corresponds to the normal vector $\hat{u}$, and where $\theta = \norm{u}$ is a parameter. \paragraph{} The generalization to spherical geometry involves simply changing the margin to reflect the appropriate geometry: \begin{eqnarray*} \lefteqn{p(y|x\,; \hat u,\theta) \propto} && \\ && \hskip-.3in \exp\left(y\,\text{sign}(\langle x,\hat u\rangle) \,\, \theta \arccos\left(\norm{x}_A\sqrt{1-\langle x|_A,\hat{u}|_A\rangle^2}\right)\right). \end{eqnarray*} Denoting $s_x=y\,\text{sign}(\langle x,\hat u\rangle)$, the log-likelihood of the example $(x,y)$ is \begin{align*} &\ell(\hat u,\theta\,; (x,y))\\ %&= \log %\frac{e^{s_x\theta\arccos\left(\|x\|_A\sqrt{1-\langle x|_A,\hat{u}|_A\rangle^2}\right)}}{ %e^{s_x\theta\arccos\left(\theta \|x\|_A\sqrt{1-\langle x|_A,\hat{u}|_A\rangle^2}\right)} %+e^{-s_x\theta\arccos\left(\norm{x}_A\sqrt{1-\langle x|_A,\hat{u}|_A\rangle^2}\right)}}\\ &=-\log \left(1+ e^{-2s_x\theta\arccos\left(\norm{x}_A\sqrt{1-\langle x|_A,\hat{u}|_A\rangle^2}\right)}\right). \end{align*} \paragraph{} We compute the derivatives of the log-likelihood in several steps, using the chain rule and the notation $z=\langle x|_A,\hat u|_A\rangle$. We have \begin{eqnarray*} {\frac{\partial\arccos\left(\norm{x}_A\sqrt{1-z^2}\right)}{\partial z} =} \frac{z\norm{x}_A}{\sqrt{1-\norm{x}_A^2(1-z^2)}\sqrt{1-z^2}} \end{eqnarray*} and hence \begin{eqnarray} \label{eq:llGrad} \frac{\partial\ell(\hat u,\theta\,; (x,y))}{\partial z}\;=\; \frac{2s_x\theta z\norm{x}_A/(1+e^{2s_x\theta\arccos\left(\norm{x}_A\sqrt{1-z^2}\right)})}{ \sqrt{1-\norm{x}_A^2(1-z^2)}\sqrt{1-z^2}}. \end{eqnarray} The log-likelihood derivative with respect to $\hat u_i$ is equation \eqref{eq:llGrad} times \begin{align*} \frac{\partial \langle x|_A,\hat u|_A\rangle}{\partial \hat u_i} &=\begin{cases} 0 & i\in A\\ \frac{ (x|_A)_i}{\norm{\hat u}_A} - \hat u_i \frac{\langle x|_A,\hat u|_A\rangle}{\norm{u}_A^2} & i\not\in A\,. \end{cases} \end{align*} The log-likelihood derivative with respect to $\theta$ is \begin{align*} \frac{\partial \ell(\hat u,\theta\,; (x,y))}{\partial\theta}& %\frac{e^{-2s_x\theta\arccos(\norm{x}_A\sqrt{1-z^2})} %}{1+e^{-2s_x\theta\arccos(\norm{x}_A\sqrt{1-z^2})}} %2s_x \arccos(\norm{x}_A\sqrt{1-z^2}) =\frac{2s_x \arccos(\norm{x}_A\sqrt{1-z^2})}{1+e^{2s_x\theta\arccos(\norm{x}_A\sqrt{1-z^2})}}. \end{align*} \paragraph{} Optimizing the log-likelihood with respect to $\hat{u}$ requires care. Following the gradient $\hat u^{(t+1)}=\hat u^{(t)}+\alpha \grad \ell(\hat u^{(t)})$ results in a non-normalized vector. Performing the above gradient descent step followed by normalization has the effect of moving along the sphere in a curve whose tangent vector at $\hat u^{(t)}$ is the projection of the gradient onto the tangent space $T_{\hat u^{(t)}}\Sn$. This is the technique used in the experiments described in Section~\ref{sec:exp}. \paragraph{} Note that the spherical logistic regression model has $n+1$ parameters in contrast to the $n+2$ parameters of Euclidean logistic regression. This is in accordance with the intuition that a hyperplane separating an $n$-dimensional manifold should have $n$ parameters. The extra parameter in the Euclidean logistic regression is an artifact of the embedding of the $n$-dimensional multinomial space, on which the data lies, into an $(n+1)$-dimensional Euclidean space. \paragraph{} The derivations and formulations above assume spherical data. If the data lies on the multinomial manifold, the isometry $\pi$ mentioned earlier has to precede these calculations. The net effect is that $x_i$ is replaced by $\sqrt{x_i}$ in the model equation, and in the log-likelihood and its derivatives. \paragraph{} Synthetic data experiments contrasting Euclidean logistic regression and spherical logistic regression on $\Snp$, as described in this section, are shown in Figure~\ref{fig:toy-data-hyperplane}. The leftmost column shows an example where both models give a similar solution. In general, however, as is the case in the other two columns, the two models yield significantly different decision boundaries. \begin{figure} \def\figsz{.37} \begin{center} \begin{tabular}{cc} \includegraphics[scale=\figsz]{hyperplaneFigures/toyDataPlots/toyData11c.eps}& \includegraphics[scale=\figsz]{hyperplaneFigures/toyDataPlots/toyData12c.eps}\\ \includegraphics[scale=\figsz]{hyperplaneFigures/toyDataPlots/toyData21c.eps}& \includegraphics[scale=\figsz]{hyperplaneFigures/toyDataPlots/toyData22c.eps}\\ \includegraphics[scale=\figsz]{hyperplaneFigures/toyDataPlots/toyData31c.eps}& \includegraphics[scale=\figsz]{hyperplaneFigures/toyDataPlots/toyData32c.eps} \end{tabular} \end{center} \caption{Experiments contrasting Euclidean logistic regression (left column) with multinomial logistic regression (right column) for several toy data sets in $\mathbb{P}^2$.}\label{fig:toy-data-hyperplane} \end{figure} %----------------------------------------------------------------------- \subsection{Hyperplanes in Riemannian Manifolds}\label{sec:general} The definition of hyperplanes in general Riemannian manifolds has two essential components. In addition to discriminating between two classes, hyperplanes should be regular in some sense with respect to the geometry. In Euclidean geometry, the two properties of discrimination and regularity coincide, as every affine subspace of dimension $n-1$ separates $\mathbb{R}^n$ into two regions. In general, however, these two properties do not necessarily coincide, and have to be considered separately. \paragraph{} The separation property implies that if $N$ is a hyperplane of $M$ then $M\setminus N$ has two connected components. Note that this property is topological and independent of the metric. The linearity property is generalized through the notion of auto-parallelism explained below. The following definitions and propositions are taken from \emcite{Spivak1975}, Volume 3. All the connections described below $\nabla$ are the metric connections inherited from the metric $g$. \begin{defn} Let $(\cM,g)$ be a Riemannian manifold and $\nabla$ the metric connection. A submanifold $N\subset M$ is {\sl auto-parallel\/} if parallel translation in $M$ along a curve $C\subset N$ takes vectors tangent to $N$ to vectors tangent to $N$. \end{defn} \begin{prop} A submanifold $N\subset M$ is auto-parallel if and only if \[X,Y\in T_pN\Rightarrow \nabla_XY\in T_pN.\] \end{prop} \paragraph{} \begin{defn} A submanifold $N$ of $M$ is {\sl totally geodesic at $p\in N$} if every geodesic $\gamma$ in $M$ with $\gamma(0)=p, \gamma'(0)\in T_pN$ remains in $N$ on some interval $(-\epsilon,\epsilon)$. The submanifold $N$ is said to be {\it totally geodesic\/} if it is totally geodesic at every point. \end{defn} As a consequence, we have that $N$ is totally geodesic if and only if every geodesic in $N$ is also a geodesic in $M$. \begin{prop} Let $N$ be a submanifold of $(M,\nabla)$. Then \begin{enumerate} \item If $N$ is auto-parallel in $M$ then $N$ is totally geodesic. \item If $M$ is totally geodesic and $\nabla$ is symmetric then $M$ is auto-parallel. \end{enumerate} \end{prop} Since the metric connection is symmetric, the last proposition gives a complete equivalence between auto-parallelism and totally geodesic submanifolds. \paragraph{} We can now define linear hyperplanes on Riemannian manifolds. \begin{defn} A {\sl linear decision boundary\/} $N$ in $M$ is an auto-parallel submanifold of $M$ such that $M\setminus N$ has two connected components. \end{defn} \paragraph{} Several observations are in order. First note that if $M$ is an $n$-dimensional manifold, the separability condition requires $N$ to be an $(n-1)$-dimensional submanifold. It is easy to see that every affine subspace of $\mathbb{R}^n$ is totally geodesic and hence auto-parallel. Conversely, since the metric connection is symmetric, every auto-parallel submanifold of Euclidean space that separates it is an affine subspace. As a result, we have that our generalization does indeed reduce to affine subspaces under Euclidean geometry. Similarly, the above definition reduces to spherical hyperplanes $H_u\cap \Sn$ and $H_u\cap\Snp$. Another example is the hyperbolic half plane $\mathbb{H}^2$ where the linear decision boundaries are half-circles whose centers lie on the $x$ axis. \paragraph{} Hyperplanes on $\Sn$ have the following additional nice properties. They are the set of equidistant points from $x,y\in\Sn$ (for some $x,y$), they are isometric to $\mathbb{S}^{n-1}$ and they are parameterized by $n$ parameters. These properties are particular to the sphere and do not hold in general \cite{Bridson1999}. %---------------------------------------------------------------------- \subsection{Experiments} \label{sec:exp} A natural embedding of text documents in the multinomial simplex is the $L_1$ normalized term-frequency (tf) representation \cite{Joachims2000} \begin{align*} \hat{\theta}(x) = \left(\frac{x_1}{\sum_i x_i},\ldots,\frac{x_{n+1}}{\sum_i x_i} \right). \end{align*} Using this embedding we compared the performance of spherical logistic regression with Euclidean logistic regression. Since Euclidean logistic regression often performs better with $L_2$ normalized tf representation, we included these results as well. \paragraph{} Experiments were conducted on both the Web-KB \cite{Craven1998} and the Reuters-21578 \cite{Lewis94} datasets. In the Web-KB dataset, the classification task that was tested was each of the classes faculty, course, project and student vs.~the rest. In the Reuters dataset, the task was each of the 8 most popular classes vs.~the rest. The test error rates as a function of randomly sampled training sets of different sizes are shown in Figures \ref{fig:hyperplaneExp1}-\ref{fig:hyperplaneExp3}. In both cases, the positive and negative example sets are equally distributed, and the results were averaged over a 20-fold cross validation with the error bars indicating one standard deviation. As mentioned in Section~\ref{sec:snp}, we assume that the boundary set of $q=\argmin_{y\in \Snp\cap E_u}d(y,x)$ is equal to $A = \{i:(p_u(x))_i \leq 0\}$. \paragraph{} The experiments show that the new linearity and margin concepts lead to more powerful classifiers than their Euclidean counterparts, which are commonly used in the literature regardless of the geometry of the data. \begin{figure} \begin{center} \begin{tabular}{cc} {\tiny Web-KB: faculty vs.~all} & {\tiny Web-KB: course vs.~all}\\ \includegraphics[scale=0.37]{hyperplaneFigures/webkb-results/webkbPlot1.eps}& \includegraphics[scale=0.37]{hyperplaneFigures/webkb-results/webkbPlot2.eps}\\ {\tiny Web-KB: project vs.~all} & {\tiny Web-KB: student vs.~all} \\ \includegraphics[scale=0.37]{hyperplaneFigures/webkb-results/webkbPlot3.eps}&%\\[8pt] \includegraphics[scale=0.37]{hyperplaneFigures/webkb-results/webkbPlot4.eps} \end{tabular} \end{center} \caption{ Test error accuracy of spherical logistic regression (solid), and linear logistic regression using tf representation with $L_1$ normalization (dashed) and $L_2$ normalization (dotted). The task is Web-KB binary ``one vs.~all'' classification, where the name of the topic is listed above the individual plots. Error bars represent one standard deviation over 20-fold cross validation for spherical logistic regression. The error bars of the other classifiers are of similar sizes and are omitted for clarity.} \label{fig:hyperplaneExp1} \end{figure} \begin{figure} \begin{center} \begin{tabular}{cc} {\tiny Reuters: earn vs.~all} & {\tiny Reuters: acq vs.~all} \\ \includegraphics[scale=0.4]{hyperplaneFigures/reuters-results/reutersPlot1.eps}& \includegraphics[scale=0.4]{hyperplaneFigures/reuters-results/reutersPlot2.eps}\\%[8pt] {\tiny Reuters: money-fx vs.~all} & {\tiny Reuters: grain vs.~all} \\ \includegraphics[scale=0.4]{hyperplaneFigures/reuters-results/reutersPlot3.eps}& \includegraphics[scale=0.4]{hyperplaneFigures/reuters-results/reutersPlot4.eps} \end{tabular} \end{center} \caption{ Test error accuracy of spherical logistic regression (solid), and linear logistic regression using tf representation with $L_1$ normalization (dashed) and $L_2$ normalization (dotted). The task is Reuters-21578 binary ``one vs.~all'' classification, where the name of the topic is listed above the individual plots. Error bars represent one standard deviation over 20-fold cross validation for spherical logistic regression. The error bars of the other classifiers are of similar sizes and are omitted for clarity.} \label{fig:hyperplaneExp2} \end{figure} \begin{figure} \begin{center} \begin{tabular}{cc} {\tiny Reuters: crude vs.~all} & {\tiny Reuters: trade vs.~all} \\ \includegraphics[scale=0.4]{hyperplaneFigures/reuters-results/reutersPlot5.eps}& %[8pt] \includegraphics[scale=0.4]{hyperplaneFigures/reuters-results/reutersPlot6.eps}\\ {\tiny Reuters: interest vs.~all} & {\tiny Reuters: ship vs.~all} \\ \includegraphics[scale=0.4]{hyperplaneFigures/reuters-results/reutersPlot7.eps}& \includegraphics[scale=0.4]{hyperplaneFigures/reuters-results/reutersPlot8.eps} \end{tabular} \end{center} \caption{ Test error accuracy of spherical logistic regression (solid), and linear logistic regression using tf representation with $L_1$ normalization (dashed) and $L_2$ normalization (dotted). The task is Reuters-21578 binary ``one vs.~all'' classification, where the name of the topic is listed above the individual plots. Error bars represent one standard deviation over 20-fold cross validation for spherical logistic regression. The error bars of the other classifiers are of similar sizes and are omitted for clarity.} \label{fig:hyperplaneExp3} \end{figure} %---------------------------------------------------------------- %---------------------------------------------------------------- %---------------------------------------------------------------- %---------------------------------------------------------------- \section{Metric Learning} \label{sec:metricLearning} Machine learning algorithms often require an embedding of data points into some space. Algorithms such as $k$-nearest neighbors and neural networks assume the embedding space to be $\mathbb{R}^n$ while SVM and other kernel methods embed the data in a Hilbert space through a kernel operation. Whatever the embedding space is, the notion of metric structure has to be carefully considered. The popular assumption of a Euclidean metric structure is often used without justification by data or modeling arguments. We argue that in the absence of direct evidence of Euclidean geometry, the metric structure should be inferred from data (if available). After obtaining the metric structure, it may be passed to a learning algorithm for use in tasks such as classification and clustering. \paragraph{} Several attempts have recently been made to learn the metric structure of the embedding space from a given data set. \emcite{Saul1997} use geometrical arguments to learn optimal paths connecting two points in a space. \emcite{Xing2003} learn a global metric structure. that is able to capture non-Euclidean geometry, but only in a restricted manner since the metric is constant throughout the space. \emcite{Lanckriet2002} learn a kernel matrix that represents similarities between all pairs of the supplied data points. While such an approach does learn the kernel structure from data, the resulting Gram matrix does not generalize to unseen points. \paragraph{} Learning a Riemannian metric is also related to finding a lower dimensional representation of a dataset. Work in this area includes linear methods such as principal component analysis and nonlinear methods such as spherical subfamily models \cite{Gous1998} or locally linear embedding \cite{Roweis2000} and curved multinomial subfamilies \cite{Hall2000}. Once such a submanifold is found, distances $d(x,y)$ may be computed as the lengths of shortest paths on the submanifold connecting $x$ and $y$. As shown in Section~\ref{sec:problem}, this approach is a limiting case of learning a Riemannian metric for the embedding high-dimensional space. \paragraph{} Lower dimensional representations are useful for visualizing high dimensional data. However, these methods assume strict conditions that are often violated in real-world, high dimensional data. The obtained submanifold is tuned to the training data and new data points will likely lie outside the submanifold due to noise. It is necessary to specify some way of projecting the off-manifold points into the manifold. There is no notion of non-Euclidean geometry outside the submanifold and if the estimated submanifold does not fit current and future data perfectly, Euclidean projections are usually used. \paragraph{} Another source of difficulty is estimating the dimension of the submanifold. The dimension of the submanifold is notoriously hard to estimate in high dimensional sparse datasets. Moreover, the data may have different lower dimensions in different locations or may lie on several disconnected submanifolds thus violating the assumptions underlying the submanifold approach. \paragraph{} We propose an alternative approach to the metric learning problem. The obtained metric is local, thus capturing local variations within the space, and is defined on the entire embedding space. A set of metric candidates is represented as a parametric family of transformations, or equivalently as a parametric family of statistical models and the obtained metric is chosen from it based on some performance criterion. \paragraph{} In Section~\ref{sec:problem} we discuss our formulation of the Riemannian metric problem. Section~\ref{sec:transformations} describes the set of metric candidates as pull-back metrics of a group of transformations followed by a discussion of the resulting generative model in Section~\ref{sec:simplexModel}. In Section~\ref{sec:experimentsMetric} we apply the framework to text classification and report experimental results on the WebKB data. %------------------------------------------------------------------------ \subsection{The Metric Learning Problem} \label{sec:problem} The metric learning problem may be formulated as follows. Given a differentiable manifold $\mathcal{M}$ and a dataset $D=\{x_1,\ldots, x_N\}\subset\mathcal{M}$, choose a Riemannian metric $g$ from a set of metric candidates $\mathcal{G}$. As in statistical inference, $\mathcal{G}$ may be a parametric family \begin{equation} \label{eq:paramMetricFam} \mathcal{G}=\{g^{\theta}:\theta\in\Theta\subset\mathbb{R}^k\} \end{equation} or as in nonparametric statistics a less constrained set of candidates. We focus on the parametric approach, as we believe it to generally perform better in high dimensional sparse data such as text documents. The reason we use a superscript $g^{\theta}$ is that the subscript of the metric is reserved for its value at a particular point of the manifold. \paragraph{} We propose to choose the metric based on maximizing the following objective function $\mathcal{O}(g,D)$ \begin{equation} \label{eq:normInvVol} \mathcal{O}(g,D)= \prod_{i=1}^{N}\frac{(\dvol g(x_i))^{-1}}{\int_{\mathcal{M}}(\dvol g(x))^{-1}\text{d}x} \end{equation} where $\dvol g(x)=\sqrt{\det G(x)}$ is the differential volume element, and $G(x)$ is the Gram matrix of the metric $g$ at the point $x$. Note that $\det G(x)>0$ since $G(x)$ is positive definite. \paragraph{} The volume element $\dvol g(x)$ summarizes the size of the metric $g$ at $x$ in one scalar. Intuitively, paths crossing areas with high volume will tend to be longer than the same paths over an area with low volume. Hence %For this reason, maximizing the inverse volume in \eqref{eq:normInvVol} will result in shorter curves across densely populated regions of $\mathcal{M}$. As a result, the geodesics will tend to pass through densely populated regions. This agrees with the intuition that distances between data points should be measured on the lower dimensional data submanifold, thus capturing the intrinsic geometrical structure of the data. \paragraph{} The normalization in \eqref{eq:normInvVol} is necessary since the problem is clearly unidentifiable without it. Metrics $cg$ with $0