\documentclass[11pt]{article}
\usepackage{epsfig}
\usepackage{theapa}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{jair}
\newcommand{\bra}[1]{\left(#1\right)}
\newcommand{\pot}[1]{\psi_{#1}}
\newcommand{\nodes}[1]{S_{#1}}
\newcommand{\nodescap}[1]{S_\mathrm{\scriptscriptstyle #1}}
\newcommand{\nodescapprime}[1]{S_\mathrm{\scriptscriptstyle #1}'}
\newcommand{\statecap}[1]{\vec s_\mathrm{\scriptscriptstyle #1}}
\newcommand{\statecapprime}[1]{{\vec s\makebox[0pt][l]{$\,'$}}_\mathrm{\scriptscriptstyle #1}}
\newcommand{\ncap}[1]{n_\mathrm{\scriptscriptstyle #1}}
\newcommand{\marnodes}{\nodescap{mar}}
\newcommand{\marnodesprime}{\nodescapprime{mar}}
\newcommand{\sepnodes}{\nodescap{sep}}
\newcommand{\othnodes}{\nodescap{oth}}
\newcommand{\encnodes}{\nodescap{enc}}
\newcommand{\clnodes}{\nodescap{cl}}
\newcommand{\marstate}{\statecap{mar}}
\newcommand{\marstateprime}{\statecapprime{mar}}
\newcommand{\sepstate}{\statecap{sep}}
\newcommand{\othstate}{\statecap{oth}}
\newcommand{\nenc}{\ncap{enc}}
\newcommand{\nsep}{\ncap{sep}}
\newcommand{\prob}[1]{p\bra{#1}}
\newcommand{\qrob}[1]{q\bra{#1}}
\newcommand{\ubound}[1]{p_+\bra{#1}}
\newcommand{\lbound}[1]{p_-\bra{#1}}
\newcommand{\norm}[1]{\left\|#1\right\|}
\newcommand{\trinorm}[1]{\left\{#1\right\}}
\newcommand{\hide}[1]{}
\newcommand{\equationref}[1]{Equation~\ref{#1}}
\newcommand{\sectionref}[1]{Section~\ref{#1}}
\newcommand{\figureref}[1]{Figure~\ref{#1}}

\ShortHeadings{Bound Propagation}{Leisink \& Kappen}
\jairheading{19}{2003}{139-154}{9/02}{8/03}
\firstpageno{139}

\begin{document}
\title{Bound Propagation}
\author{\name Martijn Leisink \email martijn@mbfys.kun.nl \\
\name Bert Kappen \email bert@mbfys.kun.nl \\
\addr University of Nijmegen, Department of Biophysics \\
Geert Grooteplein 21, 6525EZ Nijmegen, The Netherlands }

\maketitle

\begin{abstract}
In this article we present an algorithm to compute bounds on the
marginals of a graphical model.  For several small clusters of nodes upper
and lower bounds on the marginal values are computed independently of
the rest of the network.  The range of allowed probability distributions
over the surrounding nodes is restricted using earlier computed bounds.
As we will show, this can be considered as a set of constraints in
a linear programming problem of which the objective function is the
marginal probability of the center nodes.  In this way knowledge about
the maginals of neighbouring clusters is passed to other clusters thereby
tightening the bounds on their marginals.  We show that sharp bounds
can be obtained for undirected and directed graphs that are used for
practical applications, but for which exact computations are
infeasible.
\end{abstract}

\section{Introduction}

Graphical models have become a popular and powerful tool to deal with
dependencies in a complex probability model.  For small networks an
exact computation of the marginal probabilities for certain nodes is
feasible \shortcite<see, for example>{zhang94a}.  Unfortunately,
due to an exponential scaling of the computational complexity, these
exact algorithms soon fail for somewhat more complex and more realistic
networks.

To deal with this computational problem two main streams can be
distinguished.  Firstly, there is the approach of approximating
the exact solution.  Numerous algorithms have been developed
among which are mean field \shortcite{peterson87a,saul96a} and \textsc{tap}
\shortcite{thouless76a,plefka81a}.  Nowadays, approximation methods such as
belief propagation \shortcite{murphy99a,yedidia00a} using the Bethe or Kikuchi
energy are gaining popularity.  In contrast with mean field expansions,
the latter methods try to model only the local probability distributions
of several small clusters, which communicate with each other in a
way known as message passing.

The second approach is that of finding upper and lower bounds on important
quantities of the graphical model.  In contrast to approximation
techniques, these methods give reliable information.  One may find,
for instance, that a marginal probability is definitely between 0.4
and 0.5.  These methods usually focus on the partition function of
the network, which can be lower bounded \shortcite{mync} or upper bounded
\shortcite{wainwright02a}.  Little work has been done, however, on bounding
the marginal probabities for certain nodes in the network directly.
Up to now researchers have focused their attention on combining upper
and lower bounds on the partition function in which way bounds on the
marginals can be derived.  For the special case of the Boltzmann machine
this was done by \shortciteA{mynips2}.

In this article we present a method called `bound propagation' which
has certain similarities to belief propagation in the sense that
small clusters of a few nodes try to bound their marginals as tight
as possible.  The knowledge contained in these bounds is shared with
surrounding clusters.  Iteration of this procedure improves bounds on all
marginals until it is converged.  It is important to realize that bounding
and approximation techniques are not mutually exclusive.  We propose to
always run the bound propagation algorithm first.  When the result is not
satisfactory, one can always fall back upon approximation techniques as
belief propagation paying the price of being uncertain about the solution.

In \sectionref{theory} we explain the general method of bound propagation,
which is described more algorithmicly in \sectionref{algorithm}.
In \sectionref{toysec} we show the behaviour of the algorithm on a toy
problem, whereas in \sectionref{results} the results are shown for more
realistic networks.  Finally, in \sectionref{discussion}, we conclude
and discuss possible paths for future research.

\section{What Bounds Can Learn From Each Other}\label{theory}

Consider a Markov network defined by the potentials
$\pot{i}\bra{\nodes{i}}$ each dependent on a set of nodes $\nodes{i}$.
The probability to find the network in a certain state $S$ is
proportional to
\begin{equation}
\prob{S}\propto\prod_i\pot{i}\bra{\nodes{i}}
\end{equation}
Problems arise when one is interested in the marginal probability over a
small number of nodes (denoted by $\marnodes$), since in general this
requires the summation over all of the exponentially many states.

\begin{figure}
{\centering\begin{tabular}{cccc}
\epsfig{file=example_single.eps,width=2.1cm}&
\epsfig{file=example_vertical.eps,width=2.1cm}&
\epsfig{file=example_diagonal.eps,width=2.8cm}&
\epsfig{file=example_vertical_other.eps,width=2.1cm}\\
a&b&c&d
\end{tabular}\par
\caption{Examples of separator and marginal nodes in an (infinite)
Ising grid.  Light shaded nodes correspond to $\marnodes$, black nodes
to $\sepnodes$.  (a) A single node enclosed by its tightest separator.
(b), (c) Two marginal nodes enclosed by their tightest separator. (d)
Another choice for a separator for a single node marginal; the other
node enclosed by that separator belongs to $\othnodes$.}\par
\label{ising_example}}
\end{figure}

Let us define the set of separator nodes, $\sepnodes$, to be the nodes
that separate $\marnodes$ from the rest of the network.  One can take,
for instance, the union of the nodes in all potentials that contain at
least one marginal node, while excluding the marginal nodes itself in that
union (which is the Markov blanket of $\marnodes$):
\begin{equation}
\sepnodes=\bigcup_{\nodes{i}\cap\marnodes\neq\emptyset}
	\!\!\!\nodes{i}/\marnodes
\end{equation}
See for instance \figureref{ising_example}a--c.  It is not necessary,
however, that $\sepnodes$ defines a tight separation.  There can
be a small number of other nodes, $\othnodes$, inside the separated
area, but not included by $\marnodes$ (\figureref{ising_example}d).
A sufficient condition for the rest of the theory is that the size of
the total state space, $\norm{\sepnodes\cup\marnodes\cup\othnodes}$,
is small enough to do exact computations.

Suppose that we know a particular setting of the separator nodes,
$\sepstate$, then computing the conditional distribution of
$\marnodes$ given $\sepnodes$ is easy:
\begin{equation}
\prob{\marstate\,|\,\sepstate}=\sum_{\othstate\in\othnodes}
	\!\!\prob{\marstate\othstate\,|\,\sepstate}
\end{equation}
This expression is tractable to compute since $\sepnodes$ serves as the
Markov blanket for the union of $\marnodes$ and $\othnodes$.  Thus the actual numbers
are independent of the rest of the network and can be found either by
explicitly doing the summation or by any other more sophisticated method.
Similarly, if we would know the exact distribution over the separator nodes,
we can easily compute the marginal probability of $\marnodes$:
\begin{equation}
\prob{\marstate}=\sum_{\sepstate\in\sepnodes}\sum_{\othstate\in\othnodes}
	\!\!\prob{\marstate\othstate\,|\,\sepstate}\prob{\sepstate}\label{marginal}
\end{equation}

Unfortunately, in general we do not know the distribution $\prob{\sepstate}$
and therefore we cannot compute $\prob{\marstate}$.  We can, however,
still derive certain properties of the marginal distribution, namely
upper and lower bounds for each state of $\marnodes$.  For instance,
an upper bound on $\marstate$ can easily be found by computing
\begin{equation}
\ubound{\marstate}=\mathrm{max}_{\qrob{\sepstate}}\sum_{\sepstate\in\sepnodes}
	\sum_{\othstate\in\othnodes}
	\!\!\prob{\marstate\othstate\,|\,\sepstate}\qrob{\sepstate}\label{linprog}
\end{equation}
under the constraints
$\forall_{\sepstate\in\sepnodes}\;\qrob{\sepstate}\ge0$ and
$\sum_{\sepstate\in\sepnodes}\qrob{\sepstate}=1$.  The distribution
$\qrob{\sepstate}$ corresponds to $\norm{\sepnodes}$ free parameters
(under the given constraints) used to compute the worst case (here:
the maximum) marginal value for each $\marstate$.  Obviously, if
$\qrob{\sepstate}=\prob{\sepstate}$ the upper bound corresponds to the
exact value, but in general $\prob{\sepstate}$ is unknown and we have to
maximize over all possible distributions.
One may recognize \equationref{linprog} as a simple linear programming
(\textsc{lp}) problem with $\norm{\sepnodes}$ variables.  With a similar
equation we can find lower bounds $\lbound{\marstate}$.

Upto now we did not include any information about $\qrob{\sepstate}$
we might have into the equation above.  But suppose that we
already computed upper and lower bounds on the marginals of other
nodes than $\marnodes$.  How can we make use of this knowledge?
We investigate our table of all sets of marginal nodes we bounded
earlier and take out those that are subsets of $\sepnodes$.  Obviously,
whenever we find already computed bounds on $\marnodesprime$ where
$\marnodesprime\subseteq\sepnodes$\footnote{In the rest of this article
we will use the symbol $\marnodesprime$ to indicate an earlier computed
bound over $\marnodesprime$ with $\marnodesprime\subseteq\sepnodes$.},
we can be sure that
\begin{equation}
\sum_{\sepstate\in\sepnodes/\marnodesprime}\!\!\!\qrob{\sepstate}=
\qrob{\marstateprime}\le\ubound{\marstateprime}\label{aconstraint}
\end{equation}
and similarly for the lower bound.  These equations are nothing more
than extra contraints on $\qrob{\sepstate}$ which we may include in the
\textsc{lp}-problem in \equationref{linprog}.
There can be several $\marnodesprime$ that are subsets of $\sepnodes$ each
defining $2\norm{\marnodesprime}$ extra constraints.

In this way, information collected by bounds computed earlier, can
propagate to $\marnodes$, thereby finding tighter bounds for these
marginals.  This process can be repeated for all sets of marginal
nodes we have selected until convergence is reached.

\subsection{Linear Programming}

\begin{figure}
{\centering\begin{tabular}{cc}\epsfig{file=tri.eps,width=70mm}&
\epsfig{file=cut.eps,width=70mm}\\a&b\end{tabular}\par
\caption{The area in which $\qrob{s_1s_2}$ must lie is shown.
$\qrob{00}$ is implicitly given since the distribution is normalized.
a) The pyramid is the allowable space.  The darker planes show how
this pyramid can be restricted by adding earlier computed bounds as
constraints in the linear programming problem.  b) This results in a
smaller polyhedron.  The black lines show the planes where the
function $\sum_{s_1s_2}\prob{s_0\!=\!1\,|\,s_1s_2}\qrob{s_1s_2}$ is
constant for this particular problem.}\par
\label{example3d}}
\end{figure}

In terms of standard linear programming, \equationref{linprog} can be
expressed as
\begin{equation}
\mathrm{max}\,\Big\{\vec c\cdot\vec x\:\Big\vert\: \vec x\ge0
	\,\wedge\, A\vec x\le\vec b\,\Big\}
\end{equation}
where the variables are defined as
\begin{align}
\vec x&=\qrob{\sepstate}\label{x}\\
\vec c&=\sum_{\othstate\in\othnodes}\!\!\prob{\marstate\othstate\,|\,\sepstate}
=\prob{\marstate\,|\,\sepstate}
\label{c}\\
A&=\left\{\begin{array}{l}
		\delta\bra{\marstateprime,\sepstate}\\
		-\delta\bra{\marstateprime,\sepstate}
	\end{array}\right.\quad
\vec b=\left\{\begin{array}{l}
		\ubound{\marstateprime}\\
		-\lbound{\marstateprime}
	\end{array}\right.
\qquad\bra{\forall_{\marnodesprime\subseteq\sepnodes}
	\forall_{\marstateprime\in\marnodesprime}}\label{Ab}
\end{align}
where $\delta\bra{\marstateprime,\sepstate}=1$ iff the states of the
nodes both node sets have in common are equal.  Thus the columns of
the matrix $A$ correspond to $\sepstate$, the rows of $A$ and $b$ to
all the constraints (of which we have $2\norm{\marnodesprime}$ for each
$\marnodesprime\subseteq\sepnodes$).  The ordering of the rows of $A$
and $b$ is irrelevant as long as it is the same for both.  The constraint
that $\qrob{\sepstate}$ should be normalized can be incorporated in $A$
and $\vec b$ by requiring
\begin{equation}
\sum_{\sepstate}\qrob{\sepstate}\le1\qquad\textrm{and}\qquad
-\sum_{\sepstate}\qrob{\sepstate}\le-1\label{normalization}
\end{equation}
The maximum of $\vec c\cdot\vec x$ corresponds to $\ubound{\marstate}$.
The negative of $\lbound{\marstate}$ can be found by using $-\vec
c\cdot\vec x$ as the objective function.

\subsection{A Simple Example}

Imagine a single neuron, $s_0$, that is separated from the rest of
the network by two other neurons, $s_1$ and $s_2$.  We want to compute
bounds on the marginal $\prob{s_0\!=\!1}$.  Since we do not know the
distribution over the separator nodes, we have to consider all possible
distributions as in \equationref{linprog}.  Therefore we introduce
the free parameters $\qrob{s_1s_2}$.  The goal is now to compute the
minimum and maximum of the function
$\sum_{s1s2}\prob{s_0\!=\!1\,|\,s_1s_2}\qrob{s_1s_2}$.
In \figureref{example3d}a we show in three dimensions the space
(the large pyramid) in which the distribution $\qrob{s_1s_2}$ lies.
$\qrob{00}$ is implicitly given by one minus the three other values.

We can add, however, the earlier computed (single node)
bounds on $\prob{s_1}$ and $\prob{s_2}$ to the problem.
These restrict the space in \figureref{example3d}a further,
since for instance (see also \equationref{aconstraint})
\begin{equation}
\qrob{s_1\!=\!1}=\sum_{s_2}\qrob{s_1\!=\!1,s_2}
=\qrob{10}+\qrob{11}\le\ubound{s_1\!=\!1}
\end{equation}
We have four independent constraints, which are shown in
\figureref{example3d}a as planes in the pyramid.

Obviously, by adding this information the space in which $\qrob{s_1s_2}$
may lie is restricted to that shown in \figureref{example3d}b.
In the same figure we have added black lines, which correspond to
the planes where the objective function
is constant.  A standard linear programming tool will immediately return
the maximum and the minimum of this function thus bounding the marginal
$\prob{s_0\!=\!1}$.

\section{The Algorithm In Detail}\label{algorithm}

To solve an arbitrary network we first determine the set $\Omega$ of all
$\sepnodes$ we will use for the problem.  We choose
\begin{equation}
\Omega\bra{N}=\left\{\sepnodes\;\vert\;\norm{\clnodes}\le N\,
	\wedge\,\forall_{s\notin\clnodes}\norm{\clnodes\cup s}>N\right\}
\label{omega}
\end{equation}
where $\clnodes$ is the cluster $\sepnodes\cup\marnodes\cup\othnodes$
and $s$ denotes a single node.
This choice limits the state space of $\clnodes$ to $N$, which more or
less determines the computation time.  If we choose, for example, $N=2^5$
for the Ising grid in \figureref{ising_example} assuming binary nodes,
only the configuration in \figureref{ising_example}a would be included
in $\Omega$.  Choosing $N=2^8$ puts the configurations
in \figureref{ising_example}b--d and a few others not shown into
$\Omega$, but not \figureref{ising_example}a, since that separator is
already embedded in larger ones (e.g.\ \figureref{ising_example}d).
In other words:  $\Omega\bra{N}$ is the set of all $\sepnodes$ for
which $\norm{\sepnodes\cup\marnodes\cup\othnodes}$ is not larger than
$N$, but if we add a single node it would cross that boundary.

We will compute bounds for all $\marnodes$ for which there can be found
a separator in $\Omega$.  We reserve memory for $\ubound{\marnodes}$
and $\lbound{\marnodes}$ and initialize them to one and zero
respectively.  Note that if we have a bound over $\marnodes$ it is
still worthwhile to compute bounds over subsets of $\marnodes$ (e.g.\
\figureref{ising_example}b and d).  In contrast to a joint probability
table, the bounds on marginals over subsets cannot be computed simply
by summing over the marginal of $\marnodes$, since we have only bounds
and not the real values.

\subsection{The Iterative Part}

For all $\sepnodes\in\Omega$ we perform the same action.  First we set
up the $A$ and $\vec b$ for the \textsc{lp}-problem, since this
depends only on $\sepnodes$ and the already computed bounds for which
$\marnodesprime\subseteq\sepnodes$.  The number of variables in our
\textsc{lp}-problem obviously is $\norm{\sepnodes}$.  The number of
(inequality) constraints is equal to
\begin{equation}
\sum_{\marnodesprime\subseteq\sepnodes}2\norm{\marnodesprime}
\label{numberofconstraints}
\end{equation}
Once the \textsc{lp}-problem is set up this far, we try to improve all
bounds that are defined over an $\marnodes$ for which $\sepnodes$ is
its separator.  The separator in figures \ref{ising_example}b and d,
for instance, is identical, but the $\marnodes$'s differ.  For each bound we
iterate over its state space and set the $\vec c$ as in \equationref{c}
to its appropriate value.  Then we compute the new upper and lower
bounds for that $\marstate$ by solving the \textsc{lp}-problem twice:
maximize and minimize.  If the new found value improves the bound,
we store it, otherwise we abandon it.

The iterative procedure is repeated until convergence is reached.  In
our simulations we define a bound as being converged as soon as the
relative improvement after one iteration is less than one percent.  If
this holds for all bounds, the algorithms stops.

\subsection{Computational Complexity}

\textsc{lp}-problems have been thoroughly studied in the literature.
The computational complexity of such problems is shown to be
polynomial \shortcite{khachiyan}.  But more importantly, for most
practical problems the number of iterations needed is close to linear
in the number of variables or the number of constraints, whichever
is less \shortcite{todd}.  In each iteration a `pivoting action'
needs to be performed, which is some operation that roughly accesses
all the elements of the matrix $A$ once.  Therefore the expected
computational complexity for solving one \textsc{lp}-problem is ${\cal
O}\bra{m^2\norm{\sepnodes}}$, where $m$ is the number of constraints.
The actual observed computation time depends heavily on the particular
problem, but in general \textsc{lp}-problems upto tens of thousands
of variables and constraints are reasonable.  This makes the presented
method tractable for separators with $\norm{\sepnodes}$ of a similar size.

For every $\marstate\in\marnodes$ we need to solve two
\textsc{lp}-problems:  one for the upper and one for the lower bound.
Note that in this iteration, we do not need to change the matrix $A$
and vector $\vec b$ from \equationref{Ab}.  The vector $\vec c$, that
defines the objective function, obviously does vary with $\marstate$.
Therefore, we expect a total computational complexity of ${\cal
O}\bra{m^2\norm{\marnodes}\norm{\sepnodes}}$ for updating one cluster
$\marnodes$.  How quickly the algorithm converges is more difficult to
estimate.  In the next section, however, we address this topic on the
basis of a toy problem.

\hide{Obviously, the computation time is dominated by the most time comsuming
$\textsc{lp}$-problems.  That is, those $\textsc{lp}$-problems for which
$\norm{\sepnodes}$ is large.  We will now argue that in that case the
computational complexity is mainly determined by $\norm{\sepnodes}$ and
that $m^2$ should be viewed as a more or less constant multiplicative
factor.  Suppose we have an architecture with binary nodes and we try
to solve an $\textsc{lp}$-problem based upon a rather large separator
of 14 nodes.  The number of variables in such a problem is $2^{14}$.
The number of constraints is determined by adding the state spaces of all
$\marnodesprime$ we can find (see \equationref{numberofconstraints}).
First of all, the state space of the $\marnodesprime$'s found will not
be very large, since their own separator should be tractable for the
bound propagation algorithm.  Secondly, due to the fact that the nodes
in $\marnodesprime$ are always close together, we expect to find about
14 $\marnodesprime$'s with a state space $2^k$ instead of $14\choose
k$.  In worst case (speaking about computation complexity) in this
example we would expect to find about}


\hide{Let us focus for a while on the computational complexity of solving
one \textsc{lp}-problem with separator $\sepnodes$.  As mentioned
earlier the number of variables is given by $\norm{\sepnodes}$.  The
number of constraints (apart from the normalizing constraint in
\equationref{normalization}) is given in
\equationref{numberofconstraints}.

For the most time consuming \textsc{lp}-problems (for which
$\norm{\sepnodes}$ is large) it is reasonable to assume that the number
of constraints is much less than $\norm{\sepnodes}$.  This can be seen
as follows.  The only way to have a number of constraints comparable to
$\norm{\sepnodes}$ is that there exist $\marnodesprime$'s which include
most or all of the nodes in $\sepnodes$.  It is, however, very unlikely,
that the separator belonging to such an $\marnodesprime$ is small enough
to be tractable, since $\sepnodes$ itself was already close to the
largest we can handle.  Thus these $\marnodesprime$'s does not exist
in our collection of bounds.  Therefore the computational complexity
of the most time consuming \textsc{lp}-problems is mainly determined
by $\norm{\sepnodes}$.}

To conclude this section, we remark that when the \textsc{lp}-problem is
not tractable, one can always leave out as many constraints as needed,
paying the price of getting looser bounds.\hide{  Thus to summarize:  For the
most time consuming (i.e.\ large) sets of separator nodes we have to solve
$\norm{\marnodes}$ \textsc{lp}-problems for all $\marnodes$ for which
$\sepnodes$ is a separator.  The computational complexity for solving
one $\textsc{lp}$-problem is ${\cal O}\bra{m^2\norm{\sepnodes}}\approx
{\cal O}\bra{\norm{\sepnodes}^2}$.}


\section{A Toy Problem}\label{toysec}

\begin{figure}
{\centering\epsfig{file=ringtest.eps,width=10cm}\par
\caption{A Boltzmann machine ring of 200 nodes is initialized with
thresholds and weights as shown in the lower part of the figure.
Exact means are plotted as well as the results of the bound
propagation algorithm.  The number of nodes included in $\marnodes$ is
indicated in the legend.}\par
\label{ring}}
\end{figure}

To get some intuition about the performance of the algorithm we created
a ring (one large loop) of 200 binary nodes.  Nodes are connected as in
the Boltzmann machine with a weight $w$ and each node has a threshold
$\theta$.  Weights and thresholds are varying along the ring as is shown
at the bottom of \figureref{ring}.  For this simple problem we can compute
the exact single node marginals.  These are plotted as a solid line in
\figureref{ring}.  We ran the bound propagation algorithm
three times, where we varied the maximum
number of nodes included in $\marnodes$: $\left|\marnodes\right|=1$,
$\left|\marnodes\right|=2$ and $\left|\marnodes\right|=3$.  In all cases
we chose two separator nodes:  the two neighbors of $\marnodes$.  Thus these
three cases correspond to $\Omega\bra{2^3}$, $\Omega\bra{2^4}$ and
$\Omega\bra{2^5}$.  As is clear from \figureref{ring} the simplest
choice already find excellent bounds for the majority of the nodes.
For large negative weights, however, this choice is not sufficient.
Allowing larger separators clearly improves the result.  In all cases
the computation time was a few seconds\footnote{We used a computer
with a 1~GHz Athlon processor.}.

\subsection{Fixed Points and Convergence}

\begin{figure}
{\centering\begin{tabular}{cc}
\epsfig{file=bmchaintest.eps,width=7cm}&
\epsfig{file=bmchaintest_convergence.eps,width=7cm}\\
a&b
\end{tabular}\par
\caption{The result for a Boltzmann machine ring of an arbitrary number
of nodes.  All thresholds and weights have the same value.  In the
left panel the gap is shown between the upper and lower bound of the
mean of the nodes.  The algorithm converges exponentially to the final
bounding values as soon as it is close enough.  Thus the distance to its
asymptotic value is proportional to $\alpha^n$, where $n$
is the number of iterations.  The parameter $\alpha$ is shown in the
right panel.}\par
\label{bmchaintest}}
\end{figure}

It is possible to derive the fixed point of the bound propagation
algorithm analytically for a Boltzmann ring if we take a single value
for all weights ($w$) and one for all thresholds ($\theta$).  Due to this
symmetry all single node marginals are equal in such a network.  Moreover
all upper and lower bounds on the single nodes should be identical.
This implies that for the fixed point the following holds for any $i$:
\begin{equation}
\ubound{s_i}=\mathrm{max}_{\qrob{s_{i-1}s_{i+1}}}
	\prob{s_i|s_{i-1}s_{i+1}}\qrob{s_{i-1}s_{i+1}}\label{analytical}
\end{equation}
under the constraints
\begin{align}
\qrob{s_{i-1}}&\le\ubound{s_{i-1}}=\ubound{s_i}\\
\qrob{s_{i+1}}&\le\ubound{s_{i+1}}=\ubound{s_i}\\
\qrob{s_{i-1}}&\ge\lbound{s_{i-1}}=\lbound{s_i}\\
\qrob{s_{i+1}}&\ge\lbound{s_{i+1}}=\lbound{s_i}
\end{align}
and similarly for the lower bound.  The conditional probability in
\equationref{analytical} is given by
\begin{equation}
\prob{s_i|s_{i-1}s_{i+1}}\propto\exp\bra{\theta s_i+ws_{i-1}s_i+ws_is_{i+1}}
\end{equation}
From these equations one can derive the fixed point for $\ubound{s_i}$
and $\lbound{s_i}$.  It turns out, however, that it is easier to
determine the fixed point of the upper and lower bound on the mean, hence
$\ubound{s_i=1}-\lbound{s_i=-1}$ and $\lbound{s_i=1}-\ubound{s_i=-1}$.
But this has no effect on the principle results as tightness and
convergence speed.

In \figureref{bmchaintest}a the difference between the upper and
lower bound on the means is shown for various choices of the weight
and threshold.  As we saw before tight bounds can be obtained for
small weights and somewhat larger, but positive weights, whereas
negative weights result in rather poor, but still non-trivial bounds.
It should be mentioned, however, that these results correspond to
the choice of the smallest clusters ($\left|\marnodes\right|=1$
and $\left|\sepnodes\right|=2$) for the bound propagation algorithm.
The bounds can easily be improved by enlarging the clusters as we saw
in \figureref{ring}.

Close to the fixed point the bound propagation algorithm converges
exponentially.  The distance to the asymptotic value can be written
as $\alpha^n$, where $n$ is the number of iterations and $\alpha$
is a number between zero and one indicating the convergence speed.
The closer to one $\alpha$ is, the slower the algorithm converges.
In \figureref{bmchaintest}b $\alpha$ is shown for the same weights and
thresholds.  It is clear that larger weights induce a slower convergence.
That is what we see in general:  probabilities that are close to
deterministic slow down the convergence.

\section{Real World Networks}\label{results}

To assess the quality of the bound propagation algorithm for more
challenging networks, we computed bounds for a Bayesian network
(\sectionref{alarmsec}), that is commonly used in articles as a
representative of a real world problem.  Secondly, we show that the bound
propagation method can be used to find many tight bounds for a large
Ising grid (\sectionref{isingsec} and finally, in \sectionref{qmrsec}, we
show results for a bi-partite graph.  For these networks, we will briefly
discuss bound propagation in relation with the cluster variation method
\shortcite{kappen02novel} which is one of the better approximation methods
for general graphs.  Whenever possible, we show the exact marginals.

\subsection{The Alarm Network}\label{alarmsec}

\begin{figure}
{\centering\epsfig{file=alarm25000.eps,width=14cm}\par
\caption{The exact marginals (thick horizontal lines) and the bounds (rectangles)
for the alarm network.  Some bounds are so tight that the rectangle is
not visible.
We used $\Omega\bra{25,000}$.}\par
\label{alarm}}
\end{figure}

\begin{figure}[ht]
{\centering\epsfig{file=lattice.eps,width=12cm}\par
\caption{The band width (upper minus lower bound) of each node in a 40x40
Ising grid is shown by the blackness of the squares.  The exact marginal
could not be computed for such a large network.  The solid lines show the
boundary where the band width was 0.1.  We used $\Omega\bra{2500}$.}\par
\label{lattice}}
\end{figure}

The Alarm network\footnote{The data can be
downloaded from the Bayesian Network Repository at\\
\texttt{http://www.cs.huji.ac.il/labs/compbio/Repository/}} is a commonly
used network which is a representitive of a real life Bayesian network.
It was originally described by \shortciteA{beinlich} as a network for
monitoring patients in intensive care.  It consists of 37 nodes of two,
three or four states and is small enough to do exact computations.
In \figureref{alarm} the graphical model is shown.  Our task was
to determine the single node marginals in absence of any evidence.
For each node the exact marginal probabilities for its states are shown as
horizontal lines.  The top and bottom of the rectangles around these
lines indicate the upper and lower bounds respectively.
For this network we used $\Omega\bra{25,000}$\footnote{We started with
$\Omega\bra{100}$ and gradually increased this value up to 25,000.
The bounds obtained with small clusters are used as good starting points
for the more time consuming larger clusters.}.  If we make this
choice, we can treat exactly the minimal Markov blanket of all single
node marginals.  The algorithm converged in about six minutes.

We see that the bounds in the lower left corner are very
close to the exact value.  We can give some intuition why this happens.
Observe that node 7 on its own can play the role of separator for
nodes 1 to 6.  Thus all correlations between these nodes and the rest
of the network are through node 7, which means that their uncertainty
only depends on the uncertainty of node 7.  The latter happens to be
very small (consider that as a given fact) and therefore the marginal
bounds for nodes 1 to 6 are very tight.

The bounds in the upper right corner, on the other hand, are quite
poor.  This is partially due to the density of the connections there and the
number of states each of the nodes has.  In combination this leads to
large state spaces already for small clusters.  Therefore, we cannot
compute bounds on sets of multiple nodes, which generally leads to
weaker bounds.  Secondly, the local probability tables are of a
form that is difficult for the bound propagation algorithm.  That is,
a lot of probabilities are close to zero and one.

This picture is what we generally see.  There are some parts in the
network that are hard to compute (not necessarily impossible) which
usually have a more dense structure and a larger (local) state space.
Other parts can be bounded very well.  Note that this network could easily
be embedded in a very large network that is intractable for any exact
method.  Even in that case, the bound propagation algorithm could be used
to find similar results.  Since the method is based on local computations,
the existance of a large surrounding network has almost no influence.

\subsection{A Large Ising Grid}\label{isingsec}

\hide{\begin{figure}
{\centering\begin{tabular}{c|c}
\epsfig{file=100.eps,width=65mm}&\epsfig{file=2500.eps,width=65mm}\\
a&b\\
\end{tabular}\par
\caption{Each panel is the single node marginal of a 4x4 Ising grid
with nearest neighbour interaction.  The nodes took 2, 3 or 4 states
as can be seen in the figure.  The exact marginal is shown by
the thick horizontal lines.  The top and bottom of the rectangles around
these lines indicate the upper and lower bound found respectively.
(a) $\Omega=\Omega\bra{100}$ converged using 2 seconds of \textsc{cpu}
time, (b) $\Omega=\Omega\bra{2500}$ converged using 36 seconds of
\textsc{cpu} time.}\par
\label{ising}}
\end{figure}}

We created a so called two-dimensional Ising model, which is a
rectangular grid with connections between all pairs of nodes that
are direct neighbours.  We used a grid of 40 by 40 binary nodes.
The potentials were drawn from a uniform distribution between zero
and one.  In contrast with Bayesian belief networks \shortcite{neil92a}
the probability distribution over the unclamped Ising grid is not
automatically normalized, but this has no consequences for the algorithm.

For small Ising grids computing the exact marginals is tractable, but
as soon as the grid is larger than about 25 by 25 the exact algorithm
fails due to the exponential scaling of the computational complexity.
The bound propagation algorithm, on the other hand, only depends on the
local structure (i.e.\ the size of the Markov blankets of each node) and
thus scales linearly with the number of nodes.  We created an 40x40 Ising
grid with binary nodes similarly as above and ran the bound propagation
algorithm.  For this network the exact algorithm would require a storage
capacity of at least $2^{41}$ real numbers, whereas bound propagation
converged in about 71 minutes in which time it computed bounds on the
marginals over all 1600 nodes.

In \figureref{lattice} we show the 40x40 grid where the blackness of
the squares correspond to the band width of the single node marginals.
This band width is defined as the difference between the upper and
lower bound.  Due to the fact that marginal probabilities sum to one,
the two band widths for these binary neurons are identical.  We can
clearly see some regions in the lattice (the blacker area's) for which
bound propagation had some difficulties.  Most of the marginals, however,
are bounded quite well:  more than 75\% had a band width less than 0.1.

\subsection{The Bi-Partite Graph}\label{qmrsec}

\begin{figure}
{\centering\epsfig{file=qmr.eps,width=14cm,height=4cm}\par
\caption{A bi-partite network of twenty nodes in each layer.  Each child
node (bottom layer) has three parents in the top layer and each
parent node has exactly three children.  Probability tables were randomly
initialized from a uniform distribution.  The shaded nodes in the bottom
layer are clamped to a certain state.  The bound propagation algorithm
ran with $\Omega\bra{2500}$ and converged in 37 seconds.  For each node
the thick horizontal lines show their exact marginals, the top and bottom
of the rectangle are the upper and lower bounds found.}\par
\label{qmr}}
\end{figure}

A bi-partite graph is a network consisting of two layers of nodes,
where the top layer is hidden and the bottom layer visible.  The only
connections in the network are from the top layer (parent nodes)
to the bottom layer (child nodes).  Such a architecture appears
very simple, but already with several tens of nodes it is often
intractable to compute the marginals exactly when evidence is presented.
A bi-partite graph can be undirected or directed (pointing downwards).
A recent example of the first is the \textsc{p}o\textsc{e} (product of
experts) from \shortciteA{hinton99a} with discrete nodes.  The `Quick
Medical Reference' (\textsc{qmr}) network from \shortciteA{shwe91a}
is a good example for a directed bi-partite graph.  Recently there
were some approaches to use approximation techniques for this network
\shortcite{jaakkola99variational,murphy99a}.  Therefore we have chosen to use
the bi-partite network as an architecture to test the bound propagation
algorithm on.

We created bi-partite networks of an equal number of (binary valued)
nodes in both layers.  Each child node has exactly three parents and
each parent points to exactly three children.  The connections were
made randomly.  The conditional probability tables (\textsc{cpt}'s)
holding the probability of the child given its parents was initialized
with uniform random number between zero and one.  The prior probability
tables for the parents were initialized similarly.

In \figureref{qmr} we show a bi-partite network with twenty nodes in
each layer, which is still small enough to compute the exact marginals.
The arrows show the parent-child relations.  To every third child evidence
was presented (the shaded nodes).  We ran the bound propagation algorithm
with $\Omega\bra{2500}$ and reached convergence in 37 seconds.  For every
node we show the marginal over its two states.  The thick line is the
exact marginal, the top and bottom of the rectangle around it indicate
the upper and lower bound respectively.  It is clear that for most of
the nodes the computed bounds are so tight that they give almost the
exact value.  For every node (except two) the gap between upper and
lower bound is less than one third of the whole region.

\begin{figure}
{\centering\begin{tabular}{cc}
\epsfig{file=delta2.eps,width=6cm}&
\epsfig{file=cvm_vs_bound.eps,width=6cm}\\
a&b\\
\end{tabular}\par
\caption{
A bi-partite network of a thousand nodes in each layer was created and
evidence was presented on every third child node.
The left histogram shows the band widths found for the single node
marginals of all unclamped nodes.  The bound propagation
algorithm converged in about 75 minutes and used $\Omega\bra{2500}$.
The right histogram shows the relative position of the approximation
within the bounding region.  All approximations were between the upper
and lower bound, although this is not necessarily the case.}\par
\label{qmr1000}}
\end{figure}

The network shown in \figureref{qmr} is small enough to be treated
exactly.  This enables us to show the computed bounds together with the
exact marginals.  The bound propagation algorithm, however, can be used
for much larger architectures.  Therefore we created a bi-partite graph
as we did before, but this time with a thousand nodes in each layer.
Again we presented evidence to every third node in the bottom layer.
The bound propagation algorithm using $\Omega\bra{2500}$ converged in
about 75 minutes.  In \figureref{qmr1000}a we show a histograms of
the band widths (upper minus lower bound) found for the single node
marginals for all unclamped nodes.
Clearly, for the majority of the marginals very tight bounds are found.

Although exact computations are infeasible for this network, one can
still use approximation methods to find an estimate of the single node
marginals.  Therefore, we ran the cluster variation method (\textsc{cvm})
as described in \shortcite{kappen02novel}.  For the clusters needed
by \textsc{cvm}, we simply took all node sets for which a potential
is defined.  For every unclamped node we computed the relative position
of the approximation for the single node marginal within the bounding
interval, denoted by $\eta$.  Thus $\eta=0$ indicates an approximation
equal to the lower bound, and similarly $\eta=1$ for the upper bound.
Although there is no obvious reason the approximated marginals should
be inside the bounding interval, it turned out they were.  A histogram
of all the relative positions is shown in \figureref{qmr1000}b, where
we split up the full histogram into three parts, each referring to
a certain region of band widths.  It is remarkable that the shape of
the histogram does not vary much given the tightness of the bounds.
Approximations tend to be in the middle, but for all band widths there
are approximations close to the edge.

One can argue that we are comparing apples to pears, since we can
easily improve the results of both algorithms, while it is not clear
which cases are comparable.  This is certainly true, but it is not our
intention to make a comparison to determine which method is the best.
Approximation and bounding methods both have their own benefits.
We presented \figureref{qmr1000}b to give at least some intuition about
the actual numbers.  In general we found that approximations are usually
within the bounding intervals, when computation time is kept roughly
the same.  This does, however, not make the bounds irrelevant.  On the
contrary, one could use the bounding values as a validation whether
approximations are good or even use them as confidence intervals.

\hide{\begin{figure}
{\centering\begin{tabular}{cc}
\epsfig{file=single.eps,width=60mm}&\epsfig{file=multiple.eps,width=60mm}\\
a&b\\
\end{tabular}\par
\caption{The two histograms show all the band widths (upper minus
lower bound) found for (a) the single node and (b) the multiple node
marginals for a bi-partite network of a thousand nodes in each layer.
Evidence was presented on every third child node.  The bound propagation
algorithm converged in about 75 minutes and used $\Omega\bra{2500}$.
For this choice for $\Omega$ and this particular instance of the
network there were 1666 bounds with $\left|\marnodes\right|=1$, 2022
with $\left|\marnodes\right|=2$, 1343 with $\left|\marnodes\right|=3$, 300 with
$\left|\marnodes\right|=4$ and two with $\left|\marnodes\right|=5$.  Therefore the
left panel holds 3332 band widths, the right panel 23696.}\par
\label{qmr1000}}
\end{figure}}

\section{Discussion}\label{discussion}

We have shown that bound propagation is a simple algorithm with
surprisingly good results.  It performs exact computations on local parts
of the network and keeps track of the uncertainties that brings along.
In this way it is capable to compute upper and lower bounds on any
marginal probability of a small set of nodes in the intractable network.
Currently we do not understand which network properties are responsible
for the tightness of the bounds found.  In \figureref{lattice}, for
instance, we saw islands of worse results in a sea of tight bounds.
It is obvious that one bad bound influences its direct neighbourhood,
since bound propagation performs local computations in the network.
We have no clue, however, why these islands are at the location we
found them.  We tried to find a correlation with the size of the weights
(rewriting the network potentials to a Boltzmann distribution) or with
the local frustration of the network (spin glass state), but could
not explain the quality of the bounds in terms of these quantities.
Here we pose it as an open question.

\hide{Experimental results not shown here gave us the impression
that the algorithm has difficulties to find good bounds if the network
tends to be deterministic (i.e.\ all the probability tables of the belief
network have values rather close to zero or one).  This is not surprising,
since these networks usually have many long range correlations which are
difficult to cover for any algorithm that is based on local
computations.}
The computational complexity of the algorithm is mainly determined by the
state space of the separator nodes one uses, which makes the algorithm
unsuitable for heavily connected networks such as Boltzmann machines.
Nevertheless, there is a large range of architectures for which bound
propagation can easily be done.  Such networks typically have a limited
number of connections per node which makes the majority of the Markov
blankets small enough to be tractable for the bound propagation algorithm.

\hide{There are many parameters that can be chosen freely, such as which node
sets are to be bounded and which separator sets we choose.  There is
definitely a lot of research to be done to determine which node sets
are worth spending the available computational power on.}

We want to discuss one important property of the bound propagation
algorithm we did not address before.  In this article we found
our set of $\sepnodes$'s by \equationref{omega}.  Due to this
choice the separator nodes will always be in the neighbourhood of
$\marnodes$.  We have, however, much more freedom to choose $\sepnodes$.
In \sectionref{theory} we stated that a sufficient condition to make the
algorithm tractable is that $\norm{\sepnodes\cup\marnodes\cup\othnodes}$
is small enough to do exact computations.  A more general, but still
sufficient condition is that we should choose $\sepnodes$ such that
$\prob{\marnodes\,|\,\sepnodes}$ can be computed efficiently, since
this is the quantity we need (see \equationref{c}).  If the network
structure inside the area separated by $\sepnodes$ can be written as
a junction tree with a small enough tree width, we can compute these
conditional probabilities even if the number of nodes enclosed by
$\sepnodes$ is very large.  For certain architectures the separator
can even in that case be small enough.  One can think of a network
consisting of a number of tractable junction trees that are connected
to each other by a small number of nodes.  One should be aware of the
important difference between the method outlined here and the method
of conditioning \shortcite{pearl88a}, which does exact inference.
We do not require the part of the network that is outside the separator
(i.e.\ all nodes not in $\sepnodes$, $\marnodes$ or $\othnodes$) to
be tractable.  The open problem is to develop an efficient algorithm
to find suitable separators, since the nodes in such an $\sepnodes$
are generally spread widely.

To conclude this article we like to say a few more words about linear
programming.  This is a field of research that is still developing.
Any improvements on \textsc{lp}-solving algorithms directly
influence the results presented in this article.  We could imagine
that more advanced \textsc{lp}-methods could benefit from the fact that
the matrix $A$ in \equationref{Ab} is sparse.  At least half of the
entries are zero.  To be more precise:  the constraints induced by a
particular $\marnodesprime$ are exactly $2\norm{\marnodesprime}$ rows
in the matrix $A$ with exactly $2\norm{\sepnodes}$ non-zero entries,
thus having a fraction of $1-1/\norm{\marnodesprime}$ zero�s.  Moreover,
all the non-zero entries are ones.  Finally, there is a promising future
for \textsc{lp}-solvers, since the algorithms seems to be suitable for
parallel computing \shortcite{alon90parallel}.


\hide{But we see also other directions for future research.  In this article
we took the bounds on the marginals itself to restrict the uncertain
range of $\prob{\sepstate}$.  These lead to the inequality constraints
in the \textsc{lp}-problem.  There is no reason, however, why we could
not choose other constraints that do not correspond to marginals.
Suppose we have a probability distribution over the two variables $s_1$
and $s_2$.  We could try, for example, to find an upper and lower
bound on $\prob{s_1=1,s_2=1}+2\prob{s_1=1,s_2=2}-3\prob{s_1=2,s_2=2}$.
It might be the case that the upper and lower bounds on this quantity
are tighter than those on the marginals.  Once we have such bounds, we
can use them as a constraint in our \textsc{lp}-problem every time $s_1$
and $s_2$ appear together in $\sepnodes$.}

\acks{This research is supported by the Technology Foundation
\textsc{stw}, applied science devision of \textsc{nwo} and the
technology programme of the Ministry of Economic Affairs.}

\bibliography{references}
\bibliographystyle{theapa}

\end{document}

