# Formally Verified Approximate Policy Iteration

Maximilian Schäffeler
Technische Universität München
Germany &Mohammad Abdulaziz
King’s College London
United Kingdom

###### Abstract

We formally verify an algorithm for approximate policy iteration on Factored
Markov Decision Processes using the interactive theorem prover Isabelle/HOL.
Next, we show how the formalized algorithm can be refined to an executable,
verified implementation. The implementation is evaluated on benchmark problems
to show its practicability. As part of the refinement, we develop verified
software to certify Linear Programming solutions. The algorithm builds on a
diverse library of formalized mathematics and pushes existing methodologies
for interactive theorem provers to the limits. We discuss the process of the
verification project and the modifications to the algorithm needed for formal
verification.

##  1 Introduction

Markov Decision Processes (MDPs) are models of probabilistic systems, with
applications in AI, model checking, and operations research. In AI, for
instance, given a description of the world in terms of states and actions,
where actions can change those states in a randomised fashion, a solution is a
_policy_ that instructs a decision-maker as to which actions to choose in
every state, with the aim of accruing maximum possible _reward_. There is a
large number of methods to solve MDPs, most notably, value and policy
iteration, which are some of the oldest and most studied algorithms, and which
can compute policies with guarantees on the optimality of those policies.

In many safety-critical applications, e.g. in AI or autonomous systems, an
important factor of whether a system for solving MDPs can be used is the
trustworthiness of that system. One important aspect here is the assurance
that the output of the MDP solving system is correct. Such assurance can be
attained to some degree by testing and other software engineering methods.
However, the best guarantee one can hope is to mathematically prove the
correctness of the MDP solver. A successful way of mathematically proving
correctness properties of (i.e. formally verify) pieces of software is using
Interactive Theorem Provers (ITPs), which are formal mathematical systems that
one can use to devise machine-checked proofs. Indeed, ITPs have been used to
prove correctness properties of compilers [26], operating systems kernels
[24], model checkers [15], planning systems [3, 1, 2], and, most related to
the topic of this work, algorithms to solve MDPs [32]. A challenge with using
ITPs to prove algorithms correct, nonetheless, is that they require intense
human intervention. Thus for an ITP to be successfully employed in a serious
verification effort, novel ideas in the design of the software to be verified
as well as the underlying mathematical argument have to be made.

In this paper, we consider formally verifying algorithms for solving
_factored_ MDPs. A challenge to using MDPs to model realistic systems is
their, in many cases, enormous size. For such systems, MDPs are succinctly
represented as factored MDPs. The system’s state is characterised as an
assignment to a set of state variables and actions are represented in a
compact way by exploiting the structure present in the system, e.g. the fact
that each action only changes a small set of state variables. Such
representations are common in AI [16, 31, 38] and in model checking [18, 12].
Although ITPs have been used prove the correctness of multiple types of
software and algorithms, including algorithms on MDPs, algorithms on factored
MDPs are particularly challenging. The root of this difficulty is that the
succinctness of the representation comes at a cost. Naively finding a solution
for a given factored MDP could entail the construction of structures that are
exponentially bigger than the factored MDP. This necessitates the usage of
advanced data structures and heuristics, in addition to combining a large
number of computational techniques, to avoid that full exponential blow up.

Our main contribution is that we formally verify the algorithm by Guesterin et
al. [16] using the Isabelle/HOL [28] theorem prover. This algorithm computes
approximate policies, i.e. sub-optimal policies with guarantees on their
optimality, for one type of factored MDPs. For instance, the algorithm we
consider here combines scoped functions and decision lists, which are data
structures that exploit factoring of the representation, linear programming,
in addition to probabilistic reasoning and dynamic programming. The
combination of this wide range of mathematical/algorithmic concepts and
techniques is what makes this algorithm particularly hard from a verification
perspective. To get an idea of the scale, the reader is advised to look at
Fig. 2, which shows the hierarchy of concepts and definitions which we had to
develop (aka _formalize_) within Isabelle/HOL just to be able to state the
algorithm and its correctness statement. This is, of course, in addition to
the predefined notions of analysis, probabilities, and MDPs, which already
exist in Isabelle/HOL. Furthermore, to be able to prove the algorithm correct,
we had to design an architecture of the implementation that allows feasible
verification. This architecture mixes verification and certification: we
verified the entire algorithm, except for linear program solving, for which we
build a verified a certificate checker. In addition to proving the algorithm
correct, we obtain a formally verified implementation of the algorithm, the we
experimentally show to be practical. Our work, as far as we are aware, is the
first work on verifying algorithms for solving factored MDPs. We believe that
our work provides infrastructure as well as methodological insights enabling
the verification of other algorithms for factored MDPs, e.g. algorithms for
planning under uncertainty and probabilistic model checking.

##  2 Background

We introduce the interactive theorem prover Isabelle/HOL, our formal
development of factored MDPs and show how it relates to existing
formalizations of MDPs in Isabelle/HOL.

###  2.1 Isabelle/HOL

An interactive theorem prover (ITP) is a program that implements a formal
mathematical system, in which definitions and theorem statements can be
expressed, and proofs are constructed from a set of axioms. To prove a fact in
an ITP, the user only provides the high-level steps while the ITP fills in the
details at the level of axioms. Specifically, our developments use the
interactive theorem prover Isabelle/HOL [28] based on Higher-Order Logic, a
combination of functional programming with logic. Isabelle is highly
trustworthy, as the basic inference rules are implemented in a small, isolated
kernel. Outside the kernel, several tools implement proof tactics, data types,
recursive functions, etc..

For improved readability, our presentation of definitions and theorems
deviates slightly from the formalization in Isabelle/HOL. Specifically, we use
subscript notation for list indexing and some function applications. We use
parentheses for function application opposed to juxtaposition in Isabelle/HOL.
For a list $\mathit{xs}$, $\mathit{len}(\mathit{xs})$ returns the length,
while $\mathit{xs}_{:i}$ returns the first $i$ elements of $\mathit{xs}$. List
concatenation is written as $\mathit{xs}\mathit{\cdot}\mathit{ys}$,
$x~{}\mathit{::}~{}\mathit{xs}$ inserts $x$ at the front of the list
$\mathit{xs}$, $\mathit{map}(f,\mathit{xs})$ applies $f$ to every element of
$\mathit{xs}$. Finally, $\mathit{\mathbb{N}}_{<i}$ is short for the first $i$
natural numbers (including $0$).

###  2.2 Factored MDPs

Factored MDPs are compactly represented MDPs that exploit commonly found
structure in large MDPs. This can often lead to an exponential reduction in
the size of the model. Common formats to store factored MDPs are JANI [10],
the PRISM language [18], and RDDL [31]. We implement the factoring described
in [16]. In Isabelle/HOL, we define factored MDPs using _locales_ [8] (see
Listings A.2.1, A.2.1). A locale introduces a mathematical context with
constants and assumptions, in which we develop our formalization. Locales can
be instantiated with concrete constants and a proof that discharges the
assumptions of the locale, yielding all the theorems proved within the locale.
For example, we instantiate the existing MDP formalization with the factored
systems introduced in this section, and are therefore able to reuse important
definitions and theorems. Moreover, reducing factored MDPs to a tested and
reviewed library improves the trust in our definitions.

##### State Space

A state of the MDP is an assignment of values to the $\mathit{\mathit{n}}$
state variables. Each state variable $i\in\mathit{\mathbb{N}}_{<n}$ has a
finite, nonempty domain $\mathit{X}_{i}$. In Isabelle/HOL, we implement a
state $x$ of the MDP with a map, i.e. a function with an explicit domain
$\mathit{Dom}(x)$ (see Listing A.2.1). A _partial state_ is a map where only a
subset of the state variables are assigned
($\mathit{Dom}(x)\subseteq\mathit{\mathbb{N}}_{<\mathit{\mathit{n}}}$), but
all entries are valid, i.e. $x_{i}\in\mathit{X}_{i}$. The set $\mathit{X}$ of
_states_ of the MDP consists of all partial states $x$ where
$\mathit{Dom}(x)=\mathit{\mathbb{N}}_{<\mathit{\mathit{n}}}$. The domain of a
state can be restricted to a set $Y$ with the notation $x|_{Y}$. One partial
state $x$ is called _consistent_ with another partial state $t$ (written
$x\sqsubseteq t$), if and only if $x|_{\mathit{Dom}(t)}=t$.

##### Example

As a running example, we will use a model of a computer network with ring
topology from Guestrin et al.’s original paper (see Figure 1). In the ring,
each machine could be either working or broken, and their states change
stochastically. Each machine $C_{i}$’s state of operation is characterised by
a state variable $i$, s.t. all domains $\mathit{X}_{i}=\\{\text{\text{W},
\text{B}}\\}$. In a ring with four machines, one valid partial state is
$s\coloneqq[1\mapsto\text{\text{W}},3\mapsto\text{\text{B}}]$,
$\mathit{Dom}(s)=\\{1,3\\}$. It holds that
$s\sqsubseteq[1\mapsto\text{\text{W}}]$, but
$s\not\sqsubseteq[2\mapsto\text{\text{B}}]$.

##### Scoped Functions

For many MDPs, the transition behavior and the rewards can be computed from a
combination of functions that only depend on a few state variables. Such
_scoped functions_ take a partial state as input and determine the output
inspecting only dimensions that are within their scope. The algorithm we
formalize shows how to express policy evaluation and policy improvement with
scoped functions, which avoids enumerating the state space. Scoped functions
pose a challenge for formalization, as scopes could be represented implicitly
or explicitly: we can either prove that a given function has a restricted
scope, or the function can store its own scope explicitly as data. In
Isabelle/HOL we do both: proving that a function has a certain scope only
after the function definition is more flexible, as one may derive multiple
scopes for a single function (see Listing A.2.1). However, since it is
infeasible to calculate the scope of a function in general, we use explicit
scopes in the executable version of the algorithm. An important operation on
scoped functions is instantiation with a partial state $t$. The operation
$\mathit{inst}_{t}$ applied to a scoped function returns a new scoped function
with reduced scope, where all input dimensions that $t$ provides have been
fixed to the values of $t$.

##### Transitions

The MDP provides a finite set of actions $A$ with a default action $d\in A$.
For each action, transition probabilities
$\mathit{\mathcal{P}}^{a}_{i}:\mathit{X}\to\textrm{$\mathbb{P}$}(\mathit{X}_{i})$
with
$\mathit{scope}(\mathit{\mathcal{P}}^{a}_{i})\subseteq\mathit{\mathbb{N}}_{<\mathit{\mathit{n}}}$
determine the system’s evolution. The sets $\mathit{effects}_{a}$ define the
state variables where the behavior of $a$ differs from the default action
(i.e. variables $i$ s.t.
$\mathit{\mathcal{P}}^{a}_{i}\neq\mathit{\mathcal{P}}^{d}_{i}$). The combined
transition probabilities, induced by action $a$, between two states $x$ and
$x^{\prime}$ are defined as
$\mathit{\mathcal{P}}^{a}(x,x^{\prime})\coloneqq\prod_{i<\mathit{\mathit{n}}}~{}\mathit{\mathcal{P}}^{a}_{i}(x,x^{\prime}_{i})$
(see Listing A.2.1).

##### Example

In the ring topology domain, there is an action to restart each machine.
Restarting a machine guarantees that it will work in the next step. The
default action $d$ is to do nothing. The probability that a machine works in
the next step depends on its own state and the state of the predecessor. This
means that e.g. $\mathit{scope}(\mathit{\mathcal{P}}^{d}_{2})={1,2}$. The
exact conditional probability distribution for the MDP’s evolution under this
action is shown in Figure 1. Using an explicit representation to model this
factored action, we would need, in the worst-case, $16$ transitions, each of
which with a distribution over $16$ possible successor states. This is in
contrast to four tables in the factored case. A common way to model these
scopes is using dependency graphs as shown in Figure 1.

##### Rewards

Each action comes with scoped reward functions
$\mathit{R}^{a}_{i}:\mathit{X}\to\mathit{\mathbb{R}}$, for $i<r_{a}$. The
reward for taking an action is a sum of its reward functions
$\mathit{R}^{a}(x)\coloneqq\textstyle\sum_{i<\mathit{r}_{a}}\mathit{R}^{a}_{i}(x).$
We assume that the first $r_{d}$ reward functions are the same for all
actions. Now, given a policy $\mathit{\pi}:X\to A$ we are interested in the
discounted expected total reward
$\mathit{\nu}_{\mathit{\pi}}(x)\coloneqq\mathit{\mathbb{E}}_{\omega\sim\mathcal{T}(\pi,x)}\left[\textstyle\sum_{i}\mathit{\gamma}^{i}\mathit{R}^{\mathit{\pi}(\omega_{i})}(\omega_{i})\right]$
with discount factor $\mathit{\gamma}<1$ and trace space $\mathcal{T}$. Our
goal is to achieve the optimal reward
$\mathit{\nu}^{*}(x)\coloneqq\textstyle\sup_{\mathit{\pi}\in\Pi}\mathit{\nu}_{\pi}(x)$.
For a value estimate $v:\mathit{X}\to\mathit{\mathbb{R}}$ and an action $a$,
the $Q$ function is defined as

$\mathit{Q}^{a}_{v}(x)\coloneqq\mathit{R}^{a}(x)+\mathit{\gamma}\textstyle\sum_{x^{\prime}\in\mathit{X}}\mathit{\mathcal{P}}^{a}(x,x^{\prime})\cdot
v(x^{\prime}).$

For each state, the maximum lookahead w.r.t. all actions is
$\mathit{Q}^{*}_{v}(x)$. A _greedy policy_ $\pi$ is a policy where
$\mathit{Q}^{\mathit{\pi}(x)}_{v}(x)=\mathit{Q}^{*}_{v}(x)$ for all states
$x$. The Bellman error, denoted by $\|v-Q^{\pi}_{v}\|$, is the maximum
distance between $Q^{\mathit{\pi}(x)}_{v}(x)$ and $v(x)$ for any state $x\in
X$, i.e. the $L_{\infty}$ distance.

##### Example

In our example, a factored representation of the rewards of $d$,
$\mathit{R}^{d}_{i}$, is $\mathit{R}^{d}_{i}(\text{W})=1$ and
$\mathit{R}^{d}_{i}(\text{B})=0$, for all $1\leq i\leq 4$. This gives rise to
an exponentially smaller reward function in comparison to an explicitly
represented MDP, where the reward function for $d$ would have $16$ entries.

$C_{1}$$C_{2}$$C_{3}$$C_{4}$ (a)

$X_{1}$$X_{2}$$X_{3}$$X_{4}$$X_{1}^{\prime}$$X_{2}^{\prime}$$X_{3}^{\prime}$$X_{4}^{\prime}$
(b)

$X\mid_{\mathit{scope}(\mathit{\mathcal{P}}^{d}_{1})}$ | $\mathit{\mathcal{P}}^{d}_{1}(X,X_{1}^{\prime}\mapsto\text{W})$
---|---
$X_{1}\mapsto\text{W}$, $X_{4}\mapsto\text{W}$ | 0.9
$X_{1}\mapsto\text{B}$, $X_{4}\mapsto\text{W}$ | 0.2
$X_{1}\mapsto\text{W}$, $X_{4}\mapsto\text{B}$ | 0.7
$X_{1}\mapsto\text{B}$, $X_{4}\mapsto\text{B}$ | 0.1
(c)

Figure 1: (a) Ring network of 4 connected computers. (b) Variable dependencies
of the default action. (c) Probabilities of machine $C_{1}$ working in the
next step, for every state of the machines $C_{1}$ and $C_{4}$.

###  2.3 Linear Value Functions

Even if all reward and transition functions are scoped functions, the value
function $\mathit{\nu}_{\mathit{\pi}}$ may still be unstructured, i.e.
computing $\mathit{\nu}_{\mathit{\pi}}$ might require the construction of an
exponentially big mapping, in the worst case [16]. However, the value of a
policy can be approximated as the weighted sum of $\mathit{m}$ basis functions
$h_{i}:\mathit{X}\to\mathit{\mathbb{R}}$. Given a weight $w_{i}$ for each
$h_{i}$, the value of a state can be computed as a weighted sum
$\mathit{\nu}_{w}(x)\coloneqq\textstyle\sum_{i<\mathit{m}}w_{i}h_{i}(x)$. Note
that the efficiency of the algorithm we verify here crucially depends on the
fact that $h_{i}$ is a scoped function, where
$\mathit{scope}(h_{i})\subset\mathit{X}$.

For each action choice $a$ and a basis function $h_{i}$, we can compute its
expected evaluation $\mathit{g}^{a}_{i}$ in the successor state. Since
$\mathit{g}$ does not depend on the weights, it only needs to be computed once
for a set of basis functions. In Isabelle/HOL (see Listing A.2.1), we prove
that $\mathit{g}^{a}_{i}$ has a structured representation with scope
$\Gamma^{a}_{i}\coloneqq\bigcup_{j\in\mathit{scope}(h_{i})}\mathit{scope}(\mathit{\mathcal{P}}^{a}_{j})$
as follows

$\displaystyle\mathit{g}^{a}_{i}(x)\coloneqq$
$\displaystyle\textstyle\sum{x^{\prime}\in\mathit{X}}.~{}\mathit{\mathcal{P}}^{a}(x,x^{\prime})\cdot\mathit{h}_{i}(x^{\prime})=\textstyle\sum{x^{\prime}\in\mathit{X}|_{\Gamma^{a}_{i}}}.~{}\mathit{\mathcal{P}}^{a}(x|_{\Gamma^{a}_{i}},x^{\prime})\cdot\mathit{h}_{i}(x^{\prime}).$

The second equality holds since
$\sum{x^{\prime}\in\mathit{X}|_{Y}}.~{}\mathit{\mathcal{P}}^{a}(x,x^{\prime})=1$
for any scope $Y$. This also leads to an efficient computation of the $Q$
functions: $\mathit{Q}_{w}^{a}(x)\coloneqq
Q^{a}_{\mathit{\nu}_{w}}(x)=\mathit{R}^{a}(x)+\mathit{\gamma}\textstyle\sum_{i<\mathit{m}}w_{i}\mathit{g}^{a}_{i}(x)$.

##### Example

The value functions in the ring domain can be approximately represented using
the basis function $h_{0}=1$ and one function per machine: $h_{i}=1$ if
$X_{i}=\text{W}$ and $0$ otherwise, for $1\leq i\leq 4$. Note that because
these basis functions have very limited scopes, the accuracy of the best
possible approximation of the value function using these basis will be
limited. This may be improved by including more dimensions in the scopes to
reflect dependencies between variables.

##  3 Approximate Policy Iteration

Approximate Policy Iteration (API) is a variant of policy iteration that uses
basis functions to scale to large systems [16]. Each iteration of API consists
of three parts: policy evaluation, policy improvement and Bellman error
computation. The algorithm terminates when either a timeout $t_{\mathit{max}}$
is reached, the error dips below a threshold $\epsilon$, or the weights given
to the basis functions converge. In Isabelle/HOL, the algorithm is implemented
as a function $\mathit{api}(t,\mathit{\pi},w)$ (Algorithm 1, Listing A.2.2).
It takes a time $t$, weights for the basis functions $w$ and a policy
$\mathit{\pi}$. The initial call to the algorithm is
$\mathit{api}(0,\mathit{\pi}^{0},w^{0})$ where $w_{0}=0$ and
$\mathit{\pi}^{0}$ is some greedy policy w.r.t. $w^{0}$. A single iteration of
API first uses the current policy $\mathit{\pi}$ to compute updated weights
$w^{\prime}$, then computes a new greedy policy $\mathit{\pi}^{\prime}$, and
finally determines the Bellman error $\mathit{err}$ of the updated policy. If
one of the termination conditions is met, the algorithm returns the current
iteration number, weights and policy, as well as the error and whether the
weights converged. Otherwise, $\mathit{api}$ is called again recursively with
updated inputs.

######  Algorithm 1 (Approximate Policy Iteration).

$\displaystyle\mathit{api}(t,\mathit{\pi},w)\coloneqq\textbf{if}~{}t\geq
t_{\mathit{max}}~{}\textbf{or}~{}\mathit{err}\leq\epsilon~{}\textbf{or}~{}w_{=}~{}\textbf{then}~{}(t,\mathit{\pi}^{\prime},w^{\prime},\mathit{err},w_{=})~{}\textbf{else}~{}\mathit{api}(t+1,\mathit{\pi}^{\prime},w^{\prime})$
$\displaystyle\quad\begin{aligned}
\textbf{where}~{}&w^{\prime}&&\coloneqq\mathit{upd\\_w}(\mathit{\pi})\\\
&\mathit{\pi}^{\prime}&&\coloneqq\mathit{greedy\\_\pi}(w^{\prime})\\\
&\mathit{err}&&\coloneqq\mathit{factored\\_err}(\mathit{\pi}^{\prime},w^{\prime})\\\
&w_{=}&&\coloneqq w^{\prime}=w\end{aligned}$

We structure and decouple the algorithm using locales. Conceptually, this
usage of locales is similar to using modules in programming languages. Here,
we use locales to postulate the existence of three functions
($\mathit{upd\\_w}$, $\mathit{greedy\\_\pi}$, $\mathit{factored\\_err}$) along
with their specifications:

######  Specification 1 ($\mathit{upd\\_w}$).

A decision list policy a list representation of policies where each entry
(called branch) is a pair $(t,a)$ of a partial state and an action. To select
an action in a state $x$, we search the list for the first branch where $x$ is
consistent with $t$. Now fix a decision list policy $\mathit{\pi}$, let
$w^{\prime}=\mathit{upd\\_w}(\mathit{\pi})$. Then $\mathit{\nu}_{w^{\prime}}$
is the best possible estimate of $\mathit{\nu}_{\mathit{\pi}}$:
$\|\mathit{\nu}_{w^{\prime}}-\mathit{\nu}_{\mathit{\pi}}\|=\textstyle\inf_{w}\|\mathit{\nu}_{w}-\mathit{\nu}_{\mathit{\pi}}\|.$

######  Specification 2 ($\mathit{greedy\\_\pi}$).

For all weights $w$, $\mathit{greedy\\_\pi}(w)$ is a greedy decision list
policy for $\mathit{\nu}_{w}$.

######  Specification 3 ($\mathit{factored\\_err}$).

Given weights $w$ and a greedy decision list policy $\mathit{\pi}$ for
$\mathit{\nu}_{w}$, $\mathit{factored\\_err}$ determines the Bellman error:
$\mathit{factored\\_err}(\mathit{\pi},w)=\|\mathit{Q}^{*}_{w}-\mathit{\nu}_{w}\|$.

For now, we merely state the specifications for the algorithms approximate
policy iteration builds upon, only later will we show how to implement the
specifications concretely. Viewed differently, the specifications serve as
inefficient or even non-executable implementations of algorithms. This
approach keeps the assumptions on individual parts of the algorithm explicit
and permits an easier exchange of implementations, e.g. in our developments
one may swap the LP certification algorithm for a verified LP solver
implementation. It also facilitates gradual verification of software: the
correct behavior of the software system can be proved top-down starting from
assumptions on each component. Finally, the small number of axioms improves
performance of automated proof methods in ITPs. In the following sections we
show that these specifications have efficient implementations. See Figure 2
for an overview of all components.

Within the context of the locale, assuming all specifications, we can derive
the same error bounds as presented in [16] (see Listing A.2.2). One exemplary
important result is that if the weights converge during API, then in the last
step the bellman error is equal to the approximation error. This observation
leads to the following _a posteriori_ optimality bound:

######  Theorem 1.

Let
$api(w_{0},\mathit{\pi}_{0})=(t^{\prime},\mathit{\pi},w,\mathit{err},\mathit{True})$.
Then $(1-\mathit{\gamma})\|\mathit{\nu}^{*}-\mathit{\nu}_{w}\|\leq
2\mathit{\gamma}\cdot\mathit{err}$.

ptApproximate Policy IterationAlgo 1: Spec 1, 2, 3 $\Rightarrow$ Thm 1
ptDecision List PolicyAlgo 2: Spec 2 ptValue DeterminationSpec 6 $\Rightarrow$
Spec 1 ptLP BranchAlgo 6: Spec 7 $\Rightarrow$ Spec 6 ptFactored LPSpec 7 ptLP
Certification ptFactored Bellman ErrorAlgo 3: Spec 4 $\Rightarrow$ Spec 3
ptBranch ErrorAlgo 4: Spec 5 $\Rightarrow$ Spec 4 ptVariable EliminationAlgo
5: Spec 5 Figure 2:  Simplified Locale Hierarchy

##  4 Policy Improvement

Given weights $w$, the policy improvement phase determines a greedy policy
w.r.t. $\nu_{w}$. The policy takes the form of a decision list, so each
element of the list stores a partial state and an action. The main idea for an
efficient computation is that only actions better than the default action need
to be considered. This notion is made precise by the bonus function
$\delta_{a}$ (Algorithm 2, Listing A.2.3) with scope
$\mathit{\mathbf{T}}_{a}$. Originally, $\Gamma_{d}$ was not part of
$\mathit{\mathbf{T}}_{a}$[16], which we assume to be an oversight in the
definition, since the behavior of the default action does influence the bonus.
Unless components cancel out, the scope of a function difference is the union
of the scopes of both arguments. Finally, we concatenate the branches
$\mathit{\pi}_{a}$ for every action but $d$, add the default action as a
fallback and sort the decision list policy by decreasing bonus. Here, the
empty map with no entries is $\bot_{M}$, and
$\mathit{Dom}(\bot_{M})=\emptyset$. We can show that $\mathit{greedy\\_\pi}$
satisfies Specification 2, since action selection proceeds by decreasing
bonus.

######  Algorithm 2 (Decision List Policy).

$\displaystyle\mathit{greedy\\_\pi}\coloneqq\mathit{sort\\_\pi}((\bot_{M},d,0)~{}\mathit{::}~{}\mathit{concat}([\mathit{\pi}_{a}~{}|~{}a\in
A-\\{d\\}])$ $\displaystyle\quad\begin{aligned}
\textbf{where}~{}\mathit{\pi}_{a}&\coloneqq[(x,a,\mathit{\delta}_{a}(x))\mid\mathit{\delta}_{a}(x)>0,x\in
X_{\mid\mathit{\mathbf{T}}_{a}}]&\quad&\text{(Positive bonus, states with
domain $\mathit{\mathbf{T}}_{a}$)}\\\ \mathit{\delta}_{a}&\coloneqq
Q^{a}_{w}-Q^{d}_{w}&\quad&\text{(Bonus function)}\\\
\mathit{\mathbf{T}}_{a}&\coloneqq\mathit{scope}(R^{a})\cup\textstyle\bigcup_{i\in\mathbf{I}_{a}}\Gamma^{a}_{i}\cup\Gamma^{d}_{i}&\quad&\text{(Scope
of the bonus function)}\\\
\mathbf{I}_{a}&\coloneqq\\{i<m\mid\mathit{effects}_{a}\cap\mathit{scope}(h_{i})\neq\emptyset\\}&\qquad&\text{(Basis
functions relevant to the action)}\end{aligned}$

##  5 Factored Bellman Error

The Bellman error $\|\mathit{Q}^{*}_{w}-\mathit{\nu}_{w}\|$ is an indicator of
the degree of optimality of a policy. An inefficient computation would
enumerate every state, and return the maximum error. However, for a decision
list policy, we can compute the error incurred by each branch separately. The
total error then equals the maximum error of any branch (Algorithm 3, Listing
A.2.4). For now, assume that we have a function $\mathit{branch\\_err}$ that
computes the error for a single branch, i.e. the maximum error for any state
that selects the respective branch. These are all states that are consistent
with the current branch $t$, but did not match any prior branch
$t^{\prime}\in\mathit{ts}$ of the policy, formally
$\mathit{X_{(t,ts)}}\coloneqq\\{x\in\mathit{X}.~{}t\sqsubseteq x\land\forall
t^{\prime}\in\mathit{ts}.~{}t^{\prime}\not\sqsubseteq x\\}$. Hence we also
need to past the prefix of the decision list policy to the function that
computes the error of a branch. Note that if the branch is selected by no
state its error is $-\infty$. We show that if Specification 4 is met and
$\mathit{\pi}$ is a greedy policy w.r.t. $w$, then $\mathit{factored\\_err}$
satisfies Specification 3 (see Listing A.2.4).

######  Algorithm 3 (Factored Bellman Error).

$\mathit{factored\\_err}\coloneqq\textstyle\sup_{i<\mathit{len}(p)}\mathit{branch\\_err}(p_{i},\mathit{map}(\mathit{fst},p_{:i}))$

######  Specification 4 ($\mathit{branch\\_err}$).

Given a prefix of a policy $\pi$ with the current branch $(t,a)$ and a list of
partial states $\mathit{ts}$ from prior branches,
$\mathit{branch\\_err}(t,a,\mathit{ts})=\textstyle\sup_{x\in\mathit{X_{(t,ts)}}}|Q^{a}_{w}(x)-\mathit{\nu}_{w}(x)|$.

##### Branch Error

Take a branch, $(t,a)$, of the policy, for a partial state $t$ an action $a$.
The states of all prior branches form the list of partial states
$\mathit{ts}$. To find the Bellman error of the current branch, we need to
maximize $|\mathit{Q}^{a}_{w}(x)-\mathit{\nu}_{w}(x)|$ w.r.t. states
$x\in\mathit{X_{(t,ts)}}$. Note that for all states $x$

$\displaystyle Q^{a}_{w}(x)-\mathit{\nu}_{w}(x)=$
$\displaystyle\textstyle\sum_{i<r_{a}}\mathit{R}^{a}_{i}(x)+\textstyle\sum_{i<\mathit{m}}w_{i}(h_{i}-\mathit{\gamma}g_{i})(x)$
$\displaystyle{}={}$ $\displaystyle\textstyle\sum
f\in[\mathit{R}^{a}_{0},\ldots,\mathit{R}^{a}_{\mathit{r}_{a}},w_{0}(h_{0}-\mathit{\gamma}g_{0}),\ldots,w_{m}(h_{m}-\mathit{\gamma}g_{m})].~{}f(x).$
(1)

If we had an algorithm to efficiently compute the maximum sum of scoped
functions, we could determine the error of a branch. Again, we specify an
algorithm $\mathit{max_{\Sigma}}$ to do exactly that (Specification 5). In
Algorithm 4 (see Listing A.2.4), we call $\mathit{max_{\Sigma}}$ with the
functions from section 5. To restrict the maximization to the states
$\mathit{X_{(t,ts)}}$, we instantiate all functions with the partial state
$t$. Additionally, we define the functions $\mathcal{I}^{\prime}$ that
evaluate to $-\infty$ on states that select a different branch of the policy.
Hence, these states are ignored in the error computation. We also apply
$\mathit{max_{\Sigma}}$ to the negated functions to compute the absolute value
of the error. Finally, we formally prove that $\mathit{branch\\_err}$
satisfies Specification 4.

######  Specification 5 (Variable Elimination).

For scoped functions $\mathit{fs}$,
$\mathit{max_{\Sigma}}(\mathit{fs})=\sup_{x\in\mathit{X}}\textstyle\sum_{f\in\mathit{fs}}f(x)$.

######  Algorithm 4 (Branch Error).

$\displaystyle\mathit{branch\\_err}\coloneqq\max(\mathit{max_{\Sigma}}(\mathit{fs}~{}\mathit{\cdot}~{}\mathcal{I}^{\prime}),\mathit{max_{\Sigma}}(-\mathit{fs}~{}\mathit{\cdot}~{}\mathcal{I}^{\prime}))$
$\displaystyle\quad\begin{aligned}
\textbf{where}~{}&\mathit{rs}&&\coloneqq[\mathit{R}^{a}_{0},\ldots,\mathit{R}^{a}_{\mathit{r}_{a}}]\\\
&\mathit{ws}&&\coloneqq[w_{0}(\mathit{h}_{0}-\mathit{\gamma}\mathit{g}_{0}),\ldots,w_{\mathit{m}}(\mathit{h}_{\mathit{m}}-\mathit{\gamma}\mathit{g}_{\mathit{m}})]\\\
&\mathit{fs}&&\coloneqq\mathit{map}(\mathit{inst}_{t},\mathit{rs}~{}\mathit{\cdot}~{}\mathit{ws})\end{aligned}$

##### Variable Elimination

The specification for $\mathit{max_{\Sigma}}$ can be efficiently implemented
with a variable elimination algorithm (Algorithm 5, Listing A.2.4). In each
iteration, the algorithm selects a dimension of the state space, collects all
functions that depend on this dimension in a set $E$. It then creates a new
function $e$ that maximizes all functions in $E$ over that dimension of the
state space. The algorithm keeps track of a set of functions to maximize and
the number of the current iteration. Since the number of operations performed
by the algorithm varies greatly with the elimination order, the variables can
be reordered with a bijection
$\mathit{\mathit{\mathcal{O}}}:\mathit{\mathbb{N}}_{<\mathit{\mathit{n}}}\to\mathit{\mathbb{N}}_{<\mathit{\mathit{n}}}$.
We formally prove that at any point during execution, the maximum of the
current list of functions stays constant, thus $\mathit{max_{\Sigma}}$ meets
Specification 5.

######  Algorithm 5 (Variable Elimination).

$\displaystyle\mathit{max\\_step}(i,\mathit{fs})\coloneqq(i+1,e~{}\mathit{::}~{}E^{\prime})$
$\displaystyle\quad\begin{aligned}
\textbf{where}~{}~{}&(E,E^{\prime})\coloneqq\mathit{partition}(f\mapsto\mathit{\mathit{\mathcal{O}}}(i)\in\mathit{scope}(f),\mathit{fs})\
\ \text{(Partition based on
$\mathit{\mathit{\mathcal{O}}}(i)\in\mathit{scope}(f)$?)}\\\
&e\coloneqq(x\mapsto\textstyle\max_{y\in\mathit{X}_{\mathit{\mathit{\mathcal{O}}}(i)}}\textstyle\sum_{f\in
E}f(x_{\mathit{\mathit{\mathcal{O}}}(i)\mapsto y}),\textstyle\bigcup_{f\in
E}\mathit{scope}(f)-\\{\mathit{\mathit{\mathcal{O}}}(i)\\})\end{aligned}$
$\displaystyle\mathit{max_{\Sigma}}(\mathit{fs})\coloneqq\textstyle\sum_{f\in\mathit{fs^{\prime}}}f(\mathit{\bot_{M}})$
$\displaystyle\quad\begin{aligned}
\textbf{where}~{}&(\\_,\mathit{fs^{\prime}})\coloneqq\mathit{max\\_step}^{\mathit{\mathit{n}}}(0,\mathit{fs})\end{aligned}$

##  6 Value Determination

After a new candidate policy is found, it is evaluated in the value
determination phase. Since the exact value function cannot be represented in
general as a linear combination of the basis functions, the aim is to find
weights for the basis functions with the smallest approximation error (see
Specification 1). This problem can be expressed using linear programming (LP).
The structure of the algorithm that finds optimal weights is analogous to the
factored Bellman error computation. For each branch of the policy, we generate
a set of LP constraints (according to Specification 6, see Listing A.3) that
expresses the approximation error incurred by this branch. The union of all
constraints is then

$\displaystyle\mathit{weight\\_lp}\coloneqq\textstyle\bigcup_{i<\mathit{len}(p)}\mathit{branch\\_lp}(p_{i},\mathit{map}(\mathit{fst},p_{:i})).$

Variables of the LP are the approximation error $\phi$ and the weights $w$.
The LP is optimized for minimal error $\phi$, the values of the variables $w$
in an optimal solution determine the new weights. Given a set of LP
constraints $\mathit{cs}$, $\langle\mathit{cs}\rangle_{\mathit{LP}}$ denotes
the set of feasible solutions. We show that for any optimal solution
$(\phi^{*},w^{*})\in\langle\mathit{weight\\_lp}\rangle_{\mathit{LP}}$, setting
$\mathit{upd\\_w}\coloneqq w^{*}$ satisfies Specification 1. We also formally
prove that the LP always has an optimal solution, since the set of potentially
optimal solutions is compact.

######  Specification 6 ($\mathit{branch\\_lp}$).

Given a partial state $t$, an action $a$, a list of partial states
$\mathit{ts}$, $\mathit{branch\\_lp}$ constructs an LP that minimizes the
approximation error over the states $\mathit{X_{(t,ts)}}$:

$\displaystyle(\phi,w)\in\langle\mathit{branch\\_lp}(t,a,ts)\rangle_{\mathit{LP}}\longleftrightarrow\forall
x\in\mathit{X_{(t,ts)}}.~{}\phi\geq|\mathit{Q}^{a}_{w}(x)-\mathit{\nu}_{w}(x)|.$

##### LPs for Branches

For each branch of the policy, we proceed similarly to the Bellman error
computation: we create two constraint sets, for positive and negative errors
respectively (Algorithm 6, Listing A.3). We omit the scopes for brevity. At
this level, we make use of another algorithm $\mathit{min\\_lp}$. Its first
input $C$ is a list of $m$ scoped functions, the second input $b$ is another
list of scoped functions. Now $\mathit{min\\_lp}(C,b)$ creates an LP that
minimizes $Cw-b$ w.r.t. $w$ over all states (see Specification 7). We can show
that $\mathit{branch\\_lp}$ fulfills Specification 6, since the feasible
solutions are upper bounds of the smallest approximation error. The derivation
of $C$ and $b$ follows the Bellman error computation. Interested readers
should consult the formal proofs.

######  Algorithm 6 (Branch LP).

$\displaystyle\mathit{branch\\_lp}\coloneqq\mathit{min\\_lp}(C,-b~{}\mathit{\cdot}~{}\mathcal{I}^{\prime})\cup\mathit{min\\_lp}(-C,b~{}\mathit{\cdot}~{}\mathcal{I}^{\prime})$
$\displaystyle\quad\begin{aligned} \textbf{where}~{}&b&&\coloneqq
map(\mathit{inst}_{t},[\mathit{R}^{a}_{i}~{}|~{}i<\mathit{r}_{a}])&&\hfill\qquad\qquad\qquad\qquad\qquad\qquad\text{($Cw-b~{}\widehat{=}~{}\mathit{Q}^{a}_{w}-\mathit{\nu}_{w}$)}\\\
&C&&\coloneqq
map(\mathit{inst}_{t},[h_{i}-\mathit{\gamma}g_{i}^{a}~{}|~{}i<\mathit{m}])\end{aligned}$

######  Specification 7 ($\mathit{min\\_lp}$).

$\mathit{min\\_lp}(C,b)$ generates an LP that minimizes $Cw-b$:

$\displaystyle(\phi,w)\in\langle\mathit{min\\_lp}(C,b)\rangle_{\mathit{LP}}\longleftrightarrow\forall
x\in\mathit{X}.~{}\phi\geq\textstyle\sum_{i<\mathit{len}(C)}w_{i}C_{i}(x)+\textstyle\sum_{i<\mathit{len}(b)}b_{i}(x).$

##### Factored LP Construction

The algorithm $\mathit{min\\_lp}$ resembles $\mathit{max_{\Sigma}}$, so we
only point out the challenges encountered during verification. Full details
can be found in the formalization (see Listing A.3). There are two significant
modifications we made to the algorithm to make verification feasible. First,
the original algorithm may create equality constraints that constrain
variables to $-\infty$. Since these constraints are not supported by the LP
solvers we use, we formally prove that one can modify the algorithm to omit
such constraints. Second, the combination of LP constraints in the definition
of $\mathit{weight\\_lp}$ requires some care, to avoid interactions between LP
variables created in different branches. The $\mathit{min\\_lp}$ algorithm
creates new (private) LP variables, we need to make sure that the variables
have distinct names for each branch. This issue was not discussed by Guestrin
et al.. The problem can be solved by adding a tag to each generated variable.
The tags contain $t$, $a$, and a flag to differentiate the two invocations of
$\mathit{min\\_lp}$ in each branch. For distinct tags $p$, $p^{\prime}$ we can
then show: the solutions to the union of two constraint sets are equivalent to
the intersection of the solution spaces of the individual constraint sets
(concerning $\phi$ and $w$):

$\displaystyle\langle\mathit{min\\_lp}(p,C,b)\cup\mathit{min\\_lp}(p^{\prime},C^{\prime},b^{\prime})\rangle_{\mathit{LP}}=\langle\mathit{min\\_lp}(p,C,b)\rangle_{\mathit{LP}}\cap\langle\mathit{min\\_lp}(p^{\prime},C^{\prime},b^{\prime})\rangle_{\mathit{LP}}.$

We show that $\mathit{min\\_lp}$ creates an LP that is equivalent to the
explicit (potentially exponentially larger) LP that creates a constraint for
each state of the MDP. It immediately follows that $\mathit{min\\_lp}$
satisfies Specification 7. This completes the correctness proof of the
abstract version of approximate policy iteration. Due to our approach to the
verification of algorithms using loosely coupled locales, $\mathit{min\\_lp}$
and $\mathit{max_{\Sigma}}$ are not tied to MDPs and are thereby reusable
components.

##  7 Code Generation

We now discuss the process of deriving a verified efficiently executable
version from the verified abstract algorithm discussed above. To do so, we
follow the methodology of program refinement [37], where one starts with an
abstract, potentially unexecutable, version of the algorithm and verifies it.
Then one devises more optimised versions of the algorithm, and only proves the
optimizations correct in this latter step, thus separating mathematical
reasoning from implementation specific reasoning. This approach was used in
most successful algorithm verification efforts [24, 15, 23]. In this work, the
three most important stages are the initial abstract algorithm, an
implementation with abstract data structures, and finally an implementation
with concretized data structures. As a last step, we export verified code in
the programming language Scala for the approximate policy iteration algorithm.

###  7.1 Refinement using Locales

Our implementation of step-wise refinement is based on Isabelle/HOL locales.
For each locale of the abstract algorithm, we define a corresponding locale
where we define the executable version of the algorithm. Finally, we relate
the abstract version of the MDP to the concrete version. For each definition,
we then show that corresponding inputs lead to corresponding outputs, i.e. our
abstract algorithm and the implementation behave the same. At this point in
the refinement, data structures remain abstract APIs, with the concrete
implementations chosen only later. We use the data structures provided by the
Isabelle Collections Framework [25] for code generation, that we extend with
our own data structure for scoped functions, represented as a pair of a
function and a set for its scope. The data structure also provides an
operation to evaluate a function on all of its scope for memoization.

###  7.2 Certification of Linear Programming Solutions

Nodes | Constrs | Vars | $\mathbf{t(s)}$ | $\mathbf{t_{LP}(s)}$
---|---|---|---|---
1 | 74 | 41 | 0.27 | 0.02
3 | 1258 | 693 | 0.89 | 0.15
5 | 4378 | 2455 | 1.89 | 0.46
7 | 9418 | 5305 | 3.78 | 0.98
9 | 16378 | 9243 | 6.74 | 1.82
11 | 25258 | 14269 | 12.44 | 3.36
13 | 36058 | 20383 | 20.95 | 5.31
15 | 48778 | 27585 | 34.69 | 8.67
17 | 63418 | 35875 | 58.30 | 16.86
19 | 79978 | 45253 | 92.19 | 30.25
Table 1: Evaluation on the ring domain, the columns denote the number of
variables (number of states grows exponentially), the number of LP constraints
and variables generated in the last iteration. The last two columns give tlhe
total running time and the time spent in the LP solver.

An implementation of approximate policy iteration depends on efficient LP
solvers. In our verified implementation, we use precise but unverified LP
solvers and certify their results. This avoids implementing a verified,
optimized LP solver but retains formal guarantees. The tradeoff here is that
the unverified LP solver may now return solutions that cannot be certified. At
the cost of performance, it would also be possible to connect the
formalization to an existing executable simplex implementation for
Isabelle/HOL [34]. To achieve formal guarantees, the LP has to be solved
exactly, i.e. using rational numbers. Two potential candidates for precise LP
solvers are QSopt_ex [4] and SoPlex [9]. For larger LPs, the performance of
SoPlex was more reliable in our setting.

We certify optimality using the dual solution and the strong duality of linear
programming. In our formalization we also formally prove that infeasibility
and unboundedness can be certified similarly using farkas certificates or
unbounded rays. The linear program is first preprocessed to a standard form:
variable bounds and equality constraints are reduced to only inequality
constraints. The constraints of the resulting LP are in the standard form
$Ax\leq b$, with no restrictions on $x$. During benchmarking of the
certification process, the normalization operation usually applied to
rationals after each operation proved to be very costly. For certification, we
represent rational numbers as pairs that are never normalized, which leads to
faster certificate checking in our experiments.

The exported Scala program takes an arbitrary function from linear programs to
their solutions as input. If this function returns invalid solutions, they are
rejected by the certificate checker, so there are no assumptions we need to
place on the LP solver. However, there is the implicit assumption that the LP
solver is deterministic. Since we are working in a fragment of a functional
programming language, calling the LP solver twice on the same problem should
lead to the same solution. In theory, a nondeterministic LP solver could be
misused to lead to inconsistencies. As we do not compare LP solutions in our
algorithm and SoPlex is actually deterministic, this problem does not impact
our verified software. A possible more general solution to the problem could
be the use of memoization or to explicitly model the nondeterminism using
monads.

###  7.3 Experimental Evaluation

We show the practicality of our verified implementation by applying it to the
ring topology domain. The run our implementation on an Intel i7-11800H CPU and
we set the discount factor to $0.9$ in all our experiments. We do not reorder
variables. In all runs, the weights converged after 5 iterations. The results
of the experiments (Table 1, for full table see Table 2) show that the
algorithm can deal with systems of half a million states and 20 actions. For
larger systems, the precise mode of the LP solver Soplex cannot find a
rational solution. The experiment shows that our implementation of linear
programming certification can process linear programs with tens of thousands
of constraints.

##  8 Discussion

Our work combines and integrates a wide range of tools and formalization
efforts to produce a verified implementation of approximate policy iteration.
In Isabelle/HOL, we build on the formal libraries for LPs [35], MDPs [19],
linear algebra, analysis [21], probability theory [20], and the collections
framework [25]. Furthermore, we certify the results computed by precise LP
solvers. In total, our development is comprised of approximately 20,000 lines
of code, two thirds of which are concerning code generation. We show how to
facilitate locales, powerful automation in Isabelle/HOL, and certification to
develop complex formally verified software. We also make the case that the
process of algorithm verification provides a detailed understanding of the
algorithm: all hidden assumptions are made explicit, while we modularize the
algorithm into components with precisely specified behaviour.

The methodology presented here can be applied to a large number of algorithms
since the algorithm we verified is a seminal algorithm for solving factored
MDPs, combining a large number of concepts that are widely used in many
contexts, like AI [13, 30, 14] and model checking [18, 12].

Applications of interactive theorem provers in verification and formalizing
mathematics have been recently attracting a lot of attention [6, 5, 27]. In
most applications, especially in computer science [15, 24] and AI [7, 33], the
emphasis is on the difficulty of the proofs, whether that is due to many cases
or complex constructions, etc., and how theorem proving helped find mistakes
or find missing cases in the proofs. A distinct feature of this work is that
its complexity comes from the large number of concepts it combines, shown in
Fig. 2, which is a more prevalent issue in formalizing pure mathematics. Our
project has contributed a better restructuring of the algorithm and untangling
of the different concepts, leading to better understandability.

Multiple directions can be considered to extend and build on our work. An
alternative to using exact LP solvers is to use their highly optimized
floating point counterparts that do not calculate exact solutions. In this
case, one needs to certify the error bound derived from a dual solution to the
LP. Then the error bound has to be incorporated in the error analysis of the
algorithm to obtain formal guarantees. Moreover, we may initialize the
algorithm with weights computed by an unverified, floating point
implementation of API to reduce number of iterations performed by the verified
implementation. In Isabelle/HOL, some of the correctness proofs for code
generation can be automated using tools to transfer theorems. Furthermore, the
ergonomics for instantiating locales and inheriting from other locales could
be improved for situations with many locale parameters.

##### Related Work

Several formal treatments of MDPs have been developed in the theorem provers
Isabelle/HOL [19, 32, 11, 17] and Coq [36]. All developments verify algorithms
for explicitly represented MDPs, which limits their practical applicability to
solve large MDPs. We adapt and integrate the implementation of Schäffeler and
Abdulaziz with our formalization. The algorithm we verified in Isabelle/HOL
was first presented by [16]. An elementary approach to Linear Programming
certification certification with Isabelle/HOL has been done before as part of
the Flyspeck project [29], where the feasibility of a solution to a linear
program is checked by Isabelle/HOL. The work presents a method that uses dual
solutions produced by floating-point LP solvers to find bounds on the
objective value of an LP. The tool Marabou for neural network verfication uses
Farkas vectors to produce proofs for its results [22].

## References

  * Abdulaziz and Koller [2022] M. Abdulaziz and L. Koller.  Formal semantics and formally verified validation for temporal planning.  In _The 36th AAAI Conference on Artificial Intelligence (AAAI)_ , 2022.  doi: 10.1609/aaai.v36i9.21197.
  * Abdulaziz and Kurz [2023] M. Abdulaziz and F. Kurz.  Formally verified SAT-Based AI planning.  In _The 37th AAAI Conference on Artificial Intelligence (AAAI)_ , 2023.  doi: 10.1609/aaai.v37i12.26714.
  * Abdulaziz and Lammich [2018] M. Abdulaziz and P. Lammich.  A Formally Verified Validator for Classical Planning Problems and Solutions.  In _The 30th International Conference on Tools with Artificial Intelligence (ICTAI)_ , 2018.  doi: 10.1109/ICTAI.2018.00079.
  * Applegate et al. [2007] D. L. Applegate, W. Cook, S. Dash, and D. G. Espinoza.  Exact solutions to linear programming problems.  _Operations Research Letters_ , Nov. 2007.  ISSN 01676377.  doi: 10.1016/j.orl.2006.12.010.
  * Avigad [2023] J. Avigad.  Mathematics and the formal turn, Nov. 2023.
  * Avigad and Harrison [2014] J. Avigad and J. Harrison.  Formally verified mathematics.  _Commun. ACM_ , 2014.  doi: 10.1145/2591012.
  * Bagnall and Stewart [2019] A. Bagnall and G. Stewart.  Certifying the True Error: Machine Learning in Coq with Verified Generalization Guarantees.  In _The 33rd AAAI Conference on Artificial Intelligence (AAAI)_ , 2019.  doi: 10.1609/aaai.v33i01.33012662.
  * Ballarin [2014] C. Ballarin.  Locales: A Module System for Mathematical Theories.  _J. Autom. Reason._ , 2014.  doi: 10.1007/s10817-013-9284-7.
  * Bestuzheva et al. [2023] K. Bestuzheva, M. Besançon, W.-K. Chen, A. Chmiela, T. Donkiewicz, J. van Doornmalen, L. Eifler, O. Gaul, G. Gamrath, A. Gleixner, L. Gottwald, C. Graczyk, K. Halbig, A. Hoen, C. Hojny, R. van der Hulst, T. Koch, M. Lübbecke, S. J. Maher, F. Matter, E. Mühmer, B. Müller, M. E. Pfetsch, D. Rehfeldt, S. Schlein, F. Schlösser, F. Serrano, Y. Shinano, B. Sofranac, M. Turner, S. Vigerske, F. Wegscheider, P. Wellner, D. Weninger, and J. Witzig.  Enabling Research through the SCIP Optimization Suite 8.0.  _ACM Trans. Math. Softw._ , June 2023.  ISSN 0098-3500.  doi: 10.1145/3585516.
  * Budde et al. [2017] C. E. Budde, C. Dehnert, E. M. Hahn, A. Hartmanns, S. Junges, and A. Turrini.  JANI: Quantitative Model and Tool Interaction.  In _The 23rd International Conference on Tools and Algorithms for the Construction and Analysis of Systems_ , 2017\.  doi: 10.1007/978-3-662-54580-5_9.
  * Chevallier and Fleuriot [2021] M. Chevallier and J. D. Fleuriot.  Formalising the Foundations of Discrete Reinforcement Learning in Isabelle/HOL.  _CoRR_ , 2021.  URL https://arxiv.org/abs/2112.05996.
  * Dehnert et al. [2017] C. Dehnert, S. Junges, J.-P. Katoen, and M. Volk.  A Storm is Coming: A Modern Probabilistic Model Checker.  In _The 29th International Conference on Computer Aided Verification_ , 2017.  doi: 10.1007/978-3-319-63390-9_31.
  * Delgado et al. [2011] K. V. Delgado, S. Sanner, and L. N. De Barros.  Efficient solutions to factored MDPs with imprecise transition probabilities.  _Artificial Intelligence_ , 2011.  doi: 10.1016/j.artint.2011.01.001.
  * Deng et al. [2022] Z. Deng, S. Devic, and B. Juba.  Polynomial Time Reinforcement Learning in Factored State MDPs with Linear Value Functions.  In _The 25th International Conference on Artificial Intelligence and Statistics_ , May 2022.  URL https://proceedings.mlr.press/v151/deng22c.html.
  * Esparza et al. [2013] J. Esparza, P. Lammich, R. Neumann, T. Nipkow, A. Schimpf, and J.-G. Smaus.  A Fully Verified Executable LTL Model Checker.  In _25th International Conference on Computer Aided Verification (CAV)_ , 2013.  doi: 10.1007/978-3-642-39799-8_31.
  * Guestrin et al. [2003] C. Guestrin, D. Koller, R. Parr, and S. Venkataraman.  Efficient solution algorithms for factored MDPs.  _J. Artif. Intell. Res._ , 2003.  doi: 10.1613/jair.1000.
  * Hartmanns et al. [2023] A. Hartmanns, B. Kohlen, and P. Lammich.  Fast Verified SCCs for Probabilistic Model Checking.  In _The 21st International Symposium on Automated Technology for Verification and Analysis (ATVA)_ , 2023.  doi: 10.1007/978-3-031-45329-8_9.
  * Hinton et al. [2006] A. Hinton, M. Z. Kwiatkowska, G. Norman, and D. Parker.  PRISM: A Tool for Automatic Verification of Probabilistic Systems.  In _The 22nd International Conference on Tools and Algorithms for the Construction and Analysis of Systems (TACAS)_ , 2006.  doi: 10.1007/11691372_29.
  * Hölzl [2017] J. Hölzl.  Markov Chains and Markov Decision Processes in Isabelle/HOL.  _J. Autom. Reason._ , 2017.  doi: 10.1007/s10817-016-9401-5.
  * Hölzl and Heller [2011] J. Hölzl and A. Heller.  Three Chapters of Measure Theory in Isabelle/HOL.  In _2nd International Conference on Interactive Theorem Proving (ITP)_ , 2011.  doi: 10.1007/978-3-642-22863-6_12.
  * Hölzl et al. [2013] J. Hölzl, F. Immler, and B. Huffman.  Type Classes and Filters for Mathematical Analysis in Isabelle/HOL.  In _4th International Conference on Interactive Theorem Proving (ITP)_ , 2013.  doi: 10.1007/978-3-642-39634-2_21.
  * Isac et al. [2022] O. Isac, C. Barrett, M. Zhang, and G. Katz.  Neural Network Verification with Proof Production.  In _The 22nd Conference on Formal Methods in Computer-Aided Design (FMCAD)_ , Oct. 2022.  doi: 10.34727/2022/isbn.978-3-85448-053-2_9.
  * Kanav et al. [2014] S. Kanav, P. Lammich, and A. Popescu.  A Conference Management System with Verified Document Confidentiality.  In _The 26th International Conference on Computer Aided Verification (CAV)_ , 2014.  doi: 10.1007/978-3-319-08867-9_11.
  * Klein et al. [2009] G. Klein, K. Elphinstone, G. Heiser, J. Andronick, D. Cock, P. Derrin, D. Elkaduwe, K. Engelhardt, R. Kolanski, M. Norrish, T. Sewell, H. Tuch, and S. Winwood.  seL4: Formal verification of an OS kernel.  In _22nd ACM Symposium on Operating Systems Principles 2009 (SOSP)_ , 2009.  doi: 10.1145/1629575.1629596.
  * Lammich and Lochbihler [2019] P. Lammich and A. Lochbihler.  Automatic Refinement to Efficient Data Structures: A Comparison of Two Approaches.  _J Autom Reasoning_ , June 2019.  ISSN 1573-0670.  doi: 10.1007/s10817-018-9461-9.
  * Leroy [2009] X. Leroy.  Formal verification of a realistic compiler.  _Commun. ACM_ , 2009.  doi: 10.1145/1538788.1538814.
  * Massot [N/A] P. Massot.  Why formalize mathematics?  N/A.
  * Nipkow et al. [2002] T. Nipkow, L. C. Paulson, and M. Wenzel.  _Isabelle/HOL - A Proof Assistant for Higher-Order Logic_.  2002\.
  * Obua and Nipkow [2009] S. Obua and T. Nipkow.  Flyspeck II: The basic linear programs.  _Ann Math Artif Intell_ , Aug. 2009.  ISSN 1573-7470.  doi: 10.1007/s10472-009-9168-z.
  * [30] I. Osband and B. V. Roy.  Near-optimal Reinforcement Learning in Factored MDPs.
  * [31] S. Sanner.  Relational Dynamic Influence Diagram Language (RDDL): Language Description.
  * Schäffeler and Abdulaziz [2023] M. Schäffeler and M. Abdulaziz.  Formally Verified Solution Methods for Infinite-Horizon Markov Decision Processes.  In _The 37th AAAI Conference on Artificial Intelligence (AAAI)_ , 2023.  doi: 10.1609/aaai.v37i12.26759.
  * Selsam et al. [2017] D. Selsam, P. Liang, and D. L. Dill.  Developing Bug-Free Machine Learning Systems With Formal Mathematics.  In _The 34th International Conference on Machine Learning (ICML)_ , 2017.  URL http://proceedings.mlr.press/v70/selsam17a.html.
  * Spasić and Marić [2012] M. Spasić and F. Marić.  Formalization of Incremental Simplex Algorithm by Stepwise Refinement.  In _The 18th International Symposium on Formal Methods_ , 2012\.  doi: 10.1007/978-3-642-32759-9_35.
  * Thiemann [2022] R. Thiemann.  Duality of Linear Programming.  _Arch. Formal Proofs_ , 2022.  URL https://www.isa-afp.org/entries/LP_Duality.html.
  * Vajjha et al. [2021] K. Vajjha, A. Shinnar, B. M. Trager, V. Pestun, and N. Fulton.  CertRL: Formalizing convergence proofs for value and policy iteration in Coq.  In _The 10th International Conference on Certified Programs and Proofs (CPP)_ , 2021.  doi: 10.1145/3437992.3439927.
  * Wirth [1971] N. Wirth.  Program Development by Stepwise Refinement.  _Commun. ACM_ , 1971.  doi: 10.1145/362575.362577.
  * Younes and Littman [2004] H. L. Younes and M. L. Littman.  PPDDL1. 0: The language for the probabilistic part of IPC-4.  In _Proc. International Planning Competition_ , 2004.

##  Appendix A Appendix / supplemental material

###  A.1 Experiments

Computers | Constraints | Variables | Time (s) | Time in LP (s)
---|---|---|---|---
1 | 74 | 41 | 0.27 | 0.02
2 | 498 | 268 | 0.58 | 0.08
3 | 1258 | 693 | 0.89 | 0.15
4 | 2578 | 1438 | 1.34 | 0.28
5 | 4378 | 2455 | 1.89 | 0.46
6 | 6658 | 3744 | 2.61 | 0.69
7 | 9418 | 5305 | 3.78 | 0.98
8 | 12658 | 7138 | 5.08 | 1.31
9 | 16378 | 9243 | 6.74 | 1.82
10 | 20578 | 11620 | 8.89 | 2.45
11 | 25258 | 14269 | 12.44 | 3.36
12 | 30418 | 17190 | 15.75 | 4.19
13 | 36058 | 20383 | 20.95 | 5.31
14 | 42178 | 23848 | 26.92 | 6.66
15 | 48778 | 27585 | 34.69 | 8.67
16 | 55858 | 31594 | 44.42 | 10.79
17 | 63418 | 35875 | 58.30 | 16.86
18 | 71458 | 40428 | n.a. | n.a.
19 | 79978 | 45253 | 92.19 | 30.25
20 | 88978 | 50350 | n.a. | n.a.
Table 2: Evaluation on the ring domain, the columns denote the number of
variables (number of states grows exponentially), the number of LP constraints
and variables generated in the last iteration. The last two columns give the
total running time and the time spent in the LP solver.

###  A.2 Isabelle/HOL Snippets

Here we show listings of the most important parts of the formalization.

####  A.2.1 Factored MDPs

Listing 1: Definition of Factored MDP

* _locale_ Factored_MDP_Consts =

 _fixes_

dims :: nat  _and_ — Number of state variables

doms ::  _‹_ nat  _$\Rightarrow$_ ’x :: countable set _›_  _and_ — Domain of
each state variable

reward_dim ::  _‹_ ’a ::countable  _$\Rightarrow$_ nat _›_  _and_

— number of rewards for action

rewards ::  _‹_ ’a  _$\Rightarrow$_ nat  _$\Rightarrow$_ ’x state
_$\Rightarrow$_ real _›_  _and_

— rewards a i x: reward i for action a in state x

reward_scope ::  _‹_ ’a  _$\Rightarrow$_ nat  _$\Rightarrow$_ nat set _›_
_and_

— variables the reward function depends upon

actions ::  _‹_ ’a set _›_  _and_

— finite, nonempty set of actions

transitions ::  _‹_ ’a  _$\Rightarrow$_ nat  _$\Rightarrow$_ ’x state
_$\Rightarrow$_ ’x pmf _›_  _and_

— pmf transitions a i x y: probability that next state at dim i will be y, for
action x.

transitions_scope ::  _‹_ ’a  _$\Rightarrow$_ nat  _$\Rightarrow$_ nat set _›_
_and_

— scopes of the transition functions

l :: real  _and_

— discount factor < 1

h_dim ::  _‹_ nat _›_  _and_

— number of basis functions

h ::  _‹_ nat  _$\Rightarrow$_ ’x state  _$\Rightarrow$_ real _›_  _and_

— basis functions, estimate value function

h_scope ::  _‹_ nat  _$\Rightarrow$_ nat set _›_

— scopes of the basis functions

_locale_ Factored_MDP = Factored_MDP_Consts

 _where_

rewards = rewards

 _for_

rewards ::  _‹_ ’a::countable  _$\Rightarrow$_ nat  _$\Rightarrow$_
’x::countable state  _$\Rightarrow$_ real _›_ +

 _assumes_

transitions_closed:

 _‹_ _$\bigwedge$_ a i x.

a  _$\in$_ actions  _$\Longrightarrow$_

i < dims  _$\Longrightarrow$_

x  _$\in$_ partial_states_on (transitions_scope a i)  _$\Longrightarrow$_

set_pmf (transitions a i x)  _$\subseteq$_ doms i _›_  _and_ — transitions go
to domains

transitions_scope:

 _‹_ _$\bigwedge$_ a i.

a  _$\in$_ actions  _$\Longrightarrow$_

i < dims  _$\Longrightarrow$_

has_scope_on (transitions a i) partial_states (transitions_scope a i)_›_
_and_

— the transition functions have the specified scopes

transitions_scope_dims:

 _‹_ _$\bigwedge$_ a i.

a  _$\in$_ actions  _$\Longrightarrow$_

i < dims  _$\Longrightarrow$_ transitions_scope a i  _$\subseteq$_
{..<dims}_›_  _and_

— the scopes are valid

actions_fin:  _‹_ finite actions _›_  _and_

actions_ne:  _‹_ actions  _$\not=$_ {}_›_  _and_

— finite, nonempty action sets

doms_fin:  _‹_ _$\bigwedge$_ i. i < dims  _$\Longrightarrow$_ finite (doms
i)_›_  _and_

doms_ne:  _‹_ _$\bigwedge$_ i. i < dims  _$\Longrightarrow$_ doms i  _$\not=$_
{}_›_  _and_

— finite, nonempty domains

dims_pos:  _‹_ dims > 0 _›_  _and_

— there exists a domain

reward_scope:

_‹_ _$\bigwedge$_ a i.

a  _$\in$_ actions  _$\Longrightarrow$_

i < reward_dim a  _$\Longrightarrow$_

has_scope_on (rewards a i) partial_states (reward_scope a i)_›_  _and_

reward_scope_dims:

_‹_ _$\bigwedge$_ a i.

a  _$\in$_ actions  _$\Longrightarrow$_

i < reward_dim a  _$\Longrightarrow$_ reward_scope a i  _$\subseteq$_
{..<dims}_›_  _and_

— the reward functions are proper scoped functions

h_scope:

_‹_ _$\bigwedge$_ i.

i < h_dim  _$\Longrightarrow$_

has_scope_on (h i) partial_states (h_scope i)_›_  _and_

h_scope_dims:

_‹_ _$\bigwedge$_ i. i < h_dim  _$\Longrightarrow$_ h_scope i  _$\subseteq$_
{..<dims}_›_  _and_

— the basis functions are proper scoped functions

disc_lt_one:  _‹_ l < 1 _›_  _and_

disc_nonneg:  _‹_ l  _$\geq$_ 0 _›_ — valid discount factor

Listing 2: Factored MDPs with a Default Action

* _locale_ Factored_MDP_Default_Consts = Factored_MDP_Consts

 _where_

transitions = transitions

 _for_

transitions ::  _‹_ ’a::{countable, linorder}  _$\Rightarrow$_ nat
_$\Rightarrow$_ ’x::{countable, linorder} state  _$\Rightarrow$_ ’x pmf _›_ +

 _fixes_

d ::  _‹_ ’a _›_  _and_

effects ::  _‹_ ’a  _$\Rightarrow$_ nat set _›_

_locale_ Factored_MDP_Default =

Factored_MDP_Default_Consts

 _where_

transitions = transitions +

Factored_MDP

 _where_

transitions = transitions

 _for_

transitions ::  _‹_ ’a::{countable, linorder}  _$\Rightarrow$_ nat
_$\Rightarrow$_ ’x::{countable, linorder} state  _$\Rightarrow$_ ’x pmf _›_ +

 _assumes_

default_act:  _‹_ d  _$\in$_ actions _›_  _and_

effects:  _‹_ _$\bigwedge$_ i a.

i < dims  _$\Longrightarrow$_

a  _$\in$_ actions  _$\Longrightarrow$_

(_$\exists\,$_ x  _$\in$_ partial_states. transitions a i x  _$\not=$_
transitions d i x)  _$\Longrightarrow$_

i  _$\in$_ effects a _›_  _and_

effects_default:

 _‹_ effects d = {}_›_  _and_

rewards_default_dim:  _‹_ _$\bigwedge$_ a. a  _$\in$_ actions
_$\Longrightarrow$_ reward_dim a  _$\geq$_ reward_dim d _›_  _and_

rewards_eq:  _‹_ _$\bigwedge$_ a i. a  _$\in$_ actions  _$\Longrightarrow$_ i
< reward_dim d  _$\Longrightarrow$_ rewards a i = rewards d i _›_  _and_

reward_scope_eq:  _‹_ _$\bigwedge$_ a i. a  _$\in$_ actions
_$\Longrightarrow$_ i < reward_dim d  _$\Longrightarrow$_ reward_scope a i =
reward_scope d i _›_

— The first _reward_dim d_ reward functions coincide for all actions.

Listing 3: (Partial) States

* _definition_  _‹_ partial_states = {x.  _$\forall\,$_(i, j)  _$\in$_ Map.graph x. i  _$\in$_ vars  _$\wedge$_ j  _$\in$_ doms i}_›_

_definition_  _‹_ states = partial_states  _$\cap$_ {x. dom x = vars}_›_

_definition_  _‹_ partial_states_on X = partial_states  _$\cap$_ {x. X
_$\subseteq$_ dom x}_›_

_definition_  _‹_ partial_states_on_eq X = partial_states  _$\cap$_ {x. X =
dom x}_›_

_definition_  _‹_ consistent x x’  _$\longleftrightarrow$_ x |‘ dom x’ = x’
_›_

Listing 4: Scoped Functions

* _definition_  _‹_ has_scope_on f D R =

(_$\forall\,$_ d  _$\in$_ D.  _$\forall\,$_ d’  _$\in$_ D. restrict d R =
restrict d’ R  _$\longrightarrow$_ f d = f d’)_›_

_definition_  _‹_ instantiate f t x = f (x ++ t)_›_

Listing 5: Basic Definitions on Factored MDPs

* _definition_ p _a_ ::  _‹_ ’a  _$\Rightarrow$_ ’x state  _$\Rightarrow$_ ’x state  _$\Rightarrow$_ real _›_  _where_

"p _a_ a x x’ = (_$\prod\,$_ i < dims. pmf (transitions a i x) (the (x’ i)))"

_definition_  _‹_ reward a x = (_$\sum\,$_ i < reward_dim a. rewards a i x)_›_

Listing 6: Linear Value Functions

* _definition_  _‹_ _$\nu$__w_ w x = (if x  _$\in$_ partial_states then (_$\sum\,$_ i < h_dim. w i * h i x) else 0)_›_

_definition_ g ::  _‹_ nat  _$\Rightarrow$_ ’a  _$\Rightarrow$_ ’x state
_$\Rightarrow$_ real _›_  _where_

 _‹_ g i a x = (_$\sum\,$_ x’  _$\in$_ states. p _a_ a x x’ * h i x’)_›_

_definition_  _$\Gamma$__a_ ::  _‹_ ’a  _$\Rightarrow$_ nat set
_$\Rightarrow$_ nat set _›_  _where_

 _‹_ _$\Gamma$__a_ a I = (_$\bigcup\,$_ i  _$\in$_ I. transitions_scope a
i)_›_

_lemma_ scope_g:

 _assumes_ a:  _‹_ a  _$\in$_ actions _›_

 _assumes_ i:  _‹_ i < h_dim _›_

 _shows_  _‹_ has_scope_on (g i a) states (_$\Gamma$__a_ a (h_scope i))_›_

_lemma_ Q_g:

 _assumes_  _‹_ a  _$\in$_ actions _›_  _‹_ x  _$\in$_ partial_states _›_

 _shows_  _‹_ Q w x a = reward a x + l * (_$\sum\,$_ i < h_dim. w i * g i a
x)_›_

####  A.2.2 Approximate Policy Iteration

Listing 7: Approximate Policy Iteration

* _locale_ API_Interface_Consts =

Factored_MDP_Default_Consts +

 _fixes_

 _$\varepsilon$_ :: real  _and_ — desired precision

t_max :: nat  _and_ — max number of iterations

update_dec_list ::  _‹_ weights  _$\Rightarrow$_ (’a state  _$\times$_ ’b)
list _›_  _and_

update_weights ::  _‹_ nat  _$\Rightarrow$_ (’a state  _$\times$_ ’b) list
_$\Rightarrow$_ weights _›_

_begin_

_definition_  _‹_ is_greedy_dec_list_pol w p  _$\longleftrightarrow$_

is_dec_list_pol p  _$\wedge$_

(is_greedy_pol_w w (dec_list_to_pol p)) — turned into a policy, it is optimal
wrt the weights _›_

_definition_  _‹_ opt_weights_pol p w = is_arg_min (proj_err_w p) (_$\lambda$_
_. True) w _›_

_definition_  _‹_ update_weights_spec  _$\longleftrightarrow$_ (_$\forall\,$_
p i. is_dec_list_pol p  _$\longrightarrow$_ opt_weights_pol (dec_list_to_pol
p) (update_weights i p))_›_

_definition_  _‹_ dec_list_pol_spec  _$\longleftrightarrow$_ (_$\forall\,$_ w.
is_greedy_dec_list_pol w (update_dec_list w))_›_

_definition_  _‹_ proj_err_w p w = (_$\bigsqcup\,$_ x  _$\in$_ states. dist
(_$\nu$__w_ w x) (Q w x (p x)))_›_

_function_ api_aux  _where_

_‹_ api_aux pw t = (

let

(p, w) = pw;

w’ = update_weights t p;

p’ = update_dec_list w’;

err = proj_err_w (dec_list_to_pol p’) w’;

t’ = t + 1;

w_eq = (_$\forall\,$_ i < h_dim. w’ i = w i);

err_le = (err  _$\leq$_  _$\varepsilon$_);

timeout = (t’  _$\geq$_ t_max) in

(if w_eq  _$\vee$_ err_le  _$\vee$_ timeout then (w’, p’, err, t, w_eq,
err_le, timeout)

else api_aux (p’, w’) t’))_›_

_definition_  _‹_ api =

(let

w0 = (_$\lambda$_ _. 0);

p0 = update_dec_list w0

in api_aux (p0, w0) 0)_›_

_end_

_locale_ API_Interface =

API_Interface_Consts +

Factored_MDP_Default +

 _assumes_

update_weights_spec:  _‹_ update_weights_spec _›_  _and_

dec_list_pol_spec:  _‹_ dec_list_pol_spec _›_

_begin_

_theorem_ api_correct:

 _assumes_  _‹_ api = (w, p, err, t, True, err_le, timeout)_›_

 _shows_  _‹_(1 - l) * pol_err (dec_list_to_pol p)  _$\leq$_ 2 * l * err _›_

_end_

####  A.2.3 Decision List Policy

Listing 8: Decision List Policy Definitions

* _definition_ U _a_ ::  _‹_ ’a  _$\Rightarrow$_ nat set _›_  _where_

 _‹_ U _a_ a = (_$\bigcup\,$_ i  _$\in$_ {reward_dim d..<reward_dim a}.
reward_scope a i)_›_

_definition_  _‹_ I _a_ a = {i. i < h_dim  _$\wedge$_ effects a  _$\cap$_
h_scope i  _$\not=$_ {}}_›_

_definition_  _$\Gamma$__a_ ::  _‹_ ’a  _$\Rightarrow$_ nat set
_$\Rightarrow$_ nat set _›_  _where_

 _‹_ _$\Gamma$__a_ a I = (_$\bigcup\,$_ i  _$\in$_ I. transitions_scope a
i)_›_

_definition_  _$\Gamma$__a_ ’ ::  _‹_ ’a  _$\Rightarrow$_ nat set
_$\Rightarrow$_ nat set _›_  _where_

 _‹_ _$\Gamma$__a_ ’ a X =  _$\Gamma$__a_ a X  _$\cup$_  _$\Gamma$__a_ d X _›_

_definition_ T _a_ ::  _‹_ ’a  _$\Rightarrow$_ nat set _›_  _where_

 _‹_ T _a_ a = U _a_ a  _$\cup$_ (_$\bigcup\,$_ i  _$\in$_ I _a_ a.
_$\Gamma$__a_ ’ a (h_scope i))_›_

_definition_ bonus ::  _‹_ weights  _$\Rightarrow$_ ’a  _$\Rightarrow$_ ’x
state  _$\Rightarrow$_ real _›_  _where_

 _‹_ bonus w a x = Q w x a - Q w x d _›_

Computationally efficient version of _bonus_.

_definition_  _‹_ bonus’ w a x = R _a_ a x + l * (_$\sum\,$_ i  _$\in$_ I _a_
a. w i * (g’ i a x - g’ i d x))_›_

_definition_ dec_list_act ::  _‹_ weights  _$\Rightarrow$_ ’a  _$\Rightarrow$_
(’x state  _$\times$_ ’a  _$\times$_ real) list _›_  _where_

 _‹_ dec_list_act w a = (

let

ts = assignment_list (T _a_ a);

ts’ = filter (_$\lambda$_ t. bonus’ w a t > 0) ts;

ts’’ = map (_$\lambda$_ t. (t, a, bonus’ w a t)) ts’

in

ts’’

)_›_

_abbreviation_  _‹_ _$\pi$_ _unsorted w  _$\coloneqq$_ (_$\lambda$_ _. None,
d, 0) # concat (map (dec_list_act w) actions_nondef)_›_

_abbreviation_  _‹_ _$\pi$_ _sorted w  _$\coloneqq$_ sort_dec_list (_$\pi$_
_unsorted w)_›_

_definition_ dec_list_pol ::  _‹_ weights  _$\Rightarrow$_ (’x state
_$\times$_ ’a) list _›_  _where_

 _‹_ dec_list_pol w = map (_$\lambda$_(t, a, _). (t, a)) (_$\pi$_ _sorted
w)_›_

_definition_ dec_list_pol_sel ::  _‹_ weights  _$\Rightarrow$_ ’x state
_$\Rightarrow$_ ’a _›_  _where_

 _‹_ dec_list_pol_sel w = dec_list_to_pol (dec_list_pol w)_›_

Listing 9: Properties of the Decision List Policy

* The main property of _dec_list_pol_ is: Given a state, the first element of the policy it is consistent with is the greedy action in that state.

_lemma_ dec_list_pol_sel_greedy:

 _assumes_  _‹_ x  _$\in$_ states _›_

 _shows_  _‹_ is_greedy_act w x (dec_list_pol_sel w x)_›_

_lemma_ dec_list_pol_is_pol:  _‹_ is_dec_list_pol (dec_list_pol w)_›_

####  A.2.4 Bellman Error

Listing 10: Bellman Error

* _locale_ Bellman_Err_Consts =

Factored_MDP_Default_Consts

 _where_ rewards = rewards  _for_

rewards :: "’a::{countable, linorder}  _$\Rightarrow$_ nat  _$\Rightarrow$_
’x::{countable, linorder} state  _$\Rightarrow$_ real" +

_fixes_

error_branch ::  _‹_ ’x state  _$\Rightarrow$_ ’a  _$\Rightarrow$_ (’x state)
list  _$\Rightarrow$_ ereal _›_  _and_

dec_pol ::  _‹_(’x state  _$\times$_ ’a) list _›_  _and_

w ::  _‹_ weights _›_

_begin_

_definition_  _‹_ greedy_spec =

(_$\forall\,$_ x  _$\in$_ states. is_greedy_act w x (dec_list_to_pol dec_pol
x))_›_

_definition_  _‹_ dec_pol_spec  _$\longleftrightarrow$_ is_dec_list_pol
dec_pol  _$\wedge$_ greedy_spec _›_

_definition_  _‹_ error_branch_spec = (_$\forall\,$_ t a ts.

a  _$\in$_ actions  _$\longrightarrow$_ t  _$\in$_ partial_states
_$\longrightarrow$_ set ts  _$\subseteq$_ partial_states  _$\longrightarrow$_

error_branch t a ts =

(_$\bigsqcup\,$_ x  _$\in$_ {x. x  _$\in$_ states  _$\wedge$_ consistent x t
_$\wedge$_ list_all (_$\lambda$_ t’.  _$\neg$_ consistent x t’) ts}.

ereal (dist (Q w x a) (_$\nu$__w_ w x))))_›_

_definition_  _‹_ update_err_iter =

(_$\lambda$_(t, a) (ts, err). (t#ts, sup (error_branch t a ts) err))_›_

_definition_  _‹_ err_list xs = snd (fold update_err_iter xs ([], 0))_›_

_definition_  _‹_ factored_bellman_err = real_of_ereal (err_list dec_pol)_›_

_end_

Listing 11: Factored Bellmann Error Correctness

* _lemma_ bellman_err_w_eq_Q:

 _shows_  _‹_ bellman_err_w w =

(_$\bigsqcup\,$_ x _$\in$_ states. dist (Q w x (dec_list_to_pol dec_pol x))
(_$\nu$__w_ w x))_›_

Listing 12: Error of a Single Branch

* _locale_ Bellman_Err_Branch_Consts = Factored_MDP_Default_Consts

 _where_

rewards = rewards

 _for_

rewards ::  _‹_ ’a::{countable, linorder}  _$\Rightarrow$_ nat

 _$\Rightarrow$_ (’x::{countable, linorder} state)  _$\Rightarrow$_ real _›_ +

 _fixes_

variable_elim ::  _‹_((’x state  _$\Rightarrow$_ ereal)  _$\times$_ (nat set))
list  _$\Rightarrow$_ ereal _›_

— Variable elimination algorithm, maximizes sum of input functions

_begin_

Specification of the variable elimination algorithm.

_definition_  _‹_ variable_elim_correct fs  _$\longleftrightarrow$_

variable_elim fs = (MAX x  _$\in$_ states. sum_list (map (_$\lambda$_(f, s). f
x) fs))_›_

Valid input functions: don’t take (positive) infinite values + scopes are
valid.

_definition_  _‹_ scope_inv b  _$\longleftrightarrow$_ (_$\forall\,$_(b,
scope_b)  _$\in$_ set b.

(_$\forall\,$_ x. b x  _$\not=$_  _$\infty$_)  _$\wedge$_ scope_b
_$\subseteq$_ vars  _$\wedge$_ has_scope_on b partial_states scope_b)_›_

_definition_  _‹_ scope_inv’ b  _$\longleftrightarrow$_ (_$\forall\,$_(b,
scope_b)  _$\in$_ set b.

(_$\forall\,$_ x. b x  _$\not=$_ -_$\infty$_  _$\wedge$_ b x  _$\not=$_
_$\infty$_)  _$\wedge$_ scope_b  _$\subseteq$_ vars  _$\wedge$_ has_scope_on b
partial_states scope_b)_›_

_definition_  _‹_ variable_elim_spec  _$\longleftrightarrow$_ (_$\forall\,$_
fs. scope_inv fs  _$\longrightarrow$_ variable_elim_correct fs)_›_

Functions to maximize that depend on the weights.

_definition_  _‹_ hg_scope t a i = h_scope i  _$\cup$_  _$\Gamma$__a_ a
(h_scope i) - dom t _›_

_definition_  _‹_ hg_inst w t a i = instantiate (_$\lambda$_ x. w i * (h i x -
l * g’ i a x)) t _›_

_definition_  _‹_ C w t a = (map (_$\lambda$_ i. (hg_inst w t a i, hg_scope t
a i)) [0..<h_dim])_›_

_definition_  _‹_ neg_C w t a = map neg_scoped (C w t a)_›_

Functions to maximize that independent of the weights.

_definition_  _‹_ r_act_dim a = reward_dim a - reward_dim d _›_

_definition_  _‹_ r_scope t a i = reward_scope a i - dom t _›_

_definition_  _‹_ r_inst t a i = instantiate (rewards a i) t _›_

_definition_  _‹_ b t a = (map (_$\lambda$_ i. (r_inst t a i, r_scope t a i))
[0..<reward_dim a])_›_

_definition_  _‹_ neg_b t a = map neg_scoped (b t a)_›_

Indicator functions that make us ignore states that choose earlier branches.

_definition_  _‹_ _$\mathcal{I}$_ t x = (if consistent x t then -_$\infty$_
::ereal else 0)_›_

_definition_  _‹_ _$\mathcal{I}$_ ’ ts t’ = map (_$\lambda$_ t. (instantiate
(_$\mathcal{I}$_ t) t’, dom t - dom t’)) ts _›_

_definition_  _‹_ scope__$\mathcal{I}$_ ’ ts t’ =  _$\bigcup\,$_(set (map snd
(_$\mathcal{I}$_ ’ ts t’)))_›_

_definition_  _‹_ b’ t a ts = (b t a) @  _$\mathcal{I}$_ ’ ts t _›_

_definition_  _‹_ neg_b’ t a ts = (neg_b t a) @  _$\mathcal{I}$_ ’ ts t _›_

The maximum of the negative + positive error is the absolute error.

_definition_  _‹_ _$\varepsilon$_ _pos w t a ts = variable_elim (C w t a @
neg_b’ t a ts)_›_

_definition_  _‹_ _$\varepsilon$_ _neg w t a ts = variable_elim (neg_C w t a @
b’ t a ts)_›_

_definition_  _‹_ _$\varepsilon$_ _max w t a ts = max (_$\varepsilon$_ _pos w
t a ts) (_$\varepsilon$_ _neg w t a ts)_›_

_end_

Listing 13: Branch Error Correctness

* _lemma_  _$\varepsilon$_ _max_correct’’:

 _assumes_  _‹_ t  _$\in$_ partial_states _›_  _‹_ a  _$\in$_ actions _›_  _‹_
set ts  _$\subseteq$_ partial_states _›_

 _shows_  _‹_ _$\varepsilon$_ _max w t a ts =

(_$\bigsqcup\,$_ x  _$\in$_ {x. x  _$\in$_ states  _$\wedge$_ consistent x t
_$\wedge$_ (_$\forall\,$_ t’  _$\in$_ set ts.  _$\neg$_ consistent x t’)}.

ereal (dist (Q w x a) (_$\nu$__w_ w x)))_›_

Listing 14: Variable Elimination Definitions

* _locale_ Variable_Elimination = Variable_Elimination_Consts  _$\mathcal{F}$_ dims  _$\mathcal{O}$_ doms

 _for_

 _$\mathcal{F}$_ ::  _‹_(((nat  _$\rightharpoonup$_ ’a)  _$\Rightarrow$_
ereal)  _$\times$_ nat set) list _›_  _and_ — Functions to be maximized

dims ::  _‹_ nat _›_  _and_ — Number of variables

 _$\mathcal{O}$_ ::  _‹_ nat  _$\Rightarrow$_ nat _›_  _and_ — Elimination
order

doms ::  _‹_ nat  _$\Rightarrow$_ ’a set _›_ — Domains of the dimensions +

 _assumes_

 _$\mathcal{O}$_ _bij:  _‹_ bij_betw  _$\mathcal{O}$_ {..<dims} {..<dims}_›_
_and_

 _$\mathcal{F}$_ _scopes:  _‹_ _$\forall\,$_(f, s)  _$\in$_ set
_$\mathcal{F}$_. s  _$\subseteq$_ {..<dims}  _$\wedge$_

has_scope_on f {x.  _$\forall\,$_(i, j)  _$\in$_ Map.graph x. i  _$\in$_
{..<dims}  _$\wedge$_ j  _$\in$_ doms i} s _›_  _and_

domains_ne:  _‹_ _$\bigwedge$_ i. i < dims  _$\Longrightarrow$_ doms i
_$\not=$_ {}_›_  _and_

domains_fin:  _‹_ _$\bigwedge$_ i. i < dims  _$\Longrightarrow$_ finite (doms
i)_›_  _and_

 _$\mathcal{F}$_ _not_inf:  _‹_ _$\bigwedge$_ f x. f  _$\in$_ set
_$\mathcal{F}$_  _$\Longrightarrow$_ (fst f) x  _$\not=$_  _$\infty$__›_

_begin_

The body of the algorithm, in a single iteration, we first select the variable
to eliminate as __$\mathcal{O}$_ i_, then partition the work set of functions
_F_ based on whether __$\mathcal{O}$_ i_ is in their scope. Next we build a
new function that maximizes all functions that contain _@{term _‹_
_$\mathcal{O}$_ i _›_}_ over this variable, and add it to the working set.

_definition_  _‹_ elim_step iF = (

let

(i, F) = iF;

l =  _$\mathcal{O}$_ i;

E = filter (_$\lambda$_ f. l  _$\in$_ scope f) F;

e = (_$\lambda$_ x. MAX x _l_  _$\in$_ doms l.  _$\sum\,$_ f  _$\leftarrow$_
E. (fn f (x(l  _$\mapsto$_ x _l_))));

scope_e = (_$\bigcup\,$_ f _$\in$_ set E. scope f) - {l}

in

(i + 1, (e, scope_e) # filter (_$\lambda$_ f. l  _$\notin$_ scope f) F)

)_›_

The body of the algorithms needs to be iterated _dims_ times.

_definition_  _‹_ elim_aux = (elim_step ^^ dims) (0,  _$\mathcal{F}$_)_›_

Finally, we are left with only functions with empty scopes, and we return
their sum.

_definition_  _‹_ elim_max = (_$\sum\,$_(f, f_scope)  _$\leftarrow$_ snd
elim_aux. f Map.empty)_›_

_expl_max_ is the inefficient variant of the algorithm that enumerates all
states.

_definition_  _‹_ expl_max__$\mathcal{F}$_  _$\mathcal{F}$_ ’ = (MAX x
_$\in$_ full_vecs.  _$\sum\,$_ f  _$\leftarrow$_  _$\mathcal{F}$_ ’. fst f
x)_›_

_definition_  _‹_ expl_max = expl_max__$\mathcal{F}$_  _$\mathcal{F}$__›_

The invariant _invar_max_ can be used to show the algorithm _elim_max_
correct.

_definition_  _‹_ invar_max  _$\mathcal{F}$_ ’ i  _$\longleftrightarrow$_

(_$\forall\,$_ f  _$\in$_ set  _$\mathcal{F}$_ ’.  _$\forall\,$_ x. fst f x
_$\not=$_  _$\infty$_)  _$\wedge$_

(_$\forall\,$_ f  _$\in$_ set  _$\mathcal{F}$_ ’. has_scope_on (fst f)
partial_vecs (snd f))  _$\wedge$_

(_$\forall\,$_ f  _$\in$_ set  _$\mathcal{F}$_ ’. snd f  _$\subseteq$_
{..<dims} -  _$\mathcal{O}$_ ‘ {..<i})  _$\wedge$_

(MAX x _$\in$_ vecs_on ({..<dims} -  _$\mathcal{O}$_ ‘ {..<i}).  _$\sum\,$_ f
_$\leftarrow$__$\mathcal{F}$_ ’. fst f x) = expl_max__$\mathcal{F}$_
_$\mathcal{F}$__›_

_lemma_ elim_max_correct:  _‹_ elim_max = expl_max _›_

_end_

###  A.3 Policy Evaluation

Listing 15: Policy Evaluation / Weight Update

* _locale_ ValueDet_API_Consts = …

_fixes_

min_sol ::  _‹_ real  _$\times$_ (nat  _$\Rightarrow$_ real)_›_

dec_pol ::  _‹_(’x state  _$\times$_ ’a) list _›_

_begin_

_definition_  _‹_ dec_pol_spec  _$\longleftrightarrow$_ fst ‘ set dec_pol
_$\subseteq$_ partial_states  _$\wedge$_ snd ‘ set dec_pol  _$\subseteq$_
actions  _$\wedge$_ distinct dec_pol  _$\wedge$_

(_$\forall\,$_ x  _$\in$_ states.  _$\exists\,$_ y  _$\in$_ fst ‘ set dec_pol.
consistent x y)_›_

_definition_  _‹_ update_weights_iter =

(_$\lambda$_(t, a) (ts, cs). (t#ts, union_constr (gen_constr t a ts) cs))_›_

_definition_  _‹_ constrs_list xs = snd (fold update_weights_iter xs ([],
empty_constr))_›_

_definition_  _‹_ constrs = constrs_list dec_pol _›_

_definition_  _‹_ update_weights = snd min_sol _›_

_end_

Listing 16: Policy Evaluation Correctness

* _lemma_ proj_err_upd_eq_pol:  _‹_ proj_err_w (dec_list_to_pol dec_pol) update_weights =

proj_err_pol (dec_list_to_pol dec_pol)_›_

Listing 17: Linear Program for a Branch

* _locale_ Process_Branch_Consts = …

_fixes_

t ::  _‹_ ’x state _›_  _and_

a ::  _‹_ ’a _›_  _and_

ts ::  _‹_ ’x state list _›_  _and_

factored_lp ::  _‹_

bool  _$\Rightarrow$_

((’x state  _$\Rightarrow$_ real)  _$\times$_ (nat set)) list  _$\Rightarrow$_
— C

((’x state  _$\Rightarrow$_ ereal)  _$\times$_ (nat set)) list
_$\Rightarrow$_ — b

’c _›_

_begin_

Here we define the specification of the Factored LP algorithm:

_definition_  _‹_ factored_lp_constrs pos C b  _$\longleftrightarrow$_

constr_set (factored_lp pos C b) = {(_$\varphi$_ , w).  _$\forall\,$_ x
_$\in$_ states.

ereal  _$\varphi$_  _$\geq$_ ereal (_$\sum\,$_ i < length C. w i * (fst (C !
i) x))

\+ (_$\sum\,$_ i < length b. fst (b ! i) x)}_›_

_definition_  _‹_ factored_lp_inv pos C b  _$\longleftrightarrow$_

inv_constr (factored_lp pos C b)_›_

_definition_  _‹_ factored_lp_privs pos C b  _$\longleftrightarrow$_

privates (factored_lp pos C b)  _$\subseteq$_ {((t, a), pos)}_›_

_definition_  _‹_ inv_b b  _$\longleftrightarrow$_

(_$\forall\,$_(b, scope_b)  _$\in$_ set b. has_scope_on b partial_states
scope_b  _$\wedge$_ scope_b  _$\subseteq$_ vars  _$\wedge$_ (_$\forall\,$_ x.
b x  _$\not=$_  _$\infty$_))_›_

_definition_  _‹_ inv_C b  _$\longleftrightarrow$_

(_$\forall\,$_(b, scope_b)  _$\in$_ set b. has_scope_on b partial_states
scope_b  _$\wedge$_ scope_b  _$\subseteq$_ vars)_›_

_definition_  _‹_ factored_lp_spec  _$\longleftrightarrow$_

(_$\forall\,$_ pos C b. inv_C C  _$\longrightarrow$_ inv_b b
_$\longrightarrow$_

factored_lp_constrs pos C b  _$\wedge$_ factored_lp_inv pos C b  _$\wedge$_

factored_lp_privs pos C b)_›_

_definition_  _‹_ neg_scoped = (_$\lambda$_(f, s). (-f, s))_›_

_definition_  _‹_ hg_scope i = h_scope i  _$\cup$_  _$\Gamma$__a_ a (h_scope
i) - dom t _›_

_definition_  _‹_ hg_inst i = instantiate (_$\lambda$_ x. h i x - l * g’ i a
x) t _›_

_definition_  _‹_ C = (map (_$\lambda$_ i. (hg_inst i, hg_scope i))
[0..<h_dim])_›_

_definition_  _‹_ neg_C = map neg_scoped C _›_

— the function _r -_$\mathcal{I}$__ becomes _b_

_definition_  _‹_ r_act_dim = reward_dim a - reward_dim d _›_

_definition_  _‹_ r_scope i = reward_scope a i - dom t _›_

_definition_  _‹_ r_inst i = instantiate (rewards a i) t _›_

_definition_  _‹_ b = (map (_$\lambda$_ i. (r_inst i, r_scope i))
[0..<reward_dim a])_›_

_definition_  _‹_ neg_b = map neg_scoped b _›_

_definition_  _‹_ _$\mathcal{I}$_ t’ x = (if consistent x t’ then -_$\infty$_
::ereal else 0)_›_

_definition_  _‹_ _$\mathcal{I}$_ ’ = map (_$\lambda$_ t’. (instantiate
(_$\mathcal{I}$_ t’) t, dom t’ - dom t)) ts _›_

_definition_  _‹_ scope__$\mathcal{I}$_ ’ =  _$\bigcup\,$_(set (map snd
(_$\mathcal{I}$_ ’)))_›_

_definition_  _‹_ b’ = b @  _$\mathcal{I}$_ ’ _›_

_definition_  _‹_ neg_b’ = neg_b @  _$\mathcal{I}$_ ’ _›_

_definition_  _‹_ _$\Omega$_ _pos = factored_lp True C neg_b’ _›_

_definition_  _‹_ _$\Omega$_ _neg = factored_lp False neg_C b’ _›_

_definition_  _‹_ _$\Omega$_ = union_constr  _$\Omega$_ _pos  _$\Omega$_ _neg
_›_

_end_

Listing 18: Linear Program for a Branch – Correctness

* _lemma_  _$\Omega$_ _set_correct’:

 _shows_  _‹_ constr_set  _$\Omega$_ = {(_$\varphi$_ , w).

 _$\forall\,$_ x  _$\in$_ {x. x  _$\in$_ states  _$\wedge$_ consistent x t
_$\wedge$_ (_$\forall\,$_ t’  _$\in$_ set ts.  _$\neg$_ consistent x t’)}.

Listing 19: Factored Linear Program

* _locale_ factored_lp_consts =  _fixes_

C :: "(((nat  _$\rightharpoonup$_ ’a)  _$\Rightarrow$_ real)  _$\times$_ nat
set) list"  _and_

B :: "(((nat  _$\rightharpoonup$_ ’a)  _$\Rightarrow$_ ereal)  _$\times$_ nat
set) list"  _and_

doms :: "nat  _$\Rightarrow$_ ’a set"  _and_

— Each dimension has a domain

dims :: nat

— for each _i < num_c_, _c i v_ assigns to each vector _v _$\in$_ Dom num_c_ a
real value  _and_

prefix ::  _‹_ ’x _›_  _and_

— prefixes are added to private LP variables to make them distinct between
branches

order :: "nat  _$\Rightarrow$_ nat"

_begin_

Enumerate all functions in c + all states with matching scopes.

_definition_  _‹_ vars_c = {var_f’ (f_c i) z | z i. z  _$\in$_ vecs_on (scope_c i)  _$\wedge$_ i < num_c}_›_

_definition_  _‹_ constr_c z i = Eq (_$\lambda$_ v.

if v = var_f’ (f_c i) z then -1

else if v = var_w i then c_r i z

else 0) 0 _›_

_definition_ "constrs_c =

{constr_c z i | z i. z  _$\in$_ vecs_on (scope_c i)  _$\wedge$_ i < num_c}"

_definition_ "vars_b = { var_f’ (f_b j) z | z j.

z  _$\in$_ vecs_on (scope_b j)  _$\wedge$_ j < num_b}"

Ensures _f (b, j, z) = b_j z_

_definition_ "constr_b z j =

Eq (_$\lambda$_ v. if v = var_f’ (f_b j) z then 1 else 0) (real_of_ereal (b_r
j z))"

_definition_ "constrs_b =

{constr_b z j | z j. z  _$\in$_ vecs_on (scope_b j)  _$\wedge$_ j < num_b  _$\wedge$_ b_r j z  _$\not=$_ -_$\infty$_}"

_definition_ "vars_w = {var_w i| i. i < num_c}"

Initial set of function to maximize constructed from B and C.

_definition_ "_$\mathcal{F}$_ _init = {(f_c i) | i. i < num_c}  _$\cup$_ {(f_b i) | i. i < num_b}"

_definition_ "scopes_init f = (case f of f_b i  _$\Rightarrow$_ scope_b i | f_c i  _$\Rightarrow$_ scope_c i | _  _$\Rightarrow$_ {})"

_definition_ "constrs_init = constrs_c  _$\cup$_ constrs_b"

_definition_ "vars_init = vars_c  _$\cup$_ vars_b  _$\cup$_ vars_w"

Create constraints that ensure _f_e_ is at least the sum of the functions in
E: sum f in E, _f(z(l := xl) <= f(e, l)_.

_definition_ "constr_max E l scopes z xl =

Le (_$\lambda$_ v.

case v of

var_f p f z’  _$\Rightarrow$_

if p = prefix  _$\wedge$_ f = f_e l  _$\wedge$_ z’ = z then -1

else if p = prefix  _$\wedge$_ f  _$\in$_ E  _$\wedge$_ z’ = (z(l  _$\mapsto$_
xl)) |‘ (scopes f) then 1

else 0

| _  _$\Rightarrow$_ 0) 0"

_definition_ "constrs_max E l scopes scope_e =

{constr_max E l scopes z xl | xl z. z  _$\in$_ vecs_on scope_e  _$\wedge$_ xl  _$\in$_ doms l}"

Single iteration step: 1\. select variable to eliminate, 2\. partition
functions based on that, 3\. create new function (variable) with constraints
to ensure maximization.

_definition_ "elim_step  _$\Omega$_  _$\mathcal{F}$_ scopes i = (let

l = order i;

E = {e | e. e  _$\in$_  _$\mathcal{F}$_  _$\wedge$_ l  _$\in$_ scopes e};

scope_e = (_$\bigcup\,$_ e  _$\in$_ E. scopes e) - {l};

 _$\Omega$_ ’ =  _$\Omega$_  _$\cup$_ constrs_max E l scopes scope_e in

(_$\Omega$_ ’,  _$\mathcal{F}$_ \- E  _$\cup$_ {f_e l}, scopes(f_e l :=
scope_e), i+1)

)"

_definition_ "elim_vars_aux =

((_$\lambda$_(_$\Omega$_ ,  _$\mathcal{F}$_ , scopes, i). elim_step
_$\Omega$_  _$\mathcal{F}$_ scopes i) ^^ dims) (constrs_init,  _$\mathcal{F}$_
_init, scopes_init, 0)"

_definition_ "gen_constrs arg = (

case arg of (_$\Omega$_ ,  _$\mathcal{F}$_ , scopes, i)  _$\Rightarrow$_
(_$\Omega$_  _$\cup$_ {Le (_$\lambda$_ v.

case v of

var_f p f z  _$\Rightarrow$_ if p = prefix  _$\wedge$_ f  _$\in$_
_$\mathcal{F}$_  _$\wedge$_ z = Map.empty then 1 else 0

| var_phi  _$\Rightarrow$_ -1

| _  _$\Rightarrow$_ 0) 0}))"

_definition_ "elim_vars = gen_constrs elim_vars_aux"

_end_

Listing 20: Factored LP Correctness

* _lemma_ constr_set_factored_eq’:

 _‹_{(x var_phi,  _$\lambda$_ i. x (var_w i)) |x. x  _$\in$_ ({v |v. Ball
elim_vars (sat_constr v)})} =

{(_$\varphi$_ , w).  _$\forall\,$_ x _$\in$_ full_vecs.

ereal (_$\sum\,$_ i<length C. w i * fst (C ! i) x) +

(_$\sum\,$_ i<length B. fst (B ! i) x)  _$\leq$_ ereal  _$\varphi$_}_›_