\documentclass[twoside,11pt]{article}

%\usepackage{llncs}
%\usepackage{amsmath}
%\usepackage{amssymb}
\usepackage{standard}
\usepackage{graphicx}

\newcommand{\risi}[1]{\ignore{#1}}
\newcommand{\risiS}[2]{#1\ignore{#2}}

\renewcommand{\emph}[1]{\textbf{#1}}

\DeclareMathAlphabet{\mathpzc}{OT1}{pzc}{m}{it}
\newcommand{\pzc}[1]{\mathpzc{#1}}
%\newcommand{\pzc}[1]{\math{#1}}
\newcommand{\pzcx}{\pzc x}
%\newcommand{\rf}[1]{(\ref{#1})}
\newcommand{\todo}[1]{\mbox{}\\ \textbf{#1}\\  \mbox{}}

%\newtheorem

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

\usepackage{jmlr2e}

% Definitions of handy macros can go here

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

\newcommand{\Tcal}{\mathcal T}

% Heading arguments are {volume}{year}{pages}{submitted}{published}{authors}

\jmlrheading{5}{2004}{819--844}{1/04}{7/04}{Tony Jebara, Risi Kondor and Andrew Howard}
\ShortHeadings{Probability Product Kernels}{Jebara, Kondor and Howard}
\firstpageno{819}

\long\def\COMMENT#1{}
\newcommand{\step}{\mbox{\rm step}}
 \newcommand{\sign}{\mbox{\rm sign}}
\newcommand{\arginf}{\mbox{\rm arg}\inf}
\newcommand{\argsup}{\mbox{\rm arg}\sup}
\newcommand{\argmax}{\mbox{\rm arg}\max}
\newcommand{\argmin}{\mbox{\rm arg}\min}
\renewcommand{\[}{\begin{eqnarray*}}
\renewcommand{\]}{\end{eqnarray*}}


\begin{document}

\title{Probability Product Kernels}

\author{\name Tony Jebara \email jebara@cs.columbia.edu \\
        \name Risi Kondor \email risi@cs.columbia.edu \\
        \name Andrew Howard \email ahoward@cs.columbia.edu \\
       \addr Computer Science Department\\
       Columbia University\\
       New York, NY 10027, USA}
 
\editor{Kristin Bennett and Nicol{\`o} Cesa-Bianchi}

\maketitle


\begin{abstract}
The advantages of discriminative learning algorithms and kernel
machines are combined with generative modeling using a novel kernel
between distributions. In the probability product kernel, data points
in the input space are mapped to distributions over the sample space
and a general inner product is then evaluated as the integral of the
product of pairs of distributions. The kernel is straightforward to
evaluate for all exponential family models such as multinomials and
Gaussians and yields interesting nonlinear kernels. Furthermore, the
kernel is computable in closed form for latent distributions such as
mixture models, hidden Markov models and linear dynamical systems. For
intractable models, such as switching linear dynamical systems,
structured mean-field approximations can be brought to bear on the
kernel evaluation. For general distributions, even if an analytic
expression for the kernel is not feasible, we show a straightforward
sampling method to evaluate it. Thus, the kernel permits
discriminative learning methods, including support vector machines, to
exploit the properties, metrics and invariances of the generative
models we infer from each datum. Experiments are shown using
multinomial models for text, hidden Markov models for biological
data sets and linear dynamical systems for time series data.
\end{abstract}

\begin{keywords}
Kernels, support vector machines, generative models, Hellinger
divergence, Kullback-Leibler divergence, Bhattacharyya affinity, expected
likelihood, exponential family, graphical models, latent models,
hidden Markov models, dynamical systems, mean field
\end{keywords}

\section{Introduction}

Recent developments in machine learning, including the emergence of
support vector machines, have rekindled interest in kernel methods
\citep{Vap98,HastieBook}. Typically, kernels are designed and used by
the practitioner to introduce useful nonlinearities, embed prior
knowledge and handle unusual input spaces in a learning algorithm
\citep{SmolaBook}. In this article, we consider kernels between
distributions. Typically, kernels compute a generalized inner product
between two input objects $\pzc x$ and $\pzc x'$ which is equivalent
to apply a mapping function $\Phi$ to each object and then computing a
dot product between $\Phi(\pzc x)$ and $\Phi(\pzc x')$ in a Hilbert
space. We will instead consider the case where the mapping $\Phi(\pzc
x)$ is a probability distribution $p(x|\pzc x)$, thus restricting the
Hilbert space since the space of distributions is trivially embedded
in the Hilbert space of functions. However, this restriction to
positive normalized mappings is in fact not too limiting since general
$\Phi({\pzc x})$ mappings and the nonlinear decision boundaries they
generate can often be mimicked by a probabilistic mapping. However, a
probabilistic mapping facilitates and elucidates kernel design in
general. It permits kernel design to leverage the large body of tools
in statistical and generative modeling including Bayesian networks
\citep{Pearl}. For instance, maximum likelihood or Bayesian estimates
may be used to set certain aspects of the mapping to Hilbert space.

In this paper, we consider the mapping from datum $\pzc x$ to a
probability distribution $p(x|\pzc x)$ which has a straightforward
inner product kernel $k(\pzc x,\pzc x')=\int p^\rho(x|\pzc x)
p^\rho(x|\pzc x') dx$ (for a scalar $\rho$ which we typically set to
$1$ or $1/2$). This kernel is merely an inner product between two
distributions and, if these distributions are known directly, it
corresponds to a linear classifier in the space of probability
distributions. If the distributions themselves are not available and
only observed data is given, we propose using either Bayesian or
Frequentist estimators to map a single input datum (or multiple
inputs) into probability distributions and subsequently compute the
kernel. In other words, we map points as ${\pzc x} \rightarrow
p(x|\pzc x)$ and ${\pzc x'} \rightarrow p(x|\pzc x')$. For the setting
of $\rho=1/2$ the kernel is none other than the classical
Bhattacharyya similarity measure \citep{KonJeb03}, the affinity
corresponding to Hellinger divergence \citep{JebKon03}.\footnote{This has
interesting connections to Kullback-Leibler divergence
\citep{Topsoe}.}

One goal of this article is to explore a contact point between
discriminative learning (support vector machines and kernels) and
generative learning (distributions and graphical
models). Discriminative learning directly optimizes performance for a
given classification or regression task.  Meanwhile, generative
learning provides a rich palette of tools for exploring models,
accommodating unusual input spaces and handling priors. One approach
to marrying the two is to estimate parameters in generative models
with discriminative learning algorithms and optimize performance on a
particular task. Examples of such approaches include conditional
learning \citep{bengio96inputoutput} or large margin generative
modeling \citep{JaaMeiJeb99}. Another approach is to use kernels to
integrate the generative models within a discriminative learning
paradigm. The proposed probability product kernel falls into this
latter category. Previous efforts to build kernels that
accommodate probability distributions include the Fisher kernel
\citep{JaaHau98}, the heat kernel \citep{LafLeb02} and kernels arising
from exponentiating Kullback-Leibler divergences
\citep{MorHoVas04}. We discuss, compare and contrast these approaches
to the probability product kernel in Section~\ref{sec: relation}. One
compelling feature of the new kernel is that it is straightforward and
efficient to compute over a wide range of distributions and generative
models while still producing interesting nonlinear behavior in, for
example, support vector machine (SVM) classifiers. The generative
models we consider include Gaussians, multinomials, the exponential
family, mixture models, hidden Markov models, a wide range of
graphical models and (via sampling and mean-field methods) intractable
graphical models. The flexibility in choosing a distribution from the
wide range of generative models in the field permits the kernel to
readily accommodate many interesting input spaces including sequences,
counts, graphs, fields and so forth. It also inherits the properties,
priors and invariances that may be straightforward to design within a
generative modeling paradigm and propagates them into the kernel
learning machine. This article thus gives a generative modeling route
to kernel design to complement the family of other kernel engineering
methods including super-kernels \citep{superkernels}, convolutional
kernels \citep{Hau99,ColDuf02}, alignment kernels \citep{Wat00},
rational kernels \citep{CorHafMoh03} and string kernels
\citep{Leslie02,VisSmo03}.

This paper is organized as follows. In Section~\ref{sec: kernel}, we
introduce the probability product kernel and point out connections to
other probabilistic kernels such as Fisher kernels and exponentiated
Kullback-Leibler divergences. Estimation methods for mapping points
probabilistically to Hilbert space are outlined. In Section~\ref{sec:
exp fam} we elaborate the kernel for the exponential family and show
its specific form for classical distributions such as the Gaussian and
multinomial. Interestingly, the kernel can be computed for
mixture models and graphical models in general and this is elaborated
in Section~\ref{sec: graphical}. Intractable graphical models are then
handled via structured mean-field approximations in Section~\ref{sec:
mean field}. Alternatively, we describe sampling methods for
approximately computing the kernel. Section~\ref{sec: relation}
discusses the differences and similarities between our probability
product kernels and previously presented probabilistic kernels. We
then show experiments with text data sets and multinomial kernels, with
biological sequence data sets using hidden Markov model kernels, and
with continuous time series data sets using linear dynamical system
kernels. The article then concludes with a brief discussion.

\section{A Kernel Between Distributions}\label{sec: kernel}

Given a positive (semi-)definite kernel $k:\Xcal\times\Xcal\mapsto\RR$ on the input 
space $\Xcal$ and examples $\seq{\pzcx}{m}\<\in\Xcal$ with corresponding labels $\seq{y}{m}$, 
kernel based learning algorithms return a hypothesis of the form 
$h(\pzcx)=\Sum{i}{m}\alpha_i\,k(\pzcx_i,\pzcx)+b$. 
%Equivalently, regarding $h$ as a vector in the Reproducing Kernel Hilbert Space 
%$\Hcal$ associated with $k$, we can say that $h$ is in the span of 
%$\{k_{\pzcx_i}=k(\pzcx_i,\cdot)\}_{i=1}^m$. 
The role of the kernel is to capture our prior assumptions of similarity in \m{\Xcal}. 
When \m{k(\pzcx,\pzcx')} is large, we expect that the corresponding labels be similar or 
the same.%, whereas when \m{k(\pzcx,\pzcx')} is low, this is less likely. 

Classically, the inputs \m{\cbr{\pzcx_i}_{i=1}^m} are often represented 
as vectors in \m{\RR^n}, and \m{k} is chosen to be one of a small number of popular 
positive definite functions on \m{\RR^n}, such as  the Gaussian RBF
\begin{equation}\label{eq: Gaussian RBF}
  k(\pzcx,\pzcx')=e^{-\nm{\pzcx-\pzcx'}^2/(2\sigma^2)}.
\end{equation}
%Recently, a number of imaginative alternatives have been proposed for cases when the 
%inputs are not easily embeddable in \m{\RR^n}, specifically, when they are discrete objects  
%such as a graphs, tree, vertices of a graph, automata, etc.. 
%However, once an appropriate kernel has been found, the nature of the input objects 
%and their representation is of no relevance to the learning algorithm itself. 
Once we have found an appropriate kernel,  
the problem becomes accessible to 
a whole array of non-parametric discriminative kernel based learning 
algorithms, such as support vector machines, Gaussian processes, etc. \citep{SchSmo02}. 

In contrast, generative models fit probability distributions \m{p(\pzcx)} to 
\m{\seq{\pzcx}{m}} and base their predictions on the likelihood under different models. 
When faced with a    
discriminative task,  approaching it from the generative angle is 
generally suboptimal.  
On the other hand, it is often easier to capture the structure of complex objects 
with generative models than directly with kernels. 
Specific examples include 
multinomial models for text documents, 
hidden Markov models for speech, Markov random fields for images and 
linear dynamical systems for motion.   
%graphical models for .... and .... . 
In the current paper we aim to combine the advantages of both worlds by 
\begin{enumerate}
\item fitting separate probabilistic models \m{p_1(x),p_2(x),\ldots p_m(x)} 
  to \m{\seq{\pzcx}{m}};
\item defining a novel kernel \m{k^{\text{prob}}(p,p')} between probability distributions on \m{\Xcal}; 
\item finally, defining the kernel between examples to equal \m{k^{\text{prob}}} 
  between the corresponding distributions: 
\begin{equation*}
  k(\pzcx,\pzcx')=k^\text{prob}(p,p').
\end{equation*}
\end{enumerate}
We then plug this kernel into any established kernel based learning algorithm 
and proceed as usual. 
We define \m{k^\text{prob}} in a rather general form and then investigate special cases. 

\begin{definition}
  Let \m{p} and \m{p'} be probability distributions on a space \m{\Xcal} 
  %\m{q} be a bounded positive function on \m{\Xcal} 
   and \m{\rho} be a positive constant.  
  Assume that \m{p^\rho,p'^\rho\tin L_{2}(\Xcal)}, i.e. that 
  %\begin{equation*}
    \m{\int_\Xcal p(x)^{2\rho} dx} and \m{\int_\Xcal p'(x)^{2\rho} dx} 
  %\end{equation*}
  %and the analogous integral for \m{p'} 
  are well defined (not infinity). 
  The {\rm \emph{probability product kernel}} between distributions \m{p} and \m{p'} 
  % with parameters \m{q} and \m{\rho} 
  is defined 
  \begin{equation}\label{eq: prob prod kernel}
    k^\textrm{prob}(p,p')=\int_\Xcal p(x)^\rho p'(x)^\rho\,dx=\inp{p^\rho,{p'}^\rho}_{L_2}. 
  \end{equation}
\end{definition}

It is well known that \m{L_2(\Xcal)} is a Hilbert space, 
hence for any set \m{\mc{P}} of probability 
distributions over \m{\Xcal} such that  
\m{\int_\Xcal p(x)^{2\rho} dx} is finite for any \m{p\tin\mc{P}}, 
the kernel defined by \rf{eq: prob prod kernel} is positive definite. 
The probability product kernel is also a simple and 
intuitively compelling notion of similarity between distributions. 

\ignore{
%Unless noted otherwise, in most of what follows, we will simply take
%\m{q(x)\<=1} but other choices include taking $q(x)$ to be a power of
%the maximum likelihood fit of a model to the whole aggregated input
% data set. Furthermore, 
In the following we will drop the qualifying superscript from
\m{k^\text{prob}}, and by slight abuse of notation use the same symbol
\m{k} to denote both the kernel between distributions and the kernel
between examples.

It is also easy to see that \m{k} is positive definite. 
Recall that, by definition, \m{k} is positive definite if for any 
\m{\seq{x}{m}\<\tin\Xcal} and coefficients \m{\seq{\alpha}{m}} not all zero   
\m{\Sum{i}{m}\Sum{j}{m}\alpha_i \alpha_j k(x_i,x_j)>0}. 
%is greater than (or greater than or equal to) zero. 

\begin{lemma}
  Let \m{\psi} be a mapping from some space \m{\Omega} to a linear inner product 
  space \m{\mathcal Z} and \m{k:\Omega\times\Omega\mapsto\RR} be such that 
  \m{g(x,x')=\inp{\psi(x),\psi(x')}_{\mathcal Z}}. 
  Then \m{g} is positive definite on \m{\Omega}. 
\end{lemma}
\begin{proof}
  For any positive integer \m{m}, any \m{\seq{x}{m}\<\in\Xcal} and \m{\seq{\alpha}{m}} 
  not all zero, the vector \m{v=\Sum{i}{m}\alpha_i \psi(x)} is a non-zero element of 
  \m{\mathcal Z}, and by linearity 
  \begin{equation*}
    \Sum{i}{m}\Sum{j}{m}\alpha_i \alpha_j g(x_i,x_j)=\inp{v,v}. 
  \end{equation*}
  By the definition of the inner product, this quantity is strictly positive.  
\end{proof}

\begin{proposition}
  The probability product kernel defined in \rf{eq: prob prod kernel} is positive 
  definite. 
\end{proposition}
\begin{proof}
  Consider the mapping \m{\psi:L_{2\rho}(\Xcal)\mapsto L_{2\rho}(\Xcal)}
  given by \m{\psi:p\mapsto p \sqrt{q}}.  
  The above lemma proves that 
  \begin{equation}
    g(p,p')=\inp{\psi(p),\psi(p')}_{L_{2\rho}}=
    \sqb{\int_\Xcal p(x)^\rho p'(x)^\rho q(x)\,dx}^{1/\br{2\rho}}  
  \end{equation}
  is positive definite on the space of probability distributions. 
  It is well known that any positive power of a positive definite kernel is also 
  positive definite, hence \m{k^\text{prob}(p,p')=\br{g(p,p')}^{2\rho}} is positive definite.  
  %This maps all our distributions to \m{L_{2\rho}} in such a way that the kernel 
  \ignore{
  Since \m{p(x)} is a probability distribution,  
  \begin{equation*}
    \int_\Xcal p(x)^{2\rho}\,q(x)\,dx \leq 
    \br{\sup_{x\in \Xcal} q(x)}\underbrace{\int_\Xcal p(x)\,dx}_{1}, 
  \end{equation*}
  hence \m{\psi:p\mapsto p\sqrt q} defines a mapping 
  from probability distributions on \m{\Xcal} 
  to \m{L_2(\Xcal)}. In this space the kernel becomes the inner product, 
  \m{k(p,p')=\inp{\psi(p),\psi(p')}_{L_2}},  and we can apply our previous lemma. 
  }
\end{proof}
}
%The probability product kernel is also related to some well known measures of similarity 
%in statstics. 
\subsection{Special Cases and Relationship to Statistical Affinities}

For \m{\rho=1/2} 
\begin{equation*}
  k(p,p')=\int \sqrt{p(x)} \sqrt{p'(x)}\,dx,
\end{equation*}
which we shall call the \emph{Bhattacharyya kernel}, because 
in the statistics literature it is known as Bhattacharyya's affinity 
between distributions \citep{bhattacharyya}, 
related to the better-known Hellinger's distance 
\begin{equation*}
  H(p,p')=\ovr{2} \int  \br{\sqrt{p(x)}-\sqrt{p'(x)}}^2 dx
\end{equation*}
by \m{H(p,p')=\sqrt{2-2k(p,p')}}\:. 
Hellinger's distance can be seen as a symmetric approximation to the Kullback-Leibler (KL) 
divergence, and in fact is a bound on KL, as shown in \citep{Topsoe}, where relationships 
between several common divergences are discussed.  
The Bhattacharyya kernel was first introduced in \citep{KonJeb03} 
and has the important special property $k(\pzc x,\pzc x)=1$.

When $\rho=1$, the kernel takes the form of the expectation of one distribution 
under the other: 
\begin{equation}
k(\pzc x, \pzc x')=\int p(x)\, p'(x)\:dx =\E_{p}[p'(x)]=\E_{p'}[p(x)]. 
\end{equation}
We call this the \emph{expected likelihood kernel}.

It is worth noting that when dealing with distributions over discrete spaces 
\m{\Xcal=\cbr{\seqi{x}}},  
probability product kernels have a simple geometrical interpretation. 
In this case \m{p} can be represented as a vector %of probabilities 
\m{p=\br{\seqi{p}}} where \m{p_i=\Pr\br{X\es\, x_i}}. 
The probability product kernel is then the dot product between 
\m{\tilde p=\br{\seqi{p^\rho}}} and \m{\tilde p'=\br{\seqi{{p'}^\rho}}}. 

In particular, in the case of  the Bhattacharyya kernel,  
\m{\tilde p} and \m{\tilde p'} are vectors on the unit 
sphere.
This type of similarity measure is not unknown in the literature. 
In text recognition, for example, the so-called ``cosine-similarity'' 
(i.e. the dot product measure) is well entrenched, and it has been observed that 
preprocessing word counts by taking square roots often improves performance 
\citep{GolSah98,CutKarPedTuk92}. 
Lafferty and Lebanon  also arrive at a similar kernel, but from a very different angle, 
looking at the diffusion kernel on the statistical manifold of multinomial distributions 
\citep{LafLeb02}. 

\subsection{Frequentist and Bayesian Methods of Estimation}

Various statistical estimation methods can be used to fit the distributions 
\m{\seq{p}{m}} to the examples \m{\seq{\pzcx}{m}}. 
Perhaps it is most immediate to consider a parametric family \m{\cbr{p_\theta(x)}}, 
take the maximum likelihood estimators 
\m{\hat \theta_i =\arg \max_\theta p_\theta(\pzcx_i)}, 
and set \m{p_i(x)=p_{\hat\theta_i}(x)}. 

The alternative, Bayesian, strategy is to postulate a prior on \m{\theta}, 
invoke Bayes' rule 
\begin{equation*}
p(\theta|\pzc{x})=\fr{p(\pzc{x}|\theta)\,p(\theta)}
{\int p(\pzc{x}|\theta)\,p(\theta) \:d\theta}\:,
\end{equation*}
and then use either the distribution $p(x|\hat \theta_{\text{MAP}})$ based on the 
maximum a posteriori estimator 
$\hat \theta_\text{MAP}=\arg \max_\theta p(\theta|\pzc x)$, or the true posterior 
\begin{equation}
p(x|\pzc{x})=\int p(x|\theta)\,p(\theta|\pzc{x})\:d\theta.
\end{equation}

At this point the reader might be wondering to what extent it is justifiable to 
fit distributions to single data points. 
It is important to bear in mind that \m{p_i(x)} is but an intermediate step in 
forming the kernel, and is not necessarily meant to model anything in and of itself. 
The motivation is to exploit the way that probability distributions capture similarity: 
assuming that our model is appropriate, if a distribution fit to one data point gives high likelihood to 
another data point, this indicates that the two data points are in some sense similar. 
This is particularly true of intricate graphical models commonly used to model structured 
data, such as times series, images, etc.. 
Such models can capture relatively complicated relationships 
in real world data. 
Probability product kernels allow kernel methods to harness the power of such models. 

%At the opposite end of the spectrum, 
When \m{\pzcx} is just a point in \m{\RR^n} with no further structure,
probability product kernels are less compelling. \COMMENT{more difficult to}
Nevertheless, it is worth noting that even the Gaussian RBF
kernel~\rf{eq: Gaussian RBF} can be regarded as a probability product
kernel (by setting setting \m{\rho\<=1} and choosing \m{p_i(x)} to be a \m{\sigma/2} variance
Gaussian fit to \m{\pzcx_i} by maximum likelihood).  In particular
situations when the point \m{\pzcx} does not yield a reliable estimate
of $p$ (typically because the class of generative models is too
flexible relative to the datum), we may consider regularizing the
maximum likelihood estimate of each probability by fitting it to the
neighbourhood of a point or by employing other more global estimators
of the densities. Other related methods that use points to estimate
local distributions including kernel density estimation (KDE) where
the global density is composed of local density models centered on
each datum \citep{Sil86}. We next discuss a variety of distributions
(in particular, the exponential family) where maximum likelihood
estimation is well behaved and the probability product kernel is
straightforward to compute.

\section{Exponential Families}\label{sec: exp fam}


\begin{table}[t]
\begin{center}
\begin{tabular}{ | l | c | c | c |}
\hline
 & \m{\Tcal(x)} & ${\cal A}(x)$ & ${\cal K}(\theta)$ \\ %&  constraint\\
\hline
\hline
Gaussian (\,\m{\theta\<=\mu,\ \sigma^2\<=1}\,) & \m{x} &
$ -\frac{1}{2} x^T\! x - \frac{D}{2} \log(2\pi) $ &
$ \qqquad \frac{1}{2} \theta^T\! \theta  \qqquad $  \\
\hline
Gaussian (\,\m{\theta\<=\sigma^2\!,\ \mu\<=0}\,) & \m{x} &
$ - \frac{1}{2} \log(2\pi) $ &
$ -\frac{1}{2} \log \theta$  \\ % & $\theta > 0$ \\
\hline
Exponential (\,\m{\theta\<=-1/\beta}\,)& \m{x} & $ 0 $ & $ -\log(-\theta)$  \\  % & $\theta < 0$ \\
\hline
Gamma (\m{\theta\<=\alpha})& \m{\log x} & $-\log x-x$  & $\log \Gamma(\theta)$  \\ % & $\theta > 0$ \\
\hline
Poisson ($\,\theta\<=\log p\,$) & \m{x} & $-\log(x!)$ & $\exp\theta$   \\% & \\
\hline
Multinomial (\,\m{\theta_i\<=\log \alpha_i }\,)& \m{\br{x_1,x_2,\ldots}} &
$ \log\brbig{\brbig{\sum_i x_i}!/\brbig{\prod_i x_i!}}$ & $ 1 $  \\ %&
%$  \log\brbig{\brbig{\Sum{i}{D}x_i}!} \<- \Sum{i}{m}\log(x_i!) $  & $1$ & \\
%$  \eta \log( 1\<+\sum_{d=1}^D \exp(\theta_d)) $  & \\
\hline
\end{tabular}
\label{tbl:efamilies}
\end{center}
\caption{\label{table:exp fam} Some well-known exponential family distributions}
\end{table}


A family of distributions parameterized by \m{\theta\<\in\RR^D}
is said to form an exponential family if its members are of the form  
\[
p_\theta(x) &=& \exp({\cal A}(x)+\theta^T \mc{T}(x)-{\cal K}(\theta)).
\]
Here ${\cal A}$ is called the measure, \m{\cal K} is the 
the cumulant generating function and \m{\mc T} are the sufficient statistics.  
%(all of these can be scalar or vector valued).  
%Often $\mc{T}$ is of the simple form $\mc{T}\!=\theta^T\!x$. 

Some of the most common statistical distributions, including 
the Gaussian, multinomial, Poisson and 
Gamma distributions, are exponential families (Table \ref{table:exp fam}).  
Note that to ensure that \m{p_\theta(x)} is normalized, 
${\cal A}$,\ ${\cal K}$ and $\theta$ must be related through the Laplace transform 
\[{\cal K}(\theta)=\log \int\exp\br{{\cal A}(x)+\theta^T \mc{T}(x)}\:dx.\]
%to ensure that $p(x|\theta)$ is normalized. 

The Bhattacharyya kernel ($\rho\<=1/2$) 
can be computed in closed form for any exponential family:
\[
k(\pzc x,\pzc x')=k(p,p') &=& \int_x p_\theta(x)^{1/2} p_{\theta'}(x)^{1/2}\, dx \\
&=& \int_x \exp \left ({\cal A}(x)+
\left(  \hbox{$\textstyle \frac{1}{2}$}\theta+ 
\hbox{$\textstyle \frac{1}{2}$}\theta' \right )^T\! \Tcal(x)- 
\hbox{$\textstyle \frac{1}{2}$}{\cal K}(\theta)- 
\hbox{$\textstyle \frac{1}{2}$}{\cal K}(\theta') \right )dx \\
&=& \exp\left ({\cal K}\left(  \hbox{$\textstyle \frac{1}{2}$}\theta+ 
\hbox{$\textstyle \frac{1}{2}$}\theta' \right )- 
\hbox{$\textstyle \frac{1}{2}$}{\cal K}(\theta)- 
\hbox{$\textstyle \frac{1}{2}$}{\cal K}(\theta') \right ).
\]
For \m{\rho\<\neq 1/2},\ \m{k(\pzc x,\pzc x')} 
can only be written in closed form when  
\m{\mc A} and \m{\mc T} obey some special properties. 
Specifically, if \m{2\rho\mc{A}(x)=\mc{A}(\eta x)} %for some \m{\eta\tin\RR} 
for some \m{\eta} 
and \m{\Tcal} is linear in \m{x}, then 
\begin{multline*}
  \log \int\exp\br{2\rho\ts\mc{ A}(x)+\rho\br{\theta+\theta'}^T\! \mc{T}(x)}\:dx=\\
  \log \sqbP{\ovr{\eta} \int\exp\brBig{{\cal A}(\eta x)+\fr{\rho}{\eta}\br{\theta+\theta'}^T \mc{T}(\eta x)}\:d(\eta x)}=
  \mc{K}\brbig{{\textstyle \fr{\rho}{\eta}}\br{\theta+\theta'}} -\log\eta,  
\end{multline*}
hence 
\begin{equation*}
  k(p,p')=\exp\sqb{\mc{K}\brbig{{\textstyle \fr{\rho}{\eta}}\br{\theta+\theta'}} -\log\eta-
  \hbox{$\textstyle \frac{\rho}{2}$}{\cal K}(\theta)- 
  \hbox{$\textstyle \frac{\rho}{2}$}{\cal K}(\theta')}. 
\end{equation*}
%For some of the most common families we can also get closed form expressions 
%for more general probability product kernels. 
In the following we  look at some specific examples. 

\subsection{The Gaussian Distribution}

The $D$ dimensional Gaussian distribution $p(x)\!=\!\mc{N}(\mu,\Sigma)$ is of the form
\[p(x)={(2\pi)^{-D/2} \abs{\Sigma}^{-1/2}}\:
\exp\br{-{\textstyle \ovr{2}}(x\ms\mu)^T \Sigma^{-1} (x\ms\mu)}\]
%written $p(x)\!\sim\!\mc{N}(\mu,\Sigma)$, 
where $\Sigma$ is a positive definite matrix and 
$\abs{\Sigma}$ denotes its determinant. 
For a pair of Gaussians $p\!=\! \mc{N}(\mu,\Sigma)$ 
and $p'\!=\! \mc{N}(\mu',\Sigma')$, completing the square in the exponent gives the general probability product kernel 
\begin{multline}\label{eq: kernel between Gaussians}
K_\rho(\pzc x,\pzc x') = K_\rho(p,p')= %\br{\inp{p,p'}_{L_\rho}}^\rho=
\int_{\Real^D} p(x)^{\rho} p'(x)^{\rho} dx = \\
\br{2\pi}^{(1-2\rho)D/2} \rho^{-D/2}\:
{\bigl\vert\Sigma^\dagger\ts \bigr\vert ^{1/2}}
{\bigl\vert\ts\Sigma\ts\bigr\vert^{-\rho/2} 
{\bigl\vert\ts\Sigma'\ts\bigr\vert}^{-\rho/2}}
%& & \:\:\:\: 
\exp\br{ 
-\fr{\rho}{2}\br{\mu^T \Sigma^{-1} \mu
+{\mu'}^T {\Sigma'}^{-1} \mu'
-{\mu^\dag}^T{\Sigma^\dag} \mu^\dag}}, 
\end{multline}
%\todo{Check this!}
where $\Sigma^\dag=\big(\Sigma^{-1}\!\ps {\Sigma'}^{-1}\big)^{-1}$ 
and $\mu^\dag=\Sigma^{-1}\mu+{\Sigma'}^{-1}\mu'$.
When the covariance is isotropic and fixed,  
%In the special case that the of fixed and isotropic covariance 
$\Sigma=\sigma^2 I$, this simplifies to 
\[K_\rho(p,p')=(2\rho)^{-D/2}\, (2\pi \sigma^2)^{(1-2\rho)D/2}\:
e^{-\nm{\mu-\mu'}^2/(4\sigma^2\!/\rho)},
\]
which, for $\rho\<= 1$ (the expected likelihood kernel) simply gives 
%the following Gaussian (whose variance is effectively
%double the original $\Sigma=\sigma^2 I$):
\[k(p,p')={\textstyle \ovr{\br{4\pi {\sigma}^2}^{D/2}}}\:
e^{-\nm{\mu'\ms\mu}^2/(4{\sigma}^2)},\]
recovering, up to a constant factor, the Gaussian RBF kernel~\rf{eq: Gaussian RBF}. 

\subsection{The Bernoulli Distribution}

The Bernoulli distribution $p(x) = \gamma^x (1-\gamma)^{1-x}$
%models binary observations $x\tin\{0,1\}$. 
with parameter $\gamma\tin(0,1)$,
and its $D$ dimensional variant, sometimes referred to as Naive Bayes,%\vspace{-6pt}
\[
p(x) = \prod_{d=1}^D \gamma_d^{x_d} (1-\gamma_d)^{1-x_d}
\]
with $\gamma\tin(0,1)^D$, are used to model binary $x\tin\{0,1\}$ or
multidimensional binary $x\tin\{0,1\}^D$ observations. 
The probability product kernel %between a pair of such distributions
\[
K_\rho(\pzc x,\pzc x')=
K_\rho(p,p') = \sum_{x\in\{0,1\}^D}  \prod_{d=1}^D ({\gamma}_d {\gamma'_d} )^{\rho x_d}
((1-{\gamma}_d)(1-{\gamma'_d}))^{\rho (1-x_d)}
\]
factorizes as 
%trivially  (for any setting of $\rho$) as:
\[
K_\rho(p,p') = 
\prod_{d=1}^D \left [ ({\gamma}_d {\gamma'_d} )^\rho + 
(1-{\gamma}_d)^\rho (1-{\gamma'_d})^\rho \right ].\] %\: \: = \:\:  
\COMMENT{In the special case of the expected likelihood kernel ($\rho\es 1$), 
this further simplifies to \vspace{-6pt} 
\[k(p,p')=\prod_{d=1}^D  (1 -{\gamma}_d - {\gamma'_d} ). \]}

\subsection{The Multinomial Distribution}


The multinomial model 
\[
p(x)=\fr{s!}{x_1!\, x_2!\ldots x_D!}\: \alpha_1^{x_1} \alpha_2^{x_2}\ldots \alpha_D^{x_D}
\]
with parameter vector $\alpha=(\seq{\alpha}{D})$ subject to $\Sum{d}{D}\alpha_d=1$ 
is commonly used to model discrete integer counts $x=(\seq{x}{D})$ with 
\m{\Sum{i}{D} x_i=s} fixed, 
such as the number of occurrences of words in a document. 
% and multinomial coefficient defined as 
%\[{X \choose \sprod{x}{D}}=\fr{X!}{x_1! x_2!\ldots x_D!}.\]
The maximum likelihood estimate given observations $\useq{x}{n}$ is
\[\hat\alpha_d=\fr{\Sum{i}{n}x_d^{(i)}}{\Sum{i}{n}\Sum{d}{D}x_d^{(i)}}\:.\]

For fixed \m{s}, the Bhattacharyya kernel (\m{\rho\<=1/2}) %$k(\pzc x,\pzc x')=k(p,p')$
can be computed explicitly using the multinomial theorem:
\begin{eqnarray}
\label{eqn:bhatmulti}
k(p,p')=
\sum_{
\begin{array}{c}
\scriptstyle x=(\seq{x}{D})\\ 
{\scriptstyle \sum_i x_i=s}
\end{array}}
\fr{s!}{x_1!\, x_2!\ldots x_D!}\; \Prod{d}{D} (\alpha_d \alpha'_d)^{x_d/2}
\:=\sqb{\Sum{d}{D} (\alpha_d \alpha_d')^{1/2}}^s
\end{eqnarray}
which is equivalent to the homogeneous polynomial kernel of order $s$ between 
the vectors $(\sqrt{\alpha_1},\sqrt{\alpha_2},\ldots,\sqrt{\alpha_D})$ and 
$(\sqrt{\alpha'_1},\sqrt{\alpha'_2},\ldots,\sqrt{\alpha'_D})$. 
When $s$ is not constant, we can sum over all its possible values
\[
k(p,p')=\sum_{s=0}^{\infty} \sqb{\:\Sum{d}{D} (\alpha_d \alpha'_d)^{1/2}}^s =
\br{1-\Sum{d}{D} (\alpha_d \alpha'_d)^{1/2}}^{-1}
\] 
or weight each
power differently, leading to a power series expansion.  
\ignore{
For the general
case of $\rho \!\not =\! 1/2$, a general formula for discrete
multinomial events is available (if $\rho \!\not =\! 1/2$ there is no
general closed form formula for counts except for $X\es 1$ where, as
opposed to counts, we really have single mutually exclusive events). This discrete events scenario if arguably more relevant and yields the form:
%\vspace{-6pt}
\[k(p,p')=\Sum{d}{D} (\alpha_d \alpha_d')^\rho.\]
}

\subsection{The Gamma and Exponential Distributions}

The gamma distribution \m{\Gamma(\alpha, \beta)} 
parameterized by \m{\alpha\<>0} and \m{\beta\<>0} is of the form 
\begin{equation*}
  p(x)=\ovr{\Gamma(\alpha)\beta^\alpha}\:x^{\alpha-1}\, e^{-x/\beta}. 
\end{equation*}
The probability product kernel between \m{p\sim\Gamma(\alpha,\beta)} and 
\m{p'\sim\Gamma(\alpha',\beta')} will then be 
\begin{equation*}
  k_\rho(p,p')=\fr{\Gamma(\alpha^\dag)\,{\beta^\dag}^{\alpha^\dag}}
  {\sqb{\Gamma(\alpha)\,\beta^\alpha\,\Gamma(\alpha')\,{\beta'}^{\alpha'}}^\rho}
\end{equation*}
where \m{\alpha^\dag=\rho\br{\alpha\<+\alpha'\<-2}+1} and 
\m{1/\beta^\dag=\rho\br{1/\beta+1/\beta'}}. 
When \m{\rho\<=1/2} this simplifies to \m{\alpha^\dag=\br{\alpha\<+\alpha'}/2}\: and 
\m{1/\beta^\dag=\br{1/\beta+1/\beta'}/2}. 

The exponential distribution 
%\begin{equation*}
  \m{p(x)=\ovr{\beta}\:e^{-x/\beta}}
%\end{equation*}
can be regarded as a special case of the gamma family with 
\m{\alpha\<=1}. In this case 
\begin{equation*}
  k_\rho(p,p')=\fr{\rho\,(1/\beta+1/\beta')}{\br{\beta\beta'}^\rho}\:. 
\end{equation*}

\section{Latent and Graphical Models} 
\label{sec: graphical}

We next upgrade beyond the exponential family of distributions and
derive the probability product kernel for more versatile generative
distributions. This is done by considering latent variables and
structured graphical models. While introducing additional complexity
in the generative model and hence providing a more elaborate
probability product kernel, these models do have the caveat that
estimators such as maximum likelihood are not as well-behaved as in
the simple exponential family models and may involve
expectation-maximization or other only locally optimal estimates. We
first discuss the simplest latent variable models, mixture models,
which correspond to a single unobserved discrete parent of the
emission and then upgrade to more general graphical models such as
hidden Markov models.

Latent variable models identify a split between the variables in $x$
such that some are observed and are denoted using $x_{\rm o}$ (where
the lower-case 'o' is short for 'observed') while others are hidden
and denoted using $x_{\rm h}$ (where the lower-case 'h' is short for
'hidden'). In the latent case, the probability product kernel is
computed only by the inner product between two distributions $p(x_{\rm
o})$ and $p'(x_{\rm o})$ over the observed components, in other words
$k(p,p')=\int p(x_{\rm o}) p'(x_{\rm o})dx_{\rm o}$. The additional
variables $x_{\rm h}$ are incomplete observations that are not part of
the original sample space (where the datum or data points to kernelize
exist) but an augmentation of it. The desired $p(x_{\rm o})$
therefore, involves marginalizing away the hidden variables:
\[
p(x_{\rm o}) &=& \sum_{x_{\rm h}} p(x_{\rm o},x_{\rm h}) \:\: = \:\: \sum_{x_{\rm h}} p(x_{\rm h}) p(x_{\rm o}|x_{\rm h}).
\]
For instance, we may consider a mixture of exponential families model
(a direct generalization of the exponential family model) where
$x_{\rm h}$ is a single discrete variable and where $p(x_{\rm
o}|x_{\rm h})$ are exponential family distributions themselves. The
probability product kernel is then given as usual between two such
distributions $p(x_{\rm o})$ and $p'(x_{\rm o})$ which expand as
follows:
\[
k(p,p') &=& \sum_{x_{\rm o}} \left (p(x) \right )^\rho \left ( p'(x) \right )^\rho \:\: = \:\: \sum_{x_{\rm o}} \left ( \sum_{x_{\rm h}} p(x_{\rm h}) p(x_{\rm o}|x_{\rm h}) \right )^\rho
 \left ( \sum_{x_{\rm h}'} p'(x_{\rm h}') p'(x_{\rm o}|x_{\rm h}') \right )^\rho.
\]
We next note a slight modification to the kernel which makes latent
variable models tractable when $\rho \neq 1$. This alternative kernel
is denoted ${\tilde k}(p,p')$ and also satisfies Mercer's condition
using the same line of reasoning as was employed for $k(p,p')$ in the
previous sections. Essentially, for ${\tilde k}(p,p')$ we will assume
that the power operation involving $\rho$ is performed on each entry
of the joint probability distribution $p(x_{\rm o},x_{\rm h})$ instead
of on the marginalized $p(x_{\rm o})$ alone as follows:
\[
{\tilde k}(p,p') &=& 
\sum_{x_{\rm o}} \sum_{x_{\rm h}} \left ( p(x_{\rm h}) p(x_{\rm o}|x_{\rm h}) \right )^\rho
  \sum_{x_{\rm h}'} \left ( p'(x_{\rm h}') p'(x_{\rm o}|x_{\rm h}') \right )^\rho .
\]
While the original kernel $k(p,p')$ involved the mapping to Hilbert
space given by ${\pzc x} \rightarrow p(x_{\rm o})$, the above ${\tilde
k}(p,p')$ corresponds to mapping each datum to an augmented
distribution over a more complex marginalized space where probability
values are squashed (or raised) by a power of $\rho$. Clearly,
$k(p,p')$ and ${\tilde k}(p,p')$ are formally equivalent when
$\rho=1$. However, for other settings of $\rho$, we prefer handling
latent variables with ${\tilde k}(p,p')$ (and at times omit the $\sim$
symbol when it is self-evident) since it readily accommodates
efficient marginalization and propagation algorithms (such as the
junction tree algorithm \citep{JordanNewBook}).

One open issue with latent variables is that the ${\tilde k}$ kernels
that marginalize over them may produce different values if the
underlying latent distributions have different joint distributions
over hidden and observed variables even though the marginal
distribution over observed variables stays the same. For instance, we
may have a two-component mixture model for $p$ with both components
having the same identical emission distributions yet the kernel
${\tilde k}(p,p')$ evaluates to a different value if we collapse the
identical components in $p$ into one emission distribution. This is to
be expected since the different latent components imply a slightly
different mapping to Hilbert space.

We next describe the kernel for various graphical models which
essentially impose a factorization on the joint probability
distribution over the $N$ variables $x_{\rm o}$ and $x_{\rm h}$. This
article focuses on directed graphs (although undirected graphs can
also be kernelized) where this factorization implies that the joint
distribution can be written as a product of conditionals over each
node or variable $x_i \in \{ x_{\rm o},x_{\rm h} \}, i=1\ldots N$
given its parent nodes or set of parent variables $x_{pa(i)}$ as
$p(x_{\rm o},x_{\rm h}) = \prod_{i=1}^N p(x_i|x_{pa(i)})$.  We now
discuss particular cases of graphical models including mixture models,
hidden Markov models and linear dynamical systems.

\subsection{Mixture Models}

In the simplest scenario, the latent graphical model could be a
mixture model, such as a mixture of Gaussians or mixture of
exponential family distributions. Here, the hidden variable $x_{\rm
h}$ is a single discrete variable which is a parent of a single
$x_{\rm o}$ emission variable in the exponential family. To derive the
full kernel for the mixture model, we first assume it is
straightforward to compute an elementary kernel between any pair of
entries (indexed via $x_{\rm h}$ and $x_{\rm h}'$) from each mixture
model as follows:
\[
{\tilde k}(p(x_{\rm o}|x_{\rm h}),p'(x_{\rm
o}|x_{\rm h}')) &=& \sum_{x_{\rm o}} \left ( p(x_{\rm o}|x_{\rm h})
p'(x_{\rm o}|x_{\rm h}') \right )^\rho.
\]
In the above, the conditionals could be, for instance, exponential
family (emission) distributions. Then, we can easily compute the
kernel for a mixture model of such distributions using these
elementary kernel evaluations and a weighted summation of all possible
pairings of components of the mixtures:
\[
{\tilde k}(p,p') &=& 
\sum_{x_{\rm o}} \sum_{x_{\rm h}} \left ( p(x_{\rm h}) p(x_{\rm o}|x_{\rm h}) \right )^\rho
  \sum_{x_{\rm h}'} \left ( p'(x_{\rm h}') p'(x_{\rm o}|x_{\rm h}') \right )^\rho  \\
&=&  \sum_{x_{\rm h},x_{\rm h}'} \left ( p(x_{\rm h})  p'(x_{\rm h}') \right )^\rho {\tilde k}(p(x_{\rm o}|x_{\rm h}),p'(x_{\rm
o}|x_{\rm h}')).
\]
Effectively, the mixture model kernel involves enumerating all
settings of the hidden variables $x_{\rm h}$ and $x_{\rm h}'$. Taking
the notation $| .|$ to be the cardinality of a variable, we
effectively need to perform and sum a total of $| x_{\rm h} | \times |
x_{\rm h}' |$ elementary kernel computations ${\tilde k}_{x_{\rm
h},x_{\rm h}'}$ between the individual emission distributions.

\subsection{Hidden Markov Models}

\begin{figure}
  \center
  \begin{tabular}[b]{ccc} 
\scalebox{0.16}{\rotatebox{270}{\includegraphics*[-0.2in,-0.2in][3.9in,11.7in]{hmm.prn}}} &
\scalebox{0.16}{\rotatebox{270}{\includegraphics*[-0.2in,-0.2in][3.9in,11.7in]{hmmp.prn}}} &
\scalebox{0.16}{\rotatebox{270}{\includegraphics*[-0.2in,-0.2in][5.9in,11.7in]{hmmp2.prn}}} \\
\COMMENT{ \includegraphics[height=1.2in,angle=270]{hmms2.eps} &
\includegraphics[height=1.2in,angle=270]{hmmu2.eps} &
\includegraphics[height=1.2in,angle=270]{hmmus2.eps} \\}
(a) HMM for $p$. & (b) HMM for $p'$. & (c) Combined graphs.
  \end{tabular}
\caption{Two hidden Markov models and the resulting graphical model as
the kernel couples common parents for each node creating undirected
edges between them.}
\label{fig:hmm}
\end{figure}

We next upgrade beyond mixtures to graphical models such as hidden
Markov models (HMMs). Considering hidden Markov models
as the generative model $p$ allows us to build a kernel over sequences
of variable lengths. However, HMMs and many popular Bayesian networks
have a large number of (hidden and observed) variables and enumerating
all possible hidden states for two distributions $p$ and $p'$ quickly
becomes inefficient since $x_{\rm h}$ and $x_{\rm h}'$ are
multivariate and/or have large cardinalities. However, unlike the
plain mixture modeling case, HMMs have an efficient graphical model
structure implying a factorization of their distribution which we can
leverage to compute our kernel efficiently. A hidden Markov model over
sequences of length $T+1$ has observed variables $x_{\rm
o}=\{ x_0,\ldots,x_{T} \}$ and hidden states $x_{\rm
h}=\{ q_0,\ldots,q_{T} \}$. The graphical model for an HMM in
Figure~\ref{fig:hmm}(a) reflects its Markov chain assumption and leads
to the following probability density function (note here we define
$q_{-1} = \{ \}$ as a null variable for brevity or, equivalently,
$p(q_0|q_{-1})=p(q_0)$):
\[
p(x_{\rm o}) &=& \sum_{x_{\rm h}} p(x_{\rm o},x_{\rm h}) \:\: = \:\:
\sum_{q_0} \ldots \sum_{q_T} \prod_{t=0}^T p(x_t|q_t) p(q_t|q_{t-1}).
\]

To compute the kernel ${\tilde k}(p,p')$ in a brute force manner, we
need to sum over all configurations of $q_0,\ldots,q_T$ and
$q_0',\ldots,q_T'$ while computing for each each configuration its
corresponding elementary kernel ${\tilde
k}(p(x_0,\ldots,x_T|q_0,\ldots,q_T),p'(x_0,\ldots,x_T|q_0',\ldots,q_T'))$. Exploring
each setting of the state space of the HMMs (shown in
Figure~\ref{fig:hmm}(a) and (b)) is highly inefficient requiring
$|q_t|^{(T+1)} \times |q_t'|^{(T+1)}$ elementary kernel evaluations. However, we
can take advantage of the factorization of the hidden Markov model by
using graphical modeling algorithms such as the junction tree
algorithm \citep{JordanNewBook} to compute the kernel efficiently.
First, we use the factorization of the HMM in the kernel as
follows:
\[
{\tilde k}(p,p') &=& 
\sum_{x_0,\ldots,x_T}
\sum_{q_0 \ldots q_T} \prod_{t=0}^T p(x_t|q_t)^\rho p(q_t|q_{t-1})^\rho
\sum_{q_0' \ldots q_T'} \prod_{t=0}^T p'(x_t|q_t')^\rho p(q_t'|q_{t-1}')^\rho \\
&=& \sum_{q_0 \ldots q_T} \sum_{q_0' \ldots q_T'} \prod_{t=0}^T p(q_t|q_{t-1})^\rho  p(q_t'|q_{t-1}')^\rho \left (
\sum_{x_t} p(x_t|q_t)^\rho  p'(x_t|q_t')^\rho 
\right ) \\
&=& \sum_{q_0 \ldots q_T} \sum_{q_0' \ldots q_T'} \prod_{t=0}^T p(q_t|q_{t-1})^\rho  p(q_t'|q_{t-1}')^\rho \psi (q_t,q_t') \\
&=& \sum_{q_T} \sum_{q_T'} \psi(q_T,q_T') \prod_{t=1}^{T} \sum_{q_{t-1}}\sum_{q_{t-1}'} p(q_{t}|q_{t-1})^\rho  p(q_{t}'|q_{t-1}')^\rho \psi (q_{t-1},q_{t-1}') p(q_0)^\rho p'(q_0')^\rho.
\]

\begin{figure}
  \center
  \begin{tabular}[b]{c} 
\scalebox{0.4}{\rotatebox{270}{\includegraphics*[0.0in,0.0in][3.5in,11.5in]{jta2.prn}}} 
  \end{tabular}
\caption{The undirected clique graph obtained from the kernel for
efficiently summing over both sets of hidden variables $x_{\rm h}$ and
$x_{\rm h}'$.}
\label{fig:jta}
\end{figure}


In the above we see that we only need to compute the kernel for each
setting of the hidden variable for a given time step on its own.
Effectively, the kernel couples the two HMMs via their common children
as in Figure~\ref{fig:hmm}(c). The kernel evaluations, once solved,
form non-negative clique potential functions $\psi(q_t,q_t')$ with a
different setting for each value of $q_t$ and $q_t'$. These couple the
hidden variable of each HMM for each time step. The elementary kernel
evaluations are easily computed (particularly if the emission
distributions $p(x_t|q_t)$ are in the exponential family) as:
\[
{\tilde k}(p(x_t|q_t),p'(x_t|q_t')) &=& \psi (q_t,q_t') \:\: = \:\: \sum_{x_t}
p(x_t|q_t)^\rho p'(x_t|q_t')^\rho.
\]
Only a total of $(T+1)\times |q_t| \times |q_t'|$ such elementary
kernels need to be evaluated and form a limited number of such clique
potential functions.  Similarly, each conditional distribution
$p(q_t|q_{t-1})^\rho$ and $p'(q_t'|q_{t-1}')^\rho$ can also be viewed
as a non-negative clique potential function, i.e. $\phi(q_t,q_{t-1}) =
p(q_t|q_{t-1})^\rho$ and
$\phi(q_t',q_{t-1}')=p'(q_t'|q_{t-1}')^\rho$. Computing the kernel
then involves marginalizing over the hidden variables $x_{\rm h}$,
which is done by running a forward-backward or junction-tree algorithm
on the resulting undirected clique tree graph shown in
Figure~\ref{fig:jta}. Thus, to compute the kernel for two sequences
${\pzc x}$ and ${\pzc x}'$ of lengths $T_{\pzc x}+1$ and $T_{{\pzc
x}'}+1$, we train an HMM from each sequence using maximum likelihood
and then compute the kernel using the learned HMM transition matrix
and emission models for a user-specified sequence length $T+1$ (where
T can be chosen according to some heuristic, for instance, the average
of all $T_{\pzc x}$ in the training data set). The following is
pseudo-code illustrating how to compute the ${\tilde k}(p,p')$ kernel
given two hidden Markov models $(p,p')$ and user-specified parameters
$T$ and $\rho$.
\[
&&\Phi(q_0,q_0') = p(q_0)^\rho p'(q_0')^\rho  \\
&& {\rm for} \:\: t=1\ldots T \\
&&\:\:\:\: \:\:\:\: \Phi(q_t,q_t') = \sum_{q_{t-1}} \sum_{q_{t-1}'} p(q_t|q_{t-1})^\rho  p'(q_t'|q_{t-1}')^\rho \psi(q_{t-1},q_{t-1}') \Phi(q_{t-1},q_{t-1}') \\
&& {\rm end} \\
&&{\tilde k}(p,p') = \sum_{q_T} \sum_{q_T'} \Phi(q_T,q_T') \psi(q_T,q_T').
\]
Note that the hidden Markov models $p$ and $p'$ need not have the same
number of hidden states for the kernel computation.


\subsection{Bayesian Networks}

The above method can be applied to Bayesian networks in general where
hidden and observed variables factorize according to
$p(x)=\prod_{i=1}^N p(x_i|x_{pa(i)})$. The general recipe mirrors the
above derivation for the hidden Markov model case. In fact, the two
Bayesian networks need not have the same structure as long as they
share the same sample space over the emission variables $x_{\rm
o}$. For instance, consider computing the kernel between the Bayesian
network in Figure~\ref{fig:bayesnet}(a) and the hidden Markov model in
Figure~\ref{fig:hmm}(b) over the same sample space. The graphs can be
connected with their common children as in
Figure~\ref{fig:bayesnet}(b). The parents with common children are
married and form cliques of hidden variables. We then only need to
evaluate the elementary kernels over these cliques giving the
following non-negative potential functions for each observed node (for
instance $x_i \in x_{\rm o}$) under each setting of all its parents'
nodes (from both $p$ and $p'$):
\[
\psi(x_{{\rm h},{\rm pa}(i)},x_{{\rm h}',{\rm pa}(i)}') &=& {\tilde k}(p(x_i|x_{{\rm
h},{\rm pa}(i)}),p'(x_i|x_{{\rm h}',{\rm pa}(i)'})) \\
&=& \sum_{x_i} p(x_i|x_{{\rm
h},{\rm pa}(i)})^\rho p'(x_i|x_{{\rm h}',{\rm pa}(i)'})^\rho.
\]
After marrying common parents, the resulting clique graph is then
built and used with a junction tree algorithm to compute the overall
kernel from the elementary kernel evaluations.  Although the resulting
undirected clique graph may not always yield a dramatic improvement in
efficiency as was seen in the hidden Markov model case, we note that
propagation algorithms (loopy or junction tree algorithms) are still
likely to provide a more efficient evaluation of the kernel than the
brute force enumeration of all latent configurations of both Bayesian
networks in the kernel (as was described in the generic mixture model
case).  While the above computations emphasized discrete latent
variables, the kernel is also computable over continuous latent
configuration networks, where $x_{\rm h}$ contains scalars and vectors,
as is the case in linear Gaussian models, which we develop 
in the next subsections.

\begin{figure}
  \center
  \begin{tabular}[b]{ccc} 
\scalebox{0.22}{\rotatebox{270}{\includegraphics*[-0.2in,-0.2in][3.9in,11.7in]{bayesnet.prn}}} &
\scalebox{0.22}{\rotatebox{270}{\includegraphics*[-0.2in,-0.2in][4.9in,11.7in]{bayesnet2.prn}}} \\
(a) Bayesian network for $p$. & (b) Combined graphs. 
  \end{tabular}
\caption{The resulting graphical model from a Bayesian network and a hidden Markov model
as the kernel couples common parents for each node creating undirected
edges between them and a final clique graph for the junction tree algorithm.}
\label{fig:bayesnet}
\end{figure}

\subsection{Linear Gaussian Models}\label{sec: LG}

\newcommand{\xhid}{x_\text{h}}
\newcommand{\xobs}{x_\text{o}}

A linear Gaussian model on a  
directed acyclic graph \m{\mc{G}} with variables \m{\seq{x}{N}} associated to its 
vertices is of the form 
\begin{equation*}\label{eq: linear guassian} 
p\br{\seq{x}{N}}=
\prod_{i} \mc{N} \br{\:x_i\:;\: \beta_{i,0} + 
\! {\textstyle \sum_{j\in \text{pa}(i)}} \beta_{i,j}x_j,\, \Sigma_i\:}
\end{equation*}
where $\mc{N}(x;\mu,\Sigma)$ is the multivariate Gaussian distribution, and $\text{pa}(i)$ is 
the index set of parent vertices associated with the $i^\text{th}$ vertex. The unconditional joint distribution can be recovered from the conditional form and is itself a Gaussian distribution over the variables in the model \m{p\br{\seq{x}{N}}}.  Hence, the probability 
product kernel between a pair of linear Gaussian models can be computed in closed form 
by using the unconditional joint distribution and \rf{eq: kernel between Gaussians}.

%The probability product kernel 
%can be expressed in closed form for any choice of $\rho$ by 
%recovering the joint Gaussian distribution and applying the formula for the kernel 
%between Gaussian distributions:
%\begin{equation*}
%k(p,p') = \int \prod_{i} p(x_i|x_{\text{pa}(i)})^\rho %\prod_{i} 
%p'(x_i|x_{\text{pa}(i)})^\rho\, dx \\
%&=& \int p(x)^\rho p'(x)^\rho\, dx
%\end{equation*}
%can be computed 

Latent variable models require an additional marginalization step. % to compute. 
Variables are split into two sets: the observed variables \m{\xobs}
and the hidden variables \m{{\xhid}}. 
After the joint unconditional distribution $p(\xobs,\xhid)$ has been computed, 
the hidden variables can be trivially integrated out: 
%with some additional manipulation.
\begin{eqnarray*}
 % k(p,p') &=& \int p(\xobs,\xhid)^\rho\, 
  %p'(\xobs,x'_\text{h})^\rho\, d\xobs\, d\xhid\, dx'_\text{h} \\
  %&=& \int  p(\xobs)^\rho\,p(\xhid|\xobs)^\rho 
  %p'(\xobs)^\rho\,p'(x'_{\text{h}}|\xobs)^\rho\, d\xobs\, d\xhid\, dx'_{\text{h}} \\
  %&=& \sqbP{ (2\pi)^{(D+D')}\, \absbig{\Sigma_{\xhid|\xobs}}\, 
  %  \absbig{\Sigma'_{\xhid|x'_o}} }^\frac{1-\rho}{2} \rho^\frac{-(D+D')}{2} 
  %\int p(\xobs)^\rho\, p'(\xobs)^\rho\, d\xobs,
 k(p,p') &=& \int \left( \int p(\xobs,\xhid) \,d\xhid \right)^\rho\,\left ( \int p(\xobs,x'_\text{h})\, dx'_\text{h}\right )^\rho\, d\xobs \:\: = \:\:
\int p(\xobs)^\rho\, p'(\xobs)^\rho\, d\xobs,
\end{eqnarray*}
so that the kernel is computed between Gaussians over only the observed variables.
%where $|\Sigma_{\xhid|\xobs}|$ is the determinant of the covariance of the conditional 
%distribution $p(\xhid|\xobs)$ and $D$ is the number of latent variables. 
%Note that when $\rho\<=1$ this additional step can be skipped because the conditional distribution 
%over the latent variables integrates to unity.
\subsection{Linear Dynamical Systems}
A commonly used linear Gaussian model with latent variables is the
linear dynamical system (LDS) \citep{ShuSto82}, also known as the
state space model or Kalman filter model.  The LDS is the continuous
state analog of the hidden Markov model and shares the same
conditional independence graph Figure~\ref{fig:hmm}(a). Thus, it is
appropriate for modeling continuous time series data and continuous
dynamical systems. The LDS's joint probability model is of the form
%with all variables distributed as Linear Gaussians.
\begin{equation*}
p(x_0,\ldots x_T,s_0,\ldots s_T)=
\mc{N}(s_0;\mu,\Sigma)\, \mc{N}(x_0; Cs_0,R) 
\prod_{t=1}^T \mc{N}(s_t;As_{t-1},Q)\,\mc{N}(x_t; Cs_t,R),
\end{equation*}
where $x_t$ is the observed variable at time $t$, and $s_t$ is the latent state space 
variable at time $t$, obeying the Markov property.

The probability product kernel between LDS models is 
\begin{equation*}
 % k(p,p') = \sqbP{ (2\pi)^{(D+D')}\, \abs{\Sigma_{s|x}}\, 
 %   \absbig{\Sigma'_{s|x}} }^\frac{1-\rho}{2} 
 % \rho^\frac{-(D+D')}{2} \mc{N}(x;\mu_x,\Sigma_{xx})^\rho\, \int\mc{N}(x;\mu'_x,\Sigma'_{xx})^\rho\,dx
 k(p,p') = \int\mc{N}(x;\mu_x,\Sigma_{xx})^\rho\,\mc{N}(x;\mu'_x,\Sigma'_{xx})^\rho\,dx,
\end{equation*}
where $\mu_x$ and $\Sigma_{xx}$ are the unconditional mean and covariance, which can be computed from the recursions
%The relevant terms to evaluate the kernel can be computed from the recursions 
\begin{eqnarray*}
\mu_{s_t} &=& A\mu_{s_{t-1}}\\
\mu_{x_t} &=& C\mu_{s_t}\\
\Sigma_{s_ts_t} &=& A\Sigma_{s_{t-1}s_{t-1}}A' + Q\\
\Sigma_{x_tx_t} &=& C\Sigma_{s_ts_t}C' + R. 
%\Sigma_{x_ts_t} &=& C \Sigma_{s_ts_t}\\
%\Sigma_{s_ts_{t-1}} &=& A \Sigma_{s_{t-1}s_{t-1}}. 
\end{eqnarray*}
%The joint covariance is reconstructed from the smaller covariance terms 
%\begin{eqnarray*}
%\Sigma = \left[ \begin{array}{cc}
%\Sigma_{ss} & \Sigma_{sx}  \\
%\Sigma_{xs} & \Sigma_{xx}  \end{array} \right],
%\end{eqnarray*}
%although it is not necessary to do this explicitly.
%The final term to be computed is 
%\begin{equation*}
%\Sigma_{s|x} = \Sigma_{ss} - \Sigma_{sx}\Sigma_{xx}^{-1}\Sigma_{xs}. 
%\end{equation*}

As with the HMM, we need to set the number of time steps before computing the kernel. 
Note that the most algorithmically expensive calculation when computing the probability product kernel between Gaussians is taking the inverse of the covariance matrices. 
The dimensionality of the covariance matrix will grow linearly with the number 
of time steps and the running time of the matrix inverse grows cubically with the dimensionality. 
However, the required inverses can be efficiently computed because $\Sigma_{xx}$ is block diagonal. 
Each extra time step added to the kernel will simply add another block,  
and the inverse of a block diagonal matrix can be computed by inverting the 
blocks individually. 
Therefore, the running time of inverting the block diagonal covariance matrix will only grow linearly with the 
number of time steps $T$ and cubically with the (small) block size.

\section{Sampling Approximation}
%In the context of generative modelling, complicated models are often necessary to 
%sufficiently describe the data of interest. 
%However, closed form solutions for a probability product kernel for any setting of 
%$\rho$ cannot always be found for these complex models. 
To capture the structure of real data, generative models often need to be quite intricate. 
Closed form solutions for the probability product kernel are then unlikely to exist,  
forcing us to rely on approximation methods to compute \m{k(p,p')}.  

Provided that evaluating the likelihood  \m{p(x)} is computationally 
feasible, and that  we can also efficiently sample from  
the model, in the \m{\rho\<=1} case (expected likelihood kernel) we may    
employ the Monte Carlo estimate  
%We must then rely upon efficient approximation methods to compute the kernel. 
%Monte Carlo sampling methods can be used to approximate the expected likelihood kernel 
%($\rho=1$). 
%This requires that we are able to efficiently sample and evaluate the likelihood of 
%these samples for a given generative model. In most cases this will be true. 
%The approximate kernel can be computed as
\[
k(p,p')  &\approx&  
\frac{\beta}{N} 
\Sum{i}{N}
%\sum_{\begin{array}{cc} \scriptstyle x^i \sim p(x) \\ 
%\scriptstyle i=1,\ldots,N \end{array}} 
p'(x^i) \:\: + \:\:
\frac{(1-\beta)}{N'}  
\Sum{i}{N}
%\sum_{\begin{array}{cc} {\scriptstyle x'^i \sim p'(x)} \\ 
%{\scriptstyle i=1,\ldots,N'} \end{array}}
 p(x'^{\,i}),
\]
where $x^1,\ldots,x^N$ and $x'^{1},\ldots,x'^{N}$ 
are i.i.d. samples from $p$ and $p'$ respectively,  
and $\beta\tin[0,1]$ is a parameter of our choosing. 
By the law of large numbers, 
this approximation will converge to the true kernel value in the limit \m{N\to\infty} 
for any \m{\beta}. 
Other sampling techniques such as importance sampling and a 
variety of flavors of Markov chain Monte Carlo (MCMC), including Gibbs sampling, 
can be used to improve the rate of the sampling approximation's convergence to the true kernel. 

In some cases it may be possible to use the sampling approximation for an general $\rho$. 
This is possible if a normalizing factor $Z$ can be computed to renormalize $p(x)^\rho$ into a valid distribution. 
$Z$ is computable for most discrete distributions and some continuous distribution, such as the Gaussian. 
The sampling approximation then becomes
\[
k(p,p')  &\approx&  
\frac{\beta}{N} 
\Sum{i}{N}
Z\:p'(\hat{x}^i)^\rho \:\: + \:\:
\frac{(1-\beta)}{N'}  
\Sum{i}{N}
 Z'\:p(\hat{x}'^{\,i})^\rho,
\]
where $Z$ and $Z'$ are the normalizers of $p$ and $p'$ after they are taken to the power of $\rho$. In this formulation $\hat{x}^1,\ldots,\hat{x}^N$ and $\hat{x}'^{1},\ldots,\hat{x}'^{N}$ 
are i.i.d. samples from $\hat{p}$ and $\hat{p}'$ where $\hat{p}= \frac{p^\rho}{Z}$ and $\hat{p}'= \frac{p'^\rho}{Z'}$. 

In practice, for a finite number of samples, the approximate kernel may not be a valid Mercer kernel. 
The non-symmetric aspect of the approximation can be alleviated by computing only the lower 
triangular entries in the kernel matrix and replicating them in the upper half. 
Problems associated with the approximation not being positive definite can be rectified by either using more samples, 
or by using an implementation for the kernel based classifier algorithm that can converge to a local optimum such as sequential minimal optimization \citep{Platt}.  
\subsection{Switching Linear Dynamical Systems}
In some cases, part of the model can be evaluated by sampling, while the rest is  
computed in closed form. 
This is true for the switching linear dynamical system (SLDS) 
\citep{PavReh00}, which combines the discrete hidden states of the HMM and the continuous 
state space of the LDS:   
%The joint distribution is defined as:
\begin{equation*}
p(x,s,q)=p(q_0)\,p(s_0|q_0)\,p(x_0|s_0)
\prod_{t=1}^{T} p(q_t|q_{t-1})\,p(s_t|s_{t-1},q_t)\,p(x_t|s_t), 
\end{equation*}
where $q_t$ is the discrete hidden state at time \m{t}, 
$s_t$ is the continuous hidden state,  
and $x_t$ are observed emission variables (Figure~\ref{fig:slds}). 
\begin{figure}
  \center
\scalebox{.180}{\rotatebox{270}{\includegraphics*[0.7in,0.5in][8in,11in]{SLDS2.eps}}}  \\
\caption{Graph for the Switching Linear Dynamical System.}
\label{fig:slds}
\end{figure}
Sampling \m{q^1\,\ldots,q^N} according to \m{p(q)} and 
\m{q'^1,\ldots,q'^N} according to \m{p'(q')}, the remaining factors in 
the corresponding sampling approximation %to the probability product kernel 
\begin{equation*}
%k(p,p') \approx 
\frac{1}{N} 
%\!\!\!\!\!\! 
%\sum_{ \begin{array}{c} 
%\scriptstyle q^i\sim p(q)\\
%\scriptstyle q'^i\sim p'(q')\\
%\scriptstyle i=1,\dots,N  
%\end{array} }
%{\br{q_i}_{i=1}^N, (q'_i)_{i=1}^N}
%\!\!\!\!\!
\Sum{i}{N}
\int \left( \int p(s_0|q_0^i)\, p'(s'_0|q'^i_0) 
\prod_{t=1}^{T}p(s_t|s_{t-1},q_t^i)\, p'(s'_t|s'_{t-1},q'^i_t) 
\prod_{t=0}^{T}p(x_t|s_t)\, p'(x_t|s'_t) \:ds\ts ds'\right)^\rho\! dx
\end{equation*} 
form a non-stationary linear dynamical system. Once the joint distribution is recovered, it
can be evaluated in closed form by \rf{eq: kernel between Gaussians}
as was the case with a simple LDS. 
For $\rho\neq1$ this is a slightly different formulation than the sampling approximation in the previous section. 
This hybrid method has the nice property that it is always a valid kernel for any number of samples since 
it simply becomes a convex combination of kernels between LDSs.
An SLDS model is useful for
describing input sequences that have a combination of discrete and
continuous dynamics, such as a time series of human motion containing
continuous physical dynamics and kinematics variables as well as
discrete action variables \citep{PavReh00}.
%Once the $q_t$'s are known, the kernel 
%reduces to a kernel between linear Gaussians.
%which can be solved in closed form.

\section{Mean Field Approximation}
\label{sec: mean field}

Another option when the probability product kernel between graphical models cannot be computed 
in closed form is to use variational bounds, 
such as the mean field approximation \citep{Jaa00}.  
In this method, we introduce a variational distribution \m{Q(x)}, and using Jensen's inequality, lower bound the kernel:
%the kernel value is lower bounded in the form  
%In this technique the so-called potential \m{\Psi(x)=p(x)p'(x)} is approximated 
%iteratively by introducing 
 %Another approximation method useful when computing the kernel in closed form is either 
%impossible or intractable is a form of the variational inference technique known as 
%mean field \citep{Jaa00}. 
%In this method we compute a lower bound on the kernel value and transform the problem 
%into computing sufficient statistics for a simpler more tractable distribution. 
%We first define the potential $\Psi(x)$ as the product of our two models. 
%This factors over a graph in the same form as the original distributions. 
%Once we've defined the potential, 
%we then introduce a simpler variational distribution Q(x) as follows:  
\begin{eqnarray*}
k(p,p')&=&\int p(x)^\rho\ts p'(x)^\rho\,dx = \int \Psi(x)^\rho\,dx = 
\exp{\left(\log{\int \frac{Q(x)}{Q(x)}\,\Psi(x)^\rho\, dx}\right)} \\
%&=& \exp{\left(\log{\int \Psi(x)^\rho dx}\right)} \\
&\geq& \exp{\left(\int Q(x)\log{\frac{\Psi(x)^\rho}{Q(x)}}\:dx\right)} = 
\exp{\left( \rho\, \E_Q\sqb{\log{ \Psi(x)}} + H(Q)\right) } = B(Q), 
\end{eqnarray*}
where \m{\Psi(x)\<=p(x)\ts p'(x)} is called the potential, 
$\E_Q\sqb{\:\cdot\:}$ denotes the expectation with respect to $Q(x)$,  
and $H(Q)$ is the entropy. This transforms the problem from evaluating an intractable integral into that of finding the sufficient statistics of $Q(x)$ and computing the subsequent tractable expectations.
%$B(Q)$ is a valid lower bound for any choice of $Q(x)$.
Note, using the mean field formulation gives a different (and arguably
closer result) to the desired intractable kernel than simply computing
the probability product kernels directly between the approximate
simple distributions. It is interesting to note the following form of
the bound when it is expressed in terms of the Kullback-Leibler
divergence, $D(Q \|p) \: = \: \int Q(x) \log \frac{Q(x)}{p(x)} dx$,
and an entropy term. The bound on the probability product kernel is a
function of the Kullback-Leibler divergence to $p$ and $p'$:
\[
%B(Q) &=& \exp(-\frac{1}{2} D(Q \|p) - \frac{1}{2}D(Q \|p')) \:\: {\rm
%where} \:\: D(Q \|p) \:\: =  \:\: \int Q(x) \log \frac{Q(x)}{p(x)} dx
k(p,p') \:\: \geq \:\: B(Q)  \:\: = \:\: \exp\left(-\rho\: D(Q \|p) - \rho\: D(Q \|p') + (1-2\rho)\: H(Q)\right) .
\label{eq: boundkl}
\]

The goal is to define $Q(x)$ in a way that makes computing the sufficient statistics easy. 
%The graph structure of $Q(x)$ is chosen so that calculation of the sufficient 
%statistics is computationally easy. 
When the graphical structure of \m{Q(x)} has no edges, which 
corresponds to the case that the underlying variables are independent, 
this is called the mean field method.  
When $Q(x)$ is made up of a number of more elaborate independent structures, 
it is called the structured mean field method. 
In either case, the variational distribution factors in the form 
\m{Q_1(x_1) Q_2(x_2)\ldots Q_N(x_N)} with \m{\seq{x}{N}} disjoint subsets 
of the original variable set \m{x}. 

Once the structure of \m{Q(x)} has been set, its (locally) optimal parameterization  
is found by %iteratively maximizing the lower bound $B(Q)$. 
%Hence, $B(Q)$ can be viewed as a function of $Q_n(x_n)$, the distribution over the 
%variables $x_n$ associated with the $n^\text{th}$ independant substructure. 
%The bound is optimized by taking the derivative with respect to $Q(x_n)$:
cycling through the distributions \m{Q_1(x_1), Q_2(x_2),\ldots,Q_N(x_N)} and 
maximizing the lower bound $B(Q)$ with respect to each one, holding the others 
fixed:   
\begin{eqnarray*}
\frac{\partial \log{B(Q)}}{\partial\, Q_n(x_n)} 
= \frac{\partial}{\partial\, Q_n(x_n)} 
\sqbP{ \rho \int Q_n(x_n)\ts \E_Q \sqb{\log{ \Psi(x)}|x_n} \ts dx_n + H(Q_n) + \text{constant}} = 0.  \\
 %&=& \rho E_Q{\{\log{ \Psi(x)}|x_n\}} + log(Q_n) + 1 = 0 \\
\end{eqnarray*}
%Solving for $Q_n(x_n)$ we get: 
This leads to the update equations 
\begin{eqnarray*}
 Q_n(x_n) = \frac{1}{Z_n}\exp{\left(\rho\, \E_Q\sqb{\log\Psi(x)|x_n}\right)}\,,\\
\end{eqnarray*}
where the conditional expectation $\E_Q\sqb{\:\cdot\:|\,x_n}$ is computed by taking the 
expectation over all variables except \m{x_n} and 
\begin{equation*}
 Z_n = \int \exp{\left(\rho\,  \E_Q\sqb{\log\Psi(x)|x_n} \right)}\: dx   
\end{equation*}
is the normalizing factor guaranteeing that \m{Q(x)} remains a valid distribution. 
%This gives a set of fixed point equations. 
%Each $Q(x_n)$ is updated one at a time with the others held fixed. 
%The mean field equations are iteratively updated until the bound on the kernel 
%value converges.

In practice this approximation may not give positive definite (PD) kernels. 
Much current research has been focused on kernel algorithms using non-positive definite kernels. 
Algorithms such as sequential minimal optimization \citep{Platt} can converge to locally optimal solutions. 
The learning theory for non-PD kernels is currently of great interest. See \cite{Ong04} for an interesting discussion and additional references.
\subsection{Factorial Hidden Markov Models}

\begin{figure}
  \center
  \begin{tabular}[b]{ccc} 
\scalebox{0.25}{\rotatebox{270}{\includegraphics*[0.6in,2.0in][8in,9in]{FHMM2.eps}}} &
\scalebox{0.25}{\rotatebox{270}{\includegraphics*[0.6in,2.0in][8in,9in]{FHMMjoined.eps}}} &
\scalebox{0.25}{\rotatebox{270}{\includegraphics*[0.6in,2.0in][8in,9in]{FHMMmf.eps}}} \\
(a) FHMM $P(x,q)$. & (b) Potential $\Psi(x,q,q')$. & (c) Variational $Q(x,q,q')$.
  \end{tabular}
\caption{Graphs associated with the factorial hidden Markov model, its probability product kernel potential, 
and its structured mean field distribution.}
\label{fig:fhmm}
\end{figure}

One model that traditionally uses a structured mean field approximation for inference 
is the factorial hidden Markov model (FHMM)\citep{ZoubinJordan}. Such
a model is appropriate for modeling output sequences that are emitted
by multiple independent hidden processes. The FHMM is given by %defined as
%follows:
\begin{equation*}
p(x_0,\seq{x}{T})=\prod_{c=1}^C \,\sqbP{p(q_0^c)\, 
\prod_{t=1}^T\, p(q_t^c|q_{t-1}^c)  }\,
\prod_{t=0}^T \, p(x_t|q_t^1,\dots,q_t^ C), 
\end{equation*}
where $q_t^c$ is the factored discrete hidden state at time $t$ for chain $c$ and $x_t$ is 
the observed Gaussian emission variable at time $t$. 
The graphs of  the joint FHMM distribution, 
the potential $\Psi(x,q,q')$, and the variational distribution $Q(x,q,q')$ 
are depicted in Figure~\ref{fig:fhmm}. 
%The graph representing the joint FHMM distribution, the potential $\Psi(x,q,q')$ and the 
%simpler variational distribution $Q(x,q,q')$ are in Figure~\ref{fig:fhmm} 
%a, b and c respectively.
The structured mean field approximation is parameterized as %as follows:
\begin{eqnarray*}
Q(\,q_0^c \<= i\,) &\propto& \exp{\left(\rho \log \phi^c(i) - \frac{\rho}{2} \log{\abs{\Sigma_i}} - 
\frac{\rho}{2} \E_Q\sqb{(x_0\<-\mu_i)^T\Sigma_i^{-1}(x_0\<-\mu_i)} \right)} \\
Q(\,q_t^c \<= i\,|\,q_{t-1}^c \<= j\,) &\propto& 
\exp{\left(\rho \log \Phi^c(i,j) - \frac{\rho}{2} \log{\abs{\Sigma_i}} - 
\frac{\rho}{2} \E_Q\sqb{(x_t\<-\mu_i)^T\Sigma_i^{-1}(x_t\<-\mu_i)}\right)}\\
Q(x_t) &=& \mc{N}(x_t; \hat{\mu}_t, \hat{\Sigma}_t) \\
\hat{\Sigma}_t &=& \ovr{\rho}\:  
\sqbbigg{\: \sum_{i,c} \E_Q\sqb{s_t^c(i)} \Sigma_i^{-1} + 
\sum_{i,c} \E_Q\sqb{q'^c_t(i)} \Sigma'^{-1}_i\:}^{-1}\\
\hat{\mu}_t &=& \rho\, \hat{\Sigma}_t 
\sqbbigg{\: \sum_{i,c} \E_Q\sqb{q_t^c(i)} \Sigma_i^{-1} \mu_i + 
\sum_{i,c} \E_Q\sqb{q'^c_t(i)} \Sigma'^{-1}_i \mu'_i\:}.
\end{eqnarray*}

It is necessary to compute the expected sufficient statistics to both update the mean 
field equations and to evaluate the bound. 
The expectations $\E_Q\sqb{x_t}$ and $\E_Q\sqb{x_t x^T_t}$ can easily be computed 
from $Q(x_t)$, whereas 
$Q(\,s_t^c \<= i\,)$ and $Q(\,s_t^c \<= i\,,\: s_{t-1}^c \<= j\,)$ 
can be computed efficiently by means of a junction tree algorithm on each 
independent chain in the approximating distribution.

\section{Relationship to Other Probabilistic Kernels}
\label{sec: relation}

While the probability product kernel is not the only kernel to
leverage generative modeling, it does have advantages in terms of
computational feasibility as well as in nonlinear flexibility.  For
instance, heat kernels or diffusion kernels on statistical manifolds
\citep{LafLeb02} are a more elegant way to generate a kernel in a
probabilistic setting, but computations involve finding geodesics on
complicated manifolds and finding closed form expressions for the heat
kernel, which are only known for simple structures such as spheres and
hyperbolic surfaces. Figure~\ref{fig:manifold} depicts such a
statistical manifold. The heat kernel can sometimes be approximated
via the geodesic distance from $p$ which is parameterized by $\theta$
to $p'$ parameterized by $\theta'$. But, we rarely have
compact closed-form expressions for the heat kernel or for these
geodesics for general distributions. Only Gaussians with spherical
covariance and multinomials are readily accommodated. Curiously, the
probabilistic product kernel expression for both these two
distributions seems to coincide closely with the more elaborate heat
kernel form. For instance, in the multinomial case, both involve an
inner product between the square-rooted count frequencies. However,
the probability product kernel is straightforward to compute for a
much wider range of distributions providing a practical alternative to
heat kernels.

\begin{figure}
  \center
  \begin{tabular}[b]{c} 
\scalebox{0.27}{\rotatebox{270}{\includegraphics*[2in,0.0in][7in,11.5in]{manifold.prn}}}
  \end{tabular}
\caption{A statistical manifold with a geodesic path between two
generative models $\theta$ and $\theta'$ as well as a local tangent
space approximation at, for instance, the maximum likelihood model
$\theta^*$ for the aggregated data set.}
\label{fig:manifold}
\end{figure}


Another probabilistic approach is the Fisher kernel \citep{JaaHau98}
which approximates distances and affinities on the statistical
manifold by using the local metric at the maximum likelihood estimate
$\theta^*$ of the whole data set as shown in
Figure~\ref{fig:manifold}. One caveat, however, is that this
approximation can produce a kernel that is quite different from the
exact heat kernel above and may not preserve the interesting
nonlinearities implied by the generative model. For instance, in the
exponential family case, all distributions generate Fisher kernels
that are linear in the sufficient statistics. Consider the Fisher
kernel
\[
k(\pzc x,\pzc x') &=&  U_{\pzc x} I^{-1}_{\theta^*} U_{\pzc x '}
\]
where $I^{-1}_{\theta^*}$ is the inverse Fisher information matrix at
the maximum likelihood estimate and where we define
\[
U_{\pzc x} &=&\nabla_\theta \log p(\pzc x|\theta)|_{\theta^*}.
\]
In exponential families, $p_\theta(x) = \exp({\cal A}(x)+\theta^T
\mc{T}(x)-{\cal K}(\theta))$ producing the following $ U_{\pzc x} =
\mc{T}(\pzc x)-{\cal G}(\theta^*)$ where we define ${\cal G}(\theta^*)
= \nabla_\theta {\cal K}(\theta))|_{\theta^*}$. The resulting Fisher
kernel is then
\[
k(\pzc x,\pzc x') &=& \left ( \mc{T}(\pzc x) -{\cal G}(\theta^*)
\right)^T I^{-1}_{\theta^*} \left ( \mc{T}({\pzc x}')^T -{\cal G}(\theta^*) \right ),
\]
which is an only slightly modified form of the linear kernel in
$\mc{T}(\pzc x)$. Therefore, via the Fisher kernel method, a Gaussian
distribution over means generates only linear kernels. A Gaussian
distribution over means and covariances generates quadratic
kernels. Furthermore, multinomials generate log-counts unlike the
square root of the counts in the probability product kernel and the
heat kernel. In practical text classification applications, it is
known that square root squashing functions on frequencies and count
typically outperform logarithmic squashing functions \citep{GolSah98,
CutKarPedTuk92}. In \citep{TsuKinAsa02,KasTsuIno03}, another related
probabilistic kernel is put forward, the so-called marginalized
kernel. This kernel is again similar to the Fisher kernel as well as
the probability product kernel. It involves marginalizing a kernel
weighted by the posterior over hidden states given the input data for
two different points, in other words the output kernel $k(x,x')=\sum_h
\sum_{h'} p(h|x) p(h'|x') k((x,h),(x',h'))$. The probability product
kernel is similar, however, it involves an inner product over both
hidden variables and the input space $x$ with a joint distribution.

A more recently introduced kernel, the exponentiated symmetrized
Kullback-Leibler (KL) divergence \citep{MorHoVas04} is also related to
the probability product kernel. This kernel involves exponentiating the negated
symmetrized KL-divergence between two distributions
\[
k(p,p') &=& \exp( - \alpha D(p \| p') - \alpha D(p' \| p)  + \beta )
\:\: {\rm where} \:\: D(p\|p') \:\: = \:\: \int_x p(x) \log
\frac{p(x)}{p'(x)} dx
\]
where $\alpha$ and $\beta$ are user-defined scalar parameters.
However, this kernel may not always satisfy Mercer's condition and a
formal proof is not available. Another interesting connection is that
this form is very similar to our mean-field bound on the probability
product kernel \rf{eq: boundkl}. However, unlike the probability
product kernel (and its bounds), computing the KL-divergence between
complicated distributions (mixtures, hidden Markov models, Bayesian
networks, linear dynamical systems and intractable graphical models)
is intractable and must be approximated from the outset using
numerical or sampling procedures which further decreases the
possibility of having a valid Mercer kernel \citep{MorHoVas04}.

\COMMENT{
Another advantage to the probability product kernel's mappings ${\pzc
x} \rightarrow p(x)^\rho$ is that they are almost as general as the
mapping to an arbitrary Hilbert space. In other words, we can mimic
${\pzc x} \rightarrow \Phi({\pzc x})$ even though we have the
restriction $p(x) \geq 0$ and $\int_x p(x) dx = 1$. Consider
converting any $\Phi({\pzc x})$ and its corresponding kernel $k({\pzc
x},{\pzc x}')= \Phi({\pzc x})^T\Phi({\pzc x}')$ to a probability
mapping by making its features non-negative and normalized. This can
be done by adding a constant to all feature entries and dividing by
the total sum of features
\[
p(x|{\pzc x}) &=& \frac{ \Phi({\pzc x})+\Phi_c }{\sqrt{ \Phi({\pzc x})^T \Phi({\pzc x}) + 2 \Phi_c^T \Phi({\pzc x}) + \Phi_c^T \Phi_c  } }
\]
In the $\rho=1$ case, such a mapping produces the surrogate kernel
\[
k({\pzc x},{\pzc x}')^* &=& \int_x p(x|{\pzc x})p(x|{\pzc x}') dx \\
&=& \frac{k({\pzc x},{\pzc x}') + c {\vec 1}^T \Phi({\pzc x}) + c {\vec 1}^T \Phi({\pzc x}') + c^2 {\vec 1}^T {\vec 1}}
{
{\vec 1}^T\Phi({\pzc x}) {\vec 1}^T\Phi({\pzc x}')
+ c{\vec 1}^T{\vec 1}{\vec 1}^T\Phi({\pzc x})
+ c{\vec 1}^T{\vec 1}{\vec 1}^T\Phi({\pzc x}')
+ c^2{\vec 1}^T{\vec 1}{\vec 1}^T{\vec 1}
} .
\]
For a support vector machine (SVM) whose decision boundary is given by
this new kernel which has only non-negative features we have:
\[
f({\pzc x}) &=& \sum_t \alpha_t y_t k({\pzc x},{\pzc x}_t)^* + b \\
&=& \sum_t \alpha_t y_t k({\pzc x},{\pzc x}_t) + c {\vec 1}^T \Phi({\pzc x}) + b^*
\]
where the new bias term $b^*$ absorbs the extra constants.
}

\section{Experiments}

In this section we discuss three different learning problems: text
classification, biological sequence classification and time series
classification. For each problem, we select a generative model that is
compatible with the domain which then leads to a specific form for the
probability product kernel. For text, multinomial models are used, for
sequences hidden Markov models are natural and for time series data we
use linear dynamical systems. The subsequent kernel computations are
fed to a discriminative classifier (a support vector machine) for
training and testing.

\subsection{Text Classification}

For text classification, we used the standard WebKB data set which
consists of HTML documents in multiple categories. Only the text
component of each web page was preserved and HTML markup information
and hyper-links were stripped away.  No further stemming or elaborate
text pre-processing was performed on the text. Subsequently, a
bag-of-words model was used where each document is only represented by
the frequency of words that appear within it. This corresponds to a
multinomial distribution over counts. The frequencies of words are the
maximum likelihood estimate for the multinomial generative model to a
particular document which produces the vector of parameters $\hat
\alpha$ as in Section~\ref{sec: exp fam}.

\begin{figure}[htbp]
  \center
  \begin{tabular}[b]{cc}
\includegraphics[height=2.2in]{webkb77.eps} &
\includegraphics[height=2.2in]{webkb622.eps} \\
(a) Training Set Size = 77 & (b) Training Set Size = 622
  \end{tabular}
\caption[ ]{SVM error rates (and standard deviation) for the
probability product kernel (set at $\rho=1/2$ as a Bhattacharyya
kernel) for multinomial models as well as error rates for traditional
kernels on the WebKB data set. Performance for multiple settings of the
regularization parameter $C$ in the support vector machine are shown.
Two different sizes of training sets are shown (77 and 622).
All results are shown for 20-fold cross-validation.}
\label{fig:webkb}
\end{figure}

A support vector machine (which is known to have strong empirical
performance on such text classification data) was then used to build a
classifier. The SVM was fed our kernel's evaluations via a gram matrix
and was trained to discriminate between faculty and student web pages
in the WebKB data set (other categories were omitted). The probability
product kernel was computed for multinomials under the parameter
settings $\rho=1/2$ (Bhattacharyya kernel) and $s=1$ as in
Section~\ref{sec: exp fam}.  Comparisons were made with the (linear)
dot product kernel as well as Gaussian RBF kernels. The data set
contains a total of 1641 student web pages and 1124 faculty web
pages. The data for each class is further split into 4 universities
and 1 miscellaneous category and we performed the usual training and
testing split as described by \citep{LafLeb02,JoaCriSha01} where
testing is performed on a held out university. The average error was
computed from 20-fold cross-validation for the different kernels as a
function of the support vector machine regularization parameter $C$ in
Figure~\ref{fig:webkb}. The figure shows error rates for different
sizes of the training set (77 \COMMENT{, 155, 311} and 622 training
points). In addition, we show the standard deviation of the error rate
for the Bhattacharyya kernel.  Even though we only explored a single
setting of $s=1,\rho=1/2$ for the probability product kernel, it
outperforms the linear kernel as well as the RBF kernel at multiple
settings of the RBF $\sigma$ parameter (we attempted
$\sigma=\{1/4,1,4\}$). This result confirms previous intuitions in the
text classification community which suggest using squashing functions
on word frequencies such as the logarithm or the square root
\citep{GolSah98,CutKarPedTuk92}. This also related to the
power-transformation used in statistics known as the Box-Cox
transformation which may help make data seem more
Gaussian \citep{DavMac93}.

\subsection{Biological Sequence Classification}

For discrete sequences, natural generative models to consider are
Markov models and hidden Markov models. We obtained labeled gene
sequences from the $HS^3D$ data set\footnote{This is a collection of
unprocessed intron and exon class gene sequences referred to as the
{\em homo sapiens splice sites data set} and was downloaded from
www.sci.unisannio.it/docenti/rampone.} The sequences are of variable
lengths and have discrete symbol entries from the 4-letter alphabet
(G,A,T,C). We built classifiers to distinguish between gene sequences
of two different types: introns and exons using raw sequences of
varying lengths (from dozens of characters to tens of thousands of
characters). From the original (unwindowed) 4450 introns and 3752
exons extracted directly from GenBank, we selected a random subset of
500 introns and 500 exons. These 1000 sequences were then randomly
split into a 50\% training and a 50\% testing subset. In the first
experiment we used the training data to learn stationary hidden Markov models
whose number of states $M$ was related to the length $T_n$ of a given
sequence $n$ as follows: $M={\rm
floor}(\frac{1}{2}\sqrt{O^2+4(T\zeta+O+1)}-\frac{1}{2}O)+1$. Here, $O$
is the number of possible emission symbols (for gene sequences, $O=4$)
and $\zeta$ is the ratio of parameters to the number of symbols in the
sequence $T$ (we used $\zeta=1/10$ although higher values may help
avoid over-parameterizing the models). Each hidden Markov model was
built from each sequence using the standard Expectation-Maximization
algorithm (which is iterated for 400 steps to ensure convergence, this
is particularly crucial for longer sequences). Gram matrices of size
$1000 \times 1000$ were then formed by computing the probability
product kernel between all the hidden Markov models.  It was noted
that all Gram matrices were positive definite by a singular value
decomposition. We also used the following standard normalization
of the kernel (a typical pre-processing step in the literature):
\[ 
{\tilde k}(p,p') & \leftarrow & \frac{{\tilde k}(p,p')}{\sqrt{{\tilde
k}(p,p)} \sqrt{{\tilde k}(p',p')}}.
\]

\begin{figure}[htbp]
  \center
  \begin{tabular}[b]{cc}
\includegraphics[height=2.27in]{hmmintronexon.eps}  &
\includegraphics[height=2.27in]{kmersintronexon.eps} \\
(a) Hidden Markov Models (1st Order)  & (b) Markov Models (0th to 4th Order)
  \end{tabular}\caption[ ]{SVM error rates for the probability product kernel at $\rho=1$ for hidden Markov models and Markov models (equivalent to string or spectrum kernels). In (a) test error rate under various levels of regularization is shown for 5 different settings of $T$ in the probability product kernel. In (b) test error rate under various levels of regularization is shown for 5 different settings of the order $K$ of the Markov model.}
\label{fig:genes}
\end{figure}

Subsequently, we trained a support vector machine classifier on 500
example sequences and tested on the remaining
500. Figure~\ref{fig:genes}(a) depicts the resulting error rate as we
vary $C$, the regularization parameter on the SVM as well as $T$ the
number of time steps used by the kernel. Throughout, these
experiments, we only used the setting $\rho=1$.  Note that these
models were 1st order hidden Markov models where each state only
depends on the previous state. Varying $T$, the length of the hidden
Markov models used in the kernel computation has a slight effect on
performance and it seems $T=9$ at an appropriate $C$ value performed
best with an error rate as low as 10-11\%.

\begin{figure}
  \center
  \begin{tabular}[b]{c} 
\scalebox{0.16}{\rotatebox{270}{\includegraphics*[-0.2in,-0.2in][3.9in,11.7in]{mm.prn}}}
  \end{tabular}
\caption{A higher order (2nd order or $K=3$) Markov model.}
\label{fig:mm}
\end{figure}

For comparison, we also computed more standard string kernels such as
the k-mer counts or spectrum kernels on the same training and testing
gene sequences as shown in Figure~\ref{fig:genes}(b). These basically
correspond to a fully observed (non-hidden) stationary Markov model as
in Figure~\ref{fig:mm}. The zeroth order ($K=1$) Markov model does poorly
with its lowest error rate around 34\%. The first order Markov model ($K=2$)
performs better with its lowest error rate at 18\% yet is slightly worse
than first order hidden Markov models. This helps motivate the
possibility that dependence on the past in first order Markov models
could be better modeled by a latent state as opposed to just the
previous emission. This is despite the maximum likelihood estimation
issues and local minima that we must contend with when training hidden
Markov models using expectation-maximization. The second order Markov
models with $K=3$ fare best with an error rate of about 7-8\%. However,
higher orders ($K=4$ and $K=5$) seem to reduce
performance. Therefore, a natural next experiment would be to consider
higher order hidden Markov models, most notably second order ones
where the emission distributions are conditioned on the past two
states as opposed to the past two emitted symbols in the sequence.

\COMMENT{

load logQ2T1; load logQ2T2; load logQ2T3; load logQ2T4; load logQ2T5;
clf;
semilogx(logQ2T1(:,3),logQ2T1(:,5),'X-');
hold on;
semilogx(logQ2T2(:,3),logQ2T2(:,5),'X:');
semilogx(logQ2T3(:,3),logQ2T3(:,5),'X--');
semilogx(logQ2T4(:,3),logQ2T4(:,5),'O-');
semilogx(logQ2T5(:,3),logQ2T5(:,5),'O:');
h = gca;
set(h,'FontSize',15);
xlabel('Regularization C');
ylabel('Error Rate');
axis([3 500 0.05 0.4])
legend('T=3','T=4','T=5','T=6','T=7');



figure(7);
cd finalRho10.7.8.9.10.11
load logQ2T1; load logQ2T2; load logQ2T3; load logQ2T4; load logQ2T5;
clf;
plot(log(logQ2T1(:,3)),logQ2T1(:,5),'X-');
hold on;
plot(log(logQ2T2(:,3)),logQ2T2(:,5),'X:');
plot(log(logQ2T3(:,3)),logQ2T3(:,5),'X--');
plot(log(logQ2T4(:,3)),logQ2T4(:,5),'O-');
plot(log(logQ2T5(:,3)),logQ2T5(:,5),'O:');
h = gca;
set(h,'FontSize',15);
xlabel('Regularization log(C)');
ylabel('Error Rate');
axis([log(2.25) log(985) 0.05 0.4])
legend('T=7','T=8','T=9','T=10','T=11');


figure(8);
cd fullRho10.37
load logQ2T1; load logQ2T2; load logQ2T3; load logQ2T4; load logQ2T5;
clf;
plot(log(logQ2T1(:,3)),logQ2T1(:,5),'X-');
hold on;
plot(log(logQ2T2(:,3)),logQ2T2(:,5),'X:');
plot(log(logQ2T3(:,3)),logQ2T3(:,5),'X--');
plot(log(logQ2T4(:,3)),logQ2T4(:,5),'O-');
plot(log(logQ2T5(:,3)),logQ2T5(:,5),'O:');
h = gca;
set(h,'FontSize',15);
xlabel('Regularization log(C)');
ylabel('Error Rate');
axis([log(2.25) log(985) 0.05 0.4])
legend('T=3','T=4','T=5','T=6','T=7');

cd ..
figure(9);
cd kmers
load logQ2T1; load logQ2T2; load logQ2T3; load logQ2T4; load logQ2T5;
clf;
plot(log(logQ2T1(:,3)),logQ2T1(:,5),'X-');
hold on;
plot(log(logQ2T2(:,3)),logQ2T2(:,5),'X:');
plot(log(logQ2T3(:,3)),logQ2T3(:,5),'X--');
plot(log(logQ2T4(:,3)),logQ2T4(:,5),'O-');
plot(log(logQ2T5(:,3)),logQ2T5(:,5),'O:');
h = gca;
set(h,'FontSize',15);
xlabel('Regularization log(C)');
ylabel('Error Rate');
axis([log(2.25) log(985) 0.05 0.4])
legend('K=1','K=2','K=3','K=4','K=5');
cd ..

}

\subsection{Time Series Experiments}
In this experiment we compare the performance of the probability
product kernel between LDSs and the sampling approximation of the SLDS
kernel with more traditional kernels and maximum likelihood techniques
on synthetic data. We sampled from 20 exemplar distributions to
generate synthetic time series data. The sampled distributions were 5
state, 2 dimensional SLDS models generated at random. Each model was
sampled 5 times for 25 time steps with 10 models being assigned to
each class. This gave 50 time series of length 25 per class, creating
a challenging classification problem. For comparison, maximum
likelihood models for LDSs and SLDSs were fit to the data and used for
classification as well as SVMs trained using Gaussian RBF kernels. All
testing was performed using leave one out cross validation. The SLDS
kernel was approximated using sampling (50 samples) and the
probability product kernels were normalized to improve performance.

The maximum likelihood LDS had an error rate of .35, the SLDS .30, and
the Gaussian RBF .34. Exploring the setting of the parameters
$\rho$ and $T$ for the LDS and SLDS kernels, we were able to
outperform the maximum likelihood methods and the Gaussian RBF kernel
with an optimal error for the LDS and SLDS kernel of .28 and .25
respectively. Figure~\ref{fig:ldsgraph} (a) and (b) show the error
rate versus $T$ and $\rho$ respectively. It can be seen that
increasing $T$ generally improves performance which is to be
expected. Although it does appear that $T$ can be chosen too
large. Meanwhile, $\rho$, generally appears to perform better when it
is smaller (a counter example is the LDS kernel at the setting
$T=5$). Overall, the kernels performed comparably to or better than
standard methods for most settings of $\rho$ and $T$ with the exception of
extreme settings.
\begin{figure}
  \center
\setlength\tabcolsep{0pt}
  \begin{tabular}[b]{cc} 
\scalebox{0.427}{\includegraphics*[0.7in,2.5in][8in,9in]{ldsgraph1.eps}} &
\scalebox{0.427}{\includegraphics*[0.7in,2.5in][8in,9in]{ldsgraph2.eps}} \\
(a)  & (b) 
  \end{tabular}
\caption{A comparison of the choice of parameters (a) $T$ and (b) $\rho$. The $x$-axis is the parameter and the $y$-axis is the error rate.}
\label{fig:ldsgraph}
\end{figure}

\section{Conclusions}

Discriminative learning, statistical learning theory and kernel-based
algorithms have brought mathematical rigor and practical success to
the field making it is easy to overlook the advantages of generative
probabilistic modeling.  Yet generative models are intuitive and offer
flexibility for inserting structure, handling unusual input spaces,
and accommodating prior knowledge.  By proposing a kernel between the
models themselves, this paper provides a bridge between the two
schools of thought.

The form of the kernel is almost as simple as possible, yet it still
gives rise to interesting nonlinearities and properties when applied
to different generative models. It can be computed explicitly for some
of the most commonly used parametric distributions in statistics such
as the exponential family. For more elaborate models (graphical
models, hidden Markov models, linear dynamical systems, etc.),
computing the kernel reduces to using standard algorithms in the
field. Experiments show that probability product kernels hold promise
in practical domains. Ultimately, engineering kernels between
structured, discrete or unusual objects is an area of active research
and, via generative modeling, probability product kernels can bring
additional flexibility, formalism and potential.



\acks{The authors thank the anonymous reviewers for their helpful
comments. This work was funded in part by the National Science
Foundation (grants IIS-0347499 and CCR-0312690).}


\bibliography{everything}

\end{document}
