# What type of inference is planning?

Miguel Lázaro-Gredilla Li Yang Ku Kevin P. Murphy Dileep George
Google Deepmind
{lazarogredilla, liyangku, kpmurphy<EMAIL_ADDRESS>

###### Abstract

Multiple types of inference are available for probabilistic graphical models,
e.g., marginal, maximum-a-posteriori, and even marginal maximum-a-posteriori.
Which one do researchers mean when they talk about “planning as inference”?
There is no consistency in the literature, different types are used, and their
ability to do planning is further entangled with specific approximations or
additional constraints. In this work we use the variational framework to show
that all commonly used types of inference correspond to different weightings
of the entropy terms in the variational problem, and that planning corresponds
_exactly_ to a _different_ set of weights. This means that all the tricks of
variational inference are readily applicable to planning. We develop an
analogue of loopy belief propagation that allows us to perform approximate
planning in factored-state Markov decisions processes without incurring
intractability due to the exponentially large state space. The variational
perspective shows that the previous types of inference for planning are only
adequate in environments with low stochasticity, and allows us to characterize
each type by its own merits, disentangling the type of inference from the
additional approximations that its practical use requires. We validate these
results empirically on synthetic MDPs and tasks posed in the International
Planning Competition.

##  1 Introduction

There are many kinds of probabilistic inference that are used, such as
marginal, maximum-a-posteriori (MAP), or marginal MAP (MMAP). In this work we
show that planning is a distinct type of inference, and that _under stochastic
dynamics_ does not correspond exactly with any of the above methods.
Furthermore, we show how to rank the above methods in terms of quality as it
pertains to planning.

Our approach is based on variational inference. The variational perspective
that we adopt allows direct comparison between different inference types, and
allows us to develop analogues of existing approximate inference algorithms
for this “planning inference” task. Given a "flat" Markov Decision Process
(MDP), respresented as the probabilistic graphical model (PGM) in Fig 1[Left],
we show that planning inference provides the same (exact) results as value
iteration. And using an analogue of loopy belief propagation (LBP), we show
how to apply approximate planning inference to the factored MDP in Fig.
1[Right], which has an exponentially large state space, and for which exact
solutions are no longer tractable. Under moderate stochasticity in the
dynamics, we show that this approximate planning inference is superior to
other more established types of inference.

Figure 1: Factor graphs: [Left] Standard MDP [Right] Factored MDP with sparse
factor connectivity.

##  2 Background

###  2.1 Markov Decision processes (MDPs)

A finite-horizon Markov decision process (MDP) is a tuple $({\cal X},{\cal
A},p(x_{1}),{\cal P},{\cal R},T)$, where ${\cal X}$ is the state space, ${\cal
A}$ the action space, $p(x_{1})$ the starting state distribution, ${\cal P}$
the transition probabilities $P(x_{t+1}|x_{t},a_{t})$ (i.e., the dynamics of
the process), $R_{t}(x_{t},a_{t},x_{t+1})$ the reward for transitioning from
$x_{t}$ to $x_{t+1}$ under action $a_{t}$ at timestep $t$, and $T$ is the
horizon. For simplicity of exposition we will consider discrete states and
actions, but analogous results apply in the continuous case. Solving an MDP
corresponds to finding a policy $\pi_{t}(a_{t}|x_{t})$ that maximizes the
expectation of some function $g(\cdot)$ of the reward
$\mathbb{E}[g(R_{1}(x_{1},a_{1},x_{2}),\ldots,R_{T-1}(x_{T-1},a_{T-1},x_{T}))]$,
typically its sum. The optimal policy can be different at each time step,
i.e., non-stationary. We do not use a discount factor, but it can be trivially
included in the reward function, since it is also non-stationary. Policies,
states and actions at all timesteps will be noted
$\bm{\pi}\equiv\\{\pi_{t}\\}_{t=1}^{T-1},\bm{x}\equiv\\{x_{t}\\}_{t=1}^{T}$
and $\bm{a}\equiv\\{a_{t}\\}_{t=1}^{T-1}$.

In a state-factored MDP (factored MDP for short), each state $x_{t}$
factorizes into $N_{e}$ r.v. or _entities_ $\\{x^{(i)}_{t}\\}_{i=1}^{N_{e}}$,
each with cardinality $N_{s}$, so it has an exponentially large state space of
size $N_{s}^{N_{e}}$. Transitions factorize as
$P(x_{t+1}|x_{t},a_{t})=\prod_{i=1}^{N_{e}}P(x^{(i)}_{t+1}|{\bm{x}}^{\text{pa}(i)}_{t},a_{t})$,
where ${\bm{x}}^{\text{pa}(i)}_{t}$ is a small subset of $x_{t}$ on which
$x^{(i)}_{t+1}$ depends, allowing for a tabular definition of the transition
without exponential cost ($\text{pa}(i)$ stands in for “parents of $i$” ).
Similarly, the reward factorizes as
$R_{t}(x_{t})=\sum_{i=N_{e}+1}^{N_{e}+N_{r}}R_{t}({\bm{x}}^{\text{pa}(i)}_{t})$.

###  2.2 Variational inference

Let’s consider running different _types_ of inference in the
unnormalized111Just like in undirected graphical models, the different types
of inference are well-defined despite lack of normalization. Table 1 shows the
precise definitions including two reasonable definitions of marginalization.
factor graph shown in Fig. 1[Left]. Marginal inference would compute the sum
over all the ${\bm{x}},{\bm{a}}$ configurations (the partition function); MAP
inference would compute the maximizing configuration (the posterior mode); and
marginal MAP (MMAP) inference (see [23]) would find the assignment of $\bm{a}$
that maximizes the summation over $\bm{x}$ conditional on that $\bm{a}$.
Although marginal, MAP, and MMAP inference are distinct, a lot is shared. They
all target a “quantity of interest” (e.g., partition function, maximum
probability state, best conditional partition function); they all produce a
distribution as a result of inference (e.g., delta at the mode, full
posterior, best conditional posterior); naïve computation requires a number of
summations and/or maximizations exponential in the number of variables; and
finally, they can all be represented as a variational inference (VI)
problem.222 We assume familiarity with variational inference, see [18] for an
introduction.  The VI representation naturally leads to efficient message-
passing solutions and approximate inference algorithms, as we show below.

For a general factor graph $f({\bm{x}},{\bm{a}})$, the VI problem is
$\max_{q({\bm{x}},{\bm{a}})}\langle\log
f({\bm{x}},{\bm{a}})\rangle_{q({\bm{x}},{\bm{a}})}+H^{\text{type}}_{q}({\bm{x}},{\bm{a}})$
where $q({\bm{x}},{\bm{a}})$ is an arbitrary variational distribution over the
variables of the factor graph, the term $-\langle\log
f({\bm{x}},{\bm{a}})\rangle_{q({\bm{x}},{\bm{a}})}$ is known as the _energy
term_ , and $H^{\text{type}}_{q}({\bm{x}},{\bm{a}})$ is a particular entropy
choice that will determine the inference type. It is a standard result, (e.g.,
[18]) that using the Shannon entropy $H_{q}({\bm{x}},{\bm{a}})$ results in
marginal inference. Setting it to zero (aka the zero-temperature limit)
corresponds to MAP inference [37, 34, 36, 20, 26]. Setting it to the
conditional entropy of
$H_{q}({\bm{x}}|{\bm{a}})=H_{q}({\bm{x}},{\bm{a}})-H_{q}({\bm{a}})$ results in
MMAP inference [23]. Note that the only difference across VI problems is the
_weighting_ of the entropy terms.

Despite the similarities, the computational complexity of these types of
inference can differ significantly, even for tree-structured graphs. In
particular, the entropy term for MMAP is not concave, and inference is NP-hard
[23], whereas for marginal and MAP inference the entropy term is concave (the
energy term is always linear), and inference is polynomial.

##  3 Methods

Table 1: Different types of inference from a variational perspective, including a proper “planning” inference type. They all share the same energy term $E_{\lambda}(\bm{q})$ defined in Eq. (3), and differ only in the entropy term. The closed-form expressions provide the optimal value of the bound, but are generally intractable. The general tractability of the bound maximization for MDPs is marked in the Tr column. All bounds are monotonically related, and listed in descending order, see Section 4.1. Type of |  $\Big{\downarrow}$ | Closed form for quant. of interest |  $\Big{\downarrow}$ | Entropy term for variational bound | Tr
---|---|---|---|---|---
inference | $F_{\lambda}=\max_{\bm{q}}F_{\lambda}({\bm{q}})$ | $F_{\lambda}({\bm{q}})=\frac{1}{\lambda}(-E_{\lambda}(\bm{q})+H(\bm{q}))$ |
Marginal333 Two types of marginal inference are included for precision. “Marginal” refers to marginal inference directly applied on the same exact factor graph as the other types of inference. Because the factor graph lacks a prior over $\bm{a}$, it is not a properly normalized joint distribution. Adding a uniform prior over the actions resolves this and results in “Marginal${}^{\text{U}}$”. Both are a constant apart and are of independent interest. | $\frac{1}{\lambda}\log\sum_{{\bm{x}},{\bm{a}}}P({\bm{x}}|{\bm{a}})e^{\lambda R({\bm{x}},{\bm{a}})}$ | $H_{\bm{q}}(x_{1})+\sum_{t=1}^{T-1}H_{\bm{q}}(x_{t+1},a_{t}|x_{t})$ | ✓
Planning | $\frac{1}{\lambda}\max_{\bm{\pi}}\log\mathbb{E}_{\bm{\pi}}\big{[}e^{\lambda R({\bm{x}},{\bm{a}})}\big{]}$ | $H_{\bm{q}}(x_{1})+\sum_{t=1}^{T-1}H_{\bm{q}}(x_{t+1}|x_{t},a_{t})$ | ✓
M. MAP | $\frac{1}{\lambda}\max_{\bm{a}}\log\sum_{\bm{x}}P({\bm{x}}|{\bm{a}})e^{\lambda R({\bm{x}},{\bm{a}})}$ | $H_{\bm{q}}(x_{1})+\sum_{t=1}^{T-1}H_{\bm{q}}(x_{t+1},a_{t}|x_{t})-H_{\bm{q}}(a_{t})$ | ✗
MAP | $\frac{1}{\lambda}\max_{{\bm{x}},{\bm{a}}}\log P({\bm{x}}|{\bm{a}})e^{\lambda R({\bm{x}},{\bm{a}})}$ | $0$ | ✓
Marginal${}^{\text{U}}$ | $\frac{1}{\lambda}\log\sum_{{\bm{x}},{\bm{a}}}P({\bm{x}}|{\bm{a}})\frac{1}{N_{a}}e^{\lambda R({\bm{x}},{\bm{a}})}$ | $H_{\bm{q}}(x_{1})+\sum_{t=1}^{T-1}H_{\bm{q}}(x_{t+1},a_{t}|x_{t})-\log N_{a}^{\frac{T-1}{\lambda}}$ | ✓

In this section we introduce our VI framework, which we will use to derive a
novel linear programming formulation of planning problems, a novel value
belief propagation (VBP) algorithm and a novel closed form (sampling-free)
approach to determinization.

###  3.1 VI for standard MDPs

The main quantity of interest in this paper is the _best exponential utility_
, which we will refer to simply as the _utility_. Given an MDP with horizon
$T$ and a risk parameter $\lambda>0$, the utility is defined as

$\displaystyle F_{\lambda}^{\text{planning}}$
$\displaystyle=\frac{1}{\lambda}\log\max_{\bm{\pi}}\sum_{\bm{a},\bm{x}}\exp\Big{(}\lambda\sum_{t=1}^{T-1}R_{t}(x_{t},a_{t},x_{t+1})\Big{)}P(x_{1})\prod_{t=1}^{T-1}P(x_{t+1}|x_{t},a_{t})\pi_{t}(a_{t}|x_{t})$
(1)
$\displaystyle=\frac{1}{\lambda}\max_{\bm{\pi}}\log\mathbb{E}_{\bm{\pi}}\Big{[}\exp(\lambda
R({\bm{x}},{\bm{a}}))\Big{]}.$

Observe that $F_{\lambda\rightarrow
0^{+}}^{\text{planning}}=\max_{\bm{\pi}}\mathbb{E}_{\bm{\pi}}\Big{[}R({\bm{x}},{\bm{a}})\Big{]}$
[25] corresponds to the expected reward under the optimal policy, so we are
tackling a strictly more general case. On the one hand, this definition is
more suitable for a factor graph representation, since all the terms (both
dynamics and exponentiated reward) are combined in a multiplicative way. On
the other hand, this more general formulation allows us to trade-off between
risk-neutral ($\lambda\rightarrow 0^{+}$) and risk seeking ($\lambda>0$)
policies, see [25, 12, 33] for more details.

We can turn this new quantity of interest, the utility, into the solution of a
VI problem on the factor graph of Fig. 1[Left]. Crucially, the factor graph
includes the known dynamics and rewards terms, but _not_ any policy term,
since the policy is the outcome of inference.

######  Theorem 1 (Variational formulation of planning).

Given known dynamics $P(x_{t+1}|x_{t},a_{t})$, an initial distribution
$P(x_{1})$ and reward functions $R_{t}(x_{t},a_{t},x_{t+1})$, the best
exponential utility $F^{\text{planning}}_{\lambda}$ from Eq. (1) can be
expressed as the result of a concave variational optimization problem

$F_{\lambda}^{\text{planning}}=\max_{\bm{q}}F^{\text{planning}}_{\lambda}({\bm{q}});~{}~{}~{}~{}~{}~{}F^{\text{planning}}_{\lambda}({\bm{q}})=\frac{1}{\lambda}(-E_{\lambda}(\bm{q})+H^{\text{planning}}(\bm{q}))$
(2)

with energy $E_{\lambda}(\bm{q})$ and entropy $H^{\text{planning}}(\bm{q})$
terms

$\displaystyle E_{\lambda}(\bm{q})$ $\displaystyle=-\langle\log
P(x_{1})\rangle_{q(x_{1})}-\sum_{t=1}^{T-1}\langle\log
P(x_{t+1}|x_{t},a_{t})+\lambda
R_{t}(x_{t},a_{t},x_{t+1})\rangle_{q(x_{t+1},x_{t},a_{t})}$ (3) $\displaystyle
H^{\text{planning}}(\bm{q})$
$\displaystyle=H_{\bm{q}}(x_{1})+\sum_{t=1}^{T-1}H_{\bm{q}}(x_{t+1}|x_{t},a_{t})$
(4)

where ${\bm{q}}\equiv q({\bm{x}},{\bm{a}})$ is an arbitrary distribution over
the space of states and actions.

Proof is in Appendix A. This entropy thus corresponds to “planning inference”.
The optimal policy can be obtained from the variational posterior
$q(a_{t}|x_{t})$. Table 1 lists the types of inference problems and their
associated entropies. As we will discuss in Section 4, they have a monotonic
ordering.

The VI problem Eq. (2) reduces to the standard one when $\lambda=1$, and
extends VI in a meaningful way in the presence of rewards, regardless of the
type of inference used: rewards interact in an additive way when
$\lambda\rightarrow 0^{+}$, rather than the default multiplicative (or more
precisely, exponentiated summation) interaction of $\lambda=1$. Furthermore,
it turns out that it is possible to take the $\lambda\rightarrow 0^{+}$ limit
exactly, to obtain the dual LP formulation of an MDP [30].

######  Corollary 1.1 (Additive limit).

In the limit $\lambda\rightarrow 0^{+}$, the concave problem Eq. (2) becomes
the following linear program (LP):

$\displaystyle F_{\lambda\rightarrow 0^{+}}^{\text{planning}}$
$\displaystyle=\max_{\bm{q}}F^{\text{planning}}_{\lambda\rightarrow
0^{+}}({\bm{q}})=\max_{\\{q(x_{t},a_{t})\\}_{t=1}^{T-1}}\sum_{t=1}^{T-1}\langle
R_{t}(x_{t},a_{t},x_{t+1})\rangle_{P(x_{t+1}|x_{t},a_{t})q(x_{t},a_{t})}$
$\displaystyle\text{s.t.
}q(x_{1})=P(x_{1});~{}~{}~{}~{}~{}\sum_{a_{t+1}}q(x_{t+1},a_{t+1})=\sum_{x_{t},a_{t}}P(x_{t+1}|x_{t},a_{t})q(x_{t},a_{t})~{}~{}\forall_{t};$
$\displaystyle~{}~{}~{}~{}~{}~{}q(x_{t})=\sum_{a_{t}}q(x_{t},a_{t})~{}~{}\forall_{t};~{}~{}~{}~{}~{}q(x_{t},a_{t})\geq
0~{}\forall t,$

which corresponds to the maximum expected reward $F_{\lambda\rightarrow
0^{+}}^{\text{planning}}=\max_{\bm{\pi}}\mathbb{E}_{\bm{\pi}}[R({\bm{x}},{\bm{a}})]$.

See Appendix B for proof. This $\lambda$-based formulation will prove useful
next in the factored case.

###  3.2 VI LP and VBP for factored MDPs

Factored MDPs (e.g., Fig. 1[Right]) are loopy factor graphs with an
exponentially large state space, so the previous approaches cannot be applied
directly. An effective approximate marginal inference approach for this type
of problem is loopy belief propagation (LBP). Since planning is now seen as a
type of inference, we can create an analogue to LBP which we call value belief
propagation (VBP).

Following LBP, we make two approximations to Eq. (2) to make it tractable.

First, we replace the variational distribution ${\bm{q}}$ with pseudo-
marginals $\bm{\tilde{q}}$. Eqs. (3) and (4) never access to the full joint
$q({\bm{x}},{\bm{a}})$, but only the local marginals of each factor. Pseudo-
marginals are a collection of local distributions, consistent at each
variable, but not necessarily marginals of any distribution. Just like
${\bm{q}}$ is defined in a convex region called the _marginal polytope_
$\cal{M}$, the pseudo-marginals $\bm{\tilde{q}}$ are defined in an outer
convex region called the _local polytope_ $\cal{L}$ that contains $\cal{M}$
[38].

Second, we replace the entropy $H^{\text{planning}}({\bm{q}})$ with its
_Bethe_ approximation

$H_{\text{Bethe}}^{\text{planning}}({\bm{\tilde{q}}})=\sum_{i=1}^{N_{e}}H_{\bm{\tilde{q}}}(x^{(i)}_{1})+\sum_{t=1}^{T-1}\Big{(}H_{\text{Bethe}}(\tilde{q}({\bm{x}}_{t}))+\sum_{i=1}^{N_{e}}H_{\bm{\tilde{q}}}(x^{(i)}_{t+1}|{\bm{x}}^{\text{pa}(i)}_{t},a_{t})\Big{)},$

where
$H_{\text{Bethe}}(\tilde{q}({\bm{x}}_{t}))=\sum_{j=1}^{N_{e}+N_{r}}\Big{(}H_{\bm{\tilde{q}}}({\bm{x}}^{\text{pa}(j)}_{t})-\sum_{k\in\text{pa}(j)}H_{\bm{\tilde{q}}}(x^{(k)}_{t})\Big{)}$,
see Appendix C for details. This approximation is tractable as long as the
factors (transition and rewards) are tractable, typically by connecting to a
small number of variables. Note that the policy (which would connect to all
the state variables in a time slice and introduce exponential cost) is not a
factor in the graph.

The term $H_{\text{Bethe}}(q({\bm{x}}_{t}))$ is key. It is non-concave in
general and can be interpreted as a (negative) mutual information correcting
the discrepancy between the entropy of a collection of variables and the sum
of their entropies. It makes the optimization problem harder, but also more
accurate.

In the case of MDPs $H_{\text{Bethe}}(q({\bm{x}}_{t}))=0$ and
${\bm{q}}={\bm{\tilde{q}}}$, so we recover Eq. (2). For general factored MDPs
this is not true. We can still choose to ignore this correction to obtain a
concave bound

$H_{\text{concave}}^{\text{planning}}({\bm{\tilde{q}}})=\sum_{i=1}^{N_{e}}H_{\bm{\tilde{q}}}(x^{(i)}_{1})+\sum_{t=1}^{T-1}\sum_{i=1}^{N_{e}}H_{\bm{\tilde{q}}}(x^{(i)}_{t+1}|{\bm{x}}^{\text{pa}(i)}_{t},a_{t})\geq
H_{\text{Bethe}}^{\text{planning}}({\bm{\tilde{q}}}).$ (5)

Then, the planning Bethe approximation of the variational bound and its
concave upper bound are

$\displaystyle\tilde{F}_{\lambda}^{\text{planning}}=\max_{\bm{\tilde{q}}}\tilde{F}_{\lambda}^{\text{planning}}(\bm{\tilde{q}})=\max_{\bm{\tilde{q}}}\frac{1}{\lambda}(-E(\bm{\tilde{q}})+H_{\text{Bethe}}^{\text{planning}}({\bm{\tilde{q}}}))~{}~{}~{}s.t.~{}~{}\bm{\tilde{q}}\in\cal{L}$
(6)
$\displaystyle\hat{F}_{\lambda}^{\text{planning}}=\max_{\bm{\tilde{q}}}\hat{F}_{\lambda}^{\text{planning}}(\bm{\tilde{q}})=\max_{\bm{\tilde{q}}}\frac{1}{\lambda}(-E(\bm{\tilde{q}})+H_{\text{concave}}^{\text{planning}}({\bm{\tilde{q}}}))~{}~{}~{}s.t.~{}~{}\bm{\tilde{q}}\in\cal{L}.$
(7)

We see that
$\hat{F}_{\lambda}^{\text{planning}}\geq\tilde{F}_{\lambda}^{\text{planning}}$
and $\hat{F}_{\lambda}^{\text{planning}}\geq F_{\lambda}^{\text{planning}}$.
The former is trivial given the negative term removed. The latter follows from
(a) switching the optimization domain from $\cal{M}$ to
$\cal{L}\supseteq\cal{M}$, which can only increase the value of the bound, and
(b) Eq. (5) corresponds to Eq. (4), but with joint entropies over entities
replaced with sums of the entropies, which is an upper bound. The fact that
$\hat{F}_{\lambda}^{\text{planning}}(\bm{\tilde{q}})$ is concave and upper
bounds the exact utility has two advantages: it can be computed without local
minima problems, and it is an _admissible heuristic_ of the original utility,
meaning that it can be used as a heuristic for algorithms that emit a
certificate of optimality or infeasibility.

######  Lemma (Additive limit for factored MDPs).

In the limit $\lambda\rightarrow 0^{+}$, the concave problem Eq. (7) becomes
the following VI LP:

$\displaystyle\hat{F}_{\lambda\rightarrow 0^{+}}^{\text{planning}}$
$\displaystyle=\max_{\bm{\tilde{q}}}\hat{F}^{\text{planning}}_{\lambda\rightarrow
0^{+}}({\bm{\tilde{q}}})=\max_{\bm{\tilde{q}}}\sum_{t=1}^{T}\sum_{j=N_{e}}^{N_{e}+N_{r}}\langle
R_{t}({\bm{x}}^{\text{pa}(j)}_{t})\rangle_{\tilde{q}({\bm{x}}^{\text{pa}(j)}_{t})}$
$\displaystyle\text{s.t.
}~{}~{}\tilde{q}(x^{(i)}_{1})=P(x^{(i)}_{1})~{}~{}\forall
i;~{}~{}~{}~{}~{}\bm{\tilde{q}}\in\cal{L}$
$\displaystyle~{}~{}~{}~{}~{}~{}\tilde{q}(x^{(i)}_{t+1})=\sum_{{\bm{x}}^{\text{pa}(i)}_{t},a_{t}}P(x^{(i)}_{t+1}|{\bm{x}}^{\text{pa}(i)}_{t},a_{t})\tilde{q}({\bm{x}}^{\text{pa}(i)}_{t},a_{t})~{}~{}\forall
t,x_{t+1},a_{t},1\leq i\leq N_{e}$

which upper bounds the maximum expected reward $\hat{F}_{\lambda\rightarrow
0^{+}}^{\text{planning}}\geq F_{\lambda\rightarrow
0^{+}}^{\text{planning}}=\max_{\bm{\pi}}\mathbb{E}_{\bm{\pi}}[R({\bm{x}},{\bm{a}})]$.
Alternatively, the same expression can be obtained from Corollary 1.1 by
relaxing the marginal polytope into the local polytope. Since it is a
relaxation, the upper bounding is trivial.

To the best of our knowledge, this is a novel VI LP formulation and it can be
used to tractably (over) estimate the optimal expected reward in factored
MDPs. Similarities with [19, 14] are only surface level, see Section 5.

$\hat{F}_{\lambda}^{\text{planning}}(\bm{\tilde{q}})$ from Eq. (7) can be
maximized with a conic solver (or an LP solver if $\lambda\rightarrow 0^{+}$).
The non-concave $\tilde{F}_{\lambda}^{\text{planning}}(\bm{\tilde{q}})$ from
Eq. (6) is more challenging. Conveniently,
$\tilde{F}_{\lambda}^{\text{planning}}(\bm{\tilde{q}})$ looks just like the
Bethe free energy that motivates LBP, but with a different weighting of the
local entropy terms.

Multiple works consider modifying the entropy weighting in LBP, usually with
the aim of “concavifying” the overall entropy term and developing convergent
alternatives to LBP. In particular, [15] provides fixed-point message updates
for arbitrary entropy weights. For the specific weighting of
$H_{\text{Bethe}}^{\text{planning}}({\bm{\tilde{q}}})$ the message updates
approach a singularity, which we will avoid by using
$(1-\epsilon)H_{\text{Bethe}}^{\text{planning}}({\bm{\tilde{q}}}))+\epsilon
H_{\text{Bethe}}^{\text{marginal}}({\bm{\tilde{q}}})$. The resulting message
passing algorithm updates are well-defined for any $\epsilon>0$, interpolate
between “planning inference” and marginal inference, and can get arbitrarily
close to the former by making $\epsilon\rightarrow 0^{+}$. The use of
$\epsilon$ is not just a mathematical trick, but can be practically used to
continuously anneal from LBP into VBP, which can improve convergence. See [23]
for an analogous technique with the same purpose.

VBP inherits many of the properties of LBP: the message updates are not
guaranteed to converge, but if they do, they do so at a fixed point of Eq. 6.
Convergence can be improved by the use of damping and annealing. The precise
message updates for the general case are provided in Appendix D.

Computation associated with VBP and VI LP scale as ${\cal
O}(T\sum_{i=1}^{N_{e}+N_{r}}N_{a}N_{s}^{\text{pa(i)}+1})$, as expected.

###  3.3 VBP for standard MDPs

It is instructive to look at the VBP updates for a standard, non-factored MDP.
In this case, it is possible to take the limit $\epsilon\rightarrow 0^{+}$ and
get well-defined updates. For $\lambda=1$ and a single reward at $T$

Backward updates: $\displaystyle
m_{\text{b}}(x_{T})=e^{R_{T}(x_{T})};~{}~{}~{}m_{\text{b}}(x_{t})=\max_{a_{t}}Q(x_{t},a_{t});$
Forward updates: $\displaystyle
m_{\text{f}}(x_{1})=P(x_{1});~{}~{}~{}m_{\text{f}}(x_{t+1})=\sum_{x_{t},a_{t}}p(x_{t+1}|x_{t},a_{t}=\operatorname*{argmax}_{a^{\prime}_{t}}Q(x_{t},a^{\prime}_{t}))m_{\text{f}}(x_{t})$
Posterior: $\displaystyle q(x_{t+1},x_{t},a)\propto
m_{\text{b}}(x_{t+1})p(x_{t+1}|x_{t},a_{t}=\operatorname*{argmax}_{a^{\prime}_{t}}Q(x_{t},a^{\prime}_{t}))m_{\text{f}}(x_{t})$

where
$Q(x_{t},a_{t})=\sum_{x_{t+1}}m_{\text{b}}(x_{t+1})p(x_{t+1}|x_{t},a_{t})$.
Iterating these updates converges in a single backward and forward pass to the
global optimum. The backward messages correspond to the value function (hence
the name VBP), and the familiar intermediate quantity $Q(x_{t},a_{t})$ matches
the Q-function. The forward messages correspond to occupancy probabilities
under the optimal policy.

###  3.4 Determinization in hindsight

The previous presentation implies that all VI tricks are now applicable to
planning. As an example, we can show that _determinization_ [41], a technique
from the planning literature to extend deterministic planning algorithms to
stochastic domains, and that is usually computed via sampling, can be obtained
in closed form by leveraging the variational formulation.

To be more precise, we can compute $F_{\lambda=1}^{\text{det. planning}}$ (for
MDPs) and $\hat{F}_{\lambda=1}^{\text{det. planning}}$ (for factored MDPs)
exactly as a single concave optimization problem (avoiding sampling) when the
inner deterministic planning problem is solved with a standard LP MAP
relaxation [34]. Additionally we can prove that for factored MDPs
$\hat{F}_{\lambda=1}^{\text{planning}}\leq\hat{F}_{\lambda=1}^{\text{det.
planning}}$ (i.e., the superiority of the bound introduced here wrt
determinization in the case of factored MDPs). See Appendix E for further
details.

##  4 The different types of inference and their adequacy for planning

###  4.1 Ranking inference types for planning

As we show in Section 5, the term “planning as inference” has been used in the
literature to refer to different inference types, none of which corresponds
with the exact “planning inference” from this work, which is exact. Table 1
gives us direct insight as to how the different types compare when used for
planning. First, from inspection of the entropy terms it is trivial to see
that, in the exact case, and for a given variational posterior ${\bm{q}}$,
$F^{\text{marginal}^{\text{U}}}_{\lambda}({\bm{q}})\leq
F^{\text{MAP}}_{\lambda}({\bm{q}})\leq F^{\text{MMAP}}_{\lambda}({\bm{q}})\leq
F^{\textbf{planning}}_{\lambda}({\bm{q}})\leq
F^{\text{marginal}}_{\lambda}({\bm{q}})$. This in turn means that for the
optimal posterior of each type of inference we have

$F^{\text{marginal}^{\text{U}}}_{\lambda}\leq F^{\text{MAP}}_{\lambda}\leq
F^{\text{MMAP}}_{\lambda}\leq F^{\textbf{planning}}_{\lambda}\leq
F^{\text{marginal}}_{\lambda}.$ (8)

Since the tightness of a lower bound is an indication of its quality, it
follows that MMAP inference is no worse and potentially better than all other
common types of inference. However, as noted in Section 2.2, MMAP inference is
particularly hard, even in trees, meaning that in the case of a non-factorial
MDP like the one in Fig. 1[Left], the computation of
$F^{\text{MMAP}}_{\lambda}$ is intractable, even though all the other
quantities, including the one of interest, $F^{\text{planning}}_{\lambda}$,
are exactly computable. What we can tractably compute is the lower bound
$F^{\text{MMAP}}_{\lambda}({\bm{q}})\leq F^{\text{MMAP}}_{\lambda}$ and try to
maximize it wrt $\bm{q}$, but without guarantees of finding the optimal value.
Thus, among the common inference typles, MMAP is the best, but it can be
subject to bad local minima problems. This seems more acceptable in the
factored MDP case, but it is disappointing that the problem remains in the
standard (single entity) MDP case.

###  4.2 The stochasticity of the dynamics is key

The energy term Eq. (3), which is common to all inference methods, contains
subterms $\langle\log P(x_{t+1}|x_{t},a_{t})\rangle_{q(x_{t+1},x_{t},a_{t})}$.
When dynamics are deterministic, this forces the optimal variational
conditional to be $q(x_{t+1}|x_{t},a_{t})=P(x_{t+1}|x_{t},a_{t})$, since any
other choice would make those subterms, and therefore the bound, $-\infty$.
This affects the tightness of the quantities of interest

$F^{\text{marginal}^{\text{U}}}_{\lambda}\leq
F^{\text{MAP}}_{\lambda}=F^{\text{MMAP}}_{\lambda}=F^{\textbf{planning}}_{\lambda}\leq
F^{\text{marginal}}_{\lambda},$

and _justifies the use of MAP and MMAP inference as planning when dynamics are
deterministic_. When using approximate inference, if dynamics are close to
deterministic, it might make more sense to choose the type of inference based
on the quality of the approximation, rather than its tightness. If dynamics
are stochastic, the suboptimality of MMAP can be explained as a _lack of
reactivity to the environment_. Indeed, if we reduce the planning problem to a
non-reactive policy $\pi(a_{t}|x_{t})=\pi(a_{t})$ we recover MMAP inference as
optimal. We test this experimentally in Section 6 and further expand on it in
Appendix F.2. MAP has the same problem, but additionally lacks integration
over observation sequences (“trajectories”). Even with deterministic dynamics,
marginal inference might not produce good utility estimates, but its action
posterior will be proportional to the reward of the action sequence, so if we
additionally assume $\exp(\lambda R({\bm{x}}))\in\\{0,1\\}$ (i.e., pure
planning where we want to attain any of a subset of states), it will produce
optimal planning _choices_.

##  5 Related work

As stated, the meaning of “planning as inference” is uneven across the
literature. [35] introduce the policy in the MDP factor graph and maximize the
likelihood wrt to its parameters using EM. This is an exact approach, although
it is more appropriate to say that it is planning as _learning_ rather than a
type of inference. [22] is a well-known reference that uses _MAP inference_ to
plan. It only considers the multiplicative setting $\lambda=1$. Acknowledging
the limitations of that proposal under stochastic dynamics, it suggests a
structured posterior where $q(x_{t+1}|x_{t},a_{t})=P(x_{t+1}|x_{t},a_{t})$.
Eq. (2) shows that with this choice all entropy terms cancel with energy ones,
resulting in a simpler problem. Still, the rigid posterior prevents the
optimal solution to the original problem from being found in general. In
contrast, our entropy term allows for an unrestricted variational posterior.
[10] introduces ARollout, which can be seen as a single-forward-pass LBP
approximation to _marginal inference_ , and is applicable to factored MDPs. In
the follow-up works [9, 8] the authors develop conformant SOGBOFA, which
approximates _marginal-MAP_ inference by using ARollout in an inner loop and
gradient descent to optimize over the actions in an outer loop. A number of
refinements are added for superior performance. This is a strong baseline and
was the runner-up in the international probabilistic planning competition
(IPPC) 2018, which agrees with our analysis from Section 4. Most works, such
as [2] choose MAP inference for planning.

Two frameworks [27, 39] have been recently introduced to analyze planning from
a message-passing perspective. The former analyzes six update rules and their
qualitative effect on the plans; the latter focuses on disentangling the
direction in which messages are propagated from the _approximation_ type,
although keeps the approximation type entangled with the _inference_ type,
with different approximations tacitly using different inference types.

Using a rigorous variational framework [7, 6] tackle graph-based MDPs, similar
to factored MDPs, but with a factorized action space: multiple actions are
taken at each time step, each locally affecting a single entity. This locality
results in additional efficiencies, so direct application to a state-only-
factored MDPs would still result in exponential cost.

LP formulations for the solution non-stationary, finite-horizon MDPs have
received significant attention over the last decade (e.g., [21, 5, 1, 4]), but
they lack a variational perspective and do not generalize easily to handle
state-factored MDPs. The LPs in [19, 14, 24] on the other hand do handle
factored MDPs and have a closer connection to our work. The problem setup is
slightly different, infinite-horizon MDPs with a stationary policy in their
case vs our finite horizon MDPs, which allows us to use _local_ non-stationary
policies. More importantly, their solution approximately decomposes the value
function into a _linear combination_ of user provided basis functions. The
quantity and quality of these will determine the quality of the approximation,
introducing additional design decisions. In VI planning we use the structure
of the graph to determine the approximation, which is multiplicative over the
value of backward messages, rather than additive.

##  6 Empirical validation

Figure 2: Performance of different types of inference on factored MDPs as a
function of their level of stochasticity (normalized entropy). [Left]
Estimation error of the best utility. Lower is better. [Right] Advantage of
the next action prescribed by a method vs. optimal planning. Higher is better.

##### Synthetic MDPs

We generate $5,000$ synthetic factored MDPs structured as in Fig. 1[Right]
with random dynamics, all-or-nothing reward at the last timestep, and
controlled normalized entropies, defined as
$H_{\text{MDP}}=\frac{\sum_{x_{t+1},x_{t},a_{t}}H(p(x_{t+1}|x_{t},a_{t}))}{N_{e}N_{a}N_{s}\log
N_{s}}\in[0,1]$. See Appendix F.1 for more details.

They are purposefully small so that we can compute
$F^{\text{marginal}^{\text{U}}}$, $F^{\text{MAP}}$, $F^{\text{MMAP}}$,
$F^{\text{planning}}$ exactly, even though they are intractable in general. We
also compute $\tilde{F}^{\text{marginal}}$ (tractable), and
$\tilde{F}^{\text{MMAP}}$ (tractable bound, generally intractable
optimization), which correspond to ARollout [10] and optimal SOGBOFA-LC* [8],
respectively. Finally, we include VBP (tractable, imperfect optimization of
$\tilde{F}^{\text{planning}}({\bm{q}})$) and the tractable VI LP
$\hat{F}_{\lambda=0}^{\text{planning}}$ and VI CVX
$\hat{F}_{\lambda=1}^{\text{planning}}$.

Fig. 2 shows the effect of stochasticity in the estimation of the utility
[Left] and the next best action [Right]. For high stochasticity, VBP, even if
approximate, dominates all other types of inference. The concave upper bounds
$\hat{F}_{\lambda=0}^{\text{planning}}$ and
$\hat{F}_{\lambda=1}^{\text{planning}}$ also improve over exact MMAP but not
as much as VBP. For low stochasticity, exact MAP and MMAP dominate VBP. The
(intractably optimal) SOGBOFA-LC* remains close to the exact MMAP. These
results agree well with the theory and observations laid out in Section 4. We
see good correlation between the accuracy of the utility estimation and the
quality of the planning choices ([Left] and [Right] panels). ARollout and
exact marginal seem to be an exception to this; this is explained in Section
4.2: for pure planning problems, with low stochasticity both methods are a
constant away from the right utility and make good choices.

##### Reactivity avoidance

We craft a multi-entity MDP in which the agent controls the level of
reactivity (see Section 4.2) needed to solve the environment, but is penalized
for lower ones. VBP keeps the reactivity at a maximum, to achieve a reward of
1. SOGBOFA-LC* “aware” that it cannot plan reactively (despite replanning),
takes step to reduce it, getting a reward of 0.33. See Appendix F.2.

Figure 3: Cumulative rewards on 6 problem domains from the ICAPS 2011 IPPC. A
small horizontal jitter was introduced in all data points for visual clarity.
Each cumulative reward is averaged over 30 simulations per instance. Datasets
are ordered from left to right and top to bottom by increasing normalized
entropy levels. Only the last two have a significant stochasticity level >5%.

##### International probabilistic planning competition tasks (IPPC)

We follow [39] and compare on the same tasks and with the same methods. We use
the 6 different domains from IPPC2011, each with 10 instances (factored MDPs)
of increasing difficulty, with given dynamics and (stationary) rewards,
40-step episodes, and mildly stochastic dynamics. As baselines, we use MFVI-
Bwd [39], ARollout ($\tilde{F}^{\mathrm{marginal}^{\text{U}}}_{\lambda=0}$,
see [10]), and SOGBOFA-LC ($\tilde{F}^{\mathrm{MMAP}}_{\lambda=0}$, see [8]).
We provide details about these competing methods in Appendix F.3. From our
proposed variational framework, we use VI LP
($\hat{F}^{\mathrm{planning}}_{\lambda=0}$), and VBP444Setting $\lambda=0$ can
result in degeneracy problems, so we use a small value instead, see Appendix
F.3. ($\tilde{F}^{\mathrm{planning}}_{\lambda\approx 0}$).

Fig. 3 shows the average cumulative reward for all domains and methods. Four
domains are highly deterministic ($H_{\text{MDP}}<0.05$), but planning
inference manages to be competitive wrt the best baselines. The other two are
Game of Life and SysAdmin, which have an average $H_{\text{MDP}}$ of 0.18 and
0.23 respectively. We notice a significant advantage of our proposals wrt the
most sophisticated method, SOGBOFA-LC
($\tilde{F}^{\mathrm{MMAP}}_{\lambda=0}$). This is consistent with our
expectation of MMAP degrading with increased stochasticity (see Section 4).
Elevators is well known for its challenging rewards [10] and the only one for
which ARollout performs noticeably worse. In this domain, VBP manages to match
or exceed SOGBOFA-LC on most instances. Overall, we observe that VBP is more
consistent across varying stochasticities, matching the performance of the
best method for each dataset. None of these domains reach the larger
stochasticity levels shown in Fig. 2 where VBP dominates. VI LP also performs
generally well, although not as well as VBP due to the missing mutual
information mentioned in Section 3.2. See Appendix F.3 for details.

##  7 Discussion

The variational framework offers a powerful tool to analyze and understand how
different existing types of inference approximate planning, the key role of
stochasticity, what the ideal type of inference for planning is, and how to
design new approximations. We anticipate the introduced VI perspective to lead
to novel algorithms for planning and to further the understanding of existing
methods.

## References

  * [1] Eitan Altman  “Constrained Markov decision processes”  Routledge, 2021
  * [2] Hagai Attias  “Planning by probabilistic inference”  In _International workshop on artificial intelligence and statistics_ , 2003, pp. 9–16  PMLR
  * [3] Elias Bareinboim and Judea Pearl  “Causal inference and the data-fusion problem”  In _Proceedings of the National Academy of Sciences_ 113.27 National Acad Sciences, 2016, pp. 7345–7352
  * [4] Sanjay P Bhat, Veeraruna Kavitha and Nandyala Hemachandra  “Finite-Horizon Constrained MDPs With Both Additive And Multiplicative Utilities”  In _arXiv preprint arXiv:2303.07834_ , 2023
  * [5] Arnab Bhattacharya and Jeffrey P Kharoufeh  “Linear programming formulation for non-stationary, finite-horizon Markov decision process models”  In _Operations Research Letters_ 45.6 Elsevier, 2017, pp. 570–574
  * [6] Feng Chen, Qiang Cheng, Jianwu Dong, Zhaofei Yu, Guojun Wang and Wenli Xu  “Efficient approximate linear programming for factored MDPs”  In _International Journal of Approximate Reasoning_ 63 Elsevier, 2015, pp. 101–121
  * [7] Qiang Cheng, Qiang Liu, Feng Chen and Alexander T Ihler  “Variational planning for graph-based MDPs”  In _Advances in Neural Information Processing Systems_ 26, 2013
  * [8] Hao Cui, Thomas Keller and Roni Khardon  “Stochastic planning with lifted symbolic trajectory optimization”  In _Proceedings of the International Conference on Automated Planning and Scheduling_ 29, 2019, pp. 119–127
  * [9] Hao Cui and Roni Khardon  “Online Symbolic Gradient-Based Optimization for Factored Action MDPs.”  In _IJCAI_ , 2016, pp. 3075–3081
  * [10] Hao Cui, Roni Khardon, Alan Fern and Prasad Tadepalli  “Factored MCTS for large scale stochastic planning”  In _Proceedings of the AAAI Conference on Artificial Intelligence_ 29.1, 2015
  * [11] Hao Jackson Cui, Radu Marinescu and Roni Khardon  “From stochastic planning to marginal MAP”  In _Advances in Neural Information Processing Systems_ 31, 2018
  * [12] Hans Föllmer and Alexander Schied  “Stochastic finance: an introduction in discrete time”  Walter de Gruyter, 2011
  * [13] Hector Geffner and Blai Bonet  “A concise introduction to models and methods for automated planning”  Springer Nature, 2022
  * [14] Carlos Guestrin, Daphne Koller, Ronald Parr and Shobha Venkataraman  “Efficient solution algorithms for factored MDPs”  In _Journal of Artificial Intelligence Research_ 19, 2003, pp. 399–468
  * [15] Tamir Hazan and Amnon Shashua  “Norm-product belief propagation: Primal-dual message-passing for approximate inference”  In _IEEE Transactions on Information Theory_ 56.12 IEEE, 2010, pp. 6294–6316
  * [16] Malte Helmert  “The fast downward planning system”  In _Journal of Artificial Intelligence Research_ 26, 2006, pp. 191–246
  * [17] Jörg Hoffmann and Bernhard Nebel  “The FF planning system: Fast plan generation through heuristic search”  In _Journal of Artificial Intelligence Research_ 14, 2001, pp. 253–302
  * [18] Michael I Jordan, Zoubin Ghahramani, Tommi S Jaakkola and Lawrence K Saul  “An introduction to variational methods for graphical models”  In _Machine learning_ 37 Springer, 1999, pp. 183–233
  * [19] Daphne Koller and Ronald Parr  “Computing factored value functions for policies in structured MDPs”  In _IJCAI_ 99, 1999, pp. 1332–1339
  * [20] Vladimir Kolmogorov  “Convergent tree-reweighted message passing for energy minimization”  In _International Workshop on Artificial Intelligence and Statistics_ , 2005, pp. 182–189  PMLR
  * [21] Atul Kumar, Veeraruna Kavitha and Nandyala Hemachandra  “Finite horizon risk sensitive MDP and linear programming”  In _2015 54th IEEE Conference on Decision and Control (CDC)_ , 2015, pp. 7826–7831  IEEE
  * [22] Sergey Levine  “Reinforcement learning and control as probabilistic inference: Tutorial and review”  In _arXiv preprint arXiv:1805.00909_ , 2018
  * [23] Qiang Liu and Alexander Ihler  “Variational algorithms for marginal MAP”  In _Journal of Machine Learning Research_ 14, 2013, pp. 3165–3200
  * [24] Alan Malek, Yasin Abbasi-Yadkori and Peter Bartlett  “Linear programming for large-scale Markov decision problems”  In _International conference on machine learning_ , 2014, pp. 496–504  PMLR
  * [25] Alexandre Marthe, Aurélien Garivier and Claire Vernade  “Beyond Average Reward in Markov Decision Processes”  In _Sixteenth European Workshop on Reinforcement Learning_ , 2023
  * [26] André FT Martins, Mário AT Figueiredo, Pedro MQ Aguiar, Noah A Smith and Eric P Xing  “Ad3: Alternating directions dual decomposition for map inference in graphical models”  In _The Journal of Machine Learning Research_ 16.1 JMLR. org, 2015, pp. 495–545
  * [27] Francesco AN Palmieri, Krishna R Pattipati, Giovanni Di Gennaro, Giovanni Fioretti, Francesco Verolla and Amedeo Buonanno  “A unifying view of estimation and control using belief propagation with application to path planning”  In _IEEE Access_ 10 IEEE, 2022, pp. 15193–15216
  * [28] Judea Pearl  “The causal foundations of structural equation modeling”  In _Handbook of structural equation modeling_ , 2012, pp. 68–91
  * [29] Laurent Perron and Vincent Furnon  “OR-Tools”, 2024  Google  URL: https://developers.google.com/optimization/
  * [30] Martin L Puterman  “Markov decision processes: discrete stochastic dynamic programming”  John Wiley & Sons, 2014
  * [31] Paul K Rubenstein, Sebastian Weichwald, Stephan Bongers, Joris M Mooij, Dominik Janzing, Moritz Grosse-Wentrup and Bernhard Schölkopf  “Causal consistency of structural equation models”  In _arXiv preprint arXiv:1707.00819_ , 2017
  * [32] Scott Sanner  “ICAPS 2011 international probabilistic planning competition”, 2011  URL: https://users.cecs.anu.edu.au/~ssanner/IPPC_2011/
  * [33] Yun Shen, Michael J Tobia, Tobias Sommer and Klaus Obermayer  “Risk-sensitive reinforcement learning”  In _Neural computation_ 26.7 MIT Press, 2014, pp. 1298–1328
  * [34] David Sontag, Amir Globerson and Tommi Jaakkola  “Introduction to dual decomposition for inference”, 2011
  * [35] Marc Toussaint and Amos Storkey  “Probabilistic inference for solving discrete and continuous state Markov Decision Processes”  In _Proceedings of the 23rd international conference on Machine learning_ , 2006, pp. 945–952
  * [36] Martin J Wainwright, Tommi S Jaakkola and Alan S Willsky  “MAP estimation via agreement on trees: message-passing and linear programming”  In _IEEE transactions on information theory_ 51.11 IEEE, 2005, pp. 3697–3717
  * [37] Yair Weiss, Chen Yanover and Talya Meltzer  “MAP estimation, linear programming and belief propagation with convex free energies”  In _arXiv preprint arXiv:1206.5286_ , 2012
  * [38] Adrian Weller, Kui Tang, Tony Jebara and David A Sontag  “Understanding the Bethe approximation: When and how can it go wrong?”  In _UAI_ , 2014, pp. 868–877
  * [39] Zhennan Wu and Roni Khardon  “Approximate Inference for Stochastic Planning in Factored Spaces”  In _International Conference on Probabilistic Graphical Models_ , 2022, pp. 433–444  PMLR
  * [40] Jonathan S Yedidia, William T Freeman and Yair Weiss  “Constructing free-energy approximations and generalized belief propagation algorithms”  In _IEEE Transactions on information theory_ 51.7 IEEE, 2005, pp. 2282–2312
  * [41] Sung Wook Yoon, Alan Fern, Robert Givan and Subbarao Kambhampati  “Probabilistic planning via determinization in hindsight.”  In _AAAI_ , 2008, pp. 1010–1016

##  Appendix A Proof of the variational formulation of planning

In this proof we will use two identities, the first is the variational
identity [18]:

$\log\sum_{{\bm{x}},{\bm{a}}}f({\bm{x}},{\bm{a}})=\max_{q({\bm{x}},{\bm{a}})}\langle\log
f({\bm{x}},{\bm{a}})\rangle_{q({\bm{x}},{\bm{a}})}+H(q({\bm{x}},{\bm{a}}))$
(9)

and the second is

$\max_{\pi}\sum_{\bm{a}}q({\bm{a}})\log\pi({\bm{a}})=\max_{\pi({\bm{a}})}-H(q({\bm{a}}))-\operatorname{KL}(q({\bm{a}})||\pi({\bm{a}}))=\sum_{\bm{a}}q({\bm{a}})\log
q({\bm{a}}),$ (10)

which follows because the term $\pi({\bm{a}})$ is an arbitrary distribution
that can make the KL divergence exactly zero (by choosing
$\pi({\bm{a}})=q({\bm{a}})$). With this we can proceed to the main proof:

$\displaystyle F_{\lambda}^{\text{planning}}=$
$\displaystyle\frac{1}{\lambda}\max_{\bm{\pi}}\log\sum_{\bm{a},\bm{x}}\exp\Big{(}\lambda\sum_{t=1}^{T-1}R_{t}(x_{t},a_{t},x_{t+1})\Big{)}P(x_{1})\prod_{t=1}^{T-1}P(x_{t+1}|x_{t},a_{t})\pi_{t}(a_{t}|x_{t})$
$\displaystyle\overset{\text{Eq.~{}\eqref{eq:app_var}}}{=}$
$\displaystyle\frac{1}{\lambda}\max_{\bm{\pi},\bm{q}}\Big{\langle}\log\Big{(}P(x_{1})\prod_{t=1}^{T-1}\exp(\lambda
R_{t}(x_{t},a_{t},x_{t+1}))P(x_{t+1}|x_{t},a_{t})\pi_{t}(a_{t}|x_{t})\Big{)}\Big{\rangle}_{q(\bm{x},\bm{a})}+H(q(\bm{x},\bm{a}))$
$\displaystyle=$
$\displaystyle\frac{1}{\lambda}\max_{\bm{\pi},\bm{q}}\Big{(}-E_{\lambda}({\bm{q}})+H(q(\bm{x},\bm{a}))+\sum_{t=1}^{T-1}\langle\log\pi_{t}(a_{t}|x_{t})\rangle_{q(x_{t},a_{t})}\Big{)}$
$\displaystyle=$
$\displaystyle\frac{1}{\lambda}\max_{\bm{q}}\Big{(}-E_{\lambda}({\bm{q}})+H_{\bm{q}}(x_{1})+\sum_{t=1}^{T-1}H_{\bm{q}}(x_{t+1},a_{t}|x_{t})+\langle\max_{\pi_{t}}\langle\log\pi_{t}(a_{t}|x_{t})\rangle_{q(a_{t}|x_{t})}\rangle_{q(x_{t})}\Big{)}$
$\displaystyle\overset{\text{Eq.~{}\eqref{eq:max_lik}}}{=}$
$\displaystyle\frac{1}{\lambda}\max_{\bm{q}}\Big{(}-E_{\lambda}({\bm{q}})+H_{\bm{q}}(x_{1})+\sum_{t=1}^{T-1}H_{\bm{q}}(x_{t+1},a_{t}|x_{t})+\langle\langle\log
q(a_{t}|x_{t})\rangle_{q(a_{t}|x_{t})}\rangle_{q(x_{t})}\Big{)}$
$\displaystyle=$
$\displaystyle\frac{1}{\lambda}\max_{\bm{q}}\Big{(}-E_{\lambda}({\bm{q}})+H_{\bm{q}}(x_{1})+\sum_{t=1}^{T-1}H_{\bm{q}}(x_{t+1},a_{t}|x_{t})-H_{\bm{q}}(a_{t}|x_{t})\Big{)}$
$\displaystyle=$
$\displaystyle\frac{1}{\lambda}\max_{\bm{q}}\Big{(}-E_{\lambda}({\bm{q}})+H^{\text{planning}}({\bm{q}})\Big{)}=\max_{\bm{q}}F_{\lambda}^{\text{planning}}(\bm{q}),$

where we additionally see that we have chosen
$\pi_{t}(a_{t}|x_{t})=q(a_{t}|x_{t})$, therefore, the optimal policy will be
$\pi_{t}(a_{t}|x_{t})=q^{*}(a_{t}|x_{t})$ where $q^{*}$ is the value of $q$
that maximizes $F_{\lambda}^{\text{planning}}(\bm{q})$.

##  Appendix B Proof that for $\lambda\rightarrow 0^{+}$ the variational bound
turns into an LP

Let us rewrite the concave optimization functional Eq. (2) (derived in
Appendix A) by combining the energy and entropy terms into KL divergences:

$\displaystyle F_{\lambda}^{\text{planning}}(\bm{q})=$
$\displaystyle\frac{1}{\lambda}\Big{(}-E_{\lambda}({\bm{q}})+H^{\text{planning}}({\bm{q}})\Big{)}=\frac{1}{\lambda}\Big{(}-\operatorname{KL}(q(x_{1})||P(x_{1}))$
$\displaystyle-\sum_{t=1}^{T-1}\langle\operatorname{KL}(q(x_{t+1}|x_{t},a_{t})||P(x_{t+1}|x_{t},a_{t}))\rangle_{q(x_{t},a_{t})}+\langle\lambda
R_{t}(x_{t},a_{t},x_{t+1})\rangle_{q(x_{t+1},x_{t},a_{t})}\Big{)}$
$\displaystyle=$ $\displaystyle\sum_{t=1}^{T-1}\langle
R_{t}(x_{t},a_{t},x_{t+1})\rangle_{q(x_{t+1},x_{t},a_{t})}$
$\displaystyle-\frac{\operatorname{KL}(q(x_{1})||P(x_{1}))+\sum_{t=1}^{T-1}\langle\operatorname{KL}(q(x_{t+1}|x_{t},a_{t})||P(x_{t+1}|x_{t},a_{t}))\rangle_{q(x_{t},a_{t})}}{\lambda}$

It is clear that if any of the KL terms are larger than 0, then
$\lim_{\lambda\rightarrow 0^{+}}F_{\lambda}^{\text{planning}}(\bm{q})=-\infty$
(for bounded rewards), whereas if all KL terms are 0, the limit is finite.
That means that to maximize the bound wrt $\bm{q}$ in the $\lambda\rightarrow
0^{+}$ limit, we will choose $q(x_{1})=P(x_{1})$ and
$q(x_{t+1}|x_{t},a_{t})=P(x_{t+1}|x_{t},a_{t})$, which allows to remove the KL
terms and results in the (constrained) LP optimization of Corollary 1.1 (which
also explicitly includes the marginalization constraints $\bm{q}\in\cal{M}$).

##  Appendix C Derivation of the planning Bethe entropy for factored MDPs

The standard Bethe entropy of a factored MDP, such as the one in Fig. 4, is
[40]:

$\displaystyle H^{\text{marginal}}_{\text{Bethe}}({\bm{\tilde{q}}})=$
$\displaystyle\sum_{i=1}^{N_{e}}H_{\bm{\tilde{q}}}(x^{(i)}_{1})+\sum_{t=1}^{T-1}\Big{(}H_{\text{Bethe}}(q({\bm{x}}_{t},a_{t}))$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\sum_{i=1}^{N_{e}}H_{\bm{\tilde{q}}}(x^{(i)}_{t+1},{\bm{x}}^{\text{pa}(i)}_{t},a_{t})-H_{\bm{\tilde{q}}}({\bm{x}}^{\text{pa}(i)}_{t},a_{t})\Big{)},$

where we use ${\bm{x}}^{\text{pa}(j)}_{t}$ to refer to the “parent” variables.
To simplify notation, when $j\in 1,\ldots,N_{e}$,
${\bm{x}}^{\text{pa}(j)}_{t}$ is the collection variables on which the
distribution of $x^{(j)}_{t+1}$ depends according to the dynamics model. When
$j\in N_{e}+1,\ldots,N_{e}+N_{r}$, ${\bm{x}}^{\text{pa}(j)}_{t}$ is the
collection of variables on which the $j-N_{e}$-th reward depends. See also
Fig. 4 for clarification on the notation of parents of dynamics and rewards.

The Bethe entropy above was defined, for conciseness, in terms of the Bethe
entropy of a portion of a time slice

$\displaystyle H_{\text{Bethe}}(q({\bm{x}}_{t},a_{t}))=$
$\displaystyle(1-N_{e})H_{\bm{\tilde{q}}}(a_{t})+\sum_{i=1}^{N_{e}}H_{\bm{\tilde{q}}}({\bm{x}}^{\text{pa}(i)}_{t},a_{t})+\sum_{r=N_{e}+1}^{N_{e}+N_{r}}H_{\bm{\tilde{q}}}({\bm{x}}^{\text{pa}(r)}_{t})-\sum_{j=1}^{N_{e}+N_{r}}\sum_{k\in\text{pa}(j)}H_{\bm{\tilde{q}}}(x^{(k)}_{t}).$

Finally, the Bethe entropy of the states of factored MDP is

$H_{\text{Bethe}}(q({\bm{x}}_{t}))=\sum_{j=1}^{N_{e}+N_{r}}\Big{(}H_{\bm{\tilde{q}}}({\bm{x}}^{\text{pa}(j)}_{t})-\sum_{k\in\text{pa}(j)}H_{\bm{\tilde{q}}}(x^{(k)}_{t})\Big{)}.$

Note that all the Bethe entropy definitions are linear combinations of
standard Shannon entropies defined over subsets of variables in the factor
graph. The subsets are defined by the factor graph, as groups of variables
connected by the same factor. The idea of the Bethe entropy approximation is
to sum the entropies of all such subsets and then discount the “overcounted”
entropy corresponding to variables that appear in multiple subsets. The
pseudo-marginals
${\bm{\tilde{q}}}\equiv\\{q(x^{(i)}_{t+1},{\bm{x}}^{\text{pa}(i)}_{t},a_{t})\\}_{t=1,i=1}^{t=T-1,i=N_{e}}\cup\\{q({\bm{x}}^{\text{pa}(r)}_{t})\\}_{t=1,r=N_{e}+1}^{t=T,r=N_{e}+N_{r}}$
are all the local distributions that correspond to each factor. The pseudo-
marginals are locally consistent at the variables, i.e., two pseudo-marginals
that contain the same variable should provide the same marginal for that
variable. However, there is no further need for consistency, and in
particular, they do not need to correspond to the marginals of global joint
distribution. The (convex) domain of pseudo-marginals contains the (also
convex) domain of the marginals, but it is larger, so switching from marginals
to pseudo-marginals in any optimization problem relaxes it and provides an
upper bound on the original optimization problem. See [38] for more details.

Since the planning entropy is a linear combination of the standard entropy of
the factor graph $H^{\text{marginal}}({\bm{q}})$ and two local entropies per
time step, we can simply approximate each Shannon entropy by its corresponding
Bethe entropy to get the Bethe planning entropy

$\displaystyle H^{\text{planning}}({\bm{q}})$
$\displaystyle=H^{\text{marginal}}({\bm{q}})+\sum_{t=1}^{T-1}H_{\bm{q}}({\bm{x}}_{t})-H_{\bm{q}}({\bm{x}}_{t},a_{t})$
$\displaystyle\approx
H^{\text{marginal}}_{\text{Bethe}}({\bm{\tilde{q}}})+\sum_{t=1}^{T-1}H_{\text{Bethe}}(q({\bm{x}}_{t}))-H_{\text{Bethe}}(q({\bm{x}}_{t},a_{t}))$
$\displaystyle=\sum_{i=1}^{N_{e}}H_{\bm{\tilde{q}}}(x^{(i)}_{1})+\sum_{t=1}^{T-1}\Big{(}H_{\text{Bethe}}(q({\bm{x}}_{t}))+\sum_{i=1}^{N_{e}}H_{\bm{\tilde{q}}}(x^{(i)}_{t+1}|{\bm{x}}^{\text{pa}(i)}_{t},a_{t})\Big{)}=H_{\text{Bethe}}^{\text{planning}}({\bm{\tilde{q}}}).$

Note that we sometimes use the superscript marginal for consistency with the
main text, but the marginal entropy is simply the standard entropy so that
superscript can be safely dropped.

##  Appendix D Value BP message updates

We are interested in optimizing the cost function

$\frac{1}{\lambda}\max_{\bm{\tilde{q}}}\Big{(}-E_{\lambda}({\bm{\tilde{q}}})+\epsilon
H_{\text{Bethe}}^{\text{marginal}}({\bm{\tilde{q}}})+(1-\epsilon)H_{\text{Bethe}}^{\text{planning}}({\bm{\tilde{q}}})\Big{)}$

As we saw in Appendix C,
$H_{\text{Bethe}}^{\text{marginal}}({\bm{\tilde{q}}})$ and
$H_{\text{Bethe}}^{\text{planning}}({\bm{\tilde{q}}})$ are linear combinations
of local entropies over small subsets of variables defined by the factor
graph, therefore so is their linear combination. This score function can be
seen as a standard Bethe free energy with non-standard entropy weightings. We
can directly derive LBP-like message updates by using these modified weights
instead ([15]), resulting in the following VBP message updates

$\displaystyle Q({\bm{x}}^{\text{pa}(i)}_{t},a_{t})$
$\displaystyle=\sum_{x^{(i)}_{t+1}}m_{\text{b}}(x^{(i)}_{t+1})p(x^{(i)}_{t+1}|{\bm{x}}^{\text{pa}(i)}_{t},a_{t})$
$\displaystyle m_{\text{b}}({\bm{x}}^{\text{pa}(i)}_{t})$
$\displaystyle=\Big{(}\sum_{a}(Q({\bm{x}}^{\text{pa}(i)}_{t},a_{t})n^{(i)}(a_{t}))^{\frac{1}{\epsilon}}\Big{)}^{\epsilon}$
$\displaystyle m^{(i)}(a_{t})$
$\displaystyle=\Big{(}\sum_{{\bm{x}}^{\text{pa}(i)}_{t}}\left(\frac{Q({\bm{x}}^{\text{pa}(i)}_{t},a_{t})}{m_{\text{b}}({\bm{x}}^{\text{pa}(i)}_{t})}\right)^{\frac{1}{\epsilon}}m_{\text{f}}({\bm{x}}^{\text{pa}(i)}_{t})m_{\text{b}}({\bm{x}}^{\text{pa}(i)}_{t})\Big{)}^{\epsilon}$
$\displaystyle n^{(i)}(a_{t})$ $\displaystyle=\prod_{k\neq i}m^{(k)}(a_{t})$
$\displaystyle m_{\text{f}}(x^{(i)}_{t+1})$
$\displaystyle=\sum_{{\bm{x}}^{\text{pa}(i)}_{t},a}\left(\frac{Q({\bm{x}}^{\text{pa}(i)}_{t},a_{t})n^{(i)}(a_{t})}{m_{\text{b}}({\bm{x}}^{\text{pa}(i)}_{t})}\right)^{\frac{1}{\epsilon}}m_{\text{f}}({\bm{x}}^{\text{pa}(i)}_{t})m_{\text{b}}({\bm{x}}^{\text{pa}(i)}_{t})\frac{p(x^{(i)}_{t+1}|{\bm{x}}^{\text{pa}(i)}_{t},a_{t})}{Q({\bm{x}}^{\text{pa}(i)}_{t},a_{t})}$
$\displaystyle m_{\text{f}}({\bm{x}}^{\text{pa}(i)}_{t})$
$\displaystyle=\prod_{k\in\text{pa}(i)}n^{(i)}_{\text{f}}(x^{(k)}_{t})$
$\displaystyle n^{(i)}_{\text{f}}(x^{(j)}_{t})$
$\displaystyle=m_{\text{f}}(x^{(j)}_{t})\prod_{k|j\in\text{pa}(k),k\neq
i}n^{(k)}_{\text{b}}(x^{(j)}_{t})~{}~{}\text{ towards parents of entity }i$
$\displaystyle n^{(i)}_{\text{b}}(x^{(j)}_{t})$
$\displaystyle=\sum_{\\{x^{(k)}_{t}\\}_{k\neq
j}}m_{\text{b}}({\bm{x}}^{\text{pa}(i)}_{t})\prod_{k\neq
i}n^{(i)}_{\text{f}}(x^{(k)}_{t})~{}~{}\text{ from parents of entity }i$
$\displaystyle m_{\text{b}}(x^{(j)}_{t})$
$\displaystyle=\prod_{k|j\in\text{pa}(k)}n^{(k)}_{\text{b}}(x^{(j)}_{t})$

The following messages should be held constant

$\displaystyle m_{\text{b}}({\bm{x}}^{\text{pa}(i)}_{t})$
$\displaystyle=\exp(\lambda R({\bm{x}}^{\text{pa}(i)}_{t}))~{}~{}\forall
i>N_{e}$ $\displaystyle m_{\text{f}}(x^{(i)}_{1})$
$\displaystyle=P(x^{(i)}_{1})$

See Fig. 4 to track the correspondence between the above messages and the
factor graph of the factored MDP. Note that most updates correspond exactly
with standard loopy BP, and only a few (the ones involving $\epsilon$) are
specific to VBP. These messages updates should be iterated until convergence
or for a fixed number of iterations. Two tricks to improve convergence are (a)
damping the message updates in log space; (b) as the iterations progress,
anneal between LBP ($\epsilon=1$) and VBP (very small $\epsilon$). We
typically use both, with a damping of 0.5 (i.e., the mean of the old and the
new message in log space) and anneal by using as epsilon
$\max\\{0.01,1.0/\text{iter}\\}$, where “iter” is the iteration number.
Scheduling also plays a role in convergence. We propagate the messages by
alternating backward and forward schedules in an outer loop, and solving each
time slice to convergence in an inner loop. We did observe a correlation
between the quality of the solutions and VBP converging.

Figure 4: Correspondence between the message passing updates and the
factorized MDP.

##  Appendix E An application of VI for planning: determinization in hindsight

In the planning literature, many algorithms make the assumption of a
deterministic environment ([13, 17, 16]; etc). In order to extend these
algorithms to the case of stochastic dynamics, the idea of _determinization in
hindsight_ has been developed [41]. It starts by factoring out all the
stochasticity from the transition probability into the random variables
${\bm{\gamma}}=\\{\gamma_{t}\\}_{t=1}^{t=T-1}$

$P(x_{t+1}|x_{t},a_{t})=\sum_{\gamma_{t}}P_{\text{det}}(x_{t+1}|x_{t},a_{t},\gamma_{t})P(\gamma_{t})=\mathbb{E}_{\gamma_{t}}[P_{\text{det}}(x_{t+1}|x_{t},a_{t},\gamma_{t})].$

This factorization is always possible, and ${\bm{\gamma}}$ is known as the
collection of exogenous variables in structural equation modeling [28, 3, 31].
We can use it to obtain an upper bound on the exact utility

$F_{\lambda=1}^{\text{planning}}=\log\max_{\bm{\pi}}\mathbb{E}_{\bm{\pi}}[\mathbb{E}_{\bm{\gamma}}[\exp(R_{\bm{\gamma}}({\bm{x}},{\bm{a}}))]]\leq\log\mathbb{E}_{\bm{\gamma}}[\max_{\bm{\pi}}\mathbb{E}_{\bm{\pi}}[\exp(R_{\bm{\gamma}}({\bm{x}},{\bm{a}}))]]=F_{\lambda=1}^{\text{det.
planning}}$

that has a natural interpretation: the dynamics of a stochastic environment
depends on random variables $\bm{\gamma}$ that get drawn at each time step, if
we knew their value ahead of time, we could treat the environment as
deterministic, find an optimal plan and estimate its utility. Because we have
_hindsight_ , i.e., the deterministic planner can see the value of future
$\gamma_{t}$ ahead of time, it can make better choices than if these were
revealed online, so the utility estimation, if exact, upper bounds the
original quantity. See [41] for more details.

###  E.1 Closed form determinization for standard MDPs

The above process uses a sample average: it alternates between sampling
${\bm{\gamma}}\sim P(\bm{\gamma})$ and running a deterministic planner. Using
the variational framework for planning it is possible to find a closed form
expression for $F_{\lambda=1}^{\text{det. planning}}$.

First, let us expand the determinization bound

$F_{\lambda=1}^{\text{det.
planning}}\overset{\text{Eq.~{}\eqref{eq:app_var}}}{=}\max_{q(\gamma)}\langle\log\max_{\bm{\pi}}\mathbb{E}_{\bm{\pi}}[\exp(R({\bm{x}},{\bm{a}}))]\rangle_{q(\gamma)}+H_{q}(\gamma)=\max_{q(\gamma)}\langle
F_{\lambda=1,\gamma}^{\text{planning}}\rangle_{q(\gamma)}+H_{q}(\gamma)$

in terms of our exact planning utility for a given future $\bm{\gamma}$,
$F_{\lambda=1,\gamma}^{\text{planning}}$. This quantity and the corresponding
bound are defined as the usual ones, but conditioned on $\bm{\gamma}$

$F_{\lambda=1,\gamma}^{\text{planning}}=\max_{\bm{q}}F_{\lambda=1}^{\text{planning}}({\bm{q}}|\bm{\gamma});~{}~{}~{}~{}~{}~{}F_{\lambda=1}^{\text{planning}}({\bm{q}}|\bm{\gamma})=-E_{\lambda=1}(\bm{q}|\bm{\gamma})+H^{\text{planning}}(\bm{q}|\bm{\gamma}).$

To be precise, the corresponding energy $E(\bm{q}|\bm{\gamma})$ and entropy
$H^{\text{det. planning}}(\bm{q}|\bm{\gamma})$ terms are

$\displaystyle E_{\lambda=1}(\bm{q}|\bm{\gamma})=$ $\displaystyle-\langle\log
P_{\text{det}}(x_{1}|\gamma_{0})\rangle_{q(x_{1}|\gamma_{0})}$
$\displaystyle-\sum_{t=1}^{T-1}\langle\log
P_{\text{det}}(x_{t+1}|x_{t},a_{t},\gamma_{t})+R_{t}(x_{t},a_{t},x_{t+1})\rangle_{q(x_{t+1},x_{t},a_{t}|\gamma_{t})}$
(11) $\displaystyle H^{\text{det. planning}}(\bm{q}|\bm{\gamma})=$
$\displaystyle
H_{\bm{q}}(x_{1}|\gamma_{1})+\sum_{t=1}^{T-1}H_{\bm{q}}(x_{t+1}|x_{t},a_{t},\gamma_{t})$
(12)

Since the conditional dynamics are now deterministic, $\log
P_{\text{det}}(x_{t+1}|x_{t},a_{t},\gamma_{t})$ can only take the values 0 and
$-\infty$. When maximizing over $\bm{q}$, it is obvious we must choose
$q(x_{t+1}|x_{t},a_{t},\gamma_{t})=P_{\text{det}}(x_{t+1}|x_{t},a_{t},\gamma_{t})~{}\forall_{t}$
and $q(x_{1}|\gamma_{0})=P_{\text{det}}(x_{1}|\gamma_{0})$, to avoid putting
any mass on the $-\infty$ value of the previous term, which would result in
the whole score function becoming $-\infty$. For such choice, $H^{\text{det.
planning}}(\bm{q}|\bm{\gamma})=0$ and the mentioned energy term also cancels,
leaving

$F_{\lambda=1}^{\text{planning}}(\bm{q}|\bm{\gamma})=\sum_{t=1}^{T-1}\langle
R_{t}(x_{t},a_{t},x_{t+1})\rangle_{P_{\text{det}}(x_{t+1}|x_{t},a_{t},\gamma_{t})q(x_{t},a_{t}|\gamma_{t})}.$
(13)

Therefore, _determinization for multiplicative exponentiated rewards can be
exactly expressed as a concave variational bound of the form_

$\displaystyle F_{\lambda=1}^{\text{det. planning}}$
$\displaystyle=\max_{{\bm{q}}_{\gamma}}\langle
F_{\lambda=1}^{\text{planning}}(\bm{q}|\bm{\gamma})\rangle_{q({\bm{\gamma}})}+\sum_{t=0}^{T-1}\langle\log
P(\gamma_{t})\rangle_{q(\gamma_{t})}+H_{\bm{q}}(\gamma_{t})$ (14)
$\displaystyle\text{s.t.
}\sum_{a_{1},\gamma_{1}}q(x_{1},a_{1},\gamma_{1})=\sum_{\gamma_{0}}P_{\text{det}}(x_{1}|\gamma_{0})q(\gamma_{0})$
$\displaystyle\sum_{a_{t+1},\gamma_{t+1}}q(x_{t+1},a_{t+1},\gamma_{t+1})=\sum_{x_{t},a_{t},\gamma_{t}}P_{\text{det}}(x_{t+1}|x_{t},a_{t},\gamma_{t})q(x_{t},a_{t},\gamma_{t})~{}~{}\forall_{t};$
$\displaystyle
q(x_{t},\gamma_{t})=\sum_{a_{t}}q(x_{t},a_{t},\gamma_{t})~{}~{}\forall_{t};~{}~{}~{}~{}~{}q(x_{t},a_{t},\gamma_{t})>0~{}\forall
t,$

where we have defined
${{\bm{q}}_{\gamma}}\equiv{\\{q(x_{t},a_{t},\gamma_{t})\\}_{t=1}^{T-1}\cup
q(\gamma_{0})}$

###  E.2 Factored MDP

In the case of a standard MDP, it is trivial to see that determinization,
which provides an upper bound on the exact utility, cannot improve on our VI
bound $F_{\lambda=1}^{\text{planning}}$, since the latter is exact. In this
section we will show that this is also the case for factored MDPs in which the
deterministic MAP problem is solved using an LP MAP relaxation.

First, we will show that $F_{\lambda=1}^{\text{planning}}$ can be rewritten in
a more convenient way to compare with determinization, in the standard, non-
factored case. We expand

$F_{\lambda=1}^{\text{planning}}=\max_{\bm{q}}F_{\lambda=1}^{\text{planning}}({\bm{q}})=\max_{\bm{q}}-E_{\lambda=1}({\bm{q}})+H^{\text{planning}}({\bm{q}})$

and substitute

$\log P(x_{t+1}|x_{t},a_{t})=\langle\log
P_{\text{det}}(x_{t+1}|x_{t},a_{t},\gamma_{t})\rangle_{q(\gamma_{t}|x_{t+1},x_{t},a_{t})}-\sum_{t=1}^{T}\operatorname{KL}(q(\gamma_{t}|x_{t+1},x_{t},a_{t})||P(\gamma_{t}))$

inside $E_{\lambda=1}({\bm{q}})$ to get

$F_{\lambda=1}^{\text{planning}}=\max_{\bm{q}}\langle
F_{\lambda=1}^{\text{planning}}(\bm{q}|\bm{\gamma})\rangle_{q({\bm{\gamma}})}+\sum_{t=0}^{T-1}\langle\log
P(\gamma_{t})\rangle_{q(\gamma_{t})}+H_{\bm{q}}(x_{t+1},\gamma_{t}|x_{t},a_{t})$

where
$q(x_{t+1}|x_{t},a_{t},\gamma_{t})=P_{\text{det}}(x_{t+1}|x_{t},a_{t},\gamma_{t})~{}\forall_{t}$
for the same reasons as above. To compact notation, we use the convention that
$H_{\bm{q}}(x_{1},\gamma_{0}|x_{0},a_{0})=H_{\bm{q}}(x_{1},\gamma_{0})$ (since
$x_{0},a_{0}$ are not defined), and similarly
$H_{\bm{q}}(\gamma_{0}|x_{0},a_{0})=H_{\bm{q}}(\gamma_{0})$. Noting that
$H_{\bm{q}}(x_{t+1},\gamma_{t}|x_{t},a_{t})=H_{\bm{q}}(x_{t+1}|x_{t},a_{t},\gamma_{t})+H_{\bm{q}}(\gamma_{t}|x_{t},a_{t})$
and that $H_{\bm{q}}(x_{t+1}|x_{t},a_{t},\gamma_{t})=0$ because
$q(x_{t+1}|x_{t},a_{t},\gamma_{t})$ is deterministic, we get

$\displaystyle F_{\lambda=1}^{\text{planning}}$
$\displaystyle=\max_{{\bm{q}}_{\gamma}}\langle
F_{\lambda=1}^{\text{planning}}(\bm{q}|\bm{\gamma})\rangle_{q({\bm{\gamma}})}+\sum_{t=0}^{T-1}\langle\log
P(\gamma_{t})\rangle_{q(\gamma_{t})}+H_{\bm{q}}(\gamma_{t}|x_{t},a_{t})$ (15)
$\displaystyle\text{s.t.
}\sum_{a_{1},\gamma_{1}}q(x_{1},a_{1},\gamma_{1})=\sum_{\gamma_{0}}P_{\text{det}}(x_{1}|\gamma_{0})q(\gamma_{0})$
$\displaystyle\sum_{a_{t+1},\gamma_{t+1}}q(x_{t+1},a_{t+1},\gamma_{t+1})=\sum_{x_{t},a_{t},\gamma_{t}}P_{\text{det}}(x_{t+1}|x_{t},a_{t},\gamma_{t})q(x_{t},a_{t},\gamma_{t})~{}~{}\forall_{t};$
$\displaystyle
q(x_{t},\gamma_{t})=\sum_{a_{t}}q(x_{t},a_{t},\gamma_{t})~{}~{}\forall_{t};~{}~{}~{}~{}~{}q(x_{t},a_{t},\gamma_{t})>0~{}\forall
t.$

Compare Eqs. (14) and (15). Although we already knew this, with this
particular formulation it is very easy to see that
$F_{\lambda=1}^{\text{planning}}\leq F_{\lambda=1}^{\text{det. planning}}$,
since they are identical except for the terms
$H_{\bm{q}}(\gamma_{t}|x_{t},a_{t})\leq H_{\bm{q}}(\gamma_{t})~{}\forall t$.
This holds for the same ${\bm{q}}$, but also establishes a relation between
both upper bounds at their maximum wrt ${\bm{q}}$.

With this result, we can return to the factored MDP case. As in the main text,
here we will consider $R$ to be only a function of $x_{t}$. We can take Eqs.
(14) and (15) and relax $\bm{q}_{\gamma}$ from the marginal polytope $\cal{M}$
to the local polytope $\cal{L}$, using the pseudo-marginals
$\bm{\tilde{q}}_{\gamma}\equiv\\{q({\bm{x}}^{\text{pa}(i)}_{t},a_{t},\gamma^{(i)}_{t})\\}_{t=1,i=1}^{t=T-1,i=N_{e}}\cup\\{q(\gamma^{(i)}_{0})\\}_{i=1}^{N_{e}}\cup\\{q({\bm{x}}^{\text{pa}(r)}_{t})\\}_{t=1,r=N_{e}+1}^{t=T,r=N_{e}+N_{r}}$.
For Eq. (14), and substituting the value from Eq. (13), this produces a
factored MDP upper bound corresponding to determinization with $\lambda=1$

$\displaystyle\hat{F}_{\lambda=1}^{\text{det.
planning}}=\max_{\bm{\tilde{q}}_{\gamma}}\sum_{t=1}^{T}$
$\displaystyle\Big{(}\sum_{r=N_{e}+1}^{N_{e}+N_{r}}\langle
R({\bm{x}}^{\text{pa}(r)}_{t})\rangle_{q({\bm{x}}^{\text{pa}(r)}_{t})}+\sum_{i=1}^{N_{e}}\langle\log
P(\gamma^{(i)}_{t-1})\rangle_{q(\gamma^{(i)}_{t-1})}+H_{\bm{q}}(\gamma^{(i)}_{t-1})\Big{)}$
s.t.
$\displaystyle\sum_{a_{1},\gamma^{(i)}_{1}}q(x^{(i)}_{1},a_{1},\gamma^{(i)}_{1})$
$\displaystyle=\sum_{\gamma^{(i)}_{0}}P_{\text{det}}(x^{(i)}_{1}|\gamma^{(i)}_{0})q(\gamma^{(i)}_{0})~{}\forall
i$
$\displaystyle\sum_{a_{t+1},\gamma^{(i)}_{t+1}}q(x^{(i)}_{t+1},a_{t+1},\gamma^{(i)}_{t+1})$
$\displaystyle=\sum_{{\bm{x}}^{\text{pa}(i)}_{t},a_{t},\gamma^{(i)}_{t}}p_{\text{det}}(x^{(i)}_{t+1}|{\bm{x}}^{\text{pa}(i)}_{t},a_{t},\gamma^{(i)}_{t})q({\bm{x}}^{\text{pa}(i)}_{t},a_{t},\gamma^{(i)}_{t})~{}~{}\forall_{t},i\in
1,\ldots,N_{e}$ $\displaystyle\text{and pseudo-marginal constraints},$

which can be interpreted as an outer maximization over $q(\bm{\gamma})$ and an
inner MAP optimization of the (conditional on $\bm{\gamma}$) deterministic
dynamics problem, solved with an LP relaxation [34]. For Eq. (15), and also
substituting the value from Eq. (13), this results in a different formulation
of the problem Eq. (6)

$\displaystyle\hat{F}_{\lambda=1}^{\text{planning}}=\max_{\bm{\tilde{q}}_{\gamma}}\sum_{t=1}^{T}\Big{(}\sum_{r=N_{e}+1}^{N_{e}+N_{r}}$
$\displaystyle\langle
R({\bm{x}}^{\text{pa}(r)}_{t})\rangle_{q({\bm{x}}^{\text{pa}(r)}_{t})}+\sum_{i=1}^{N_{e}}\langle\log
P(\gamma^{(i)}_{t-1})\rangle_{q(\gamma^{(i)}_{t-1})}+H_{\bm{q}}(\gamma^{(i)}_{t-1}|x_{t-1},a_{t-1})\Big{)}$
s.t.
$\displaystyle\sum_{a_{1},\gamma^{(i)}_{1}}q(x^{(i)}_{1},a_{1},\gamma^{(i)}_{1})$
$\displaystyle=\sum_{\gamma^{(i)}_{0}}P_{\text{det}}(x^{(i)}_{1}|\gamma^{(i)}_{0})q(\gamma^{(i)}_{0})~{}\forall
i$
$\displaystyle\sum_{a_{t+1},\gamma^{(i)}_{t+1}}q(x^{(i)}_{t+1},a_{t+1},\gamma^{(i)}_{t+1})$
$\displaystyle=\sum_{{\bm{x}}^{\text{pa}(i)}_{t},a_{t},\gamma^{(i)}_{t}}p_{\text{det}}(x^{(i)}_{t+1}|{\bm{x}}^{\text{pa}(i)}_{t},a_{t},\gamma^{(i)}_{t})q({\bm{x}}^{\text{pa}(i)}_{t},a_{t},\gamma^{(i)}_{t})~{}~{}\forall_{t},i\in
1,\ldots,N_{e}$ $\displaystyle\text{and pseudo-marginal constraints},$

In the same way as for the non-factorized MDP, it is easy to see that
$\hat{F}_{\lambda=1}^{\text{planning}}\leq\hat{F}^{\text{det. planning}}$,
while both upper bound the exact $F_{\lambda=1}^{\text{planning}}$. This means
that _in the case of (exponentiated) multiplicative rewards ( $\lambda=1$) our
proposed bound should be no worse than determinization_. Also,
$\hat{F}_{\lambda=1}^{\text{det. planning}}$ allows us to compute the
determinization bound as a single, solvable, concave optimization, rather than
having to resort to sampling and solving a theoretically infinite number of
deterministic problems (this holds when each problem would be solved using an
LP relaxation of the MAP problem).

##  Appendix F Empirical validation: details

###  F.1 Synthetic MDPs

The synthetic MDPs use binary state and actionss, $T=4$ timesteps and follow
the connectivity scheme of Fig. 1[Right]. The stochasticity level of the MDPs
is controlled by generating their dynamics according to
$P(x_{t+1}|x_{t},a_{t})=\bar{P}_{x_{t+1},x_{t},a_{t}}^{s}/Z_{x_{t},a_{t}}$
where $\bar{P}_{x_{t+1},x_{t},a_{t}}\sim U[0,1]$ and choosing the exponent $s$
and the divisor $Z_{x_{t},a_{t}}$ to obtain normalized distributions of the
desired total normalized entropy $H_{\text{MDP}}$. Each entity has two parent
entities and we provide a single reward of 1 for reaching the state 0 of
entity 1 at the last step $T$. Having a single reward allows comparing methods
with arbitrary $\lambda$ and create a pure planning problem (see Section 4.2).

###  F.2 Reactivity avoidance of MAP and MMAP inference

As discussed in Section 4.2, a common flaw among all inference types (other
than planning inference) is that they cannot anticipate reacting to the
environment. Even the tightest one, MMAP, makes plans assuming that it will
only be able to take a predefined sequence of actions. This can severely
underestimate the value of an action if, to extract future reward, reacting to
the environment state is necessary. In other words, MMAP is optimal only when
the best policy $\pi_{t}(a_{t}|x_{t})$ ignores the state and can be
represented as $\pi_{t}(a_{t})$. MAP inference shares that limitation, and
additionally lacks path integration. I.e., MAP inference makes decisions based
on the score of a single action-state trajectory, rather than a combination of
these.

In order to analyze the behavior of different inference types of inference
with reactivity, we create an MDP with two entities, each with categorical
states $0,\ldots,5$. We use a horizon of $T=7$ timesteps, placing reward only
on the last timestep. There are a total of 8 actions.

The first entity describes the location of the agent, and has a special goal
state “0” that needs to be reached at the last timestep. The second entity
describes the dynamics of the first entity and the reward achieved at the
goal. Actions $0,\ldots,5$ allow the agent to move to any location in a single
step. Actions $6,7$ allow the agent to modify the dynamics of the first
entity, respectively decreasing or increasing the needed level of reactivity
of a planner.

In more detail, the second entity has deterministic dynamics and acts as a
knob that the agent can freely turn up or down at each timestep:

$x^{(2)}_{t+1}=\begin{cases}x^{(2)}_{t}-1,&\text{if }a=6\text{ and
}0<x^{(2)}_{t}\\\ x^{(2)}_{t}+1,&\text{if }a=7\text{ and }x^{(2)}_{t}<5\\\
x^{(2)}_{t},&\text{otherwise}.\end{cases}$

The values of that knob $x^{(2)}_{t}$ alter the dynamics of the first entity,
requiring more or less reactivity:

$x^{(1)}_{t+1}=\begin{cases}\text{if }x^{(1)}_{t}=0\text{ or }6\leq a&\sim
U[1,5]\\\ \text{if }x^{(1)}_{t}\neq 0\text{ and }0\leq
a<6&\begin{cases}x^{(1)}_{t}+a\mod 6&\text{ with probability }x^{(2)}_{t}/5\\\
0&\text{ with probability }1-x^{(2)}_{t}/5\\\ \end{cases}\end{cases}$

In words, if the agent is at the goal state 0 or the reactivity of the
environment is modified, it will jump to a random non-goal state. If the agent
is not at the goal state and $x^{(2)}_{t}=5$, it can use the actions
$0,\ldots,5$ to jump to any desired state with absolute certainty. As the
value of the knob $x^{(2)}_{t}$ is reduced, the agent starts losing control of
where it is going, and it is instead more likely to jump to the goal. Now, the
final ingredient, the reward, is:

$R_{t}(x^{(1)}_{t},x^{(2)}_{t})=\begin{cases}1.0,&\text{if
}x^{(1)}_{t}=0\text{ and }x^{(2)}_{t}=5\text{ and }t=T\\\ 0.33,&\text{if
}x^{(1)}_{t}=0\text{ and }x^{(2)}_{t}<5\text{ and }t=T\\\
0.0,&\text{otherwise}.\end{cases}$

When we put all these pieces together, we have the following situation: to
achieve the maximum reward, we want to keep the the knob $x^{(2)}_{t}$ at its
“maximum reactivity” value of 5. But when we do that, _the action that will
take the agent to the goal depends on the state it is located at_. Therefore
the agent needs to _react_ to the location that it is at. And when planning
ahead, to properly choose the best first action to take, the agent needs to
score the options according to this future ability to react to actions. By
turning the knob $x^{(2)}_{t}$ all the way down to 0, we can jump to the goal
by taking any action $0,\ldots,6$, regardless of where we are, i.e., we do not
need any reactivity.

We present this factored MDP, initialized at $x^{(1)}_{0}=0,x^{(2)}_{0}=5$ to
both VBP and SOGBOFA-LC*. Their behavior is very different. VBP is aware that
it can jump to the goal state $x^{(1)}_{T}=0$ at any time as long as it is not
there yet. So it just idles in the other states, keeping the reactivity at the
maximum level, until the last timestep, in which it jumps to the goal,
capturing a reward of 1.0. SOGBOFA-LC*, even though it is replanning at each
state, can only evaluate the possible actions by considering a single best
sequence of non-reactive actions at a time. Since the reactive plans that VBP
prefers are not “visible” to SOGBOFA-LC*, it decides to use the first 5
actions to reduce the reactivity of the environment all the way down to
$x^{(2)}_{T-1}=0$. In that way, it is guaranteed to be able to jump to the
goal in the last instant without having to be reactive, even if only capturing
a reward 0.33.

This example highlights how SOGBOFA-LC* (and MMAP planning in general)
struggle with reactive environments in which the best future actions cannot be
known at the current timestep. Essentially, the agent is myopic to the fact
that it will be able to replan and chooses to be conservative and _reactivity
avoidant_ , proposing sequences of actions that will be guaranteed to work in
any scenario, even if this reduces the obtained reward. We can make the
difference between VBP and SOGBOFA-LC*’s performance arbitrarily large simply
by increasing the number of states and planning horizon.

MAP inference presents the same non-reactivity problems as MMAP, combined with
the lack of integration across multiple paths. On this problem, the average
reward achieved by a perfect MAP inference agent is $\sim 0.13$.

###  F.3 Experimental Details on ICAPS 2011 International Probabilistic
Planning Competition Problems

Figure 5: Cumulative rewards of different types of inference on 6 problem
domains from the ICAPS 2011 International Probabilistic Planning Competition.
A small horizontal jitter was introduced in all data points for visual
clarity. This is the same plot as Figure 3 but with an additional CSVI
approach. Each cumulative reward is averaged over 30 simulations per instance
except for the CSVI approach that is averaged over 12 simulations.

####  F.3.1 Experimental Settings

Five different inference approaches, MFVI-Bwd, ARollout, SOGBOFA-LC, VI LP,
VBP, and a random agent are evaluated on 6 different domains (Crossing
traffic, Elevators, Game of life, Skill teaching, Sysadmin, and Traffic) used
in the ICAPS 2011 International Probabilistic Planning Competition ([32]).
Figure 5 shows the average cumulative reward of all domains plus the result of
CSVI in [39]. Table 2 shows the corresponding objective and reference for each
approach. There are 10 instances for each domain. All instances have a horizon
of 40 and a discount factor of 1. The cumulative reward is averaged over 30
simulations for each instance and the plotted error bar shows its standard
error of the mean. Given the current state $x_{t}$, the transition probability
$P(x_{t+1}|x_{t},a_{t})$, and the reward function $R$, each approach infers
the best next action $a_{t}$ at each time step. The reward function for all 6
domains only depend on the state $x$. The cumulative reward is the sum of the
rewards of the initial state and the following 39 states.

For all inference approaches, we run with a look ahead horizon of both 4 and
9. We follow the settings in [39] where the look ahead horizon is truncated if
it extends beyond the remaining time steps. The horizon with the higher
average cumulative reward is then selected for each instance. The experiment
evaluates over two look ahead horizons due to the observation that for some
domain instances a longer horizon may lead to worse estimates. For each
simulated run, the same set of random numbers are applied to the simulated
environment to reduce noise when comparing different inference approaches. All
experiments were ran on CPU machines in the cloud.

Table 2: Five different inference approaches used for comparison in the ICAPS 2011 International Probabilistic Planning Competition problems and their corresponding objective and reference. Inference Approach | Objective | Reference
---|---|---
MFVI-Bwd | — | [39]
CSVI-Bwd | — | [39]
ARollout | $\tilde{F}^{\mathrm{marginal}^{\text{U}}}_{\lambda=0}$ |  [10], [9]
SOGBOFA-LC | $\tilde{F}^{\mathrm{MMAP}}_{\lambda=0}$ | [8]
VI LP | $\hat{F}^{\mathrm{planning}}_{\lambda=0}$ | Ours
VBP | $\tilde{F}^{\mathrm{planning}}_{\lambda\approx 0}$ | Ours

####  F.3.2 MFVI-Bwd

The backward mean field variational inference apporach (MFVI-Bwd) is based on
the implementation in [39] (https://github.com/Zhennan-Wu/AISPFS), which the
mean field approximation $q_{\phi}$ is the product of independent factors.
$q_{\phi}$ is first obtained by maximizing the ELBO of the true posterior
under the condition that the maximum accumulated reward is reached using the
EM algorithm. The action is then selected based on the marginal
$q_{\phi}(\bm{a})$. In the experiment, the maximum number of iterations is set
to 100 and the convergence threshold is set to 0.1 for the EM algorithm. MFVI-
Bwd experiments were ran on a single CPU machine with 32 virtual cores. All 30
simulations were ran in parallel for each task instance.

####  F.3.3 CSVI-Bwd

The backward collapsed variational inference approach (CSVI-Bwd) is proposed
in [39]. It uses collapsed variational inference to effectively marginalizes
out variables other than the actions to achieve a tighter ELBO. The authors
have shown that this approach out performs MFVI-Bwd. The results of CSVI-Bwd
in Figure 5 is taken directly from [39]. They are averaged over 12 simulations
instead of 30 and are only based on a look ahead horizon of 9. The random
numbers are also different between the 12 simulations and the 30 simulations.
The CSVI approach took a lot longer to compute and we were not able to
complete the experiment with 30 simulations in time.

####  F.3.4 ARollout

The Algebraic Rollout Algorithm (ARollout) introduced in [10] is equivalent to
belief propagation when conditioned on actions as shown in [11]. In our
experimental setting which the search depth is fixed with no computation time
limit, ARollout is equivalent to the original SOGBOFA (symbolic online
gradient based optimization for factored actions) approach introduced in [9].
While ARollout performs approximate aggregate rollout simulations to evaluate
each action, SOGBOFA uses the gradient of the accumulated reward to update
actions for exploration and can be advantageous when the action space is large
given a computation time constraint. In this experiment, we use the results of
forward belief propagation to represent ARollout.

For each time step, the action with the highest estimated accumulated reward
is selected. The estimated accumulated reward is calculated by running a
forward pass on the factored MDP representing the problem conditioned on the
next action. ARollout experiments were ran on CPU machines with 2 virtual
cores. All simulations were ran in parallel.

####  F.3.5 SOGBOFA-LC

SOGBOFA-LC is based on the Lifted Conformant SOGBOFA implemented in [8]
(https://github.com/hcui01/SOGBOFA). SOGBOFA-LC has two differences over the
original SOGBOFA. First, it uses a lifted graph that saves computation by
identifying same operations. This improvement in computation speed doesn’t
have an effect in our experiment result since we provide enough computation
time for the given maximum search depth. Second, it uses conformant solutions,
which the evaluation of the next action is based on a linear rollout plan that
best supports the action. This is achieved by calculating the gradient with
respect to all actions within the search depth.

In our experiment, the search depth is set to 9 or 4 based on the look ahead
horizon. The number of gradient updates is set to 500 following the
experimental setting in [39]. The allowed time is set to 50000 per iteration,
which is sufficient for 500 updates given the assigned depth. SOGBOFA-LC
experiments were ran on CPU machines with 4 virtual cores sequentially.

Note that the RDDL files in the repository (https://github.com/hcui01/SOGBOFA
or https://github.com/Zhennan-Wu/AISPFS/tree/master/SOGBOFA) have different
initial states from the original competition for some instances in the
elevator and skill teaching domains. These differences were corrected to match
the original competition settings in our experiments. Modifications were also
made to use the same random numbers in the environment as other experiments
and to measure the standard error of the mean instead of the standard
deviation.

####  F.3.6 VI LP

The Variational Inference Linear Programming (VI LP) approach uses the GLOP
solver in Google’s OR-Tools ([29]) to solve the linear programming (LP)
problem derived from each task instance with the target of maximizing the
expected accumulated reward. Constraints on states that are specified in the
original RDDL problem is added to the LP problem. Among the six domains only
the elevator and crossing traffic domains have such constraints. These
constraints specifies that the elevator/robot can’t be at different
floors/locations at the same time step.

The solver is ran for each next action at each time step. The next action that
has the highest estimated expected accumulated reward based on the solver is
selected. See Section 3.1 for more details of this approach. VI LP experiments
were ran on CPU machines with 16 virtual cores. All simulations were ran in
parallel. For each iteration, LP solvers for each next action were ran
concurrently with multiprocessing.

####  F.3.7 VBP

Note that setting a value for $\lambda$ is meaningless if the reward can have
arbitrary scaling (and for the case of IPPC in which rewards are additive, the
scaling is indeed arbitrary). Therefore, we first normalize all rewards to
have a maximum point-to-point variation of 1.0. $\lambda$ is chosen as a
reward multiplier on top of that normalization. The closer to 0 that we set
$\lambda$ then, the closer we are to the additive limit. Although outside the
scope of this work, it is actually possible to set $\lambda$ to exactly zero
in VBP (the problem remains non-convex, but the message updates within each
time step become an LP problem). However, at least in a straighforward
implementation, we observed such message updates to not result in a good
optimization of our score function, and often not converging. We attribute
this to the degeneracy of the solutions in the limit. Using instead a small
value of $\lambda=0.3$ had a favorable effect on convergence, while remaining
close enough to the additive limit. Empirically, we found that most problems
were not very sensitive to the choice of $\lambda$.

For each time step, VBP messages are propagated concurrently for a maximum of
150K iterations with 0.1 damping. The $\epsilon$ value is annealed every 300
iterations from a value of 1 to 0.01 based on the formula described in
Appendix D. The next action with the highest expected accumulated reward is
then selected. See Section 3.1 and Appendix D for more details of this
approach. VBP experiments were ran on CPU machines with 2 virtual cores. All
simulations were ran in parallel.