Generating Origin-Destination Matrices in Neural Spatial Interaction Models

Authors: ioannis  zachos
mark  girolami
theodoros  damoulas
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 12:53:39

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.

Main paper: 
 introduction: 
 High-resolution complex simulators such as agent-based models (ABMs) are increasingly deployed to assist policymaking in transportation [10, 20] , social sciences [3, 8, 14, 33] , and epidemiology [16, 21] . 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 [24] . Considering ABMs for the COVID-19 pandemic [16] as an example, the continuous mean field process corresponds to 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 [10] .The formal object of interest that describes the discrete count of these spatial interactions, e.g. agent trips between locations, is the origin-destination matrix (ODM). It is an I × J (two-way) contingency table T with elements T i,j ∈ N counting the interactions of two spatial categorical variables i, j ∈ N >0 , see Fig. 1 . 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 at 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 a continuous representation Λ coupled with the Harris-Wilson SDE (SIM rectangle). In the latter, the ground truth can be reconstructed by sampling in the discrete combinatorial space of constrained ODMs conditioned on Λ (GENSIT rectangle). ABM simulations scale with O(M log(M )) compared to GENSIT which scales with O(IJ), where M ≫ I + J is the size of the agent interaction graph. 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 as 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-hoc discretisation required when working with continuous approximations of the underlying discrete ODM T * .Traditional ABMs, Fig. 1 (Left), simulate individual-level and spatially granular discrete ODMs at a high computational cost [2] , which scales at least with O(M log(M )), where M ≫ I + J is the size of the agent interaction graph. These ODMs are then aggregated using a sum pooling operation to regional ones, T (s) , whose summary statistics C correspond to observed data. However, the lower-dimensional subspace of contingency tables T satisfying constraints C (e.g. row and column marginals), denoted by T C , is known to be combinatorially large [11] and therefore sampling and optimisation in that space is notoriously hard since it requires enumerating all elements of T C . This underlying challenge is the reason why prior art [34, 13, 40, 31, 26, 17] has been imposing a continuous relaxation of the ODM estimation problem, where the continuous approximation of the discrete contingency table restricts inference at the agent trip intensity level Λ. This leads to quantisations and inefficient rejection sampling, see Fig. 2 . Such intensity-level Spatial Interaction Models (SIMs) [44, 30] , Fig. 1 (Right), are derived from entropy maximisation arguments, Sec. 2, with summary statistics constraints. These models are embedded in the Harris-Wilson (HW) system of differential equations [19] 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 facilitates a mean-field ABM approximation. This approximation effectively acts as a cheap ABM surrogate or emulator, facilitating 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 a fraction of the computational cost. However, if 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 [44] 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.Competing approaches that operate directly on the discrete ODM space and are motivated by econometric arguments are discrete choice models [43] . However, these cannot encode summary statistic constraints without introducing large rejection rates. The work of [9] 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 [13] and the associated discrete ODM table over its entire support [46] . Such routines incur a computational overhead in the order of at least I or J due to the intractability of the Harris-Wilson 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 [17] by using a 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 origin-destination pairs IJ, which is at least I or J times faster than MCMC [13, 46] . We offer faster ODM reconstruction and enhanced uncertainty quantification compared to continuous approaches in seminal works [13, 17] and hybrid (discrete and continuous) approaches [46] in terms of the number of iterations required. It is the first framework that jointly explores the constrained continuous and discrete combinatorial ODM spaces in 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 statistics extends to other fields. Contingency tables have been widely studied in multiple instance learning [32, 12, 47] and ecological inference [36, 35, 38] . 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 [37] , [18] . Epidemiology also investigates social contact matrices quantifying the number of contacts (T ij ) between types of individuals (i, j) stratified by demographics, such as age [28] .The target distributions µ T | Λ, C T we sample from in Alg. 1 are listed below in Tab. 4.We note that the quantity ω ij = Λij Λ++ Λi+Λ+j 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 β = 0 in (3), (4). This translates to the travel cost having no effect on the agents' choice of destination. Fisher's non-central multivariate hypergeometric reduces to its central version if and only if ω ij = 1 [1] . The central version corresponds to a uniform distribution over T C . Additionally, the normalising constant of Fisher's non-central multivariate hypergeometric is a partition function involving a sum over all elements of T C satisfying the conditioned margins. By virtue of the extension to Chu-Vandermonde's theorem [5] proved in [46] , this normalising constant can be computed in constant time usingEQUATION) ODM CT µ T | Λ, CT DISTRIBUTION TRACTABLE Unconstrained ∅ I,J i,j exp(-Λ ij )Λ T ij ij T ij ! Poisson ✓ Totally Constrained T++ I,J i,j T ++ ! T ij ! Λ ij Λ ++ T ij Multinomial ✓ Singly Constrained T•+ I,J i,j T i+ ! T ij ! Λ ij Λ i+ T ij Product Multinomial ✓ Doubly Constrained T•+, T+• I i=1 T i+ ! J j=1 T +j ! T ++ ! I,J i,j=1 T ij ! I,J i,j=1 Λ ij Λ ++ Λ i+ Λ +j T ij Fisher's non-central multivariate hypergeometric × Doubly Constrained with 10% cells fixed T•+, T+•, TX 1 -"- -"- × Doubly Constrained with 20% cells fixed T•+, T+•, TX 2 -"- -"- ×Table 4 : List of target T distributions with their associated constraints C T .Theoretically analysing the full joint framework is challenging future work, however we offer some additional theoretical insight. We encourage the reader to 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 a polynomial ring. By virtue of the Hilbert basis theorem [39] , 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.Moreover, an MB sampler converges to its stationary distribution in at most A 2 steps, where A is the total number of agents [45] . Although this result is not directly applicable to our framework in its given form, we conjecture that it can be extended to η = ±1 and Fisher's central hypergeometric 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.e. the longest path between two elements of the fiber graph.Incorporating additional structural constraints remains an open problem. The case for sparse matrices can be handled through the cell constraints T X ′ for X ′ ⊆ X . Symmetric ODMs (e.g. an adjacency matrix of a bipartite graph) can be explored as follows. The subspace of T C consists of square matrices I × I with equal row and column sums. We construct an MB on T C by applying the following modification to equation ( 14):f l (x) =      +η if x ∈ {(i 1 , j 1 ), (i 2 , j 2 ), (j 1 , i 1 ), (j 2 , i 2 )} -η if x ∈ {(i 1 , j 2 ), (i 2 , j 1 ), (j 2 , i 1 ), (j 1 , i 2 )} 0 otherwise.First, this guarantees that every 2 × 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, j) is modified then so is (j, i). Second, without loss of generality assume that a move modifies a 2 × 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.The Disjoint scheme consists only of loss terms that depend directly on fully observed data D (either log-destination attraction y or total distance agents travelled by origin location D •+ ). In contrast, the Joint scheme consists of loss terms that either depend on 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 a Gibbs sampler on the full posterior marginals θ|(x, T, C, D), . Therefore, we use the Joint scheme when we have reason to believe that the covariance between θ, x and 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 θ, x and 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.In this section, we compare our method's computational complexity to competitive approaches and list the validation metrics employed for empirical comparisons in the main paper.Let N be the number of iterations, and I, J be the number of origins, destinations.1. GENSIT: O(N E(τ J + IJ)), where τ is the number of time steps in the Euler-Maruyama solver, E is the ensemble size. We set τ = 1, and E = 1.2. SIM-MCMC [13] : O(N J(LI + J 2 )) (low σ regime) and O(N IJLKn p n t ) (high σ 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 ∈ [10, 100], n t ∈ [10, 30] , and L ∈ [1, 10] and E ∈ [1, 10]. 
spatial interaction intensities and contingency tables: 
 Consider A agents travelling from I residences (origins) to J workplaces (destinations). The expected number of trips (intensity) between origin i and destination j is Λ ij and is unobserved. The average number of agents starting (ending) their journey from each origin i (to each destination j) is:EQUATIONThe expected total number of agents travelling is assumed conserved:EQUATIONThe family of models for intensities Λ that assimilate any collection of the above constraints are called Spatial Interaction Models [44] . The demand for each destination depends on its attractiveness denoted by z := (z 1 , . . . , z J ) ∈ R J >0 . In our example, this is the number of jobs available at each destination. Let the log-attraction be x := log(z). Between two destinations of similar attractiveness, agents are assumed to prefer nearby destinations. Therefore, a cost matrix C = (c i,j ) I,J i,j=1 is introduced to reflect travel impedance. These assumptions are justified by economic arguments [30] and establish the basis for the agents' utility function. The maximum entropy distribution of agent trips subject to Λ ++ = A yields a totally constrained SIM intensity:EQUATIONwhere C Λ = {Λ ++ } is the set of summary statistic constraints on Λ. Henceforth, we set Λ ++ ∈ C Λ unless otherwise stated. The vector θ = (α, β) contains the agents' two utility parameters controlling the effects of attractiveness and deterrence on the expected number of trips Λ. If α 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 Λ •+ , then we obtain the following singly (also known as production) constrained SIM intensity:EQUATIONwhere C Λ is expanded to include Λ •+ in this case. Moreover, settingEQUATIONwhere O(i), D(j) are called balancing factors that ensure that C Λ are satisfied. The balancing factors introduce I + J unknown parameters, rendering the intensity model unidentifiable. Alternatively, these factors are approximately recursively using iterative proportional fitting [44] , which is sensitive to initialisation. Including individual cell constraints in C Λ breaks the continuity of Λ 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.1.1 for more information on SIMs as a modelling choice.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 Λ functional forms in equations ( 3), (4), and (5). We note that the SIM's Λ is equivalent to the multinomial logit [29] , which generalises our Λ 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 families of SIMs is governed by the Harris-Wilson system of J coupled stochastic differential equations (SDEs) [19, 13] :EQUATIONwhere ϵ > 0 is a responsiveness parameter, κ > 0 is the number of agents competing for one job, δ ≥ 0 is a parameter related to the smallest number of jobs a destination can have, σ > 0 is the standard deviation of SDE's noise, and B t is a J-dimensional Wiener process. The term Λ +jκz j + δ in (6) reflects the net job capacity at 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 a 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 Λ 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:EQUATIONT is the I × J discrete contingency table summarising the number of agents travelling from i to j, and T ij ⊥ T i ′ j ′ |Λ ij , Λ i ′ j ′ ∀ i ̸ = i ′ , j ̸ = j ′ . 
neural calibration of spatial interaction models: 
 We now introduce our framework1 (GENSIT). We estimate the parameters of the continuous SIM intensity by employing an ensemble of Neural Networks ψ N N : R J → R 2 . This allows us to bypass Figure 3 : GENSIT: (a) successive iterations of Alg. 1 for a given ensemble member, (b) plate diagram for every iteration, ensemble member. We propose two schemes: a Joint and a Disjoint (see App. B.3.1 for details). Contrary to the latter, the former passes table T information to the loss L (see in (b)). We perform an optimisation step in the intensity Λ space and a sampling step in T space, with associated complexities O(τ J + IJ) and O(IJ). The Λ arises by the well-known family of SIMs (3),(4) coupled with the HW-SDE (6). The T sampling step generates discrete C T -constrained ODMs contrary to [13, 17] , which only operate on the continuous mean-field level Λ. the computational challenges of solving the HW-SDE inverse problem within a Bayesian framework [13, 46] . Conditioned on those estimates and for every ensemble member e = 1, . . . , E, we solve the HW-SDE to get estimates of the time-evolved log destination attraction x after τ time steps using an Euler-Maruyama numerical solver ϕ HW : R J → R J [27] . This allows us to incorporate the HW physics model into our parameter estimation without sampling from the SDE's intractable steady-state Boltzmann-Gibbs distribution, which was the case in [13, 46] . Instead of a log-destination attraction data (y) likelihood, we compute a loss operator L(• ; D, ν) that can assimilate data D from multiple sources on any transformation of T, Λ, x. We note that ν is used to denote the loss hyperparameters (see Tab. 5). Then, the NN parameters (weights and biases) W are updated using back-propagation using a suite of optimisation algorithms [22] . This step requires derivatives of the loss with respect to the NN parameters ∇ W L(• ; D, ν) to be computed, which is achieved using off-the-shelf auto-differentiation libraries [25] . This gradient is informed by x estimates and therefore by the dynamics of the HW-SDE in (6) . Hence, the resulting SIM intensity is both stochastic and physics-driven. These steps are depicted in Fig. 3b by following the arrows from right to left and in steps 6 to 12 of Alg. 1. We note that our framework is divided into two sampling schemes: a Joint and a Disjoint. The former passes T information to the loss operator L, whereas the latter does not. The Joint scheme corresponds to a Gibbs sampler on the full posterior marginals θ|(x, T, C, D), x|(θ, T, C, D) and T|(θ, x, C, D). The Disjoint scheme corresponds to a collapsed Gibbs sampler where we sample from θ|(x, C, D), x|(θ, C, D) and then from T|(θ, x, C, D) by integrating out T (see App. B.3.1).The described optimisation routine yields a point estimate for Λ for any given ensemble member e and iteration n. We either increase the ensemble size E to obtain a distribution over Λ, or we can treat Λ as a deterministic mapping of realisations of random variables x, θ whose generalised posterior [7, 23] is:EQUATIONwhere ν are loss-related hyperparameters. Each loss evaluation can be treated as a sample from a generalised likelihood (first term), while samples from the x prior (second term) are obtained by forward solving ϕ HW . Finally, we can enforce a prior over θ by appropriately initialising the NN parameters W.The continuous agent trip intensity Λ (continuous ODM) samples need to be mapped to discrete agent trips T (discrete ODM), which are required for agent population synthesis. We proceed by introducing the necessary machinery to achieve this. We wish to sample from the target measure µ evaluated as µ(T|Λ, C).Let the table cells be indexed byX = {(i, j) : 1 ≤ i ≤ I, 1 ≤ j ≤ J} such that T (x) = T ij is the table value of cell x = (i, j) ∈ 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 k : X k → {0, 1} I+J is defined to beEQUATIONHence, we can define summary statistic operators S k : T → N I+J as linear combinations of the basis operator, that isS k (T) = x∈X k T(x)1 k (x).A collection of such summary statistic operators is abbreviated by S(T).Algorithm 1 : Generating Neural Spatial Interaction Tables. O(N E(τ J + IJ))1: Inputs: evidence: C, C, D, funcs.: L, ψ N N , ϕ HW , M, µ, hyperparams.: N , E, τ , ν, κ, δ, σ, ϵ. 2: Outputs: x (1:N ) , θ (1:N ) , Λ (1:N ) , T (1:N ) . 3:for each ensemble member e = 1, . . . , E do 4:Initialise W (0) , T (0) .for each iteration n = 1, . . . , N do 6:θ (n) ← ψ N N (y; W (n-1)), with y ∈ D. ▷ Forward solving Neural Net.7:EQUATIONx(n) t ← x (n-1) t ∀ t = 1, . . . , τ . 9:for each time step t = 1, . . . , τ -1 do 10:x(n) t+1 ← ϕ HW (x (n) t , θ (n) ).▷ Solving Harris-Wilson SDE.11:x (n) ← x (n) τ .12: 3) or (4). ▷ Computing SIM intensity.Λ (n) ← Λ C (x (n) , θ (n) , C) using (13:if Joint sampling scheme is used then 14: 1) , Λ (n) ; D, ν).L (n) ← L(x (n) , T (n-else 16:L (n) ← L(x (n) ; {y}, ν). 17: Compute ∇ (n) W L. 18:Update W (n) using back-propagation. ▷ Updating Neural Net Weights with SGD.19:if µ is not tractable then 20:Sample f l uniformly at random from M.Find supp{η} such that T (n-1) + ηf l ≥ 0.22:T (n) ∼ µ(T (n-1) + ηf l |Λ (n) , C T ).else 24:T (n) ∼ µ(• |Λ (n) , C T ).Computing Neural Net loss.Closed-form sampling.Constraints on T can be expressed as fixed summary statistics C T = {s 1 , . . . , s K }, where each s k is a fixed evaluation of S k (T) with respect to the basis operator. For example, constraint {T •+ , T +• } can be expressed in terms of the basis operator (9) over the entire cell set X . Definition 3.1. Let C T = {s 1 , . . . , s K } be a set of constraints based on summary statistics operators S.A table T is C T -admissible if and only if S(T) = {s k } K k=1 ∈ C T . The subspace of T containing all C T -admissible I × J contingency tables is T C = {T ∈ T : S(T) = C T }.This space contains all discrete ODMs consistent with the aggregate summary statistics C T . Our goal is to efficiently sample from a measure µ on T C subject to arbitrary C T .The target distribution is tractable, i.e. the entire table can be 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 whereC T tractable ⊆ P({T ++ , T •+ , {T X l |X l ⊆ X , l ∈ N}}) ∪ P({T ++ , T +• , {T X l |X l ⊆ X , l ∈ N}}).Note that the unconstrained case of C T = ∅ is handled by the construction in (7) . This facilitates sampling T (1:N ) in closed-form and in parallel as shown in step 22 of Alg. 1. The most notable cases of tractable distributions are (see App. B.2.1):EQUATIONEQUATIONEQUATIONEach of the distributions above can assimilate cell constraints of the formC T = {T X l |X l ⊆ X , l ∈ N} without violating µ's tractability by limiting the support of T|Λ, C. If the constraint set C T contains at least both marginals, that is C T intractable ⊆ P({T •+ , T +• , {T X l |X l ⊆ X , l ∈ N}}) \ C T tractable, then the target distribution becomes Fisher's non-central multivariate hypergeometric [1, 42] :EQUATION)Direct sampling without rejection from this distribution for arbitrary I, J is infeasible [1] .Therefore, we devise an MCMC proposal on T C . Using a suite of greedy deterministic Algorithms [6] , we can initialise our MCMC with a T (0) ∈ T C . By virtue of definition 3.2 we guarantee that no proposed moves modify the summary statistics in C k (Condition 1) and that there exists a path between any two tables such that any table member of the path is C k -admissible (Condition 2). The collection C of K constraints generates K sets of Markov bases M 1 , . . . , M K . Our proposal mechanism consists of the universe M = K k=1 M k . See App. B.2.2 for info on Markov Bases. Definition 3.2. : A Markov basis M k is a set of moves f 1 , . . . , f L : X → Z satisfying:EQUATIONA Gibbs Markov Basis (GMB) sampler can now be constructed (see steps 20 -22 in Alg. 1). Proposition 3.1. (Adapted from [11] ): Let µ be a probability measure on T C . Given a Markov basis M that satisfies 3.2, generate a Markov chain in T C by sampling l uniformly at random from {1, . . . , L}. Let η ∈ Z. If the chain is at T ∈ T C , determine supp (η) such that T + ηf l ≥ 0. ChooseP(η) ∝ x∈X :f l (x)̸ =0 (µ(T (x) + ηf l (x))) -1and move to T ′ = T + ηf l for the choice of η. An aperiodic, reversible, connected Markov chain in T C is constructed with stationary distribution proportional to µ(T). Proof is found in [11] .We note that when individual cells are fixed (T X l ∈ C), the summary statistics operator is applied over X l 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 state space, leading to better mixing times [41] .: O(N J(LI + J 2 )) (low noise regime) O(N IJLKn p n t ) (high noise regime). See items 1,2 for details.4. SIM-NN [17] : O(N E(τ J)). See item 1 for details.5. GMEL [26] : O(n n n l ((I + J)D + IJ)) where Y ∈ R (I+J)×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 ∈ [1, 10] andn n ∈ [1, D].We leverage the following three validations metrics employed throughout the literature [26, 46, 34, 40, 31] . The Standardised Root Mean Square Error (SRMSE) between the target ODM M (table T or intensity Λ) samples and the ground truth T * is equal toSRMSE(M (1:N ) , T * ) = I,J i,j=1 E[M (1:N ) ij |C, D] -T * ij 2 IJ   IJ i,j=1 E[M (1:N ) ij |C, D] IJ   -1, and the Sorensen Similarity Index (SSI) is equal toSSI(M (1:N ) , T * ) = 1 IJ I,J i,j=1 2 min E[M (1:N ) ij |C, D], T * ij E[M (1:N ) ij |C, D] + T * ij ,whose domain is [0, 1] and values closer to 1 imply a better ODM fit. Finally, the coverage probability of the ground truth table T * is calculated by first identifying the lower L q T(1:N ) ij and upperU q T (1:N ) ijboundaries of the high probability region containing q% of the total mass for each table cell x = (i, j) ∈ X . Then, the q% cell coverage probability is equal toCP q (M (1:N ) , T * ) = 1 IJ IJ i,j=1 1 L q M (1:N ) ij |C, D ≤ T * ij ≤ U q M (1:N ) ij |C, D .In this section, we detail the implementation details and experimental protocols used to obtain the results in the main paper.Regarding the HW-SDE, an Euler-Maruyama numerical solver ϕ HW is employed throughout the paper with a time discretisation step of ∆t = 0.01 and number of steps τ = 1. The low and high SDE noise levels correspond to σ = 0.014 and σ = 0.141, respectively, to establish a comparison against [13, 17, 46] . Following [13, 17, 46] , we set the responsiveness parameter ϵ = 1, and the parameter δ 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 Λ+++δJ The set of observation data D may include the observed log destination attraction y ∈ R J and the total distance travelled from each origin D •+ = (T * ⊙ C) ∈ R I >0 . For Cambridge, both y, 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 y, D •+ , respectively. For DC, we only have access to a feature matrix Y ∈ R (I+J)×D from which we extract a column to use as y.Our NN is a multi-layer perceptron with one hidden layer, implemented in PyTorch [25] and depicted in Fig. 7 . The input layer is set to the observed log-destination attractions y ∈ R J since we are learning the θ that generates the observed physics y. The output layer is two-dimensional due to the parameter vector θ ∈ R 2 . For both datasets we set the number of hidden layers to one and number of nodes to 20. The hidden, output layers have a linear and absolute activation functions, respectively. The latter is important in guaranteeing that our θ estimates are always positive. The NN parameters W are initialised by sampling uniformly over the region [0, 4] (J+1)×20+21×2 . We use the Adam optimizer [22] with 0.002 learning rate. Bias is initialised uniformly at [0, 4]. We follow [13, 46] in fixing σ d = 0.03 and σ T , σ Λ to 0.07 to Figures based on synthetic experiments are produced as follows. We generate ground truth tables (100 × 100), (200 × 200), (300 × 300), (400 × 400), (500 × 500), (600 × 600), (700 × 700), (800 × 800), (900 × 900), (1000 × 1000) as explained in the main paper. For each of these tables, Algorithm 1 and SIM-NN are run for N = 10 We follow the setup of [13, 46] to initialise SIM-MCMC , SIT-MCMC and ensure comparability with our models. The computational budget in Tab. 1 and Figs. 10,8 is set to N = 10 5 . SIM model parameters are fixed to ϵ = 1, κ = 1.025, δ = 0.0128 and σ d = 3% × log(J). For x updates the acceptance is monitored to be at least 90%, while the θ-update'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 normalizing constant, the number of particles is set to 100 and a uniform temperature schedule of 50 inverse temperatures is employed. The percentage of positive signs is maintained around 75% and above.In the case of the Washington DC data, we employ the same train/test/validation test split as in [26] . In terms of hyperparameter optimisation for GENSIT , we leverage the same architecture as in the Cambridge data and optimise the learning rate on the validation set. To ensure a meaningful comparison, we also optimise the learning rate and multitask weights of the GMEL framework on the validation set. We run each method appearing in Tab. 2 for an ensemble of size E = 10 obtain error bars of this ensemble.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 approximately 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.In this section we provide auxiliary experimental results for the Cambridge and Washington datasets.ODM's reconstruction sensitivity to different ensembles sizes E of Alg. 1 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 a larger number of steps N . Despite the fact that a larger ensemble size E facilitates exploration of the non-convex loss landscape, this is not materialised in T-space. The SRMSE's invariance under different N, E configurations detected in 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 in the estimation of both Λ and T.The information captured in different L and D is portrayed in Fig. 10 . Our choice of L(x, T, Λ ; D, ν) in the Joint scheme (see 3rd x value from the left) achieves one of the lowest SRMSEs and highest CPs across σ, C. Tuning of ν becomes essential when loss terms assimilating multiple data sources may yield conflicting metrics (lower SRMSE and low CP, and vice versa). -: This case is not handled by the approach mentioned in the column.MATRIX C σ TARGET M ↓ SRMSE (M (1:N ) , T * ) ↑ 99% CP (M (1:N ) , T * ) GENSIT [SIT- MCMC] [SIM- NN] [SIM- MCMC] GENSIT [SIT- MCMC] [SIM- NN] [SIM- MCMC] DISJOINT JOINT DISJOINT JOINT Λ 2.(b) ODMs with intractable T distribution (13) , where conditioning Λ on C is problematic. Our codebase and the real-world data we used are accessible from the Supplementary Material. For a detailed explanation of how to 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.Codebase found at https://github.com/YannisZa/GeNSIT 
experimental results: 
 We empirically test our framework on both synthetic and real-world data from Cambridge, UK and Washington, DC, USA. We compare GENSIT against SIM-MCMC [13], SIM-NN [17] , SIT-MCMC [46] and the Geo-contextual Multitask Embedding Learner (GMEL) [26] . See App. C.1 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 Λ 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,E,F for experimental protocols, implementation details, a comprehensive sensitivity study on our framework, and reproducibility using our package gensit. We begin by offering synthetic experiments for varying ODM dimensions (I × J) and number of agents A, comparing computation time and ground truth reconstruction error (SRMSE) in Figs. 4 and 5 . For a fixed (I × J) dimension we sample an intensity uniformly at random and subsequently generate a ground truth table T * with fixed A = T ++ by sampling from the Multinomial distribution in (10) . For Figs. 4 and 5b we set T ++ = 10 6 while for Fig. 5a we vary this for 0.01, 0.1, 0.25, 0.5 and 1 × 10 6 . Our framework achieves a speed-up of O(J 2 ) compared to previous methods such as SIM-MCMC, SIT-MCMC. Reconstruction errors grow linearly in dimension (I × J) and exponentially in number of agents A.C T = {T •+ , T +• , T X 50% }. Constraint T XWe now assess our framework's capacity to estimate agent home-to-work trips in real-world ODM prediction problems in Cambridge, UK [46] and Washington DC, USA [26] .In the Cambridge dataset, the ground truth ODM is a 69 × 13 contingency table with 33, 704 agents.Tab. 1 shows 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-MCMCand SIM-NN, which only operate on the Λ level. We also outperform the competitive SIT-MCMC approach [46], which can operate in the discrete ODM level, since in the limit of C T our method can move to high T probability regions much faster than MCMC due to the optimisation of the SIM Λ parameters. This effectively deflates the Λ estimator variance relative to the variance of the Λ samples in MCMC. Only in the absence of rich C T data does SIT-MCMC outperform our method (total constrained case) as it limits the support of θ to [0, 2] 2 which acts as regularisation. Letting the support of θ span the entire R 2 allows information carried by stronger Λ constraints (such as Λ +• ) to permeate to the T space much faster in GeNSIT compared to SIT-MCMC. The plethora of sources of uncertainty ranging from the HW-SDE ( 6) to the combinatorial nature of T C suggest that GENSIT is faster in reconstructing ground truth ODMs compared to a fully Bayesian approach such as SIT-MCMC. -: This case is not handled by the approach mentioned in the column.ORIGIN- DESTINATION MATRIX C TARGET M ↓ SRMSE (M (1:N ) , T * ) ↑ 99% CP (M (1:N ) , T * ) GENSIT [SIT- MCMC] [SIM- NN] [SIM- MCMC] GENSIT [SIT- MCMC] [SIM- NN] [SIM- MCMC] DISJOINT JOINT DISJOINT JOINT Λ 2.(a) ODMs with closed-form (tractable) T distributions (10) . (b) ODMs with intractable T distribution (13) , where conditioning Λ on C is problematic.ORIGIN- DESTINATION MATRIX C TARGET M ↓ SRMSE (M (1:N ) , T * ) ↑ 99% CP (M (1:N ) , T * ) GENSIT [SIT- MCMC] GENSIT [SIT- MCMC] DISJOINT JOINT DISJOINT JOINT Λ 2.Table 1 : Ground truth T * validation metrics comparing our method against [13, 17, 46] in the Λ and T levels across constraint sets C and σ = 0.141 (best) for the Cambridge dataset. On a C basis the best metric in T,Λ spaces is emphasised 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, 2] 2 instead of R 2 , which has a regularisation effect. In the absence of substantial C data, this effect is more pronounced. See Tab. 6 (App. E.1) for full table across multiple σ regimes.Fig. 6 supports this claim by shedding light on the convergence rate of running 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 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 = 10 4 . In the Disjoint GeNSIT the information encoded in larger C T is not propagated to the Λ updates. As a result, no SRMSE or CP improvements are detected in the course of N .We apply our method to the Washington dataset, where the ground truth ODM is a 179 × 179 contingency table with 200, 029 agents. Tab. 2 reports the reconstruction error (SRMSE) and the % coverage of the ground truth cells (CP) in Λ and T. Comparisons against SIM-MCMCand SIT-MCMCwere infeasible due to their high computational complexity (300 hours to obtain 500 samples on a 32-core NVIDIA GPU). Instead, we leveraged GMEL [26] which operates only in the [26, 17] in the Λ and T levels for C = {T ++ , T X train } and σ = 0.141 (best) for the Washington dataset. Arrow ↑ indicates higher values are better, and vice versa. We achieve the best error (SRMSE) and ground truth coverage (CP) overall (see bold cells). The observed data D leveraged to train GENSIT, SIM-NN is a small subset of the data required to train GMEL (y is a column vector of Y). See Tab. 7 (App. E.2) for full table across multiple σ regimes.M ↓ SRMSE (M (1:N ) , T * ) ↑ SSI (M (1:N ) , T * ) ↑ 99% CP (M (1:N ) , T * ) [GMEL] Y, C Λ 2.continuous Λ space by learning a mapping between a large feature space D and T X train . Tab. 2 shows that we outperform both GMEL and SIM-NN in terms of reconstruction error and coverage. 
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 Λ. We surmount the limitations of methods which operate strictly on Λ space [13, 17, 26] and of methods that incur a large computational cost [46] . We accomplish this by introducing the first framework operating on both T, Λ that scales linearly with the number of origin-destination pairs IJ. We offer enhanced reconstruction error and coverage of the ground truth ODMs in Cambridge, UK and Washington, DC. Although NNs require a much larger number of internal parameters to be calibrated relative to MCMC, their embedding of physics models regularises this parameter space and prevents over-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-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.[46] Ioannis Zachos, Theodoros Damoulas, and Mark Girolami. B.1 Continuous ODMA fundamental assumption is that the continuous agent trip intensity Λ 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 Λ. 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 Λ 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 [11] and empirically shown in Tabs. 1,2 of our Experimental Section. In poorly constrained settings, we rely more on Λ for spatial imputation, which increases reconstruction error. 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.In the ODMs considered we have constructed agent utility functions that leverage the log destination attraction at 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 to derive the functional forms of Λ in equations ( 3), (4), and ( 5). The SIM intensity is equivalent to the multinomial logit [29] , 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 o ∈ R I in addition to the existing x, C in the totally constrained SIM (total = Λ ++ ). Then, we need to maximise:-I,J i,j Λ ij (log(Λ ij ) -o i -αx j + βc ij ) -µ( I,J ij Λ ij -Λ ++ ),where µ is the Lagrange multiplier. This yieldsΛ ij = Λ ++ (exp(o i + αx j -βc ij )) I,J i,j exp(o i + αx j -βc ij ) . B.2 Discrete ODM B.2. 
 
 Appendix:  


Decision: 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.

Reviews: review: paper_summary: paper_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.
strength_weakness: strength_weakness: 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.) 

The 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.

The 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.

The 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?

The 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.
questions: 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.

Would an ablation study be feasible? Or is there a reason that all the parts of the algorithm must be included together?
limitations: limitations: The paper mentions that the problem of more complex cell structure constraints for population synthesis remains unsolved.

There 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).

The 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.
paper_summary: paper_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.
strength_weakness: strength_weakness: I would flag two main concerns:
1) 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
2) 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.
questions: 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?
2) 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: limitations: NA
paper_summary: paper_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.
strength_weakness: strength_weakness: 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.

To 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+}$?

The authors aim for precision and rigor, but the text falls short and is confusing. Here are some examples:

- 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?
- 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.
- 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)?
- In practice, how can $\mathcal{C}_{\Lambda}$ be known to be used as input in the algorithm?
- What happens if the inputs $\mathcal{C}_{\Lambda}$ are completely inconsistent with the inputs $\mathcal{C}_T$? For example, $C_{++} = 20000$ and $\Lambda_{++} = 100$?
- 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.
- "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?
- "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})$).
- "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.
questions: questions: The real data sets are relatively small (69x13 and 179x179). Can you run large synthetic examples to study the scalability of GeNSIT?
limitations: limitations: They addressed adequately the limitations.
paper_summary: paper_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.
strength_weakness: strength_weakness: - 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.
- 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.
questions: 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?
- 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: limitations: NA
paper_summary: paper_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.
strength_weakness: strength_weakness: 1.	The 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.
2.	The 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.
questions: questions: 1.	In 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: limitations: The major limitations are discussed in the experiment section, and no significant negative social impact can be observed in this paper.
paper_summary: paper_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.
strength_weakness: strength_weakness: 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. 

2. Computational cost should be further discussed. The computational complexity is given by $O(NEJ(\tau|W|+I))$, what is the scale of each parameter?
questions: 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?

2. Is this method robust to outliers?

3. This method seems very complicated to me. Can you summarize the intuition behind it in a short paragraph?
limitations: limitations: See comments above.

Comments: title: Thanks for the response
Acknowledgement
Acknowledgement of rebuttal and thanks
response to Authors' rebuttal
comment: Thanks for the response. I will keep my original score.
Thanks for responding to my questions. Based on the reply, I have increased my score.
Dear Reviewers,

Thank you for your reviews and taking the time to respond to our rebuttal. If you have any additional comments or requests we would be happy to address them in the remaining reviewer-author discussion period. If you believe we have already addressed your corrections/questions we would appreciate if this could be reflected in your updated score.

Many Thanks,
Authors
I thank the authors for providing this detailed rebuttal that addresses all my points. In the global rebuttal it seems you have made many improvements to the paper that cover all my suggested improvements. I think this definitely improves the paper on the dimensions of clarity and presentation in particular.

I have updated my 'presentation' score to 3: good
I thank the authors for their responses. Accordingly, I have updated my score.
I appreciate the response, and everything generally makes sense. I will maintain my score for now.