\documentclass[twoside,11pt]{article}


\usepackage{jair}
\usepackage{theapa}
\usepackage{latexsym}


\input{epsf.tex}
\input{macros}

\def\definedas{{\:\:\:\:\stackrel{\it def}{=}\:\:\:\:\:\:\:}}


\def\Q{{\bf Q}}
\def\q{{\bf q}}
\def\E{{\bf E}}
\def\e{{\bf e}}


\def\MAXSNP{{\it MAXSNP}}

\def\dpr{{\bf D-PR}}
\def\dmpe{{\bf D-MPE}}
\def\dmap{{\bf D-MAP}}
\def\sat{{\bf SAT}}
\def\msat{{\bf MAJSAT}}
\def\emsat{{\bf E-MAJSAT}}
\def\maxsat{{\bf MAXSAT}}

\jairheading{21}{2004}{101-133}{3/03}{2/04}
\ShortHeadings{Complexity Results and Approximation Strategies for
  MAP Explanations}{Park \& Darwiche}
\firstpageno{101}
\begin{document}
\title{Complexity Results and Approximation Strategies for
MAP Explanations}

\author{\name James D. Park \email jd@cs.ucla.edu \\
\name Adnan Darwiche \email darwiche@cs.ucla.edu \\
\addr Computer Science Department  \\
University of California \\
Los Angeles, CA 90095 \\
}

\maketitle
\begin{abstract}
MAP is the problem of finding a most probable instantiation of a
set of variables given evidence. MAP has always been perceived to
be significantly harder than the related problems of computing the
probability of a variable instantiation ($\Pr$), or the problem of
computing the most probable explanation (MPE). This paper
investigates the complexity of MAP in Bayesian networks.
Specifically, we show that MAP is complete for $\NP^\PP$ and
provide further negative complexity results for algorithms based
on variable elimination. We also show that MAP remains hard even
when MPE and $\Pr$ become easy. For example, we show that MAP is
$\NP$-complete when the networks are restricted to polytrees, and
even then can not be effectively approximated.

Given the difficulty of computing MAP exactly, and the difficulty
of approximating MAP while providing useful guarantees on the
resulting approximation, we investigate best effort
approximations. We introduce a generic
 MAP approximation framework.  We provide two instantiations of the
 framework; one for networks which are amenable to exact inference ($\Pr$),
 and one for networks for which even exact inference is too hard.
 This allows MAP approximation on networks that are too complex to
 even exactly solve the easier problems, $\Pr$ and MPE.  Experimental
 results indicate
 that using these approximation algorithms provides much better
 solutions than standard techniques, and provide accurate MAP
 estimates in many cases.
\end{abstract}

\section{Introduction}
The task of computing the Maximum a Posteriori Hypothesis (MAP) is
to find the most likely configuration of a set of variables given
{\em partial} evidence about the complement of that set. The focus
of this paper is on the complexity of computing MAP in Bayesian
networks, and on a class of best effort methods for approximating
MAP.

One specialization of MAP which has received a lot of attention is
the Most Probable Explanation (MPE).  MPE is the problem of
finding the most likely configuration of a set of variables given
{\em complete} evidence about the complement of that set. The
primary reason for this attention is that MPE seems to be a much
simpler problem than its MAP generalization. Unfortunately, MPE is
not always suitable for the task of providing explanations.
Consider for example the problem of system diagnosis, where each
component has an associated variable representing its health.
Given some evidence about the system behavior, one is usually
interested in computing the most probable configuration of health
variables. This is a MAP problem since the available evidence does
not usually specify the value of each non--health variable. It is
common to approximate this problem using MPE, in which case one
is finding the most likely configuration of every unknown
variable, including health variables and some other variables of
no particular interest, such as the inputs and outputs of system
components. However, the projection of an MPE solution on health
variables is usually not a most likely configuration. Neither is the
configuration obtained by choosing the most likely state of each
health variable separately.

MAP turns out to be a very difficult problem, even when compared
to the MPE problem, or to the \(\Pr\) problem for computing the
probability of evidence. Specifically, we provide in
Section~\ref{complexity section} some complexity results which
indicate that neither exact nor approximate solutions can be
guaranteed for MAP, even under very restricted circumstances. Yet,
MAP remains an important problem for which we would like to
generate solutions. Therefore, we propose in
Section~\ref{approximation section I} a general framework based on
local search for best--effort approximation of MAP. We also
provide two specific instances of the proposed framework, one is
applicable to networks that are amenable to exact computation of
\(\Pr\) and is given in Section~\ref{approximation section I},
while the other is applicable to networks that are not even
amenable to \(\pr\) and is given in Section~\ref{approximation
section II}. We report on experimental results for each method
using real--world and randomly generated Bayesian networks, which
illustrate the effectiveness of our proposed framework on a wide
range of networks. We close our paper by some concluding
remarks in Section~\ref{sec:conclusion}.


\section{MAP Complexity}\label{complexity section}
We begin this section by reviewing some complexity theory classes
and terminology that pertain to the complexity of MAP. We then
examine the complexity of MAP in the general case, followed by
examining the complexity when the number of MAP variables is
constrained.  We then consider the complexity of MAP algorithms based
on variable elimination.  We conclude the complexity section by
examining the complexity of MAP on polytrees.


\subsection{Complexity Review}

We assume that the reader is familiar with the basic notions of
complexity theory like the hardness and completeness of languages, as
well as the complexity class $\NP$.  

In addition to $\NP$, we will also be interested in the class $\PP$
and a derivative of it.  Informally, $\PP$ is the class which contains
the languages for which there exists a nondeterministic Turing machine
where the majority of the nondeterministic computations accept if and
only if the string is in the language.  $\PP$ can be thought of as the
decision version of the functional class $\SharpP$.  As such,  $\PP$
is a powerful language.  In fact $\NP\subseteq\PP$, and the inequality
is strict unless the polynomial hierarchy collapses to the second
level.\footnote{This is a direct result of Toda's theorem
  \cite{toda91}.  From Toda's theorem $\P^{\PP}$ contains the entire
  polynomial hierarchy ($\PH$), so if $\NP=\PP$, then $\PH \subseteq
  \P^{\PP}=\P^{\NP}$.}

Another idea we will need is the concept of an oracle.
Sometimes it is useful to ask questions about what could be done if an
operation were free.  In complexity theory this is modeled as a Turing
machine with an oracle.  An oracle Turing machine is a Turing machine
with the additional capability of being able to obtain answers to
certain queries in a single time step.  For example, we may want to
designate the class of languages that could be recognized in
nondeterministic polynomial time if any $\PP$ query could be answered
for free.  The class of languages would be $\NP$ with a $\PP$ oracle,
which is denoted $\NPPP$.

Consider now a Boolean expression \(\phi\) over variables
\(X_1,\ldots,X_n\). The following three classical problems are
complete for the above complexity classes:
\begin{description}
\item[\sat:] Is there a truth assignment (world) that satisfies \(\phi\)?
This problem is \(\NP\)--complete.

\item[\msat:] Do the majority of worlds satisfy \(\phi\)?
This problem is \(\PP\)--complete.

\item[\emsat:] Is there an instantiation of variables
\(X_1\),\ldots,\(X_k\), \(1 \leq k \leq n\), under which the
majority of worlds satisfy \(\phi\)? This problem is
\(\NPPP\)--complete.
\end{description}
Intuitively, to solve an \(\NP\)--complete problem we have to {\em
search} for a solution among an exponential number of candidates,
where it is easy to decide whether a given candidate constitutes a
solution. For example, in \sat, we are searching for a world that
satisfies a sentence (testing whether a world satisfies a sentence
can be done in time linear in the sentence size). To solve a
\(\PP\)--complete problem, we have to {\em add up} the weights of
solutions, where it is easy to decide whether a particular
candidate constitutes a solution and it is also easy to compute
the weight of a solution. For example, in \msat, a solution is a
world that satisfies the sentence and the weight of a solution is
\(1\). Finally, to solve an \(\NPPP\)--complete problem, we have
to search for a solution among an exponential number of
candidates, but we also need to solve a \(\PP\)--complete problem
in order to decide whether a particular candidate constitutes a
solution. For example, in \emsat, we are searching for an
instantiation \(x_1,\ldots,x_k\), but to test whether an
instantiation satisfies the condition we want, we must solve a
\msat\ problem.

\subsection{Decision Problems}

We will be dealing with the decision versions of Bayesian network
problems in this paper, which we define formally in this section.

A Bayesian network is a pair \((G,\Theta)\), where \(G\) is a
directed acyclic graph (DAG) over variables \(\X\), and \(\Theta\)
defines a conditional probability table (CPT) \(\Theta_{X|\U}\)
for each variable \(X\) and its parents \(\U\) in the DAG \(G\).
That is, for each value \(x\) of variable \(X\) and each
instantiation \(\u\) of parents \(\U\), the CPT \(\Theta_{X|\U}\)
assigns a number in \([0,1]\), denoted by \(\theta_{x|\u}\), to
represent the probability of \(x\) given \(\u\).\footnote{Hence,
we must have \(\sum_x \theta_{x|\u}=1\).} The probability
distribution \(\Pr\) induced by a Bayesian network \((G,\Theta)\)
is given as follows. For a complete instantiation \(\x\) of the
network variables \(\X\), the probability of \(\x\) is given by
\[ \pr(\x) \definedas \prod_{x\u \sim \x} \theta_{x|\u}, \]
where \(x\u\) is the instantiation of a family (a variable and its
parents) and \(\sim\) represents the compatibility relation among
instantiations. That is, the probability assigned to a complete
variable instantiation \(\x\) is the product of all parameters
that are consistent with that instantiation.

The following decision problems assume that we are given a
Bayesian network \((G,\Theta)\) that has rational parameters and
induces the probability distribution \(\pr\). Moreover, by
evidence \(\e\), we mean an instantiation of variables \(\E\).
\begin{description}
\item[\dmpe:] Given a rational number \(p\), evidence \(\e\), and the set
of network variables \(\X\), is there an instantiation \(\x\) such
that \(\pr(\x,\e) > p\)?

\item[\dpr:] Given a rational number \(p\) and evidence \(\e\), is
\(\pr(\e) > p\)?

\item[\dmap:] Given a rational number \(p\), evidence \(\e\), and some
set of variables \(\Q\), is there an instantiation \(\q\) such
that \(\pr(\q,\e) > p\)? Variables \(\Q\) are called the MAP
variables in this case.
\end{description}
%It is well known that \dmpe\ is \(\NP\)--complete\footnote{The
%NP-hardness of the functional version of
%  MPE was shown in \cite{shimony94}. We are not aware of any published
%  proof of the completeness of the decision problem, so we sketch it
%  here.  Membership is immediate, since the score for a purported
%  solution can be tested in linear time.  Hardness is based on using
%  the standard Bayesian network simulation of a Boolean expression
%  (c.f. Theorem~\ref{general theorem}) to solve SAT.
%  MPE$(V_{\phi}=T)>0$ if and only if the expression is satisfiable.}
%and that \dpr\ is \(\PP\)--complete \cite{littman01}.
%We show in the
%following section that \dmap\ is \(\NPPP\)--complete.
While decision problems are useful for examining the complexity of
finding an exact solution, what we are really interested in is the
functional problem of actually computing the solution.  When we can't
solve the problem exactly, we would also like to know how close we
can get efficiently.  For that we consider
approximation algorithms.
We now define the notion of an {\em approximation factor} which we
will use when discussing the complexity of approximation
algorithms. Specifically, we will say that an approximate solution
$M'$ is within the approximation factor $\epsilon >1$ of the true
solution $M$ in case $\frac{M}{\epsilon} \leq M' \leq \epsilon M$.
Moreover, we will say that an algorithm provides an $f(n)$--factor
approximation in case for all problems of size $n$, the
approximate solutions returned by the algorithm are within the
approximation factor $f(n)$.

\subsection{MAP Complexity for the General Case}
Computing MPE, $\Pr$, and MAP are all \(\NP\)--Hard, but there
still appears to be significant differences in their complexity.
MPE is basically a combinatorial optimization problem. Computing
the probability of a complete instantiation is trivial, so the only
real difficulty is
determining which instantiation to choose. \dmpe\ is $\NP$-complete
\cite{shimony94}. \(\Pr\) is a completely different type of
problem, characterized by counting instead of optimization, as we
need to add up the probability of network instantiations. \dpr\
is $\PP$-complete \cite{littman01}---notice that this is the
complexity of the decision version, not the functional version
which is $\SharpP$-complete \cite{roth96}. MAP combines both the
counting and optimization paradigms.  In order to compute the
probability of a particular instantiation, a \(\pr\) query is
needed. Optimization is also required, in order to be able to
decide between the many possible instantiations. This is reflected
in the complexity of MAP.
\begin{theorem}\label{general theorem}
\dmap\ is $\NPPP$-complete.\footnote{This result was stated
without proof
  by \citeA{littman99}.}
\end{theorem}
{\em Proof:} Membership in $\NPPP$ is immediate.  Given any
instantiation $\q$ of the MAP variables, we can verify if it is a
solution by querying the $\PP$ oracle if $\Pr(\q,\e)>k$.

To show hardness, we reduce \emsat\ \cite{littman98} to \dmap\ by
first creating a Bayesian network that models a Boolean formula
\(\phi\). For each variable \(X_i\) in the formula \(\phi\), we
create an analogous variable in the network with values
\(\{T,F\}\) and a uniform prior probability. Then, for each
logical operator, we create a variable with values \(\{T,F\}\)
whose parents are the variables corresponding to its operands, and
whose CPT encodes the truth table for that operator (see
Figure~\ref{circuit figure} for a simple example). Let $V_\phi$ be
the network variable corresponding to the top level operand.

For a complete instantiation $\x$ of all of the variables \(\X\)
appearing in the Boolean expression \(\phi\), with evidence $\V_\phi=T$, we have:
\[\Pr(\x,V_\phi=T)=
\left\{
\begin{array}{ll}
\frac{1}{2^n} & \mbox{$\x$ satisfies $\phi$} \\
0 & \mbox{otherwise}
\end{array}
\right.
\]
For a particular instantiation $\q$ of MAP variables
$X_1,...,X_k$, and evidence $V_\phi=T$, we have:
\begin{eqnarray*}
\Pr(\q,V_\phi=T)& = &\sum_{x_{k+1},...,x_n} \Pr(\q,x_{k+1},..x_{n})\\
& = & \frac{\#_\q}{2^n}
\end{eqnarray*}
where $\#_{\q}$ is the number of complete variable instantiations
compatible with $\q$ that satisfies $\phi$. Since there are
$2^{n-k}$ possible instantiations of $X_{k+1},...,X_n$, the
fraction $f_\q$ satisfied is $\#_\q/2^{n-k}$, so
\[
\Pr(\q,V_\phi=T)=\frac{f_\q}{2^k}
\]
Thus, an instantiation $\q$ of the MAP variables is compatible
with more than half of the complete, satisfying instantiations if
$\Pr(\q,V_\phi=T)> 1/2^{k+1}$.  So the MAP query over
variables $X_1,...,X_k$ with evidence $V_\phi=T$ and threshold $1/2^{k+1}$ is true if
and only if the \emsat\ query is also true. $\Box$

\begin{figure}[tb]
\centerline{\epsfxsize=50mm \leavevmode \epsffile{circuit.eps}}
\caption{The Bayesian network produced using the reduction in
  Theorem~\ref{general theorem} for Boolean formula $\lnot(x_1 \lor
  x_2)\land \lnot x_3$.}\label{circuit figure}
\end{figure}

In fact, the above theorem can be strengthened.
\begin{theorem}\label{big theorem}
\dmap\ remains complete for $\NPPP$ even when (1) the network has
depth 2, (2)  there is no
evidence, (3) all variables are Boolean, and
(4) all network parameters lie in the interval
$[\frac{1}{2}-\epsilon,\frac{1}{2}+\epsilon]$ for any fixed $\epsilon>0$.
\end{theorem}
The proof appears in Appendix~\ref{big theorem proof}. Unlike
computing probabilities,  which becomes easy when the number of
evidence nodes is bounded by a constant and the parameters are
bounded away from 0 (it falls into $\RP$ as described by \citeR{dagum97}), MAP
retains its $\NPPP$ complexity even under these restrictive
circumstances.\footnote{This is not altogether surprising since
when evaluating the score of a possible solution, the MAP
variables act as evidence variables.}

$\NPPP$ is a powerful class, even compared to $\NP$ and $\PP$.
$\NPPP$ contains other important AI problems, such as
probabilistic planning problems \cite{littman98}. The three
classes are related by $\NP\subseteq\PP\subseteq\NPPP$, where the
equalities are considered very unlikely. In fact, $\NPPP$ contains
the entire polynomial hierarchy \cite{toda91}. Additionally,
because MAP generalized $\Pr$, MAP inherits the wild
non--approximability of $\Pr$.  From the ``Bayesian network
simulates SAT'' reduction we get:
\begin{corollary}
Approximating MAP within any approximation factor $f(n)$ is NP--hard.
\end{corollary}
{\bf Proof:} Using the evidence $V_\phi=T$, the exact MAP
solution is the number of satisfying instantiations divided by $2^n$,
which is 0 if it is unsatisfiable, and positive if it is
satisfiable.  If the formula is unsatisfiable, then the approximate
solution must be 0 because $0/\epsilon=0\leq M'\leq \epsilon 0=0$,
where $\epsilon=f(n)$.  If the formula is satisfiable, then the
approximate solution must be positive since $M/\epsilon>0$.  Thus we can test
satisfiability by testing if the approximate MAP solution is zero or
not.$\Box$
\subsection{Complexity Parameterized by the Number of Maximization
  Variables}
We now examine the complexity when we restrict the number of
maximization variables.  Let $n$ be the number of non--evidence
variables, and $k$ be the number of maximization variables.  In the
extreme case of $k=0$, this is simply \dpr, so it is
$\PP$--complete.  At the other extreme, when $k=n$, it becomes \dmpe,
and is $\NP$--complete.  So constraining the number of maximization
variables can have a dramatic impact on the complexity.  We now
examine this issue in detail.  Let $\mbox{\dmap}_{m}$ be the subset of
\dmap\ problems where $k=O(m)$, and let $\mbox{\dmap}^{m}$ be the
subset of \dmap\ problems where $n-k=O(m)$.  We can then consider the
complexity of these parameterized classes of problems.  The primary
results are the following:

\begin{theorem}
$\mbox{\dmap}_{\log n}$ is in $\P^{\PP}$, and
$\mbox{\dmap}^{\log n}$ is in $\NP$.  However, for any $\epsilon>0$,
  both $\mbox{\dmap}_{n^\epsilon}$ and $\mbox{\dmap}^{n^\epsilon}$
  remain $\NPPP$--complete. 
\end{theorem}
{\em Proof:}  First, if $k=O(\log n)$, then the number of possible
instantiations of the maximization variables is bounded by a
polynomial.  Thus, given a $\PP$ oracle, it is possible to decide the
problem in polynomial time by asking for each instantiation $\q$ of
the maximization 
variables whether $\Pr(\q,\e)$ exceeds the threshold.  Similarly, if
$n-k=O(\log n)$, then for any instantiation $\q$ of the maximization
variables, we can test to see if $\Pr(\q,\e)$ exceeds the threshold
by summing over the polynomial number of compatible instantiations.

For $k=O(n^\epsilon)$ we can provide a
simple reduction to solve any \dmap\ problem by creating a polynomially
larger one satisfying the constraint on the number of maximization
variables.  From the unconstrained problem, we simply create a new 
problem by adding a polynomial number of irrelevant variables, with
no parents or children.  Similarly, we can provide a reduction of the
general \dmap\ problem to one constrained to have $n-k=O(n^\epsilon)$,
by adding a polynomial number of maximization variables with no
parents, no children, and deterministic priors.
$\Box$
\subsection{Results for Elimination Algorithms}
Solution to the general MAP problem seems out of reach, but what
about
for ``easier'' networks?  State--of--the--art exact inference
algorithms (variable elimination \cite{dechterUAI96}, join trees
\cite{lauritzen88,shenoy86,jensen90}, recursive conditioning
\cite{darwicheAIJ-RC}) can compute $Pr$ and MPE in space and time
complexity that is exponential only in the width of a given elimination
order.  This allows many networks to be solved using reasonable
resources even though the general problems are very difficult.
Similarly, state--of--the--art MAP algorithms can also solve MAP with time
and space complexity that is exponential only in width of used
elimination order but, for MAP, not all orders can be
used.  In this section, we investigate the complexity of variable
elimination for MAP.


Before analyzing the complexity of variable elimination for MAP,
we review variable elimination. First, we need the concept of a
potential. A potential is simply a function over some subset of
the variables, which maps each instantiation of its variables to a
real number.  The size of a potential is parameterized by the
number of instantiations of its variables, and so is exponential
in the number of variables. Notice that CPTs  are potentials.  In
order to use variable elimination for $\Pr$, MPE and MAP, we need
three simple operations: multiplication, summing--out, and maximization.  Multiplication of two potentials $\phi_1$ and
$\phi_2$ with variables $\X\Y$ and $\Y\Z$ respectively (where $\Y$
is the set of variables they have in common), is defined as
$(\phi_1\phi_2)(\x\y\z)=\phi_1(\x\y)\phi_2(\y\z)$.  Notice that if
both $\X$ and $\Z$ are nonempty, then the size of $\phi_1\phi_2$
is greater than the size of either $\phi_1$ or $\phi_2$.
SumOut($\phi$,$\Y$) where $\phi$ is over variables $\X\Y$ is
defined as
\[\mbox{SumOut}(\phi,\Y)(\x)=\sum_{\y} \phi(\x\y),\]
where $\y$ ranges over the instantiations of $\Y$. Maximization is
similar to summing out but maximizes out the unwanted variables:
\[\mbox{Maximize}(\phi,\Y)(\x)=\max_{\y} \phi(\x\y).\]

In order to handle evidence, we need the concept of an evidence
indicator.  The evidence indicator $\lambda_E$ associated with evidence
$E=e$ is a potential over variable $E$ where $\lambda_E(e)=1$, and is
zero for all other values.  

Given a variable ordering \(\pi\), variable elimination can be
used to compute the probability of evidence \(\e\) as follows:
\begin{enumerate}
\item Initialize \(P\) to contain the evidence indicators for $\e$
  and all of the conditional probability tables.
\item For each variable $X$, according to order \(\pi\),
  \begin{enumerate}
  \item Remove from $P$ all potentials mentioning $X$.
  \item Let $M_X$ be the product of all of those potentials.\label{ve
    multiplication step}
  \item add SumOut($M_X,X$) to $P$.\label{ve project step}
  \end{enumerate}
\item Return the product of all potentials in $P$.\label{ve last step}
\end{enumerate}
In each iteration, a variable $X$ is eliminated which leads to
removing all mention of $X$ from $P$. By step~\ref{ve last step},
all variables have been removed, so the potentials remaining are
just constants and the resulting product is a single number
representing the probability of evidence \(\e\). MPE is computed
in the same way, except the projection in step~\ref{ve project
step} is replaced by maximization. The complexity of variable
elimination is linear in the number of variables and linear in the
size of the largest potential $M_X$ produced in step~\ref{ve
multiplication step}. The size of the largest potential varies
significantly based on the elimination order. The {\em width} of
an elimination order is simply $\log_2 (\mbox{size}(M_X))-1$ where
$M_X$ is the largest potential produced using that elimination
order.\footnote{The -1 in the definition is to preserve
compatibility
  with the previously defined notion of treewidth in graph theory.}
  The {\em treewidth} of a Bayesian network is the minimum
width of all elimination orders. For $\Pr$ and MPE, any elimination
order can be used, so the complexity is linear in the number of
variables and exponential in the treewidth.  The same is not true for
MAP.  Variable elimination for MAP is similar to the other methods,
but with an extra constraint.  In step~\ref{ve project step}, if $X$
is a MAP variable the projection is replaced with maximization.  If
it is not a MAP variable projection is used.  The extra constraint is
that not all orders are valid.  Maximization and projection do not
commute, and maximization must be performed last.  This means that
for an elimination order to be valid for performing MAP, when $X$ is
a MAP variable, the potential in step~\ref{ve multiplication step}
must not mention any non-MAP variables.
In practice this is ensured by requiring that
the elimination order eliminate all of the MAP variables last.
This tends to produce elimination orders with widths much larger
than those available for $\Pr$ and MPE, often placing exact MAP
solutions out of reach.

%In order to assess the magnitude of increase in width caused by
%restricting elimination orders, we generated 1000 Bayesian
%networks randomly as given in Appendix~\ref{random networks} and
%then computed the constrained and unconstrained elimination orders
%for these networks using the min--fill heuristic
% \cite{kjaerulff90,darwicheIJAR96}. For constrained orders, all non--MAP variables
%were eliminated first. Each network had 100 nodes and the set of
%MAP variables consisted of 10--25 root nodes. For all of the networks
%with a particular unconstrained width, we measured the
%minimum, maximum, average and weighted average of the constrained
%width.  The average was computed as $\sum_{i=1}^{k}
%w_i/k$  where $w_1,\ldots,w_k$ are the constrained widths for the $k$
%networks having the same unconstrained width.  Since the complexity is exponential in
%the width, a weighted sum gives a better representation of the
%average complexity. It was computed as $\log_2(\sum_{i=1}^k
%2^{w_i}/k)$. Figure \ref{width change} summarizes the results.

%\begin{figure*}[tb]
%\centerline{\epsfxsize=100mm \leavevmode \epsffile{mapwidths.ps}}
%\caption{The minimum, maximum, average and weighted average of the
%  constrained widths. The X axis is the unconstrained width, and the Y
%  axis is the constrained width.  Notice that the constrained
%width can grow to be unmanageable even for networks with small
%unconstrained width. \label{width change}}
%\end{figure*}

%In many cases, the constrained width was much larger than the
%unconstrained width, often making the MAP problem unreasonably
%expensive, even when the MPE problem could be solved exactly with
%reasonable resources. For example, the weighted average width
%increased from about 13 to about about 27 due to MAP constraints.

In order to assess the magnitude of increase in width caused by
restricting the elimination order, we randomly generated 100 Bayesian
networks, each containing 100 variables, 
according to the first method in Appendix B.  For each
network, we computed the width using the min--fill heuristic
\cite{kjaerulff90,darwicheIJAR96}. Then, we repeatedly added a single
variable to the set of MAP variables, and computed the constrained
width, again using min--fill, but eliminating the MAP variables last.
This process was repeated until all variables were in the MAP variable set.
Figure~\ref{width experiments figure} contains statistics on these
experiments.  The $X$--axis corresponds to the number of MAP
variables (thus $X=0$ corresponds to the unconstrained width).  The
$Y$--axis corresponds to the width found.  The graph details the
minimum, maximum, mean, and weighted mean for each of the 100
networks.  The weighted mean takes into account that the complexity
is exponential in the width, and so provides a better representation
of the average complexity.  It was computed as
$\log_2(\frac{1}{n}\sum_{i=1}^n 2^{w_i})$.   Notice that the
unconstrained widths range from 11 to 18, and that as the number of
MAP variables increases, the width increases dramatically.  For
example, even when only a quarter of the variables are MAP variables ($X=25$)
the widths range between 22 and 34, (which corresponds roughly from
difficult but doable to well beyond what today's inference algorithms
can handle on today's computers) with a weighted average above 30. 
Notice also, that as we would expect from our complexity analysis,
problems with very few or very many MAP variables are significantly
easier than those in the middle range.
\begin{figure}[tb!]
\centerline{\epsfxsize=5in \leavevmode
\epsffile{widths.eps}} 
\caption{Statistics on the constrained width for different numbers of
  MAP variables.  The $X$--axis is the number of MAP variables, and the
  $Y$--axis is the width.  Notice that the widths required for
  the general MAP 
problem are significantly larger than for  Pr  or MPE, which correspond
to $X=0$ and $X=100$ respectively.}\label{width experiments figure}
\end{figure}

We now consider the question of whether there are less stringent
conditions for valid elimination orders, that may allow for orders
with smaller widths.

As described earlier, given an ordering, elimination
algorithms work by stepping through the ordering, collecting the
potentials mentioning the current variable, multiplying them, then
replacing them with the potential formed by summing out (or
maximizing) the current variable from the product. This process
can be thought to induce an {\em evaluation tree;} see
Figure~\ref{elimination tree figure}. The evaluation tree for an
elimination order \(\pi\) is described as follows. The leaves
correspond to the CPTs of given Bayesian network, and the internal
nodes correspond to potentials created during the elimination
process. The children of a potential \(\phi\) represent other
potentials that had to be multiplied together when constructing
\(\phi\). Note that each internal node in the elimination tree
corresponds to variable in the order \(\pi\), whose elimination
leads to constructing the node; Figure~\ref{elimination tree
figure}(b). Therefore, an evaluation tree can be viewed as
inducing a partial elimination order; see Figure~\ref{elimination
tree figure}(c).

\begin{figure*}[tb]
\centerline{\epsfxsize=120mm \leavevmode
\epsffile{eliminationtree.eps}} \caption{(a) A Bayesian network,
(b) an evaluation tree for the elimination order \(A, E, B, D,
E\), and (c) a partial elimination order induced by the evaluation
tree.}\label{elimination tree figure}
\end{figure*}

The standard way of constructing a valid elimination order for MAP
is to eliminate the MAP variables \(\Q\) last. Two questions present
themselves. First, are there valid orderings in which variables
\(\Q\) are not eliminated last? And second, if so, can they
produce widths smaller than those generated by eliminating
variables \(\Q\) last?

The answer to the first question is yes, there are other valid
elimination orders in which variables \(\Q\) are not eliminated
last. To see that, suppose we have a variable order \(\pi\) which
induces a particular evaluation tree \(T\), and let \(\sigma\) be
the partial elimination order corresponding to \(T\). Any variable
order \(\pi^\prime\) which is consistent with the partial order
\(\sigma\) will also induce the tree \(T\). Hence, if order
\(\pi\) was valid, then order \(\pi^\prime\) will also be valid.
Figure~\ref{elimination tree figure} shows the evaluation tree
induced by using the order \(A,E,B,D,C\) for computing MAP over
variables \(\Q = C, D\). The order \(A, B, D, E, C\) is consistent
with this evaluation tree and, hence, is also valid for computing
MAP over variables \(C, D\). Yet, variables \(C, D\) do not appear
last in the order.

Unfortunately, orders in which variables \(\Q\) are not eliminated
last do not help.

\def\nQ{\overline{\bf Q}}

\begin{theorem}\label{order theorem}
For any elimination order \(\pi\) which is valid for computing MAP
over variables \(\Q\), there is an ordering of the same width in
which variables \(\Q\) are eliminated last.
\end{theorem}
{\em Proof:} Consider the evaluation tree induced by any valid
elimination order, and the corresponding partial order it induces.
No variable in \(\nQ\) is the parent of any variable in \(\Q\). To
prove this, suppose that \(X\) is a parent of \(Y\) in the
evaluation tree, where \(\X \in \nQ\) and \(\Y \in \Q\). This
means that the potential which results from eliminating variable
\(Y\) includes variable \(X\), which also means that \(X\) must
have appeared in some of the potentials we multiplied when
elimination variable \(Y\). But this is a contradiction since the
evaluation tree and its underlying order are valid. Since no
variable in \(\nQ\) is a parent of a variable in \(\Q\), all
variables in \(\nQ\) can be eliminated first in any order
consistent with the partial order defined by the evaluation tree.
Then, all variables in \(\Q\) can be eliminated, again obeying the
partial ordering defined by the evaluation tree. Because the
produced order has the same elimination tree as the original
order, they have the same width. $\Box$

\subsection{MAP on Polytrees}
Theorem~\ref{order theorem} has significant complexity implications
for elimination algorithms even on polytrees.
\begin{theorem}\label{polytree order}
Elimination algorithms require exponential resources to perform MAP,
even on some polytrees.
\end{theorem}
{\em Proof:} Consider computing MAP over variables
\(X_1,\ldots,X_n\) given evidence \(S_n=T\) for the network shown
in Figure~\ref{reduction figure}. By Theorem~\ref{order theorem},
there is no order whose width is smaller than that of an order
\(\pi\) in which we eliminate variables \(S_0,\ldots,S_n\) first,
and then variables \(X_1,\ldots,X_n\) last. It is easy to show
though that any such order \(\pi\) will have width \(n\). Hence,
variable elimination will require exponential resources using such
an order. $\Box$

The set of MAP variables makes a crucial difference in the
complexity of MAP computations.  For example, if the MAP variables
were $X_1, \ldots, X_{n/2},S_0, \ldots, S_{n/2}$ instead of $X_1,
\ldots, X_n$ can be solved in linear time.


The above negative findings are specific to variable elimination
algorithms.  The question then is whether this difficulty is an
idiosyncrasy of variable elimination which can be avoided if we
were to use some other method for computing MAP. The following
result, however, shows that finding a good general algorithm for
MAP on polytrees is unlikely.

\begin{theorem}\label{polytree theorem}
MAP is \(NP\)--complete when restricted to polytrees.
\end{theorem}

{\em Proof:} Membership is immediate.  Given a purported solution
instantiation $\q$, we can compute $\Pr(\q,\e)$ in linear time and
test it against the bound. To show hardness, we reduce \maxsat\ to
MAP on a polytree.\footnote{Actually, we only need to reduce it to
  SAT, but the \maxsat\ result will be used in Theorem~\ref{big time
    polytree approximation theorem}.}  Similar reductions were used by
 \citeA{papadimitriou87} and \citeA{littman98} relating to partially
 observable Markov decision problems, and probabilistic planning respectively.
The \maxsat\ problem is
defined as follows:
\begin{quote}
 Given a set of clauses $C_1,...,C_m$ over variables
$Y_1,...,Y_n$
  and an integer bound $k$, is there an assignment of the variables,
  such that more than $k$ clauses are satisfied.
\end{quote}
The idea behind the reduction is to model the random selection of
a clause, and then successively checking whether the instantiation
of each variable satisfies the selected clause. The clause
selector variable $S_0$ with possible values $1,2,...,m$ has a
uniform prior.  Each propositional variable $Y_i$ induces two
network variables $X_i$ and $S_i$, where  $X_i$ represents the
value of $Y_i$, and has a uniform prior, and $S_i$ represents
whether any of $Y_1,...,Y_i$ satisfy the selected clause.  $S_i=0$
indicates that the selected clause was satisfied by one of
$Y_1,...,Y_i$. $S_i=c>0$ indicates that the selected clause $C_c$
was not satisfied by $Y_1,...,Y_i$.  The parents of $S_i$ are
$X_i$ and $S_{i-1}$ (the topology is shown in
Figure~\ref{reduction figure}).  The CPT for $S_i$, for $i\ge 1$
is defined as
\[Pr(s_i|x_i,s_{i-1})=\left\{ \begin{array}{ll}
1 & \mbox{if $s_i=s_{i-1}=0$} \\
1 & \mbox{if $s_i=0$ and $s_{i-1}=j$, and}\\
  & \mbox{$x_i$ satisfies $c_j$}\\
1 & \mbox{if $s_i=s_{i-1}=j$ and $x_i$ does }\\
  & \mbox{not satisfy $c_j$}\\
0 & \mbox{otherwise}
\end{array}\right. \]
In words, if the selected clause was not satisfied by the first
$i-1$ variables ($s_{i-1}\neq 0$), and $x_i$ satisfies it, then
$S_i$ becomes satisfied ($s_i=0$). Otherwise, we have
$s_i=s_{i-1}$. Now, for a particular instantiation $s_0$ of $S_0$,
and instantiation $\x$ of variables $X_1,...,X_n$,
\[
\Pr(s_0,\x,S_n=0) =
 \left\{
\begin{array}{ll}
1/(m 2^n) & \mbox{if $\x$ satisfies clause $C_{s_0}$;}\\
0 & \mbox{otherwise.}
\end{array}
\right.
\]

Summing over $S_0$ yields
$Pr(\x,S_n=0) = \#_C/(m 2^n)$
where $\#_C$ is the number of clauses that $\x$ satisfies.
Thus MAP over $X_1,...,X_n$ with evidence $S_n=0$ and bound $k/(m
2^n)$ solves the \maxsat\ problem as well. $\Box$
\begin{figure}[tb]

\centerline{\epsfxsize=50mm \leavevmode \epsffile{reduction.eps}}

\caption{The network used in the reduction of Theorem ~\ref{polytree theorem}.}\label{reduction figure}

\end{figure}


Since \maxsat\ has no polynomial time approximation scheme (unless $\P=\NP$), 
no polynomial time approximation scheme exists for MAP on
polytrees. In fact, approximating MAP on polytrees
appears to be much harder than approximating MAXSAT.

\begin{theorem}\label{big time polytree approximation theorem}
Approximating MAP on polytrees within a factor of $2^{n^\epsilon}$ is
$\NP$-hard for any fixed $\epsilon$, $0\leq \epsilon < 1$, where $n$ is the size of the problem.
\end{theorem}
The proof appears in Appendix~\ref{proof appendix}.
So, not only is it hard to approximate within a constant factor, it is
hard to approximate within a polynomial factor, or even a subexponential factor.

We close this section by a summary of the complexity results in
this section:
\begin{itemize}
\item MAP is \(\NPPP\)--complete for arbitrary Bayesian networks,
even if we have no evidence, every variable is binary, and the
parameters are arbitrarily close to \(1/2\).

\item It is $\NP$--hard to approximate MAP within any factor $f(n)$.

\item Variable elimination for MAP requires exponential time, even on
  some polytrees.

\item MAP is \(\NP\)--complete for networks with polytree
structure.

\item Approximating MAP on polytrees within a factor of
  $2^{n^\epsilon}$ is \(\NP\)--hard for any fixed $\epsilon\in[0,1)$.

\end{itemize}


\section{Approximating MAP when Inference is Easy}
\label{approximation section I}

Since exact MAP computation is
often intractable, approximation techniques are needed.
Unfortunately, in spite of MAP's utility, relatively little work
has been done in approximating it.  In fact, there are only two
previous methods for approximating MAP which we are aware of.  The
first \cite{dechter98} uses the mini--bucket technique.  The other
\cite{campos99}, uses genetic algorithms to approximate the best k
configurations of the MAP variables (this problem is known as
partial abduction). Practitioners typically resort to one of two
simple approximation methods. One common approximation technique
is to compute an MPE instantiation and then project the result on
the MAP variables. That is, if we want to compute MAP for
variables \(\S\) given evidence \(\e\), and if \(\SS^\prime\) is
the complement of variables \(\SS \cup \E\), we compute an
instantiation \(\s,\s^\prime\) that maximizes
\(\pr(\s,\s^\prime\mid\e)\) and then return \(\s\). The other
method computes posterior marginals for the MAP variables,
\(\pr(S\mid\e), S \in \SS\), and then choose the most likely state
\(s\) of each variable given \(\e\).

We propose a general framework for approximating MAP.  MAP consists of
  two problems that are hard in general---optimization and inference.
  A MAP approximation algorithm can be produced by substituting
  approximate versions of either the optimization or inference
  component (or both).  The optimization problem is defined over the
  MAP variables, and the score for each solution candidate
  instantiation $\s$ of the MAP variables is the (possibly
  approximate) probability $\Pr(\s,\e)$ produced by the inference
  method.
  This allows solutions tailored to the specific problem.  For
  networks whose treewidth is manageable, but contains a hard
  optimization component (e.g. the polytree examples discussed
  previously), exact structural inference can be used, coupled with an
  approximate optimization algorithm.  Alternatively, if the
  optimization problem is easy (e.g. there are few MAP variables) but
  the network isn't amenable to exact inference, an exact optimization
  method could be coupled with an approximate inference routine.  If
  both components are hard, both the optimization and inference
  components need to be approximated.

  We investigate in this section a family of approximation algorithms
  based on local search.  We first consider the case when inference is
  tractable, then develop an extension to handle the case when
  inference is infeasible. The local search algorithms work basically as follows:
\begin{enumerate}
\item Start from an initial guess \(\s\) at the solution.

\item Iteratively try to improve the solution by moving
to a better neighbor \(\s^\prime\): \(\pr(\s^\prime\mid\e)
> \pr(\s\mid\e)\), or equivalently \(\pr(\s^\prime,\e) >
\pr(\s,\e)\).
\end{enumerate}
A {\em neighbor} of instantiation \(\s\) is defined as an
instantiation which results from changing the value of a single
variable \(X\) in \(\s\). If the new value of \(X\) is \(x\), we
will denote the resulting neighbor by \(\s-X,x\). In order to
perform local search efficiently, we need to compute the scores
for all of the neighbors \(\s-X,x\) efficiently. That is, we need
to compute \(\pr(\s-X,x,\e)\) for each \(X\in\SS\) and each of its
values \(x\) not in \(\s\). If variables have binary values, we
will have \(\mid\SS\mid\) neighbors in this case.

Local search has been proposed as a method for approximating MPE
 \cite{kask99,mengshoel00}. For MPE, the MAP variables \(\SS\)
contain all variables which are not in \(\E\) (the evidence
variables).  Therefore, the score of a neighbor,
\(\pr(\s-X,x,\e)\), can be computed easily since \(\s-X,x,\e\) is
a complete instantiation. In fact, given that we have computed
\(\pr(\s,\e)\), the score \(\pr(\s-X,x,\e)\) can be computed in
constant time.\footnote{This assumes that none of entries in the CPTs
  are 0.  If there are 0 entries in the CPTs, it may take time linear
  in the number of network variables to compute the score.
  \(\pr(\s,\e)\) is the product of the single entry of each CPT that
  is compatible with \(\s,\e\).  When changing the state of variable
  $X$ from $x$ to $x'$, the only values in the product that change are
  those from the CPTs of $X$ and its children.  If none of the CPT
  entries are 0, \(Pr(\s-X,x',\e)\) can be computed by dividing
  \(\pr(\s,\e)\) by the old and multiplying by the new entry for the
  CPTs for $X$ and its children.  This can be done in constant time if
  the number of children is bounded by a constant.}

Unlike MPE, computing the score of a neighbor, \(\pr(\s-X,x,\e)\),
in MAP requires a global computation since \(\s-X,x,\e\) may not
be a complete instantiation. One of the main observations
underlying our approach, however, is that the score
\(\pr(\s-X,x,\e)\) can be computed in $O(n \exp(w))$ time and
space, where $n$ is the number of network variables and
 $w$ is the width of an arbitrary elimination order, i.e., we can
 use any elimination order for this purpose, no need for any
 constraints. In fact, we can even do better
 than this by computing the scores of all neighbors,
 \(\pr(\s-X,x,\e)\)  for all \(X \in \SS\) and every value
 \(x\) of \(X\), in $O(n \exp(w))$  time and space. 
Thus, if we have an elimination order of width \(w\) for the given
Bayesian network, then we can perform each search step in $O(n
\exp(w))$ time and space. As we shall see later, it takes a small
number of search steps to obtain a good MAP solution. Hence, the
overall runtime is often $O(n \exp(w))$ too. Therefore, we can
produce good quality MAP approximations in time and space which are
exponential in the unconstrained width instead of the constrained
one, which is typically much larger.



The local search method proposed in this section differs from
the local search methods used for MPE in that the unconstrained width
must be small enough so that a search step can be performed relatively
efficiently.  It is pointless to use this method to approximate MPE
since in the time to take one step, the MPE could be computed exactly.
{\em This method is applicable when the unconstrained width is
  reasonable but the constrained width is not (see Figure~\ref{width
    experiments figure}).}

\subsection{Computing Neighbor Scores Efficiently}
The key to computing the neighbor scores efficiently is to express
the inference problem as a function over the evidence indicators.
For each state $x$ of variable $X$, the evidence indicator
$\lambda_x$ is one if it is compatible with the evidence, and zero
otherwise. This is a common technique that is typically used to allow
the modeling of a wider range of evidence assertions.  For example,
this allows evidence 
assertions such as $X\ne x$ by setting $\lambda_x=0$, and the
remaining indicators for $X$ to one.  We will use them for a
different purpose however.  When the inference problem is cast as a
function $f$ of the evidence indicators ($\f(\lambda_\e)=\Pr(\e)$,
where $\lambda_\e$ consists of all of the evidence indicators, set so
that they are compatible with $\e$), then $\frac{\partial f}{\partial
  \lambda_x}(\lambda_\e)=\Pr(\e-X,x)$.  When our we add the current state to the evidence,
this partial derivative yields $\Pr(\s-X,x,\e)$, which is precisely
the score for one of the neighbors. 

We can use the jointree algorithm
\cite{parkNIPS02}, 
 or the differential inference approach
 \cite{darwicheJACM03} to compute all of the partial derivatives
 efficiently.  In the differential approach, these values are immediate,
 as the entire approach is based on evaluating and differentiating
 the expression $f$ above.  It can also be computed using jointrees
 by using the Shenoy--Shafer propagation scheme.  Specifically, an
 evidence indicator table is added for each variable, and evidence
 about that variable is entered by setting the appropriate indicator
 entries.  The partial 
 derivatives of all of the indicators associated with a variable are
 obtained by multiplying all other tables assigned to the same
 cluster, and all messages into the cluster, then projecting that product
 onto the variable.  In either case, the partial derivatives for all
 indicators, and thus the score for all neighbors, can be computed in
 $O(n \exp(w))$ time, which is the same complexity as simply
 computing the score for the current state.

\subsection{Search Methods} We tested two common local search
methods, {\em stochastic hill climbing} and {\em taboo search.}
Stochastic hill climbing proceeds by repeatedly either changing
the the state of the variable that creates the maximum probability
change, or changing a variable at random. Figure~\ref{hill figure}
gives the algorithm explicitly.

\begin{figure}[tbp!]
\begin{center}
\fbox{
\begin{minipage}{10in}
\begin{small}
\begin{tabbing}
\underline{Given:} Probability distribution \(\pr\), evidence
\(\e\), MAP
variables \(\SS\), \\
the probability of taking a random flip $P_f$, and initial state $s_0$. \\
\underline{Compute:} An instantiation \(\s\) which
(approximately) maximizes \(\pr(\s\mid\e)\). \\ \\
Initialize current state $\s$ to $s_0$. \\
$\s_\best=\s$\\
Repe\=at m\=any\= times:\\
\>With probability $P_f$ do\\
\>\>$\s=\s^\prime$, where $\s^\prime$ is a randomly selected neighbor
of $s$.\\
\>Otherwise\\
\>\>Compute the score $Pr(\s-X,x,\e)$ for each neighbor $\s-X,x$.\\
\>\>If no neighbor has a higher score that the score for $\s$ then\\
\>\>\>$\s=\s^\prime$, where $\s^\prime$ is a randomly selected
neighbor
of $s$.\\
\>\>Else\\
\>\>\>$\s=\s^\prime$ where $\s^\prime$ is the neighbor with the highest score.\\
\>If $Pr(\s,\e)>Pr(\s_\best,\e)$ then \\
\>\> $\s_\best=\s$\\
Return $\s_\best$
\end{tabbing}
\end{small}
\end{minipage}
 }
\caption{Stochastic hill climbing algorithm. \label{hill figure}}
\end{center}
\end{figure}

Taboo search is similar to hill climbing except that the next
state is chosen as the best state that hasn't been visited
recently. Because the number of iterations is relatively small we
save all of the previous states so that at each iteration a
unique point is chosen. Pseudocode for taboo search appears in
Figure~\ref{taboo figure}.

\begin{figure}[tbp!]
\begin{center}
\fbox{
\begin{minipage}{10in}
\begin{small}
\begin{tabbing}
\underline{Given:} Probability distribution \(\pr\), evidence
\(\e\), MAP
variables \(\SS\). \\
\underline{Compute:} An instantiation \(\s\) which
(approximately) maximizes \(\pr(\s\mid\e)\). \\ \\
Initialize current state $\s$. \\
$\s_\best=\s$\\
Repe\=at m\=any \=times\\
\>Add $\s$ to $\mathit{visited}$\\
\>Compute the score $Pr(\s-X,x,\e)$ for each neighbor $\s-X,x$.\\
\>$\s=\s^\prime$ where $\s^\prime$ is a neighbor with the highest score not in $\mathit{visited}$.\\
\>If no such neighbor exists  (this rarely occurs)\\
\>\>Repeat for several times\\
\>\>\>$\s=\s^\prime$ where $\s^\prime$ is a randomly selected neighbor of $\s$.\\
\>If $Pr(\s,\e)>Pr(\s_\best,\e)$ then \\
\>\> $\s_\best=\s$\\
Return $\s_\best$
\end{tabbing}
 \end{small}
\end{minipage}
}

\caption{Taboo search.  Notice that the action taken is to
 choose the best neighbor that hasn't been visited.  This
leads to moves that decrease the score after a peak is discovered.
\label{taboo figure}}
\end{center}
\end{figure}

\subsection{Initialization} The quality of the solution returned
by a local search routine depends to a large extent on which part
of the search space it is given to explore.  We implemented
several algorithms to compare the solution quality with different
initialization schemes. Suppose that \(n\) is the number of
network variables, \(w\) is the width of a given elimination
order, and \(m\) is the number of MAP variables.

\begin{enumerate}
\item{\em Random initialization (Rand).}  For each MAP variable, we select
a value uniformly from its set of states.  This method takes
\(O(m)\) time.

\item{\em MPE based initialization (MPE).}  We compute the MPE solution
given the evidence.  Then, for each MAP variable, we set its
value to the value that the variable takes on in the MPE solution.
This method takes \(O(n \exp(w))\) time.

\item{\em Maximum likelihood initialization (ML).}  For each MAP
variable $X$, we set its value to the instance $x$ that maximizes
$Pr(x\mid\e)$.  This method takes \(O(n \exp(w))\) time.

\item{\em Sequential initialization (Seq).}  This method considers the
MAP variables \(X_1,\ldots,X_m\), choosing each time a variable
\(X_i\) that has the highest probability \(\pr(x_i\mid\e,\y)\) for
one of its values \(x_i\), where \(\y\) is the instantiation of
MAP variables considered so far. This method takes \(O(m n
\exp(w))\) time.
\end{enumerate}

\subsection{Experimental Results}\label{exact experiment section}
Two search methods ({\it Hill} and {\it Taboo}) and four
initialization methods ({\it Rand}, {\it MPE}, {\it ML}, {\it
Seq}) lead to 8 possible algorithms. Each of the initialization
methods can also be viewed as an approximation algorithm since one
can simply return the computed initialization. This leads to a
total of 12 different algorithms. We experimentally evaluated and
compared 11 of these algorithms, leaving out the algorithm
corresponding to random initialization.

\begin{table}[tb]
\begin{center}
Data Set 1 Solution Quality\\
\begin{tabular}{|r|l|l|l|l|l|} \hline
    & 0 & .125 & .250 & .375 & .5 \\ \hline
Rand-Hill & 147 & 805 & 917 & 946 & 966\\
Rand-Taboo & 181 & 969 & 985 & 993 & 995\\ \hline
ML & 526 & 497 & 676 & 766 & 817\\
ML-Hill & 920 & 947 & 989 & 993 & 997\\
ML-Taboo & 942 & 988 & 999 & 999 & \bld{1000}\\ \hline
MPE & 999 & 333 & 160 & 127 & 100\\
MPE-Hill & 999 & 875 & 923 & 952 & 973\\
MPE-Taboo & \bld{1000} & 986 & 992 & 990 & 998\\ \hline
Seq & 930 & 965 & 990 & 999 & 997\\
Seq-Hill & 941 & 971 & 992 & 999 & 997\\
Seq-Taboo & 962 & \bld{998} & \bld{1000} & \bld{1000} & \bld{1000}\\ \hline
\end{tabular}
\end{center}
\caption{The solution quality of each method for the first data set.  The number associated with each method and bias is the number of instances solved correctly out of 1000.  The best scores for each bias are shown in bold.\label{data1}}
\end{table}

\begin{table}[tb]
\begin{center}
Data Set 2 Solution Quality\\
\begin{tabular}{|r|l|l|l|l|l|} \hline
    & 0 & .125 & .250 & .375 & .5 \\ \hline
Rand-Hill & 20 & 634 & 713 & 799 & 845\\
Rand-Taboo & 20 & 851& 907 & 943 & 965\\ \hline
ML & 749 & 453 & 495 & 519 & 514\\
ML-Hill & 966 & 922 & 947 & 963 & 962\\
ML-Taboo & 973 & 960 & 986 & 987 & 990\\ \hline
MPE & 858 & 505 & 365 & 275 & 206\\
MPE-Hill & 961 & 853 & 850 & 874 & 891\\
MPE-Taboo & 978 & 952 & 962 & 977 & 980\\ \hline
Seq & 988 & 955 & 964 & 985 & 972\\
Seq-Hill & 988 & 960 & 966 & 986 & 976\\
Seq-Taboo & \bld{994} & \bld{977} & \bld{990} & \bld{994} & \bld{994}\\ \hline
\end{tabular}
\end{center}
\caption{The solution quality of each method for the second data set. The number associated with each method and bias is the number of instances solved correctly out of 1000.  The best scores for each bias are shown in bold.\label{data2}}
\end{table}

We tested the algorithms on various synthetically generated data sets as well as real world networks.
For the synthetic networks, we generated random
network structures using two generation methods (see
Appendix~\ref{random networks}).  For each structure, we
quantified the CPTs for different bias coefficients from 0
(deterministic except the roots), to .5 (values chosen uniformly)
so we could evaluate the influence of CPT quantification on the
solution quality. Each network consisted of 100 variables, with
some of the root variables chosen as the MAP
variables. If there were more than 25 root variables, we
randomly selected 25 of them for the MAP variables.  Otherwise we
used all of the root variables.  We chose root nodes for MAP
variables because typically some subset of the root nodes are the
variables of interest in diagnostic applications.  Evidence was
set by instantiating leaf nodes. Care was taken to insure that
the instantiation had a non zero probability.  Each algorithm was
allowed 150 network evaluations.\footnote{An evaluation takes
\(O(n\exp(w))\) time and space, where \(n\) is the number of
network variables and \(w\) is the width of given elimination
order.} We computed the true MAP and compared it to the solutions
found by each algorithm. Additionally, we measured the number of
network evaluations needed to find the solution each algorithm
subsequently returned, and the number of peaks discovered before
that solution was discovered.  The hill climbing method used in these
data sets is pure hill climbing with random walk restart.  That is, it
hill climbs until it reaches a peak, then randomly flips some of the values to move to a new location.

We generated 1000 random network structures for each of the
two structural generation methods.  For each random structure
generated, and each quantification method, we quantified the
network, computed the exact MAP, and applied each of the
approximation algorithms.  Tables~\ref{data1} and~\ref{data2}
show the solution quality of each of the methods by reporting the
fraction of networks that were solved correctly; that is, the
approximate answer had the same value as the exact answer.

One can draw a number of observations based on these experiments:
\begin{itemize}
\item In each case, taboo search performed slightly better
than hill climbing with random restarts.

\item The search methods were typically able to perform much
better than the initialization alone.

\item Even from a random start, the search methods were able to
find the optimal solution in the majority of the cases.

\item Overall, taboo search with sequential initialization performed
the best, but required the most network evaluations.
\end{itemize}

Table~\ref{evaluations required} contains some statistics on the
number of network evaluations (including those used for
initialization)
 needed to achieve the value that the method finally returned.
 The mean number of evaluations is quite small for all of the methods.
 Surprisingly, for the hill climbing methods, the maximum is also
 quite small.  In fact, after analyzing the results we discovered
 that the hill climbing methods never improved over the first peak
 they discovered.\footnote{It appears that the random walk used in
 restarting does not make eventually selecting a better region
 very likely when using so few search steps.  Often, when a sub
 optimal hill was encountered, the optimal hill was just 2 or 3
 moves away.  In those cases, the taboo search was usually able
 to find it (because its search was more guided), while random walking
 was not.}
 This suggests that one viable method for quick approximation is to
 simply climb to the first peak and return the result.  Taboo search
  on the other hand was able to improve on the first peak in some cases.


\begin{table}[tb]
\begin{center}
Evaluations Required\\
\begin{tabular}{|r|r|r|r|}\hline
Method & Mean & Stdev & Max \\\hline
Rand Hill & 12.5 & 2.5 & 21\\
Rand Taboo & 14.3 & 11.0 & 144 \\\hline
MPE & 1 & 0& 1\\
MPE Hill & 2.6 & 1.3 & 8\\
MPE Taboo & 4.0 & 8.3 & 137\\\hline
ML & 1 & 0 & 1\\
ML Hill & 1.6 & .74 & 4 \\
ML Taboo & 1.9 & 3.3 & 62 \\\hline
Seq & 25 & 0 & 25 \\
Seq Hill & 25.0 & .04 & 26\\
Seq Taboo & 25.0 & .9 & 45 \\\hline
\end{tabular}
\end{center}
\caption{Statistics on the number of evaluations each method
required before achieving the value it eventually returned. These
are based on the random method 2, bias .5 data set. The
statistics for the other data sets are similar.
\label{evaluations required}}
\end{table}

\begin{table}[!tbp]
\begin{center}
\begin{tabular}{|r|r|r|r|r|r|r|r|r|r|r|r|r|}
\hline
& \multicolumn{3}{|c|}{MPE} & \multicolumn{3}{|c|}{ML} & \multicolumn{3}{|c|}{Seq} & \multicolumn{3}{|c|}{MB}\\\hline
& I & H & T & I & H & T & I & H & T & 14 & 16 & 18\\\hline
Barley & 3 & 9 & 8 & 3 & 10 & 9 & 7 & 10 & 10 & 1 & 3 & 5\\\hline
Mildew & 6 & 10 & 10 & 8 & 10 & 10 & 8 & 10 & 10 & 4 & 4 & 7\\\hline
Munin2 & 6 & 10 & 10 & 10 & 10 & 10 & 10 & 10 & 10 & 4 & 5 & 7\\\hline
Munin3 & 9 & 10 & 10 & 10 & 10 & 10 & 10 & 10 & 10 & 4 & 6 & 2\\\hline
Pigs & 0 & 0 & 0 & 5 & 9 & 9 & 8 & 8 & 8 & 3 & 1 & 6\\\hline
Water & 9 & 10 & 10 & 8 & 10 & 10 & 10 & 10 & 10 & 6 & 6 & 9\\\hline
\end{tabular}
\end{center}
\caption{Number of times out of ten that each algorithm found the
  instantiation that yielded the highest score.  I, H, and T refer to
  initialization only, hill climbing, and taboo search respectively.}\label{best exact table}
\end{table}



\begin{table}[!tbp]
\begin{center}
\begin{small}
\mbox{
\begin{tabular}{|r|r|r|r|r|r|}
\multicolumn{6}{c}{{\bf Barley network results}}\\\hline
Method & Best & $> .9$ & $> .5$ & $> .01$ & $ \leq .01 $ \\\hline
MPE & 3 & 2 & 3 & 2 & 0\\
Hill & 9 & 0 & 0 & 1 & 0\\
Taboo & 8 & 0 & 1 & 1 & 0\\\hline
ML & 3 & 2 & 2 & 0 & 3\\
Hill & 10 & 0 & 0 & 0 & 0\\
Taboo & 9 & 0 & 1 & 0 & 0\\\hline
Seq& 7 & 3 & 0 & 0 & 0\\
Hill & 10 & 0 & 0 & 0 & 0\\
Taboo& 10 & 0 & 0 & 0 & 0\\\hline
MB 14 & 1 & 2 & 1 & 3 & 3\\
MB 16 & 3 & 3 & 1 & 2 & 1\\
MB 18 & 5 & 2 & 0 & 3 & 0\\\hline
%\multicolumn{6}{c}{ }
\end{tabular}
\begin{tabular}{|r|r|r|r|r|r|}
\multicolumn{6}{c}{{\bf Mildew network results}}\\\hline
Method & Best & $> .9$ & $> .5$ & $> .01$ & $ \leq .01 $ \\\hline
MPE & 6 & 1 & 3 & 0 & 0\\
Hill & 10 & 0 & 0 & 0 & 0\\
Taboo & 10 & 0 & 0 & 0 & 0\\\hline
ML & 8 & 1 & 1 & 0 & 0\\
Hill & 10 & 0 & 0 & 0 & 0\\
Taboo & 10 & 0 & 0 & 0 & 0\\\hline
Seq & 9 & 1 & 0 & 0 & 0\\
Hill & 10 & 0 & 0 & 0 & 0\\
Taboo& 10 & 0 & 0 & 0 & 0\\\hline
MB 14 & 4 & 1 & 0 & 1 & 4\\
MB 16 & 4 & 0 & 0 & 1 & 5\\
MB 18 & 7 & 0 & 1 & 0 & 2\\\hline
%\multicolumn{6}{c}{ }
\end{tabular}}

\mbox{
\begin{tabular}{|r|r|r|r|r|r|}
\multicolumn{6}{c}{{\bf Munin2 network results}}\\\hline
Method & Best & $> .9$ & $> .5$ & $> .01$ & $ \leq .01 $ \\\hline
MPE & 6 & 0 & 4 & 0 & 0\\
Hill & 10 & 0 & 0 & 0 & 0\\
Taboo & 10 & 0 & 0 & 0 & 0\\\hline
ML & 10 & 0 & 0 & 0 & 0\\
Hill & 10 & 0 & 0 & 0 & 0\\
Taboo & 10 & 0 & 0 & 0 & 0\\\hline
Seq & 10 & 0 & 0 & 0 & 0\\
Hill & 10 & 0 & 0 & 0 & 0\\
Taboo & 10 & 0 & 0 & 0 & 0\\\hline
MB 14 & 4 & 0 & 1 & 2 & 3\\
MB 16 & 5 & 0 & 1 & 2 & 2\\
MB 18 & 7 & 0 & 0 & 1 & 2\\\hline
%\multicolumn{6}{c}{ }
\end{tabular}
\begin{tabular}{|r|r|r|r|r|r|}
\multicolumn{6}{c}{{\bf Munin3 network results}}\\\hline
Method & Best & $> .9$ & $> .5$ & $> .01$ & $ \leq .01 $ \\\hline
MPE & 9 & 0 & 0 & 1 & 0\\
Hill & 10 & 0 & 0 & 0 & 0\\
Taboo & 10 & 0 & 0 & 0 & 0\\\hline
ML & 10 & 0 & 0 & 0 & 0\\
Hill & 10 & 0 & 0 & 0 & 0\\
Taboo & 10 & 0 & 0 & 0 & 0\\\hline
Seq & 10 & 0 & 0 & 0 & 0\\
Hill & 10 & 0 & 0 & 0 & 0\\
Taboo & 10 & 0 & 0 & 0 & 0\\\hline
MB 14 & 4 & 0 & 2 & 0 & 4\\
MB 16 & 6 & 0 & 1 & 0 & 3\\
MB 18 & 2 & 0 & 1 & 0 & 7\\\hline
%\multicolumn{6}{c}{ }
\end{tabular}
}

\mbox{
\begin{tabular}{|r|r|r|r|r|r|}
\multicolumn{6}{c}{{\bf Pigs network results}}\\\hline
Method & Best & $> .9$ & $> .5$ & $> .01$ & $ \leq .01 $ \\\hline
MPE & 0 & 0 & 0 & 0 & 10\\
Hill & 0 & 0 & 0 & 0 & 10\\
Taboo & 0 & 0 & 2 & 3 & 5\\\hline
ML & 5 & 1 & 3 & 1 & 0\\
Hill & 9 & 0 & 1 & 0 & 0\\
Taboo & 9 & 0 & 1 & 0 & 0\\\hline
Seq & 8 & 0 & 2 & 0 & 0\\
Hill & 8 & 0 & 2 & 0 & 0\\
Taboo & 8 & 0 & 2 & 0 & 0\\\hline
MB 14 & 3 & 0 & 3 & 4 & 0\\
MB 16 & 1 & 1 & 4 & 3 & 1\\
MB 18 & 6 & 0 & 2 & 2 & 0\\\hline
\end{tabular}
\begin{tabular}{|r|r|r|r|r|r|}
\multicolumn{6}{c}{{\bf Water network results}}\\\hline
Method & Best & $> .9$ & $> .5$ & $> .01$ & $ \leq .01 $ \\\hline
MPE & 9 & 0 & 1 & 0 & 0\\
Hill & 10 & 0 & 0 & 0 & 0\\
Taboo & 10 & 0 & 0 & 0 & 0\\\hline
ML & 8 & 1 & 1 & 0 & 0\\
Hill & 10 & 0 & 0 & 0 & 0\\
Taboo & 10 & 0 & 0 & 0 & 0\\\hline
Seq & 10 & 0 & 0 & 0 & 0\\
Hill & 10 & 0 & 0 & 0 & 0\\
Taboo & 10 & 0 & 0 & 0 & 0\\\hline
MB 14 & 6 & 1 & 2 & 0 & 1\\
MB 16 & 6 & 1 & 2 & 1 & 0\\
MB 18 & 9 & 0 & 0 & 1 & 0\\\hline
\end{tabular}}
\end{small}
\end{center}
\caption{Detailed performance measures on the real world networks.
Each column contains the number of times out of 10 that each
algorithm was able to achieve the given performance relative to the
best solution found. \label{exact detail table}}
\end{table}

We ran ten MAP queries for each real world network we tested. For
each query we randomly selected one fourth of the nodes to be the
variables of interest, and selected one fourth of the nodes to be
evidence nodes.
 The evidence values were chosen uniformly among the nonzero configurations.
  As our previous experiments demonstrated that a large number of iterations rarely helps, we reduced the number of iterations to 30.
 Also, we moved away from hill climbing with random restart to stochastic hill climbing (performing a random move with probability .35) since in our previous experiments the random restart never helped.
 Also, we ran a mini--bucket approximation algorithm (the only other MAP
 approximation algorithm we are aware of that is not subsumed by our
 technique) to compare its performance to our algorithms.  Since exact
 MAP computations on these networks is too hard for current algorithms
 to handle, we compare the algorithms based on relative performance
 only.

Table~\ref{best exact table} shows the number of times (out of
ten) that each algorithm was able to produce the highest
probability configuration discovered.  The search based methods
again performed much better than the other algorithms.  Note that
each of them outperformed the mini--bucket approximations on each
network.  Table~\ref{exact detail table} provides more specific
details about the relative performance for each network.  Each
block contains the  count of the number of times that each method
produced solutions within some range of the best found solution.
So for example, in the Barley group, in the MPE row, for the
column labeled ``$> .5$'' there is a 3, indicating that in 3 of
the 10 cases the solution found was between .5 and .9 times the
best solution found for that query.

Qualitatively, these results are very similar to those obtained
for the random networks. Again the search methods outperformed the
static initialization methods.  Note that for different networks,
different initializations perform better.  Notice also, that the
search methods significantly outperformed the mini--bucket
approximations in every network.



\section{Approximating MAP when Inference is Hard}
\label{approximation section II}

The techniques developed thus far depend on the ability to perform
exact inference.  For many networks, even inference is intractable.
In these cases, approximate inference can be substituted in order to
produce MAP approximations.

We investigate using belief propagation as the approximate
inference scheme, and local search for the optimization scheme.
Iterative belief propagation is a useful approximate inference algorithm for
approximating MAP for a number of reasons and has proven to be a
very effective and efficient approximation method for a variety of
domains.  It has the ability to approximate MPE,
posterior marginals, and probability of evidence, allowing for the
same initialization schemes as we used for exact inference.
Additionally, as we will show in section~\ref{bp
  retracted}, after a single inference call,
  the scores of neighbors in the search space
 can be computed
locally, allowing us to obtain the same linear speed up that we
obtained using a similar approach in the exact inference case.
Thus belief propagation allows all of the techniques for
approximating MAP for inference tractable networks to be applied
approximately when inference is not tractable.


\subsection{Belief Propagation Review} Belief propagation was
introduced as an exact inference method on polytrees
\cite{pearl88}. It is a message passing algorithm in which each
node in the network sends a message to its neighbors.  These
messages, along with the CPTs and the evidence can be used to
compute posterior marginals for all of the variables.  In networks
with loops, belief propagation is no longer guaranteed to be
exact, and successive iterations generally produce different
results, so belief propagation is typically run until the message
values converge.  This has been shown to provide very good
approximations for a variety of networks
\cite{mceliece95,murphy99}, and has recently received a
theoretical explanation \cite{yedidia00}.

Belief propagation works as follows. Each node $X$, has an
evidence indicator $\lambda_X$ where evidence can be entered. If
the evidence sets $X=x$, then $\lambda_X(x)=1$, and is 0
otherwise. If no evidence is set for $X$, then $\lambda_X(x)=1$
for all $x$. After evidence is entered, each node $X$ sends a
message to each of its neighbors. The  message a node $X$ with
parents $\U$ sends to child $Y$ is computed as \[M_{XY}=\alpha
\sum_{\U} \lambda_{X}\Pr(X|\U)\prod_{Z\neq Y} M_{ZX}\] where $Z$
ranges over the neighbors of $X$ and $\alpha$ is a normalizing
constant.\footnote{We use potential notation more common to join
trees than the standard descriptions of belief propagation because
we believe the many indices required in standard presentations
mask the simplicity of the algorithm.} Similarly, the message $X$
sends to a parent $U$ is \[M_{XU}=\alpha \sum_{X\U- \{U\}}
\lambda_{X}\Pr(X|\U)\prod_{Z\neq U} M_{ZX}.\]
Message passing
continues until the message values converge. The posterior of $X$
is then approximated as \[\Prp(X|\e)=\alpha \sum_{\U}
\lambda_{X}\Pr(X|\U)\prod_{Z} M_{ZX}.\]

The messages are all initialized to 1.  There are two main schemes
for ordering the messages.  In the first scheme, all of the messages
are computed simultaneously, based on the previous set of messages.
In the other scheme, messages are updated incrementally, and in two
phases, consistent with some ordering of the variables.  In the first
phase, in reverse order, each variable sends a message to its neighbors that precede it
in the order.  In the second phase, in order, each variable sends a
message to its neighbors that come after it in the order.
We implemented the second scheme since empirically it seems to
converge faster than the first scheme \cite{murphy99}.

\subsection{Approximating Neighbors' Scores}\label{bp retracted}
\begin{figure}
\centerline{\epsfxsize=75mm \leavevmode \epsffile{sct.ps}}
\caption{A scatter plot of the exact versus approximate retracted
  values of 30 variables of the Barley network. The x-axis is the
  true probability, and the y-axis the approximate probability.}\label{scatterplot
  figure}
\end{figure}

Belief propagation allows us to approximate the scores of
neighbors in the local search space efficiently, similar to what
we have done in the case of exact inference. The key as we shall
show next is to be able to compute the quantity $\Pr(x|\e-X)$ for
each variable \(X\) efficiently, as we can use this quantity to
rank the neighbors according to the desired score.

Specifically, in polytrees, the incoming messages are independent
of the value of the local CPT or any evidence entered. Hence,
leaving the evidence out of the product yields
\[\Pr(X|\e-X)=\alpha \sum_{\U} \Pr(X|\U)\prod_{Z} M_{ZX}.\]
Therefore, we can compute the above quantity for each variable
after a single belief propagation. In networks that are not
polytrees, the incoming messages are not necessarily independent
of the evidence or the local CPT, but as is done with other BP
methods, {\em we ignore that and hope that it is nearly
independent.} Empirically, the approximation seems to be quite
accurate. Figure~\ref{scatterplot figure} shows a representative
example, comparing the correspondence between the approximate and
exact retracted probabilities for 30 variables in the Barley
network. The x axis corresponds to the true retracted probability,
and the y axis to the approximation produced using belief
propagation.

Still, $\Pr(x|\e-X)$ is not quite what we want to score neighbors
in the local search space. But this quantity can be used to
compute the ratio of the neighboring score to the current score
which allows such comparisons.  Specifically, simple algebra shows
that:
\[\frac{\Pr(x,\s-X,\e)}{\Pr(\s,\e)}=\frac{\Pr(x|\s-X,\e)}{\Pr(x_{\s}|\s-X,\e)}\]
where $x_{\s}$ is the value that $X$ takes on in the current
instantiation $\s$.  Thus, we can find the neighbor with the best
score after a single belief propagation.

\subsection{Experimental Results}

For the first experiment, we consider the improvement possible
over what is typically done (MPE or ML) using it as a starting
point and hill climbing from there. For the first experiment, we
generated 100 synthetic networks with 100 variables each using the
first 
method described in Appendix~\ref{random networks} with bias
parameter $0.25$ and width parameter of 13.  We generated the
networks to be small enough that we could often compute the exact
MAP value, but large enough to make the problem challenging.  We
chose the MAP variables as the roots (typically between 20 and 25
variables), and the evidence values were chosen randomly from 10
of the leaves.  We computed the true MAP for the ones which memory
constraints (512 MB of RAM) allowed.  We computed the true
probability of the instantiations produced by the two standard
methods.  For both initialization methods we also computed the true probability of the
instantiations returned by pure hill climbing\footnote{We compare
pure and stochastic hill climbing to evaluate what can be gained
by stochastic methods.  The initial hill climb usually requires
very few evaluations, so if stochastic methods make little
difference, efficiency considerations would dictate that pure hill
climbing be used.} (i.e. only greedy steps were taken), and
stochastic hill climbing with 100 steps, where random moves were taken
with probability $p_{f}=.3$.  Of the 100
networks, we were able to compute the exact MAP in 59 of them.
Table~\ref{random table} shows the number exactly solved for each
method, as well as the worst instantiation produced, measured as
the ratio of the probabilities of the found instantiation to the
true MAP instantiation.  All of the hill climbing methods improved
significantly over their initializations in general, although for
2 of the networks, the hill climbing versions were slightly worse
than the initial value (the worst was a ratio of .835), because of
a slight mismatch in the true vs. approximate probabilities.  Over
all, the stochastic hill climbing routines outperformed the other
methods.

\begin{table}[tb]
\begin{center}
\begin{small}
\begin{tabular}{|l|ll|}
\hline
 & \# solved exactly & worst \\ \hline
MPE & 9 & .015 \\
MPE-Hill & 41 & .06 \\
MPE-SHill & 43 & .21 \\
ML & 31 & .34 \\
ML-Hill & 38 & .46 \\
ML-SHill & 42 & .72 \\
\hline
\end{tabular}
\end{small}
\end{center}
\caption{Solution quality for the random networks. Shows the
number solved exactly of the 59 for which we could compute the
true MAP value.  Worst is the ratio of the probabilities of the
found instantiation to the true MAP instantiation.  Each hill
climbing method improved significantly over the
initializations.}\label{random table}
\end{table}

In the second experiment, we generated 25 random MAP problems for
the Barley network, each with 25 randomly chosen MAP variables,
and 10 randomly chosen evidence assignments.  We use the same
parameters as in the previous experiment.  The problems were too
hard to compute the exact MAP, so we report only on the relative
improvements over the initialization methods.  Table~\ref{barley
table} summarizes the results.  Again, the stochastic hill
climbing methods were able to significantly improve the quality of
the instantiations created.
\begin{table}[tb]
\begin{center}
\begin{small}
\begin{tabular}{|l|llll|}
\hline
& min & median & mean & max \\ \hline
MPE-Hill & $1.0$ & $8.4$ & $1.3$x$10^{11}$ & $3.1$x$10^{12}$ \\
MPE-SHill & $1.0$ & $8.4$ & $1.3$x$10^{11}$ & $3.1$x$10^{12}$ \\
ML-Hill & $1.0$x$10^4$ & $3.6$x$10^7$ & $3.4$x$10^{15}$ & $8.4$x$10^{16}$\\
ML-SHill &$7.7$x$10^3$ & $3.6$x$10^7$ & $3.4$x$10^{15}$ & $8.4$x$10^{16}$\\
\hline
\end{tabular}
\end{small}
\end{center}
\caption{The statistics on the improvement over just the
initialization method for each search method on the data set
generated from the Barley network.  Improvement is measured as the
ratio of the found probability to the probability of the
initialization instantiation.}\label{barley table}
\end{table}

In the third experiment, we performed the same type of experiment
on the Pigs network.  None of the search methods were able to
improve on ML initialization.  We concluded that the problem was
too easy.  Pigs has over 400 variables, and it seemed that the
evidence didn't force enough dependence among the variables.  We
ran another experiment with Pigs, this time using 200 MAP
variables and 20 evidence values to make it more difficult.
Table~\ref{pigs table} summarizes the results.  Again, the
stochastic methods were able to improve significantly over the
initialization methods.
\begin{table}[tb]
\begin{center}
\begin{small}
\begin{tabular}{|l|llll|}
\hline
& min & median & mean & max \\ \hline
MPE-Hill & $1.0$ & $1.7$x$10^5$ & $1.5$x$10^7$ & $3.3$x$10^8$\\
MPE-SHill & $1.0$ & $2.5$x$10^5$ & $4.5$x$10^{11}$ & $1.1$x$10^{13}$\\
ML-Hill & $13.0$ & $2.0$x$10^3$ & $3.3$x$10^5$ & $4.5$x$10^6$\\
ML-SHill & $13.0$ & $1.2$x$10^4$ & $8.2$x$10^5$ & $8.2$x$10^6$\\
\hline
\end{tabular}
\end{small}
\end{center}
\caption{The statistics on the improvement over just the
initialization method alone for each search method on the data set
generated from the Pigs network.  Improvement is measured as the
ratio of the found probability to the initialization
probability.}\label{pigs table}
\end{table}


\begin{table}[!tbp]
\begin{center}
\begin{tabular}{|r|r|r|r|r|r|r|r|r|r|r|r|r|}
\hline
& \multicolumn{3}{|c|}{MPE} & \multicolumn{3}{|c|}{ML} & \multicolumn{3}{|c|}{Seq} & \multicolumn{3}{|c|}{MB}\\\hline
& I & H & T & I & H & T & I & H & T & 14 & 16 & 18\\\hline
Barley & 2 & 6 & 7 & 2 & 7 & 7 & 7 & 9 & 9 & 1 & 3 & 5\\\hline
Mildew & 5 & 1 & 10 & 5 & 0 & 10 & 9 & 0 & 10 & 4 & 4 & 7\\\hline
Munin2 & 5 & 0 & 9 & 5 & 0 & 9 & 10 & 0 & 10 & 4 & 5 & 7\\\hline
Munin3 & 8 & 0 & 1 & 8 & 0 & 1 & 10 & 0 & 1 & 4 & 6 & 2\\\hline
Pigs & 0 & 0 & 1 & 0 & 0 & 1 & 8 & 0 & 8 & 3 & 1 & 6\\\hline
Water & 9 & 0 & 0 & 9 & 0 & 0 & 10 & 0 & 0 & 6 & 6 & 9\\\hline
\end{tabular}
\end{center}
\caption{Number of times out of ten that each algorithm found the
instantiation that yielded the highest score.  The I, H and T
entries stand for initial, hill climbing, and taboo
respectively. MB stands for mini--buckets, and 14, 16, and 18 are the
width bounds.}\label{best approx table}
\end{table}



\begin{table}[!tbp]
\begin{center}
\begin{small}
\mbox{
\begin{tabular}{|r|r|r|r|r|r|}
\multicolumn{6}{c}{{\bf Barley network results}}\\\hline
Method & Best & $> .9$ & $> .5$ & $> .01$ & $ \leq .01 $ \\\hline
MPE & 2 & 3 & 1 & 3 & 1\\
Hill & 6 & 1 & 1 & 2 & 0\\
Taboo & 7 & 0 & 1 & 2 & 0\\\hline
ML & 2 & 3 & 1 & 3 & 1\\
Hill & 7 & 1 & 1 & 1 & 0\\
Taboo & 7 & 0 & 1 & 2 & 0\\\hline
SEQ  & 7 & 1 & 1 & 1 & 0\\
Hill & 9 & 0 & 1 & 0 & 0\\
Taboo & 9 & 0 & 1 & 0 & 0\\\hline
MB 14 & 1 & 2 & 1 & 3 & 3\\
MB 16 & 3 & 3 & 1 & 2 & 1\\
MB 18 & 5 & 2 & 0 & 3 & 0\\\hline
%\multicolumn{6}{c}{ }
\end{tabular}
\begin{tabular}{|r|r|r|r|r|r|}
\multicolumn{6}{c}{{\bf Mildew network results}}\\\hline
Method & Best & $> .9$ & $> .5$ & $> .01$ & $ \leq .01 $ \\\hline
MPE & 5 & 2 & 3 & 0 & 0\\
Hill & 1 & 0 & 0 & 0 & 9\\
Taboo & 10 & 0 & 0 & 0 & 0\\\hline
ML & 5 & 2 & 3 & 0 & 0\\
Hill & 0 & 0 & 0 & 0 & 10\\
Taboo & 10 & 0 & 0 & 0 & 0\\\hline
Seq & 9 & 1 & 0 & 0 & 0\\
Hill & 0 & 0 & 0 & 0 & 10\\
Taboo & 10 & 0 & 0 & 0 & 0\\\hline
MB 14 & 4 & 1 & 0 & 1 & 4\\
MB 16 & 4 & 0 & 0 & 1 & 5\\
MB 18 & 7 & 0 & 1 & 0 & 2\\\hline
%\multicolumn{6}{c}{ }
\end{tabular}}

\mbox{
\begin{tabular}{|r|r|r|r|r|r|}
\multicolumn{6}{c}{{\bf Munin2 network results}}\\\hline
Method & Best & $> .9$ & $> .5$ & $> .01$ & $ \leq .01 $ \\\hline
MPE & 5 & 0 & 3 & 1 & 1\\
Hill & 0 & 0 & 0 & 0 & 10\\
Taboo & 9 & 0 & 0 & 0 & 1\\\hline
ML & 5 & 0 & 3 & 1 & 1\\
Hill & 0 & 0 & 0 & 0 & 10\\
Taboo & 9 & 0 & 0 & 0 & 1\\\hline
Seq& 10 & 0 & 0 & 0 & 0\\
Hill & 0 & 0 & 0 & 0 & 10\\
Taboo & 10 & 0 & 0 & 0 & 0\\\hline
MB 14 & 4 & 0 & 1 & 2 & 3\\
MB 16 & 5 & 0 & 1 & 2 & 2\\
MB 18 & 7 & 0 & 0 & 1 & 2\\\hline
%\multicolumn{6}{c}{ }
\end{tabular}
\begin{tabular}{|r|r|r|r|r|r|}
\multicolumn{6}{c}{{\bf Munin3 network results}}\\\hline
Method & Best & $> .9$ & $> .5$ & $> .01$ & $ \leq .01 $ \\\hline
MPE & 8 & 0 & 0 & 1 & 1\\
Hill & 0 & 0 & 0 & 0 & 10\\
Taboo & 1 & 0 & 0 & 0 & 9\\\hline
ML  & 8 & 0 & 0 & 1 & 1\\
Hill & 0 & 0 & 0 & 0 & 10\\
Taboo & 1 & 0 & 0 & 0 & 9\\\hline
Seq & 10 & 0 & 0 & 0 & 0\\
Hill & 0 & 0 & 0 & 0 & 10\\
Taboo & 1 & 0 & 0 & 0 & 9\\\hline
MB 14 & 4 & 0 & 2 & 0 & 4\\
MB 16 & 6 & 0 & 1 & 0 & 3\\
MB 18 & 2 & 0 & 1 & 0 & 7\\\hline
%\multicolumn{6}{c}{ }
\end{tabular}
}

\mbox{
\begin{tabular}{|r|r|r|r|r|r|}
\multicolumn{6}{c}{{\bf Pigs network results}}\\\hline
Method & Best & $> .9$ & $> .5$ & $> .01$ & $ \leq .01 $ \\\hline
MPE & 0 & 0 & 0 & 0 & 10\\
Hill & 0 & 0 & 0 & 0 & 10\\
Taboo & 1 & 0 & 0 & 4 & 5\\\hline
ML & 0 & 0 & 0 & 0 & 10\\
Hill & 0 & 0 & 0 & 0 & 10\\
Taboo & 1 & 0 & 0 & 4 & 5\\\hline
Seq & 8 & 0 & 2 & 0 & 0\\
Hill & 0 & 0 & 0 & 0 & 10\\
Taboo & 8 & 0 & 2 & 0 & 0\\\hline
MB 14 & 3 & 0 & 3 & 4 & 0\\
MB 16 & 1 & 1 & 4 & 3 & 1\\
MB 18 & 6 & 0 & 2 & 2 & 0\\\hline
\end{tabular}
\begin{tabular}{|r|r|r|r|r|r|}
\multicolumn{6}{c}{{\bf Water network results}}\\\hline
Method & Best & $> .9$ & $> .5$ & $> .01$ & $ \leq .01 $ \\\hline
MPE & 9 & 0 & 1 & 0 & 0\\
Hill & 0 & 0 & 0 & 0 & 10\\
Taboo & 0 & 0 & 0 & 0 & 10\\\hline
ML  & 9 & 0 & 1 & 0 & 0\\
Hill & 0 & 0 & 0 & 0 & 10\\
Taboo & 0 & 0 & 0 & 0 & 10\\\hline
Seq & 10 & 0 & 0 & 0 & 0\\
Hill & 0 & 0 & 0 & 0 & 10\\
Taboo & 0 & 0 & 0 & 0 & 10\\\hline
MB 14 & 6 & 1 & 2 & 0 & 1\\
MB 16 & 6 & 1 & 2 & 1 & 0\\
MB 18 & 9 & 0 & 0 & 1 & 0\\\hline
\end{tabular}}
\end{small}
\end{center}
\caption{Detailed performance measures on the real world networks
  using Belief propagation approximation methods.
Each column contains the number of times out of 10 that each
algorithm was able to achieve the given performance relative to the
best solution found. \label{last approx table}}
\end{table}



We also ran these algorithms on the same queries on the real world
networks that were used in Section~\ref{exact experiment section}
to be able to compare performance between the methods.
Table~\ref{best approx table} shows how they performed and compares
their performance to the
mini--bucket algorithms.  Table~\ref{last approx table} gives a more
detailed exposition of their performance.  There are a couple of
interesting items about this data set.  One is the surprising
performance of simple sequential initialization.  Over all, it
performed the best of the approximate algorithms.  Another
interesting thing to note is that hill climbing often negatively
impacted performance.  This suggests that marginal computations
are often more accurate than probability of evidence computations.
This problem is especially acute in networks with significant
determinism.  While belief propagation believes a configuration
has significant probability, it may actually have 0 probability
because one of its constraints is violated.  These experiments
suggest that it is possible to improve on the standard approaches
used when inference is intractable (approximating MPE, or ML or
using a mini--bucket scheme) by using belief propagation to estimate
the joint, and successively moving to states with higher
approximate scores.


\section{Conclusion}
\label{sec:conclusion}

 MAP is a computationally very hard problem
which is not in general amenable to exact solution even for very
restricted classes (ex. polytrees).  Even approximation is
difficult.  Still, we can produce approximations that are much
better than those currently used by practitioners (MPE, ML)
through using approximate optimization and inference methods.  We
showed one method based on belief propagation and stochastic hill
climbing that produced significant improvements over those
methods, extending the realm for which MAP can be approximated to
networks that work well with belief propagation.


\section*{Acknowledgement}
This work has been partially supported by MURI grant N00014-00-1-0617



\appendix
\section{Proofs of Theorems}\label{proof appendix}

\subsection*{Proof of Theorem~\ref{big theorem}}\label{big theorem proof}
We want to show that MAP remains $\NPPP$-complete even when restricted
to networks of depth 2, with no
evidence, only binary variables, and parameters that are arbitrarily
close to $1/2$.  Membership in $\NPPP$ was established in
Theorem~\ref{general theorem}.  We  show hardness by providing a reduction from
\emsat.

The flow of the proof is as follows.
First, we construct a depth 2 Bayesian network from the \emsat\ problem.  Then, we show that
by asserting some evidence, we can overcome the constraint that all
of the parameters lie within $[1/2-\epsilon,1/2+\epsilon]$, and use
MAP to obtain the \emsat\ solution.  Finally,
we show that by including the evidence variables as MAP
variables instead, no evidence is needed.

\begin{figure}[tb]
\begin{center}
\centerline{\epsfxsize=60mm \leavevmode \epsffile{hardMapReduction.eps}}
\caption{The network produced using the construction in the proof of
  Theorem~\ref{big theorem} for the formula $(x_1 \land x_2) \lor
  x_3$.}\label{hard map reduction figure}
\end{center}
\end{figure}

The network is constructed as follows.  Each logical variable $x_i$ induces a
network variable $X_i$ with uniform prior.  Each operand $y_i$
induces a network variable $Y_i$ with a uniform prior.  Notice that
they are not connected, so unlike the reduction in
Theorem~\ref{general theorem}, the CPT entries do not enforce that an
operator variable take on a value that is consistent its operands with
respect to the
to the logic of the formula.  For example, the network will assign
positive probability to an ``and'' node being true, and both of its
operand variables being false.
We say that a variable $Y_i$ is consistent with the variables $\P_i$
associated with its operands, if the logical function of operator
$y_i$ yields the value of $Y_i$ on input $\p_i$.  Consistency,
instead of being enforced rigidly, is weighted by introducing $r$ weight
variables $W_{i1}...W_{ir}$ (the actual value of $r$ will be
discussed subsequently) associated with each $Y_i$.  The parents of
$W_{ij}$ are the operator
variable $Y_i$ and the variables corresponding to its operands.
The CPT  of $W_{ij}$ is defined as
\[
\Pr(W_{ij}=T | Y_i,\P_i)=\left\{
\begin{array}{ll}
\frac{1}{2}+\epsilon & \mbox{$Y_i$ is consistent with $\P_i$}\\
\frac{1}{2} & \mbox{otherwise}
\end{array}
\right.
\]

where $\P_i$ are the variables associated with the operands of $y_i$.
Finally,
as children of $Y_m$ (which corresponds to the top level operator) we add $r$ additional
binary variables
$W_{1}...W_{r}$, where
\[
\Pr(W_{j}=T | Y_m)=\left\{
\begin{array}{ll}
\frac{1}{2}+\epsilon & \mbox{if $Y_m=T$}\\
\frac{1}{2} & \mbox{otherwise}
\end{array}
\right.
\]

for the purpose of weighting states in which the formula is
satisfied.  See Figure~\ref{hard map reduction figure} for an example
network construction.


Now, consider the probability of a complete instantiation of the
variables, where all of the weight variables (which
includes both the
consistency weighting variables $W_{ij}$, and satisfiability weighting
variables $W_i$) are set
to true, which we denote as $\W=\mathbf{T}$.
\[
\Pr(\x,\y,\W=\mathbf{T})=\left(\frac{1}{2}\right)^{m+n}\left(\frac{1}{2}+\epsilon\right)^{kr}\left(\frac{1}{2}\right)^{(m-k)r}
\left(\frac{1}{2}+\epsilon\right)^{sr}\left(\frac{1}{2}\right)^{(1-s)r}
\]

where $\x$ is an instantiation of $X_1...X_n$, $\y$ is an
instantiation of $\Y_1...Y_m$, $k$ is the number of operator variables
that are consistent with their operands variables and s=1 if $Y_m=T$,
0 otherwise.  For a consistent satisfying assignment $\x\y$,

\[
\Pr(\x,\y,\W=\mathbf{T})=\left(\frac{1}{2}\right)^{m+n}\left(\frac{1}{2}+\epsilon\right)^{(m+1)r}
\]
while for an inconsistent, or unsatisfying assignment $\x\y$,
\[
\Pr(\x,\y,\W=\mathbf{T})\leq \left(\frac{1}{2}\right)^{m+n+r}(\frac{1}{2}+\epsilon)^{mr}.
\]

We want to choose $r$ such that the probability of a single
consistent satisfying instance is greater than twice the sum of all of the
probabilities of inconsistent or unsatisfying instances.  The
number of inconsistent or unsatisfying instances is bounded by
$2^{n+m}$, so we want an $r$ where
\[
\left(\frac{1}{2}\right)^{m+n}\left(\frac{1}{2}+\epsilon\right)^{(m+1)r}
> 2^{n+m+1}\left(\frac{1}{2}\right)^{m+n+r}(\frac{1}{2}+\epsilon)^{mr}.
\]
Solving for $r$ yields
\[
r > \frac{(m+n+1)}{1+\log_2\left(\frac{1}{2}+\epsilon\right)}
\]
which is linear in the size of the formula, so the size of the
reduction remains polynomially bounded.

Let
$C=\Pr(\x,\y,\W=\mathbf{T})
=\left(\frac{1}{2}\right)^{m+n}\left(\frac{1}{2}+\epsilon\right)^{(m+1)r}$
where $\x\y$ is a consistent satisfying instance.
Then, for a particular instantiation $\q$ of $X_1...X_k$,
\begin{eqnarray*}
\Pr(\q,\W=\mathbf{T}) & = & \sum_{x_{k+1},...,x_n,y_1,...,y_m}
\Pr(\q,x_{k+1},...x_n,y_1,...y_m,\W=\mathbf{T})\\
&=& \#_{\q} C + \sum_{\x\y} \Pr(\x\y,W=\mathbf{T})
\end{eqnarray*}
where $\#_{\q}$ is the number of complete variable instantiations
compatible with $\q$ that satisfies $\phi$ and $\x\y$ ranges over the
inconsistent and unsatisfying assignments compatible with
$x_{1}...x_{k}$.   Since for any instantiation of $\x$ there is only
one compatible instantiation, $\#_{\q}$ also corresponds to the
number of satisfiable instantiations of $\phi$ consistent with $\q$.
The choice of $r$ ensures that the sum of the unsatisfying or
inconsistent instantiations is less than $C/2$, and is always greater
than 0 assuming there is at least one operator (since there is some
instantiation where the operator and its operands are not consistent).  Thus $\#_q C<
\Pr(\q,\W=\mathbf{T}) < (\#_q+1/2) C$.  There are $2^{n-k}$ possible
instantiations of $X_{k+1}...X_n$, so if half or less are satisfied
then $\Pr(\q,\W=T)<(2^{n-k-1}+1/2)C$, while if more than half are
satisfied then
$\Pr(\q,\W=T)>(2^{n-k-1}+1)C$.  Thus the \dmap\ query, using MAP variables
$X_1...X_k$, evidence $\W=\mathbf{T}$, and threshold $(2^{n-k-1}+1)C$
is true if and only if the \emsat\ query is also true.

Now, notice that in every table that contains a weight variable, the
value of the configuration where it takes on true is greater or equal to
the value when it takes on false.  Thus
$\Pr(\q,\W=\mathbf{T})\ge\Pr(\q,\W=\w)$, for all $\q$ and $\w$.  It
then follows that
MAP$(X_1...X_k,\W=\mathbf{T})$=MAP$(X_1...X_k\W,\emptyset)$.
Therefore, the \dmap\ query, using MAP variables $X_1...X_k$,$\W$ and no
evidence, with threshold $(2^{n-k-1}+1)C$ is true if and only if the
\emsat\ query is also true. $\Box$

\subsection*{Proof of Theorem~\ref{big time polytree approximation theorem}}\label{big time polytree theorem proof}

\begin{figure}[tb]
\begin{center}
\centerline{\epsfxsize=60mm \leavevmode \epsffile{app2.eps}}
\caption{The network used in the reduction of Theorem~\ref{big time polytree
    approximation theorem}.}\label{polytree approx figure}
\end{center}
\end{figure}


As part of the proof of the theorem, we will use the following lemma.
\begin{lemma}\label{bound lemma}
For all $x\ge 1$, $4 x+\frac{1}{2}>\frac{1}{\ln(1+\frac{1}{4x})}$.
\end{lemma}
{\bf Proof:}
First, we show that $f(x)=\ln(1+\frac{1}{4 x})-\frac{1}{4
  x+\frac{1}{2}}$ is monotonically decreasing for $x\ge 1$.
\begin{eqnarray*}
\frac{df}{dx} & = & \frac{-1}{4 x^2+x}+\frac{4}{(4 x+\frac{1}{2})^2} \\
& = & \frac{16 x^2+4 x-(4 x+\frac{1}{2})^2}{(4 x+\frac{1}{2})^2 (4 x^2 +x)} \\
& = & \frac{-1}{4 (4 x+\frac{1}{2})^2 (4 x^2 +x)}
\end{eqnarray*}
which is always negative for $x\ge 1$, and hence $f(x)$ is
monotonically decreasing.

Now, since $f(x)$ is monotonically decreasing, and
$\displaystyle \lim_{x \rightarrow \infty} f(x)=0$, $f(x)$ must be
strictly positive.  Thus, for all $x\ge 1$,
$\ln(1+\frac{1}{4 x}) > \frac{1}{4 x+\frac{1}{2}}$ which implies
$4 x+\frac{1}{2}>\frac{1}{\ln(1+\frac{1}{4 x})}$. $\Box$


The basic idea of the proof is to show by repeating the construction of
Theorem~\ref{polytree theorem} a polynomial number of times,  that if we can approximate MAP on
polytrees within relative error $2^{\mbox{size}^\epsilon}$ for any
$\epsilon\in[0,1)$, where the size of the network is parameterized by
  the number of conditional probability parameters, then we can solve SAT in polynomial time.

Given a SAT problem instance with $n$ variables, and $m$ clauses, we
create a Bayesian network by replicating
the construction from Theorem~\ref{polytree theorem} q times, and
connecting them to form a polytree.   Specifically, to each copy
$i$ of the construction, we add a variable $B^i$ (we use
superscripts to denote variables associated with a particular copy
of the construction), with parents $S^i_n$, and if $i>1$, parent
$B^{i-1}$ (see Figure~\ref{polytree approx figure}). The
conditional probability of $B^i$ is uniform for all parent
instantiations.  We choose $q$ to satisfy
$(1+\frac{1}{4m})^q>2^{\mbox{size}^\epsilon}$.  We now show that $q$
can be chosen so that the network size remains polynomial in the size of
the logical formula.
The resulting network has $q(2n+2)$ variables, and each conditional
probability table
has at most $2(m+1)^2$ parameters, so the total size of the reduction is
bounded by $q(m+1)^2(4n+4)$.  Replacing size with the size bound
places the constraint  \[\left(1+\frac{1}{4m}\right)^q>2^{\left(q(m+1)^2(4n+4)\right)^\epsilon}\]
on $q$. Since $0\leq\epsilon <1$, solving for $q$ yields
\begin{eqnarray*}
\left(1+\frac{1}{4m}\right)^q & > & 2^{\left(q(m+1)^2(4n+4)\right)^\epsilon} \\
q \ln\left(1+\frac{1}{4m}\right)& > & q^\epsilon (m+1)^{2\epsilon}(4n+4)^\epsilon \ln
2\\
q^{1-\epsilon} & > & \frac{(m+1)^{2\epsilon}(4n+4)^\epsilon\ln 2}{\ln\left(1+\frac{1}{4m}\right)}\\
q & > & \left(\frac{(m+1)^{2\epsilon}(4n+4)^\epsilon\ln 2}{\ln\left(1+\frac{1}{4m}\right)}\right)^{\frac{1}{1-\epsilon}}
\end{eqnarray*}
Now, from Lemma~\ref{bound lemma}, $4m+1/2>1/\ln(1+1/4m)$, so
substitution yields a stronger bound,
\begin{eqnarray*}
q & > & \left(\left(4m+\frac{1}{4}\right)(m+1)^{2\epsilon}(4n+4)^\epsilon\ln
2\right)^{\frac{1}{1-\epsilon}}\\
\end{eqnarray*} which is polynomially bounded.
Thus the network can be constructed in time polynomial in the size of
the formula.

Then, for a particular instantiation $\x$ of all $X$ variables
$X_1^1...X_n^q$, and evidence $\s$ asserting $S^i_n=0$ for each $i$,
the probability is
\begin{eqnarray*}
\Pr(\x,\s) &=&\prod_i \Pr(\x^i,S^i_n=0) \\
&=&\prod_i  \frac{\mbox{\# clauses satisfied by $\x^i$}}{m 2^n}
\end{eqnarray*}
 because each subnetwork is
independent.  Thus the solution $M$ to MAP over $X_1^1...X_n^q$ with
evidence $S^1_n=....=S^q_n=0$ is $\left(\frac{k}{m 2^n}\right)^q$ where $k$ is the
maximum number of clauses that can be simultaneously satisfied in the
original SAT problem.  If the problem is satisfiable then $k=m$, and
so the approximate solution $M'$ obeys

\[M' \geq \frac{M}{2^{\mbox{size}^\epsilon}} > \left(\frac{4 m}{4
  m+1}\right)^q\left(\frac{m}{m 2^n}\right)^q >
\left(\frac{m-\frac{1}{4}}{m}\right)^q
\left(\frac{1}{2^n}\right)^q=\left(\frac{m-\frac{1}{4}}{m2^n}\right)^q
\]

  On the other hand, if it
isn't satisfiable, then $k\leq m-1$, so

\[ M' \leq 2^{\mbox{size}^\epsilon} M < \left(\frac{4m+1}{4m}\right)^q
\left(\frac{m-1}{m2^n}\right)^q =
\left(\frac{(4m+1)(m-1)}{4m}\right)^q\left(\frac{1}{m2^n}\right)^q < \left(\frac{m-\frac{3}{4}}{m2^n}\right)^q\]
The upper bound of $M'$ if the SAT problem is unsatisfiable is bounded
below the lower bound of $M'$ if it is satisfiable.  Because the
network construction and the bound tests can be accomplished in polynomial time,
if the MAP problem itself can be approximated within a factor of
$2^{\mbox{size}^\epsilon}$ in polynomial time then SAT
can be decided in polynomial time.

\section{Generating Random Networks}\label{random networks}
We generated several types of networks to perform our experiments.
We used two methods for generating the structure, and a single parametric
method for generating the quantification.

\subsection{Generating the Network Structure}
The first method is parameterized by the number of variables
 $N$ and the connectivity $c$.  This method tends to produce
 structures with widths that are close to $c$. \citeA{darwicheAIJ-RC}
 provides an algorithmic description.

The second method is parameterized by the number of variables $N$,
 and the probability $p$ of an edge being present. We generate an
 ordered list of $N$ variables, and add an edge between variables
 $X$ and $Y$ with probability $p$.  The edges added are directed
 toward the variable that appears later in the order.




%For the experiments in Figures~\ref{quality1} and~\ref{quality2},
%we used $N=100$, $c$ between 6 and 12, and $p=.025$. These
%numbers were chosen so that the MAP width would be small enough
%that we could compute the exact value to measure the solution
%quality.


\subsection{Quantifying the Dependencies}
The quantification method is parameterized by a bias parameter
$b$. The values of the CPTs for the roots were chosen uniformly.
The values for the rest of the nodes were based on a bias,
where one of the values $v$ was chosen uniformly in $[0,b)$,
and the other as $1-v$.  For example, for $b=.1$, each non root
variable given its parents has one value in $[0,.1)$, and the
other in $(.9,1]$.  Special cases $b=0$, and $b=.5$ produce
deterministic, and uniformly random quantifications respectively.

\bibliography{references}
\bibliographystyle{theapa}

\end{document}
