\documentstyle[fleqn,psfig]{weigend}

\newcommand{\ggg}{ }
\newcommand{\lll}{ }

\include{/nfs/home/dayan/pl2/co/include/inct}


\begin{document}

\bibliographystyle{weigbib}

\chapter{Curved Gaussian Models}

\author{Juan K Lin \& Peter Dayan}

\abstract
  In principal component analysis, data are modeled with arbitrary
  linear transformations of a spherical Gaussian.  In this paper, we
  present a simple extension to generalised linear transformations,
  which allow for an efficient local description of curvature. In a
  predetermined coordinate system, direct analytic computation of the
  curvature parameters is possible; otherwise, the coordinate system
  can also be inferred.
\endabstract

\section{Introduction}

Gaussian models for data, as found in standard density models,
principal components analysis, factor analysis, are key {\it building
  blocks\/} for more complicated density estimation models. They have
numerous advantages, notably including the simplicity with which they
can be fit to data. There is a plethora of more complicated density
estimation edifices built using these blocks, including simple
mixtures of either radial Gaussians or Gaussians with arbitary
covariance matrices ({\it eg\/} Nowlan, 1991); and, in more recent
work, mixtures of principal component or factor analysers (Bregler \&
Omohundro, 1995; Kambhatala \& Leen, 1994; Hinton, Dayan \& Revow,
1997; Roweis, 1998; Tipping \& Bishop, 1997). 

Figure~\ref{fig:fit} illustrates one rationale underlying the use of
mixture models. Figure~\ref{fig:fit}A shows $1769$ two dimensional
data points from a particular financial time series. One reasonable
way to characterise these points is as lying near to a one-dimensional
curved manifold. In this representation, two dimensional Gaussian
models are ellipses (assuming that the full covariance matrices are
modeled). Figures~\ref{fig:fit}B;C show two alternative fits of four
two-dimensional Gaussians to the data (using standard expectation
maximisation algorithm, Dempster, Laird \& Rubin, 1977, to fit
mixtures of Gaussian). The ellipses show $\sqrt{2}\sigma$ contours for
the Gaussians; the widths of the outlines are proportional to the
prior probabilities (mixing proportions) of the Gaussians. These
Gaussian construct the manifold by `bolting' pieces together -- each
piece is strongly elliptical, matching the one-dimensional structure
of the data.

\begin{figure}
\setlength{\unitlength}{.1in}
\begin{picture}(49, 15)(0, 0)
\put(0,0) {\psfig{figure=figs/basedata.eps,width=1.5in}}
\put(17,0) {\psfig{figure=figs/mog1.eps,width=1.5in}}
\put(34,0) {\psfig{figure=figs/mog2.eps,width=1.5in}}
\put(0,13) {\Large A}
\put(17,13) {\Large B}
\put(34,13) {\Large C}
\end{picture}
\caption{Mixture of Gaussians models for curved data. A) Daily values
  of the Japanese Yen {\it versus\/} the Korean Won measured in US
  Dollar units from 23rd October 1993 to 27th August 1998 normalised
  relative to their values at the end. B;C) Two different mixtures of
  Gaussian fits to these data using four full covariance
  two-dimensional Gaussians. The centers of the Gaussians are marked
  by a cross, $\sqrt{2}\sigma$ values by the ellipses, and the mixing
  proportions by the widths of the outlines of the ellipses.
\label{fig:fit}
}
\end{figure}

These Gaussian models capture part of the essence of the data, and
indeed, had we used more components, would have fit them even better.
However, they miss the simple observation that the data seem to lie on
a curved manifold. In this paper, we consider ways of changing the
Gaussian building block so that it can capture curvatures of this
form, whilst retaining as much of its computational attractiveness as
possible. 

Section~\ref{sec:bas} discusses the basic curved Gaussian model;
section~\ref{sec:cord} considers the problem of inferring appropriate
coordinate systems in which to describe the curvature. Relationships
with methods such as principal curves (Hastie \& Tibshirani, 19??) are
described in section~\ref{sec:conc}.

\section{THE BASIC BUILDLING BLOCK\label{sec:bas}}

\subsection{Two-dimensions}

One way to write a Gaussian distribution $\Nn[0,\Sigma]$ over two
variables $\xx=(x_1,x_2)$ is as:
\[
p(x_1,x_2) \propto \sqrt{|\Omega|} e^{-\xx^T \Omega
  \xx/2}
= \sqrt{w_1}e^{-w_1\phi_1^2/2}
 \sqrt{w_2}e^{-w_2\varphi_2^2/2}, \qquad \mbox{where}
\]
$\Omega=\Sigma^{-1}$, and 
\begin{eqnarray}
\varphi_1(\xx) & = & x_1 \ ,\\
\varphi_2(\xx) & = & x_2 + k_{21} x_1  \ ,\label{eq:key}\\
w_1 & = & \Omega_{11} - \Omega_{12}^2/\Omega_{22} \qquad , \qquad
w_2 \  = \  \Omega_{22} \qquad , \qquad
k_{21} \  = \  \Omega_{12}/\Omega_{22} \ .
\end{eqnarray}
The key facet of the change of coordinates from $x_1,x_2$ to
$\varphi_1,\varphi_2$ is the linear dependence of the mean of
$\varphi_2$ on $x_1$ (and hence $\varphi_1$). Figure~\ref{fig:gau}A
shows an example of the contours of such a distribution (the
ellipses), and the underlying linear dependence (the dashed line).
Note that the principal axis of the covariance matrix is not
necessarily aligned with this linear relation.

\begin{figure}
\setlength{\unitlength}{.1in}
\begin{picture}(49, 19)(0, 0)
\put(0,0) {\psfig{figure=figs/gasel.eps,width=2in}}
\put(25,0) {\psfig{figure=figs/smil.eps,width=2in}}
\put(0,17) {\Large A}
\put(25,17) {\Large B}
\end{picture}
\caption{$x_1,x_2$ plots of conventional and curved Gaussians. A) Conventional
  two-dimensional Gaussian (in $x_1,x_2$ space) with covariance matrix
  $\Sigma=\left( \begin{array}{cc} 1 & .8 \\ .8 & 2
\end{array}
\right)$. The ellipses
  are $\sigma, \sqrt{2}\sigma$ and $\sqrt{5}\sigma$ contours; the
  dashed line is $x_2 = 0.8x_1$ as in equation~\ref{eq:key};
  $3000$ points were drawn from the distribution. B) Curved
  two-dimensional Gaussian with $w_1=1$, $w_2=3$ and $k_{21}=1.5$. The ellipses
  are $\sigma, \sqrt{2}\sigma$ and $\sqrt{5}\sigma$ contours; the
  dashed line is the skeleton $x_2 = k_{21}x_1^2$; $3000$ points were
  drawn from the distribution. 
\label{fig:gau}
}
\end{figure}

This linear change in coordinates also suggests a natural
generalisation:
\begin{equation}
\varphi_2 = x_2 - \sum_{j} k_{2j} \psi_{2j}(x_1)
\label{eq:ckey}
\end{equation}
where $\psi_{2j}(x)$ are any suitably well-behaved (everywhere finite and
differentiable) and {\it fixed\/} basis functions of $x_1$, and
$k_{2j}$ are parameters. This is a {\it generalised linear model\/}
(Hastie \& Tibshirani, 19??) for $x_2$, since the mean depends
linearly on the underlying variable parameters ($k_{2j}$). 

One of the simplest possible generalised linear models is the quadratic:
\begin{equation}
\varphi_2 = x_2 k_{21} x_1^2
\label{eq:cckey}
\end{equation}
Consider therefore the full Gaussian model:
\[
p(\xx;\ww;k_{21}) \propto 
(w_1w_2)^{1/2} e^{-w_1 \varphi_1^2 - w_2 \varphi_2^2}
\]
Since the Jacobian of the transformation from $x_1,x_2$ to
$\varphi_1,\varphi_2$ using equation~\ref{eq:cckey} in place of
equation~\ref{eq:key} is lower diagonal everywhere, with the diagonal
elements consisting of ones, the determinant of the Jacobian is always
one. Therefore, this is the probability density function in the new
curved coordinate system.  More generally, the multivariate
statistical models being considered are subsets of volume preserving
non--linear transformations of multivariate Gaussian models.
Figure~\ref{fig:gau}B shows a contour plot of this Gaussian---the
`scaffolding' for $x_2$ provided by the value of $x_1$ (the dashed
line) is now quadratic rather than linear.

\subsection{$n$ dimensions\label{sec:nd}}

More generally, in the case that the input data
$\xx=(x_1,x_2,\ldots,x_n)$ are $n$\-dimensional, consider the following
non-linear change of variables:
\begin{eqnarray*}
\varphi_1 &=& x_1 + k_{10} \\
\varphi_2 &=& x_2 + k_{20} + \sum_j k_{2j} \psi_{2j}(x_1)  \\
\varphi_3 &=& x_3 + k_{30} + \sum_j k_{3j} \psi_{3j}(x_1,x_2)\\
\vdots \\
\varphi_n &=& x_n + k_{n0} + \sum_j k_{3j}
\psi_{nj}(x_1,x_2,\ldots,x_{n-1}) \ ,
\end{eqnarray*}
together with a likelihood model
\begin{equation}
p(\mathbf{x};\mathbf{w},k_1,\mathbf{k_2},...,\mathbf{k_n}) \propto
(w_1 ... w_n)^{1/2} \exp 
[-w_1 \varphi_1^2 ... -w_{n-1} \varphi_{n-1}^2 -w_n \varphi^2],
\label{eq:like}
\end{equation}
where $\phi_{ij}$ are all differentiable functions of their arguments
and $k_{ij}$ are parameters. Again, the determinant of the Jacobian of
the transformation is $1$, and so the model is well defined. Natural
choices for $\phi_{ij}$ are polynomials of non-negative orders in
their parameters, although other choices are possible.

\subsection{Parameter estimation}

Consider fitting a model of the form of section~\ref{sec:nd} to data
$\Dd=\left\{\xx^m\right\}$. Maximising the log likelihood with respect to
the model parameters gives optimising equations
\begin{equation}
k_{ij} = \left[\left(\Phi^i\right)^{-1} \cc^i\right]_j
\label{eq:invp}
\end{equation}
where $\Phi^i_{ab}=\langle \phi_{ia}\phi_{ib}\rangle$, $c^i_a =
\langle \phi_{ia} x_i\rangle$, where $\langle \cdot \rangle$ indicates
taking an average over $\Dd$, and, for convenience, we define
$\phi_{l0}=1$ for all $l$. Equation~\ref{eq:invp} is only valid if
matrix $\Phi^i$ is invertible; otherwise, the maximum likelihood
solution is degenerate with whole manifolds of solutions.  Equating to
zero the partial derivatives of the log likelihood with respect to the
inverse variance parameters, we find
\begin{equation}
\frac{\partial L}{\partial w_j}=
\frac{1}{2 w_j} - \langle \varphi_j^2 \rangle.
\end{equation}
Thus the maximum likelihood solution is 
\begin{equation}
w_j = \frac{1}{2 \langle \varphi_j^2 \rangle}.
\label{eq:wfit}
\end{equation}
using the optimal values of $k_{ij}$ given in equation~\ref{eq:invp}. 

\begin{figure}
\setlength{\unitlength}{.1in}
\begin{picture}(49, 19)(0, 0)
\put(0,0) {\psfig{figure=figs/smifit.eps,width=2in}}
\put(25,0) {\psfig{figure=figs/smimfit.eps,width=2in}}
\put(0,17) {\Large A}
\put(25,17) {\Large B}
\end{picture}
\caption{Curved Gaussian fits using
  equations~\ref{eq:invp}~and~\ref{eq:wfit} to the data in
  figure~\ref{fig:gau}B with $\phi_{21}=x_1^2$. The dashed line shows
  the quadratic skeleton; the solid line the $\sqrt{1.5}\sigma$
  contour. A) Fit using the correct ordering of the
  variables $x_1,x_2$. B) Fit using the opposite ordering $x_2, x_1$.
\label{fig:cvfit}
}
\end{figure}

Figure~\ref{fig:cvfit}A shows the maximum likelihood fit of a curved
Gaussian to the data in figure~\ref{fig:gau}B. The curved model
captures the structure of the data points well, as should be expected,
given the way that they were generated.

We have therefore suggested a natural curved generalisation of a
conventional Gaussian model. Given the coordinate systems ({\it ie\/}
the {\it ordering\/}) of the variables, computing the maximum
likelihood values of the parameters of the generalised linear model
requires nothing more than matrix inversion---in particular, it is
{\it not\/} an iterative process. One could use even more general
parameterised versions of the functions $\phi_{ij}$, and determine the
values of the parameters by gradient ascent in the log likelihood.


\section{COORDINATE TRANSFORMATIONS\label{sec:cord}}

Figure~\ref{fig:cvfit}B shows an extreme example of the consequence of
using a different coordinate system to fit the curved model -- in this
case, one that has the direction of the non-linear dependence exactly
backwards. In this case, the model clearly mis-fits the data -- doing
just as badly, or perhaps even worse, as a non-curved Gaussian
model. For fitting general datasets, such as that in
figure~\ref{fig:fit}, it is therefore necessary to be able to change
coordinate systems. 

Drawing inspiration from Independent Component Analysis, one way to do
this is to consider a curved, source datapoint $\xx$ being used to
generate an observed datapoint $\yy$ through a linear transformation:
\begin{equation}
\mathbf{y} = \Wws^{-1} \mathbf{x} + \bb
\end{equation}
where $\mathbf{b}$ is an offset, and $\Wws^{-1}$ is an invertible
`mixing' transformation. Given a dataset $\Dd=\left\{\yy^m\right\}$,
the task is to find the unmixing transformations $\bb$ and $\Wws$,
together with the parameters $\ww$ and $k_{ij}$ of the curved source
model so as to fit the data as tightly as possible, by maximising the
log likelihood, the contribution to which from $\yy^m$ is:
\[
\log |\Wws^{-1}| + \log
p(\Wws^{-1}(\yy^m-\bb),\kk_1,\mathbf{k_2},...,\mathbf{k_n})
\]
using the expression in equation~\ref{eq:like}.  It does not seem
possible, in general, to derive closed-form optima for the likelihood
as a function of all the paramters. Instead, we consider a form of
coordinate-wise ascent of the likelihood, finding the global optimum
in $\ww$ and $k_{ij}$ for values of $\Wws$ and $\bb$, and then
changing the latter variables by gradient ascent.

The derivatives of the log likelihood with respect to the mixing
parameters are straightforward to calculate, albeit rather messy. Note
that there is substantial redundancy in the model -- $\Wws$ can be
restricted to only $n(n-1)/2$ parameters corresponding to an arbitrary
orthogonal matrix with unit determinant and there are also
interdependencies between $\bb$ and $k_{i0}$ for $i$.  However, the
learning rule is considerably simpler with the model parameterization
given above. Just as for ICA, we can use the natural gradient ascent
algorithm (Amari \& Cichocki \& Yang 1995, Cardoso \& Laheld 1996) to
update the mixing transformation along the direction of steepest
ascent by right multiplying the partial derivatives of the log
likelihood with respect to $\Wws$ by $\Wws^T \Wws$. One way to choose
initial conditions for the mixing transformation is to apply standard
principal components analysis -- although this can be problematical
depending on the scaling of the variables.

\begin{figure}
\setlength{\unitlength}{.1in}
\begin{picture}(49, 15)(0, 0)
\put(0,0) {\psfig{figure=figs/fulfit.eps,width=1.5in}}
\put(17,0) {\psfig{figure=figs/fulfit2.eps,width=1.5in}}
\put(34,0) {\psfig{figure=figs/fulfail.eps,width=1.5in}}
\put(0,13) {\Large A}
\put(17,13) {\Large B}
\put(34,13) {\Large C}
\end{picture}
\caption{Inference of the mixing matrix. All graphs show the starting
  and ending $\sqrt{1.5}\sigma$ contours, the ending skeleton (dashed
  line), and intermediate contours every $10$ iterations. A;B)
  Examples of inference to the global optimal coordinate
  transformation. C) Example of inference to the poor local optimum.
\label{fig:fulfit}
}
\end{figure}

Figures~\ref{fig:fulfit}A-C show examples of the way the algorithm fits
the curved Gaussian data. In each case, the initial values of the
mixing transformation are arbitrary, and the shape of the
$\sqrt{1.5}\sigma$ contour is plotted every $10$ iterations (using a
rather arbitrary step-size). Figures~\ref{fig:fulfit}A;B show examples
of successful inference of the transformation and the curved Gaussian;
figure~\ref{fig:fulfit}C shows an example of failure, in which only the local
maximum of the likelihood (shown in figure~\ref{fig:cvfit}B) is
attained. The absolute values of the likelihood can be used to select
between alternative solutions.

Finally, we can apply the full algorithm to the currency data from
figure~\ref{fig:fit}. Figure~\ref{fig:fulcur} shows an example of a
full curved Gaussian model for the data. It clearly captures the
structure of the data more competently than the mixture of non-curved
Gaussians shown in figures~\ref{fig:fit}B;C.

\begin{figure}
\centerline{
\psfig{figure=figs/fulcur.eps,width=2in}
}
\caption{Full curved model fit to the currency data from
  figure~\ref{fig:fit}A. The dashed line is the quadratic skeleton;
  the solid line the $\sqrt{2}\sigma$ contour for the model.
\label{fig:fulcur}
}
\end{figure}

\section{DISCUSSION\label{sec:conc}}

We have presented a simple extension to a Gaussian density model that
captures weakly curved manifolds. We used a form of generalised
linear model for the mean of one latent variable given the values of
the other latent variables that precede it in some particular ordering
of the latent variables. The form of the non-linear variable
transformation chosen is convenient since the Jacobian is guaranteed
to be $1$. The advantage of such generalised linear models is that
there is a simple procedure for fitting them, although more general
models also possible. We also showed how to use ideas from independent
components analysis to infer the appropriate coordinate system in
which to fit the curved model, although this leads to an algorithm
that has to be fit iteratively.

There are various natural extensions of the work. Particularly
important is the idea of using {\it mixtures\/} of the curved Gaussian
distributions, in the same way that one uses mixtures of standard
Gaussian distributions. This presents no conceptual or computational
hurdle -- and can be done simply using the expectation-maximisation
algorithm. 

Another extension is that the local curvature information could be
used for things other than fitting a Gaussian model -- for instance,
it could be used to enhance local kernel methods in something like a
curved form of tangent distance (Simard, Le~Cun \& Denker, 1993).
Also, more general curved models could be used in place of the simple
quadratic form we used to fit the data of figure~\ref{fig:fit}A. As
ever, there are trade-offs of bias for variance in assuming more
flexible parameterisations, and more sophisticated cross-validation
or Bayesian methods could be considered to choose the orders of the
polynomials or the forms of the curvature.

A further possibility is to consider reduced-order, curved, Gaussian
models of data (which is closely related to the idea of performing ICA
when there are fewer sources than inputs, or in the presence of
assumed noise in the transformation from the sources to the inputs,
Attias, 19??). Inference of the coordinate system is more difficult in
this case, since the equivalent of matrix $\Wws$ is singular, but
standard methods can be applied.

\end{document}

%for higher order polynomial fit with $m$ parameters, just need to
%invert an mxm matrix, still analytic solution for global maximum.
%DECOMPOSITION OF ARBITRARY MATRIX INTO COMPOSITION OF ROTATION AND SHEER
%SVD DECOMPOSITION

%\bibliographystyle{refsty}
%\bibliography{juan}
%\bibstyle{refsty}
%\bibdata{juan}

\begin{thebibliography}{99}
%\fontsize{8}{8pt}\selectfont
\parskip 0pt
\bibitem{amari}
Amari, S \& Cichocki, A \& Yang, H (1996).
``A new learning algorithm for blind signal separation'',
in {\em Advances in Neural Information Processing Systems 8},
757--763.
%%%%%%%%%%%%%%%%%  add in editors

\bibitem{bo}
Bregler, C \& Omohundro, SM (1995).
``Nonlinear image interpolation using manifold learning'',
in {\em Advances in Neural Information Processing Systems 7},
  G.~Tesauro, D.~S. Touretzky, and T.~K. Leen, Eds., pp. 971--980. MIT Press,

\bibitem{cardoso}
Cardoso, J-F \& Laheld, B (1996).
``Equivariant adaptive source separation'',
in {\it IEEE Transactions on Signal Processing}, {\bf 45:2} 434--444.

\bibitem{dlr}
Dempster, AP, Laird, NM \& Rubin, DB (1977) Maximum likelihood from
incomplete data via the {EM} algorithm.\ {\it Journal of the Royal Statistical Society, Series B} {\bf 39}:1-38.

\bibitem{hdfn}
Hinton, GE, Dayan, P, Frey, BJ \& Neal, RM (1995). The wake-sleep
algorithm for unsupervised neural networks. {\it Science,\/} {\bf
268}, 1158-1160.

\bibitem{hdr}
Hinton, GE, Dayan, P \& Revow, M (1997). Modeling the manifolds of
images of handwritten digits. {\it IEEE Transactions on Neural
Networks,\/} {\bf 8}, 65-74.

%\bibitem{kl}
%Kambhatla, N and Leen, TK (1994).
%Fast non-linear dimension reduction.
%In JD Cowan, G Tesauro \& J Alspector, editors, {\em Advances in Neural
%  Information Processing Systems 6,\/}, San Mateo, CA: Morgan
%Kaufmann, 152-159.

\bibitem{karh}
Karhunen, J and Joutsensalk, J (1995).
Generalizations of principle component analysis, optimization
problems and neural networks'', {\it Neural Networks} vol 8,
no 4, 549-562

\bibitem{kl2}
Kambhatla, N and Leen, TK (1997).
Dimension reduction by local principal component analysis.
{\it Neural Computation} {\bf 9} 1493--1516.

\bibitem{n}
Nowlan, SJ (1991) {\it Soft Competitive Adaptation: Neural Network
Learning Algorithms Based on Fitting Statistical Mixtures}.\ PhD
Thesis, Deparment of Computer Science, Carnegie-Mellon University.

\bibitem{oja}
Oja, E. (1995)
The nonlinear PCA learning rule and signal separation - mathematical
analysis. Helsinki University of Technology, Laboratory of Computer and
Information Science, Report A26. 

\bibitem{r}
Roweis, S (1998). EM algorithms for PCA and SPCA. In M Jordan, M
Kearns \& S Solla, editors, {\it Advances in Neural Information
  Processing, 10.\/} Cambridge, MA: MIT Press.

\bibitem{sld}
Simard, P, LeCun, Y \& Denker, J (1993).
Efficient pattern recognition using a new transformation
  distance.
In JD Cowan, SJ hanson \& CL Giles, editors, {\em Advances in Neural
  Information Processing Systems 5,\/}, San Mateo, CA: Morgan
Kaufmann, 50-58.

\bibitem{tb}
Tipping, ME \& Bishop, CM (1997). 
     Mixtures of principal component analyzers.
In {\it Fifth International Conference on Artificial Neural
     Networks.\/} London, UK: IEE.

\end{thebibliography}

\end{document}