\documentclass[twoside,11pt]{article}
%\usepackage{pierre}
\usepackage[fleqn]{amstex}
\usepackage{epsfig}
%\usepackage{float}
\usepackage{subfigure}
%\usepackage{hhline}
%\usepackage{rotate}
\usepackage{xspace}
%\usepackage{fancyheadings}
%\usepackage{html}
\usepackage{jair}
\usepackage{theapa}

\newcommand{\home}{/home/snn/davidb}
\newcommand{\eg}{e.g.,\xspace}
\newcommand{\ie}{i.e.,\xspace}
\newcommand{\eqn}[1]{equation (#1)}
\newcommand{\Eqn}[1]{Equation (#1)}	%At beginning of sentence
\newcommand{\eqns}[2]{equations (#1) and~(#2)}
\newcommand{\fig}[1]{Figure #1}
\newcommand{\figs}[2]{Figures #1 and #2}
\newcommand{\subfig}[2]{Figure #1(#2)}
\newcommand{\chap}[1]{chapter #1}		%
\newcommand{\Chap}[1]{Chapter #1}		%At beginning of sentence
\newcommand{\sect}[1]{section (#1)}
\newcommand{\Sect}[1]{Section (#1)}		%At beginning of sentence
\newcommand{\subsect}[1]{subsection (#1)}
\newcommand{\Subsect}[1]{Subsection (#1)}	%At beginning of sentence
\newcommand{\tab}[1]{table #1}
\newcommand{\Tab}[1]{Table #1}			%At beginning of sentence
%\newcommand{\cut}[1]{}			%At beginning of sentence
%\renewcommand{\secref}[1]{section #1}			%At beginning of sentence
\newcommand{\partialdd}[2]{\frac{\partial #1}{\partial #2}}

\newcommand{\qzero}{Q_0}
\newcommand{\qone}{Q_1}
\newcommand{\qtwo}{{\tilde{Q}}_0}
\newcommand{\sm}[1]{{s}_{ \backslash #1}}
\newcommand{\INP}{X}
\newcommand{\hid}{h}
\newcommand{\HID}{H}
\newcommand{\out}{o}
\newcommand{\OUT}{O}
\newcommand{\TARGET}{T}
\newcommand{\learningparameter}{\ensuremath{\eta}}

\newcommand{\nrofinput}{N}
\newcommand{\nrofoutput}{K}

\newcommand{\order}{{\cal{O}}}

\newcommand{\weight}{\ensuremath{w}}
%\newcommand{\probability}{{\cal Q}}
\newcommand{\probability}{{Q}}
\newcommand{\normalconst}{Z}
\newcommand{\entropy}{{\cal S}}
\newcommand{\meanfield}{m}

\newcommand{\spin}{s}

% NIET THESIS DEFINITIONS
\newcommand{\KL}{\text{KL}}
%\newcommand{\myvector}[1]{\ensuremath{\mathbf #1}}
\newcommand{\myvector}[1]{\mbox{\boldmath $#1$}}

\input{myutils}
\renewcommand{\var}{\text{var}}		%variance

\newcommand{\newweight}{\theta}		%NEW WEIGHTS (TRACTABLE)
\newcommand{\decweight}{\theta'}	%NEW WEIGHTS (DECIMATION&TRACTABL)
\newcommand{\hu}{S_{H}}
\newcommand{\vu}{S_{V}}
\newcommand{\teacher}[1]{p(#1)}
%\newcommand{\pdiff}[2]{\frac{\partial{#1}}{\partial{#2}}}
\newcommand{\pdiffo}[1]{\frac{\partial}{\partial{#1}}}
\newcommand{\J}{J}
\newcommand{\s}{s}

\newcommand{\bfw}{\myvector{\weight}}
%\newcommand{\bfgamma}{\mbox{\boldmath \gamma}}

%
% hyphenations (how to break words) in alfabetic order
%
\hyphenation{Hum-phreys}
\hyphenation{knowl-edge}
\hyphenation{lat-er-al}
\hyphenation{Nijs-kens}
\hyphenation{par-a-digm}


\jairheading{10}{1999}{435-455}{8/98}{6/99}
\ShortHeadings{Variational Cumulant Expansions for Intractable Distributions}
{Barber \& van de Laar}
\firstpageno{435}


\begin{document}
\title{Variational Cumulant Expansions\\ for Intractable Distributions}

\author{\name David Barber \email davidb@@mbfys.kun.nl\\
		\name Pi\"erre van de Laar \email pierre@@mbfys.kun.nl\\
        \addr RWCP (Real World Computing Partnership) \\
		Theoretical Foundation
        SNN (Foundation for Neural Networks)
		\\
%        Department of Medical Physics and Biophysics\\
        University of Nijmegen, 
        CPK1 231, Geert Grooteplein 21\\
        Nijmegen, The Netherlands\\
%        telephone: +31 24 3614235\\
%        fax: +31 24 3541435\\
    }

\maketitle


\begin{abstract}
Intractable distributions present a common difficulty in inference
within the probabilistic knowledge representation framework and
variational methods have recently been popular in providing an
approximate solution. In this article, we describe a perturbational
approach in the form of a cumulant expansion which, to lowest order,
recovers the standard Kullback-Leibler variational bound.
Higher-order terms describe corrections on the variational approach
without incurring much further computational cost.  The relationship
to other perturbational approaches such as TAP is also elucidated.  We
demonstrate the method on a particular class of undirected graphical
models, Boltzmann machines, for which our simulation results confirm
improved accuracy and enhanced stability during learning.
\end{abstract}

\section{Introduction}

In recent years, interest has steadily grown over many associated
fields in representing and processing information in a probabilistic
framework.
%
%This has been particularly the case in the artificial
%intelligence community, where it is generally accepted that
%incorporating prior beliefs of the problem domain can greatly
%facilitate the accuracy and plausibility of
%prediction~\cite{Bishop95,MacKay95,Jordan98b}.
%
%
%
%
Arguably, the reason for this is that prior knowledge of a problem
domain is rarely absolute and the probabilistic framework therefore
arises naturally as a candidate for dealing with this
uncertainty. Better known examples of models are Belief Networks
\cite{Pearl:book88} and probabilistic neural
networks \cite{Bishop95,MacKay95}.  However, incorporating
many prior beliefs can make models of domains so ambitious that
dealing with them in a probabilistic framework is intractable, and
some form of approximation is inevitable.  One well-known approximate
%technique is Monte Carlo sampling \citeA<see, for example>{Neal93,
%gilks-book-96}], which has the benefit of universal applicability.
technique is Monte Carlo sampling \cite<see \eg>{Neal93,gilks-book-96}, which has the benefit of universal
applicability.  However, not only can sampling be prohibitively slow,
but also the lack of a suitable convergence criterion can lead to low
confidence in the results.  Recently, variational methods have
provided a popular alternative, since they not only approximate but
can also bound quantities of interest giving, potentially, greater
confidence in the results~\cite{jordan-in-book-98,barber-wiegerinck-98b,Barber98b}.
The accuracy of variational methods, however, is limited by the
flexibility of the variational distribution employed.  Whilst, in
principle, the flexibility, and therefore, the accuracy of the model,
can be increased at will
%(see \eg \citeA{Jaakkola98}, \citeA{Lawrence98})
\cite<see \eg>{Jaakkola98,Lawrence98}, this generally is at the expense of incurring
even greater computational difficulties and optimization problems.

In this article, we describe an alternative approach which is a
perturbation around standard variational methods. This has the
property of potentially improving the performance at only a small
extra computational cost.  Since most quantities of interest are
derivable from knowledge of the normalizing constant of a
distribution, we present an approximation of this quantity which, to
lowest order, recovers the standard Kullback-Leibler variational
bound.  Higher-order corrections in this expansion are expected to
improve on the accuracy of the variational solution, despite the loss
of a strict bound.

In \sect{\ref{sec: Probability distributions}} we will briefly discuss
why the normalizing constant of probability distributions is of such
importance.  We briefly review the Kullback-Leibler variational bound
in \sect{\ref{sec: bound}} before introducing the perturbational
approach in \sect{\ref{sec: Approximation}}.  This approach is
illustrated on a well-known class of undirected graphical models,
Boltzmann machines, in \sect{\ref{sec: BM}}, in which we also explain
the relation of this approach to other perturbational methods such as
TAP~\cite{Thouless77,Kappen98}.  We conclude in \sect{\ref{sec:
Discussion}} with a discussion of the advantages and drawbacks of this
approach compared to better known techniques.

\section{Normalizing Constants and Generating Functions}
\label{sec: Probability distributions}

We consider a family of probability distributions $\probability$ over
a vector of random variables\\ $\myvector{\spin} = \left \{\spin_{i} |
i \in [1,\ldots,N] \right \}$, parameterized by
$\myvector{\weight}$,\footnote{We assume that the distribution is
strictly positive over the domain.}
%
\begin{equation}
    \probability(\myvector{\spin};\myvector{\weight}) = \frac{ \exp(H (\myvector{\spin};\myvector{\weight})) }{ \normalconst (\myvector{\weight})}
    \thickspace \text{.}
\end{equation}
With a slight abuse of notation, we write 
the normalizing constant as
\begin{equation}
	\label{eq: normalizing constant}
    \normalconst(\myvector{\weight}) = \int \negthinspace {d \myvector{\spin}} \thinspace \exp(H (\myvector{\spin};\myvector{\weight}))
    \thickspace \text{.}
\end{equation}
The `potential' $H (\myvector{\spin};\myvector{\weight})$ thus
uniquely determines the distribution $\probability$.  We have assumed
here that the random variable $\myvector{\spin}$ is continuous---if
%not, the corresponding integral [see \eqn{\ref{eq: normalizing
not, the corresponding integral (see equation \ref{eq: normalizing
constant}) should be replaced by a summation over all possible
discrete states.


One approach in statistics to obtain marginals and other quantities of
interest is given by considering generating functions
\cite{Grimmett92}, which take the form
%
\begin{equation}
G(\boldlambda;\bfw) = \int d \myvector{\spin}   \exp\br{\log Q(\myvector{\spin};\bfw) + \ip{\boldlambda}{\myvector{\spin}}}
\end{equation}
%
Averages of variables with respect to $Q$ are given by derivatives of
$G(\boldlambda;\bfw)$ with respect to $\boldlambda$ for fixed $\bfw$,
so that, for example, $\av{\spin_{1} \spin_{2}} = {\partial^2
G(\boldlambda;\bfw)}/{\partial \lambda_1 \partial\lambda_2}$,
evaluated at $\boldlambda=0$. We can equally well consider the
function
%
\begin{equation}
Z(\boldlambda;\bfw) = \int d \myvector{\spin}   \exp\br{ H(\myvector{\spin};\bfw) + \ip{\boldlambda}{\myvector{\spin}}}
\label{eq:modZ}
\end{equation}
%
so that
%
\begin{equation}
	\left < \spin_{1} \spin_{2} \right > =  {\left . \frac{1}{\normalconst(\lambda_1,\lambda_2;\bfw)} \frac{\partial^2 \normalconst(\lambda_1,\lambda_2;\bfw)}{\partial \lambda_1 \partial\lambda_2} \right |}_{\lambda_1,\lambda_2 = 0}
	\thickspace \text{.}
\label{eq:linresp}
\end{equation}
Formally, all quantities of interest can be calculated from
``normalizing constants'' of distributions, possibly modified in a
suitable manner, such as in \eqref{eq:modZ} above.
%
\cut{
For example, the expectation
of $\spin_{1} \spin_{2}$ with respect to $\probability$ can be
calculated using the following normalizing constant
$\normalconst(\lambda) = \int \negthinspace {d \myvector{\spin}}
\thinspace \exp(H + \lambda \spin_{1} \spin_{2})$, from which
\begin{equation}
	\left < \spin_{1} \spin_{2} \right > =  {\left . \frac{1}{\normalconst(\lambda)} \frac{\partial \normalconst(\lambda)}{\partial \lambda} \right |}_{\lambda = 0}
	\thickspace \text{.}
\end{equation}
}
\cut{
For example, the expectation
of $\spin_{1} \spin_{2}$ with respect to $\probability$ can be
calculated using the following normalizing constant
$\normalconst(\lambda_1,\lambda_2) = \int \negthinspace {d \myvector{\spin}}
\thinspace \exp(H + \lambda_1 \spin_{1} + \lambda_2 \spin_{2})$, from which
\begin{equation}
	\left < \spin_{1} \spin_{2} \right > =  {\left . \frac{1}{\normalconst(\lambda_1,\lambda_2)} \frac{\partial^2 \normalconst(\lambda_1,\lambda_2)}{\partial \lambda_1 \partial\lambda_2} \right |}_{\lambda_1,\lambda_2 = 0}
	\thickspace \text{.}
\label{eq:linresp}
\end{equation}
%
}
%
It is clear that all moments of $\probability$ can be generated by
differentiating \eqref{eq:modZ} in a similar manner, and can be as equally
easily obtained from the generating function approach, as from the
normalising constant approach. We prefer here the normalising constant
approach since this enables us  more easily to make contact with
results from statistical physics.
%In some cases, however, it may be
%more natural to consider the generating function, as in directed
%graphical models (??).
The normalising constant approach is also more
natural for the examples of undirected networks that we shall consider
in later sections.

In the following
sections, we layout how to approximate the normalising constant,
assuming that any additional terms in the normalising constant, of the
form $\ip{\boldlambda}{\myvector{\spin}}$ in equation \eqref{eq:modZ},
are absorbed into the definition of the potential $H$.

As we will see in \secref{sec:learning}, the
normalising constant also plays a central role in learning.  Unfortunately,
for many probability distributions that arise in applications,
calculation of the normalizing constant is intractable due to the high
dimensionality of the random variable $\myvector{\spin}$.  When
$\myvector{\spin}$ is a $N$--dimensional binary vector, naively at
least, a summation over $2^N$ states has to be performed to calculate
the normalizing constant.

However, sometimes it is possible to exploit the structure of the
distribution in order to reduce the computational expense.  A trivial
example is that of factorised models, $\probability(\myvector{\spin})
\propto \prod_{i} \Psi(\spin_{i})$, in which $\normalconst = \int
\negthinspace {d \myvector{\spin}} \thinspace \prod_{i}
\Psi(\spin_{i}) = \prod_{i} \int \negthinspace {d \spin_{i}}
\thinspace \Psi(\spin_{i})$, so that the computation scales only
linearly with $N$.  Similarly, in the case of a Gaussian distribution,
computation scales (roughly) with $N^{3}$.  These tractable
distributions have recently been exploited in the variational method
to approximate intractable distributions, see for
example~\cite{Saul96,Jaakkola97,barber-wiegerinck-98b,jordan-in-book-98}.
In the following section we will briefly review one of the most common
variational methods which exploits the Kullback-Leibler bound.

\section{The Kullback-Leibler Variational Bound}
\label{sec: bound}

Our aim in this section is to briefly describe the current
state-of-the-art in approximating the normalizing constant of an
intractable probability distribution. We denote the intractable
distribution of interest as $\probability_{1}$ to distinguish it from
an auxiliary distribution $\probability_{0}$, we will introduce.
Without loss of generality, we write the distribution over a vector of
random variables $\myvector{\spin}$, parameterized by a vector
$\myvector{\weight}$ as
%
\begin{equation}
    \probability_{1} (\myvector{\spin}; \myvector{\weight}) = \frac{ e^{H_{1} (\myvector{\spin}; \myvector{\weight}) }}{\normalconst_{1} (\myvector{\weight})}
	\thickspace \text{,}
\end{equation}
%
where the potential $H_{1} (\myvector{\spin}; \myvector{\weight})$ is
a known function of $\myvector{\spin}$ and $\myvector{\weight}$, \eg
in the case of a factorised model this would have the form of $H_{1} =
\sum_{k} \weight_{k} \spin_{k}$.  The corresponding (assumed
intractable) normalizing constant is given by
%
\begin{equation}
    \normalconst_{1} (\myvector{\weight}) = \int \negthinspace {d \myvector{\spin}} \thinspace e^{H_{1}(\myvector{\spin}; \myvector{\weight})}
	\thickspace \text{.}
\end{equation}
We will use a family of tractable distributions, $\probability_{0}$,
parameterized by $\myvector{\newweight}$ which we write as
\begin{equation}
    \probability_{0} (\myvector{\spin}; \myvector{\newweight}) = \frac{ e^{H_{0} (\myvector{\spin}; \myvector{\newweight}) }}{\normalconst_{0} (\myvector{\newweight})}
	\thickspace
	\text{ with }
	\thickspace
    \normalconst_{0} (\myvector{\newweight}) = \int \negthinspace {d \myvector{\spin}} \thinspace e^{H_{0}(\myvector{\spin}; \myvector{\newweight})}
\label{eq:deftractq}
\end{equation} 
where, by assumption, the normalizing constant $\normalconst_{0}$ is tractably computable.

Standard variational methods attempt to find the best approximating
distribution $\probability_{0} (\myvector{\spin};
\myvector{\newweight^{*}})$ that matches $\probability_{1}
(\myvector{\spin}; \myvector{\weight})$ by minimizing the
Kullback-Leibler divergence,
%
\begin{equation}
	\KL (\probability_{0},\probability_{1}) = \negthickspace \int \negthinspace {d \myvector{\spin}} \thinspace \probability_{0} (\myvector{\spin}; \myvector{\newweight}) \log \frac{\probability_{0} (\myvector{\spin}; \myvector{\newweight})}{\probability_{1}(\myvector{\spin}; \myvector{\weight})}
	\thickspace \text{.}
\label{eq:kldiv}
\end{equation}
%
Since, using Jensen's inequality, $\KL
(\probability_{0},\probability_{1}) \geq 0$, we immediately obtain the
lower bound
%
\begin{equation}
\log \normalconst_{1} \geq - \int \negthinspace {d \myvector{\spin}} \thinspace \probability_{0} (\myvector{\spin}; \myvector{\newweight}) \log \probability_{0} (\myvector{\spin}; \myvector{\newweight}) +
\int \negthinspace {d \myvector{\spin}} \thinspace \probability_{0} (\myvector{\spin}; \myvector{\newweight}) H_{1} (\myvector{\spin}; \myvector{\weight})
	\thickspace \text{.}
\label{eq: variational bound}
\end{equation}
%
Using the definition \eqref{eq:deftractq} we obtain the bound in the more intuitive form,
%
\begin{equation}
\log  \normalconst_{1} \geq \log  \normalconst_{0} + \int \negthinspace {d \myvector{\spin}} \thinspace \probability_{0} (\myvector{\spin}; \myvector{\newweight}) \br{H_{1} (\myvector{\spin};\bfw) - H_{0} (\myvector{\spin};\myvector{\newweight})}
\end{equation}
%
The lower bound is made as tight as possible by maximizing the
right-hand side with respect to the parameters $\myvector{\newweight}$
of $\probability_{0}$, which corresponds to minimising the
Kullback-Leibler divergence \eqref{eq:kldiv}.  This approach has
recently received  attention in the artificial intelligence
community and has been demonstrated to be a useful
technique~\cite{Saul96,Jaakkola97,barber-wiegerinck-98b}. Note that whilst this
lower bound is useful in providing an objective comparison of two
approximations to the normalising constant (the approximation with the
higher bound value is preferred), it does not translate directly into
a bound on conditional probabilities. An upper bound on the
normalising constant is also required in this case. Whilst, in some cases, this may be feasible, in general this tends to be a rather more difficult task, and we restrict our attention to the lower bound \cite{Jaakkola97}.

In the next section, we describe another approach that enables us to
exploit further such tractable distributions.

\section{The Variational Cumulant Expansion\label{sec: Approximation}}

In order to extend the Kullback-Leibler variational lower bound with a view to
improving its accuracy without much greater computational expense, we introduce
the following family of probability distributions $\probability_{\alpha}$, 
parameterized by $\alpha$, $\myvector{\weight}$ and $\myvector{\newweight}$:
\begin{equation}
    \probability_{\alpha} (\myvector{\spin}; \myvector{\weight}, \myvector{\newweight}) = \frac{ e^{H_{\alpha} (\myvector{\spin}; \myvector{\weight}, \myvector{\newweight}) }}{\normalconst_{\alpha} (\myvector{\weight}, \myvector{\newweight})}
	\thickspace
\end{equation}
where the potential is given by
\begin{equation}
    H_{\alpha}(\myvector{\spin}; \myvector{\weight}, \myvector{\newweight}) = \alpha H_{1} (\myvector{\spin}; \myvector{\weight}) + (1-\alpha) H_{0} (\myvector{\spin}; \myvector{\newweight})
	\thickspace \text{.}
\end{equation}
The probability distributions in this family 
interpolate between a tractable distribution, $\probability_{0}$
which is obtained for $\alpha = 0$, and an intractable distribution, $\probability_{1}$
which is obtained for $\alpha = 1$.
For intermediate values of $\alpha \in (0,1)$, the distributions
remain intractable.
The normalizing constant of a distribution from this family is given by
\begin{eqnarray}
	\log \normalconst_{\alpha} (\myvector{\weight}, \myvector{\newweight}) &=&
	\log \int \negthinspace {d \myvector{\spin}} \thinspace e^{H_{0} (\myvector{\spin}; \myvector{\newweight})
		+ \alpha \left ( H_{1} (\myvector{\spin}; \myvector{\weight}) - H_{0} (\myvector{\spin}; \myvector{\newweight}) \right )}\nonumber\\
	&=&
	\log \normalconst_{0} + \log {\left < e^{\alpha \left ( H_{1} (\myvector{\spin}; \myvector{\weight}) - H_{0} (\myvector{\spin}; \myvector{\newweight}) \right )} \right >}_{0}
\label{eq: cumulant generating function}
	\thickspace \text{,}
\end{eqnarray}
where ${\left < . \right>}_{0}$ denotes expectation with respect to
the tractable distribution $\probability_{0}$.  The final term in
\eqref{eq: cumulant generating function} is intractable for any
$\alpha\neq 0$, and we shall therefore develop a Taylor series
expansion for it around $\alpha=0$.  That is,
%
\begin{multline}
\log {\left < e^{\alpha \left ( H_{1} (\myvector{\spin}; \myvector{\weight}) - H_{0} (\myvector{\spin}; \myvector{\newweight}) \right )} \right >}_{0} \approx \alpha
\av{H_{1} (\myvector{\spin}; \myvector{\weight}) - H_{0} (\myvector{\spin})}_0 \\
+\frac{\alpha^2}{2}\br{{H_{1} (\myvector{\spin}; \myvector{\weight}) - H_{0} (\myvector{\spin})} - {\av{H_{1} (\myvector{\spin}; \myvector{\weight})} + \av{H_{0} (\myvector{\spin})}_0}}^2 + O(\alpha^3)
\label{eq:taylorexpone}
\end{multline}
%
In terms of the random variable $\Delta H \equiv H_{1}
(\myvector{\spin}; \myvector{\weight}) - H_{0} (\myvector{\spin};
\myvector{\newweight})$, we see that \eqref{eq:taylorexpone} contains
the first terms in a power series in $\Delta H$, the first being the
mean, and the second the variance of $\Delta H$.  More formally, we
recognize the final term in \eqref{eq: cumulant generating function}
as the cumulant generating function $K_{\Delta H} (\alpha) \equiv \log
{\left < e^{\alpha \Delta H} \right >}_{0}$, with respect to  the random variable
%$\Delta H \equiv H_{1} (\myvector{\spin}; \myvector{\weight}) - H_{0}
%(\myvector{\spin}; \myvector{\newweight})$~\cite{Grimmett92}.  If the
$\Delta H$ ~\cite{Grimmett92}.  If the
logarithm of the moment generating function ${\left < e^{\alpha \Delta
H} \right >}_{0}$ is finite in the neighborhood of the origin then the
cumulant generating function $K_{\Delta H} (\alpha)$ has a convergent
Taylor series,
\begin{equation}
	K_{\Delta H} (\alpha) = \sum_{n=1}^{\infty} \frac{1}{n!} k^{0}_{n} (\Delta H) \alpha^{n}
	\thickspace \text{,}
\end{equation}
where $k^{0}_{n} (\Delta H)$ is the $n$th cumulant of $\Delta H$ with
respect to the distribution $\probability_{0}$. (The
definition of the $n$th cumulant is $k^{0}_{n} (\Delta H) =
\frac{\partial^n}{\partial \alpha^n} \log {\left < e^{\alpha \Delta H}
\right >}_{0} \vline_{\alpha=0}$). For convenience we have no longer
denoted the explicit dependence on the parameters
$\myvector{\newweight}$ and $\myvector{\weight}$.
%and the random
%variables $\myvector{\spin}$.

Since the intractable normalizing constant $\normalconst_{1}$ corresponds to setting $\alpha = 1$
in \eqn{\ref{eq: cumulant generating function}}, we can write the Taylor series as
\begin{equation}
	\log \normalconst_{1} = \log \normalconst_{0} + 
			\sum_{n=1}^{l} \frac{1}{n!} k^{0}_{n} (\Delta H)
			+ 
			\frac{1}{(l+1)!} k^{\xi}_{l+1} (\Delta H)
			\thickspace \text{,}
\label{eq: exact}
\end{equation}
where the remainder follows from the
the Mean Value Theorem (see, for example \citeA{Spivak67}),
and $k^{\xi}_{l+1} (\Delta H)$ 
is the $(l+1)$th cumulant of $\Delta H$ with respect to $\probability_{\xi}$, 
$0 \leq \xi \leq 1$.
An approximation is obtained by simply neglecting this (intractable) remainder and
using the
truncated expansion:
\begin{equation}
	\log \normalconst_{1} \approx \log \normalconst_{0} + 
			\sum_{n=1}^{l} \frac{1}{n!} k^{0}_{n} (\Delta H)
	\thickspace \text{,}
\label{eq: approx}
\end{equation}
Note that the right-hand side of \eqn{\ref{eq: approx}} depends on the free parameters $\myvector{\newweight}$,
and we shall discuss later how these parameters can be set to improve the accuracy of the approximation.

In the following two subsections, we will look in more detail at the first
and second-order approximation of this Taylor series, since these illuminate the
variational method and the difficulties in retaining a lower bound.

\subsection{First Order and the Variational Bound}
\label{sec: first order}

	The Taylor series (\eqn{\ref{eq: exact}}) up to first order yields
	\begin{equation}
		\label{eq: first order}
		\log \normalconst_{1} = \log \normalconst_{0} + {\left < \Delta H \right>}_{0}
								 + \frac{1}{2} k^{\xi}_{2} (\Delta H)
		\thickspace \text{,}
	\end{equation}
	where the remainder $k^{\xi}_{2} (\Delta H)$ is equal to the variance 
	(the second cumulant) of $\Delta H$ with respect to $\probability_{\xi}$ with
	$0 \leq \xi \leq 1$,
	\ie
	$k^{\xi}_{2} = \var_{\xi} ( \Delta H ) =
	 {\left < {\Delta H}^{2} \right >}_{\xi} 
	- {{\left < \Delta H \right>}_{\xi}}^{2}$,
	and
	${\left < . \right>}_{\xi}$ denotes expectation with respect to
    the distribution $\probability_{\xi}$. 
	Since the variance is nonnegative for any value of $\xi$, we have
	\begin{equation}
\label{eq: bound}
		\log \normalconst_{1} \geq \log \normalconst_{0} + {\left < \Delta H \right>}_{0}
		\thickspace \text{,}
	\end{equation}
	so that $\log \normalconst_{1}$ is not only approximated but also bounded from
	below by the right hand side of \eqn{\ref{eq: bound}}.
	This bound is the standard mean field bound (see also
	\eqn{\ref{eq: variational bound}}),
	and is employed in many variational methods
	\cite<see, for example>{Saul96,Jaakkola97,barber-wiegerinck-98b}.
	This bound can then be made as tight as possible by optimizing
	with respect to the free parameters $\myvector{\newweight}$
	of the tractable distribution $\probability_{0}$.

\subsection{Second and Higher Order}
\label{sec: second order}

	To second order \eqn{\ref{eq: exact}} yields \begin{equation}
	\label{eq: second order} \log \normalconst_{1} = \log
	\normalconst_{0} + {\left < \Delta H \right>}_{0} +
	\frac{1}{2} \var_{0} ( \Delta H) + \frac{1}{6} k^{\xi}_{3}
	(\Delta H) \thickspace \text{,} \end{equation} with the
	remainder $k^{\xi}_{3} \negthinspace = \negthinspace {\left <
	{\Delta H}^{3} \right>}_{\xi} \negmedspace - 3 {\left <
	{\Delta H}^{2} \right>}_{\xi} {\left < \Delta H \right>}_{\xi}
	\negmedspace + 2 {\left < \Delta H \right>}_{\xi}^{3}$ and $0
	\negthinspace \leq \negthinspace \xi \negthinspace \leq
	\negthinspace 1$.  Since this remainder cannot be bounded in a
	tractable manner for general distributions $\probability_{0}$
	and $\probability_{1}$, we have no obvious criterion to prefer
	one set of parameters $\myvector{\newweight}$ of the tractable
	distribution over another.  Similarly, for higher-order
	expansions it is unclear, within this expansion, how to obtain
	a tractable bound and, consequently, how to select a set of
	parameters $\myvector{\newweight}$.  We stress that any
	criterion is essentially a heuristic since the error made by
	the resulting approximation is, by assumption, not tractably
	computable. Whilst, in principle, a whole range of criteria
	could be considered in this framework, we mention here only
	two.  The first, the bound criterion, is arguably the simplest
	and most natural choice.  The second criterion we briefly
	mention is the independence method, which we discuss mainly
	for its relation to the TAP method of statistical
	physics~\cite{Thouless77}.

	\subsubsection{Bound Criterion}

	The first criterion assigns the free parameters $\myvector{\newweight}$ 
	to maximize the (first-order) bound, \eqn{\ref{eq: bound}}.
	The resulting distribution
	$\probability_{0} (\myvector{\spin}; \myvector{\newweight^{*}})$ is used
	in the truncated approximation, \eqn{\ref{eq: approx}}.
	To second order this is explicitly given by
	\begin{equation}
        \log \normalconst_{1} \approx  \log \normalconst_{0}^{*} + {\left < \Delta H \right>}^{*}_{0}
								 + \frac{1}{2} \var_{0}^{*} (\Delta H)
\label{eq: MI}
		\thickspace \text{,}
	\end{equation}
	where a `$*$' denotes that the distribution
	$\probability_{0} (\myvector{\spin}; \myvector{\newweight^{*}})$ is to be used.
	Note that the computational cost involved over that of the variational
	method of \sect{\ref{sec: first order}} is small.
	No further optimization is required, only the calculation of the
	second cumulant $\var_{0}^{*} (\Delta H)$ which by assumption is tractable.
	%Similarly, including higher-order cumulants incurs only a small computational cost.

	\subsubsection{Independence Criterion}
	
	The second criterion we consider is based on the fact that 
	the exact value of $\normalconst_{1}$ is independent of the free parameters
	$\myvector{\newweight}$, \ie
	\begin{equation}
\label{eq: tap condition}
		\frac{\partial \log \normalconst_{1}}{\partial \newweight_{i}} \equiv 0
		\thickspace \text{.}
	\end{equation}
	Since $\log \normalconst_{1}$ is intractable, we substitute the 
	truncated cumulant expansion into 
	\eqn{\ref{eq: tap condition}}, 
	\begin{equation}
		\frac{\partial}{\partial \newweight_{i}}
		 \left ( \log \normalconst_{0} + \sum_{n=1}^{l} \frac{1}{n!} k^{0}_{n} (\Delta H) \right )
		\equiv 0
		\thickspace \text{.}
	\end{equation}
	In general, there are many solutions resulting from the
	independence criterion and alone it is too weak to provide
	reliable solutions.
	Under certain restrictions, however, this approach is closely
	related to the TAP approach of statistical
	mechanics~\cite{Thouless77,Plefka82,Kappen98,Kappen98b},
	as will also be described in \sect{\ref{sec:TAP}}.

\section{Boltzmann Machines}
\label{sec: BM}

To illustrate the foregoing theory, we will consider a class of
discrete undirected graphical models, Boltzmann
machines~\cite{Ackley85}.  These have application in artificial
intelligence as stochastic connectionist
%models~\cite{Russell95,Jordan98}, in image restoration~\cite{Geman84},
models~\cite{jordan-in-book-98}, in image restoration~\cite{Geman84},
and in statistical physics~\cite{Itzykson89}.  The potential of a
Boltzmann machine, with binary random variables $\spin_{i} \in
\{0,1\}$, is given by
%
\begin{equation}
\label{eq: intractable}
H(\myvector{\spin}; \myvector{\weight}) = \sum_{i} \weight_{i} \spin_{i} + \frac{1}{2} \sum_{i,j} \weight_{ij} \spin_{i} \spin_{j}
\end{equation}
%
with $\weight_{ij} \equiv \weight_{ji}$ and $\weight_{ii} \equiv 0$.
Unfortunately, Boltzmann machines are in general intractable since
calculation of the normalizing constant involves a summation over an
exponential number of states.
%Formally, this can be proved to be $NP$-hard~\cite{Cooper90}.
In an early effort to overcome this
intractability, \citeA{Peterson87} proposed the variational
approximation using a factorised model as a fast alternative to
stochastic sampling.  However, \citeA{Galland93} pointed out that this
approach often fails since it inadequately captures second
order statistics of the intractable distribution.  Using the
potentially more accurate TAP approximation did not overcome these
difficulties and often lead to unstable solutions
\cite{Galland93}. Furthermore, it is not clear how to extend the TAP
approach to deal with using more complex, non-factorised approximating
distributions. We examine the relationship of our approach to the TAP
method in \secref{sec:TAP}. Another approach which aims to improve the
accuracy of the correlations, though not the normalising constant, is
linear response \cite{parisi-book-88} which was applied to Boltzmann
machines by \citeA{Kappen98} for the case of using a factorised model,
see \secref{sec:Factorised Boltzmann machines}. The approach we take
here is a little different in that we wish to find a better
approximation primarily to the normalising constant. Approaches such
as linear response can then be applied to this approximation to derive
approximations for the correlations, if desired.

More flexible approximating distributions have also been considered in
the variational approach, for example, mixtures of factorised
distributions~\cite{Jaakkola98,Lawrence98}, and decimatable
structures~\cite{Saul94,barber-wiegerinck-98b}.  We show how
to combine the use of such more flexible approximating distributions
with higher-order corrections to the variational procedure.  For
clarity and simplicity, we will initially approximate the normalizing
constant of a general intractable Boltzmann machine
$\probability_{1}$, see \eqn{\ref{eq: intractable}}, using only a
factorised Boltzmann machine as our tractable model
$\probability_{0}$.  Subsequently, in our simulations, we will use a
more flexible tractable model, a decimatable Boltzmann
machine.

\subsection{Using Factorised Boltzmann Machines}
\label{sec:Factorised Boltzmann machines}

The potential of a factorised Boltzmann machine is given by
%
\begin{equation}
	H_{0} = \sum_{i} \newweight_{i} \spin_{i}
	\label{eq: factorised}
        \thickspace \text{.}
\end{equation}
%
\noindent
For this simple model, calculating higher-order cumulants is straightforward once the value of the free
parameters $\myvector{\newweight}$ are known since,
%
\begin{equation}
		{\left < \spin_{i_{1}} \spin_{i_{2}} \ldots \spin_{i_{l}} \right >}_{0}
			= 
			{\left < \spin_{i_{1}} \right >}_{0}
			{\left < \spin_{i_{2}} \right >}_{0}
			\ldots
			{\left < \spin_{i_{l}} \right >}_{0}
			\thickspace
			\text{ with }
			\thickspace
			i_{1} \neq i_{2} \neq \ldots \neq i_{l}
		\thickspace \text{.}
	\end{equation}
	For notational convenience, we define $\meanfield_{i}$ to be
	the expectation value of $\spin_{i}$,
\begin{equation}
\meanfield_{i} = {\left < \spin_{i} \right >}_{0} = {\left (1+\exp\left(-\newweight_{i}\right) \right )}^{-1}.
\label{eq:meandef}
\end{equation}

\subsubsection{Bound Criterion}

Using a factorised distribution to approximate the intractable
normalizing constant corresponding to \eqn{\ref{eq: intractable}}
gives, from \eqn{\ref{eq: bound}}, the variational bound
%
\begin{equation}
\label{eq: bound fact}
    \log \normalconst_{1} \geq  \entropy(\myvector{\meanfield})
					   + \sum_{i} \weight_{i} \meanfield_{i} 
						+ \frac{1}{2} \sum_{i,j} \weight_{ij} \meanfield_{i} \meanfield_{j}
	    \thickspace \text{,}
\label{eq:varboundone}
\end{equation}
	where, for convenience, we have defined the entropy
	\begin{equation}
	    \entropy(\myvector{\meanfield}) = -\sum_{i} \left \{ \meanfield_{i} \log \meanfield_{i} + \left (1-\meanfield_{i} \right ) \log \left ( 1-\meanfield_{i} \right ) \right \}
	\thickspace \text{.}
	\end{equation}
	One can either maximize the bound, \ie the right-hand side of
	\eqn{\ref{eq: bound fact}}, directly, or use
	the fact that at the maximum, 
	\begin{equation}
\label{eq: derivative zero}
		\frac{\partial}{\partial \newweight_{k}} 
		\left [ \entropy(\myvector{\meanfield})
                       + \sum_{i} \weight_{i} \meanfield_{i}
                        + \frac{1}{2} \sum_{i,j} \weight_{ij} \meanfield_{i} \meanfield_{j}
			  \right ]
			= 0
		\thickspace \text{,}
	\end{equation}
%
which leads to the fixed point equation (with the understanding that
the means $\myvector{m}$ are related to the parameters
$\myvector{\newweight}$ by equation \eqref{eq:meandef})
%
\begin{equation}
\label{eq: factorised mf equation}
		\newweight_{i} = \weight_{i} + \sum_{j} \weight_{ij} \meanfield_{j}.
%		\thickspace \text{,}
\label{eq:mfeqone}
\end{equation}
The well-known ``mean field'' equations \cite{parisi-book-88} are
a re-expression of \eqref{eq:mfeqone} in terms of the means
only, given by \eqref{eq:meandef}.  The optimal parameters of the
factorised model might therefore also be obtained by iterating, either
synchronously or asynchronously~\cite{Peterson87}, these fixed point
equations.  Note that \eqn{\ref{eq: derivative zero}} also holds for
minima and saddle-points, so in general additional checks are needed
to assert that the parameters correspond to a maximum of \eqref{eq:varboundone}.
%For this
%particular case, however, the fixed point corresponds to a
%maximum~\cite{Kurita96}.
Insertion of the fixed point solution into
the second-order expansion for the bound criterion, \eqn{\ref{eq:
MI}}, yields (up to second order)
%
\begin{equation}
		\label{eq: MI second order}
		\log \normalconst_{1} \approx  \entropy(\myvector{\meanfield})    
				 \sum_{i} \weight_{i} \meanfield_{i}
					\frac{1}{2} \sum_{i,j} \weight_{ij} \meanfield_{i} \meanfield_{j}
					\frac{1}{4} \sum_{i,j} {\weight_{ij}}^{2}
				\meanfield_{i} (1-\meanfield_{i}) \meanfield_{j} (1-\meanfield_{j})
%		\thickspace 
		\text{,}
\label{eq:holnZ}
\end{equation}
%
where we emphasize that the $\meanfield_{i}$ are given by the
variational bound solution, \ie \eqn{\ref{eq: factorised mf
equation}}.

\subsubsection{Independence Criterion}
\label{sec:TAP}

Using a factorised distribution to approximate the intractable
normalizing constant, see \eqn{\ref{eq: intractable}}, gives up to
second order
%
%\begin{equation}
\begin{multline}
%\hspace*{-8.5mm}
\log \normalconst_{1} \approx
\entropy(\myvector{\meanfield}) \negthinspace + \negthinspace \sum_{i}
\weight_{i} \meanfield_{i} \negthinspace + \negthinspace \frac{1}{2}
\sum_{i,j} \weight_{ij} \meanfield_{i} \meanfield_{j} \negthinspace +
\negthinspace \frac{1}{4} \sum_{i,j} {\weight_{ij}}^{2} \meanfield_{i}
\left ( 1 - \meanfield_{i} \right ) \meanfield_{j} \left ( 1 -
\meanfield_{j} \right )\\
+ \frac{1}{2} \sum_{i}
{\left(\newweight_{i} - \weight_{i} - \sum_{j} \weight_{ij}
\meanfield_{j} \right )}^{2} \meanfield_{i} \left ( 1 - \meanfield_{i}
\right ) \thickspace \text{.}
\label{eq: TAP+}
\end{multline}
%\end{equation}
%
The independence criterion leads to the equations
%
\begin{multline}
{\left ( \sum_{j} {\weight_{ij}}^{2} \meanfield_{j} \left ( 1 -
\meanfield_{j} \right ) + {\left(\newweight_{i} - \weight_{i} -
\sum_{j} \weight_{ij} \meanfield_{j} \right )}^{2} \right )} \left (
\frac{1}{2} - \meanfield_{i} \right )\\ = \sum_{j}
{\left(\newweight_{j} - \weight_{j} - \sum_{k} \weight_{jk}
\meanfield_{k} \right )} \weight_{ij} \meanfield_{j} \left ( 1 -
\meanfield_{j} \right ) \thickspace \text{.}
\end{multline}
%
In general, there are many solutions to these equations and the search
can be limited by inserting in \eqn{\ref{eq: TAP+}} the constraint
that the second-order solution is close to the first-order
solution%\footnote{ The implicit assumption being that the parameters
%$\myvector{\weight}$ are small.  }
, \ie $\newweight_{i} = \weight_{i}
+ \sum_{j} \weight_{ij} \meanfield_{j} + \order (\weight^{2})$, and by
neglecting terms of $\order (\weight^{4})$ and higher.  This
simplifies \eqn{\ref{eq: TAP+}} to $-F_{\text{TAP}}$, the negative TAP
free energy\footnote{ In \citeA{Thouless77}, \citeA{Plefka82} and
\citeA{Kappen98} the random variables $\spin_{i} \in
\{-1,1\}$.}~\cite{Thouless77,Plefka82,Kappen98}.  The independence
criterion, \ie $\partial F_{\text{TAP}}/{\partial \newweight} = 0$,
leads now to the fixed point condition \begin{equation} \newweight_{i}
= \weight_{i} + \sum_{j} \weight_{ij} \meanfield_{j} + \frac{1}{2}
\sum_{j} {\weight_{ij}}^{2} (1-2\meanfield_{i}) \meanfield_{j}
(1-\meanfield_{j}) \thickspace \text{,} \end{equation} These TAP
equations can have poor convergence properties and often produce a
poor solution~\cite{Bray79,Nemoto85,Galland93}.  The additional
constraint that the TAP solution should correspond to a minimum of
$F_{\text{TAP}}$ improves the solution.  However, since TAP is
therefore essentially an expansion around the first-order solution, we
expect the numerical difference between the bound criterion and
convergent TAP results to be small.  For these reasons, we prefer the
straightforward bound criterion and leave other criteria for separate
study.
%\footnote{Note that in certain cases, in particular, the SK model,
%the TAP approach can be shown to yield exact results, whereas the naive mean field theory provides
%only an approximation.}


\subsection{Numerical Results}
\label{sec: Simulations}

In this section we present results of computer simulations to validate
our approach.  In \sect{\ref{sec: approx normalizing constant}} we
will compare the different methods of \sect{\ref{sec: first order}}
and \sect{\ref{sec: second order}} on approximating the normalizing
constant $\normalconst$. A similar comparison is made for
approximating the correlations in \secref{sec:hocorr}. In
\sect{\ref{sec:learning}} we show the results of a learning problem in
which a Boltzmann machine is used to learn a set of sample patterns, a
typical machine learning problem.  In all cases, the numbers of random
variables in the Boltzmann machines are chosen to be small to
facilitate comparisons with the exact results.

In some simulations we will not only use factorised Boltzmann machines
but also more flexible models that can make the variational bound (see
\eqn{\ref{eq: bound}}) tighter, namely decimatable Boltzmann
machines~\cite{Saul94,rueger-97,barber-wiegerinck-98b}.  Decimation is a technique
that eliminates random variables so that the normalizing constant for
the distribution on the remaining variables remains unchanged up to a
constant known factor.  For completeness, we briefly describe this
technique in the current context.

	\begin{figure}
		\begin{center}
		\subfigure[With $\spin_{1}$]
		{
		\setlength{\unitlength}{1pt}
		\begin{picture}(120,90)(-60,-20)
			%neurons
			\put(-42,  0){\circle{14}}
			\put( 42,  0){\circle{14}}
			\put(  0, 42){\circle{14}}
			%labels
			\put(  0, 42){\makebox(0,0){$\spin_1$}}
			\put(-42,  0){\makebox(0,0){$\spin_2$}}
			\put( 42,  0){\makebox(0,0){$\spin_3$}}
			\put(  0, 57){\makebox(0,0){$\newweight_{1}$}}
			\put(-26, 21){\makebox(0,0){$\newweight_{12}$}}
			\put( 26, 21){\makebox(0,0){$\newweight_{13}$}}
			\put(-42,-15){\makebox(0,0){$\newweight_{2}$}}
			\put( 42,-15){\makebox(0,0){$\newweight_{3}$}}
			\put(  0,- 8){\makebox(0,0){$\newweight_{23}$}}
			%connections
			\put(-49,  0){\line( -1, 1){10}}
			\put(-49,  0){\line( -1,-1){10}}
			\put(-35,  0){\line( 1, 0){70}}
			\put(-35,  0){\line( 1, 1){35}}
			\put( 35,  0){\line(-1, 1){35}}
			\put( 49,  0){\line(  1, 1){10}}
			\put( 49,  0){\line(  1,-1){10}}
		\end{picture}
		}
		\quad
		\quad
		\subfigure[Without $\spin_{1}$]
		{
		\setlength{\unitlength}{1pt}
		\begin{picture}(120,90)(-60,-20)
			%neurons
			\put(-42,  0){\circle{14}}
			\put( 42,  0){\circle{14}}
			%labels
			\put(-42,  0){\makebox(0,0){$\spin_2$}}
			\put( 42,  0){\makebox(0,0){$\spin_3$}}
			\put(-42,-15){\makebox(0,0){$\decweight_{2}$}}
			\put( 42,-15){\makebox(0,0){$\decweight_{3}$}}
			\put(  0,- 8){\makebox(0,0){$\decweight_{23}$}}
			%connections
			\put(-49,  0){\line( -1, 1){10}}
			\put(-49,  0){\line( -1,-1){10}}
			\put(-35,  0){\line( 1, 0){70}}
			\put( 49,  0){\line(  1, 1){10}}
			\put( 49,  0){\line(  1,-1){10}}
		\end{picture}
		}
		\end{center}
		\caption{
			The Boltzmann machine with the random variable $\spin_1$ can be decimated
			to a Boltzmann machine without $\spin_1$.
		}
\label{fig: decimation}
	\end{figure}

	Suppose we have a Boltzmann machine with many random variables,
	for which a three-variable subgraph is depicted
	in \subfig{\ref{fig: decimation}}{a},
	so that random variable $\spin_1$ is connected only to random variables
	$\spin_2$ and $\spin_3$.
	Mathematically, the condition for invariance of $\normalconst$ after removing 
	random variable $\spin_1$, \subfig{\ref{fig: decimation}}{b}, becomes
	\begin{equation}
		\sum_{\myvector{\spin} \backslash \spin_{1}}
		 e^{
		 \decweight +
	  	 \decweight_{2} \spin_{2} + 
		 \decweight_{3} \spin_{3} + 
		 \decweight_{23} \spin_{2} \spin_{3}}
		 = 
		\sum_{\myvector{\spin}}
		e^{
			\newweight +
			\newweight_{1} \spin_{1} + 
			\newweight_{2} \spin_{2} + 
			\newweight_{3} \spin_{3} + 
			\newweight_{12} \spin_{1} \spin_{2} + 
			\newweight_{13} \spin_{1} \spin_{3} +
			\newweight_{23} \spin_{2} \spin_{3}
		}
		\thickspace 
		\text{,}
	\end{equation}
	where we include the constant $\newweight$ for convenience.
	Invariance is fulfilled if
	\begin{equation}
	\begin{split}
	\decweight_{2} &= \newweight_{2} + \log \frac{1+e^{\newweight_{1} + \newweight_{12}}}{1 + e^{\newweight_{1}}}
	\thickspace \text{,} \thickspace
	\decweight_{3} = \newweight_{3} + \log \frac{1+e^{\newweight_{1} + \newweight_{13}}}{1 + e^{\newweight_{1}}}
	\thickspace \text{,}\\
	\decweight_{23} &= \newweight_{23} + \log \frac{\left (1 + e^{\newweight_{1}} \right ) \left (1+ e^{\newweight_{1}+ \newweight_{12} + \newweight_{13}}\right)}{\left( 1+e^{\newweight_{1} + \newweight_{12}} \right ) \left (1+e^{\newweight_{1} + \newweight_{13}} \right)}
	\thickspace \text{ and } \thickspace
	\decweight = \newweight + \log \left (1 + e^{\newweight_{1}}\right)
	\thickspace \text{.}
	\end{split}
	\end{equation}

For decimatable Boltzmann machines, one can repeat the decimation
process until, finally, one ends up with a Boltzmann machine which
consists of a single random variable whose normalizing constant is
trivial to compute.  This means that for decimatable Boltzmann
machines the normalizing constant can be computed in time linear in
the number of variables. Decimation in undirected models is similar in
spirit to variable elimination schemes in graphical
models. For example, in directed belief networks, ``bucket elimination''
enables variables to be eliminated by passing compensating messages to
other buckets (collections of nodes) in such that the desired marginal
on the remaining nodes remains the same \cite{dechter-98}.

\subsubsection{Approximating the Normalizing Constant}
\label{sec: approx normalizing constant}

\begin{figure}
	\begin{center}
		\subfigure[factorised model]
		{
		\setlength{\unitlength}{0.6pt}
		\begin{picture}(150,150)(-75,-75)
			%neurons
			\put(-50,-50){\circle{13}}
			\put(-50,-50){\circle{14}}
			\put(-50,-50){\circle{15}}
			\put(-50, 50){\circle{13}}
			\put(-50, 50){\circle{14}}
			\put(-50, 50){\circle{15}}
			\put( 50,-50){\circle{13}}
			\put( 50,-50){\circle{14}}
			\put( 50,-50){\circle{15}}
			\put( 50, 50){\circle{13}}
			\put( 50, 50){\circle{14}}
			\put( 50, 50){\circle{15}}
			\put(-71,  0){\circle{13}}
			\put(-71,  0){\circle{14}}
			\put(-71,  0){\circle{15}}
			\put(  0,-71){\circle{13}}
			\put(  0,-71){\circle{14}}
			\put(  0,-71){\circle{15}}
			\put( 71,  0){\circle{13}}
			\put( 71,  0){\circle{14}}
			\put( 71,  0){\circle{15}}
			\put(  0, 71){\circle{13}}
			\put(  0, 71){\circle{14}}
			\put(  0, 71){\circle{15}}
		\end{picture}
		}
		\quad
		\subfigure[decimatable]
		{
		\setlength{\unitlength}{0.6pt}
		\begin{picture}(150,150)(-75,-75)
			%neurons
			\put(-50,-50){\circle{13}}
			\put(-50,-50){\circle{14}}
			\put(-50,-50){\circle{15}}
			\put(-50, 50){\circle{13}}
			\put(-50, 50){\circle{14}}
			\put(-50, 50){\circle{15}}
			\put( 50,-50){\circle{13}}
			\put( 50,-50){\circle{14}}
			\put( 50,-50){\circle{15}}
			\put( 50, 50){\circle{13}}
			\put( 50, 50){\circle{14}}
			\put( 50, 50){\circle{15}}
			\put(-71,  0){\circle{13}}
			\put(-71,  0){\circle{14}}
			\put(-71,  0){\circle{15}}
			\put(  0,-71){\circle{13}}
			\put(  0,-71){\circle{14}}
			\put(  0,-71){\circle{15}}
			\put( 71,  0){\circle{13}}
			\put( 71,  0){\circle{14}}
			\put( 71,  0){\circle{15}}
			\put(  0, 71){\circle{13}}
			\put(  0, 71){\circle{14}}
			\put(  0, 71){\circle{15}}
			%connections
			\put(-43,-50){\line( 1, 0){86}}
			\put(-50,-43){\line( 0, 1){86}}
			\put(-43, 50){\line( 1, 0){86}}
			\put( 50,-43){\line( 0, 1){86}}
			\put(-45,-45){\line( 1, 1){90}}
		%	\put(-45, 45){\line( 1,-1){90}}
		%	\put(-64,  0){\line( 1, 0){128}}
		%	\put(  0,-64){\line( 0, 1){128}}
		%	\put(-66,  5){\line( 1, 1){61}}
		%	\put(-66, -5){\line( 1,-1){61}}
			\put(-44,-54){\line( 3,-1){38}}
			\put(-54,-44){\line(-1, 3){13}}
			\put(-44, 54){\line( 3, 1){38}}
			\put(-54, 44){\line(-1,-3){13}}
			\put( 44,-54){\line(-3,-1){38}}
			\put( 54,-44){\line( 1, 3){13}}
			\put( 44, 54){\line(-3, 1){38}}
			\put( 54, 44){\line( 1,-3){13}}
		%	\put(-44,-46){\line( 5, 2){108}}
		%	\put(-46,-44){\line( 2, 5){43}}
		%	\put( 44,-46){\line(-5, 2){108}}
		%	\put( 46,-44){\line(-2, 5){43}}
		%	\put( 44, 46){\line(-5,-2){108}}
		%	\put( 46, 44){\line(-2,-5){43}}
		%	\put(-44, 46){\line( 5,-2){108}}
		%	\put(-46, 44){\line( 2,-5){43}}
		\end{picture}
		}
		\quad
		\subfigure[fully connected]
		{
		\setlength{\unitlength}{0.6pt}
		\begin{picture}(150,150)(-75,-75)
			%neurons
			\put(-50,-50){\circle{13}}
			\put(-50,-50){\circle{14}}
			\put(-50,-50){\circle{15}}
			\put(-50, 50){\circle{13}}
			\put(-50, 50){\circle{14}}
			\put(-50, 50){\circle{15}}
			\put( 50,-50){\circle{13}}
			\put( 50,-50){\circle{14}}
			\put( 50,-50){\circle{15}}
			\put( 50, 50){\circle{13}}
			\put( 50, 50){\circle{14}}
			\put( 50, 50){\circle{15}}
			\put(-71,  0){\circle{13}}
			\put(-71,  0){\circle{14}}
			\put(-71,  0){\circle{15}}
			\put(  0,-71){\circle{13}}
			\put(  0,-71){\circle{14}}
			\put(  0,-71){\circle{15}}
			\put( 71,  0){\circle{13}}
			\put( 71,  0){\circle{14}}
			\put( 71,  0){\circle{15}}
			\put(  0, 71){\circle{13}}
			\put(  0, 71){\circle{14}}
			\put(  0, 71){\circle{15}}
			%connections
			\put(-43,-50){\line( 1, 0){86}}
			\put(-50,-43){\line( 0, 1){86}}
			\put(-43, 50){\line( 1, 0){86}}
			\put( 50,-43){\line( 0, 1){86}}
			\put(-45,-45){\line( 1, 1){90}}
			\put(-45, 45){\line( 1,-1){90}}
			\put(-64,  0){\line( 1, 0){128}}
			\put(  0,-64){\line( 0, 1){128}}
			\put(-66,  5){\line( 1, 1){61}}
			\put(-66, -5){\line( 1,-1){61}}
			\put( 66,  5){\line(-1, 1){61}}
			\put( 66, -5){\line(-1,-1){61}}
			\put(-44,-54){\line( 3,-1){38}}
			\put(-54,-44){\line(-1, 3){13}}
			\put(-44, 54){\line( 3, 1){38}}
			\put(-54, 44){\line(-1,-3){13}}
			\put( 44,-54){\line(-3,-1){38}}
			\put( 54,-44){\line( 1, 3){13}}
			\put( 44, 54){\line(-3, 1){38}}
			\put( 54, 44){\line( 1,-3){13}}
			\put(-44,-46){\line( 5, 2){108}}
			\put(-46,-44){\line( 2, 5){43}}
			\put( 44,-46){\line(-5, 2){108}}
			\put( 46,-44){\line(-2, 5){43}}
			\put( 44, 46){\line(-5,-2){108}}
			\put( 46, 44){\line(-2,-5){43}}
			\put(-44, 46){\line( 5,-2){108}}
			\put(-46, 44){\line( 2,-5){43}}
		\end{picture}
		}
	\end{center}
	\caption{
		Boltzmann machines of eight random variables	
	}
	\label{fig: eight neurons}
\end{figure}

As our `intractable' distribution $\probability_{1}$, we took a fully
connected, eight-node Boltzmann machine, which is not decimatable (see
\subfig{\ref{fig: eight neurons}}{c}).   As our tractable Boltzmann
machine, $\probability_{0}$, we took both the factorised model of
\subfig{\ref{fig: eight neurons}}{a} and the decimatable model of
\subfig{\ref{fig: eight neurons}}{b}.  We then tried to approximate
the normalizing constant for the fully connected machine in which the
parameters were randomly drawn from a zero-mean, Gaussian distribution
with unit variance.  See the appendix for additional theoretical and
experimental details regarding the optimization scheme used.  The
relative error of the approximation,
%
\begin{equation}
\label{eq: relative error}
E = \frac{\log \normalconst^{\text{exact}} - \log \normalconst^{\text{approx}}}{\log \normalconst^{\text{exact}}}
	\thickspace \text{,}
\end{equation}
was then determined for both the first and second-order approach using
both the factorised and decimatable model.  This was repeated 550
times for different random drawings of the parameters of the fully
connected Boltzmann machine.  The resulting histograms are depicted in
\fig{\ref{fig: error}}.  In these histograms one can clearly see that
the bound, while present in the first-order approach, is lost in the
second-order approach since some normalising constant estimates were
above the true value, violating the lower bound.  In five out of 1100
cases (two times 550) the synchronous mean field iterations did not
converge and these outliers were neglected in the plots.  It is,
however, possible to make the number of non-convergent cases even
lower by using asynchronous instead of synchronous updates of the mean
field parameters~\cite{Peterson87,Ansari95,Saul96}.  Note that, in
both cases, the second-order method improves on average on the
standard (first order) variational results.  Indeed, using only a
factorised model with the second-order correction improves on the
standard variational decimatable result.

\begin{figure}[t]
	\begin{center}
		\subfigure[Factorised first order. 
		The mean (absolute) error is $0.036$ ($0.034$ without outliers).]
		{\epsfig{figure=factorised_first.ps,width=0.35\linewidth}}
		\subfigure[Factorised second order.
		The mean absolute error is $0.0079$ ($0.0074$ without outliers).]
		{\epsfig{figure=factorised_second.ps,width=0.35\linewidth}}
		\subfigure[Decimatable first order.
		The mean (absolute) error is $0.0186$ ($0.0181$ without outlier).]
		{\epsfig{figure=decimatible_first.ps,width=0.35\linewidth}}
		\subfigure[Decimatable second order.
		The mean absolute error is $0.0031$ ($0.0027$ without outlier).]
		{\epsfig{figure=decimatible_second.ps,width=0.35\linewidth}}
	\end{center}
	\caption{
		Approximation of the normalizing constant of a fully connected, eight-node
		Boltzmann machine (see \subfig{\protect\ref{fig: eight neurons}}{c}),
		whose parameters are randomly drawn from 
		a zero-mean, Gaussian distribution with unit variance.
		A factorised model 
		(see \subfig{\protect\ref{fig: eight neurons}}{a})
		is used as the approximating distribution in (a) and (b)
		and a decimatable model
		(see \subfig{\protect\ref{fig: eight neurons}}{b})
		is used in (c) and (d).
		The histograms are based on 550 different random drawings
		of the parameters. 	
		For the readability of the histogram, we have excluded five outliers
		(four in (a) and (b), one in (c) and (d)),
		for which the synchronous mean field iteration did not converge. 
		The mean error is the mean of the absolute value of \eqn{\ref{eq: relative error}}.
	}
\label{fig: error}
\end{figure}

\begin{figure}[tbp]
	\begin{center}
		\subfigure[Factorised model.]
		{ \epsfig{figure=factorised_paired.eps,width=0.35\linewidth} }
		\subfigure[Decimatable model.]
		{ \epsfig{figure=decimatable_paired.eps,width=0.35\linewidth} }
	\end{center}
	\caption{
		Difference between first and second-order approximation
		of the normalizing constant of a fully connected, eight-node
		Boltzmann machine using (a) a factorised and (b) 
		a decimatable model.
		The histograms are based on the same 550 random drawings
		as used in \fig{\protect\ref{fig: error}}.
		All 1100 differences, including the five non-convergent runs,
		were larger than zero.
		Again for readability of the histogram, 
		the non-convergent runs are not included.
		The mean difference was
		0.0281 (0.0269 without the non-convergent runs)
		and 0.0155 (0.0154 without the non-convergent run)
		for the factorised model
		and the decimatable model, respectively.
	}
\label{fig: paired}
\end{figure}

We also determined whether the second-order approach improves the accuracy
of the prediction in each individual run.
In order to determine this we used the paired difference, 
\begin{equation}
\Delta = \left | \frac{\log \normalconst^{\text{exact}} - \log \normalconst^{\text{first}}}{\log \normalconst^{\text{exact}}} \right |
- \left | \frac{\log \normalconst^{\text{exact}} - \log \normalconst^{\text{second}}}{\log \normalconst^{\text{exact}}} \right |
\thickspace \text{.}
\end{equation}
If $\Delta$ is positive the second order improves the first order approach.
The corresponding results of our computer simulations are depicted 
in \fig{\ref{fig: paired}}.
In all runs $\Delta$ was positive, corresponding to a consistent improvement
in accuracy.

\subsubsection{Approximating the Means and Correlations}
\label{sec:hocorr}
There are many techniques that can be called upon to approximate
moments and correlations. Primarily, the cumulant expansion, as we
have described it, is intended to improve the accuracy in estimating
the normalising constant, and not directly the moments of the
distribution.  Arguably the simplest way to approximate correlations
is to use the standard variational approach to find the best
approximating distribution to $\probability_{1}$ and then to
approximate the moments of $\probability_{1}$ by the moments of
$\probability_{0}$. In this approach, however, certain correlations
that are not present in the structure of the approximating
distribution $\probability_{0}$, will be trivially approximated. For
example, approximating $\av{s_i s_j}_{\probability_{1}}$ using the
factorised distribution gives
$\av{s_is_j}_{\probability_{1}}\approx\av{s_i}_{\probability_{0}}\av{s_j}_{\probability_{0}}$. We
will examine in this section some ways of improving the estimation of
the correlations.

One procedure that can be used to approximate correlations such as
$\av{s_i s_j}_{\probability_{1}}$ is given by the so-called linear
response theory \cite{parisi-book-88}. This approach uses the
relationship \eqref{eq:linresp} in which the intractable normalising
constant is replaced with its approximation \eqref{eq:varboundone}.
%The derivatives in equation \eqref{eq:linresp}
%can be calculated is a self-consistent manner by analytically
%differentiating the mean field solution \eqref{eq: bound fact}.
\cut{
Similarly, one could use the higher-order approximation to the
normalising constant \eqref{eq:holnZ} in the linear response
formalism.
%% final corrections


calculate by defining another potential
\begin{equation}
{\tilde{H}}
_{\alpha}(\s;\weight,\theta,\boldgamma)
= \alpha \left[ H_1(\s;\weight) +
\sum_i \gamma_i s_i \right]
+ (1-\alpha)\left[ H_0(\s;\theta) +
\sum_i \gamma_i s_i \right]
\end{equation}
With this definition, the mean of a variable with respect to the
distribution defined by $H_1(\s;\weight)$ is given by
$\av{s_i}_{\probability_{1}} =
d\log {\tilde{Z}} / d\gamma_i$
evaluated at $\boldgamma = 0$. Using the
variational cumulant expansion, it is straightforward to show that
this gives the correction,
\begin{equation}
\av{s_i}_{\probability{1}} \approx \av{s_i}_{\probability_{0}} + \av{\s_i ( H_1 - H_0)}_{\probability_{0}} -
\av{\s_i}_{\probability_{0}} \av{ H_1 - H_0}_{\probability_{0}}
\end{equation}
where the parameters of $H_0$ have been chosen to minimize the KL
divergence between the distributions $\probability_{0}$ and
$\probability_1$.  The use of such correction terms to the first and
second moments of the fitted distribution is computationally
attractive since only one fit of the distribution defined by $H_0$ is
required. We lay aside here the study of such techniques, noting only
that the variational cumulant expansion readily enables such higher-order corrections to be evaluated.

%%



One often-made assumption, and which forms the basis for some approximation
techniques, is that the local fields are approximately Gaussian
distributed  \cite<see, for example>{saul-jordan-98}. In  the Boltzmann machine
case, this means $\sum_{j \neq i}J_{ij}s_j \sim {\cal
N}(\mu,\sigma)$. Such relations can be used, for example, in the
following manner:
%
In 0/1 notation, the mean of a variable is given by $m_i = Q(\s_i =
1)$. We can re-express this quantity using the identity,
%
\begin{align}
Q(\s_i = 1) & = \sum_{\sm{i}}Q(\s_i = 1 | \sm{i})Q(\sm{i})\\
& = \sum_{\sm{i}}\br{Q(\sm{i}|\s_i =0)Q(\s_i=0) + Q(\sm{i}|\s_i = 1)Q(\s_i=1)}Q(\s_i = 1 | \sm{i})
\end{align}
%
where $\sm{i}$ represents all the variables except variable $s_i$.
Using the definition of the mean, this gives the following exact
result\footnote{Using Bayes' rule, it is straightforward to re-express
equation \eqref{eq:cavity} in terms of one distribution alone, {\eg}
$Q(\sm{i}|\s_i =0)$. This is related to the cavity approach of
statistical physics \cite{mezard-book-87}.},
%
\begin{equation}
m_i = \frac{\sum_{\sm{i}} Q(\s_i = 1 | \sm{i})Q(\sm{i}|\s_i =0)}{1+\sum_{\sm{i}} Q(\s_i = 1 | \sm
{i})Q(\sm{i}|\s_i =0) - \sum_{\sm{i}} Q(\s_i = 1 | \sm{i})Q(\sm{i}|\s_i =1)}
\label{eq:cavity}
\end{equation}
%
For Boltzmann machines, $Q(\s_i = 1 | \sm{i}) = \sigma(\sum_{j \neq
i}J_{ij}s_j)$, from which we recognise that $m_i$ is determined by the
distribution of the variable $f_i = \sum_{j \neq i}J_{ij}s_j$,
conditional on the variable $s_i$ being 1 or 0, as in
\eqref{eq:cavity}. Assuming that $Q(\sm{i}|\s_i = 0/1)$ is well
approximated by a factorised distribution (or a mixture thereof),
approximating the mean and variance of the field distribution $f_i =
\sum_{j \neq i}J_{ij}s_j$ is straightforward, given the first two
moments of $s_j$. Assuming that the $f_i$ are Gaussian distributed,
the averages in equation \eqref{eq:cavity} can be computed
straightforwardly. Whilst this approach can produce accurate results
when the connections $J_{ij}$ are suitably random, this argument
generally breaks down in the case that we are most interested in,
namely learning, as in \secref{sec:learning} \cite<see
also>{barber-sollich-98a}. For this reason, we neglect such Gaussian
field approximations in favour of less constrained
assumptions\footnote{Whilst more advanced procedures can be considered
to improve on the Gaussian field assumption, it is not clear that this
would provide a valuable approach given limited computational
resources.}.
}
\cut{
Alternatively, one can use the approach outlined
above in which the correlations are given by ratios of normalising constant
functions, where at each stage, the normalising constants can be
approximated either by the standard mean field theory or by using the
higher-order corrections.

The question is, how are we going to approximate these correlations?
Lets look first at using a simple factorised model. We use the fact that
}

The approach that we adopt here, due to it's
straightforward relationship to normalising constants, is given by
considering the following relationship:
%
\begin{equation}
\av{s_i s_j}_{\probability_{1}} = \probability_{1}(s_i= 1, s_j =1) = e^{h_i + h_j + 2J_{ij}}\frac{Z_{-[i,j]}}{Z}
\label{eq:ratioZ}
\end{equation}
where $Z_{-[i,j]}$ is the normalising constant of the original network
in which nodes $i$ and $j$ have been removed and the biases are transformed to
$h_k' = h_k + J_{ik} + J_{jk}$.  This means that we need to
evaluate $O(n^2)$ normalising constants to approximate the
correlations. To approximate the normalising constants, we use the higher-order expansion \eqref{eq:holnZ}.
\cut{
Note that in the case of
using a factorial model, this is not really any more work since
everything reduces to single variable averages. Indeed, this is given
already by formula .. However, in more general cases, finding the
higher-order correlations that are required for the higher-order terms
may not be so computationally attractive. In the case of using mixtures
of factorised models, however, this is straightforward.
Here we are simply going to use factorial models.
Here are some results about approximating the correlations:
}
%
Similarly, we attempt to obtain a more accurate approximation to the
means $\av{s_i}$, based on the identity,
\begin{equation}
\av{s_i}_{\probability_{1}} = \probability_{1}(s_i= 1) = e^{h_i}\frac{Z_{-[i]}}{Z}
\label{eq:ratioZmean}
\end{equation}
where $Z_{-[i]}$ is the normalising constant of the original network
in which node $i$ has been removed and the biases are transformed to
$h_k' = h_k + J_{ik}$.  The ``intractable'' Boltzmann machine in our
simulations has 8 fully connected nodes with biases and weights drawn
from a standard zero-mean, unit-variance Gaussian distribution. In
figures 5 and 6 we plot results for approximating the means $\av{s_i}$
and correlations $\av{s_i s_j}$ respectively.  As can be seen, the approach we take to approximate the
correlations is an improvement on the standard variational factorised
model. We emphasize that, since we are primarily concerned with
approximating the normalising constant, there may be more suitable
methods dedicated to the task of approximating the correlations
themselves.  Indeed, one of the drawbacks of these perturbative
approaches is that physical constraints, such as that moments should
be bounded between 0 and 1, can be violated.
%Nevertheless, in the
%experiments considered here, using the higher-order approximations
%improves upon the standard approach of using the factorized
%variational approximation.

\begin{figure}[tbp]
\begin{center}
\subfigure[The mean absolute error in approximating the means $\av{s_i}$ is 0.0350.]
{\epsfig{figure=correl_mf1.eps,width=0.35\linewidth}}
\subfigure[The mean absolute error in approximating the means $\av{s_i}$ is 0.0129.]
{\epsfig{figure=correl_mf2.eps,width=0.35\linewidth}}
\end{center}
\label{fig:results_mf1}
\caption{Approximating the means $\av{s_i}$ for 1000 randomly chosen 8
node fully connected Boltzmann machines. The weights were drawn from
the standard normal distribution. (a) The standard variational
approach in which the means are estimated by the means of the best
approximating factorised distribution. (b) Using the ratio of
normalising constants \eqref{eq:ratioZ}, in which second-order
corrections are included in approximating the normalising constants. Using the ratio of normalising constants without the higher order corrections gave a mean error of 0.046, slightly worse than the standard result variational result.}
\end{figure}
%
\begin{figure}[tbp]
\begin{center}
\subfigure[Mean absolute error = 0.0329]
{\epsfig{figure=correl_mf3.eps,width=0.35\linewidth}}
\subfigure[Mean absolute error = 0.0103]
{\epsfig{figure=correl_mf4.eps,width=0.35\linewidth}}
\end{center}
\label{fig:results_mf2}
\caption{Approximating the correlations $\av{s_i s_j}$ for 1000 randomly chosen 8 node fully connected Boltzmann machines. The weights were drawn from the standard normal distribution. (a) The standard variational approach in which the means are estimated by the means of the best approximating factorised distribution, and (b) using the ratio of normalising constants, including second-order corrections. Without second-order corrections, this gives a mean error of 0.0428.}
\end{figure}


\subsubsection{Learning}
\label{sec:learning}

The results of the previous section show that, for randomly chosen
connection matrices $J$, the higher-order terms can significantly
improve the approximation to the normalising constant of the
distribution. However, in a learning scenario, such results may be of
little interest since the connection matrix will become intimately
related to the patterns to be learned---that is, $J$ will take on a
form that is appropriate for storing the patterns in the training
set as learning progresses.
%In this case, Gaussian field
%assumptions (see \secref{sec:hocorr}) are inappropriate
%\cite{barber-sollich-98a}.

%
\cut{
\subsection{No Hidden Units}
In our learning simulation, we used the teacher-pupil paradigm: we
trained our pupil, a fully connected Boltzmann machine of eight random
variables, on data generated by a fully connected Boltzmann machine of
the same size, the teacher.  The parameters of the teacher Boltzmann
machine were randomly drawn from a Gaussian distribution with zero
mean and a standard deviation of $0.5$.  The parameters of the pupil
were initialized using a Gaussian distribution with zero mean and a
standard deviation of $0.05$.  Naturally, in a real-world application
the training set is finite.  In our toy problem, however, we assume
that an infinite set of training patterns is available, drawn
according to the teacher distribution.  The advantage of this approach
is that the Kullback-Leibler divergence between the teacher and pupil
distributions should then tend to zero under the learning procedure.
The Kullback-Leibler divergence between the teacher and pupil is given
by
\begin{equation}
	\KL =  \sum_{\myvector{\spin}}
        \probability(\myvector{\spin};\myvector{\weight}_{\text{teacher}})
        \log
        \frac
		{\probability(\myvector{\spin};\myvector{\weight}_{\text{teacher}})}
		{\probability(\myvector{\spin};\myvector{\weight}_{\text{pupil}})}
	\thickspace \text{.}
	\label{eq: KL}
\end{equation}
Gradient descent yields the learning rules~\cite{Ackley85}, \ie
\begin{equation}
	\Delta \weight_{i} (t) = - \learningparameter \left (
						    {\left < \spin_{i} \right >}_{t} -
							{\left < \spin_{i} \right >}_{\text{teacher}}
						\right )
\end{equation}
and
\begin{equation}
	\Delta \weight_{ij} (t) = - \learningparameter \left (
						    {\left < \spin_{i} \spin_{j} \right >}_{t} -
							{\left < \spin_{i} \spin_{j} \right >}_{\text{teacher}}
						\right )
	\thickspace \text{,}
\end{equation}
where the learning rate $\learningparameter$ was chosen such that the
total change in parameters
was equal to $\frac{1/4}{1 + \log(t)}$, where $t$ denotes the time, \ie
the number of parameter updates~\cite{Heskes93}.
${\left < . \right >}_{\text{teacher}}$ is the expectation with respect to the teacher distribution
and ${\left < . \right >}_{t}$ is the expectation with respect to 
the pupil's parameters at time $t$.

Clearly the difficulty is in evaluating the expectations, 
${\left < \spin_{i}  \right >}_{t}$ and
${\left < \spin_{i} \spin_{j} \right >}_{t}$ 
(For this toy problem the teacher expectations are assumed given---in real-world problems these would need to be estimated from the data).
These are readily related to normalizing constant evaluations.
For example,
${\left < \spin_{7} \spin_{8} \right >}_{t}$ is given by
$e^{\weight_{7} + \weight_{8} + \weight_{78}} 
\frac{\normalconst_{78}}{\normalconst}$ with $\normalconst$
the normalizing constant of the eight neuron Boltzmann machine
and $\normalconst_{78}$ the normalizing constant of a six neuron fully connected 
Boltzmann machine whose parameters are given by
$\newweight_{i} = \weight_{i} + \weight_{i7} + \weight_{i8}$
and $\newweight_{ij} = \weight_{ij}$, with $i,j \in \{1,\ldots,6\}$.
$\normalconst$ and $\normalconst_{78}$ are then approximated using
the first and second-order approximations described in 
\sect{\ref{sec: Approximation}} using a factorised model. 
In \fig{\ref{fig: learning curve}} we have depicted the Kullback-Leibler
divergence (see \eqn{\ref{eq: KL}})
as a function of the number of parameter updates.

\begin{figure}
	\begin{center}
	\subfigure[]{ \epsfig{file=sim_learn_250_focus.eps,width=0.45\linewidth} }
	\subfigure[]{ \epsfig{file=sim_learn_250.eps,width=0.45\linewidth} }
	\end{center}
	\caption{
		The learning curves of a fully connected Boltzmann machine of
		eight random variables trained on data generated by
		a fully connected Boltzmann machine of the same size, using
		both the first and second-order approximation with a factorised model.
		In (a) the learning curves between parameter updates 20 and 50
		are shown. 
		In (b) the learning curves are shown for 250 parameter updates.
	}
\label{fig: learning curve}
\end{figure}

In agreement with previous observations~\cite{Lawrence98}, we find that
the first-order approximation is unstable during learning.
This instability is due to the fact that a Boltzmann distribution is,
in general, multimodal whereas the first-order approximation is unimodal.
To illustrate the `mode hopping' which is exhibited by the 
first-order approximation, we have depicted in 
\subfig{\ref{fig: evolution}}{a}
the evolution of the parameters, projected onto the first two
principal components, between updates 100 and 250.
The second-order approximation is, on the contrary, rather stable
during learning. 
As can be seen in \subfig{\ref{fig: evolution}}{b}
the parameters oscillate around a mean value,
with decreasing variance
due to the decrease of the learning rate over time.
This simply corresponds to asymptotic convergence of the parameters
and no mode hopping is observed.

\begin{figure}
	\begin{center}
	\subfigure[First order.]{
	\epsfig{file=evolution_first_100-250.eps,width=0.45\linewidth}}
	\subfigure[Second order.]{
	\epsfig{file=evolution_second_100-250.eps,width=0.45\linewidth}}
	\end{center}
	\caption{
		The evolution between parameter updates 100 and 250, of the parameters
		during training of a fully connected Boltzmann machine of eight random
		variables on data generated by a fully connected Boltzmann machine
		of the same size using the (a) first and (b) second-order
		approximation with a factorised model.
		For visualization, we reduced the dimensionality of the parameters
		by projecting onto the first two principal components,
		determined by the parameter values between updates 100 and 250. 
	}
\label{fig: evolution}
\end{figure}
}
%\subsection{Learning with Hidden Units}

\cut{
There are several approaches to the issue of dealing with hidden units.

\subsubsection{Direct Approach}

The most direct approach is given by considering the maximizing the
likelihood of a visible pattern directly.
%
\begin{align}
\pdiffo{\J} \ln \teacher{\vu}
 = \pdiffo{\J}\ln \int d\hu \teacher{\hu,\vu}\\
& = \int d\hu \teacher{\hu|\vu} \pdiffo{\J}\phi(\hu|\vu) - \int d\hu\vu \teacher{\hu,\vu} \pdiffo{\J}\phi(\hu,\vu)
\end{align}
%
In the case of BMs, since the energies are quadratic functions of the
random variables and linear in the parameters $\J$, then we need to
evaluate quadratic correlations in with respect to clamped and free
distributions.

The question is, how are we going to approximate these correlations?
Lets look first at using a simple factorised model. We use the fact that

\[
\av{s_i s_j} = p(s_i= 1, s_j =1) = e^{h_i + h_j + 2J_{ij}}\frac{Z_{-[i,j]}}{Z}
\]
where $Z_{-[i,j]}$ has node i,j removed and has transformed biases
$h_k' = h_k + 2(J_{ik} + J_{jk})$.  This means that we need to
evaluate (roughly) $n^2$ normalising constants to get the statistics. We
are going to use the higher-order terms to approximate the statistics.
Note that in the case of using a factorial model, this is not really
any more work since everything reduces to single variable
averages. Indeed, this is given already by formula .. However, in more
general cases, finding the higher-order correlations that are required
for the higher-order terms may not be so computationally attractive. In
the case of using mixtures of factorised models, however, this is
straightforward.

Here we are simply going to use factorial models.
}

%\subsection{Variational Learning}

%In the variational learning method, we do:

We are therefore interested here in training networks on a set of
visible patterns. That is, we split the nodes of the network into a set of
visible $\vu$ and a set of  hidden units $\hu$ with the aim of learning a set of
patterns $\vu^1\ldots \vu^p$ with a Boltzmann machine  $\qone(\vu,\hu)$.
By using the KL divergence between an approximating distribution
$\qtwo(\hu|\vu)$ on the hidden units and the conditional distribution
$\qone(\hu|\vu)$ we obtain the following bound
%
\begin{equation}
\ln \qone(\vu)\geq -\int \qtwo(\hu|\vu)\ln \qtwo(\hu|\vu) + \int \qtwo(\hu|\vu)\ln \qone(\hu,\vu)
\label{eq:vlbound}
\end{equation}
%
For the case of Boltzmann machines, $\ln \qone(\hu,\vu) = H(\hu,\vu) -
\ln Z$ in which $\ln Z$ is (assumed) intractable. In order to obtain
an approximation to this quantity, we therefore introduce a further
variational distribution, $\qzero(\vu,\hu)$ (not the same as
$\qtwo(\hu|\vu)$) which can be used as in \eqref{eq: variational
bound} to obtain a lower bound on $\ln Z$.
%%
%
%Using the standard
%variational approach gives a lower bound on $\ln Z$ which,
Unfortunately, used in \eqref{eq:vlbound}, this lower bound on $\ln Z$
does not give a lower bound on the likelihood of the visible units
$\qone(\vu)$. Nevertheless, we may hope that the approximation to the
likelihood gradient is sufficiently accurate such that ascent of the
likelihood can be made.
%
Taking the derivative of \eqref{eq:vlbound} with respect to the
parameters of the Boltzmann machine $\qone(\hu,\vu)$, we arrive at the following
learning rule for gradient ascent given  a pattern $\vu$:
%
\begin{equation}
\Delta J = \eta\br{\av{\s_i\s_j}_{\qtwo(\hu|\vu)} - \av{\s_i\s_j}_{\qone(\hu,\vu)}}
\label{eq:deltaJ}
\end{equation}
%
where $\eta$ is a chosen learning rate.  The correlations in the first
term of \eqref{eq:deltaJ} are straightforward to compute in the case
of using a factorised distribution $\qtwo$. The ``free'' expectations
in the second term are more troublesome and need to be approximated in
some manner. We examine using the standard factorised model
approximations for the free correlations and more accurate
``higher-order'' approximations as described in
\secref{sec:hocorr}. Here, we are primarily concerned with monitoring
the likelihood bound \eqref{eq:vlbound} under such gradient dynamics,
so we will look at the exact bound value, it's approximation using a
factorised model for $\ln Z$ \eqref{eq: bound fact}, and it's
second-order correction \eqref{eq: MI second order}.


%\subsection{Results}

We consider a small learning problem with 3 hidden units and 4 visible
units. Ten visible training patterns were formed in which the
elements were chosen to be on (value=1) with probability 0.4. In
figure 7(a) we demonstrate variational learning in
which the free statistics, required for the likelihood gradient ascent
rule \eqref{eq:deltaJ}, are approximated with a simple factorised
assumption, $\av{s_i s_j}_{\qone(\hu,\vu)}\approx m_i m_j (i\neq j)$.
%
We display the time evolution of the exact training pattern likelihood
bound\footnote{In practice, the likelihood on a test set is more
appropriate. However, the small number of possible patterns here (16)
makes it rather difficult to form a representative test set.}
\eqref{eq:vlbound} (solid line), and two approximations to it: the
first approximates the intractable $\ln Z$ term by using the standard
variational approach with a factorised model (dashed line). The second
approximation uses the second-order correction to the factorised model
value for $\ln Z$ (dot-dash line).
%
The learning rate was
fixed at 0.05.
%As can be seen, the true
%likelihood bound displays an instability after around 40 gradient
%updates, and this can be attributed to the emergence of
%non-unimodality in the learned distribution $\qone(\hu,\vu|J)$.
%However, i
In monitoring the progress of learning, the higher-order correction to
the likelihood bound can be seen to be a more accurate approximation
of the likelihood compared to that given by the use of the standard
variational approximation alone. Interestingly, using the second-order
correction to the likelihood, a maximum in the likelihood is detected
at a point close to saturation of the exact bound. This is not the case
using the first-order approximation alone and, with the standard
approximation to the likelihood, the dynamics of the learning process,
in terms of the training set likelihood, are poorly monitored.

In figure 7(b) the gradient dynamics are provided again by equation
\eqref{eq:deltaJ} but now with the free correlations $\av{s_i
s_j}_{\qone(\hu,\vu)}$ approximated by the more accurate
ratio-of-normalising-constants method, as described in
\secref{sec:hocorr}. Again, we see an improvement in the accuracy of
monitoring the progress of the learning dynamics by the simple
inclusion of the higher-order term and, again, the higher-order
approximation displays a maximum at roughly the correct position. The
likelihood is higher than in figure 7(a) since the gradient of the
likelihood bound is more accurately approximated through the more
accurate correlation estimates.  Note that the reason that the exact
likelihood lower bound can be above 1 is due to pattern repetitions in
the training set. The instabilities in learning after approximately
120 updates are due to the emergence of multimodality in the learned
distribution $\qone(\hu,\vu|J)$, indicating that a more powerful
approximation method, able to capture such multimodal effects, should
be used to obtain further improvements in the likelihood.


\begin{figure}[tbp]
\begin{center}
\subfigure[Learning using dynamics in which the free correlations are approximated using a factorised model.]
{\epsfig{figure=results_mf1_final.eps,height=0.3\linewidth}}
\subfigure[Learning using dynamics in which the free correlations are given by the ratio of normalising constants approach.]
{\epsfig{figure=results_mf2.eps,height=0.3\linewidth}}
\end{center}
\label{fig:results_mf3}
\caption{Learning a set of 10 random, 4-visible-unit patterns with a
Boltzmann machine with 3 hidden units. The solid line is the exact
value for the total likelihood bound \eqref{eq:vlbound} of all the
patterns in the training set. The dashed line is the likelihood bound
approximation based on using a factorised model to approximate $\ln
Z$, equation \eqref{eq: bound fact}. The dash-dot line uses the second
order correction to $\ln Z$, equation \eqref{eq: MI second order}.}
\end{figure}




\section{Discussion}
\label{sec: Discussion}

In this article we have described perturbational approximations of
intractable probability distributions.  The approximations are based
on a Taylor series using a family of probability distributions that
interpolate between the intractable distribution and a tractable
approximation thereto.  These approximations can be seen as an
extension of the variational method, although they no longer bound the
quantities of interest.  We have illustrated our approach both
theoretically and by computer simulations for Boltzmann machines.
These simulations showed that the approximation can be improved beyond
the variational bound by including higher-order terms of the
corresponding cumulant expansion.
%
\cut{ Surprisingly, we found a consistent improvement in accuracy
in every single simulation run, although a decrease in accuracy cannot
theoretically be excluded.  Furthermore, the computer simulations
showed that learning also becomes more stable when higher-order terms
are included.  Interestingly, the mode hopping behavior previously
observed in learning experiments~\cite{Lawrence98} was eliminated.
}
%

Simulations showed that the accuracy in monitoring the training set
likelihood during the learning process can be improved by including
higher order corrections. However, these perturbational approximations
cannot be expected to consistently improve on zeroth order
(variational) solutions. For instance, if the distribution is strongly
multimodal, then using a unimodal variational distribution cannot be
expected to be improved much by the inclusion of higher-order
perturbational terms. On the other hand, higher-order corrections to
multimodal approximations may well improve the solution considerably.

We elucidated the relationship of our approach to TAP, arguing that
our approach is expected to offer a more stable solution. Furthermore,
the application of our approach to models other than the Boltzmann
machine is transparent.  Indeed, these techniques are readily
applicable to other variational methods in a variety of contexts both
for discrete and continuous systems~\cite{Barber98b,Wiegerinck98}.

One drawback of the perturbational approach that we have described is
that known ``physical'' constraints, for example that moments must be
bounded in a certain range, are not necessarily adhered to. It would
be useful to develop a perturbation method that ensures that solutions
are at least physical.

We hope to apply these methods to variational techniques other than
the standard Kullback-Leibler bound, and also to study the evaluation
of other suitable criteria for the setting of the variational
parameters.

\section*{Acknowledgments}

We would like to thank Tom Heskes, Bert Kappen, Martijn Leisink, Peter
Sollich and Wim Wiegerinck for stimulating and helpful discussions and
the referees for excellent comments on an earlier version of this
paper.

\appendix
\section{}
\label{sec: appendix}

In this appendix we describe how to optimize the variational
bound when the normalization constant of an intractable
Boltzmann machine is approximated using another tractable
Boltzmann machine.  For convenience, we denote the potential
of an intractable Boltzmann Machine as \begin{equation} H_{1}
\equiv \sum_{I} \weight_{I} \spin_{I} \thickspace \text{,}
\end{equation} where the ``extended'' parameters are given by
$\newweight_{I}\in \{\newweight_{i}, \newweight_{ij} | i \in
[1,\ldots,N], j = [i+1,\ldots,N]\}$ and $\spin_{I}\in
\{\spin_{i}, \spin_{i} \spin_{j} | i \in [1,\ldots,N], j =
[i+1,\ldots,N]\}$.  The potential of the other tractable
Boltzmann Machine is denoted as \begin{equation} H_{0} \equiv
\sum_{I \in \Psi} \newweight_{I} \spin_{I} \label{eq:
tractable} \thickspace \text{,} \end{equation} where $\Psi$
denotes the set of parameters of the tractable Boltzmann
Machine.

We want to optimize the bound of \eqn{\ref{eq: bound}} with respect to
the free parameters $\newweight_{J}$ ($J \in \Psi$), \ie
%
\begin{equation} \frac{\partial}{\partial \newweight_{J}}
\left [\log \normalconst_{0} + {\left < \Delta H \right>}_{0} \right ]
= 0 \thickspace \text{,} \end{equation}
%
which leads to
%
\begin{equation}
{\left < \Delta H \spin_{J} \right >}_{0} - {\left < \Delta H \right
>}_{0} {\left < \spin_{J} \right >}_{0} = 0
\end{equation}
%
and thus to the fixed point equation\footnote{ This
fixed point equation is the generalized mean field equation.  This
form of iteration is not restricted to Boltzmann machines but is a
general characteristic of exponential models.}
%See, for example,
%\citeA{Laar98} in the context of multilayer perceptrons.  }
%
\begin{equation}
\label{eq: MF}
\newweight_{I} = \sum_{J \in \Psi}
\sum_{K} [{F_{\Psi \Psi}}^{-1}]_{IJ} F_{JK} \weight_{K} \thickspace
\text{,}
\end{equation}
%
where the Fisher matrix is given by $F_{IJ} = {\left < \spin_{I}
\spin_{J} \right >}_{0} - {\left < \spin_{I} \right >}_{0} {\left <
\spin_{J} \right >}_{0}$ which depends on the free parameters
$\myvector{\newweight}$, and $F_{\Psi \Psi}$ is the Fisher matrix of
the tractable Boltzmann machine.  For a factorised approximating
model, $F_{\Psi \Psi}$ is a diagonal matrix and \eqn{\ref{eq: MF}}
simplifies to \eqn{\ref{eq: factorised mf equation}}.  For a
decimatable Boltzmann machine, the elements of the Fisher matrix are
tractable.  The iteration can be performed either synchronously or
asynchronously~\cite{Peterson87}.  We prefer here the synchronous case
since for all $N$ parameter-updates we only need to calculate and
invert a single Fisher matrix instead of $N$ matrices in the
asynchronous case. In most applications, however, the asynchronous
method seems to be advantageous with respect to
convergence~\cite{Peterson87,Ansari95,Saul96}.


	\bibliographystyle{theapa}
	\bibliography{/home/snn/davidb/tex/bib/authors,/home/snn/davidb/tex/bib/myfile,/home/snn/davidb/tex/bib/bishop}

\end{document}
