1 Introduction↩︎

Origin–destination (OD) flow describes movements of people or goods between areas, reflecting social and economic activity at a macro level and supporting applications such as urban planning and behavioral analysis [1]–[3]. However, large-scale data collection faces statistical and privacy challenges [4], [5], and generalizable models are needed to inform urban development, especially in data-scarce or early-planning regions [6]. Consequently, OD flow prediction and generation have become central research topics [7]. Prediction methods estimate flows from area attributes (e.g., socio-demographics), while generation methods further employ random seeds and generative models to capture diverse mobility patterns [8], [9]. Despite their conceptual differences, both tasks share common modeling principles, and methods are often evaluated jointly [6], [10]. Early approaches, including gravity [11] and radiation models [12], were followed by machine learning techniques such as Random Forest [13], [14] and SVR [15], [16]. With the rise of deep learning [17], models such as DGM [4] and GMEL [8] achieved notable progress, and recent graph-based methods further advanced prediction and generation by modeling areas as attributed nodes and OD flows as directed weighted edges [6], [8], [18].

Although existing state-of-the-art methods represent area networks as graph structures to exploit properties such as translational and rotational equivalent [8], [19], this representation often overlooks easily accessible geographic information, thereby complicating the learning process [20]. In practice, relationships between areas are encoded in matrix form through pairwise adjacency and straight-line distances [6], [21]. While this formulation is theoretically lossless and supports network reconstruction, it fails to directly preserve fundamental geometric properties such as collinearity, which can be readily obtained from coordinate systems, and geodesic distances, which are essential but difficult to infer from pairwise matrices [22]. Recovering such properties is challenging, and their absence undermines the model’s ability to capture long-range and multi-area dependencies that are crucial for OD flow modeling [23]. Consequently, current feature representation and encoding strategies remain limited in representing global area-to-area relationships.

This limitation highlights a central challenge in OD flow modeling: how to represent and encode complex inter-area relationships. Graph structures offer desirable invariances and computational advantages, yet they often fail to capture long-range and multi-area geographic interactions [24]. Coordinate-based representations encode geometric attributes directly but are sensitive to the choice of reference origin, scale, and orientation [25]. Besides, some geospatial attributes (e.g., geodesic distance) remain nontrivial to recover from matrix representations [22]. Accordingly, integrating complementary representation methods and developing encoding–decoding strategies that preserve geometric attributes while remaining insensitive to reference choices are essential for improving OD flow prediction and generation.

This paper introduces GeoFlow, a novel framework for OD flow prediction and generation that systematically integrates geospatial attributes and encoder–decoder design. We first perform geospatial attribute augmentation by incorporating \(k\)-hop distance, geodesic distance, and relative position alongside straight-line distance, allowing the model to access key spatial relationships directly rather than extracting them from matrices. In the encoder, we propose geometry-intrinsic fusion encoder, combining Cartesian coordinates with graph-based representations: a graph attention module encodes intrinsic area attributes (e.g., points of interest) and aggregates influences from neighboring areas, while carefully chosen coordinate centers and orientations preserve relative positions, enabling the model to capture long-range and multi-area dependencies. In the decoder, we introduce an axial-global attention mechanism that models competitive relationships between areas while reducing computational costs, making joint attention across hundreds of areas feasible. Finally, for flow generation, we adopt flow matching models, which provide higher stability, efficiency, and generalization than diffusion models [26], [27]. To the best of our knowledge, this is the first work to apply flow matching in the OD flow generation task.

Extensive experiments demonstrate that GeoFlow markedly improves both OD flow prediction and generation. For the prediction task, it consistently achieves substantially lower errors than existing baselines, reflecting a clear advantage. In the generation task, GeoFlow attains superior performance in reconstruction accuracy and produces more authentic and diverse samples. Ablation studies on geospatial representations show that even lightweight attribute augmentations provide notable gains, highlighting the practical value of these previously overlooked features, and ablation experiments on the encoder and decoder architectures confirm their effectiveness. Furthermore, analytical experiments demonstrate the rationale behind the overall design and offer deeper insight into the role of geospatial attributes.

Our contributions can be summarized in three main parts:

1. We propose GeoFlow, a novel method for OD flow prediction and generation that systematically enhances feature representation and model architecture. Each area is augmented with crucial geospatial attributes to capture long-range and multi-area relationships. A geometric-intrinsic fusion encoder module encodes area properties, while an axial-global attention mechanism models competitive dependencies among areas.

2. We evaluate GeoFlow on both OD flow prediction and generation tasks. GeoFlow achieves substantial improvements in prediction accuracy and generates more realistic and diverse samples.

3. We conduct ablation experiments on model components, showing that the proposed approaches substantially improve model performance. Additionally, analytical experiments on geospatial features provide deeper insights into the dynamics of OD flow modeling.

2 Related Work↩︎

Existing work on origin–destination (OD) flow prediction or generation falls into two broad categories: principle-driven methods and data-driven methods. Principle-driven methods derive flows from explicit mechanistic physical models and statistical formulations. [11] employs the gravity model to explain intercity human mobility. [28] introduces the notion of competitive opportunities and incorporates probability theory to estimate OD flows. [12] propose a parameter-free radiation model and enhance the accuracy of flow modeling. These formulations provide transparent explanations of aggregate flows but lack flexibility for complex empirical patterns. Data-driven methods, which learn flow patterns directly from observed mobility data, have emerged and revitalized the field in recent years. Many studies apply classical machine learning algorithms to OD flow generation. Random Forest (RF) demonstrates strong potential in the task [14]. Gradient Boosting Regression Trees (GBRT) enhance predictive accuracy through boosting and have been effectively applied to urban OD prediction [29]. Support Vector Regression (SVR) employs kernel methods to estimate flows between regions based on their urban attributes [16]. Despite their demonstrated success, classical machine learning methods have limitations in fully capturing the complex dynamics and high-dimensional interactions inherent in urban OD flows, paving the way for deep learning-based techniques. The Deep Gravity Model (DGM), inspired by traditional gravity models, uses multi-layer perceptrons to estimate flow distributions [4]. Graph-based methods, including NetGAN [18] and Geographical Multi-task Embedding Learning (GMEL) [8], leverage random walks or graph neural networks to capture spatial and structural information, enabling more accurate OD flow generation. More recent approaches, such as DiffODGen [10] and WEDAN [6], employ diffusion models to model graph topology and edge weights conditioned on urban attributes. Despite their ability to capture complex spatial and structural patterns, these methods typically abstract areas as attributed nodes, overlooking spatial proximity and inter-area interactions and thereby hindering further advances. TransFlower [21] acknowledges the importance of relative location in capturing directional relations between areas; however, its reliance on pairwise attention and limited integration of topological structure restricts its capacity to model OD flows.

3 Method↩︎

3.1 Problem Definition, Formulation, and Notation↩︎

We consider OD flow modeling as a unified paradigm: the core problem is to characterize flows among areas under given attributes. The paradigm can be instantiated in two ways: (1) When the goal is to estimate missing or future entries, the formulation reduces to prediction; (2) When the goal is to capture the full stochastic behavior of flows, it becomes generation. Both tasks share the same underlying setting but differ in the learning objectives.

Formally, let \(R = \{A_i\}_{i=1}^n\) denote a region partitioned into \(n\) non-overlapping areas; each area \(A_i\) is associated with a feature vector \({\boldsymbol{x}}_i\in\mathbb{R}^{d_\mathrm{x}}\), and the collection is denoted by a feature matrix \({\boldsymbol{X}}=[{\boldsymbol{x}}_1; {\boldsymbol{x}}_2; \dots; {\boldsymbol{x}}_n]\in\mathbb{R}^{n\times d_\mathrm{x}}\). OD flows are represented as a non-negative matrix \({\boldsymbol{F}}\in\mathbb{R}_{\ge0}^{n\times n}\) where entry \({F}_{i,j}\) specifies the flow from \(A_i\) to \(A_j\), while inter-area geographic or topological information is captured by a tensor \({\boldsymbol{G}}\in {\mathbb{R}}^{n\times n\times d_\mathrm{g}}\), where \(d_\mathrm{g}\) is the dimension of the per-pair geographic descriptor, including relationships such as adjacency and distance. A single data instance is represented as \(D=\{R,{\boldsymbol{X}},{\boldsymbol{G}},{\boldsymbol{F}}\}\).

Prediction as supervised estimation. Given a training index set \(\Omega\subseteq\{1,\dots,n\}^2\), the prediction task seeks to estimate \(\hat{{\boldsymbol{F}}}=f({\boldsymbol{X}},{\boldsymbol{G}};{\boldsymbol{\theta}})\) by minimizing the discrepancy between observed and predicted flows: \[\min_{{\boldsymbol{\theta}}}\; \mathcal{J}_{\mathrm{pred}}({\boldsymbol{\theta}}) = \frac{1}{|\Omega|}\sum_{(i,j)\in\Omega}\ell\big({F}_{i,j}, \hat{F}_{i,j}\big),\] where \(\ell(\cdot,\cdot)\) is a per-entry loss (e.g., squared error).

Generation as distributional modeling. Beyond point estimation, the generation task requires capturing the full distributional structure of OD flows. This can be formulated via a conditional generator \(\hat{{\boldsymbol{F}}} = g({\boldsymbol{X}}, {\boldsymbol{G}}, z; {\boldsymbol{\theta}})\), with latent seed \(z\sim p_z\), where \(p_z\) is the seed distribution, optimized under \[\begin{align} \min_{{\boldsymbol{\theta}}} &\mathcal{J}_{\mathrm{gen}}({\boldsymbol{\theta}}) = \\ & \mathbb{E}_{z\sim p_z} \big[ \ell_{\mathrm{rec}}({\boldsymbol{F}}, \hat{{\boldsymbol{F}}}) + \alpha\ell_{\mathrm{dist}}({\boldsymbol{F}},\hat{{\boldsymbol{F}}}) + \beta\ell_{\mathrm{div}}(\hat{{\boldsymbol{F}}})\big], \end{align}\] where \(\ell_{\mathrm{rec}}\) enforces reconstruction fidelity, \(\ell_{\mathrm{dist}}\) measures distributional discrepancy between realistic data and generated ones, and \(\ell_{\mathrm{div}}\) promotes generation diversity. Hyperparameters \(\alpha\) and \(\beta\) balance reconstruction fidelity, distributional realism, and sample diversity.

3.2 GeoFlow↩︎

Unlike trajectory forecasting, which requires reconstructing detailed travel paths, OD flow tasks focus on aggregate movements between areas. A straightforward baseline concatenates origin and destination features to predict or generate flows from pairwise representations. However, this formulation treats areas as isolated nodes and neglects critical geographic and multi-area dependencies.

To address these limitations, we partition area attributes into two categories and design dedicated encoders. The first, geometric attributes, includes coordinates, pairwise distances, and spatial descriptors such as geodesic proximity, which support global and long-range interaction modeling. The second, intrinsic attributes, covers local covariates such as points of interest (POIs) and socio-demographics, which characterize production and attraction capacity and primarily affect nearby areas. These two categories differ in semantics, scale, and propagation behavior: geometric attributes impose structural and long-range constraints, whereas intrinsic attributes capture localized demand and short-range interactions. Accordingly, we encode them separately and fuse their representations downstream to integrate global structure with local context. For decoding, we introduce an axial attention mechanism that models competitive dependencies in human mobility: OD pairs sharing the same origin or destination exhibit stronger correlations than unrelated pairs. This mechanism is further combined with global attention, resulting in a decoder that enhances training efficiency while achieving superior performance within comparable computational budgets. Fig. 1 illustrates the overall architecture of the GeoFlow framework.

The following sections elaborate on these components: Sec. 3.2.1 describes the augmented geospatial attributes, Sec. 3.2.2 presents the encoding strategies, and Sec. 3.2.3 details the axial-global attention mechanism.

For notational simplicity, we use \({A}_{i,j}\) (or \({A}_{i,j,k}\)) to denote the \((i,j)\)-th entry of a two-dimensional matrix (or the \((i,j,k)\)-th entry of a three-dimensional tensor) \({\boldsymbol{A}}\). We write \({\boldsymbol{A}}_{i,:}\) for the \(i\)-th row of \({\boldsymbol{A}}\) and \({\boldsymbol{A}}_{:,j}\) for its \(j\)-th column. Analogous notation similarly applies to tensors of order three and higher.

3.2.1 Geospatial Feature Augmentation↩︎

Since graph-based approaches typically store geographical relationships as adjacency or distance matrices, recovering relative positions or geodesic distances from such representations is difficult, limiting the model’s ability to capture essential factors that shape travel behavior [20], [22]. Rather than relying solely on deeper or wider networks, explicitly providing these salient features offers a direct and cost-effective alternative, while also supplying inductive signals that stabilize and accelerate training.

In this work, we augment each area pair with three additional geographic relation descriptors: relative position, \(\boldsymbol{k}\)-hop distance, and geodesic distance. The relative position is defined as the coordinate difference between two areas under a Cartesian system assigned to the region. To ensure translational and rotational invariance while retaining meaningful orientation, we place the origin at the region’s geometric center, align the \(x\)-axis with the principal orientation (the direction of greatest variance), and normalize coordinates to \([-1,1]\), thereby mitigating scale differences across regions. The \(k\)-hop distance is the length of the shortest unweighted path between two areas, reflecting topological proximity. The geodesic distance, or free-flow distance, is the minimal travel length along the network, capturing realistic accessibility. Fig. 1 upper left illustrates an example of the augmented geospatial features. By incorporating these attributes, GeoFlow provides richer geographic context, improving both predictive accuracy and training stability in OD flow modeling.

3.2.2 Encoder: Geometric–Intrinsic Fusion↩︎

To effectively leverage the enriched geospatial features, we employ a dedicated encoder that fuses geometric and intrinsic area attributes. As shown in Fig. 1 upper middle, for the geometric attributes, which include relative position, straight-line distance, \(k\)-hop distance, and geodesic distance, we encode each pairwise attribute vector using a multilayer perceptron (MLP). These geometric attributes are normalized to remove scale-induced biases. By adopting a consistent regional coordinate system and applying systematic normalization, the model can more readily infer derived geometric cues (e.g., relative bearing) within a fixed domain. Normalized relative positions also reduce biases introduced by absolute location and orientation, obviating the need for explicit trigonometric encodings. For a region with \(n\) areas, let the geometric information be denoted by the tensor \({\boldsymbol{G}}\in\mathbb{R}^{n\times n\times d_\mathrm{g}}\). We employ an MLP \(\phi_{\mathrm{g}}:\mathbb{R}^{d_\mathrm{g}}\to\mathbb{R}^{d_\mathrm{e}}\) elementwise to obtain an encoded geographic embedding \[{\boldsymbol{E}}_{\mathrm{g}} = \phi_{\mathrm{g}}({\boldsymbol{G}})\in\mathbb{R}^{n\times n\times d_\mathrm{e}},\] where \(d_\mathrm{e}\) is the dimension of embedding \({\boldsymbol{E}}_{\mathrm{g}}\).

For intrinsic attributes, given their dependence on neighboring areas, we adopt a graph attention-based encoding [30]. Unlike prior approaches that rely solely on local message passing and are therefore constrained by network depth, we introduce an axial-attention design to enhance information propagation efficiency. An area’s influence on others depends on both its intrinsic attributes and geographic relations. Accordingly, GeoFlow conditions attention on the geographic tensor when computing inter-area interactions. Formally, let \({\boldsymbol{X}}\in\mathbb{R}^{n\times d_\mathrm{x}}\) denote the intrinsic attribute matrix. We first project \({\boldsymbol{X}}\) as intrinsic embedding \({\boldsymbol{E}}_\mathrm{x}\) via an MLP \(\phi_\mathrm{x}:\mathbb{R}^{d_\mathrm{x}}\rightarrow\mathbb{R}^{d_\mathrm{e}}\): \[{\boldsymbol{E}}_\mathrm{x} = \phi_\mathrm{x}({\boldsymbol{X}}),\] and then update \({\boldsymbol{E}}_\mathrm{x}\) via an attribute aggregation. Specifically, we first use a learnable weighted vector \({\boldsymbol{w}}_\mathrm{geo} \in {\mathbb{R}}^{d_\mathrm{g}}\) to produce a score matrix \({\boldsymbol{S}}\in {\mathbb{R}}^{n\times n}\): \[\label{eq:32w32geo} {\boldsymbol{S}}= \operatorname{softmax}({\boldsymbol{G}}{\boldsymbol{w}}_{\mathrm{geo}}).\tag{1}\] Intuitively, an area’s intrinsic attributes are influenced by information that propagates from other areas through the spatial structure. The score matrix \({\boldsymbol{S}}\) is designed to capture this influence pattern, where \({S}_{i,j}\) denotes the influence weight from area \(A_j\) to area \(A_i\), inferred from their geographic descriptors. Next, we use \({\boldsymbol{S}}\) to modulate query formation from intrinsic attributes and compute pairwise aggregation weights \({\boldsymbol{A}}\) via scaled dot-product: \[{\boldsymbol{A}}= \operatorname{softmax}\Big( \frac{({\boldsymbol{S}}{\boldsymbol{X}}{\boldsymbol{W}}_\mathrm{qx})({\boldsymbol{X}}{\boldsymbol{W}}_\mathrm{kx})^\top}{\sqrt{d_\mathrm{h}}} \Big),\] \[(\widetilde{{\boldsymbol{E}}}_\mathrm{x})_{i,:} = \sum_{k=1}^{n} {A}_{i,k}\, ({\boldsymbol{E}}_\mathrm{x})_{k,:}.\] where \({\boldsymbol{W}}_\mathrm{qx}, {\boldsymbol{W}}_\mathrm{kx} \in {\mathbb{R}}^{d_\mathrm{x}\times d_\mathrm{h}}\) are learnable projection matrices, and \(d_\mathrm{h}\) represents the hidden dimension.

Then, we employ an MLP to integrate the geographic embedding \({\boldsymbol{E}}_\mathrm{g}\) and the intrinsic embedding \(\widetilde{{\boldsymbol{E}}}_\mathrm{x}\), producing the joint representation \({\boldsymbol{Z}}_{0} \in {\mathbb{R}}^{n\times n\times d_\mathrm{z}}\). Formally, let \(\phi_{\mathrm{f}}:\mathbb{R}^{3d_\mathrm{e}}\to\mathbb{R}^{d_\mathrm{z}}\) denote the fusion MLP. For each OD pair \((i,j)\), the fused embedding is obtained as \[\label{eq:32z} ({\boldsymbol{Z}}_{0})_{i,j,:} = \phi_\mathrm{f}\big(\,\big[({\boldsymbol{E}}_\mathrm{g})_{i,j,:} \;\Vert\; (\widetilde{{\boldsymbol{E}}}_\mathrm{x})_{i,:} \;\Vert\; (\widetilde{{\boldsymbol{E}}}_\mathrm{x})_{j,:}\big]\,\big),\tag{2}\] where \(\Vert\) denotes concatenation. The joint embedding \({\boldsymbol{Z}}_{0}\) serves as the unified representation for subsequent OD flow prediction and generation.

3.2.3 Decoder: Axial–Global Attention↩︎

Naive full attention can capture interactions among OD pairs [31]. Still, it is computationally prohibitive: for a region with \(n\) areas, the \(n^2\) OD pairs incur an \(O(n^4)\) cost for pairwise attention, limiting the depth of dynamics the model can afford under practical budgets [32]. Empirically, OD pairs that share an origin or share a destination are more strongly coupled. Fixing an origin induces comparisons among candidate destinations, while a destination’s finite capacity differentiates travelers from different origins, underscoring the relevance of origin- and destination-conditioned interactions. Accordingly, we adopt an axial attention scheme that alternately applies origin-wise and destination-wise attention for stacked blocks. This design efficiently captures origin-conditioned and destination-conditioned relational patterns without evaluating full pairwise attention, reducing computational cost to \(O(n^3)\). To preserve the rich interactions, we introduce a redesigned compact global-attention mechanism to capture long-range and multi-area dependencies with \(O(n^2)\) computational complexity. The decoder takes the unified OD representation \({\boldsymbol{Z}}_{0}\) from Eq. 2 as input and refines it with a stack of \(L_\mathrm{a}\) axial-attention layers and \(L_\mathrm{g}\) compact global-attention layers.

Axial-attention mechanism. Each axial-attention layer consists of an origin-wise sublayer followed by a destination-wise sublayer, and thus performs two successive updates of the embedding. For the \((k{+}1)\)-th axial-attention layer, where \(k \in \{0,1,\dots,L_\mathrm{a}-1\}\), we use \({\boldsymbol{Z}}_{2k}\) and \({\boldsymbol{Z}}_{2k{+}1}\) to represent the inputs to the corresponding origin-wise and destination-wise attention sublayers, respectively. Specifically, \({\boldsymbol{Z}}_{0}\) is fed into the first origin-wise attention layer as input.

We first introduce the standard scaled dot-product attention [31] as \[\operatorname{Attn}({\boldsymbol{Q}},{\boldsymbol{K}},{\boldsymbol{V}}) = \operatorname{softmax}\Big(\frac{({\boldsymbol{Q}}{\boldsymbol{W}}_\mathrm{q})({\boldsymbol{K}}{\boldsymbol{W}}_\mathrm{k})^\top}{\sqrt{d_\mathrm{k}}}\Big) {\boldsymbol{V}}{\boldsymbol{W}}_\mathrm{v},\] where \({\boldsymbol{Q}}\), \({\boldsymbol{K}}\), \({\boldsymbol{V}}\) are query, key, value, respectively. \({\boldsymbol{W}}_\mathrm{q}\), \({\boldsymbol{W}}_\mathrm{k}\), \({\boldsymbol{W}}_\mathrm{v}\) are learnable matrices. \(d_\mathrm{k}\) is the feature dimension. The formulation can be readily extended to the multi-head one; here, we omit it for brevity. The self-attention can be written as \[\operatorname{SelfAttn}({\boldsymbol{M}}) = \operatorname{Attn}({\boldsymbol{M}}, {\boldsymbol{M}}, {\boldsymbol{M}}).\] The origin-wise and destination-wise attention for the \((k{+}1)\)-th axial-attention layer is then formulated as \[\big(\mathcal{A}^{\mathrm{orig}}_{k{+}1}({\boldsymbol{Z}}_{2k})\big)_{i,:,:} = \operatorname{SelfAttn}\big( ({\boldsymbol{Z}}_{2k})_{i,:,:} \big)\] \[\big(\mathcal{A}^{\mathrm{dest}}_{k{+}1}({\boldsymbol{Z}}_{2k{+}1})\big)_{:,j,:} = \operatorname{SelfAttn}\big( ({\boldsymbol{Z}}_{2k{+}1})_{:,j,:} \big).\] We define the axial-attention mechanism for the \((k{+}1)\)-th layer as \[\mathcal{A}^{\mathrm{axial}}_{k+1}({\boldsymbol{Z}}_{2k}) = \Big(\mathcal{A}^{\mathrm{dest}}_{k+1} \bullet \mathcal{A}^{\mathrm{orig}}_{k+1}\Big)({\boldsymbol{Z}}_{2k}),\] where \(\bullet\) denotes function composition. After \(L_\mathrm{a}\) axial-attention layers, we obtain \({\boldsymbol{Z}}_{2L_\mathrm{a}}\), the output representation of the axial-attention block. It explicitly models interactions among areas in direct competition (e.g., sharing the same origin or destination context), while implicitly aggregating global information through iterative propagation across layers.

Compact global-attention mechanism. To avoid losing the complex, multi-area interactions among OD pairs, we introduce a computation-efficient compact global-attention pathway that complements the axial-attention blocks. Instead of applying full attention to all \(n^2\) OD pairs, we first aggregate OD information at the area level, then further apply global-attention layers to provide an additional low-resolution channel for information exchange across the entire region. Given \({\boldsymbol{Z}}_{2L_\mathrm{a}} \in \mathbb{R}^{n\times n\times d_\mathrm{z}}\), we compute origin-wise and destination-wise pooled representations by average pooling: \[({\boldsymbol{H}}_{\mathrm{orig}})_{i,:} = \frac{1}{n}\sum_{j=1}^{n} ({\boldsymbol{Z}}_{2L_\mathrm{a}})_{i,j,:},\] \[({\boldsymbol{H}}_{\mathrm{dest}})_{j,:} = \frac{1}{n}\sum_{i=1}^{n} ({\boldsymbol{Z}}_{2L_\mathrm{a}})_{i,j,:},\] where \({\boldsymbol{H}}_{\mathrm{orig}}, {\boldsymbol{H}}_{\mathrm{dest}} \in \mathbb{R}^{n\times d_\mathrm{z}}\) summarize OD information for each origin and destination, respectively. Then we concatenate them and pass the result through a projection: \[{\boldsymbol{H}}_0 = \phi_\mathrm{h}\big([{\boldsymbol{H}}_{\mathrm{orig}} \;\Vert\; {\boldsymbol{H}}_{\mathrm{dest}}]\big),\] where \(\phi_\mathrm{h}: \mathbb{R}^{2 d_\mathrm{z}} \rightarrow \mathbb{R}^{d_\mathrm{z}}\) is a projection MLP. Here, we interpret \({\boldsymbol{H}}_0 \in {\mathbb{R}}^{n\times d_\mathrm{z}}\) as an area-level latent state that consolidates OD context from the perspectives of serving as the origin and the destination, which is subsequently refined by global self-attention \[{\boldsymbol{H}}_1 = \operatorname{SelfAttn}({\boldsymbol{H}}_0).\] We repeat this operation \(L_\mathrm{g}\) times, yielding the output \({\boldsymbol{H}}_{L_\mathrm{g}}\), which encodes explicit global interactions across all areas via a compact set of area-level embeddings. Through this global-attention pathway, all areas interact explicitly via a compact set of representations at the aggregated level, enabling GeoFlow to capture high-order, long-range, and origin– and destination-crossing dependencies while keeping the computational cost manageable.

Then, we unsqueeze the \({\boldsymbol{H}}_{L_\mathrm{g}}\) along the origin and destination axes to form origin-aligned and destination-aligned views \(\widetilde{{\boldsymbol{H}}}_{\mathrm{orig}}\) and \(\widetilde{{\boldsymbol{H}}}_{\mathrm{dest}}\in\mathbb{R}^{n\times n\times d_\mathrm{z}}\) by \[(\widetilde{{\boldsymbol{H}}}_{\mathrm{orig}})_{i,j,:} = ({\boldsymbol{H}}_{L_\mathrm{g}})_{i,:}, \quad \forall\;j\] \[(\widetilde{{\boldsymbol{H}}}_{\mathrm{dest}})_{i,j,:} = ({\boldsymbol{H}}_{L_\mathrm{g}})_{j,:}, \quad \forall\;i,\] and expand them back to the \(n\times n\) OD grid via a Hadamard product so that each OD pair receives the corresponding origin- and destination-conditioned context, \[\widetilde{{\boldsymbol{H}}}_\mathrm{OD} = \widetilde{{\boldsymbol{H}}}_{\mathrm{orig}} \circ \widetilde{{\boldsymbol{H}}}_{\mathrm{dest}}.\] Finally, we pass the resulting OD-level embedding through an MLP to obtain the final decoder output, \[{\boldsymbol{Z}}_\mathrm{dec} = \phi_{\mathrm{out}}(\widetilde{{\boldsymbol{H}}}_\mathrm{OD}).\]

3.3 OD Flow Prediction and Generation↩︎

GeoFlow provides an expressive, geometry-aware representation of intra-region area relationships that preserves multi-area spatial and relational cues necessary for accurate OD modeling. We apply GeoFlow to both OD flow prediction and generation tasks.

Prediction. For the prediction task, we append an MLP to the decoder output to produce the predicted OD flow matrix. Let \({\boldsymbol{Z}}_{\mathrm{dec}}\in\mathbb{R}^{n\times n\times d_\mathrm{z}}\) denote the decoder output; we project it via MLP \(\phi_\mathrm{p}: \mathbb{R}^{d_\mathrm{z}} \to \mathbb{R}\), as \[\label{eq:32pred32mlp} \hat{{\boldsymbol{F}}} = \phi_\mathrm{p}({\boldsymbol{Z}}_{\mathrm{dec}})\tag{3}\] where \(\hat{{\boldsymbol{F}}}\in\mathbb{R}^{n\times n}\) is the predicted OD flow matrix. The prediction loss is the average element-wise squared error \[\mathcal{L}_{\mathrm{pred}} \;=\; \frac{1}{n^2}\,\big\|\hat{{\boldsymbol{F}}}-{\boldsymbol{F}}\big\|_F^2,\] where \({\boldsymbol{F}}\) is the ground-truth OD flow matrix and \(\|\cdot\|_F\) denotes the Frobenius norm.

Generation. For the generation task, we adopt continuous-time flow matching [26] conditioned on GeoFlow’s output embedding. We embed the intermediate instantaneous state at timestamp \(t\), denoted by \({\boldsymbol{F}}_t\), via an extra MLP encoder \(\psi\) and fuse it with geometry and intrinsic embeddings. Formally, in Eq. 2 , the fused embedding is turned to obtain as \[({\boldsymbol{Z}}_{t})_{i,j,:} = \phi^\mathrm{gen}_\mathrm{f}\big(\,\big[({\boldsymbol{E}}_\mathrm{g})_{i,j,:} \Vert (\widetilde{{\boldsymbol{E}}}_\mathrm{x})_{i,:} \Vert (\widetilde{{\boldsymbol{E}}}_\mathrm{x})_{j,:} \Vert \psi({\boldsymbol{F}}_t)_{i,j,:} \big]\,\big).\] The decoder then maps \({\boldsymbol{Z}}_t\) to a time-conditioned OD representation \({\boldsymbol{Z}}_{\mathrm{dec},t}\), which is followed by an MLP to predict the velocity field \(\hat{\boldsymbol{V}}({\boldsymbol{F}}_t, t, {\boldsymbol{Z}}_{\mathrm{dec},t}; {\boldsymbol{\theta}}) \in \mathbb{R}^{n\times n}\). The model is trained by mean squared error against the analytic target velocity induced by a perturbation schedule. The loss function can be written as \[\mathcal{L}_{\mathrm{gen}}=\mathbb{E}_{D,t,\xi}\big\|\hat{\boldsymbol{V}}-{\boldsymbol{V}}_{\mathrm{tgt}}({\boldsymbol{F}}_t\mid {\boldsymbol{F}})\big\|_F^2.\]

For a common linear-Gaussian schedule \({\boldsymbol{F}}_t=\alpha(t)\boldsymbol{\xi}+\sigma(t){\boldsymbol{F}}\), where \(\alpha(t)\) and \(\sigma(t)\) are predefined scalar schedule functions that control the interpolation from Gaussian noise \(\boldsymbol{\xi}\sim\mathcal{N}(\boldsymbol{0},\boldsymbol{I})\) to data \({\boldsymbol{F}}\). In particular, \(t=0\) corresponds to the noise distribution and \(t=1\) corresponds to the data distribution. \({\boldsymbol{V}}_\mathrm{tgt}\) represents the target velocity. At inference, we integrate the learned ordinary differential equation \(\mathrm{d}\hat{{\boldsymbol{F}}}_t/\mathrm{d}t=\hat{\boldsymbol{V}}\) from \(\hat{{\boldsymbol{F}}}_0=\xi\) to obtain \(\hat{{\boldsymbol{F}}}_1\). This formulation keeps training stable, conditions sampling on geometry and area identity, and yields faster inference than diffusion methods. For more details, please refer to Appendix 6.1.

4 Experiments↩︎

3.2pt

4.1 Setup↩︎

Objectives. We aim to address three core questions through our experiments: (1) How does GeoFlow perform compared to existing methods? To this end, we design a suite of comprehensive metrics to systematically evaluate GeoFlow against competitive baselines in Sec. 4.2. (2) Are the individual components of GeoFlow effective? We conduct targeted ablation studies to verify the contribution of each component in Sec. 4.3. (3) Is the design motivation of GeoFlow well-founded? We perform a series of analysis experiments to demonstrate that our observations and motivations are not merely intuitive but also supported by empirical evidence in Sec. 4.4.

Dataset. We evaluate three OD flow regimes, covering commuting, freight logistics, and air-travel mobility. For commuting flows, we use the CommutingODGen dataset, a large-scale open-source benchmark for OD flow-related tasks that contains commuting flows across diverse urban environments in the United States [6]. For freight logistics, we use the Freight Analysis Framework (FAF) database¹ released by the Bureau of Transportation Statistics in the U.S. Department of Transportation, which provides U.S. freight flows with annual estimates from 2018 to 2024. For air-travel mobility, we use the Tourism Volumes and Flows series² published by New Zealand’s Ministry of Business, Innovation and Employment, which reports destination-level visitor volumes and origin-specific visitor flows. More details on dataset construction are provided in Appendix 7.1.

Metrics. Since GeoFlow is evaluated on both OD flow prediction and generation, we adopt a set of metrics from multiple perspectives. (1) Numerical accuracy. For the prediction task, we employ Common Part of Commuting (CPC) and Normalized Root Mean Square Error (NRMSE) to quantify the alignment between predicted flows and ground truth values. (2) Distributional similarity. Following prior work [6], [10], we use Jensen–Shannon Divergence (JSD) to measure the discrepancy between the generated and true distributions of inflow, outflow, and OD flows. (3) Diversity. For the generation task, in addition to reconstruction accuracy and distributional similarity, we adopt the sample-level dissimilarity to assess diversity (Div.) across generated samples, which helps prevent potential mode collapse. Please refer to Appendix 7.2 for more detailed information on metrics.

4.2 Comparative Evaluation↩︎

Baselines. We compare GeoFlow against a broad and representative set of baselines: the classical gravity model; machine-learning regressors including RF, GBRT, and SVR; deep-learning approaches DGM, GMEL, and TransFlower; and recent graph-based methods NetGAN, DiffODGen, and WEDAN. More details and discussion about the baselines are provided in Appendix 12.

Results and analysis. Quantitative results are summarized in Tab. 1. Overall, GeoFlow achieves the best overall performance across both prediction and generation tasks on the three OD regimes. For prediction, deep models generally outperform classical baselines, and GeoFlow further improves over the strongest transformer baseline TransFlower, with notable gains in CPC and reductions in NRMSE and JSD, e.g., on Commuting (CPC \(0.604\) vs.\(0.486\), NRMSE \(1.727\) vs.\(2.146\), JSD \(0.171\) vs.\(0.281\)), Freight (CPC \(0.645\) vs.\(0.593\), NRMSE \(0.606\) vs.\(1.042\), JSD \(0.048\) vs.\(0.098\)), and Air Travel (CPC \(0.808\) vs.\(0.773\), NRMSE \(0.360\) vs.\(0.463\), JSD \(0.013\) vs.\(0.016\)). For generation, recent graph-based generative baselines DiffODGen and WEDAN improve upon NetGAN, and GeoFlow achieves the best overall quality–diversity trade-off, yielding the highest diversity (Div.\(0.069\)) and consistently lower NRMSE and JSD across three datasets. In particular, GeoFlow attains the best CPC on Commuting and Freight (CPC \(0.566\) and \(0.485\)), while on Air Travel it remains competitive in CPC (\(0.457\) vs.\(0.486\)) and substantially improves reconstruction and distribution matching (NRMSE \(1.213\) vs.\(1.485\), JSD \(0.031\) vs.\(0.036\)). These results suggest that explicitly incorporating geospatial structure with OD-aware decoding improves both predictive accuracy and generative fidelity. Further analyses on region size effects and generalizability are provided in Appendix 8.2.

4.3 Ablation Studies↩︎

We conduct ablation studies on the Commuting dataset under both prediction and generation settings, using CPC, NRMSE, and JSD for evaluation. Additional ablation results on the Freight dataset are provided in Appendix 8.1.

3.5pt

Ablation on geospatial feature augmentation. Tab. 2 shows that augmenting straight-line distance with relative coordinates, \(k\)-hop distance, and free-flow distance consistently improves generation quality over the SL-only baseline, and the full feature set achieves the best overall performance (ID 5), increasing generation CPC from \(0.462\) to \(0.527\) and reducing generation NRMSE/JSD from \(6.162/0.252\) to \(4.330/0.194\). These features complement each other by encoding orientation-aware location cues and network-based accessibility, leading to the strongest overall feature combination.

Ablation on encoder and decoder architecture. Replacing the proposed encoder with a conventional graph convolution module or replacing the axial-global attention decoder with global-only attention degrades performance, as shown in Tab. 2. The full encoder-decoder architecture achieves the best results (ID 9). For prediction, it improves CPC/NRMSE/JSD from \(0.567/2.134/0.201\) to \(0.597/1.730/0.171\) over ID 6; for generation, it improves the corresponding metrics from \(0.510/6.780/0.209\) to \(0.527/4.330/0.194\).

4.4 Analysis Experiments↩︎

Complementarity of straight-line and network-based distances. We compare straight-line (SL) distance with \(k\)-hop and free-flow (FF) distances using row-normalized matrices. Tab. ¿tbl:tab:32sim32anal? shows that SL and FF are relatively close in aggregate (Frobenius norm \(0.741\), MAE \(0.021\)), yet still exhibit noticeable worst-case deviations (MAD \(0.128\)). In contrast, \(k\)-hop differs substantially from both SL and FF (Frobenius norm \(\approx 8.0\), MAE \(\approx 0.15\)). This indicates that SL alone cannot capture network-induced distortions, motivating the inclusion of network-based descriptors in modeling geospatial dependencies.

Difficulty of inferring higher-order geospatial attributes. We test whether higher-order cues can be recovered from adjacency and SL distance by training an MLP to predict relative angles, \(k\)-hop, and FF distances, which shows large reconstruction errors, including a \(51.1\%\) average relative error for angles and substantial Frobenius-norm discrepancies for \(k\)-hop and FF. This suggests that shortest-path and angular structure are hard to infer implicitly, supporting our choice to provide these features explicitly via preprocessing.

5 Conclusion↩︎

In this work, we introduce GeoFlow, a novel framework for modeling origin–destination (OD) flows. GeoFlow augments critical yet often overlooked geospatial relationships as representations and employs a geometric–intrinsic fusion encoder coupled with an axial–global attention decoder, achieving state-of-the-art performance on both prediction and generation tasks. Extensive experiments validate the effectiveness of the proposed approach. Ablation and analytical studies further highlight the significance of geospatial attributes and the encoder–decoder design, offering new insights for OD flow modeling and related applications.

Acknowledgements↩︎

This work was sponsored by 2025 Shanghai Key Technology Research and Development Program, Next- Generation Information Technology Project under Grant No. 25511103800.

Impact Statement↩︎

This work aims to advance machine learning methods for origin-destination flow modeling. Such models may support urban planning, mobility analysis, and infrastructure assessment, but their outputs should not be used as the sole basis for high-stakes planning or policy decisions without domain validation. Like other black-box neural network models, GeoFlow may produce plausible but unreliable outputs under distribution shifts or unexpected input conditions, and its predictions cannot be fully explained by explicit physical mobility mechanisms. The attention-based case studies in this paper provide only partial qualitative interpretability and should not be interpreted as complete causal or physical explanations of mobility behavior. In addition, GeoFlow assumes access to reliable geographic signals, such as network connectivity and geodesic distances. When these signals are unavailable, outdated, or biased, the predicted or generated flows may inherit corresponding errors. Care should therefore be taken when applying the model in data-sparse regions or in settings involving sensitive mobility patterns.

Appendix↩︎

6 Method Details↩︎

6.1 Flow Matching Models for Origin-Destination Flow Generation↩︎

To generate high-fidelity origin-destination (OD) flow matrices, we adopt Flow Matching (FM), an efficient and expressive framework for generative modeling [26], [33]. In contrast to diffusion models, which typically require simulating stochastic or ordinary differential equations (ODEs) during training, FM directly learns the vector field that transports a simple prior distribution to the complex data distribution. This simulation-free paradigm substantially improves training efficiency. In the following, we use \({\boldsymbol{x}}\) to denote the vectorized form of the OD flow matrix \({\boldsymbol{F}}\), following the conventional notation in the flow matching literature.

6.1.1 Formulation↩︎

The core principle of FM is to define a continuous-time flow from a simple prior distribution \(p_0\), typically a standard normal distribution \(\mathcal{N}(\mathbf{0}, \mathbf{I})\), to the target data distribution \(p_1\), which in our case corresponds to OD matrices. This flow is governed by an ODE: \[\frac{d{\boldsymbol{x}}_t}{dt} = {\boldsymbol{v}}_t({\boldsymbol{x}}_t),\] where \(t \in [0, 1]\) is a pseudo-time variable, \({\boldsymbol{x}}_t\) denotes the state at time \(t\), and \({\boldsymbol{v}}_t\) is a time-dependent vector field. The goal is to train a neural network \({\boldsymbol{v}}_{\boldsymbol{\theta}}({\boldsymbol{x}}_t, t)\) to approximate this vector field.

To derive a tractable training objective, we first define a conditional probability path \(p_t({\boldsymbol{x}}| {\boldsymbol{x}}_1)\) that connects a prior sample \({\boldsymbol{x}}_0 \sim p_0\) to a data sample \(vx_1 \sim p_1\). Following [26], we adopt a simple linear interpolation path: \[{\boldsymbol{x}}_t = (1-t){\boldsymbol{x}}_0 + t{\boldsymbol{x}}_1,\] with the corresponding conditional vector field: \[{\boldsymbol{u}}_t({\boldsymbol{x}}_t | {\boldsymbol{x}}_0, {\boldsymbol{x}}_1) = {\boldsymbol{x}}_1 - {\boldsymbol{x}}_0.\] For this linear path, the target vector field is constant and independent of both \({\boldsymbol{x}}_t\) and \(t\).

The FM objective is to minimize the expected mean squared error between the predicted vector field \({\boldsymbol{v}}_{\boldsymbol{\theta}}\) and the target vector field \({\boldsymbol{u}}_t\): \[\mathcal{L}_{\mathrm{FM}} = \mathbb{E}_{t \sim \mathcal{U}(0,1), {\boldsymbol{x}}_0 \sim p_0, {\boldsymbol{x}}_1 \sim p_1} \left[ \| {\boldsymbol{v}}_{\boldsymbol{\theta}}((1-t){\boldsymbol{x}}_0 + t{\boldsymbol{x}}_1, t) - ({\boldsymbol{x}}_1 - {\boldsymbol{x}}_0) \|_2^2 \right].\]

6.1.2 Training↩︎

In practice, optimization proceeds by sampling mini-batches. For each training step, a data point \({\boldsymbol{x}}_1\) (an OD matrix), a noise sample \({\boldsymbol{x}}_0 \sim \mathcal{N}(\boldsymbol{0}, \boldsymbol{I})\), and a time step \(t \sim \mathcal{U}(0,1)\) are sampled. The interpolated state \({\boldsymbol{x}}_t\) and target velocity \({\boldsymbol{v}}_{\mathrm{tgt}} = {\boldsymbol{x}}_1 - {\boldsymbol{x}}_0\) are then computed. The model \({\boldsymbol{v}}_{\boldsymbol{\theta}}\) takes \({\boldsymbol{x}}_t\), \(t\), and associated context as input to predict a velocity \({\boldsymbol{v}}_{\mathrm{pred}}\), and the mean squared error between \({\boldsymbol{v}}_{\mathrm{pred}}\) and \({\boldsymbol{v}}_{\mathrm{tgt}}\) is minimized via backpropagation.

6.1.3 Sampling↩︎

Once trained, \({\boldsymbol{v}}_{\boldsymbol{\theta}}\) enables the generation of new OD matrices by solving the initial value problem defined by the learned ODE. Starting from a prior sample \({\boldsymbol{x}}_0 \sim \mathcal{N}(\boldsymbol{0}, \boldsymbol{I})\), we integrate the ODE from \(t=0\) to \(t=1\) using a numerical solver. In our implementation, the sampling process employs the forward Euler scheme: \[{\boldsymbol{x}}_{t+\Delta t} = {\boldsymbol{x}}_t + {\boldsymbol{v}}_{\boldsymbol{\theta}}({\boldsymbol{x}}_t, t)\Delta t,\] where \(\Delta t\) is the integration step size. After \(n_{\mathrm{steps}} = 1/\Delta t\) iterations, the resulting state \({\boldsymbol{x}}_1\) constitutes a novel sample drawn from the learned distribution. This procedure is executed independently for each graph in a batch, enabling efficient generation of multiple OD matrices.

6.2 Preprocessing↩︎

In this work, GeoFlow is trained and evaluated on the three different OD datasets, in which each region is represented as a graph in matrix form. The nodes correspond to geographic areas (e.g., census tracts), and the edges represent road network adjacency and Euclidean distance. For each region, the dataset provides node attributes, spatial coordinates, a road network adjacency matrix, a Euclidean distance matrix, and the ground-truth commuting OD flow matrix. A standardized preprocessing pipeline is applied to harmonize these inputs and derive multi-view relational structures.

The preprocessing for each urban region consists of the following steps:

1. Area Feature Scaling. For each region, we construct a raw feature vector by concatenating demographic attributes with counts of points of interest (POIs). To ensure comparability across features, we apply column-wise max normalization, dividing each feature value by the maximum observed within the same region. This rescales all features to the range \([0,1]\).

2. Coordinate Normalization. The raw geographic coordinates (zone centroids) are standardized to obtain a canonical representation invariant to translation and rotation. Specifically, we first subtract the mean to center the coordinates. We then apply Principal Component Analysis (PCA) [34] to align the coordinate axes with directions of maximal variance. Finally, the transformed coordinates are rescaled to the range \([-1,1]\) by dividing by the maximum absolute value along each axis.

3. OD Flow Normalization. The raw OD matrix, containing integer-valued flows, is normalized to form a probability distribution by dividing by the total commuting volume of the region. This total volume is stored and later used during inference to rescale generated probability matrices back to absolute flows.

4. Multi-view Affinity Construction. To capture diverse spatial relationships between zones, we construct three affinity matrices from different distance measures. Given a distance matrix \({\boldsymbol{D}}\), we compute an affinity matrix \({\boldsymbol{A}}\) using an exponential kernel: \[A_{i,j} = \exp(-D_{i,j}/\tau),\] where \(\tau\) is set to the median of all positive, finite distances in \(\mathbf{D}\). This produces a dense, weighted graph where affinity decreases with distance. The three distance notions are:

   • Euclidean Affinity: Based on straight-line distance between zone centroids, capturing pure spatial proximity.

   • Topological Affinity: Based on unweighted shortest-path distance (number of hops) on the road network, reflecting connectivity in the transport infrastructure.

   • Network-Geodesic Affinity: Based on shortest-path distance on the road network weighted by edge lengths, representing the most efficient travel distance along the network.

   For each affinity matrix, we also compute its symmetrically normalized form, \[{\boldsymbol{A}}_{\mathrm{norm}} = {\boldsymbol{D}}_{\mathrm{deg}}^{-1/2} {\boldsymbol{A}}{\boldsymbol{D}}_{\mathrm{deg}}^{-1/2},\] where \({\boldsymbol{D}}_{\mathrm{deg}}\) is the weighted generalized degree matrix with diagonal entries \(({\boldsymbol{D}}_{\mathrm{deg}})_{i,i}=\sum_j A_{i,j}\), i.e., the total connection strength of node \(i\). Symmetric normalization mitigates scale differences introduced by heterogeneous node degrees and yields affinity matrices that are numerically stable and comparable across regions. We include all normalized affinity matrices in the geometric relation descriptor \({\boldsymbol{G}}\), which supplies multi-relational contexts to downstream modules.

7 Experiment Details↩︎

7.1 Dataset Construction↩︎

The datasets cover distinct OD regimes, including commuting, freight logistics, and air travel. Tab. 3 summarizes the scale and sparsity of the three datasets.

For the Commuting dataset, we use the CommutingODGen benchmark [6], where each sample corresponds to a region-level OD matrix with provided area attributes, geographic relationships, and observed commuting flows.

The Freight dataset is built from FAF5 2022 freight flows. We aggregate OD freight tonnage between FAF zones and derive distances from FAF distance-band midpoints. We then generate 5,000 freight subregions, with region sizes sampled from a Beta\((2,3.21)\) distribution scaled to \([10,132]\), resulting in an average of approximately 56.9 areas per region. Each subregion starts from a FAF area sampled according to its total inbound and outbound freight volume, and is expanded through freight-flow neighbors weighted by tonnage. This procedure produces localized freight-flow subnetworks used as region-level OD samples.

The Air Travel dataset is constructed from the raw visitor-days OD table by grouping records according to temporal resolution, population segment, and date. For each group, we create one full origin–destination bipartite network and 21 additional random subnetworks by sampling 70% of origin markets and 70% of destinations, yielding 1,980 region-level samples in total. Since the data are intrinsically bipartite, only origin-to-destination entries are valid observations. Therefore, the observation mask is set to valid for the origin-destination entries and invalid for origin-origin, destination-destination, and destination-to-origin entries. OD flows are computed only on valid entries accordingly.

7.2 Metric Details↩︎

We provide the formal definitions and additional details of the evaluation metrics as follows.

Numerical Accuracy. To evaluate OD flow prediction and reconstruction, we employ three standard metrics:

• Common Part of Commuting (CPC). CPC measures the similarity between two OD matrices by computing the fraction of overlapping flow volume: \[\text{CPC}({\boldsymbol{X}}, {\boldsymbol{Y}}) = \frac{2 \sum_{i,j} \min(X_{i,j}, Y_{i,j})}{\sum_{i,j} X_{i,j} + \sum_{i,j} Y_{i,j}},\] where \({\boldsymbol{X}}\) and \({\boldsymbol{Y}}\) denote the predicted and ground-truth OD matrices, respectively. Higher values indicate better alignment.

• Root Mean Square Error (RMSE). RMSE quantifies the average magnitude of error in flow counts: \[\text{RMSE}({\boldsymbol{X}}, {\boldsymbol{Y}}) = \sqrt{\frac{1}{n^2} \sum_{i,j} (X_{i,j} - Y_{i,j})^2}.\]

• Normalized Root Mean Square Error (NRMSE). NRMSE normalizes RMSE by the standard deviation of the ground-truth matrix, enabling comparison across different scales: \[\text{NRMSE}({\boldsymbol{X}}, {\boldsymbol{Y}}) = \frac{\text{RMSE}({\boldsymbol{X}}, {\boldsymbol{Y}})}{\sigma_{{\boldsymbol{Y}}}}, \quad \sigma_{{\boldsymbol{Y}}} = \sqrt{\frac{1}{n^2} \sum_{i,j} (Y_{i,j} - \mu_{{\boldsymbol{Y}}})^2}, \quad \mu_{{\boldsymbol{Y}}} = \frac{1}{n^2} \sum_{i,j} Y_{i,j}.\]

• Mean Absolute Error (MAE). MAE provides a scale-independent measure of average deviation: \[\text{MAE}({\boldsymbol{X}}, {\boldsymbol{Y}}) = \frac{1}{n^2} \sum_{i,j} |X_{i,j} - Y_{i,j}|.\]

Distributional Similarity. To assess whether generated flows capture realistic distributional patterns, we compute the Jensen–Shannon divergence (JSD) between the distributions of inflows, outflows, and OD entries: \[\text{JSD}(P \| Q) = \frac{1}{2}\text{KL}(P \| M) + \frac{1}{2}\text{KL}(Q \| M),\] where \(M = \frac{1}{2}(P+Q)\) and \(\text{KL}(\cdot \| \cdot)\) denotes the Kullback–Leibler divergence. JSD is symmetric and bounded between 0 and 1, with lower values indicating closer distributions.

Diversity. For the generation task, we evaluate sample diversity using the metric defined as \[\label{supp32eq:32diversity} \text{Diversity}(\{\tilde{{\boldsymbol{F}}}^k\}_{k=1}^K) = \frac{2}{K(K-1)} \sum_{k < \ell} d(\tilde{{\boldsymbol{F}}}^k, \tilde{{\boldsymbol{F}}}^\ell),\tag{4}\] where \(\{\tilde{{\boldsymbol{F}}}^k\}_{k=1}^K\) are independently generated OD matrices under the same conditioning, and \(d(\cdot,\cdot)\) is a distance measure (e.g., in this work, RMSE). This metric computes the average pairwise dissimilarity among generated samples, with higher values indicating greater diversity. Intuitively, a model affected by mode collapse produces more similar or even identical samples with uniformly sampled random seeds, leading to high similarity (low dissimilarity and diversity), whereas a well-generalized model generates sufficiently distinct OD matrices, resulting in lower similarity and higher diversity.

7.3 Training Details↩︎

GeoFlow is trained end-to-end using the AdamW optimizer with an initial learning rate of \(1 \times 10^{-3}\) and a batch size of 64. The training objective is the Mean Squared Error (MSE) between the predicted and ground-truth normalized OD matrices. A ReduceLROnPlateau learning rate scheduler is employed, which decreases the learning rate by a factor of 0.5 if the validation metric does not improve for 10 consecutive epochs. To mitigate overfitting, early stopping is applied with a patience of 100 epochs. Model performance is evaluated on the validation set after each epoch, and the checkpoint achieving the best validation performance is retained. All experiments are implemented in PyTorch and conducted on a single NVIDIA GeForce RTX 4090 GPU.

7.4 Hyperparameters↩︎

The model architecture and training configuration are controlled by a set of specified hyperparameters. The hidden dimension for geometric and intrinsic embeddings is set to 32, which we found to provide a favorable balance between computational efficiency and predictive accuracy. The encoder consists of two multilayer perceptron (MLP) layers for encoding geometric attributes and fusing geometric and intrinsic attributes, followed by two graph attention layers. The decoder is composed of two layers of axial self-attention and two layers of global self-attention. All attention modules employ four heads, and dropout is disabled throughout training.

To ensure robustness, we conducted sensitivity analyses over different hidden dimensions, numbers of encoder/decoder layers, and attention heads. The reported configuration consistently achieved the best trade-off between performance and efficiency across validation runs, and is therefore adopted in the main experiments.

8 Additional Experiments↩︎

8.1 Additional Ablation Studies↩︎

3.5pt

Tab. 4 reports ablation results on the Freight dataset. For geospatial feature augmentation, the full feature set achieves the strongest overall balance: it obtains the best prediction CPC and JSD, the second-best prediction NRMSE, and the best generation JSD. Although the best generation CPC and NRMSE are achieved by partially ablated feature variants, the full feature set remains competitive on both metrics, suggesting that different geographic cues contribute complementary benefits under different evaluation criteria. For architecture ablation, the full encoder-decoder model achieves the best prediction performance across all metrics and the best generation JSD. The decoder-only variant obtains the best generation CPC and NRMSE, indicating that the axial-global attention decoder is particularly important for reconstruction quality, while combining the proposed encoder and decoder gives the strongest distributional fidelity. Overall, the Freight results support the usefulness of both geospatial feature augmentation and the proposed architecture, with some metric-specific trade-offs in the generation setting on the Freight dataset.

8.2 Generalizability↩︎

Beyond the main-text evaluations, we further probe the generalizability of GeoFlow along three complementary axes on the Commuting dataset [6]. First, in the prediction setting, we vary the spatial scale of training and evaluation regions to quantify how well the model transfers across omitted area scales. Second, in the generation setting, we construct longitudinal partitions of the study area and test whether a model trained on one part of the space can faithfully reproduce patterns in another. Finally, we perform a fine-grained performance analysis by conditioning on OD-flow magnitude and OD distance, which reveals how error is distributed across practically relevant regimes.

8.2.1 Generalizability Across Area Scales↩︎

To better understand how predictive performance and transferability vary with the spatial scale of regions, we adopt a scale-aware evaluation protocol that explicitly varies both the region sizes present during training and the region sizes on which models are evaluated.

Setup. For the experiments, regions are partitioned by the number of constituent areas into three size classes: small, medium, and large, as summarized in Tab. ¿tbl:supp32tab:32split32area32by32size?. We train distinct models under a set of controlled training compositions formed by varying the inclusion and exclusion of these size classes (for example, ID 4 is trained on the union of small- and medium-sized regions). Each trained model is evaluated separately on small, medium, large, and on the union of all regions (All). This design isolates how particular training mixtures influence performance both within and across scales and enables a principled comparison of generalizability under differing class imbalances and distributional shifts.

Analysis. When the evaluation scale is absent from the training pool (IDs 6, 11, and 16), performance on the held-out class deteriorates consistently. RMSE and MAE increase, and CPC declines, which confirms that excluding a scale impairs direct transfer. Crucially, however, the magnitude of this degradation is modest relative to a model trained on the full corpus (IDs 5, 9, and 13). In other words, although removing the target scale produces a measurable loss, the model nevertheless retains a substantial portion of its predictive ability across scales, indicating that it learns components of OD structure that generalize beyond the observed size class.

8.2.2 Generalizability Across Longitudinal Partitions↩︎

To complement the above scale-based study, we evaluate generalizability under a different kind of spatial shift in the generation setting. We partition the study areas into eastern and western subsets based on longitude and additionally construct a random partition with the same cardinality as a baseline.

Setup. For the “All - West” configuration, the generative model is trained on all regions except those in the western subset and then evaluated exclusively on the held-out West partition. “All - East” and “All - Random” are defined analogously. We report CPC, RMSE, and JSD between the generated and ground-truth OD distributions on the corresponding held-out partitions in Tab. ¿tbl:supp32tab:32generaliability322?.

Analysis. Across all three configurations, performance remains stable. CPC stays at a similar level across splits, RMSE differences are moderate, and JSD exhibits minor fluctuations. The random partition does not confer a clear advantage over structured East/West splits, suggesting that the model does not overfit to idiosyncratic geographic patterns of any particular longitudinal band. Instead, it captures generative regularities that transfer across distinct spatial subdomains.

8.3 Generalization Across City Morphologies↩︎

Setup. To examine whether GeoFlow generalizes across heterogeneous urban layouts, we conduct a morphology-based analysis on the Commuting dataset. We first characterize each region using a set of morphology-related descriptors that summarize its spatial scale, geometric configuration, network connectivity, distance distribution, and OD-flow structure. These descriptors capture broad differences in regional size, spatial anisotropy, compactness, graph density, and flow sparsity or concentration. After standardization, we apply K-means clustering to partition the regions into four morphology groups. Based on the aggregate profiles of the resulting clusters, we assign them interpretable names: Rural-Compact, Elongated-Sparse, Suburban, and Metropolitan. The morphology labels are used only to construct evaluation partitions and are not provided as additional supervision to the model.

Analysis. Tab. 5 reports cross-morphology transfer results, where each model is trained on a single morphology group and evaluated across all groups. Within the single-morphology models, the morphology-matched setting usually achieves the best or near-best performance on its corresponding target group, suggesting that shared spatial layout, graph structure, and OD sparsity patterns facilitate transfer. However, models trained on a single morphology often generalize poorly to structurally different groups, especially to Metropolitan regions. In contrast, training on the full dataset yields the best overall performance, indicating that exposure to diverse urban forms improves aggregate robustness. These results suggest that, by using region-centered, PCA-aligned, normalized relative coordinates together with network-based geographic descriptors, GeoFlow relies on transferable relative spatial relations rather than a fixed city-specific coordinate system.

8.3.1 Performance Analysis↩︎

Scale-dependent predictive difficulty. This partitioning makes explicit two salient characteristics of the dataset: (1) small and medium regions dominate numerically, whereas large regions are comparatively scarce; and (2) the intrinsic complexity of OD structure grows with region size. Consistent with these observations, as shown in Tab. ¿tbl:supp32tab:32generalizability?, predictive difficulty exhibits a stable ordering: models attain the best results on small regions, intermediate results on medium regions, and the worst results on large regions. Large-region evaluations display substantially higher RMSE and lower CPC, indicating both greater pointwise error and poorer fidelity of the learned OD distributions. This pattern is interpretable as the consequence of increased intra-region heterogeneity and multi-modal travel structure in larger spatial extents.

Performance profile across flow scales. Tab. ¿tbl:supp32tab:32analysis32on32flow? conditions performance on OD-flow magnitude. The distribution is highly skewed. More than \(83\%\) of OD pairs have flows below \(5\), whereas flows above \(100\) account for less than \(1\%\) of the pairs. As flow increases, absolute RMSE naturally grows, but normalized error (NRMSE) drops and CPC improves from \(0.105\) in the smallest bin to around \(0.62\)–\(0.64\) for flows in \([100, 1000)\). This indicates that, in relative terms, the model is substantially more reliable on high-volume OD pairs that dominate aggregate mobility patterns, while low-flow interactions, which are both noisier and less constrained by data, remain the most challenging. Moreover, when baseline flows are very small, CPC becomes highly sensitive to minor absolute fluctuations: a handful of misallocated trips can induce large apparent drops in correlation, making the low-flow bins look disproportionately poorly modeled. This is consistent with the interpretation that OD pairs with very low counts contain a substantial random component, and such randomness, combined with a small base volume, translates into amplified variability of distributional metrics. The slight CPC decrease in the \(\ge 1000\) bin is plausibly attributable to extreme data sparsity (only \(0.01\%\) of pairs, involving \(8.55\%\) of areas).

Performance profile across distance scales. We further examine performance as a function of normalized OD distance in Tab. ¿tbl:supp32tab:32analysis32on32distance?. Most OD pairs correspond to short- to medium-range trips: distances in \([0.2, 1)\) already account for over \(70\%\) of the mass. CPC is highest for the shortest-range bin \([0, 0.1)\) and remains relatively stable for distances up to \(1\), indicating that the model captures local and intra-region mobility patterns well. For very long-range interactions (\(\ge 2\)), CPC drops markedly to \(\approx 0.25\) despite a lower NRMSE, reflecting that these flows are extremely sparse (only \(0.05\%\) of OD pairs, involving \(39.82\%\) of areas) and concentrated along a small number of long-distance corridors. Overall, the distance-conditioned analysis shows that the model is most accurate in the distance regimes that dominate real-world demand, while rare long-distance trips and very short but noisy interactions remain more difficult to predict perfectly.

8.4 Model Capacity Scaling↩︎

To assess the scalability of our architecture, we systematically increase model capacity along two axes: width (hidden dimension) and depth (number of encoder/decoder attention layers and attention heads) with all other optimization hyperparameters and training schedules kept fixed.

Tab. ¿tbl:supp32tab:32scalability? reports the results. Increasing the hidden dimension from 16 to 64 (IDs 1–3) yields consistent gains. CPC improves from \(0.576\) to \(0.593\), RMSE drops from \(107.5\) to \(92.43\), and JSD decreases from \(0.212\) to \(0.194\). This indicates that the model can effectively leverage additional representational capacity, and that the architecture behaves in a stable, roughly monotone way under width scaling. Adding moderate depth and more heads on top of a 32-dimensional backbone (ID 4) further improves all metrics, reaching higher CPC and lower RMSE, and lower JSD. The incremental improvements suggest diminishing returns at higher capacities, consistent with a regime where performance is increasingly data-limited rather than model-limited. Overall, these results show that our model scales gracefully. Larger configurations provide better accuracy, but even the smallest variant (ID 1) attains a substantial fraction of the performance of the largest one. This makes it possible to trade off compute and memory for accuracy depending on deployment constraints, while retaining the same architectural design and training procedure.

8.5 Scalability Profiling↩︎

We further profile the computational scalability of GeoFlow with respect to the structural scale of a region. Here, structural scale refers to the number of areas in a region, and therefore to the number of OD pairs to be modeled. This is distinct from the numerical scale of OD flows, which is mitigated by the OD-flow normalization step in preprocessing. For each scale group, we report the average per-step wall-clock latency for both prediction and generation, including training and inference.

As shown in Tab. 6, the latency remains nearly flat for most regions below 200 areas. Prediction training stays around 18–19 ms per step, while prediction inference remains around 4 ms. Generation training has a similar cost to prediction training, since flow matching uses a single sampled time step during each training update. In contrast, generation inference is substantially slower because sampling requires multiple numerical integration steps rather than a single forward pass. The latency increases moderately for regions with 200–500 areas and rises sharply for very large regions with more than 500 areas. This pattern indicates that GeoFlow is computationally efficient for the majority of regions in the benchmark, while very large regional graphs remain the main computational bottleneck due to the rapid growth of OD-pair interactions with the number of areas.

8.6 Comparison of Diffusion Models and Flow Matching Models↩︎

To further assess the effectiveness of our framework, we conduct a comparative study between flow-matching and diffusion models. In this setting, only the flow-matching components are replaced with their diffusion counterparts, while all other configurations remain unchanged. Since diffusion models follow a different generative paradigm and training procedure, we perform extensive hyperparameter tuning to ensure comparable model size and computational cost. As shown in Fig. 3, although the diffusion model initially exhibits normal metric progression, it undergoes severe oscillations and sustained degradation. This instability arises from the inherently slower and noisier training dynamics of diffusion, which hinder efficient convergence. In contrast, the flow-matching model exhibits consistent and monotonic improvements throughout training, indicating its stability and efficiency in capturing the underlying OD flow structure.

9 Computational Complexity and Resource Usage↩︎

We profile the computational footprint of our proposed GeoFlow architecture against the transformer-based baseline TransFlower. Table ¿tbl:supp32tab:32analysis32on32complexity? reports average training and inference wall-clock time, average and peak FLOPs, and peak GPU memory usage, computed using the same configuration as that used to obtain its main performance results reported in the main experiments. All measurements are taken from steady-state runs after warm-up, so the numbers reflect typical rather than best-case performance.

GeoFlow achieves noticeably shorter wall-clock times than TransFlower: training is roughly \(3.6\times\) faster (\(0.0160\) s vs.\(0.0580\) s) and inference is about \(4\times\) faster (\(0.1104\) s vs.\(0.4468\) s). At the same time, GeoFlow requires around \(4\times\) higher average FLOPs (\(6.12\times10^{8}\) vs.\(1.49\times10^{8}\)) and a corresponding increase in peak FLOPs (\(3.58\times10^{11}\) vs.\(8.82\times10^{10}\)), while peak memory consumption remains comparable between the two models. This contrast highlights that GeoFlow adopts a more compute-intensive yet highly parallelizable design that leads to better hardware utilization and throughput, whereas TransFlower offers lower FLOP counts at the cost of higher latency. Taken together with the accuracy improvements of GeoFlow, these results indicate that GeoFlow provides a favorable accuracy–efficiency trade-off, while remaining tractable at city scale on a single general-purpose GPU.

10 Recommendations for training-set composition↩︎

The generalizability experiments in Appendix 8.2.1 suggest several practical guidelines for choosing which size classes to include when constructing the training set. The selection of region sizes has systematic and interpretable consequences for both within-scale accuracy and cross-scale transfer.

Models trained on the combination of small and medium areas tend to produce the strongest results for small-scale and aggregate evaluations (IDs 4, 8, 12 and 16 in Tab. ¿tbl:supp32tab:32generalizability?) and, in aggregate metrics, can even outperform the model trained on the full corpus (ID 4 vs.ID 1). This apparent advantage is primarily attributable to the numerical predominance of small and medium areas: their abundance supplies stronger and more stable gradient signals during optimization, which reduces average error on the dominant scales and therefore improves aggregate statistics. However, the same training regimen is weaker on large-area evaluation than the model trained on all scales (ID 16 vs.ID 13). Including the relatively scarce but distributionally distinct large-area examples in the training set increases the model’s ability to capture large-area specific structure and yields better performance on large areas and, importantly, stronger overall generalizability across all scales.

By contrast, training on medium and large areas is most effective when the evaluation target is medium or large (IDs 10 and 14), which suggests that medium areas provide transitional patterns that facilitate transfer toward large regimes. Naively pairing small and large areas without an intermediate representation often degrades performance (IDs 7, 11 and 15), likely because the omission of medium examples amplifies the domain shift between the two extremes.

In summary, the results point to a practical trade-off. Emphasizing small and medium areas during training tends to improve aggregate metrics and accuracy on the majority of areas in our corpus, whereas incorporating the relatively scarce large-area examples improves fidelity on large areas and enhances cross-scale robustness. This pattern appears consistently in our dataset but should be interpreted with caution rather than as a universal prescription. The balance between aggregate performance and coverage of data-scarce large areas is inherently dataset dependent and therefore calls for pragmatic, case-by-case choices when assembling training data in practical applications.

11 Qualitative Evaluation↩︎

11.1 Case Study↩︎

To validate that GeoFlow captures meaningful spatial dependencies rather than merely fitting statistical correlations, we visualize the learned attention weights for a representative region, Barbour County, Alabama, United States. We examine the spatial interaction of both global and axial attention mechanisms.

In the second panel of Fig. 4, the global attention heatmap shows, for each pair of areas, how much the area in column \(j\) contributes to the representation of the area in row \(i\). A prominent vertical band of high intensity appears at area 6, which corresponds to Eufaula, the largest city and economic center of Barbour County. This indicates that many other areas place relatively high attention on area 6, consistent with a hub role in the local mobility system. In contrast, area 3, which is located near the county boundary, attracts noticeably less attention from the rest of the county. The weaker column associated with area 3 reflects its more peripheral role and suggests that changes in this edge area have a smaller impact on the representations of other areas.

At a finer scale, the axial attention heatmaps in the third and fourth panels of Fig. 4 show the distribution of attention over destinations for two different origins. The third panel corresponds to an urban center, origin area 6, and exhibits a relatively diffuse pattern with a maximum attention value of about \(0.137\) spread across several neighboring destinations such as indices 5, 7 and 8. This suggests that flows originating from the center are influenced by a combination of multiple nearby areas, and the model distributes attention over a broader neighborhood that matches the multifaceted role of the city core. The fourth panel corresponds to a peripheral origin, area 0, and shows a sharper and more self focused pattern, with a maximum attention value of about \(0.148\) concentrated on index 0. For this edge area, outflow is therefore governed mainly by the intrinsic characteristics of the origin rather than by strong interactions with alternative destinations. The contrast between diffuse urban attention and focused peripheral attention indicates that GeoFlow reflects the spatial heterogeneity of commuting flows and adjusts its effective receptive field to the functional role of each area, rather than relying only on distance as in classical gravity models.

11.2 Visualization↩︎

To qualitatively assess the performance of GeoFlow, we randomly select several regions for visualization. As illustrated in Fig. 5, for each region we provide its geographic map, the ground truth OD flow, the model prediction, and three independently generated samples. The prediction results are designed to closely approximate the most likely OD flows, while the generated samples reflect the overall flow patterns while preserving key structural characteristics. This comparison demonstrates that GeoFlow not only accurately predicts dominant flows but also effectively captures the variability and essential features of the underlying OD distribution.

12 Discussion of Baseline Methods↩︎

To provide a comprehensive evaluation of GeoFlow, we compare against a representative set of baselines spanning principle-driven, machine learning, and deep learning paradigms. In this section, we summarize their main characteristics, strengths, and limitations.

Gravity Model. The gravity model is among the earliest and most widely used approaches to OD flow modeling. It posits that flows are positively associated with the "mass" (e.g., population or economic activity) of the origin and destination and inversely associated with the distance between them, typically through a power-law or exponential decay. Its transparent mathematical form provides strong interpretability and has long served as a benchmark in mobility studies. However, its reliance on simple parametric functions limits flexibility in capturing heterogeneous and nonlinear mobility patterns observed in real data.

Classical Machine Learning Regressors. Tree-based models such as Random Forest (RF) and Gradient Boosted Regression Trees (GBRT), as well as kernel-based methods like Support Vector Regression (SVR), represent the next stage in OD flow modeling. These approaches directly learn mappings from urban attributes to flow values. RF and GBRT benefit from their ability to capture nonlinear feature interactions and to provide feature importance measures, while SVR leverages kernel functions to model similarity in high-dimensional feature spaces. Although they improve predictive accuracy over principle-driven models, their reliance on hand-crafted features and limited ability to capture spatial dependencies constrain their performance on large, complex urban systems.

Deep Learning Models. The Deep Gravity Model (DGM) builds on the intuition of gravity models but employs neural networks to learn flexible production, attraction, and impedance functions. This extension enables DGM to capture nonlinearities beyond traditional formulations while maintaining interpretability to some degree. The Geographical Multi-task Embedding Learning model (GMEL) incorporates graph attention networks to integrate spatial and structural signals, offering improved representation learning. However, these models generally treat areas as independent entities with limited consideration of long-range spatial interactions, which can restrict their expressiveness.

Graph-based Generative Models. Recent advances have explored graph neural networks and generative architectures to model OD flows as attributed graphs. NetGAN employs random walks to learn latent graph structures and can synthesize plausible OD matrices, but struggles with capturing fine-grained attributes. DiffODGen introduces diffusion processes to jointly model topology and edge weights, while WEDAN applies diffusion-based denoising to learn conditional flow generation. These methods capture complex structural dependencies and achieve strong generative performance, yet they typically abstract areas as nodes without explicit treatment of spatial proximity, limiting their ability to represent geographic relations that influence mobility.

Transformer-based Methods. TransFlower adapts transformer architectures to OD flow modeling by introducing relative positional embeddings, highlighting the importance of spatial relations between origins and destinations. While it demonstrates strong performance, its reliance on pairwise attention incurs high computational cost and can be difficult to scale to large urban systems. Furthermore, the integration of topological information remains limited, leaving room for improvement in modeling multi-relational geographic contexts.

In summary, principle-driven models emphasize interpretability but are constrained in flexibility; machine learning regressors improve predictive power but fail to capture structural dependencies; deep learning models introduce nonlinear representation learning but often overlook geographic context; and graph generative and transformer-based methods achieve strong expressiveness at the expense of interpretability or efficiency. GeoFlow is designed to address these gaps by explicitly incorporating multi-relational geographic structures while maintaining scalability, enabling more faithful modeling of OD flows across diverse urban contexts.

References↩︎

---

1. www.bts.gov/faf↩︎

2. www.mbie.govt.nz/immigration-and-tourism/tourism-research-and-data/tourism-data-releases↩︎