1 Introduction
Agent-based models (ABMs) have been widely adopted as an intuitive way to model systems that consist of a heterogeneous collection of autonomous, interacting agents. ABMs are sometimes used to show that an unexpected collective behaviour can emerge from a known agent behaviour. For example Schelling (1971) showed that agents with only a very slight racial bias can quickly form a highly racially-segregated population. This use of ABMs requires no assimilation of real-world data but generates results that cannot be used to assert anything about the real world without also providing a convincing reason to believe that the model is an abstraction of the real-world. In this paper, we’ll deal with an alternative use case in which we start with the following information:
A set of real-world observations, Ω, that may contain aggregated or macro-scale quantities and may be subject to noise during the measurement process.
An agent-based model expressed as a set of agent behaviours. These will generally be stochastic, in order to account for our uncertainty in agent behaviour, and will be a function of some set of unknown parameters, θ, in order to account for our uncertainty over which behavioural model best represents the real-world target entities. Note that the parameters can, and generally should, contain values that control the size of the model stochasticity in order to account for "unknown unknowns". The agents may also interact with a world outside of the model. These "boundary conditions" necessarily consist of agents being injected into the model at certain times and/or external influences on modelled agents’ behaviour1. The start state of the model can also be treated as a boundary condition consisting of the injection of agents into the model at time t = 0. Although the boundary conditions are semantically distinct from the parameters, their mathematical treatment will be identical so we assume our prior beliefs about boundary conditions is contained in θ. We’ll show later how to express an ABM as a "forecast" distribution, P(τ|θ), which is the probability that a set of agents would exhibit behaviours τ on a model execution given parameters/boundary conditions θ. We’ll call τ a model trajectory.
A set of probabilistic prior beliefs, P(θ), about the model parameters and boundary conditions (after accounting for any micro-calibration or other relevant direct observations we may have).
Our aim will be to calculate the expectation of some (possibly vector-valued) unobserved quantity, Λ, over the posterior distribution P(τ, θ|Ω)
Inserting Bayes’ rule gives us
where P(Ω|τ) is the observation likelihood, which we assume is easy to calculate given a trajectory that covers the observed period, and P(Ω) is just a normalising constant which is not usually explicitly calculated when sampling. Note that any information about P(τ, θ|Ω) can be expressed without loss of generality in this way.
There are various approaches in the literature to solving this equation for ABMs. One strategy is to first transform the problem into a simpler form, such as a set of partial differential equations (Lloyd et al., 2016) or a graphical model (Liao & Barooah, 2010), for which there exists a well known method of solving the problem. Tang (2019) expressed the problem in terms of creation and annihilation operators then used symbolic computation to evaluate the integral in (2). However, the most popular strategy in the literature is to sample from the posterior, P(τ, θ|Ω), and use the samples to approximate the integral using Monte Carlo integration. This is the technique we employ here.
1.1 Importance sampling from the prior
The simplest Monte Carlo approach would be to approximate 𝔼P(τ,θ|Ω)(Λ(τ, θ)) using importance sampling as follows:
1. Taking a sample from the prior parameters/boundary conditions θi ~ P(θ).
2. Executing the model forward with the given θi to generate a trajectory τi ~ P(τ|θi).
3. Calculating a weight wi = P(Ω|τi) to create a weighted sample 〈τi, θi, wi〉.
4. Repeating N times from step 1 to get a set of weighted samples {〈τ1, θ1, w1〉 . . . 〈τN, θN, wN〉}.
5. Approximating
.
However, Chatterjee & Diaconis (2018) showed that the information contained in the weighted samples decreases exponentially with the KL-divergence from the prior to the posterior, DKL (P(τ, θ|Ω) || P(τ|θ)P(θ)). This divergence can be thought of as a measure of the information contained in the observations Ω. For reasons we’ll see later, in the case of ABMs, the observations usually contain enough information to make the set of weighted samples less informative than a single unweighted sample, even when N is astronomically high. Therefore, this algorithm will not work unless the observations contain very little information.
1.2 Particle filtering
In the literature, the most popular way of solving (2) for ABMs is to use some form of particle filtering, also known as "sequential Monte Carlo". This splits the time period of interest into smaller, contiguous time windows with (not necessarily equidistant) end times 〈t1 . . . tN〉 where tm < tm+1. The trajectory and observations are similarly split based on which window they occur in, τ = 〈τ1 . . . τN〉 and Ω = 〈Ω1 . . . ΩN〉. The posterior distribution of the whole trajectory can then be split into the product of contributions from each window given the distribution over the previous window, leading to a recursion relation
where τ1:t = 〈τ1 . . . τt〉 and P(τt+1|τ1:t, θ) is the probability that an execution of the ABM that starts with trajectory τ1:t would go on to produce τt+1 given the parameters/boundary conditions θ. The final posterior P(τ1:N, θ|Ω1:N) can now be built up in steps, starting with the first window and moving to the last, using the recursion as we go and approximating each P(τ1:t, θ|Ω1:t) as a set of samples, or "particles". This makes the problem easier in two respects: firstly the dimension of the distributions we need to deal with at each recursion step are reduced, since if we consider each particle separately, we’re only concerned with the trajectory in one window at a time; secondly, we split the information in the observations into smaller chunks, so at each recursion the information in Ωt is more likely to be small enough to make it practical to perform importance sampling from the prior.
If the unobserved quantity of interest, Λ, depends only on the states of the agents, σN, at time tN and/or the parameters/boundary conditions, θ, rather than on the full trajectory, τ1:N, then we can reduce the dimensionality even further by performing the recursion on model states rather than trajectories (where a model state tells us how many agents there are in each agent state). First, we arrange the windows so that each observation lies at the end of a window, then marginalise equation (3) over τ1:t+1 for a fixed σt+1 to get
where we’ve assumed that the observation likelihoods depend only on the model state and made explicit in P(σt+1|σt, θ) that the probability of the state σt+1 depends only on σt and θ2. Note that the definition of "model state" can always be re-defined in such a way as to make these assumptions true. In fact, we can see that equation (3) can be thought of as a special case of equation (4) where the model state, σt, is replaced by the trajectory, τ1:t, in which case the integral over σt is non-zero for only a single trajectory, since P(
|τ1:t, θ) = δτ1:t()P(
|τ1:t, θ), where δ here is the multivariate delta distribution.
So, we can consider a recursion on generalised trajectories, Xt,
so that equations (3) and (4) would correspond to Xt = 〈τ1:t, θ〉 and Xt = 〈σt, θ〉 respectively.
A particle filter solves this recursion by taking a set of samples drawn from P(Xt|Ω1:t), and generating a set of samples drawn from P(Xt+1|Ω1:t+1). The simplest algorithm to do this is known as sequential importance resampling, which consists of the following steps:
1. Generate a set of samples from the prior P(X0). Set t = 0.
2. For each sample , generate a sample from the forecast ~ P(Xt+1|) by executing the ABM from the end state of to give a set of forecast samples .
3. Calculate a weight for each forecast sample
to give a new set of unweighted samples . There are a few ways of performing this resampling (Douc & Cappé, 2005), the simplest being to draw an integer i ∈ [1, R] with probability wi, returning the sample and repeating R times.
5. If t < N, increment t and repeat from step 2.
However, this algorithm suffers from the problem of "sample impoverishment" (Li et al., 2014) or "sample deprivation" which describes the situation when the resampling step leaves many particles in the same state, so the number of distinct samples in our sample set is much smaller than the total number of samples. This is because, as we’ve seen, after each importance sampling step the information in the observations causes the effective sample size to decrease and therefore the resampling step is likely to generate repeated samples. This loss of effective sample size can accumulate and in practice as the number of windows increases, the number of particles required to maintain a desired accuracy often increases exponentially. In the worst case, we encounter an observation that is impossible for all particles and the algorithm fails completely. Malleson et al. (2020) showed that sample impoverishment is an issue when applied to data assimilation with an ABM of crowd movement. Khan et al. (2003) encountered the same problem in a model of ant movement.
To get a better understanding of impoverishment, suppose we use the above algorithm to generate a set of samples {〈
, θ1〉 . . . 〈
, θS〉} of the full trajectory P(τ1:N, θ|Ω1:N), from which we generate marginalisations over each window Wt = {〈, θ1〉 . . . 〈, θS〉}. On the one hand, the repeated impoverishment at each step causes the number of distinct values in WN −L to decrease exponentially as the lag, L, increases (i.e. as we look further back in time), meaning we lose information about the true distribution of P(τN −L, θ|Ω1:N). However, on the other hand, the model dynamics and observations make P(τN, θ|Ω1:N) increasingly independent of P(τN −L, θ|Ω1:N) as L increases. So if all we care about is the final window, WN , (which is often the case) impoverishment only becomes a problem when the first effect (the rate of loss of information) dominates the second (the rate of increase in independence from the past).
Unfortunately every window’s trajectory depends on the parameters, θ, but these are sampled in the first window and never subsequently change, so after assimilating only a few windows, all samples are likely to have collapsed to the same value of θ, giving us little information about the true distribution of the marginalised posterior P(θ|Ω1:N). This is fine if the parameters are known at t = 0, but if there is prior uncertainty in θ and the observations don’t reduce that uncertainty faster than the collapse in information contained in the samples, then impoverishment will be a problem (Andrieu et al., 2004; Liu & West, 2001).
A simple way of dealing with sample impoverishment is to add noise to the particles at each step. This is known as roughening (Gordon et al., 1993; Li et al., 2014) and is equivalent to replacing the delta function in equation (6) with some other function with wider support so we end up with a weighted kernel density estimation. Kieu et al. (2020) used this technique in the context of an ABM of public transport. However, Liu & West (2001) showed that this leads to over-dispersion unless we "shrink" the samples towards their mean value. More generally, the consequence of roughening is that a finite number of samples from a roughened particle filter is no longer a draw from the true posterior, when averaged over all possible draws of the same size, so we need to convince ourselves that the error introduced by the roughening process is small enough for our requirements; not always an easy task for an ABM where distributions can be highly discontinuous.
Wang & Hu (2015) dealt with sample impoverishment in the context of ABM by using the samples at time t to approximate the marginal probabilities of the state of each individual agent at the start of the next window, P(aj), then drawing new samples from the product of the marginals P(a1 . . . an) = P(aj), i.e. sampling the state of the jth agent from the marginal P(aj). This improves sample diversity at the price of losing all correlations between agents.
In the context of ABMs, sample impoverishment can manifest in its most extreme form by generating particles that have zero weight and so contribute no information to the resampling step. This happens because a draw from the forecast of an ABM P(Xt+1|Xt) often has an extremely high probability of having a likelihood P(Ωt|Xt+1) of zero. For example, suppose we’re tracking mobile phones by IMEI and receive intermittent signals from the phones identifying which cell they’re in. Suppose for simplicity that if a phone is in a known cell at time t then it could be in any of eight cells at time t + 1, each with equal prior probability. Suppose at time t + 1 we get signals from 10 phones. Even if we know the correct model state at time t, the probability that a sample from the forecast at time t + 1 has a non-zero likelihood is 8−10, about one in a billion. So sequential importance resampling will fail in a single step in this case.
1.3 Kalman filtering
One potential solution to the problem of impoverishment is to assume that some or all of the distributions are Gaussian. This is the basis of the data assimilation algorithms that have had a great deal of success in geophysical models (Carrassi et al., 2018; Kalnay, 2003; Lewis et al., 2006; Talagrand, 1997). If the forecast P(Xt+1|Ω1:t) and likelihood function P(Ωt+1|Xt+1) are multivariate Gaussians, the recursion in (5) can be solved analytically and evaluated even when the observations Ωt are highly informative. This idea leads to the Ensemble Kalman filter (Evensen, 2003) and the Unscented Kalman filter (Wan et al., 2001), both of which have been applied to data assimilation in ABM (Clay et al., 2020; Ward et al., 2016).
However, the Gaussian assumption is not usually a good one for ABMs. For example, suppose there is a sports hall containing N people and we take the model state to be the 2N dimensional vector consisting of the (x, y) coordinates of each person. Even if the initial state is distributed according to a 2N dimensional Gaussian distribution, the either/or decisions of the people can quickly transform this into a multimodal distribution. Suppose two exits are opened on opposite sides of the hall and each person moves towards their nearest exit. Each person makes an either/or decision which exit to move towards based on their position, so we can think of the prior Gaussian as being partitioned into 2N regions corresponding to all possible sets of decisions made by the people. In the forecast distribution, the probability mass in each region will move towards a different model state (one of the 2N states where all agents are at one or other of the exits), creating a highly non-Gaussian distribution with 2N modes. In addition, the likelihood function is often highly non-Gaussian in this state space. For example, suppose there is a sensor at one of the exits that trips if any agents are close. Suppose for example that there are just two agents so the model state space is (x1, y1, x2, y2) and the origin of the agents’ coordinate system is at the sensor. If the sensor trips, the likelihood function is zero everywhere except close to the perpendicular planes (0, 0, x2, y2) and (x1, y1, 0, 0). This clearly isn’t well approximated by any Gaussian. This is a consequence of the fact that we don’t know which agent caused the sensor to trip, a problem known as the "data association problem" (Lueck et al., 2019). This problem becomes superexponentially worse if we have multiple sensors. If there are N agents and M ≤ N sensors the likelihood function is the union of hyperplanes.
An alternative representation of the model state space is the occupation number representation where each agent state becomes a dimension in the state space and the coordinate values give the number of agents in each state. In this representation, Gaussians may be appropriate in some models if there are a large number of agents in every state. However, as the occupation numbers get smaller, an increasing proportion of a Gaussian approximation lies on states that contain negative occupation numbers, which is problematic as a negative occupation number has no natural interpretation. In very high dimensions, the vast majority of the probability mass will lie in a region with at least one negative occupation number. Likelihood functions in this space can also be highly non-Gaussian. For example, the sensor in (a discretised version of) the sports hall example above would have a likelihood function that is a step function on the dimension that corresponds to the state close to the door.
1.4 Metropolis-Hastings
In the following sections we present a solution to these problems by using the Metropolis-Hastings algorithm to construct a Markov process whose stationary state is P(τ1:N, θ|Ω1:N). To our knowledge this is the first time this has been done in a way that is applicable to a wide range of ABMs and observations. For small N we can sample directly from this Markov process to get samples from the posterior. For larger N, or if we wish to perform online data assimilation where the observations form a continuous, effectively endless stream, we propose using the algorithm as part of a data assimilation cycle or MCMC Particle Filter algorithm (Finke et al., 2020; Septier et al., 2009). There are many specific ways of implementing this; unfortunately the most appropriate algorithm to use depends on the nature of the ABM dynamics, the likelihood function and the observation operator, Λ, so there’s no silver bullet. We discuss this further in section 5 and present some algorithms that are likely to be appropriate in the context of ABMs.
Since the nature of an ABM’s parameters is highly model-specific, in the following we focus on the treatment of the trajectory, although in a way that is fully integrated with a model-specific treatment of the parameters for which we give more general design advice. In the case of state estimation (where the parameters are already known) the parameter-specific part can be left out completely.
2 Methods
2.1 Defining the problem
2.1.1 Definition of an ABM. Suppose we have a timestepping ABM that consists of agents with a finite number of possible internal states and a finite number of mutually exclusive ways of acting on their world. Given this, we formally define an ABM as:
An ordered list of agent states 𝒮 = 〈σ0...σS−1〉.
An ordered list of agent actions 𝒜 = 〈α0...αA−1〉.
An agent timestep, π : ℤ × ℤS × ℤ → ℝ, which defines the probability that an agent will act in a particular way such that π(𝜓, Ψ, a) gives the probability that an agent in state σ𝜓, surrounded by agents, Ψ, will perform action αa (where Ψ is an S dimensional vector whose ith element is the number of agents in state σi at the start of the timestep).
An action function, F : ℤ × ℤ → ℤS, which defines the effect of an action on the world such that F(𝜓, a) is an S dimensional vector whose ith element gives the number of agents in state σi that result from an agent in state σ𝜓 performing act αa (including the final state of the acting agent).
As a simple illustration, we define the "cat and mouse" ABM as follows:
Agent states An agent can be either a cat or a mouse and can be on one of two grid-squares, left or right, so
Agent actions In any given timestep, an agent can either move to the other grid-square or stay still, so
Agent timestep In a timestep, a cat will move or stay still with probability 0.5, a mouse will move if there are any cats on the same gridsquare or stay still otherwise. This behaviour can be expressed as the agent timestep:
Action function This gives the result of an agent’s actions. For example, F(1, 0) is the result of an agent in state 𝜓 = 1 (right cat) performing the action a = 0 (move), which is one cat in state 𝜓 = 0 (left cat). So F(1, 0) = {1, 0, 0, 0}. This model is simple enough that we can write down F explicitly for all 𝜓 and a:
2.1.2 A note on tensor notation. In the rest of this paper we’ll make use of multidimensional arrays. These are just arrays of numbers, like vectors or matrices, but arranged in any number of dimensions. Each array has a shape which tells us how many dimensions it has and the number of elements in each dimension. For notational convenience, we’ll distinguish between subscript dimensions and superscript dimensions3 and will refer to the set of all arrays of a given shape with the symbol ℝ adorned with the size of each subscript and superscript dimension. So, for example, describes the set of all three-dimensional arrays that have one superscript dimension of size N and two subscript dimensions of sizes S and A.
An element of an array is referred to by specifying the sub- and superscript coordinates of the element. For example if T ∈ , then refers to the element of T that has coordinate t in the superscript dimension and coordinates (𝜓) in the subscript dimensions. By convention, coordinates begin at 0.
We’ll also borrow from tensor notation by using the Einstein summation convention, meaning that if the same index symbol is repeated in sub- and superscript positions in a term, then a summation over that dimension is implied. So, for example if X ∈ ℝ8 and Y ∈ ℝ8 then
When the same symbol is repeated in the same position with no implied summation, this implies universal quantification. For example if Y ∈ R8 then
is equivalent to
We refer to "slices" of an array by using the * symbol in an index position. For example if T ∈ then refers to the one-dimensional array in ℝA comprised of the elements 〈...〉.
The symbols 0 and 1 represent arrays whose elements are all 0 or 1 respectively. Their shape should be unambiguous from the context.
2.1.3 Trajectories. Let a model timestep be an array E ∈ ℝSA whose elements E𝜓a give the number of agents in state σ𝜓 that perform act αa in this timestep. For example, the timestep shown in Figure 1 for the cat and mouse example would be
where one agent in state σ1 (right cat) performs action α0 (move) and one agent in state σ2 (left mouse) performs action α1 (stay still).
Let a model trajectory, T, be an array in that represents N timesteps of a model with S agent states and A actions, so that denotes the (𝜓, a)th element of the tth timestep matrix.
Figure 1. A simple cat and mouse model.
An array must satisfy a number of constraints in order to be a valid trajectory of an ABM. Since the elements of a trajectory are counts of agents, they must be non-negative integers. We’ll denote the set of all non-negative integer arrays
A trajectory must also be continuous, by which we mean that any agent that is the result of an action at timestep t − 1 must be the cause of an action at timestep t, and any agents injected by the boundary conditions at timestep t must also act in timestep t4. We call this the "continuity constraint" and define the set of continuous arrays, with respect to an action function F:
where F ∈ is the action array = F (𝜓, a) and is the number of boundary condition injections into state ϕ at time t. Note that this assumes that a timestep is performed by updating all agents in parallel, i.e. agents must act with no information about the actions of other agents in the same timestep. This is different from serial update, where agents are updated in a particular order within a timestep and have access to information about the actions of agents that come before them in the ordering. Note also that this assumes that agents in the same state are indistinguishable, i.e. if we swap two agents in the same state, the trajectory is unchanged.
So, we define the set of trajectories, (F), as a set of arrays that satisfy (7) and (8).
2.1.4 The posterior of an ABM. Our goal is to generate samples from the Bayesian posterior over the trajectories, boundary conditions and model parameters, given a set of observations, Ω
The model forecast, P (T |θ), can be decomposed into the product of conditionals for each timestep and agent state
where ∈ ℕS is the model state at the start of timestep t and P(|θ) is defined by the boundary conditions5. The term P( |, 𝜓, θ) is the probability that agents starting in state 𝜓 will collectively perform the vector of actions . The probability that a single agent performs action a is given by the agent timestep πθ (𝜓, Ψt, a) (where we’ve added the subscript θ to emphasise its possible dependence on the parameters). Since agent actions are simultaneous, the collective behaviour is given by the multinomial distribution
If the trajectory is continuous then so the model forecast is given by
The likelihood function, P (Ω|T), can also usually be decomposed. Without loss of generality, we take Ω to consist of some number of observations that are independent of each other given the trajectory, so that Ω is a set of pairs (ω, v) that consist of a stochastic observation operator ω and an observed value v (which may be a vector). We write P(ω(T) = v) to denote the probability of observation operator ω making observation v on trajectory T. So
and the posterior can be written as
2.2 Approximating the support of the posterior
As we saw in the introduction, sampling from equation (12) is difficult because it is often very hard to find proposal trajectories that have non-zero posterior probability. Our strategy to solve this problem is to first derive an expression for a volume of trajectory/parameter space that bounds the support of the posterior, supp(P(T, θ|Ω)) (i.e. the set of 〈T, θ〉pairs that have non-zero posterior probability). If we choose a bounding volume that it is easy to sample from and is a relatively tight approximation of supp(P(T|Ω)) then there’s a good chance that a sample drawn from this volume will have non-zero probability. The bounding volume doesn’t need to be exact, it just needs to be tight enough to give us a good chance of finding a non-zero probability sample.
From equation (12)
i.e. in order for 〈T, θ〉 to have non-zero posterior probability, θ must have non-zero prior probability and T must be a trajectory of the ABM that has a possible start state, all observations must have non-zero likelihood and each non-zero element of T must denote an agent action with non-zero probability.
2.2.1 Convex mixed-integer polyhedra. We now approximate equation (13) with a convex, mixed-integer polyhedron, 𝒫 (see, for example, Conforti et al. (2010)), which we define to be a set of vectors, some of whose elements are integer and some real-valued, that satisfy a set of linear inequalities:
for some matrix C and some vectors L and U.
Any linear inequality on T and θ can be expressed in this form by letting X consist of the elements of T and θ flattened into one dimension so that
and W and G have inverses
and
such that
and
where δij and δij are 1 if i = j and 0 otherwise. In this way we can switch unambiguously between X and 〈T, θ〉.
We now consider how to express each term in equation (13) as a mixed integer polyhedron. From equation (9) we can see immediately that the set of all trajectories, , is already defined as a set of linear constraints on T. The support of the prior P(θ) can usually be easily bounded by a hyper-rectangle giving the range of each variable. The supports of the start state, , the observations, P(ω(T) = v), and the agent actions, , can often be easily expressed as linear constraints. However, if this is not the case, each of the probability distributions can be expressed as a computer program. In practice, these computer programs will be simple and it will be possible to use a technique known as "abstract interpretation" (Cousot & Cousot, 1977) using the domain of convex polyhedra (Becchi & Zaffanella, 2018; Cousot & Halbwachs, 1978; Fukuda, 2020) to efficiently calculate a convex polyhedron that bounds the support of the program. Software to perform abstract interpretation using convex polyhedra already exists (Bagnara et al., 2008; Gurfinkel & Navas, 2021; Henry et al., 2012; Jeannet & Miné, 2009) and the technique has been used in applications such as compiler optimization (Sjödin et al., 2009) and verification of safety-critical systems (Halbwachs et al., 1997). The programs may contain calls to a random number generator Random() that returns a random floating-point number 0 ≤ r < 1. This can be represented in the polyhedral domain by introducing an real valued auxiliary variable, r, that satisfies 0 ≤ r < 1 for each call to Random(). These auxiliary variables can be removed immediately using the abstract interpretation software to project the polyhedron back into the space of X. This relies on converting the polyhedron to its set of vertices (Motzkin et al., 1953) and projecting the vertices, so the number of vertices must be of a tractable size. Alternatively, the auxiliary variables can be added to X and marginalised out when we take the samples.
If the number of agents is very much smaller than the number of agent states (which is often the case with agent-based models) then we may be willing to make the assumption that in any timestep there is at most one agent performing a given action from a given start state (i.e. ∀𝜓, a, t : ∈ {0, 1}). Under this assumption, which we’ll call the "Fermionic assumption", the "Fermionic trajectories", all lie on the corners of the unit hypercube. So the intersection of with any set is a convex polyhedron and supp(P(T|Ω)) can always be exactly represented as a mixed-integer polyhedron (finding that polyhedron isn’t always tractable, though, in which case abstract interpretation software can be used to find a larger polyhedron that contains it).
So, if we let 𝒫(f) be a mixed integer polyhedron that contains supp(f) then from equation (13) we can approximate the support of the posterior as
The intersection of two mixed-integer polyhedra is easy to express as another mixed integer polyhedron by just concatenating the constraints
so the only difficulty in calculating 𝒫(P(X|Ω)) is the union in the final term of equation (14). We transform this into an intersection by introducing an auxiliary variable z and using the following identity: if the hyper-rectangle 0 ≤ X ≤ H bounds L ≤ CX ≤ U for some vector H ∈ ℤN × ℝM, then for any given i
where B are upper bounds on the values of CX in the hyper-rectangle, defined as
where |C| here denotes an element-wise absolute value of a matrix and B are lower bounds on the values of CX defined as
To see why this identity holds, note first that the constraints make z into an indicator variable that is 0 if Xi = 0 or 1 otherwise. When z = 1 the first constraint on the right hand side of (16) is equal to CX ≤ U and the second is equal to L ≤ CX so their intersection is L ≤ CX ≤ U as required, whereas when z = 0 the constraints are equal to B ≤ CX ≤ B. But B and B are lower and upper bounds on the values of CX in the hyper-rectangle so this is satisfied for all values 0 ≤ X ≤ H as required.
Two things are worth noting here. Firstly if we make the Fermionic assumption then z = Xi = so the auxiliary indicator variable becomes unnecessary. Secondly, we must bound by a hyper-rectangle with its lower corner at the origin. In practice, this is not a problem as we can just mak a change of variable, imposing upper and lower bounds on each variable such that the probability of being outside those bounds is negligible, then shifting if the lower bound isn’t on zero.
So, the support of the posterior can be approximated by a mixed-integer polyhedron, 𝒫(P(X|Ω)), by using equations (14), (15) and (16).
As a simple illustration, consider a two-timestep trajectory of the cat and mouse model described in section 2.1.1. Suppose we flip a fair coin to decide whether each agent state is occupied or empty at t = 0. Suppose also that we observe a cat in the left grid-square at time t = 1. Our aim is to construct a polyhedron, 𝒫(P(X|Ω)), that bounds the support of the posterior.
Working through (14) term by term, the term constrains to non-negative, continuous trajectories as defined in equation (9), which is already in linear form. The second term is the support of the prior over the parameters which we ignore here since there are no parameters. The next term is the support of the start state, since we’re flipping a coin for each agent state at t = 0, each state occupation must be at most 1, which can be expressed as
for each 𝜓. The fourth term is the support of the observation. Since we observe a cat in the left grid-square at time t = 1 we add the constraint
The final term is the constraint due to agent interactions. The impossible interactions are a mouse staying put when there is a cat on the same gridsquare or moving when there are no cats, so we constrain against these
for all t. If, for simplicity, we make the Fermionic assumption by adding the constraints
then using the identity in (16) to take the union of each constraint in (17) with {T : = 0, 0 ≤ T ≤ 1} gives the four constraints
for each timestep t = 0 and t = 1.
Taken together, these constraints define a polyhedron that bounds the set of (Fermionic) trajectories for the cat and mouse ABM with the given start state distribution and observation.
2.3 Sampling from the posterior
Having shown how to calculate 𝒫(P(X|Ω)), we now show how to use this to sample from P(X|Ω).
Various methods of sampling from a mixed-integer polyhedron exist in the literature. Baumert et al. (2009) presents a very general algorithm for sampling from any weighted subset of integer points within a bounding hyper-rectangle. The algorithm relies on sampling points on a random walk. However, in our case the number of valid trajectories could be a vanishingly small proportion of the number of points in the smallest bounding hyper-rectangle, so a random walk would be unlikely to come across any valid trajectories. Universal hashing (Meel et al., 2016) provides a promising way of sampling uniformly from the integer points in a convex polyhedron, but we have found that this technique doesn’t scale to the number of dimensions needed for our application.
The Metropolis-Hastings algorithm is another well known and widely used sampling algorithm and this is the approach we take here. To do this we need to define
a set of Markov states, ℳ, along with a mapping, E : ℳ → ℤN × ℝM, which maps Markov states to sample values, X,
a stochastic proposal function f : ℳ → ℳ from which we can generate transitions between Markov states,
a probability distribution P : ℳ → ℝ which gives the probability of each Markov state (this need not be normalised as the Metropolis-Hastings algorithm only ever needs probability ratios).
In order to be of use in practice, the proposal function, f, must have the following properties:
For any two Markov states there should exist a sequence of transitions which forms a path between those states and has a non-zero probability of being proposed and accepted.
Given a current Markov state, there should be a computationally efficient procedure to generate a proposal and calculate its acceptance probability.
2.3.1 The set of Markov states. Given the polyhedron 𝒫(P(X|Ω)), we split the constraints into two sets: equalities (i.e. those whose lower and upper bound have the same value) and inequalities (i.e. those whose lower and upper bounds differ) so that
Our aim will be to use the equality constraints to reduce the dimension of the problem and so improve the efficiency of the algorithm. Note that D has at least (N − 1)S rows since the continuity constraints (8) are equality constraints.
Suppose there are Ne equality constraints and we partition the elements of X into ‘basic’ and ‘non-basic’ elements, XB and XN respectively, so that there are exactly Ne basic elements.
If we let Q be a permutation matrix that separates the basic from non-basic elements so that
then we can write
so
Now, since there are Ne basic variables and Ne equality constraints, B is square so if we choose the basic variables to ensure that B has an inverse then
so letting
gives
In section 2.3.5 we’ll show how to choose XB to ensure that B has an inverse and the integer elements of XB remain integer as long as the integer elements of XN are integer. Given this, then the original polyhedron is equivalent to the reduced dimension polyhedron
where and are the number of integer and real-valued elements remaining in XN.
Since Q is just a permutation matrix, XN and XB can be bound within hyper-rectangles by transforming the 𝒫-bounding hyper-rectangle, H, from section 2.2.1
So, XN is bound by 0 ≤ XN ≤ HN. Given this, we define the set of Markov states to be
and equation (21) maps each Markov state, XN ∈ ℳ, to a unique sample, X. Although this set includes values outside of the polyhedron (22), we’ll show a way to deal with that later.
2.3.2 The proposal function. In order to propose a new Markov state given a current state, XN, first choose an element of XN and a direction in which to perturb that element. This defines a rounded perturbed state, , which is equal to XN but with one of its elements rounded to the nearest integer and perturbed by ±1. Only rounded perturbed states that remain inside the hyper-rectangle HN are considered. Real valued variables should be scaled so that a perturbation of ±1 usually changes P(X|Ω) by a small amount. The boundaries, HN, should be chosen to lie on an integer value plus 0.5 for real-valued variables. So, we can define the set of all valid perturbations
where ei is the unit vector with a 1 in the ith element, υ is the set of indices of XN and we use the notation ⌊.⌉i to denote the rounding of the ith element of a vector to the nearest integer.
Assume for now that we have an approximation of the target distribution (XN) ≈ P(XN|Ω). Given a current state, XN, the probability of choosing perturbation 〈i, δ〉 ∈ 𝒟(XN) is given by
where
Given a perturbation 〈i, δ〉, the final proposed state is given by
where 𝜖 is a fractional perturbation drawn with uniform probability from [−0.5 : 0.5] if the ith element is real-valued or 𝜖 = 0 if it is integer-valued.
2.3.3 Approximating the posterior. The proposal function requires an approximation of the target distribution, P(X|Ω). For this, we use a function of the form
where Zi are row vectors and are univariate factors of any form. We call this a "linearly factorized distribution", where the linearity lies in the ZiX terms while the can be non-linear.
The first thing to note is that we can express 𝒫(P(X|Ω)) as a linearly factorized distribution in the following way: first define the infeasibility function ι(l, x, u) to be equal to 0 if l < x < u or equal to the distance to the nearest bound otherwise
Now, given a set of constraints L ≤ CX ≤ U, if we define a factor for each constraint
then
is a linearly factorized distribution that is 1 inside the polyhedron and that decays exponentially outside at a rate set by τ. So as τ → 0 equation (27) tends to an indicator function for the set of points in 𝒫(P(X|Ω)).
A linearly factorized distribution in X can easily be transformed into a linearly factorized distribution in XN by using the linear transform in equation (21). Transformation to XN ensures that the resulting distribution satisfies the equality constraints D(V + MXN) = E, while the factors in equation (27) can be used to ensure that 0 ≤ XB ≤ HB and that the constraints in (22) are satisfied.
Turning our attention now to the distribution inside the polyhedron, if we look back at equation (12) we can see that the posterior of an ABM is already highly factorized and many of the factors are already linearly factorised.
However, when considering the efficiency of the final algorithm, there is a trade off between the closeness of to P(XN|Ω) and the number of factors involved in its calculation. In our experiments, we have found we can attain very good acceptance probabilities with relatively simple renderings of The start state, likelihood function and prior over θ are often already uncorrelated functions of the elements of θ and the state occupation numbers, so are already in the correct form. This leaves only the multinomial term in equation (12) which we approximate with the linearly factorized
where
This integral can be calculated by drawing samples from the prior P(T, θ) and using Monte Carlo integration. This approximation marginalises over θ, but if πθ has a particularly strong dependence on θ then it may be worth capturing some of that dependence by, for example, replacing the term in (28) with + Rθ for some row vector R. The can then be calculated by minimising the expectation of the difference between the log probabilities of the normalised approximation and the true multinomial.
So, the final approximation, (XN), is given by the product of the factors of the bounding convex polyhedron (27), the multinomial approximation (28) and the approximations of the observations, the start state (boundary conditions) and P(θ).
Since some of the Markov states are infeasible, calculating for these states may require the evaluation of a factor at a value that is outside its domain (for example, equation (28) is undefined if is negative). In this case we choose the value at the nearest point that is within the domain.
2.3.4 The probability of a Markov state. If a Markov state satisfies all the constraints in (22) then its probability is defined to be equal to the true posterior given by (12)
where
If the constraints are not all satisfied, then the probability is defined to be the product of the factors in (12) and the infeasibility factors in (27). If any factor in (12) becomes zero as a consequence of the infeasibility, then it is replaced by its linearly factorized approximation (which must be non-zero otherwise the state wouldn’t have been proposed).
This means that, as long as τ is finite, the Markov chain will pass through some infeasible states. However, if we just ignore these and only take samples from feasible states then we end up with samples from the true posterior. Allowing the Markov chain to pass through infeasible states has the advantage of ensuring that there exists a path between any two feasible states and improves mixing. When used in the context of a particle filter, it also solves the problem of unexpected observations that can’t be generated by any particle, as we can allow the particles to become temporarily infeasible. The price we pay for this is the computational cost of moving through infeasible states that don’t generate useful samples.
It’s worth noting that equation (27) is a Boltzmann distribution which describes a thermodynamic system whose "energy" is equal to the infeasibility of X that is at equilibrium with a heat bath of "temperature" τ. It has been shown that simulations of thermodynamic systems of this type are able to solve large sets of linear integer inequalities like the one we’re dealing with here (Kirkpatrick et al., 1983). Here, we won’t be changing the temperature during the sampling process, but we do need to choose a value for τ. Higher temperatures will increase mixing in the Markov chain, but will also increase the proportion of time spent in infeasible states so we need to find a temperature that is high enough to ensure good mixing of the chain and low enough to ensure that a reasonable proportion of samples are feasible. We have found that these two competing effects roughly cancel each other out when measuring the overall efficiency of the algorithm in effective samples per second, so the value of τ is not critical as long as it’s within a certain range. We have found that a good rule of thumb is to set the temperature so that 50–80% of the samples in a chain are infeasible.
From the Metropolis-Hastings algorithm, the probability of accepting a proposal = ⌊XN⌉i + (δ + ∈)ei is given by
Since (XN) approximates P(XN), should be close to 1. We’ll see in the next section that if we choose the basis carefully then S(XN) will not change very much between transitions so the ratio should also be close to 1 meaning that a good proportion of proposals should be accepted. In our experiments we achieved acceptance rates of around 85%.
2.3.5 Choosing a basis. Our definition of a Markov state depends on a partition of X into basic and non-basic variables such that B in equation (20) has an inverse, B−1, and if XN ∈ ℤN′
× ℝM′ then X ∈ ℤN × ℝM (i.e. if XN has all integers in the correct elements then so does X). In general there exist many partitions that satisfy these requirements so it remains to define a method of choosing one. The choice of basis affects the efficiency of the Markov chain via:
the computational cost of updating the proposal probabilities P(XN, i, δ) and the sum S(XN) after each transition.
the effect of the ratio of sums S(XN)/ on the acceptance probability
the mixing rate of the Markov chain, since the choice of basis affects the available transitions in X space.
From (23) and (24) we can see that in order to calculate P(XN, i, δ) and S(XN) we need to know all 〈i, δ〉 ∈ 𝒟(XN). Since is linearly factorised
Let
be the un-normalised proposal probability. If we accept a state transition 〈k, 〉, (XN, i, δ) needs to be updated to (⌊XN⌉k + ( + 𝜖)ek, i, δ). However, this update only affects factors j such that {j : Zjek ≠ 0}, so a transition 〈k, 〉 only changes (XN, i, δ) if {j : Zjei ≠ 0} intersects with {j : Zjek ≠ 0}. So, if we let Z be the matrix formed from the row vectors Zi and choose a basis such that Z is sparse, (XN, i, δ) can be stored for each value of 〈i, δ〉 ∈ 𝒟(XN), and only values for a few values of i will need to be recalculated on transition. Tang (2022) shows that a binary sum-tree data structure can be used to efficiently update and draw from the un-normalised probability distribution defined by (XN, i, δ) while also updating S(XN).
The effect of on the acceptance probability is also improved if Z is sparse because only a small proportion of the terms in the sum change between S(XN) and so we should expect their ratio to be close to 1.
The Markov chain will also mix better if the columns of Z are sparse. Mihelich et al. (2018) shows that the Kolmogorov-Sinai entropy of the Markov process is a good proxy for its mixing speed. In the context of our algorithm the entropy is higher if the neighbours of a Markov state have similar probability. If Z is sparse then adjacent Markov states are more likely to have similar probabilities, for a fixed mean element magnitude.
So, our aim will be to find a basis that keeps Z sparse, while maintaining the correct integer signature of X. To do this we use the following algorithm: Start with a linearly factorized distribution on X (rather than XN)
and a set of equality constraints that need to be satisfied for all Markov states XN
From these, form a mixed-integer polyhedron
We assume, without loss of generality, that all constraints involving only integer elements of X have integer coefficients. C is usually sparse to start with, so our aim is to maintain that sparsity while using the equality constraints to reduce X to XN . If we let
and
then
Our aim is to mark variables from X as basic one at a time. Once we’re done, the remaining variables will make up the non-basic variables, XN.
Given a set of constraints in the form (29), we mark the jth element of X as basic and the ith equality constraint as “reduced” by choosing an element to be a pivot point and performing Gaussian elimination on the jth column of W so that is the only non-zero element in the column Wj. This can be done by pre-multiplying by G
where G is the matrix
Letting = GW and = GF recovers the canonical form
but now the jth element of X is marked as basic and the ith row of is marked as "reduced". We now repeat the process, choosing only pivot points that are non-basic and unreduced, until no more equality constraints can be reduced.
The use of Gaussian elimination here ensures that B−1 exists. In order to ensure that XB has the correct integer structure, given that XN does, we add the requirement that we can only mark integer elements of X as basic on a pivot point such that the row vector Wi is zero on all non-integer columns and divides for all k (so for rows that have a non-zero element on a real-valued column, we always eliminate on a real-valued column). This requirement ensures that for integer elements, , after reduction = Fi − WiXN. But if XN has the correct integer structure then WiXN is a sum of products of integers, which is an integer. So, if has any integer solution then Fi must also be an integer. So, if XN has the correct integer structure, then is also an integer.
This defines the set of valid pivot points at each elimination step. In order to decide which point to choose, we note that when we pivot on , GW differs from W in at most
elements, where |·|0 is the L0-norm (i.e. the number of non-zero elements in a vector or matrix). So, this gives an upper bound on |WG|0 − |G|0, the change in L0-norm (or reduction in sparsity) of W when we pivot on . It has been found that choosing to be the element that minimises µ (known as the Markowitz criterion) is a computationally efficient way of maintaining the sparsity of a matrix while performing Gaussian elimination (Markowitz, 1957; Maros, 2002; Suhl & Suhl, 1990). So, in order to find a sparse basis, we calculate the set of valid pivot points, choose the point that minimises µ, perform the pivot and repeat until no more valid pivot points remain. Suhl & Suhl (1990) presents a computationally efficient algorithm for doing this.
3 Spatial predator-prey demonstration
We demonstrate the above techniques on a 32x32 grid of "predator" and "prey" agents. The agent acts consist of moving to an adjacent grid-square (i.e. up, down, left or right), giving birth to another agent on the an adjacent grid-square or dying. If any agent moves off the edge of the grid, it reappears at the opposite edge (i.e. the grid is topologically a torus). Predators only give birth when there is at least one prey on an adjacent grid-square (i.e. if there is a source of food close by). Prey may be eaten by predators so their probability of dying is much greater when predators are present on adjacent grid-squares. The probability of each agent’s behaviour is shown in Table 1. These values were chosen to minimise the probability that either species becomes extinct.
Table 1. The probabilities of each behaviour in the predator-prey model.
| Aget type | Adjacent agents | Behaviour | Probability |
|---|
| Prey | No predators | die | 0.100 |
| | give birth | 0.156 |
| | move | 0.744 |
| Prey | Predators > 0 | die | 0.400 |
| | give birth | 0.156 |
| | move | 0.444 |
| Predator | No prey | die | 0.100 |
| | give birth | 0.000 |
| | move | 0.900 |
| Predator | Prey > 0 | die | 0.100 |
| | give birth | 0.300 |
| | move | 0.600 |
In order to generate observations to assimilate, we used the following procedure:
2. Simulate the ABM for 16 timesteps from the start state. If the resulting trajectory isn’t Fermionic, repeat from step 1. The resulting trajectory, Treal, is considered to be the "real-world" trajectory from which we take observations.
3. For each timestep and gridsquare of Treal, take two draws from a Bernoulli distribution with probability 0.05. If the first draw is 1, then observe the number of predators in that gridsquare at that time. If the second draw is 1 then observe the number of prey. The observations were noiseless so the observed count was the true count.
Given the (Fermionic) model, the start state P() and the observations, an approximation of the posterior was generated as described in section 2.3.3 and used to generate a sparse basis as described in section 2.3.5. The resulting basis had 22.4 non-zero elements per million.
Four initial feasible solutions were generated in the following way:
1. Draw a start state from P().
2. Simulate for 16 timesteps to get a (non-Fermionic) prior trajectory.
3. Starting with the prior trajectory, generate proposals as described in section 2.3.2 and accept with probability 1.0 until a feasible state is reached.
Starting with these four initial solutions, four separate Markov chains were generated using the Matropolis-Hastings algorithm as described in section 2.3. 1,000,000 samples were then taken from each chain and discarded in order to burn them in. Finally, 10,000,000 samples were taken from each chain and split into first and last halves to give eight sample sequences of 5,000,000 samples each. For all samples, the temperature, τ, was set to 0.1.
Feasible samples were generated at a rate of one every 42µs on a single core of a 1.6GHz Intel i5-8250U. 85% of the proposals were accepted and 70% of samples were infeasible.
In order to assess the convergence of the chains, a set of summary statistics were calculated for each sample of each sequence. The summary statistics consisted of the number of agents within each shaded square shown in Figure 2 measured at the end of the last timestep. We’ll refer to these as statistic 1 to 4 going from largest to smallest area. This arrangement of regions was chosen in order to capture the convergence at different spatial scales.
Figure 2. Summary statistics consist of the total number of agents in each of the four shaded regions.
For each statistic, let xij to refer to the value of the statistic for the ith sample of the jth sequence and let
be the mean of the jth sequence and
be the mean over all sequences.
Let W be the within-sequence variance
and B be the between-sequence variance
Following Gelman et al. (2013), we calculated an over-approximation of the true variance of the statistic as
and the Gelman-Rubin diagnostic (Gelman & Rubin, 1992) was calculated using
This gives a measure of the uncertainty in the standard deviation of the statistic due to the finite length of the sample sequences. If this number is close to 1, then we have some justification in believing that we have taken enough samples. A value of less than 1.1 is often used as a criterion for convergence.
Also following Gelman et al. (2013), we approximated the autocorrelation using
where
at various time lags for each sample sequence. The autocorrelation of a convergent chain should decline to zero long before the lag reaches the total number of samples.
t can also be used to calculate the effective number of samples (i.e. the number of independent and identically distributed draws that would give the same uncertainty in the mean of the statistic) using
where the sum runs from t = 0 until the first t such that t ≤ 0.
5 Discussion
5.1 Sequential MCMC for ABM
The predator-prey example demonstrated assimilation over a relatively short time period. In practice it’s likely that assimilation would be required over a longer period, or perhaps over a continuous stream of observations. As discussed earlier, the way to deal with this is to split the longer time period into windows as in equation (5) and treat each window separately. We can make use of the MCMC sampler to implement this broad idea in many different ways.
One way would be the following: given N samples from the posterior of a window ~ P(Xt|Ω1:t), approximate the posterior from some family of simpler, analytical distributions, PΘ(Xt), by finding the Θ that minimises the KL-divergence from PΘ(Xt) to P(Xt|Ω1:t). The KL-divergence can itself be approximated from the samples
with
which is minimised when
The approximation, PΘ(Xt), can then be inserted into the recursion equation (5) to generate samples from the next window ~ P(Xt+1|Ω1:t+1), from where we can repeat the cycle indefinitely. The exact form of PΘ is model dependent. The simplest form would be a product of univariate Poisson distributions, one for each of the integer variables, and Gaussians for the real-valued variables. If we wish to capture correlations between variables, a multivariate Poisson distribution could be used (see for example Muñoz-Pichardo et al. (2021), where PΘ is a product of univariate Poisson distributions whose rate parameters are given by Λ = BXR + AXI , where XR are the real-valued variables and XI are the integer valued variables, B and A are matrices, and A is triangular).
An alternative approach is to use the Resample-Move algorithm presented in Gilks & Berzuini (2001). This is a particle filter based on the sequential importance resampling algorithm described earlier. The innovation comes in the form of an additional MCMC move step after each resample step. At each move step each particle, = 〈, θi〉, is "moved" by initialising the Markov chain with and drawing a single new sample ~ P(Xt|Ω1:t) by taking one or more Markov transitions. This solves the problem of particle impoverishment because even if multiple samples collapse to the same state during the resample step, the subsequent move step is likely to move them apart again. The algorithm remains exact because the samples coming out of the mth resample step are already distributed according to the target distribution, P(Xt|Ω1:t), but this is also the stationary distribution of the Markov transition so the moved samples are also distributed according to the target distribution.
The Resample-Move algorithm still uses an importance sampling step with a draw from the prior which, as we have seen, is often a problem for ABMs because drawing from the prior is likely to result in a sample with zero weight. If this is the case, we can replace the draw from the prior with an draw from P(Xt+1|Ω1:t+1, ) using our MCMC algorithm for each particle from the previous timestep. The correct weight of the particle is then proportional to P(Ωt+1|) (Doucet et al., 2009). Unfortunately we don’t immediately know this value. If we’re doing an assimilation step every timestep, or if the ABM is simple enough, it may be possible to calculate this value analytically, however, often this is not possible. Han & Carlin (2001); Newton & Raftery (1994); Štefankovič et al. (2009) give a number of ways of approximating these values, at least up to a multiplicative constant (which is all we need). An alternative is to maintain a population of feasible and infeasible particles. At each window, draw Xt from the prior, extract XN to give a Markov state then weight by the Markov probability as defined in 2.3.4 divided by the prior probability of drawing XN. The set of feasible particles are then draws from the posterior.
As the number of assimilated windows increases, there may be a lag L where we are happy to assume that the new observations Ωt+1 do not change our beliefs about the model state before the end of the (t − L)th window, i.e. P(τ1:t−L|Ω1:t, θ) ≈ P(τ1:t−L|Ω1:t+1, θ). In this case, it is sufficient to perturb the trajectories of only a finite number of windows into the past when performing the MCMC moves, so we perturb τt−L+1:t+1 and θ while holding τ1:t−L fixed. This means that the dimension of the MCMC problem does not increase as we assimilate more windows. The early parts of the trajectory are then subject to possible particle impoverishment again if we’re using an algorithm with a resampling step, but this time it doesn’t matter because it is, by assumption, not affecting the distribution of the current state.
5.2 Limitations and further work
This algorithm presented here is intended as the proof-of-concept of a novel approach to MCMC sampling from an ABM rather than a production-ready tool for immediate use by practitioners. At present, the bounding of the supports of the agent time-step, boundary conditions and observation likelihood functions by a convex polyhedron must either be done by hand or by third party abstract interpretation tools. Code for commonly encountered scenarios is available on the GitHub repository but this process could be automated for user supplied functions by taking as input a computer program and outputting a set of linear constraints for the support. Similarly, the temperature setting, τ, at present needs to be set by trial-and-error, whereas an automated warm-up or adaptive scheme would improve the algorithm’s efficiency and usability.
As presented, the algorithm is only applicable if the internal state of an agent can be represented in a few bytes. This is because it depends on explicitly representing the probability of each Markov transition at each iteration. But the number of possible transitions is of the order of the number of non-basic variables, and this expands exponentially with the size of the internal state of an agent. So, if agents have a large amount of internal state it will become impractical to represent this distribution explicitly. This problem can be overcome by considering only a subset of non-basic variables for update at each iteration, in a similar way to partial pricing in the simplex algorithm (Maros, 2002). For example, given the current Markov state, only consider updating the non-basic variables that correspond to alternative actions of some agent in the current trajectory at some timestep. In this situation we also can’t explicitly store the coefficients of the linear constraints, C, so in order to calculate the effect of a perturbation on the basic variables we’d calculate this on the fly from equations (7) and (8), the agent timestep function and the observation operators. Finding a good basis of non-basic variables could also be done on the fly for some subset of variables, or alternatively equality constraints could be left as-is as inequalities with equal lower and upper bounds. More work needs to be done exploring the relative efficiency of these options.
