Generating Origin-Destination Matrices in Neural Spatial Interaction Models
Authors: Ioannis Zachos
Mark Girolami
Theodoros Damoulas
Venue: NeurIPS 2024 poster (2024)
Keywords: multiagent systems
neural differential equations
contingency tables
agent-based modelling
markov bases
origin-destination matrix
Primary area: machine_learning_for_social_sciences

Abstract: Agent-based models (ABMs) are proliferating as decision-making tools across policy areas in transportation, economics, and epidemiology. In these models, a central object of interest is the discrete origin-destination matrix which captures spatial interactions and agent trip counts between locations. Existing approaches resort to continuous approximations of this matrix and subsequent ad-hoc discretisations in order to perform ABM simulation and calibration. This impedes conditioning on partially observed summary statistics, fails to explore the multimodal matrix distribution over a discrete combinatorial support, and incurs discretisation errors. To address these challenges, we introduce a computationally efficient framework that scales linearly with the number of origin-destination pairs, operates directly on the discrete combinatorial space, and learns the agents' trip intensity through a neural differential equation that embeds spatial interactions. Our approach outperforms the prior art in terms of reconstruction error and ground truth matrix coverage, at a fraction of the computational cost. We demonstrate these benefits in two large-scale spatial mobility ABMs in Washington, DC and Cambridge, UK.

Generating Origin-Destination Matrices in Neural Spatial Interaction Models
Ioannis Zachos
Mark Girolami 1.2 Theodoros Damoulas 2.3
Department of Engineering, Cambridge University, Cambridge, CB2 1PZ.
"The Alan Turing Institute, London, NW1 2DB, Departments of Statisties & Computer Science, University of Warwick, Coventry, CV4 7AL. iz2300cam.c ac uk
Abstract
Agent-based models (ABMs) are proliferating as decision-making tools across policy areas in transportation, economics, and epidemiology. In these models, a central object of interest is the discrete origin-destination matrix which captures spa- agent trip counts between locations. Existing approaches resort to continuous approximations of this matrix and subsequent ad-hoc discretisations in order to perform ABM simulation and calibration. This impedes conditioning on partially observed summary statistics, fails to explore the multimodal matrix distri- bution over discrete combinatorial support, and incurs discretisation errors. To address these challenges, introduce computationally efficient framework that scales linearly with the number of origin-destination pairs, operates directly on the discrete combinatorial space, and leams the agents' trip intensity through neural differential equation that embeds spatial interactions. Our approach outperforms the prior art in terms of reconstruction error and ground truth matrix coverage, fraction the computational cost. We demonstrate these benefits in large-scale spatial mobility ABMs in Cambridge, UK and Washington, DC, USA.
tial interactions
Introduction
High-resolution complex simulators such agent-based models (ABMs) increasingly deployed to assist policymaking in transportation [10 21], social sciences B g 14 34, and epidemiology
They simulate individual agent interactions governed by stochastic dynamic systems, giving rise to an aggregate, in a mean field sense, continuous emergent structure. This is achieved by computationally expensive forward simulations, which hinders ABM parameter calibration and large-scale testing of multiple policy scenarios 25 Considering ABMs for the COVID-19 pandemic
16 22
[6] as an example, the continuous mean field process corresponds the spatial intensity of the infections which is noisily observed at some spatial aggregation level, while the individual and discrete human contact interactions that give rise to that intensity are at best partially observed or fully latent. In transportation and mobility, running examples in this work, the continuous mean field process corresponds to the spatial intensity of trips arising from unobserved individual agent trips between discrete sets of origin and destination locations [19 [0].
The formal object of interest that describes the discrete count of these spatial interactions, e.g. agent trips between locations, is the iigiinddstinattinn matrix (ODM). It is an I x J (two-way) contingency table T with elements Tig E N counting the interactions of two spatial categorical variables i,j Noo. see Fig. [ It is typically sparse due to infeasible links between origin and destination locations, and partially observed through summary statistics such as table row and/or column marginals due to privacy concerns, data availability, and data collection costs. Operating
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

This work: GeNSIT (n" sample)
ABM (s" simulation)
SIM
Ground truth
T
T
T
A (n)
HW SDE
experasive
Markov
numerical solution
simuuattoon
Basis
A
1 n
4 4 6
sampling
4
4
R
R
2
262
2 2
2
R. Rogion
C 'uummary statistic coustraint D:Discrete Tip d:Coniimoous Aug
C:j: Travel cost from i to
Observed log destination altraction V5 (x marker size)
Figure 1: The ground truth discrete ODM (two-way contingency table) can be reconstructed through either multiple expensive ABM simulations 2 (ABM rectangle) or approximated by continuous representation A coupled with the Harris-Wilke truth can be reconstructed by sampling in the discrete combinatorial space of constrained ODMs conditioned on A (GENSIT rectangle). ABM simulations scale with OMM log(M)) compared to GENSIT which scales with O(JJ), where M > + J is the size of the agent interaction graph.
SDE (SIM rectangle). In the latter, the ground
the discrete ODM level and learning this latent contingency table from summary statistics is vital for handling high-resolution spatial constraints and partial observations such the total number of agents interacting between a pair of locations. It is also necessary for population synthesis in ABMs [15,, which is performed prior to simulation in order to reduce the size of the ABM's parameter space. Moreover, it avoids errors and biases due to ad-hoe discretisation required when working with continuous approximations of the underlying discrete ODM T
Traditional ABMs, Fig. I (Left), simulate individual-leve and spatially granular discrete ODMs at high computational cost which scales least with OM [M)( where M >> I + J is the size of the agent interaction graph. These ODMs are then aggregated using a sum pool- ing operation to regional ones, T(,), whose summary statisties C correspond to observed data. However, the lower-dimensional subspace of contingency tables T satisfying constraints C (c.g. row and column marginals), denoted by Tc. is known to be combinatorially large [1]] and therefore sampling and optimisation in that space is notoriously hard since it requires enumerating all elements of Te. This underlying challenge is the reason why prior art 35 13 32
6 2
3
0 2 O
3
2
55
2 0
27 has been imposing continuous relaxation of the ODM estimation problem, where the continuous approximation of the discrete contingency table restricts inference at the agent trip intensity level A. This leads to quantisations and inefficient rejection sampling, see Fig. 2
0.22.70.2 70.2
1. 41.
1.0.2
A: A T E Te
X Rejected samples Quanitisation
Such intensity-level Spatial Interaction Models (SIMs) [1), Fig. (Right), are derived from entropy maximisation argu- ments, Sec. _ with summary statistics constraints. These models are embedded in the Wilrisn (HW) system of differential equations (20) that describe the time evolution of location attraction and reflects the utility gained from reaching a destination. Coupling SIMs with the HW model introduces an inductive bias that regularises the continuous ODM and fa- cilitates a mean-field ABM approximation. This approximation effectively acts as a cheap ABM surrogate or emulator, facilit- ating faster forward simulations to be run by practitioners [4]. This has tangible benefits to ABM calibration allowing faster exploration of the parameter space. Our enhanced ODM reconstruction error demonstrates GeNSIT's ability to sufficiently approximate the ABM simulator at fraction of the computational cost. However, both the continuous ODM row and column marginals are fixed, the HW model becomes redundant and the SIM can either be greedily approximated through iterative proportional fitting 45 which is sensitive to initialisation or become unidentifiable as the number of parameters scales with both the number of origins and destinations. Furthermore, conditioning on individual continuous ODM cells creates discontinuities as the SIM becomes a piecewise function, also preventing its coupling into the HW model.
Quantisation field
The space Te of 3 x 3 discrete ODMs with summary stat- istics Cr. Sampling on the continu- ous relaxation of Te (A level) with quantisation can lead to either large rejection rates, or poor exploration of the distribution over Tc.
Figure 2:

---

Competing approaches that operate directly on the discrete ODM space and are motivated by econometric arguments are discrete choice models [4]] However, these cannot encode summary statistic constraints without introducing large rejection rates. The work of [ leverages SIMs to sample discrete ODMs but removes intensity constraints through log-linearity assumptions and does not exploit the physics structure, effectively stripping SIMs of their advantages over other choice models. Markov Chain Monte Carlo routines have been devised to address these issues by learning the SIM parameters [33 and the associated discrete ODM table over its entire support [47]. Such routines incur a computational overhead in the order of at least or due to the intractability of the Harris-Wilsee model prior, rendering them prohibitive for large-scale applications. Neural Network (NN) parameter calibration of SIMs has been empirically shown to achieve up to ten-fold speed-ups during training [77 by using numerical discretisation scheme for the Harris- Wilson model as a forward solver. Despite the significant advantages offered by NNs, they operate strictly in the continuous intensity level and cannot generate discrete agent-level ODMs.
In this paper, we introduce a computationally scalable framework named Generating Neural Spatial Interaction Tables (GENSIT) for exploring the constrained discrete ODM space using closed-form or Gibbs Markov Basis sampling while neurally calibrating the underlying physics-driven SIM parameters of its continuous representation, as shown in Fig. 3 Our framework scales linearly with the number of orgiinddstinattinn pairs IJ, which is at least or times faster than MCMC [3] 47 We offer faster ODM reconstruction and enhanced uncertainty quantification compared continuous approaches in seminal works [13 and hybrid (discrete and continuous) approaches [47 in terms of the number of iterations required. It is the first framework that jointly explores the constrained continuous and discrete combinatorial ODM spaces linear in the number of origin-destination pairs time. It does so while outperforming the prior art in terms of reconstruction error and ground truth table coverage while enabling the integration of a broader set of constraints, if these are available.
Our framework has merit beyond ODM sampling in ABMs. The challenge of learning discrete contingency tables constrained by their summary statisties extends to other fields. Contingency tables have been widely studied in multiple instance learning B3
48 and ecological inference 37 36 39]. In Neuroscience one estimates the efficiency, cost and resilience (equivalent to T,y) of neural pathways between pairs of brain regions (i. to understand communication and information processing 38 [18]. Epidemiology also investigates social contact matrices quantifying the number of contacts (Tiy) between types of individuals (i, j) stratified by demographies, such as age 29
2 Spatial Interaction Intensities and Contingency Tables
Consider A agents travelling from / residences (origins) to J workplaces (destinations). The expected number of trips (intensity) between origin and destination j Ai and is unobserved. The average number of agents starting (ending) their journey from each origin i (to each destination j) is:
A+t
Agg. i= I.
A+j
Σ Aij. j= L
(1)
The expected total number of agents travelling is assumed conserved:
A
A+ = A.
(2)
The family of models for intensities A that assimilate any collection of the above constraints called Spatial Interaction Models [45]. The demand for each destination depends on its attractiveness denoted by = (21...., . 2,) R In example. this is the number of jobs available each destination. Let the log-attraction be x == log(z). Between two destinations of similar attractiveness,
J.J agents are assumed to prefer nearby destinations. Therefore, a cost matrix C = (c. is introduced to reflect travel impedance. These assumptions are justified by economic arguments [3] and establish the basis for the agents' utility function. The maximum entropy distribution of agent trips subject to A++ = A yields a totally constrained SIM intensity:
A4+ exp(ax; -By)
Aij
ER. exp(aar.. BCk.)
IJJ
3

---

where CA = (A++} } is the set of summary statistic constraints on A. Henceforth, we set A++ E CA unless otherwise stated. The vector (a, 3) contains the agents' two utility parameters controlling the effects of attractiveness and deterrence on the expected number of trips A. If a grows larger relative to & then agents gravitate towards destinations with higher job availability regardless of the travel cost incurred, and vice versa. Further, if we also fix the expected origin demand A. then we obtain the following singly (also known as production) constrained SIM intensity:
Ait exp(oar -B0x,)
Ar
(4)
Σn exp(ertm - ..8mm)
where CA is expanded to include A. in this case. Moreover, setting CA a doubly constrained SIM
A. A+] yields
Ai = A++ exp(axx - exp(az;--es)000 (i)))()
(5)
where Oi),, D(j) are called balancing factors that ensure that CA are satisfied. The balancing factors introduce 1 + unknown parameters, rendering the intensity model unidentifiable. Alternatively, these factors are approximately recursively using iterative proportional fitting [45], which sensitive to initialisation. Including individual cell constraints in CA breaks the continuity of A as a function of its parameters. For these reasons, the doubly and/or cell-constrained SIMs are prohibitive for use in statistical inference. See App. B.I.I for more information on SIMs as a modelling choice.
We note that additional data at the origin, destination and origin- destinadstnnnnonnn level can be assimilated into SIMs. This can be achieved by incorporating them as terms in the maximum entropy argument used to derive the A functional forms in equations (3 and (5). We note that the SIM's A is equivalent to the multinomial logit (30), which generalises our A construction to accommodate for more data. See App. B.1.2 for a guide on eliciting agent utility functions.
It has been shown that the destination attractiveness z = exp(x) in familics of SIMs is governed by the Harris Willson system of // coupled stochastic differential equations (SDEs) 22 13
dzz dt Ea Kæj + 6) + 0z2 Bj,t, z(0)= =z'
(6)
where e >> 0 is a responsiveness parameter, x > O is the number of agents competing for one job, 8 2 0 is a parameter related the smallest number of jobs a destination can have, > 0 is the standard deviation of SDE's noise, and B, is a J-dimensional Wiener process. The term A+j-Ka + in (6) reflects the net job capacity destination j. If more agents travel to j than there are jobs there (positive capacity), this may signify a boost in j's local economy, which would trigger rise in job availability, and vice versa. The diffusion term stochastically perturbs this trend to account for unobserved events, such as local government interventions affecting employment. By the HW-SDE, the SIM intensity is a stochastic and physics-driven quantity reflecting the agents' average number of trips A between their residences and workplaces. However, the SIM intensity differs from the realised number of trips T agents make. The two notions are connected as follows:
(7) T is the I x J discrete contingency table summarising the number of agents travelling from to j. and Tij 1 Tey ( Asg, Avy V i + Y.j + jr. Any choice of model for TiylAn Au, (say Poisson or Binomial) becomes equivalent upon conditioning on sufficient summary statistics Cr M. We note that the conditional independence of the Tiy 's given the A.y's 's and that T inherits all constraints from A such that every summary statistic constraint in A-space is also applied in T-space, i.e. we always set A+ = [Cr++ICr] and As4 = E[T,+ [Cr) The hard-coded constraints Cr are realisations of the Poisson random variables T AA T A, T45 AA, and therefore no longer random. They can be thought of as noise-free data on the discrete table space. Following the notational convention for A, we define T,+ Σj= Tij. T+j = Σ Tij and T++ Σijj=1 Tij, respectively. The aforementioned summary statistics become Dirac random variables upon conditioning on T. The union of table and intensity constraints is summarised by C. We sometimes drop subscripts from C for clarity. See App. for more information on the notation used.
TiglAi ~ Poisson(Aig),
IJJ
i=
3 Neural Calibration of Spatial Interaction Models
We now introduce our framework (GENSIT). We estimate the parameters of the continuous SIM intensity by employing an ensemble Neural Networks VNN RJ + IR, This allows us to bypass
Codebase found at httpp://gthhb.com//aaaii///////

---

Neural SIM Calibration
A (n-1)
x(n-
W
(())
Tin
W
W
T(2)
Markov Basis
7
32
Patsogool operation with food
Vig 1:
Depeedinne proseet onls in pest bome
D= (D+.y) C)++
Deversintnti
covermatttti
peoduct
(a)
(b)
Figure 3: GENSIT: (a) successive iterations of Alg. for a given ensemble member, (b) plate diagram for every iteration, ensemble member. We propose two schemes: a Joint and Disjoint (see App. B.3 for details). Contrary to the latter, the former passes table T information to the loss (see (b)). We perform optimisation step in the intensity A space and sampling step in T space, with associated complexities OT IJ) and O(JJ). The A arises by the well- well-known family of SIMs 3 coupled with the HW-SDE (6) The T sampling step generates discrete Cr-constrained ODMs contrary to 13 which only operate on the continuous mean-field level A.
the computational challenges of solving the HW-SDE inverse problem within Bayesian framework 13 [77 Conditioned those estimates and for every ensemble member 1,.... E, we solve the HW-SDE to get estimates of the time-evolved log destination attraction & after T time steps using an Euler-Maruyan numerical solver Óиw R - R (28). This allows to incorporate the HW physics model into our parameter estimation without sampling from the SDE's intractable steady-state Boltzmann-Gible distribution, which was the case in 13 47
Instead of log-destination attraction data (y) likelihood, we compute loss operator C(:; DD v) that can assimilate data D from multiple sources on any transformation of T, A, x. We note that v is used to denote the loss hyperparameters (see Tab. 5). Then, the NN parameters (weights and biases) W are updated using pack-propagation using suite of optimisation algorithms (23]. This step requires derivatives of the loss with respect the NN parameters V w L( D, be computed, which
achieved using off-the-shelf auto-differentiation libraries [26]. This gradient is informed by & estimates and therefore by the dynamics the HW-SDE in (6) Hence. the resulting SIM intensity both stochastic and physics-driven. These steps are depicted in Fig. by following the arrows from right left and in steps to 12 Alg. [ We note that our framework is divided into two sampling schemes: a Joint and Disjoint. The former passes T information to the loss operator L, whereas the latter does not. The Joint scheme corresponds to Gibbs sampler on the full posterior marginals (x. T, T.C, D), (e T, C, D) and ( x, C, D). The Disjoint scheme corresponds to collapsed Gibbs sampler where we sample from 0 C)D), x(,,C, integrating out T (see App. B.
C, D) and then from (e, x, C, D) by
The described optimisation routine yields point estimate for A for any given ensemble member and iteration n. We either increase the ensemble size E' to obtain a distribution over A, or we can treat A as a deterministic mapping of realisations of random variables x, whose generalised posterior is: p(x, O exp(-L(x, 0; D, ())(())))))
(8)
where v are loss-related hyperparameters. Each loss evaluation can be treated as a sample from a generalised likelihood (first term), while samples from the prior (second term) are obtained by forward solving @nw. Finally, we can enforce a prior over by appropriately initialising the NN parameters W.
The continuous agent trip intensity A (continuous ODM) samples need to be mapped to discrete agent trips T (discrete ODM), which are required for agent population synthesis. We proceed by
5

---

introducing the necessary machinery to achieve this. We wish to sample from the target measure evaluated as C)).
3.1 Constrained Table Sampling
Let the table cells be indexed X' = ( table value of cell x = (4,j) E X. All non-negative two-way contingency tables T are assumed to be members of a discrete space T. The k-th basis operator 1 Nr + {0, 1] is defined to be
such that T(x) Tij is the
[x(z) (O..... 0,
O.
0,
0).
entry
entry I+J
I+J
entries
Hence, we can define summary statistic operators S T - N linear combinations of the basis operator, that is S&(T) Exxxx (x)) A collection of such summary statistic operators abbreviated by S(T).
Algorithm 1: Generating Neural Spatial Interaction Tables.
O(NE(TJ + IJ))
1: Inputs: evidence: C, C. D, funcs.: C, NNNN, ΦHw, M, μ. hyperparams.: N, E, T, v, K,
2: Outputs: x (1:N), (1:N), (I:N) ((I:N)
3: for each ensemble member E,...,
4: Initialise W T
5: for each iteration 1....,, N do
6:
((m) < Vyw((;; W(n-1)), with y E ED.
Forward solving Neural Net.
(()) + (0(m), Ki K,,,o
7:
x(") -x(("-1 (( Vt= .....,
8:
9:
for each time step 1 l...., do
Incorporating physics.
10:
x < dHwll (x.") (m)
(n)
(())).
Solving Harris-Wilson SDE.
(n) < (m)
11:
12:
A (n) Ac(x(1), O(4) C) using (3)
Computing SIM intensity.
13:
if Joint sampling scheme is used then
14:
I(m) + ((x(1), T(n-1), (m)) D.
else
Computing
16:
L + -(x(x)) ()) Compute vin L.
Neural Net loss.
Update (")) using back pack-propagation. if is not tractable then
Updating Neural Net Weights with SGD.
19:
20:
Sample f uniformly random from M. Find supp(n} such that T(=11) rnfi 2
Gibbs Markov Basis sampling.
Sampling discrete ODM.
22:
(1) .Cr). T(a) ~ pl. A (a) Cr).
T(m) (T((-1) + A
23:
else
}
Closed-form sampling.
24:
Constraints T can be expressed as fixed summary statistics C7 {s sK). where each S is a fixed evaluation Sg (T) with respect to the basis operator. For example, constraint (T. T4.3 can be expressed in terms of the basis operator ( over the entire cell set X.
Definition 3.1. Let Cr ..,,,..,, ,8K } be set of constraints based on summary statistics operators S. A table T Cr-udmissible and only S(T) = (Sx J E Cr.
The subspace of T containing all C2 -dmiissibee x J contingency tables is Te = (T T S(T) Cr). This space contains all discrete ODMs consistent with the aggregate summary statistics Cr. Our goal is to efficiently sample from measure y To subject arbitrary Cr.
The target distribution tractable, i.e. the entire table can sampled directly in closed-form if and only if the universe of summary statistic constraints contains at most one table marginal

---

(row or col. sum). This covers the cases where Cr Tractabbble C TT,, S X, e N)}) U + T_.. (Tx, [ _ X,I N])). Note that the unconstrained case of Cr is handled by the construction in 77) This facilitates sampling T in closed-form and in parallel as shown in step 22 of Alg. [ The most notable cases of tractable distributions are (see App. B.2.
AJA, T++ Multinomial (T++. A/A++);
(10)
TJA, T II
Multinomial A/Act);
(11)
==1
A T
II Multinomial(4g..///
(12)
Each of the distributions above can assimilate cell constraints the form Cr = (Tx, 2 C E
without violating ja's tractability by limiting the support of T A, C. If the constraint set Cr contains
at least both marginals, that is CT intraciable C T {Tx | S X, E NJ)) C7 unttable then the target distribution becomes Fisher's non-central multivariate hypergeometric
T
1.J
AAyA+
Ti
II
(13)
T4 jJ== Tyy
A14
+,==
Direct sampling without rejection from this distribution for arbitrary I, is infeasible .
3.1.1 Markov Basis MCMC
Therefore, we devise an MCMC proposal on Te. Using a suite of greedy deterministic Algorithms
(6]. we can initialise our MCMC with a T E Tc. By virtue of definition 3.2 we guarantee that proposed moves modify the summary statistics in Cg (Condition 1) and that there exists path between any two tables such that any table member of the path is Ca-admissible (Condition 2). The collection C' of K constraints generates K sets of Markov bases MM....,
Our proposal mechanism consists of the universe M 0== Ma- See App. B.2.2 for info Markov Bases. Definition 3.2. A Markov basis M₁ is set of moves fi...., fr. X + Z satisfying:
S& (G) (0) V1 1I Ma and
2. for any two C4-admissible T, T there are fl..... , fl with q E Zyo such that T T + Em=t mf. and T Σm=i mf.. 2 0 for 1 i a A.
In the case of x J ODMs constrained by C, inracttalle .M consists of functions f,.... fr such that x = (11,j1), = (12. E with i # 12, J1 + j2:
7 n if x (21 or I (12.j1); otherwise.
if x (1,, ji) (12,52);
fx)
(14)
A Gibbs Markov Basis (GMB) sampler can now be constructed (see steps 20
22 in Alg.
Proposition 3.1. (Adapted from [II]): : Let μ be a probability measure on Te. Given a Markov basis M that satisfies 3.2 generate a Markov chain in Te by sampling / uniformly random from (1...., L... Let 1) E Z. If the chain is at T E Te, determine supp (1) such that + nff 2 O. Choose
P(n))
II (0))
(m(T(a) + (x)))
and move to T = T + 1if4 for the choice of n, An aperiodic, reversible, connected Markov chain in Te is constructed with stationary distribution proportional to (T)). Proof is found in [1].
We note that when individual cells are fixed (T C), the summary statistics operator is applied over X2 instead of X. Therefore, the size of M is reduced to comply with Definition 3.2 This may shorten the diameter of the Markov Chain's 's state space, leading to better mixing times 42
7

---

Experimental Results
We empirically test framework both synthetic and real-world data from Cambridge, UK and Washington, DC, USA. We compare GENSIT against SIM-MCMC [13), SIM-NN [[7], SIT-MCMC 47 and the Geo-contextual Multitask Embedding Learner (GMEL) (7). See App.
for the computational complexities of these methods. GMEL is trained on a larger set of data that includes, besides the cost matrix C, destination-level urban indicators Y, where the log destination attraction y is a column vector of Y. We validate the generated T and A samples from all methods against the ground truth ODM T using the Standardised Root Mean Square Error (SRMSE) and the 99% high probability region cell coverage probability (CP) (see App. C.2 for definitions), See App. D for experimental protocols, implementation details, comprehensive sensitivity study on our Tramework, and reproducibility using our package gensit
C.1
4.1
Synthetic ODM Scalability
10
Total Compute Time (
12.5
3.2
Intersity Learning Time
(( GENSITPeepn)
Table Sampling Time
(Degan)
GeNSITUont
GeNSITToor
GENBITUJnn))
10.0
SIM-NN
3.0
7.5
2.8
5.0
2.6
0.00
0.25
0.50
0.75
1.00
0.00
0.25
0.50
0.75
1.00
0.00
0.25
0.50
0.75
1.00
x10
X10
x10
Figure 4: Total computation time (left) of GENSIT and SIM-NN
versus the number of
origin-destination pairs IJ with (I x J) equal to 100 x 100, 200 x 200,
1000 x 1000. The two algorithms are run for N 10 iterations with CT = (T.,, T.,, Tx Constraint T Tx means that 50% of table cells chosen uniformly at random are fixed. Total computation time is the sum of the T sampling (middle) and A learning (right) times. Intensity learning and table sampling times are computed for lines 6 12 and lines 19|24 of Alg. respectively. Our framework scales linearly with IJ, which is much faster than SIM-MCMC and SIT-MCMC. SIM-NN does not operate at all in the discrete table space, which explains its faster computational speed.
We begin by offering synthetic exper- iments for varying ODM dimensions x J) and number of agents A, com- paring computation time and ground truth reconstruction error (SRMSE) in
0.0025
(RMSE (anl) (")
(")
0.802
0.015
(FMSE (R
0.005
Figs. and For fixed (I J)
0.000
dimension we sample an intensity uni- formly at random and subsequently generate ground truth table T' with fixed A T4+ by sampling from the Multinomial distribution in (OP. For Figs. and 5b we set T++ 106 while for Fig. 5a we vary this for 0.1, 0.5 and x 1
0.00111
0.35
Number of Agents
1.0
CDM Dimension
(a)
(b)
Figure 5: SRMSE by total number of agents A (left) and
ODM dimension (I x J) (right) of GENSIT's discrete ODM sampling (line 24 of Alg. 1 for N = 10a iterations, fixed intensity and C (I x J) = 150 150. The reconstruction error (SRMSE) scales linearly in IJ and exponentially A T++.
+. Tx On the left we set
Our framework achieves speed-up of (JJ²) compared to previous methods such as SIM-MCMC, SIT-MCMC.
Reconstruction errors grow linearly in dimension (I J) and exponentially number of agents A.
4.2 Real-world ODM Reconstructions and Predictions
We now assess our framework's 's capacity to estimate agent home-to-wan trips in real-world ODM prediction problems in Cambridge, UK 47 and Washington DC, USA (27].

---

4.2.1 Cambridge, UK
In the Cambridge dataset, the ground truth ODM is a 69 x 13 contingency table with 33, 704 agents. Tab. [soows that reconstruction error (SRMSE) and the % of ground truth cells covered by the 99% high probability region of the ODM samples (CP) are significantly improved when operating in T level using our Joint scheme compared to SIM-MCMCan A level. We also outperform the competitive SIT-MCMC approach 77 which can operate in the discrete ODM level, since in the limit of Cr our method can move to high T probability regions much faster than MCMC due to the optimisation of the SIM A parameters. This effectively deflates the Á estimator variance relative to the variance of the A samples in MCMC. Only in the absence of rich Cr data does SIT-MCMC outperform our method (total constrained case) as limits the support of to [0, 2] which acts as regularisation. Letting the support of 0 span the entire I2 allows information carried by stronger A constraints (such as A+-) to permeate to the T space much faster in GeNSIT compared to SIT- SIT-MCMC. The plethora of sources of uncertainty ranging from the HW-SDE (6) to the combinatorial nature of Te suggest that GENSIT is faster in reconstructing ground truth ODMs compared to a fully Bayesian approach such SIT-MCMC.
SIM-NN, which only operate on the
ORIGIN-
TARGET M
SRMSE ENN) T
99%
CP D:W T*
C
DESTINATION
GENSIT
STT.
SIM-
SIM-
GENSIT
STT-
SIM-
SIM-
MATRIX
DISIOINT
TOINT
MCMC
NN
MCCIC
DISTOINT
JOINT MCMC
NN
MCMC
Total
2.43
2,43
0.70
1.08
0.700
23%
89%
42%
constrained
4+
T
2.43
2.43
0.70
30%
31%
68%
Singly
A
0.43
0.43
0.62
0.43
0.63
90%
90%
42%
93%
42%
constrained
A.
0.43
0.43
0.61
93%
93%
72%
This case is not handled by the approach mentioned in the column.
(a) ODMs with closed-form (tractable) T distributions
10
ORIGIN- DESTINATION
TARGET
SRMSE (M
T*
-99% CP (M (M
GENSIT
SIT-
GENSIT
SIT-
MATRIX
DISIOINT
JOINT
MCMC
DISIOINT
JOINT
MCMC
Doubly constrained
Tr.
A
2.30
0.68
1.22
18%
78%
20%
1.15
0.55
0.59
68%
89%
86%
Doubly and 10% coll constrained
T
A
2.50
0.87
1.10
18%
70%
32%
A+4
1.06
0.43
0.56
71%
92%
88%
Doubly and 20%
2.50
0.92
1.06
18%
78%
32%
cell constrained
T A++
T
1.02
0.38
0.51
74%
94%
90%
(b) ODMs with intractable T distribution
13 where conditioning A on C is problematic.
Table I: Ground truth T validation metrics comparing our method against [I3 17 47 in the A and T levels across constraint sets C and a = 0.141 (best) for the Cambridge dataset. On C basis the best metric in T.A spaces is emphusised for each of the two and highlighted between the two. Inference on the discrete table space offers lower SRMSE and higher CP compared to inference in the continuous space. On an ODM basis we obtain the best reconstruction error (SRMSE) and ground truth coverage (CP) in T-space in all but the totally constrained ODM. This is due to SIT-MCMC forcing 0 E [0, 2] instead of R2, which has a regularisation effect. In the absence substantial data, this effect is more pronounced. See Tab. 6 (App. E. for full table across multiple regimes.
Fig. 6 supports this claim by shedding light on the convergence rate of rumning mean estimates of the ground truth T The Joint GeNSIT scheme converges to a mean T estimate much earlier than SIT-MCMC across all C regimes. Mean T estimates are also improved in the Joint GeNSIT compared to SIT-MCMC in confined spaces (doubly, doubly and 10% cell, doubly and 20% cell constrained ODMs). Under the same C regimes CP does not improve significantly as N grows large, which suggests that the variance of samples appears stable as early as N 10 In the Disjoint GeNSIT the information encoded in larger Cr is not propagated to the A updates. As result, no SRMSE or CP improvements are detected the course of N.
4.2.2 Washington, DC, USA
We apply our method to the Washington dataset, where the ground truth ODM is 179 x 179 contingency table with 200, 029 agents. Tab. reports the reconstruction enror (SRMSE) and the % coverage of the ground truth cells (CP) in A and T. Comparisons against SIM-MCMCand SIT-MCMCWere infeasible due to their high computational complexity (300 hours to obtain 500 samples on 32-core NVIDIA GPU). Instead, we leveraged [2 which operates only in the

---

T-,, Tay, SIT-MGMC
T...
T...T...Tyy Joint (GENSIT)
90
T... Joint (GINTIT)
0.7
85
SRMSE
T..T... SIT-MOMC
8
SIT-MOMG
80
T. TT, Joint (GENSIT)
T.,. Joint (GENSIT)
0.5
75
70
20000
40000
60000
80000
100000
20000
40000
60000
80000
100000
N
N
Figure 6: SRMSE and CP for C C= (T (T. T +- Tx } computed cumulatively along
(T.,,
N iterations of Alg. for TiD,C samples for GENSIT (Joint) and SIT-MCMC for the Cambridge dataset. The Joint GENSIT converges to a lower SRMSE faster than SIT-MCMC while achieving better T coverage. The singly constrained ODM (C (T.+)) } has very large support Te given that CP is decreasing with N. See App. for all cases of C and frameworks.
DATA D
TARGET
SRMSE
SSI
199% CP
PRAMEWORK
M
(M T
(::) T")
(M T
GMEL
Y.C
A
2.43 = 0.15
0.38 0.02
5% 1%
SIM-NN
A
2.47 0.00
0.51 + 0.00
24% + 3%
A
2.17 0.00
0.51 0.00
19% + 3%
GENSIT
(Disjoint)
J.C
T
2.39 0.09
0I.13 + 0.01
47% + 1%
A
2.45 0.00
0.50
2% 0%
GENSIT
(Joint)
T
2.37
0.08
0.45 0.01
44% +: 1%
Table 2: Ground truth T validation metrics (mean + std. for E = 10 ensemble size) comparing
our method against [2 17 in the A and T levels for C = (T++,, x and a = 0.141 (best) for better, and vice versa. We achieve the
the Washington dataset. Arrow indicates higher values best error (SRMSE) and ground truth coverage (CP) overall (see bold cells). The observed data D leveraged train GENSIT, SIM-NN is small subset of the data required to train GMEL (y is column vector of Y). See Tab. Appp. E.2 for full table across multiple a regimes.
continuous space by learning a mapping between a large feature space D and T Tab. shows that we outperform both GMEL and SIM-NN in terms of reconstruction error and coverage.
5 Concluding Remarks
We introduced GeNSIT, an efficient framework for jointly sampling the discrete combinatorial space of agent trips (T) subject to summary statistic data and its continuous mean-field limit A. We surmount the limitations of methods which operate strictly on A space 13 17 27
and of methods that ineur large computational cost 47 We accomplish this by introducing the first framework operating on both T, A that scales linearly with the number of origin-destination pairs TJ. We offer enhanced reconstruction error and coverage of the ground truth ODMs in Cambridge, UK and Washington, DC. Although NNs require much larger number of internal parameters to be calibrated relative to MCMC, their embedding physics models regularises this parameter space and prevents over- ever-fitting. A remaining open problem on this front is the assimilation of more complex C structures in agent population synthesis and simulation (see App. B.2.3). since ground truth data is typically partially observed. Our work also relies on the SIM's assumptions about the agents' decision decision-making process, which in practise is unobserved. An examination benefit the applicability of our framework. In terms of our work's 's social impact, policy decisions made from ABMs of social systems could negatively affect individuals, necessitating expert review and ethics oversight.
different agent utility models could
10

---

Acknowledgments and Disclosure of Funding
IZ was supported by UK Research and Innovation (UKRI) Research Council and Arup Group PhD studentship. MG was supported by Royal Academy of Engineering Research Chair grant RC- SRF1718/6/34 and EPSRC grants EP/T0004 4/1, EP/W0058 16/1, EP/V056441/1, EP/V056522/1, EP/R018413/2, EP/R0347 10/1, EP/R004889/1. TD was supported by the UKRI Turing AI accelera- tion Fellowship EP/VO2678X/1. For the purpose of open access, the authors have applied a Creative Commons Attribution (CC-BY) licence any Author Accepted Manuscript version arising from this submission.
References
[1] Alan Agresti. Categorical Data Analysis. Wiley Series in Probability and Statisties. John Wiley & Sons, Inc., Hoboken, NJ, USA, July 2002. ISBN 978-0-471-36093-5 978-0-441-24968-9. doi: 10.1002/0471249688.
[2] Rushil Anirudh, Jayaraman Thiagarajan, Peer-Timo Bremer, Timothy Germann, Sara Del Valle, and Frederick Streitz. Accurate Calibration of Agent-based Epidemiological Models with Neural Network Surrogates. Proceedings of the 1st Workshop on Healthcare Al and COVID-19, ICML 2022, pages 54-62. PMLR, July 2022.
[3] Steven C. Bankes. Agent-based modeling: A revolution? Proceedings of the National Academy of Sciences, (9(suppl_3)))) 3):7199977977 7200, May 2002. doi: 10.10733pnas.077 07208 1299..
[4] David Banks and Mevin B. Hooten, Statistical Challenges in Agent-Based Modeling. The American Statistician, 75(3):235-24) 2222,,2222222222 April 2021. doi: /10.
[5] Hacène Belbachir. A combinatorial contribution to the multinomial Chu-Vandermonde convolu- tion. Les Annales RECITS, 1:27-32, January 2014.
[6] Yvonne M. M. Bishop, Stephen E. Fienberg, and Paul W. Holland. Discrete Multivariate Analysis. Springer, New York, 2007. ISBN 978-0-387-72805-6.
[7] Pier Giovanni Bissiri, Chris Holmes, and Stephen Walker. A General Framework for Updating Belief Distributions. Journal of the Royal Statistical Society Series B: Statistical Methodology, 78(5):1103 130,, November 2016. ISSN 1369-7412, 1467-9868. doi: 10.1111/rssb.. /tll1//srsb..1588
[8] Eric Bonabeau. Agent-based modeling: Methods and techniques for simulating human systems. Proceedings of the National Academy of Sciences, (tsssuppp3)7788072255 May 2002. ISSN 0027-8424, 1091-6490. doi: 10. 10.10733/nass.08208000
[9] Luis Carvalho. A Bayesian Statistical Approach for Inference on Static Origin- -eetinnttninnoio Matrices in Transportation Studies. Technometrics, 56(2):225-237, April 2014. ISSN 0040- 1 1706, 1537-2723. doi: 10. 11010000001110662113338888
[10] Andrew Crooks, Christian Castle, and Michael Batty. Key challenges agent-based modelling for geo-spatial simulation. Computers, Environment and Urban Systems, 32(6):417-430, November 2008. ISSN 01989715. doi: 10. 10.10116/..ommm htt.::..ommmeauuuusy.2000..00.
[11] Persi Diaconis and Bernd Sturmfels. Algebraic algorithms for sampling from conditional distributions. The Annals of Statistics, 26(1), February 1998. ISSN 0090-5364. doi: 10.1214/ aos/1030563990.
[12] Daniel R. Dooly, Qi Zhang. Sally Goldman, and Robert A. Amar. -ustiple-Instance Learning of Valued Data. Journal of Machine Learning Research, 3(Dec):651-66) 2002. ISSN ISSN 1533-7928.
[13] L.. Ellam, M. Girolami, G. A. Pavliotis, and A. Wilson. Stochastic modelling of urban structure. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences. 474 (013))20170700 May 2018. ISSN 1364-5021, 1471-2946. doi: 10.:1001100 00.1098/rspa.2017.000
[14] Doyne Farmer and Dunean Foley. The economy needs agent- modelling. Nature, 460 (7256):685 68--666, August 2009, ISSN 1476-4687. doi: 10. 10.1038/460685..

---

[15] Bilal Farooq, Michel Bierlaire, Ricardo Hurtubia, and Gunnar Flötteröd. Simulation based population synthesis. Transportation Research Part B: Methoddoogical,, 58:243-22 263,, December 2013. ISSN 0191-2615. doi: (..10166..bb.2013.0000
[16] Neil M. Ferguson, Derek A. T. Cummings, Christophe Fraser, James C. Cajka, Philip C. Cooley, and Donald S. Burke. Strategies for mitigating an influenza pandemic. Nature, 442(7101): 448 452, July 2006. ISSN 1476-4687. doi: 10. 1038/nauuur044795.
[17] Thomas Gaskin, Grigorios A. Pavliotis, and Mark Girolami. Neural parameter calibration for large-scale multiagent models. Proceedings of the National Academy of Sciences, 120(7): doi: 10. 1..1073/pnas.2211441.
e2216415120, February 2023. ISSN 0027-8424, 1091-6490.
[18] Joaquín Goñi, Martijn van den Heuvel, Andrea Koniigboon,,,,,reer, Nieves Velez de Mendizabal, Richard F. Betzel, Alessandra Griffa, Patric Hagmann, Bernat Corominas- Lean-Philippe Thiran, and Olaf Sporns. Resting- brain functional connectivity predicted by analytic measures of network communication. Proceedings of the National Academy of Sciences, 111(2):833-83) January 2014, doi: 10. 10.1003/pnas.. 13 15529111.
[19] Sam Griesemer, Defu Cao, Zijun Cui, Carolina Osorio, and Yan Liu. Active sequen- tial posterior estimation for sample-efficient simulation-based inference. In The Thirty- eighth Annual Conference on Neural Information Processing Systems, 2024. URL https
//openeeiiew nww. wtw.foornnndduuusuuaoo
[20] B Harris and A G Wilson. Equilibrium Values and Dynamics of Attractiveness Terms in IoaductionnCCnstrainng Spatial-Interaction Models. Environment and Planning A: Economy and Space, 10(4):371-388, April 1978. ISSN 0308-518X, 1472-3409.
doi: 10. 1068//
[21] Andreas Horni, Kai Nagel, and Kay W. Axhausen. The Multi-Agent Transport Simulation MATSim. Ubiquity Press, August 2016. ISBN 188-1-909188-75--0 doi: 10.5334/baw.
[22] Cliff C. Kerr, Robyn M. Stuart, Dina Mistry, Romesh G. Abeysuriya, Katherine Rosenfeld, Gregory R. Hart, Rafael C. Núñez, Jamie Cohen, Prashanth Selvaraj. Brittany Hagedorn, Lauren George, Michal Jastrzebski, Amanda S. Izzo, Greer Fowler, Anna Palmer, Dominic Delport, Nick Scott, Sherrie L. Kelly, Caroline S. Bennette, Bradley G. Wagner, Stewart Chang, Assaf P. Oron, Edward A. Wenger, Jasmina Panovska-Griifiths, Michael Famulare, and Daniel J. Klein. Covasim: An agent-based model COVID-19 dynamics and interventions. PLoS Computational Biology, 1707)=e1009110 149, July 2021. ISSN 1553-734X. doi: 10.1371/ journal. pcb....ibii.14.4
[23] Diederik P. Kingma and Jimmy Ba. Adam: A Method for Stochastic Optimization, January 2017.
[24] Jeremias Knoblauch, Jack Jewson, and Theodoros Damoulas. An Optimization-centric View on Bayes" Rule: Reviewing and Generalizing Variational Inference. Journal of Machine Learning Research, 23(132):1-109 2022. ISSN 1533-7928.
[25] Robert Lempert. Agent-based modeling as organizational and public policy simulators. Pro- ceedings of the National Academy of Sciences, 99(supp__3:71155-77) May 2002. doi: 10.1073/pmas 072079399.
[26] Shen Li, Yanli Zhao, Rohan Varma, Omkar Salpekar, Pieter Noordhuis, Teng Li, Adam Paszke, Jeff Smith, Brian Vaughan, Pritam Damania, and Soumith Chintala. PyTorch distributed. Proceedings of the VLDB Endowment, 13(12):00553311 August 2020. ISSN 2150-8097. doi: 11144778/3455778.411111
[27] Zhicheng Liu, Fabio Miranda, Weiting Xiong, Junyan Yang, Qiao Wang, and Claudio Silva. Learning Geo-Contextur Embeddings for Commuting Flow Prediction. Proceedings of the AAAI Conference on Artificial Intelligence, 34(01))808-811)) April 2020. ISSN 2374-3468. doi: 10. v.00//aai..34001...
[28] Gisiro Maruyama. Continuous Markov processes and stochastic equations. Rendiconti del Circolo Matematico di Palermo, 4(1):48-94 January 1955. ISSN 1973-4409, doi: 10. 10.1 BF02846028.
12

---

[29] Dina Mistry, Maria Litvinova, Ana Pastore y Piontti. Matteo Chinazzi, Laura Fumanelli, Marcelo C. Gomes, Syed A. Haque, Quan-Hui Liu, Kumpeng Mu, Xinyue Xiong. M. Elizabeth Halloran, Ira M. Longini, Stefano Merler, Marco Ajelli, and Alessandro Vespignani. Inferring high-resolution human mixing patterns for disease modeling. Nature Communications, 12(1):
323, January 2021. ISSN 2041-1723. doi: 10. (0.1038/$14467 1467-020-205444
[30] Peter Nijkamp and Aura Reggiani. Entropy, spatial interaction models and discrete choice analysis: Static and dynamic analogies. European Journal of Operational Research, 36(2):
186-196, August 1988. ISSN 03772217. doi: 11101660377-2221(88899998
[31] Jim Pooler. An extended family of spatial interaction models. Progress in Human Geography, doi: 10.1177/03091135594 800102.
18(1):17-39, March 1994. ISSN 0309-1325, 1477-0288.
[32] Nastaran Pourebrahim, Selima Sultana, Amirreza Niakanlahiji, and Jean-Claude Thill. Trip distribution modeling with Twitter data. Computers, Environment and Urban Systems, 77: 101354, September 2019. ISSN 0198-9715. doi:
[33] Soumya Ray and David Page. Multiple Instance Regression. In Proceedings of the Eighteenth International Conference on Machine Learning, ICML '01, pages 425 4355,,42,550 San Francisco, CA, USA, June 2001. Morgan Kaufmann Publishers Inc. ISBN 978-1-55860-778-1.
[34] Carl Orge Retzlaff, Martina Ziefle, and André Calero Valdez. The History of Agent-Based Modeling in the Social Sciences. In Vincent G. Duffy, editor, Digital Human Modeling and Applications in Health, Safety, Ergonomics and Risk Management. Human Body, Motion and Behavior, pages 304-319, Cham, 2021. Springer International Publishing. ISBN 978-3-030- 77817-0. doi: 1111007/978-3-00--7788
[35] Caleb Robinson and Bistra Dilkina. A machine learning approach to modeling human migration. In Proceedings of the 1sl ACM SIGCAS Conference on Computing and Sustainable Societies, Compass 18, New York, NY, USA, 2018. Association for Computing Machinery. ISBN 978-1-4503-5816-3 doi: 10.1145/32098 11.229988.
[36] Ori Rosen, Martin Tanner, Gary King, Ori Rosen, and Martin Tanner. Ecological Inference: New Methodological Strategies. Cambridge University Press, New York, 2004.
[37] Alexander A. Schuessler. Ecological inference. Proceedings of the National Academy of Sciences, 96(19):10578-1051 September 1999. doi: 10. 10.10733pnas.9..11 19.10577.
[38] Caio Seguin, Olaf Sporns, and Andrew Zalesky. Brain network communication: Concepts, models and applications. Nature Reviews Neuroscience, 24(9)5557-574, September 2023. ISSN 1471-0048. doi: 10. 00.10338441583-033-007114
[39] Daniel R Sheldon and Thomas Dietterich. Collective Graphical Models. In Advances in Neural Information Processing Systems, volume 24. Curran Associates, Inc., 2011.
[40] Stephen G. Simpson. Ordinal Numbers and the Hilbert Basis Theorem. The Journal of Symbolic Logic, 53(3):961-974, 1988. ISSN 0022-4812. 12. doi: 10.2307/2274585.
[41] Gabriel Spadon, Andre C. P. L. F. de Carvalho, Jose F. Rodrigues-Jr. and Luiz G. A. Alves. Reconstructing commuters network using machine learning and urban indicators. Scientific Reports, 9(1):11801 August 2019. ISSN 2045-2322. doi: 10. (038//41198-011)) -98-991-4
[42] Caprice Stanley and Tobias Windisch. beat-bath random walks with Markov bases. Advances in Applied Mathematics, 92: 143, January 2018. ISSN 01968858. doi: 10. 10.1016/j.amm.2117 08.002.
[43] Thomas M. Sutter, Laura Manduchi, Alain Ryser, and Julia E. Vogt. Continuous Relaxation For The Multivariate Noncentral Hypergeometric Distribution. ICIR2022 Workshop on the Elements of Reasoning: Objects, Structure and Causality, March 2022.
[44] Kenneth Train. Discrete Choice Methods with Simulation. Cambridge University Press, Cambridge, 2 edition, 2009, ISBN 978-0-521-74738-7,
13

---

[45] A G Wilson. A Family of Spatial Interaction Models, and Associated Developments. En- vironment and Planning Economy and Space, 3(1):1-32, March 1971. ISSN 0308-518X, 1472-3409. doi: 10. 10.1068/a33001.
[46] Tobias Windisch. Rapid Mixing and Markov Bases. SIAM Journal on Discrete Mathematics, (144):2330-24433 January 2016. ISSN 0895-4801, 1095-7146. doi: 10.1137/15M1022045.
[47] loannis Zachos, Theodoros Damoulas, and Mark Girolami. Table inference for combinatorial origin-destination choices in agent-based population synthesis. Star, 13(1))e656, 2024. ISSN
2049-1573. doi: 10. 10.1002/sta4.656.
[48] Yivan Zhang, Nontawat Charoenphakdee, Zhenguo Wu, and Masashi Sugiyama. Learning from Aggregate Observations, In Advances in Neural Information Processing Systems, volume 33, pages 9993-8055. Curran Associates, Inc., 2020.
14

---

Appendices
Table of Contents
Nomenclature
16
Theory and Methodology Foundations
17
B. Continuous ODM
17
B.1.1
Spatial Interaction Modelling Choice
17
B.1.2
Eliciting Agent Utility Functions
17
B.2
Discrete ODM
17
B.2.1
Target Table Distributions
17
B.2.2
Markov Bases Preliminaries
18
B.2.3
Incorporating Additional Constraints
18
B.3
GENSIT
19
B.3.1
Disjoint versus Joint scheme
19
Method Comparisons
19
Computational Complexities
19
Validation Metries
19
Experimental Settings
20
D.I Implementation Details
20
D.2 Experimental Protocols
21
Auxiliary Experimental Results
22
E.
Cambridge, UK
22
Washington, DC, USA
23
Reproducibility
23
Cambridge. UK
24
Washington, DC. USA
26
15

---

A
Nomenclature
We provide Tab. of notation (see Tab.
used throughout the paper.
Notation
Definition
Domain
R
Real numbers.
N
Natural numbers including O.
2
Integer numbers including O.
Number of origin locations
Rgg-
J
Number of destination locations.
Nya
A
Number of agents modelled.
Nya
T
Discrete destineetnn matrix (two-way contingency table).
NTXJ
T
Ground truth discrete orggin-destination matrix.
NTXJ
Continuous origin-de drgiindestinattin matrix.
R
X
Set of ODM cells whose values are fixed/conditioned.
(1,J)))
Spatial interaction model parameter controlling the eflect of log destination attract- iveness on the expected number of agent trips.
Ryo
Spatial interaction model parameter controlling the effect travel impedance on the expected number agent trips.
8
R>0
Harris Wilson (HW) SDE responsiveness parameter.
Ry0
HW SDE parameter controlling the smallest number of jobs a destination can have.
Ry
8
HW SDE parameter controlling the numher agents competing for one joh.
Ry0
0
HW SDE parameter controlling the diffusion term's s noise variance.
R>0
B,
Wiener process.
R
0
Parameter vector (a, 3) learned by the Neural Network.
R5
Destination attraction.
R5
x
Log destination attraction.
R"
y
Observed log destination attraction (log number of jobs).
R
Y
x I' matrix of destination location features including x as a column, where F is the total number of destination features.
RJxP
C
Ix J cost matrix reflecting agent travel impedance.
R
D4+
Ohserved total travel cost by origin location.
R
Observation data used the Neural Network's loss function.
P,y.
D.y.
T T.+
Number of agents travelling from each origin (random, deterministic).
N
T T+-
Number of agents travelling to each destination (random, deterministic).
N
T++.
Total number agents (random, deterministic).
N
T++
Tx. Tx
Set of discrete ODM cell values that fixed/conditioned (random, deterministic).
NI
P(T T
Cr
Set of noise-free data constraints in discrete T space.
T.,Txx
A+
Expected number agents travelling from each origin.
R30
A+
Expected number of agents travelling to each destination.
R50
A++
Expected total number of agents.
Rs0
PHA
CA
Set of noise-free data constraints in continuous A space.
or
P((A+ A+
C
Set of constraints in discrete T and continuous A spaces Cr U U..
T
Space of all non-negative discrete ODMs.
Te
Space of all non-negative discrete ODMs constrained by Cr.
P
Powerset operator.
S
Summary statistic operator.
0
Time/space complexity operator.
P
Probability mass function. Expectation.
E
Hadamard product.
Loss operator using the Neural Network.
Loss operator hyperparameters (as GA).
ONN
Neural Network forward solver.
W
Neural Network weights.
R 1)
DRW
Buler-Maruyama numerical solver for the Harris- Wilson SDE.
M
Markov Basis moves used in Gibbs sampling of discrete ODMs.
M
Target ODM T or A.
Table 3: Notation used in the paper.
16

---

B
Theory and Methodology Foundations
B.1
Continuous ODM
B.1.1 Spatial Interaction Modelling Choice
A fundamental assumption is that the continuous agent trip intensity A is governed by SIM. Therefore, agents are discouraged from travelling long distances due high travelling costs incurred and encouraged to travel to locations that offer high utility, such as areas with high job availability. This means that highly irregular or non-uniform trip data might not be captured well by A. However, operating at the discrete ODM level allows us to assimilate trip data as constraints. This shrinks the ODM support and enhances spatial imputation of missing data much better than doing so the continuous A level. In the limit of trip data constraints, we are guaranteed to improve our estimated of missing trip data regardless of their spatial distribution and presence of outliers. This is theoretically founded [I] and empirically shown in Tabs. 12 of our Experimental Section. In poorly constrained settings, we rely more on A for spatial imputation, which increases reconstruction error. Outlier presence will be reflected in the row/column sum constraints facilitating table sampling in that region of the ODM support. If no row/column sum constraints are imposed, exploring such regions of the support would be significantly hindered.
B.1.2 Eliciting Agent Utility Functions
In the ODMs considered we have constructed agent utility functions that leverage the log destination attraction each destination location x and the total distance (cost) of travelling between any origin and destination C. Additional observations at the origin, destination and origin-destination level can be assimilated by incorporating them as terms in the maximum entropy argument used derive the functional forms of A in equations (3) (4) and (5). The SIM intensity is equivalent to the multinomial logit [30), which allows us define an arbitrary utility model inside the exp in the numerator of these equations. For example, one might want to include data E R' in addition to the
existing x, C in the totally constrained SIM (total A4+). Then, we need maximise:
I Nig (log(Ai)) - 0; oay + B6,g) - pl
I.J
Σ
2 Ais A4+),
ij
where μ is the Lagrange multiplier. This yields
Agj = A (explo; + ozz - Bctz))
Σnj exp(o; + oatz - BCiz)
B.2
Discrete ODM
B.2.1 Target Table Distributions
The target distributions y (T A, Cr) we sample from in Alg. are listed below in Tab.
We note that the quantity wiy Ar in the last three rows of Tab. 4 is called the odds ratio and encodes the dependence between rows (origins) and columns (destinations). SIMs incorporate this dependence spatially in the cost matrix C. The case for spatial independence is achieved when 3 destination. Fisher's non-central multivariate hypergeometric reduces to its central version if and only if wiy = 1 [T] The central version corresponds to a uniform distribution over Te. Additionally, the normalising constant of Fisher's 's non-central multivariate hypergeometric is partition function involving a sum over all elements of Te satisfying the conditioned margins. By virtue of the extension to Chu-Vandennn Vandermonde's theorem 5 proved in [4], this normalising constant can be computed in constant time using
A45
33. ④ This translates to the travel cost having no effect on the agents' choice of
Ti
II
T T
Wij W4j
(15)
17

---

ODM
Cr
p (T) Cr
DISTRIBUTION
TRACTABLE
Unconstrained
Totally
Constrained
T++
11
Agg
Multinomial
Singly
T
Product Multinomial
Construined
Fisher's s non-central multivariate
Doubly Constrained
T,,, T-.
[T. Ru R. Tes!
24 -T-- Tay!
hypergoometrie
Doubly Constrained with 10% cells fixed Doubly Constrained with 20% cells fixed
T++, T+., Tx;
Tie Txa
Table 4: List of target distributions with their associated constraints Cr.
B.2.2 Markoy Bases Preliminaries
Theoretically analysing the full joint framework is challenging future work, however we offer some additional theoretical insight. We encourage the reader follow the references of this section for more technical details.
The fundamental theorem of Markov Bases introduced in the seminal work of Diaconis and Sturmfels 11 establishes an equivalence between a Markov Basis and the generator of an ideal of polynomial ring. By virtue of the Hilbert basis theorem [40, any ideal in polynomial ring has finite generating set. Therefore, there exists a finite MB for the class of inference problems we are dealing with. This MB connects the fiber, i.c. the support of all discrete ODMs satisfying summary statistics in the form of row, column sums and/or cell constraints. By Proposition 3.1 we can construct Gibbs MB sampler that will converge to Fisher's non-central hypergeometric distribution on the fiber finite time.
Moreover, MB sampler converges its stationary distribution in most A steps, where A is the total number of agents [46]. Although this result is not directly applicable to our framework in its given form, conjecture that distribution as evidenced by our experimental results. An MB connects tables (elements) of the fiber graph (discrete state-space of tables). The fastest convergence of an MB sampler is achieved by an MB that has the smallest possible diameter, i.c. the longest path between two elements of the fiber graph.
be extended
+1 and Fisher's central hypergeometric
B.2.3 Incorporating Additional Constraints
Incorporating additional structural constraints remains open problem. The case for sparse matrices can be handled through the cell constraints Txx for X' S X. Symmetric ODMs (e.g. an adjacency matrix of bipartite graph) be explored as follows. The subspace of To consists of square matrices with equal row and column sums. We construct an MB on Te by applying the following modification to equation (14):
+n if E ((11,)), (12, (fi h)), (j2, 12)) -1 if E (2), (12,)), (fr, i1), (2,d2))
(x)
otherwise.
First, this guarantees that every 2 move in the original M that modifies a cell along the diagonal is by definition symmetric as to the other cell it modifies, meaning if cell (i, is modified then so is i). Second, without loss of generality assume that a move modifies a 2 section in the upper triangular section of the ODM. This ensures that the same move is applied symmetrically to the lower triangular section of the ODM. In both cases, any move will guarantee that the ODM will be symmetric after removing any duplicate Markov Bases generated in the process.
18

---

B.3 GENSIT
B.3.1 Disjoint versus Joint scheme
The Disjoint scheme consists only of loss terms that depend directly fully observed data D (either log-destination attraction y or total distance agents travelled by origin location D.+). In contrast, Joint scheme consists of loss terms that either depend the same fully observed data and on the partially observed table T through the table marginal T x, C, D.
The Joint scheme is an instance of Gibbs sampler on the full posterior marginals (x. T, C,C, D), x T.C, D) and (0, x, C, D). The Disjoint scheme is an instance collapsed Gibbs sampler where we sample from 0 (x. C. D). x ( D) and then from ( C, D). This means we integrate out T by means of p [0)x C, D) En P(0 x, PT Er T, P) D). Therefore, we use the Joint scheme when we have reason to believe that the covariance between 0, x and is small. This would be the case when agent trip intensity is influenced by both the Harris-Wilson SDE and the realised number of agent trips. In contrast, we use the Disjoint scheme to accelerate convergence by reducing the covariance between 0, and T. This would be the when the agent trip intensity is governed only by the Harris Wilson SDE and not by the realised number of agent trips.
D) and p(x(0, =
C Method Comparisons
In this section, we compare our method's 's computational complexity to competitive approaches and list the validation metrics employed for empirical comparisons in the main paper.
C.1 Computational Complexities
Let N be the number of iterations, and 7, be the number of origins, destinations.
GENSIT: O(N + IJ)), where is the number of time steps in the Euler-Maruyama solver, E is the ensemble size. We set 1, and E =
2. SIM-MCMC [3]: O(N + (2)) (low a regime) and O(NIJLKn,,)) (high o
regime), where T. is the number of leapfrog steps in Hamiltonian Monte Carlo, 1 < K N is the number of stopping times, and np, are the number of particles and temperatures used in Annealed Importance Sampling, respectively. Typical ranges are np [10, 100), E 30), and L E 10) and E E 10).
3. SIT-MCMC 47
O(N J²)) (low noise regime) O(NIJLKn,m)) (high noise regime). See items 1,2 for details.
4. SIM-NN I17 O(NE(FJ)) See item 1 for details.
5. GMEL 27 (i +ID IJ)) where Y ((++J) the feature matrix for all locations, D is the number of features per location, n. < D is the number of nodes per layer, and ny the number of hidden layers. In the DC dataset, we have + 179 since every origin is also destination (we do not count them twice). Typical ranges include ny E [1, 10) and n, E [1, D].
C.2 Validation Metrics
We leverage the following three validations metrics employed throughout the literature 227 47 35 32). The Standardised Root Mean Square Error (SRMSE) between the target ODM M (table T or intensity A) samples and the ground truth T is equal to
Eij
(1:N) Ic, (E
-T;) Tz²)²
E[M Ej
(1:N) [C,
-1
SRMSE M (1:N) T =
1J
1J
19

---

and the Sorensen Similarity Index (SSI) is equal to
I.J 2 min
E[M (1:N)
Cc, D), Ti3
(SN( TT) TJ
(1:N BLM, IC, DJ Tú
whose domain is [0, 1] and values closer to imply a better ODM fit. Finally, the coverage probability
of the ground truth table T is caleulated by first identifying the lower by (T (1:N)) and upper
Ug T ) boundaries of the high probability region containing 9% of the total mass for each table by cell (i,j) X. Then, the 9% cell coverage probability is equal to
1 Lug (M (M (C, S T U. (M.l }
CPg(M "T) Σ
D Experimental Settings
In this section, we detail the implementation details and experimental protocols used to obtain the results in the main paper.
D.1 Implementation Details
Regarding the HW-SDE, an Euler- Maruyaaa numerical solver Φиw is employed throughout the paper with time discretisation step of At 0.01 and number of steps 1. The low and high SDE noise levels correspond to σ = 0.014 and a = 0.141, respectively, to establish a comparison against set the responsiveness parameter 5 relating to the job availability of a destination where no agents travel to 0, as in [17] The job competition parameter и is uniquely determined by Σ/=, exp(a;)) Moreover, the same transportation
[3 17 (7). Following [B 47,
1, and the parameter
46+
network distance-based cost matrix C as [7 is used.
The set of observation data D may include the ob- served log destination attraction y E IR and the total distance travelled from each origin D++ (T* @ C) E TR For Cambridge, both y, D.+ have been sourced from the UK's ' population census data- set provided by the Office of National Statistics. Their spatial resolution is regional: middle and lower super output areas for y, D.4. respectively. For DC, we only have access to a feature matrix Y E (R++J)xXD from which we extract a column
Input layer
Hidden layer
Output
layer
Linear
y/1
Abe.
* a
as y.
Our NN is a multi-layer perceptron with one hid- den layer, implemented in PyTorch (26) and depicted in Fig. 7 The input layer set to the observed log-destination attractions learning the that generates the observed physics y. The output layer is two-dimensional due to the parameter vector E IR?. For both datasets we set the number of hidden layers nodes to 20. The hidden, output layers have a lin- and absolute activation functions, respectively. The latter is important in guaranteeing that our es- timates are always positive. The NN parameters W are initialised by sampling uniformly over the region [o, 4] (J= 1)x 20+21x2
B
E Rdd since we are
yJ
Liuesu
and number of
Blas
Blas
21
2 # Nodes
(J+1)
x 20
(20
x
# Weights
Figure 7: Visual depiction of the Neural Net- work architecture used in GENSIT for both weal-world datasets. The size of the weights W scales linearly in the number of destina- tions J.
We use the Adam optimizer 23 with 0.002 learning rate, Bias initialised uniformly at [0, 4). We follow
13 47 in fixing σa 0.03 and or, aA to 0.07 to
20

---

reflect 1% and 3% noise levels, i.e. o/ (Jg = 3%. We assume are observations from the SDE's stationary distribution. hence our batch size is one. We initialise the Buler-Maruyann solver at y and run for r = 1 and step size Al = 0.001. At equilibrium the only change in the log-destination attractions is attributed to the SDE's diffusion.
The SDE is encoded in the NN as follows. The Buler-Maruyama solver @Hw provides forward solutions from
different time steps. The loss function computes the discrepancy between the forward solution x after steps and the observed data at the stationary equilibrium of the SDE, allowing us to encode the physics directly into the NN. more complicated architecture with larger number of weights would potentially lead to overfitting the SDE parameters 0 to compromising our framework's 's generalisability.
NAMB
OPERATOR
EVALUATION
D
SCHEME
Negative log destination attraction
[ D,
log(as) + Eim (y7x) log (T/TT) | (A/A
y
Od
Disjoint
Negative log table
((T,A)
Joint
Total T-based distance travelled by origin
log(os) + (Tan
C(T; D,M)
D
ov Joint
ELL
log(Oa))
Total A- A-based distance travelled by destination
C(A; D..)
D
a Disjoint
Joint destination attraction. table loss Joint destination attraction, total T based distance loss
T,A; (D,)
(x:; ((T,A)
y
dd
Joint
L(x,T; D,y)
4(x: D. C D,v)
y,
0a,0r
Joint
Joint destination attraction, total
A-based distance loss
C(x, A; D,v)
[(x D,v) C(A; D,
y, D.
Ca.
Disjoint
Joint tuble, total T-based distance loss Joiat table, total A-based distance loss Joint destination attraction, table, total T-based distance loss
(TT D.Y)
(TT.)) D,,)
D
dy
Joint
CTT,A; D,v)
(T, LA; D.v)
D
ax
Joint
D(x; (TT,A)
Cx,, A, D.
L(T D,v)
y, D.
da,
Joint
Joint destimation attraction, table, total A hased distance loss
L(x, T,A D.v)
L(x: D.v) (TT.A) L(T: D,v)
y.D...
0a.dx
Joint
Table 5: List of loss function operators C and their observation data requirements D.
A list of all loss functions used to train the NN is provided in Tab. 5 The Joint and Disjoint GeNSIT schemes in Figs. and 9 use the loss operators L (x; D, v). and 2 xx T, D,v). The observed total distance travelled Is assumed be more noisy than the observed employment data since the latter is a more volatile quantity than the latter and is more prone to errors. We note that SIM-NN in [17 was under the same NN configuration except for the destination attraction loss L (x: D2 for which the authors used the L-2 error. This is equivalent to our loss Tab. 5 in the limit of 02 -> 0. Regarding, table sampling, the observed cells in the Doubly Constrained with 10% cells fixed and Doubly Constrained with 20% cells fixed cases are chosen uniformly at random over the discrete cell support [1, I] x J]. All intractable distributions in Tab. 4 are sampled using Gibbs Markov Basis described in Proposition 3.1 A variant of maximum entropy iterative proportional fitting 66] is employed to initialise tables in Alg. T
D.2 Experimental Protocols
Figures based on synthetic experiments are produced as follows. We generate ground truth tables (100 x 100), (200 x 200), (300 x 300), (400 x 400), (500 x 500), (600 x 600), (700 x 700), (800 x 800), (900 x 900), (1000 x 1000) as explained in the main paper. For each of these tables, Algorithm and SIM-NN are run for N 10 A 105, and CT { T4., Tx } while monitoring The computation times to produce Fig. 4 Additionally, we generate ground truth tables of sizes
(50 50), (100 100), (150 150), (200 200), (250 x 250) with total number of agents A 10%
Cr {T T Asos and run line 24 of Alg. D to produce Fig. 5b Finally, we generate 150 x 150
ground truth tables with A = 0.01, 0.1. 0.25, 0.5. x 10c and run the same line of Alg. Jusing
Cr (T T X and N 10a create Fig.
5a
Alg. run for N 10s iterations and E = ensemble size to produce the data in Figs. 8 10 In
Fig. computational budget of N F. = 10a samples is fixed and the following schedule (N, F) used: (10, 1000), (50, 200), (100, 100), (500, 20), (1000, 10), (5000, 2), (10000, 1). In this case,
21

---

estimators of the likes of E [- | are computed over both N and E by coupling all samples across members of the ensemble. The Standardised Root Mean Square Error and R% High probability region cell Coverage Probability metries are computed in the same fashion as in [47]. Finally, in Tab.
1 the latent samples A (I:N), ( have been trimmed by applying burning and thinning of 100 and 100 samples across every member of the ensemble.
We follow the setup of 13 to initialise SIM- MCMC SIT-MCMC and ensure comparabil- ity with our models. The computational budget in Tab. | and Figs. |0 8 is set to N = 10 SIM 1.025, 5 = 0.0128 and ad 3% x log(J). For up- dates the acceptance is monitored to be at least 90%, while the 0-update's 's acceptance ranges from 30% to 70% depending on the size of the constraint data C. In the high-noise sampling scheme's importance sampling of the normal- izing constant, the number of particles is set to 100 and uniform temperature schedule of 50 inverse temperatures is employed. The percent- age of positive signs is maintained around 75% and above.
V Te., af 0,002
T. ... 0,014
T.., T.,,T.,Ta, 2001
Te., T.., Try o -014
T-14-
T,,. 0.141
Te, Tu.
T..,T.,, Tr,, o =014
model parameters are fixed to e
Ti,
T,.. T.,, loarnd
Tr,Ta,, learned
T.., Ta,,T, kammdd
1.8-
1.6-
1.4
1.2-
SRMSE
1.0-
0.8-
In the case of the Washington DC data, we em- ploy the same train/test/validation test split as in (27). In terms of hyperparameter optimisa- tion for GENSIT, we leverage the same archi- tecture as in the Cambridge data and optimise the learning rate on the validation set. To en- sure a meaningful comparison, we also optimise the learning rate and multitask weights of the GMEL framework on the validation set. We each method appearing in Tab. for an semble of size E- 10 obtain error bars of this ensemble.
0.6
0.4
0.2
0.0
LEnsero mooe
Gawone woy
(1isweg) moy
pwoo Dow--ss moden
(iswoo)
(Lisaoo)) may.
Cissero)
Gtnseee sow
Sangee woy Covoee
Lts
way
10000
20000
40000
60003
50000
100000
Method, N
All experiments were run using a 32-core CPU machine with 128GB memory. SIM-MCMC and SIT-MCMC took approximately 0.5-1 second per iteration and were for 100 hours in the Cambridge data. SIM-NN took approxim- ately 0.1 seconds per iteration and was run for 50 and 100 hours in the Cambridge, and DC data, respectively. Finally, our framework (GENSIT) took on average 0.2 seconds per iteration and was run for 70 and 120 hours in the Cambridge, and DC data, respectively.
Figure SRMSE (y-axis) and CP (xx marker
size) across C (marker colour) computed cumu- latively along N iterations of Alg. for T C samples for GENSIT (Joint and Disjoint) and SIT- MCMC for the Cambridge dataset. The Joint GENSIT converges a lower SRMSE faster than SIT-MCMC while achieving better T coverage in all but the totally constrained C cases.
E Auxiliary Experimental Results
In this section we provide auxiliary experimental results for the Cambridge and Washington datasets.
E.1 Cambridge, UK
ODM's reconstruction sensitivity to different ensembles sizes E of Alg. I and number of steps N is examined in Fig. 9 A strong preference towards N >> E is evidenced by significant enhancements in both SRMSE and CP except for the Totally and Singly constrained ODMs (see Tab. 4. All other ODMs employ the GMB sampler(3.1] since the target distribution μ is intractable. GMB mixes slower than closed-form sampling and therefore needs to be run for larger number of steps N. Despite the fact that a larger ensemble size E facilitates exploration of the non-convex loss landscape, this not materialised in T-space. The SRMSE's invariance under different N, E' configurations detected in
22

---

the Totally and Singly ODMs suggests that the conditional independence of successive T samples compensates for potentially inadequately navigating the loss landscape. Besides, any type of T sampler used is invariant under initialisation. Therefore, leveraging both optimisation and sampling offers complementary benefits the estimation of both A and T.
The information captured in different L and D is portrayed in Fig. 10 Our choice of C(x, T, A D,v) in the Joint scheme (see 3rd x value from the left) achieves one of the lowest SRMSES and highest CPs across a, C. Tuning of y becomes essential when loss terms assimilating multiple data sources may yield conflicting metrics (lower SRMSE and low CP, and vice versa).
ORIGIN DESTINATION
TARGET M
SRMSE
(MUA,
T
199
CP (MAN)) T
SIT- MCMMC
IIMM
SIME MCMAC
MATRIX
D DISOIRT
JOINT
SN
DISTOINT
JOINT
MCCAC
NN
MCCMC
0.014
T
2.52
232
2.52
2.52
0.73
11%
4%
11%
21%
67%
79%
21%
Ti
2.43
23%
4T%
89%
constrained
A++
0.141
T
2.43
2.43
0.70
30%
31%
68%
learned
2.52
88%
T
2.52
2.52
14%
15%
0.38
57%
23%
89%
23%
T
0.58
0.58
0.69
82%
81%
68%
Singly constraincd
0.43
0.43
0.62
0.43
0.63
90%
90%
42%
93%
42%
A
T
0.43
0.43
0.61
93%
93%
72%
0.58
0.58
0.44
73%
72%
88%
learned
T
0.58
0.58
86%
87%
(a) ODMs with closed-form (tractable) T distributions
10)
ORIGIN- DESTINATION
TARGET M
SRMSE MM
99% CP (M (M) T")
GENSIT
GENSIT
DISIOINT
SIT-MCMC
DISJOIST
SIT-CMC
0.014
A
T
2.52
1.19
0.67
0.59
0.59
0.71
61%
4%
87%
37%
69%
21%
Doubly constrained
T...
0.141
A
T
2.30
1.15
0.68
0.55
0.59
1.22
18%
68%
89%
78%
20%
86%
learned
T
1.17
0.59
66%
86%
Doubly and 10% cell
0.014
A
2,52
0.89
0.42
0.71
476
49%
92%
22%
0.141
A
2.50
U.37
1.10
18%
79%
32%
0.43
0.56
92%
leamed
A
0.99
31%
T
1.08
0.43
70%
92%
T...
0.014
1.04
0.38
0.51
70%
92%
89%
20% cell constrained
Doubly and
T...
Ta Ap
0.161
A
2.50
0.92
0.38
1,06
18%
78%
314
leansed
A
T
2.52
1.04
0.37
0.96
72%
8%
47%
93%
This case not handled by the approach mentioned in the column.
(b) ODMs with intractable T distribution
where conditioning on problematic.
Table 6: Table expanded multiple SDE noise regimes.
Fig. sheds light on the effect of such information propagation on the convergence rate of running estimates of [P[- 'D] to the ground truth T". The Joint GeNSIT scheme converges to mean T estimate much earlier than SIT-MCMC across all regimes. Mean T estimates are improved the Joint GeNSIT compared to SIT-MCMC in confined T spaces (Doubly, Doubly and 10% cell, Doubly and 20% cell constrained ODMs). Under the same C regimes CP does not improve significantly as N grows large, which suggests that the variance of T samples appears stable as early as N 10a. In the Disjoint GeNSIT the information encoded in larger Cp is not propagated to the A updates. As result, SRMSE or CP improvements are detected in the course of N.
E.2 Washington, DC, USA
We append the expanded Tab. of results for the Washington dataset in Tab.
F
Reproducibility
Our codebase and the real-world data we used are accessible from the Supplementary Material. Trip and employment data for Cambridge, UK are obtained from the Office of National Statistics and used
23

---

Tue
Jontel, T2-
-- 0014
0.10, I-
-- Jaunol, T-
Ti
0024 Tx
I..
Junnel T..
learned T
0.01. T..
0014 T..T. Ro
0.0 T. T.,. To
lomed Ty T.
0.411 T-
0,044, T...
0.141 T...
0.014 T.,, T.
3.2
0.41, T... T..
loamed, T.....
0.04, Tu. T Tx
0.341 T.., T..
2.5
louned, T-..
1.4
24
1.24
2.
SEMSE
1.5-
1.0
1.1
SRMSE
0.5-
1.2-
0.5
8.6
0.4
0.2
Desgoott
Dlsiviut.
Dhijoin
Aoind
Disooint
Joint
Desgouat
Disgoou
Jonu
Degiviu.
C.S.A. 29
F.K.T.: 20
ZlA plag
4/ or om
<(F.A. mla
10,1000 30,200
100,100
500,20 1000,10 30002 100001
Method, (N, E)
cix
Cix.t.A.
C(k.T.A
Figure 9: Televel SRMSE and CP (x marker size) under different combinations of the num- ber of iterations N and ensemble size E used in Alg. for a fixed budget of x E = 10a. When N > E, local A, T information is greedily leveraged compared to global information accu- mulated over the NN ensemble of (!!) initial-
Loss operator
Figure 10: T-space SRMSE and 99% HMR cell CP (x marker size) achieved for different choices of L (see Tab. [5] for definitions). Both low SRMSE and a high CP cannot be achieved for these C choices. This pattern is more prevalent in L's that assimilate both and D.+. suggesting the possibility of contradictory information in the two datasets. Losses that include L (T, A) achieve lower SRMSEs and are not significantly reduced in light of more data.
isations. Closed-form T sampling eliminates the dependence between successive T samples. contrast, Gibbs Markov Basis sampling requires a larger number of steps N to converge to the stationary target distribution μ the expense of higher sample autocorrelation.
The entire feature space for the Washington DC data accessible through this repository All data assets leveraged are under the CC-BY 4.0 licence. Once the repository is cloned and a Python virtual environment needs to be created the gensit package can be installed. Depending on the machine capabilities, the number of workers and threads per worker can be set as follows:
47
export $N WORKERS kinput_herel export $N_THREADS [input_hers
F.1
Cambridge, UK
export
SCAMBRIDGE cambridge_worke wtwpwoooooootereeeeesseecoce to_msoas
The experimental protocols used in the paper can be executed as follows. Experiment can be reproduced by running

---

DATA
TARGET M
SRMSE
SSI
199 CP
FRAMEWORK
D
(MOEN) T
(MO:N) T
(M T
GMEL
Y.C
A
2.43 0.15
01.38 0.02
5% T%
2.47 + 0.00 (low)
0.50 + 0.00 (low)
7% 1% (low)
2.47 = 0.00 (high)
0.51 + 0.00 (high)
SIM-NN
A
24% 3% (high)
2.48 0.00
0.50 + 0.00
0% = 0% (leamed)
(learned)
(learned)
2.47 = 0.00 (low)
0.50 0.00 (low)
6% + 1% (low)
2.17 + 0.00 (high)
0.51 0.00 (high)
A
19% + 3% (high)
2.48
0.50 0.00 (learned)
y. C
3% = 1% (leamed)
GENSIT
(Disjoint)
(learned)
2.40 = 0.08 (low)
0.43 0.01 (low)
35% + 1% (low)
2.39 + 0.09 (high)
0.43 0.43 - 0.01 (high)
47% + 1% (high)
T
2.40 = 0.08
0.43 + 0.01
36% + 1% (learned)
(learned)
(learned)
2.45 = 0.00 (low)
0.50 = 0.00 (low)
+ 0% (low)
2.45 + 0.00 (high)
0.50 + 01.00 (high)
A
2% - 0% (high)
2.45 0.00
0.50
2% - 0% (leamed)
GENSIT (Joint)
(learned)
(learned)
0.08 (low)
0.01 (low)
44% + 1% (low)
2.37 0.08 (high)
0.45 (high)
44% 1% (high)
T
2.37 0.09
0.15 0.01
44% 1%
(learned)
(learned)
(learned)
Table Table expanded multiple SDE noise regimes.
gensit run http://wwppppaacaaaaaaaaaecaaata toml -sm \
-et
-et SIM_NN -et NonJointtaabeeSIM.N -nt TN_THREADS -nw $N_WORKERS
and
gensit run http:///ppptsccccca//ccccllll/ar .oml
-et JointTablesIM_NW \
$N_ THREADS -nw WORKERS
The exploration-exploitation experiment can be run through commands
gensit run http:////aanunttttttt/c////////e
sm -et SIM_MCMC -et SIM_NN -et NonJointTabbleSIM.n
$1 THREADS -nw $N_ WORKERS
and
gensit httpaa/eeeanttttaaaaaaaaaaaaa/er toml -sm -et JointTableSIM.NW -nt $N_THREADS -nw $N_WORKERS
Finally, the comparison of loss functions in experiment 3 can be performed using
run http://sszhppppttoortrt////////p -et JointTableSIM.NV -nt $N_ THREADS -nw $N_WORKERS
toml
25

---

F.2 Washington, DC, USA
The GENSIT Disjoint scheme run using
gensit run http://ppptpppoootts//////////ep wtw.ttiinnniiininnnoooommmo toml
-sm -et NonJointTableeIIM.M -nt $N_THREADS -nw $N_WORKERS
The GENSIT Joint scheme is using
gensit http////.ppptt////////////////p/ toml
-am -et JointTableSIM
-nt SN_EADSADS -nw $N WORKERS
The SIM-NN comparison is rum using
gensit /ttp:///iiiitttsoo/////////////r httprrrmmmnnnnmmmnsooootoo..... toml \
-sm -et SIM_NN -nt $N_THREADS $H_WORKERS
The GMEL comparison run using
gensit run http://wwppttteeoott///////////a toml \
-sm -et htapl/ttttttioeeeeoooooeeee_e_po -nt $N_THREADS -nw $N_WORKERS
For a detailed explanation of how reproduce the Tabs. and Figs. we refer the reader to the README. md file in root directory of our codebase, which can be found in the Supplementary Material.
26

---

1. Claims
Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Question:
Answer: [Yes]
Justification: Scalability claims are backed by computational complexities, two large empir- ical studies, and convergence studies. Performance claims backed by empirical studies and appropriate metrics.
2. Limitations
Question:
Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification: Limitations in terms of handling structured cell constraints discussed in concluding remarks. Limitations in terms of performance under low-constraint regimes discussed in main experimental section.
3. Theory Assumptions and Proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [NA]
Justification: Citation of proven results provided.
4. Experimental Result Reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main ex- perimental results the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Answer: [Yes]
Justification: All necessary information needed to reproduce experimental findings is provided in the main text and Appendix. Code to reproduce computational experiments, as well the real-world datasets provided in the supplementary material.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, described in supplemental material?
Answer: [Yes]
Justification: As mentioned previously, all code and data are provided in the supplementary material and will become open access at the end of the reviewing period.
6. Experimental Setting/Details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: Hyperparameter tuning and train/test/validattin splits are provided either in the main text or in the appendix.
7. Experiment Statistical Significance
Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Answer: [Yes]
Justification: We cite all error bars where computationally feasible and discuss the statistical significance our empirical claims in the main text.
8. Experiments Compute Resources
Question: For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed reproduce the experiments?
27

---

Answer: [Yes]
Justification: All computational complexities are cited the main paper, while compute resources are referenced in both the main text and Appendix.
9. Code Of Ethics
Question: Does the research conducted the paper conform, in every respect, with the NeurlPS Code of Ethics http://:///.oripp...ccccp///a//a
Answer: [Yes]
Justification: All datasets and frameworks used conform to the Code of Ethics. Spatial datasets are provided at regional levels where individual anonymity is preserved.
10. Broader Impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [Yes]
Justification: This work can assist practitioners working with multi-agent systems in the development of agent-based models. There are potential positive impacts of applying this work to agent-based modelling allowing decision-makers to test different policy scenariors.
Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification: No risks posed.
12. Licenses for existing assets
Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Question:
Answer: [Yes]
Justification: All data assets leveraged are under the CC-BY 4.0 license
13. New Assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [Yes]
All new assets are properly documented and released under the CC-BY 4.0 license.
Justification:
14. Crowdsourcing and Research with Human Subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?
Answer: [No]
Justification: All datasets used are publicly accessible and no primary data collection has been undertaken.
15. Institutional Review Board (IRB) Approvals or Equivalent for Research with Human Subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: No primary data collection has been undertaken.
28


Reviews:
[{"review_id": "1oUPACEoe3", "reviewer": "Reviewer_sr5H", "summary": "The paper introduces a new framework, GENSIT, for calculating origin-destination matrices for Agent-Based Models (ABMs) representing the movement of individual agents, from partially observed summary statistics. The method uses a neural differential equation as a physics-based way of embedding spatial interactions. This allows benefits over prior methods such as scaling linearly with the number of origin-destination pairs IJ. The authors demonstrate improved compute performance and accuracy over prior methods on two datasets of Cambridge, UK and Washington D.C, US.", "questions": "The algorithm is quite complex and there are many different parts to the algorithm described (Algorithm 1 in the paper). It is unclear which parts are the most important for the improved performance claimed by the paper, since there is no ablation study of the individual components of the algorithm.\n\nWould an ablation study be feasible? Or is there a reason that all the parts of the algorithm must be included together?", "limitations": "The paper mentions that the problem of more complex cell structure constraints for population synthesis remains unsolved.\n\nThere is a potential limitation that it may be hard to get accurate ground truth of people movement data in real life scenarios to validate methods such as this. This is not mentioned in the paper (as far as I could see).\n\nThe paper could mention the potential risks when simulating social behaviour that the behaviour is oversimplified and doesn't capture the nuances of real human behaviour, or that decisions made on the basis of agent-based models of social systems could negatively affect some individuals.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The paper provides a novel contribution to methods for estimating origin-destination matrices, using neural differential equations. It provides a conceptual motivation for this; incorporating a physics based model which matches the domain, and also shows experimental evidence of improvement over prior techniques.\n\nThe paper appears high quality. It is well written with no typos that I noticed. There are clear diagrams to illustrate the method, and there is comprehensive explanation of the mathematics, with equations provided in detail.\n\nThe paper is significant to reducing the computational cost of estimating origin-destination matrices for Agent-based models, which is a sub-problem of Agent-based modelling.", "weaknesses": "Despite being well written, the paper is often dense and a bit tricky to understand. It is heavy with mathematics and often uses technical terms without explaining what they mean in this context (e.g. _\"explore the multimodal matrix distribution over a discrete combinatorial support, and incurs discretisation errors\"_ in the paper abstract.) \n\nThe experimental results seem good, however the table they are presented in was complicated, and hard to understand at first glance what the improvement is over existing techniques, and what each of the columns of the table actually mean. It might be worth having a summary table which only shows the high level results in a clear way.\n\nThe authors could perhaps give more context of how this method fits into Agent-based Modelling overall, and how it could be used by ABM practitioners.\n\nThe paper is generally well referenced, however there is no mention or explanation of the connection to the wider field of ABM calibration (i.e. estimating unknown parameters for ABMs), which is broader than just Origin-destination matrix estimation. It would be good to add some mention of this broader context of ABM calibration, the challenges involved, and how this technique helps.  The paper claims to reduce the cost of estimating origin-destination matrices, however it's not clear whether that is a computationally limiting step for ABMs, is this a large fraction of the overall compute required for ABMs? What practical benefit is there to having more accurate and faster calculation of Origin-destination matrices?\n\nThe paper does not give much information on the implementation details of the neural network for the neural differential equation. There is a reference for PyTorch, but it is not stated explicitly whether PyTorch is used. There are no values for hyperparameters and the architecture of the neural network is not described. It is not specified which optimisation algorithm is used for training.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their feedback. \n> Despite being well written, the paper is often dense and ...\n\nPlease see global rebuttal (**Theme 2**), where we list improvements we now make to increase clarity throughout, improve notation, and provide a more gentle introduction to technical concepts. We welcome additional suggestions.\n> The experimental results seem good, however...\n\nWe refer the reviewer to the global rebuttal (**Theme 2**) where we address all reviewers' comments on improving the presentation of experimental results. \n> The authors could perhaps give more context of how this method fits into Agent-based Modelling overall ... The paper is generally well referenced, however...\n\nWe refer the reader to the global rebuttal (**Theme 1**) where we propose adding two paragraphs in the introduction of the main paper to clarify the motivation, link to ABMs and the impact of our framework. \n> The paper does not include much information on the implementation details...\n\nAlthough we include these details in our codebase and briefly described them in our App. (lines 434-445), we refer the reviewer to the global rebuttal (**Theme 2**) where we propose updating the App. to include an expanded description of the neural network implementation details including the architecture, optimisation algorithm, and hyper-parameters. We also improve the link to that information in the main paper and provide a diagram of the architecture in Fig.2 of the pdf in the global rebuttal.\n> The algorithm is quite complex and ... Would an ablation study be feasible?...\n\nWe argue that we are effectively performing a component-based ablation study through our comparison with SIM-MCMC [R2], SIT-MCMC [R3], SIM-NN [R1] and propose that we explicitly mention this in our experimental section. SIM-MCMC and SIT-MCMC both use a Markov Chain Monte Carlo (MCMC) sampler to obtain $\\boldsymbol{\\theta}$ estimates from the Harris-Wilson SDE. This corresponds to replacing lines 6-17 of Algorithm 1 with an MCMC sampler. Tables 1 and 2 clearly indicate that our method and SIM-NN [R1], which both leverage a Neural Network to learn $\\boldsymbol{\\theta}$ and a Euler-Maruyama solver to sample $\\mathbf{x}\\vert\\boldsymbol{\\theta}$, are superior to the MCMC-based methods in terms of reconstruction error and coverage of the ground truth ODM.  Moreover, neither SIM-NN nor SIM-MCMC sample discrete ODMs as GeNSIT does in lines 18-23 of Algorithm 1. Therefore, comparing against these methods allows us to test the impact that operating only in the continuous intensity $\\boldsymbol{\\Lambda}$ space has in the overall reconstruction error and coverage of the ground truth ODM. Tables 1 and 2 demonstrate a significant improvement in these two metrics when ODM estimation is performed in the discrete $\\mathbf{T}$ space. Therefore, we can conclude that the components of GeNSIT are superior to those found in the prior art because of the computational savings and improved $\\boldsymbol{\\theta}$ estimates brought by the Neural Network (lines 6-17 of Algorithm 1) and the enhanced reconstruction and coverage brought by discrete ODM sampling (lines 18-23 of Algorithm 1). \n\nApart from this comparison, we have now run additional synthetic experiments provided in the global rebuttal (**Theme 3**) where we perform a component-based ablation study of the computation time and SRMSE of each component along increasing ODM dimensions $(I,J)$ (Fig. 2) as well as a study of the discrete ODM sampling component for varying ODM dimensions (Fig. 3) and number of agents $M$ (Fig. 4). We believe these new synthetic experiments will complement the ablation study provided in the experimental section of the paper and propose including Fig 2 in the introduction of the experimental section and Figs. 3 and 4 in Section C of the App.\n> There is a potential limitation that...\n\nIndeed, such data would not be generally accessible. This is precisely why we need a robust framework for estimating origin-destination matrices when fully observed ground truth data is not available, such as GeNSIT. We have leveraged two datasets where ground truth data exists in order to be able to validate GeNSIT. Moreover, it is often the case when summary statistic data are readily available (e.g. agent population at a regional level) from population censuses. Our method can flexibly assimilate these types of data as noise-free constraints in the discrete ODM space. We propose identifying this limitation in our concluding remarks. \n> The paper could mention the potential risks...\n\nWe propose to re-emphasize the point on the potential negative impact of policy interventions elicited from ABM situations, and refer the reviewer to the global rebuttal (**Theme 1**) where we propose extending our introduction to discuss amongst other things the social impact of ABMs and GeNSIT's limitations.\n\n# References \n\n[R1]: Gaskin, T., Pavliotis, G. A., & Girolami, M. (2023). Neural Parameter Calibration for Large-Scale Multiagent Models. Proceedings of the National Academy of Sciences.\n\n[R2]: Ellam, L., Girolami, M., Pavliotis, G. A., & Wilson, A. (2018). Stochastic Modelling of Urban Structure. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences.\n\n[R3]: Zachos, I., Damoulas, T., & Girolami, M. (2023). Table Inference for Combinatorial Origin-Destination Choices in Agent-Based Population Synthesis. Stat.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "Gv0FORGvf1", "reviewer": "Reviewer_C3cp", "summary": "The paper considers the task of estimating discrete Origin-Destination Matrices (ODM) which will be useful for generating synthetic agent populations. Rather than the computationally inefficient approach of searching the huge space of such matrices using expensive Agent-based Simulations, the study takes the approach of generating continuous approximations of ODMs (intensity matrices) using a Spatial Interaction Model, whose parameters are estimated using a neural network. Then the continuous ODMs are used to estimate the discrete ODMs through Gibbs Markov Basis Sampling.", "questions": "1) What are the observations that we have, based on which we are calibrating the SIM using neural network? Is it only the attractiveness of each location? Can there be other versions of the problem where more observations are available?\n2) I understand that the task of estimating the discrete ODM is a combinatorial one, but can't we significantly reduce the space by imposing realistic constraints on the matrix structure?", "limitations": "NA", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "The main strengths of the paper are as follows:\n1) It links the existing literature on spatial interaction models and the Harris-Wilson system with agent-based mobility simulation\n2) The work provides a solid approach to estimating the OD intensity matrices from observed data using the SIM\n3) Further, it provides a new approach to solve the combinatorial problem of finding the discrete ODM, based on Gibbs Markov Basis Sampling\n4) They validate their approach with good experimental results", "weaknesses": "I would flag two main concerns:\n1) The way this paper has been written, its contributions seem to be too specific to the task of generating discrete OD matrices, which may be useful to a specific application (agent-based population synthesis). I think the general approach (especially estimating parameters of the SIM, sampling discrete matrices using Gibbs Markov Basis Sampling) etc have merit beyond this task, but somehow that is not coming through\n2) I did not find the paper very easy to read. I think a table with notations will be useful, also a bit of more clarity on Spatial Interaction Models will be useful. Also, the nature of the observations (y, D) is not very clear.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their feedback.\n> ..contributions seem to be too specific to the task of generating discrete OD matrices...\n\nWe point to the global rebuttal (**Theme 1**) where we elaborate on our framework's wider applicability, importance, connection to ABMs, and societal impact. \n> ...I think a table with notations will be useful...\n\nPlease refer to **Theme 2** in the global rebuttal where we include a notation Tab. \n\nTo improve the exposition of SIMs, we add the following before line 119 to link SIMs to the multinomial logit:\n\n\"We note that additional data at the origin, destination and origin-destination level can be assimilated into SIMs. This can be achieved by incorporating them as terms in the maximum entropy argument used to derive the $\\boldsymbol{\\Lambda}$ functional forms in equations (3), (4), and (5). We note that the SIM's $\\boldsymbol{\\Lambda}$ is equivalent to the multinomial logit [5], which generalises our $\\boldsymbol{\\Lambda}$ construction to accommodate for more data.\"\n\nWe explain this in the App.: \n\nThe set of observation data $\\mathcal{D}$ may include the observed log destination attraction $\\mathbf{y}\\in \\mathbb{R}^{J}$ and the total distance travelled from each origin $\\mathbf{D}\\_{\\cdot+} = (\\mathbf{T}^{*} \\odot \\mathbf{C}) \\in \\mathbb{R}\\_{>0}^{I}$. For Cambridge, both $\\mathbf{y}$, $\\mathbf{D}$ have been sourced from the UK's population census dataset provided by the Office of National Statistics. Their spatial resolution is regional: middle and lower super output areas for $\\mathbf{y}$, $\\mathcal{D}$, respectively. For DC, we only have access to a feature matrix $\\mathbf{Y} \\in \\mathbb{R}^{(I+J)\\times D}$ from which we extract a column to use as $\\mathbf{y}$. \n> What are the observations that we have, based on which we are calibrating the SIM using neural network?...\n\nWe leverage the observed log destination attraction at each destination location. Our sensitivity study in Fig. 6 p.14 also leverages the total distance travelled by agents from each origin location $\\mathbf{D}\\_{\\cdot+} = (\\mathbf{T}^{*} \\odot \\mathbf{C})$. Also, our Joint scheme's loss on $\\mathbf{T}\\vert \\boldsymbol{\\Lambda},\\mathcal{C}_T,\\mathcal{D}$ propagates the information from the summary statistic data $\\mathcal{C}_T$ to $\\boldsymbol{\\Lambda}$. The SDE solution is expressed in terms of the log origin/destination attraction we need observations on the solution to calibrate the SIM parameters using the Neural Net. \n\nAdditional observations at the origin, destination and origin-destination level can assimilated by incorporating them as terms in the maximum entropy argument used to derive the functional forms of $\\boldsymbol{\\Lambda}$ in equations (3), (4), and (5). The SIM intensity is equivalent to the multinomial logit [R1], which allows us to define an arbitrary utility model inside the $\\exp$ in the numerator of these equations. For example, one might want to include data $\\mathbf{o} \\in \\mathbb{R}^{I}$ in addition to the existing $\\mathbf{x}, \\mathbf{C}$ in the totally constrained SIM (total $=M$). Then, we need to maximise:\n$$-\\sum_{i,j}^{I,J} \\Lambda_{ij}(\\log(\\Lambda_{ij}) -o_i-\\alpha x_j+\\beta c_{ij}) - \\mu(\\sum_{ij}^{I,J} \\Lambda_{ij} - M),$$\nwhere $\\mu$ is the Lagrange multiplier. This yields $$\\Lambda_{ij} = \\frac{M(\\exp(o_i+\\alpha x_j-\\beta c_{ij}))}{\\sum_{i,j}^{I,J} \\exp(o_i+\\alpha x_j - \\beta c_{ij})}.$$\nWe propose adding a paragraph in a section in the App. as guidance to practitioners who wish to construct new utility functions.\n> ...can we significantly reduce the space by imposing realistic constraints on the matrix structure?\n\nWe are exploring the space of discrete ODMs $\\mathbf{T}$ with summary statistics $\\mathcal{C}\\_T$ known as the fiber $\\mathcal{T}\\_{\\mathcal{C}}$ using a finite Markov Basis (MB). Additional structural constraints on the ODM is equivalent to finding a subspace of $\\mathcal{T}\\_{\\mathcal{C}}$ where these structural constraints are satisfied. For this, we update our MB to explore this subspace without breaking the connectivity of the fiber graph (i.e. irreducibility).\n\nSparse matrices can be handled through the cell constraints $\\mathbf{T}\\_{\\mathcal{X}'}$ for $\\mathcal{X}'\\subseteq\\mathcal{X}$. Symmetric ODMs (e.g. an adjacency matrix of a bipartite graph) can be explored as follows. The subspace of $\\mathcal{T}\\_{\\mathcal{C}}$ consists of square matrices $I\\times I$ with equal row and column sums. We construct a MB on $\\mathcal{T}\\_{\\mathcal{C}}$ by applying the following modification to equation (14) (see **Note**):\n$$\\mathbf{f}_l(x) = \\begin{cases}\n    +\\eta & \\text{if } x \\in \\{(i_1,j_1),(i_2,j_2),(j_1,i_1),(j_2,i_2)\\} \\\\ \n    -\\eta  & \\text{if } x \\in \\{(i_1,j_2),(i_2,j_1),(j_2,i_1),(j_1,i_2)\\} \\\\\n    0 & \\text{otherwise.}\n\\end{cases}$$\nFirst, this guarantees that every $2\\times 2$ move in the original $\\mathcal{M}$ that modifies a cell along the diagonal is by definition symmetric as to the other cell it modifies, meaning if cell $(i,j)$ is modified then so is $(j,i)$. Second, without loss of generality assume that a move modifies a $2\\times2$ section in the upper triangular section of the ODM. This ensures that the same move is applied symmetrically to the lower triangular section of the ODM. In both cases, any move will guarantee that the ODM will be symmetric after removing any duplicate Markov Bases generated in the process.\n\nWe propose detailing the case for symmetric ODMs in the App. In general, encoding structural constraints on the ODM remains an open challenge, and propose stating this in our conclusion.\n\n**Note**: Equation (14) should read: \n$$\\mathbf{f}_l(x) = \\begin{cases}\n    +\\eta & \\text{if } x \\in \\{(i_1,j_1),(i_2,j_2)\\} \\\\\\\n    -\\eta  & \\text{if } x \\in \\{(i_1,j_2),(i_2,j_1)\\} \\\\\\\n    0 & \\text{otherwise.}\n\\end{cases}$$\n# References\n[R1]: Nijkamp, P., & Reggiani, A. (1988). Entropy, Spatial Interaction Models and Discrete Choice Analysis: Static and Dynamic Analogies. European Journal of Operational Research.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "jvLKFc6uUd", "reviewer": "Reviewer_qQL1", "summary": "The authors propose a generative model for origin-destination matrices $I \\times J$ with known marginals. The model assumes that the counts in the matrix are realizations of Poisson random variables conditioned on parameters determined by an unknown intensity $\\Lambda$. The matrix entries must also respect the marginal constraints of this unknown intensity. Additionally, they must satisfy a model parameterized by two parameters, $\\alpha$ and $\\beta$, which control the attractiveness of each destination and the costs between pairs of origin-destination locations.", "questions": "The real data sets are relatively small (69x13 and 179x179). Can you run large synthetic examples to study the scalability of GeNSIT?", "limitations": "They addressed adequately the limitations.", "rating": 5, "confidence": 1, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "GeNSIT introduces an efficient framework for jointly sampling both discrete combinatorial space of agent trips (T) and its continuous mean-field limit ($\\Lambda$). This dual approach is a novel contribution. Another strength is that the framework scales linearly with the number of origin-destination pairs (I+J). GeNSIT regularizes the parameter space of neural networks by embedding physics models. This integration of physics models regularizes the model adding reliability and robustness to the framework.", "weaknesses": "I found this article very difficult to understand. One reason is that it seems to have been written for experts in origin-destination matrices, which is a very small group of people at this conference. A more important reason is that the paper is unnecessarily theoretical, overusing mathematical notation and jargon, and leaving several points unclear, which hinders comprehension. Despite the authors' rigor in definitions and notation, there appear to be errors or inaccuracies in several parts that left me very confused.\n\nTo be more specific, consider the constraints discussed in (3), (4), and (5). These pertain to the marginals $\\Lambda_{++}, \\Lambda_{i+}, \\lambda_{+j}$  of the intensity $Lambda$. Since this intensity is unknown, these models cannot be trained. Later on, the authors seem to state that they will replace the marginals based on $Lambda$ with the sufficient statistics $T_{++}, T_{i+}, T_{+j}$. If this is the case, why isn't the problem presented directly in terms of these statistics? Or will the marginals of the function $\\Lambda$ to be learned satisfy, for example, $\\Lambda_{i+} = T_{i+}$?\n\nThe authors aim for precision and rigor, but the text falls short and is confusing. Here are some examples:\n\n- The attractiveness of destination $z_j$ is governed by a stochastic differential equation (SDE)? But your algorithm outputs a single $ Z = (z_1, ..., z_J) $ array. Is this single array an asymptotic result of the SDE?\n- In equation (7), make it explicit that this is the distribution conditioned on $\\Lambda_{ij}$ and that you assume the counts $ T_{ij} $ are independent conditionally. This is mentioned in passing a few lines below.\n- Why is $\\mathcal{C}_{\\Lambda} \\subset \\mathcal{C}_T$? We have that $\\mathcal{C}_{\\Lambda}$ is the set of distributions with a constraint on $\\Lambda$ and  $\\mathcal{C}_T$ is the set of tables realized with fixed margins. How can one be contained within the other if the sets pertain to different elements (distributions indexed by intensity in one case, and count tables in the other)?\n- In practice, how can $\\mathcal{C}_{\\Lambda}$ be known to be used as input in the algorithm?\n- What happens if the inputs $\\mathcal{C}_{\\Lambda}$ are completely inconsistent with the inputs $\\mathcal{C}_T$? For example, $C_{++} = 20000$ and $\\Lambda_{++} = 100$?\n- In equation (9), what is $ k $? Did I miss the definition of this notation? I keep going back through the text to check if I missed a definition.\n- \"Upon conditioning on sufficient summary statistics $\\mathcal{C}_T$, any choice of probability model in (7) becomes equivalent [1].\" I don't understand what \"equivalent\" means here. Do you mean that, conditioned on the sufficient statistics, the data becomes independent of the parameters and there is no way to choose among the different realizations ${T}$ that meet the constraints?\n- \"The hard-coded constraints $\\mathcal{C}_T$ are treated as noise-free data on the discrete table space.\" How can they be free of random noise? They are sums of Poisson distributions and therefore have a Poisson distribution (for example, $T_{i+} \\sim \\text{Poisson}(\\Lambda_{+i})$).\n- \"The aforementioned random variables are summary statistics uniquely determined upon conditioning on $T$.\" This doesn’t make sense. Conditioning on the entire matrix $T_{ij}$, the marginals become degenerate random variables with a single value with probability 1.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their feedback.\n> I found this...\n\nPlease see **Themes 1,3** in global rebuttal.\n> To be more specific, consider...\n\nThe SIM's intensities in (3),(4) are embedded in the Harris-Wilson SDE. The goal is to solve the SDE's inverse problem by learning the parameters $\\boldsymbol{\\theta} = (\\alpha,\\beta)$ appearing in (3),(4). The intensity $\\boldsymbol{\\Lambda}$ marginals will indeed be fixed to satisfy $\\Lambda_{++} = \\mathbb{E}[T\\_{++}\\vert \\mathcal{C}\\_T]$ and $\\boldsymbol{\\Lambda}\\_{\\cdot+} = \\mathbb{E}[\\mathbf{T}\\_{\\cdot+}\\vert \\mathcal{C}\\_T]$ when $T\\_{++},\\mathbf{T}\\_{\\cdot+} \\in \\mathcal{C}\\_T$. So, $\\mathbf{T}$-level constraints will match $\\boldsymbol{\\Lambda}$-level constraints except for the case in (5). This is because the doubly constrained SIM is non-identifiable as argued in line 114. \n> The attractiveness of destination $z_j$ is...\n\nThe $\\mathbf{z}$ is the numerical solution of the SDE after $\\tau$ steps, which we aim to explicitly state it in the exposition of SIMs. \n> In equation (7), make it explicit that...\n\nWe replace equation (7) with $T_{ij} \\vert \\Lambda_{ij} \\sim \\text{Poisson}(\\Lambda_{ij})$ and extend line 131 to \"... number of agents travelling from $i$ to $j$, and $T_{ij} \\perp T_{i'j'} \\vert \\Lambda_{ij},\\Lambda_{i'j'} \\;\\forall\\; i\\neq i', j\\neq j'$\".\n> Why is $\\mathcal{C}_{\\Lambda} \\subseteq \\mathcal{C}_T$?...\n\nWe replace $\\mathcal{C}\\_{\\Lambda} \\subseteq \\mathcal{C}\\_T$ with the explanation that \"every summary statistic constraint in $\\Lambda$-space is also applied in $\\mathbf{T}$-space, i.e. we always set $\\Lambda\\_{++} = \\mathbb{E}[T\\_{++}\\vert \\mathcal{C}\\_T]$ and $\\Lambda_{i+} = \\mathbb{E}[T\\_{i+}\\vert \\mathcal{C}\\_T]$\".\n> ...how can $\\mathcal{C}_{\\Lambda}$ be...\n\nIn practise, we choose $\\mathcal{C}\\_{\\Lambda}$ between $\\{\\Lambda\\_{++}\\}$ and $\\{ \\Lambda\\_{++},\\boldsymbol{\\Lambda}\\_{\\cdot+} \\}$. When $\\mathcal{C}\\_T = \\{T\\_{++},\\mathbf{T}\\_{\\cdot+}\\}$, we set $\\mathcal{C}\\_{\\Lambda} = \\{\\Lambda\\_{++},\\boldsymbol{\\Lambda}\\_{\\cdot+}\\}$. In all other cases we set $\\mathcal{C}\\_{\\Lambda}= \\{\\Lambda\\_{++}\\}$. When $\\mathcal{C}\\_{T}$ contains at least $T\\_{++},\\mathbf{T}_{\\cdot+},\\mathbf{T}\\_{+\\cdot}$ and $\\mathcal{C}\\_{\\Lambda} = \\{\\Lambda\\_{++},\\boldsymbol{\\Lambda}\\_{\\cdot+}\\}$, the dependence between $\\mathbf{T}\\vert \\boldsymbol{\\Lambda},\\mathcal{C}$ and $\\mathbf{y}\\vert \\mathbf{x}, \\mathcal{C}$. As a result, constraints are implicitly weighted (hard $\\mathcal{C}\\_{\\Lambda}$ and soft $\\mathcal{C}\\_T$ constraints), which inflicts identifiability issues in $\\boldsymbol{\\Lambda}$.\n> What happens if the inputs $\\mathcal{C}\\_{\\Lambda}$ are...\n\nWe always fix $\\Lambda_{++} = \\mathbb{E}[T\\_{++}\\vert \\mathcal{C}\\_T]$ and $\\Lambda\\_{i+} = \\mathbb{E}[T\\_{i+}\\vert \\mathcal{C}\\_T]$, so we would not encounter such a case. Despite this, the $\\boldsymbol{\\Lambda}$ information is propagated to and from $\\mathbf{T}$ via $\\frac{\\boldsymbol{\\Lambda}}{\\Lambda\\_{++}}$ in $\\mathbf{T}\\vert \\boldsymbol{\\Lambda},\\mathcal{C}_T$.\n> In equation (9), what is $k$?...\n\nMany thanks for picking up our omission. Index $k$ refers to the $k$-th basis operator applied over the cell subset $\\mathcal{X}\\_k$. We will explicitly state this in line 173. \n> Upon conditioning on sufficient summary statistics $\\mathcal{C}\\_T$...\n\nWe recognise that our writing needs improvement here. The Poisson model for $T\\_{ij}\\vert \\Lambda\\_{ij}$ is a modelling choice of the practitioner. Another potential choice would have been $T_{ij} \\sim \\text{Binomial}(\\Lambda_{++},\\frac{\\Lambda\\_{ij}}{\\Lambda\\_{++}})$. By conditioning on the summary statistics $T\\_{++},T\\_{i+},T\\_{+j}$ these two modelling choices would both yield the same conditional distributions $\\mathbf{T}\\vert \\boldsymbol{\\Lambda},\\mathcal{C}\\_T$. In that sense the two modelling choices are equivalent under conditioning on $\\mathcal{C}\\_T$. We will clarify this in the main text.\n> \"The hard-coded constraints $\\mathcal{C}_T$ are treated as...\n\nWe apologise for the confusion here. $\\mathcal{C}\\_T$ are realisations of the Poisson random variables $T\\_{++}\\vert\\boldsymbol{\\Lambda},T\\_{i+}\\vert\\boldsymbol{\\Lambda},T\\_{+j}\\vert\\boldsymbol{\\Lambda}$. By conditioning on $\\mathcal{C}\\_T$, these variables are no longer random. We propose explicitly stating this in line 135.\n> \"The aforementioned random variables...\n\nWe propose replacing our sentence with \"The aforementioned summary statistics become Dirac random variables upon conditioning on $\\mathbf{T}$.\"\n> The real data sets are relatively small...\n\nWe argue that modelling between 900 (69x13) to 32k (179x179) cells and approximately 40k-100k agents are sufficiently large datasets given that they are used to model entire cities such as Cambridge, UK and DC, USA.\n\nWe also point the reviewer to the global rebuttal (**Theme 3**), where we introduce new experimental results on GeNSIT's scalability and propose to include them in the main and Appendix. \n> This is a theoretical paper with no immediate connection to societal impact\n\nWe'd like to respectfully argue against this statement and provide evidence to the contrary. This is not just a theoretical paper: it introduces methodologies and algorithms that are applied to real-world datasets and ABMs of up to 32,000 cell dimensions and 100,000 agents. We also refer the reviewer to the global rebuttal (**Theme 1**) where we propose elaborating on the benefits of ODM estimation to ABM calibration as well as the social impact of our work.\n\n # References\n\n[R1]: Ferguson, N. et. al (2006). Strategies for mitigating an influenza pandemic. Nature.\n\n[R2]: Ferguson, Neil et. al (2020). Report 9: Impact of non-pharmaceutical interventions (NPIs) to reduce COVID19 mortality and healthcare demand. Imperial College London.\n\n[R3]: W Axhausen, Kay, Andreas Horni, and Kai Nagel (2016). The multi-agent transport simulation MATSim. Ubiquity Press.\n\n[R4]: Crooks, A., et. al (2021). Agent-Based Modeling and the City: A Gallery of Applications. Urban Informatics.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "hnkBt0sX8Q", "reviewer": "Reviewer_dVmA", "summary": "This paper introduces a novel framework named Generating Neural Spatial Interaction Tables (GENSIT) for efficiently generating origin-destination (OD) matrices in neural spatial interaction models. The primary objective is to address the challenges of existing methods, such as continuous approximations and ad-hoc discretizations, by operating directly on the discrete combinatorial space of OD matrices. The proposed method utilizes neural differential equations to embed spatial interactions and learns the agents' trip intensity, thereby providing a more accurate and computationally efficient solution. The framework is validated through large-scale spatial mobility agent-based models (ABMs) on multiple real scenarios.", "questions": "- Could you provide more details on why generating OD matrices efficiently is crucial in practice? Specifically, how does this efficiency impact the usability and effectiveness of ABMs in policy-making and other applications?\n- The paper mentions that GENSIT outperforms prior art in terms of reconstruction error and computational cost. It would be helpful to include a detailed comparison with these methods, highlighting the specific scenarios where GENSIT provides the most significant advantages.", "limitations": "NA", "rating": 5, "confidence": 5, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- The paper presents a novel approach to generating OD matrices by directly operating in the discrete combinatorial space, which addresses the limitations of continuous approximations. The integration of neural differential equations for embedding spatial interactions is a significant advancement in this domain.\n- The framework is thoroughly evaluated using real-world datasets from Cambridge, UK, and Washington, DC, USA. The experimental results demonstrate that GENSIT outperforms existing methods in terms of reconstruction error and ground truth matrix coverage, while also being computationally efficient.\n- The paper is well-structured and clearly presents the motivation, methodology, and results. Figures and tables are effectively used to illustrate the performance improvements of the proposed method.", "weaknesses": "- The paper does not sufficiently explain the practical importance of generating OD matrices efficiently. For readers unfamiliar with this line of literature, it would be beneficial to provide more context on why efficiency in this process matters and how it impacts the overall utility of ABMs in real-world applications.\n- While the methodology is appealing, some parts of the explanation, particularly the notations used in the paper, could be simplified for better understanding. Including more intuitive explanations or examples might help in making the content more accessible.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their comments and questions. These will improve the paper's clarity and accessibility to a wider audience.\n\n> The paper does not sufficiently explain the practical importance of generating OD matrices efficiently. For readers unfamiliar with this line of literature, it would be beneficial to provide more context on why efficiency in this process matters and how it impacts the overall utility of ABMs in real-world applications.\n\nWe appreciate that the paper would significantly benefit from a more elaborate explanation of the importance of generating ODMs efficiently and the impact this has on ABM simulation and calibration. We refer the reviewer to the global rebuttal (**Theme 1**) where we propose elaborating on exactly this.  \n\n> While the methodology is appealing, some parts of the explanation, particularly the notations used in the paper, could be simplified for better understanding. Including more intuitive explanations or examples might help in making the content more accessible.\n\nThank you. We refer the reviewer to the global rebuttal (**Theme 2**) where we suggest clarity modifications relevant to the reviewer's comment. We also offer additional experimental results and analysis in the global rebuttal (**Theme 3**) to make the content more accessible. If the reviewer has additional quantities in mind in need of further intuition please let us know. \n\n> Could you provide more details on why generating OD matrices efficiently is crucial in practice? Specifically, how does this efficiency impact the usability and effectiveness of ABMs in policy-making and other applications?\n\nWe refer the reviewer to the global rebuttal (**Theme 1**) where we propose elaborating on exactly this in our introduction.\n\n> The paper mentions that GENSIT outperforms prior art in terms of reconstruction error and computational cost. It would be helpful to include a detailed comparison with these methods, highlighting the specific scenarios where GENSIT provides the most significant advantages.\n\nWe now better highlight the specific scenarios of most significant advantage of GenSIT by explicitly describing these settings in the start of the experimental section (see global rebuttal **Theme 2**), in the analysis of the experiments and also in the concluding remarks where we also highlight settings where GenSIT advantages are limited (totally constrained ODM). The computational cost comparison is detailed by additional experiments listed in the global rebuttal (**Theme 3**) as well as a new section in the Appendix listing all methods' computational complexities (see rebuttal to reviewer C8eo).", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "Yrndl97beD", "reviewer": "Reviewer_h7J4", "summary": "The authors proposed a novel framework that can effectively generate the discrete origin-destination matrices (ODMs) in agent-based models. By using neural differential equations to embed the spatial interactions and operating directly on the discrete combinatorial space, the method overcomes the limitations of traditional continuous models. The superiority of the proposed method is demonstrated in experiments over large-scale real-world datasets.", "questions": "1.\tIn this work, an ensemble of neural networks is used to predict the vector that governs the utility parameters. I would like to know what networks are used in this part and how different networks affect the reconstruction error.", "limitations": "The major limitations are discussed in the experiment section, and no significant negative social impact can be observed in this paper.", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1.\tThis paper contains many mathematical terms and notations. The authors introduce these notations by giving examples and explanations of their practical meaning, such as the description of the Harris-Wilson system in Eq.6, which improves the quality of the presentation.\n2.\tIn the experiment section, the authors conducted experiments on large-scale datasets that were obtained from real-world cases (Cambridge, UK, and Washington DC, USA). This significantly improves the soundness of this paper.", "weaknesses": "1.\tThe author introduced two sampling schemes in computing the loss for the neural network: the joint and disjoint schemes. This is also described in lines 13 to 15 in Alg.1. However, the description of these two functions is not adequate in the main text; only the difference in the information of T is mentioned. This could be improved if a more detailed explanation of the practical meaning of these two schemes and how they are related to the loss functions is added.\n2.\tThe author claims that the proposed method scales linearly with the number of origin-destination pairs. This major contribution should be mentioned in the experiment section in the main text. Right now, although real-world datasets are used, the experiment part is not easy to understand, and some of the results are not clearly highlighted. The experiment section could be improved if the author gives more results with different scales of data.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their helpful comments and important questions that will significantly improve the paper's clarity. \n\n> The author introduced two sampling schemes...\n\nWe propose adding the following to the Appendix:\n\n\"The Disjoint scheme consists only of loss terms that depend directly on fully observed data $\\mathcal{D}$ (either log-destination attraction $\\\\mathbf{y}$ or total distance agents travelled by origin location $\\mathbf{D}\\_{\\cdot+}$). In contrast, the Joint scheme consists of loss terms that either depend on the same fully observed data and on the partially observed table $\\mathbf{T}$ through the table marginal $\\mathbf{T}\\vert \\boldsymbol{\\theta}, \\mathbf{x},\\mathcal{C},\\mathcal{D}$. \nThe Joint scheme is an instance of a Gibbs sampler on the full posterior marginals $\\boldsymbol{\\theta}\\vert (\\mathbf{x},\\mathbf{T},\\mathcal{C},\\mathcal{D})$, $\\mathbf{x}\\vert (\\boldsymbol{\\theta}, \\mathbf{T},\\mathcal{C},\\mathcal{D})$ and $\\mathbf{T}\\vert (\\boldsymbol{\\theta}, \\mathbf{x}, \\mathcal{C}, \\mathcal{D})$. The Disjoint scheme is an instance of a collapsed Gibbs sampler where we sample from $\\boldsymbol{\\theta}\\vert (\\mathbf{x},\\mathcal{C},\\mathcal{D})$, $\\mathbf{x}\\vert (\\boldsymbol{\\theta}, \\mathcal{C},\\mathcal{D})$ and then from $\\mathbf{T}\\vert (\\boldsymbol{\\theta}, \\mathbf{x}, \\mathcal{C}, \\mathcal{D})$. This means we integrate out $\\mathbf{T}$ by means of $p(\\boldsymbol{\\theta}\\vert \\mathbf{x},\\mathcal{C},\\mathcal{D}) = \\sum_{\\mathbf{T}} p(\\boldsymbol{\\theta}\\vert \\mathbf{T},\\mathbf{x},\\mathcal{C},\\mathcal{D})P(\\mathbf{T}\\vert\\mathbf{x},\\mathcal{C},\\mathcal{D})$ and $p(\\mathbf{x}\\vert \\boldsymbol{\\theta},\\mathcal{C},\\mathcal{D}) = \\sum_{\\mathbf{T}} p(\\mathbf{x}\\vert \\mathbf{T},\\boldsymbol{\\theta},\\mathcal{C},\\mathcal{D})P(\\mathbf{T}\\vert\\boldsymbol{\\theta},\\mathcal{C},\\mathcal{D})$. \nTherefore, we use the Joint scheme when we have reason to believe that the covariance between $\\boldsymbol{\\theta},\\mathbf{x}$ and $\\mathbf{T}$ is small. This would be the case when the agent trip intensity is influenced by both the Harris-Wilson SDE and the realised number of agent trips. In contrast, we use the Disjoint scheme to accelerate convergence by reducing the covariance between $\\boldsymbol{\\theta},\\mathbf{x}$ and $\\mathbf{T}$. This would be the case when the agent trip intensity is governed only by the Harris-Wilson SDE and not by the realised number of agent trips.\"\n\nWe propose expanding the explanation in line 157:\n\"The Joint scheme corresponds to a Gibbs sampler on the full posterior marginals $\\boldsymbol{\\theta}\\vert (\\mathbf{x},\\mathbf{T},\\mathcal{C},\\mathcal{D})$, $\\mathbf{x}\\vert (\\boldsymbol{\\theta}, \\mathbf{T},\\mathcal{C},\\mathcal{D})$ and $\\mathbf{T}\\vert (\\boldsymbol{\\theta}, \\mathbf{x}, \\mathcal{C}, \\mathcal{D})$. The Disjoint scheme corresponds to a collapsed Gibbs sampler where we sample from $\\boldsymbol{\\theta}\\vert (\\mathbf{x},\\mathcal{C},\\mathcal{D})$, $\\mathbf{x}\\vert (\\boldsymbol{\\theta}, \\mathcal{C},\\mathcal{D})$ and then from $\\mathbf{T}\\vert (\\boldsymbol{\\theta}, \\mathbf{x}, \\mathcal{C}, \\mathcal{D})$ by integrating out $\\mathbf{T}$.\"\n\nTo further improve the clarity, we propose adding a column to Table 4 in the Appendix clearly stating whether each loss operator is used in the Joint or Disjoint scheme. \n> The author claims that the proposed method scales linearly...\n\nWe appreciate that a lot of results have been presented in the experimental section of the paper. We refer the reviewer to the global rebuttal where we propose ways to improve the clarity of the experimental section (**Theme 2**) as well as include new results on the scalability of GeNSIT across ODM dimensions and number of agents (**Theme 3**).We hope that these Figures will complement the experimental section and highlight that our method scales linearly with the number of origin-destination pairs $IJ$.\n> In this work, an ensemble of neural networks...\n\nAlthough we include information about the Neural Network (NN) we used in our attached codebase and in lines 434-445 of the Appendix, we recognise that this could be made clearer and more detailed. We refer the reviewer to the global rebuttal (**Theme 2**) where we propose adding two paragraphs in Section C of the Appendix to detail the types of networks and hyper-parameters used in our experiments. A further clarification will be added to the Appendix:\n\nWe leverage a multi-layer perceptron with one hidden layer to build a map from the observed log destination attraction $\\mathbf{y}$ to the Harris-Wilson SDE parameters $\\boldsymbol{\\theta}$. The SDE is encoded in the NN as follows. The Euler-Maruyama solver provides forward solutions from $\\boldsymbol{\\theta}$ to $\\mathbf{x}$ at different time steps. The loss function computes the discrepancy between the forward solution $\\mathbf{x}$ after $\\tau$ steps and the observed data $\\mathbf{y}$ at the stationary equilibrium of the SDE, allowing us to encode the physics directly into the NN. A more complicated architecture with a larger number of weights would potentially lead to overfitting the SDE parameters $\\boldsymbol{\\theta}$ to $\\mathbf{y}$ compromising our framework's generalisability.\n> The major limitations are discussed in the experiment section, and no significant negative social impact can be observed in this paper.\n\nWe refer the reader to the global rebuttal (**Theme 2**), where this is addressed.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "PfkXuRoR01", "reviewer": "Reviewer_C8eo", "summary": "The paper introduces the GENSIT framework to generate ODMs in agent-based models. The framework addresses challenges in traditional methods that rely on continuous approximations and ad-hoc discretizations, proposing a more efficient approach that operates directly on the discrete combinatorial space. The method leverages neural differential equations to model spatial interactions and demonstrates superior performance in terms of reconstruction error and computational efficiency. The proposed method is illustrated using real-world data from Cambridge, UK, and Washington, DC, USA.", "questions": "1. What is the model assumption? Does the GENSIT framework perform well in scenarios where the spatial distribution of the data is highly irregular or non-uniform?\n\n2. Is this method robust to outliers?\n\n3. This method seems very complicated to me. Can you summarize the intuition behind it in a short paragraph?", "limitations": "See comments above.", "rating": 6, "confidence": 1, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "1. Combining neural differential equations, the GENSIT framework introduces a novel approach to generating ODMs by directly operating on the discrete space.\n\n2. The framework's effectiveness is demonstrated through empirical validation on serval real-world datasets, showing superior performance in terms of reconstruction error and computational efficiency.\n\n3. The proposed method is computationally efficient and suitable for large-scale datasets.", "weaknesses": "1. No theoretical guarantee. I understand that the theoretical aspects of GENSIT may be challenging to analyze. However, I hope the author could provide some theorems on consistency and robustness, which will help the reader better understand this work. \n\n2. Computational cost should be further discussed. The computational complexity is given by $O(NEJ(\\tau|W|+I))$, what is the scale of each parameter?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their feedback.\n> No theoretical guarantee...\n\nTheoretically analysing the full joint framework is challenging future work, however we offer some existing theoretical insight that we will briefly discuss in the main and expand in the App.\n\nWe leverage the universal approximation theorem for analysis the SIM parameter learning using a Neural Net. We also leverage results from the literature on discrete contingency tables complement Proposition 3.1 for the analysis of discrete ODM sampling. \n\nFirst, the fundamental theorem of Markov Bases [R1] establishes an equivalence between a Markov Basis (MB) and the generator of an ideal of a polynomial ring. By virtue of the Hilbert basis theorem [R3], any ideal in a polynomial ring has a finite generating set. Therefore, there exists a _finite_ MB for the class of inference problems we are dealing with. This MB connects the fiber, i.e. the support of all discrete ODMs satisfying summary statistics in the form of row, column sums and/or cell constraints. By Proposition 3.1, we can construct a Gibbs MB sampler that will converge to Fisher's non-central hypergeometric distribution on the fiber in finite time. \n\nMoreover, we state an important result on the convergence rate of the MB sampler [1,2] in the App.: A MB sampler converges to its stationary distribution in at most $M^2$ steps ($M$: number of agents). Although this result is not directly applicable to us in its given form ($\\eta=\\pm 1$ on Fisher's central hypergeometric distribution), we conjecture that it can be extended to $\\vert\\eta\\vert>1$ and Fisher's _non-central_ hypergeometric distribution as evidenced by our experimental results. A MB connects tables (elements) of the fiber graph (discrete state-space of tables). The fastest convergence of an MB sampler is achieved by a MB that has the smallest possible diameter (longest path between two elements of the fiber graph).\n> Computational cost should be further discussed...\n\nWe point to the global rebuttal (**Theme 3**) where we propose adding a Fig. comparing GeNSIT's computation time with increasing ODM dimensions. We also include this section in the App.: \n\n$N$: number of iterations, $I,J$: number of origins,destinations.\n- GeNSIT: $\\mathcal{O}(NE(\\tau J\\vert+IJ))$, where $\\tau$ is the number of time steps in the Euler-Maruyama solver, $E$ is the ensemble size. We set $\\tau=1$, and $E=1$. \n- SIM-MCMC [2]: $\\mathcal{O}(NJ(LI + J^2))$ (low $\\sigma$ regime) and $\\mathcal{O}(NIJLKn_pn_t)$ (high $\\sigma$ regime), where $L$ is the number of leapfrog steps in Hamiltonian Monte Carlo, $1<K<N$ is the number of stopping times, and $n_p$,$n_t$ are the number of particles and temperatures used in Annealed Importance Sampling, respectively. Typical ranges are $n_p \\in [10,100]$, $n_t \\in [10,30]$, and $L\\in[1,10]$ and $E\\in[1,10]$.\n- SIT-MCMC [3]: $\\mathcal{O}(NJ(LI + J^2))$ (low noise regime) $\\mathcal{O}(NIJLKn_pn_t)$ (high noise regime). See 1,2 for details.\n- SIM-NN [1]: $\\mathcal{O}(NE(\\tau J))$. See 1. for details.\n- GMEL [4]: $\\mathcal{O}(n_nn_l((I+J)D + IJ))$ where $\\mathbf{Y} \\in \\mathbb{R}^{(I+J) \\times D}$ is the feature matrix for all locations, $D$ is the number of features per location, $n_n < D$ is the number of nodes per layer, and $n_l$ is the number of hidden layers. In the DC dataset, we have $I+J=179$ since every origin is also a destination (we do not count them twice). Typical ranges include $n_l \\in [1,10]$ and $n_n \\in [1,D]$.\n> What is the model assumption?...\n\nWe briefly mention the following in the main and expand in the App.:\nA fundamental assumption is that the continuous agent trip intensity $\\boldsymbol{\\Lambda}$ is governed by a SIM. Therefore, agents are discouraged from travelling long distances due to high travelling costs incurred and encouraged to travel to locations that offer high utility, such as areas with high job availability. This means that highly irregular or non-uniform trip data might not be captured well by $\\boldsymbol{\\Lambda}$. However, operating at the discrete ODM level allows us to assimilate trip data as constraints. This shrinks the ODM support and enhances spatial imputation of missing data much better than doing so at the continuous $\\boldsymbol{\\Lambda}$ level. In the limit of trip data constraints, we are guaranteed to improve our estimated of missing trip data regardless of their spatial distribution and presence of outliers. This is theoretically founded [R1] and empirically shown in Tabs. 1,2 of our experimental section. In poorly constrained settings, we rely more on $\\boldsymbol{\\Lambda}$ for spatial imputation, which increases reconstruction error. Outlier presence will be reflected in the row/column sum constraints facilitating table sampling in that region of the ODM support. If no row/column sum constraints are imposed, exploring such regions of the support would be significantly hindered.\n> Can you summarize the intuition behind the method ...?\n\nGeNSIT estimates ODMs summarising agent trip counts between origin and destination locations. It learns how probable it is for agents to travel to a specific destination captured in the continuous trip intensity. This probability is governed by the trade-off between the cost and utility gained from travelling. Conditioned on the intensity, we sample discrete trip counts (ODMs) that satisfy the observed trip summary statistics, such as the agent population at each origin and/or agent demand for a particular destination. Therefore, our framework jointly learns the continuous intensity and discrete ODM until convergence.\n\n# References\n\n[R1]: Diaconis, P., & Sturmfels, B. (1998). Algebraic Algorithms for Sampling from Conditional Distributions. The Annals of Statistics.\n\n[R2]: Windisch, T. (2016). Rapid Mixing and Markov Bases. SIAM Journal on Discrete Mathematics.\n\n[R3]: Simpson, S. G. (1988). Ordinal numbers and the Hilbert basis theorem. The Journal of Symbolic Logic.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:
We thank all reviewers. All agree on the paper's contributions: "novel contribution" [sr5H,qQL1,dVmA,h7J4,C8eo], "new approach" [C3cp]. The "paper appears high quality ... is significant" [sr5H], "provides a solid approach" [C3cp], is "well-structured and clearly presents the motivation" [dVmA]. An "efficient" [qQL1,dVmA,C8eo] framework adds "reliability", "robustness" [qQL1] and "soundness" [h7J4] with a "thorough evaluation" [dVmA] and "good experimental results" [C3cp], which demonstrate "superior performance" [C8eo], "outperforming existing methods" [dVmA]. 

We now address feedback on motivation (**Theme 1**), clarity (**Theme 2**), and scalability (**Theme 3**).
# Theme 1: Motivation (sr5H, dVmA, C3cp, qQL1)

We add at our intro:
> "Our framework has merit beyond ODM sampling in ABMs. The challenge of learning discrete contingency tables constrained by their summary statistics extends to other fields. Contingency tables have been widely studied in multiple instance learning [27, 11, 39] and ecological inference [31, 30, 32]. In Neuroscience one estimates the efficiency, cost and resilience (equivalent to $T_{ij}$) of neural pathways between pairs of brain regions $(i,j)$ to understand communication and information processing [R2,R3]. Epidemiology also investigates social contact matrices quantifying the number of contacts ($T_{ij}$) between types of individuals $(i,j)$ stratified by demographic characteristics, such as age [R4]."

We append to line 20: 
> "This is achieved by computationally expensive forward simulations, which hinders ABM parameter calibration and large-scale testing of multiple policy scenarios [R6]." 

We expand line 35: 
> "It is also necessary for population synthesis in ABMs [14], which is performed prior to simulation in order to reduce the size of the ABM's parameter space.  Moreover, it avoids ...". 

We append to line 61: 
> "This approximation effectively acts as a cheap ABM surrogate or emulator, facilitating faster forward simulations to be run by practitioners [R5]. This has tanglible benefits to ABM calibration allowing faster exploration of the parameter space. Our enhanced ODM reconstruction error demonstrates GeNSIT's ability to sufficiently approximate the ABM simulator at a fraction of the computational cost.".  

We discuss limitations and social impact in conclusions:
> "Our work also relies on the SIM's assumptions about the agents' decision-making process, which in practise is unobserved. An examination of different agent utility models could benefit the applicability of our framework. In terms of our work's social impact, policy decisions made from ABMs of social systems could negatively affect individuals, necessitating expert review and ethics oversight.".
# Theme 2: Clarity (sr5H, qQL1, C3cp, dVmA)
We now provide a Table of Notation, attached pdf (Tab. 1). 

App. Section C now gives full details on the NN. The architecture in Fig. 1 in pdf and lines 434-445 read:
> "Our NN is a multi-layer perceptron with one hidden layer, implemented in PyTorch [R1]. The input layer is set to the observed log-destination attractions $\mathbf{y}\in \mathbb{R}^{J}$ since we are learning the $\boldsymbol{\theta}$ that generates the observed physics $\mathbf{y}$. The output layer is two-dimensional due to the parameter vector $\boldsymbol{\theta}\in\mathbb{R}^2$. For both datasets we set the number of hidden layers to one and number of nodes to 20. The hidden/output layer has a linear/absolute activation function. 
>
> We use the Adam optimizer with $0.002$ learning rate. Bias is initialised uniformly at $[0,4]$. We follow [1,2] in fixing $\sigma_d=0.03$ and $\sigma_{T},\sigma_{\Lambda}$ to $0.07$ to reflect a $1\%$ and $3\%$ noise levels, i.e. $\sigma / \log(J) \approx 3\%$. We assume $\mathbf{y}$ are observations from the SDE's stationary distribution, hence our batch size is one. We initialise the Euler-Maruyama solver at $\mathbf{y}$ and run for $\tau=1$ and step size $\Delta t=0.001$. At equilibrium the only change in the log-destination attractions is attributed to the SDE's diffusion." 

In Tabs. 1 and 2 we rename columns to SRMSE$(\mathbf{M}^{(1:N)},\mathbf{T}^{\*})$, CP$(\mathbf{M}^{(1:N)},\mathbf{T}^{\*})$, and SSI$(\mathbf{M}^{(1:N)},\mathbf{T}^{\*})$ and link to their exact formulation in the App.

For Tab.  1, we keep the best SRMSE, CP for each method and ODM across all noise regimes, and move the rest to App., Section B. 

For Tab.  2, we keep the high-noise results and move the rest to App., Section B. We explain the second column (FEATs. $\mathcal{D}$) in the caption.
# Theme 3: Scalability (qQL1, h7J4, C8eo, dVmA)
We now offer synthetic experiments (attached pdf) for varying ODM dimensions $(I,J)$ and number of agents $M$, comparing computation time and ground truth reconstruction error (SRMSE). 

Our complexity is linear in $IJ$. We rederive the complexity fully in terms of $J$ and without $\vert\mathbf{W}\vert$. Our complexity is $\mathcal{O}(NE(\tau J + IJ))$ and we include Fig.2 (pdf), in the experimental section showing how our method's table sampling, intensity learning and total computation times scale with the ODM dimensions. Figs. 3-4  showing the variation of SRMSE are discussed in the main and App.

# References 
[R1]: Paszke, A. et.al. (2019). Pytorch: An imperative style, high-performance deep learning library. NeurIPS.

[R2]: Seguin, C., et.al. (2023). Brain Network Communication: Concepts, Models and Applications. Nature Reviews Neuroscience.

[R3]: Goñi, J., et. al (2014). Resting-Brain Functional Connectivity Predicted by Analytic Measures of Network Communication. PNAS.

[R4]: Mistry, D., et. al (2021). Inferring High-Resolution Human Mixing Patterns for Disease Modeling. Nature Communications.

[R5]: Banks, D. L., & Hooten, M. B. (2021). Statistical Challenges in Agent-Based Modeling. The American Statistician.

[R6]: Lempert, R. (2002). Agent-based modeling as organizational and public policy simulators. PNAS.

Meta-review: The paper discusses the estimation of so-called origin-destination matrices in agent-based models. 

The paper is generally well received by the reviewers and is written to a high standard, providing a lot of mathematical details. Reviewers also positively commented about the application to real data.  However, the method described is relevant only to a relatively small subarea of research. 

A minor issue is that the referencing could generally be improved. As some reviewers have pointed out, the connection to the field of calibrating ABMs is not really discussed. Similarly the field of Bayesian calibration of computer models, there is a broad literature discussing with posteriors of the form of eq (8), see e.g. [1 on page 4, 2, equation (3), 3 just to give some examples].

[1] Ridgway, J. (2017). Probably approximate Bayesian computation: nonasymptotic convergence of ABC under misspecification. arXiv preprint arXiv:1707.05987. 

[2] Schmon, S. M., Cannon, P. W., & Knoblauch, J. (2020). Generalized posteriors in approximate Bayesian computation. arXiv preprint arXiv:2011.08644.

[3] Park, M., Jitkrittum, W., & Sejdinovic, D. (2016, May). K2-ABC: Approximate Bayesian computation with kernel embeddings. In Artificial intelligence and statistics (pp. 398-407). PMLR.
Decision: Accept (poster)