Complexity-aware Large Scale Origin-Destination Network Generation via Diffusion Model

As shown in Figure 1, our method includes two stages, which can be treated as separating the city-wide OD matrix generation task into two steps. In the first stage, as we can see in the upper part of Figure 1, the method determines whether there is a mobility flow between regions $P_{\theta_{M}}(M|X)$, i.e., topology structure, with the topology diffusion model. Then, the region pairs without flow will be directly set to zero flow and discarded in the rest of the process. According to the sparsity of the OD matrix, this will greatly reduce the scale of city-wide OD matrix generation and thus solve the first challenge. In the second stage, we use a continuous denoising diffusion model to learn the joint probability distribution of OD flows between the remaining region pairs $P_{\theta_{F}}(F|X)$ for generation.

The diffusion models in both the two stages are in the form of classifier-free conditional diffusion models (Ho and Salimans, 2022), where pairwise regional attributes are element-wise assigned to the OD flows between each two regions in the OD matrix. This modeling can make the generated OD matrix not only globally conform to the overall properties, but also retain the advantages of the pairwise predictive model by fully considering the features of OD pairs so that OD flows are also locally reasonable.

A noteworthy issue is that combining the two diffusion models directly has the problem of cascading errors, i.e., the errors generated in the topology generation stage are completely transferred to the flow generation stage, which leads to an amplification of the overall instability. To solve this issue, we adopt the strategy of teacher-force learning to collaborate the training of networks $\theta_{M}$ and $\theta_{F}$. The training process is based on data from the source city, where the true value of the OD matrix is known. Not the adjacency matrix $\hat{M}$ generated is applied for training the network parameters in the flow diffusion model, but the union $\widetilde{M}$ of generated $\hat{M}$ and the ground truth $M$ instead, which is computed as follow,

where the $\vee$ denotes the logic operator of OR. $\widetilde{M}$ then serves as a mask in the flow diffusion model, so that it focuses only on the elements in the OD matrix whose value are not zero. In this way, the flow diffusion model has the capability of tolerating errors arising in the topology diffusion model and remedy them by generating the missing zero flows during the topology generation, thus improving the overall robustness of the method.

To avoid destroying the sparsity of the generated OD matrix, we utilize the reverse denoising process with discrete state spaces to generate the adjacency matrix in the topology diffusion model. We adopt the d3pm proposed by Austin et al (Austin et al., 2021), which works for graph generation task (Vignac et al., 2022), to approximate the data distribution $p_{M}(M)$ with $p_{\theta_{M}}(M^{0})$. In our approach, data $M_{ij}\in\mathbb{R}^{d}$, which follows a binary distribution, is encoded with one-hot encoding, where $d$ equals two, i.e., the number of classes (zero/non-zero flow). For the convenience of representation, we briefly refer to the element in adjacency matrix $M_{ij}$ as $m$ in the following. The noise-addition in the diffusion process is implemented through the matrix product with a series of transition matrix $(\mathbf{Q}^{1},...,\mathbf{Q}^{T})$,

such that $\mathbf{Q}^{t}_{ij}$ denotes the probability of state transiting from class $i$ to class $j$. It is worth noting that the noisy data state $m^{t}$ at $t$ diffusion step can be directly obtained by closed-form (Austin et al., 2021) calculations,

where $\bar{\mathbf{Q}}^{t}=\mathbf{Q}^{1}\mathbf{Q}^{2}...\mathbf{Q}^{t}$. The transition matrix $\textbf{Q}^{t}$ applied in our method introduce the uniform noise perturbation, which is calculated as follows,

where $\top$ means the matrix transposition and $\alpha$ controls the degree of adding noise. Given $m^{0}$, the reverse denoising process can also be computed in closed-form based on Bayes rule (Austin et al., 2021),

where $\odot$ means the element-wise product. Since $m^{0}$ is unknown, we construct neural networks to fit the probability of $p_{\theta}(\hat{m}^{0}|M^{t})$ given noisy topology matrix $M^{t}$ and realize the denoising process with predicted $\hat{m}^{0}$ when doing generation,

For generating the adjacency matrix based on specific city characteristics, we then extend the above denoising process into a conditional denoising process to model conditional probability $p_{\theta_{M}}(M^{0}|\mathbf{X}_{r\in\mathcal{R}})$. The conditional denoising process is as follows,

In the procedure of generation, a pure noise $M^{T}\in\mathbb{R}^{N\times N\times d}$, which follows the uniform distribution, is first sampled. The data samples that follow the original data distribution are then recovered gradually from the $M^{T}$ by the iterative calculation of reverse denoising process described by Eq. 9, where $\theta_{M}$ is the trained neural networks. It should be emphasized that the whole process of recovering data through latent states is performed on a discrete space, which can fully preserve the sparsity of the OD matrix (Vignac et al., 2022).

For better considering the adjacency matrix, i.e., topology structure and fully modeling the network properties (Saberi et al., 2017), the network parameterization of topology diffusion model is designed to build with the graph transformer proposed by Dwivedi et al (Dwivedi and Bresson, 2020). Moreover, we make two improvements to the network, adapting it to classifier-free element-wise conditional schema and integrating a node properties augmentation module to enhance the ability of modeling the network properties of mobility flows (Saberi et al., 2017).

The backbone of the graph transformer receives two inputs, node features and edge features, as shown in Figure 2. The graph transformer consists of several layers. In each layer, all nodes form a position-insensitive sequence. The attention weights between each pair of nodes will be calculated through self-attention. It is worth noting that the calculation of attention weights used for information propagation between nodes also incorporates the influence of edge features. Each node updates its own features by weighted aggregation of neighbor information. Each edge updates itself by fusing its own features with the attention information between its end nodes. The calculation process of each layer in the graph transformer is as follows,

where $h_{i}^{l}$ means the hidden states of node $i$ after $l$ layers, $e_{ij}^{l}$ means the hidden states of the edge from node $i$ to node $j$ after $l$ layers, $O_{h}^{l}$, $O_{e}^{l}$, $Q^{k,l}$, $K^{k,l}$, $V^{k,l}$ and $W^{k,l}$ are the learnable parameters, $\|$ denotes the concatenation, $k=1$ to $K$ means the number of attention heads. The output of each layer will be input to the next layer. After the last layer, the output of node features $H^{L}$ is dropped and the output of edge features $E^{L}$ is converted into a tensor of the same shape as $M^{t}$ by a fully-connected layer to predict $\hat{M}^{0}$, as shown in Figure 2.

Message-passing graph neural networks have a limitation of unable to capture some specific features (Xu et al., 2018; Morris et al., 2019), such as the weighted aggregation schema can hardly model the degree of nodes in a graph. However, the network properties of mobility flow are largely reflected in the region level, i.e., node level in mobility flow network (Saberi et al., 2017). Therefore, we augment the capability of our graph transformer with a node properties augmentation to model the node level flow characteristics by statistically counting the properties of nodes into node features based on the adjacency matrix. The statistic properties of nodes include, but are not limited to, the degree of the node and the centrality, etc. With the help of node properties augmentation, the network can comprehensively model the noisy topology structure from both node and edge level in each diffusion step.

As shown in Figure 2, we make an improvement by introduce a CondEmb module on the neural network backbone to make it classifier-free element-wise conditional denoising networks, so that diffusion models can determine whether there is a mobility flow based on its characteristics of origin and destination. First, urban attributes for each region $\{\mathbf{X}_{r}|r\in\mathcal{R}\}$ are processed by MLPs to become regional condition embeddings, which have the same dimension as the node features that will be input to the graph transformer, as shown in Figure 2. Then, regional condition embeddings are fused with node features and edge features respectively through cross attention $\otimes$ as the final input of the transformer, where the node features are directly fused with the embeddings region by region, while the edge features are fused with the concatenation of the embeddings of origin and destination as well as the distance between them. In this way, each row and column in the adjacency matrix, and each element has corresponding spatial urban features as its unique conditions that affect its generation in the denoising process.

The network $\theta_{M}$ introduced above is trained by reducing the cross entropy between the prediction $\hat{M}_{ij}^{0}$ and true value $M_{ij}^{0}$ of each element to empower the denoising process represented by Eq. 9. The overall cross entropy loss of the neural networks is calculated as follows,

where LCE means cross entropy loss. The training algorithm of networks is shown as Algorithm 1. Once neural networks are trained, they can be employed to generate topology structure of mobility flow networks for better city-wide OD matrix generation.

Different from topology generation, we adopt Gaussian noise to construct the diffusion process since OD flows are continuous values. Since DDPM (Ho et al., 2020) was proposed, diffusion models in continuous latent spaces have been explored to a significant extent (Nichol and Dhariwal, 2021). In this part, the process of noise-addition diffusion of our diffusion models is performed on the image-like tensor, i.e., the OD matrix $F\in\mathcal{R}^{N\times N\times d}$, where $d$ equals 1. Differently, the reverse denoising process is performed only at the position, where the elements equal to 1, in the mask $\widetilde{M}$, as shown in Figure1. Thanks to the sparsity of the OD matrix, although the OD matrix is large, our denoising process only focuses on the nonzero flow part, so we can effectively model it.

The diffusion process of the flow diffusion model is similar to the classical continuous diffusion models and has no parameters, which is shown in Eq. 1 from right to left. So we will not cover it repeatedly, but rather introduce the denoising process in detail. The computation of each denoising step is shown as follows,

where $\mathcal{N}$ means the Gaussian distribution, $\alpha$ is the noise schedule and $\bar{\alpha}^{t}=\alpha^{1}\alpha^{2}...\alpha^{t}$.

Similar to the topology diffusion model, we construct a conditional diffusion model by applying city characteristics as conditions to generate OD flows between regions. The network parameterization is also based on the graph transformer, except that we use the geo-contextual embedding learning methods proposed in GMEL (Liu et al., 2020) as the CondEmb module shown in Figure 2 to enhance our condition features representational capabilities. To be specific, urban attributes $\{\mathbf{X}_{r}|r\in\mathcal{R}\}$ for each region are processed by two GATs (Veličković et al., 2017) to obtain the embeddings for each region as the origin and the destination, respectively. Then, region embeddings that are fused with the node features input into the graph transformer are composed of the origin embedding and the destination embedding of the region. The edge embeddings consist of the origin embedding of the origin and the destination embedding of the destination, as well as the distance between them. In addition, node properties augmentation is also used for flow generation to enhance the graph modeling capability from both node and edge level.

Ho et al. (Ho et al., 2020) show that the denoising process can be trained by predicting the added noise in the diffusion process at the corresponding time step, i.e., optimizing the neural networks through the following objective,

where $\mathbf{\epsilon}$ denotes the noise added during the step $t$ in diffusion process and $||\cdot||_{2}$ means the $L-2$ norm. The training algorithm is referred to Algorithm 2. The $\widetilde{M}$ is used in collaborative training schema to reduce the effect of cascade errors of combining the two stages, whcih has been introduced in Sec. 4.1.2. Once the neural network is trained, it can be used to generate the OD flows for given topology of mobility flow network. First, a pure noise is sampled from the standard Gaussian noise $\mathcal{N}(0,\mathbf{I})$. The noise $\mathbf{\epsilon}_{\theta_{F}}$ is then sequentially predicted and removed from $F^{t}$ to obtain $F^{t-1}$. $F^{0}$ is finally generated via $T$ denoising step iteratively. It is worth emphasizing that elements with a mask of 0 in $\widetilde{M}$ are not involved in the training and generation process because of the mask $\widetilde{M}$.

We choose three major US metropolises for our experiments. The cities are divided into regions at the census tract level. The basic statistic information on the cities is presented in Table 1. The data for each city includes two parts, city characteristics and the OD matrix, which are introduced as follows.

City Characteristics. City characteristics depict the geographic structure of an entire city and the spatially heterogeneous distribution of city functions. In detail, each region carries attributes about demographics and urban functions, which are portrayed through American Community Survey Data collected by the U.S. Census Bureau and POIs distribution from OpenStreetMap (OpenStreetMap contributors, 2017) respectively. Each region aggregates all the POIs therein to obtain a distribution over different categories. The distance between regions is determined in terms of the planar Euclidean distance between the centroids of the regions. Based on the above data processing, we extract features for each region, as a feature vector including demographic information by age, gender, income, etc. and the number of different categories of POIs, which are normally used for OD matrix generation as in previous works (Robinson and Dilkina, 2018; Pourebrahim et al., 2019; Liu et al., 2020; Simini et al., 2021).

Commuting OD Matrix. The OD matrices are constructed based on the commuting flow collected in Longitudinal Employer-Household Dynamics Origin-Destination Employment Statistics (LODES) dataset in 2018. The commuting flows are aggregated into the regions. Each element in the OD matrices records the number of workers who are resident in one region and employed in another.

We compare our proposed DiffODGen with the following six baselines that can be categorized into three categories.

The first category contains two traditional methods.

Gravity Model (GM) (Barbosa et al., 2018). Motivated by Newton’s law of Gravitation, GM assumes that the mobility flow is positively related to the populations of origin and destination, and negatively related to the distance between them.

Random Forest (RF) (Robinson and Dilkina, 2018; Pourebrahim et al., 2019). RF is a widely adopted tree-based machine learning model thanks to its robustness, achieving rather competitive performance in the task of OD flow generation.

The second category contains two state-of-the-art (SOTA) deep-learning methods developed for OD flow generation.

Deep Gravity Model (DGM) (Simini et al., 2021). DGM introduces deep neural networks into the modeling of gravity models to obtain flows by predicting the probability distribution to different regions when the outflow is given. We adapt the model to predict volumes of OD flow directly, making it to generate the OD matrix for new cities.

Geo-contextual Multitask Embedding Layer (GMEL) (Liu et al., 2020). GMEL uses graph neural networks (GNNs) to aggregate information from neighbors of each region (node) to capture its spatial features in a city (graph), which helps learn better-quality regional embeddings to achieve better prediction accuracy.

Finally, we additionally compare two deep generative models to check the validity of our design in terms of generating a realistic mobility flow network, i.e., the OD matrix.

NetGAN (Bojchevski et al., 2018). NetGAN is a GAN-based model that mimics real networks by generating random walking sequences with the same distribution as those sampled from the real networks. We adapt it to generate directed weighted networks, i.e., OD matrices.

Denoising Diffusion Probabilistic Model (DDPM) (Ho et al., 2020). We also adopt the naive diffusion model that directly learns to generate the OD matrix as a baseline. In this baseline, no distinguishing treatment is assigned to zero elements and non-zero elements in the OD matrix.

We use two kinds of evaluation metrics to verify the generation realism, including three metrics related to flow value error and three metrics related to network statistics similarity. Specifically, the three error-based metrics are Root Mean Square Error (RMSE), Normalized Root Mean Square Error (NRMSE) and the commonly adopted Common Part of Commuting (CPC).

where the $\bar{\mathbf{F}}$ denotes the mean. To calculate the statistical similarity between generated OD matrices and real OD matrices, we adopt Jensen-Shannon Divergence (JSD) and use it to measure the distance between distributions of generated data and real data, with respect to three typical network statistics, i.e., inflow, outflow and OD flow. The computation of JSD is shown as follow,

where KL means Kullback–Leibler divergence and $\mathbf{P}$ denotes the empirical probability distribution. The inflow and outflow are calculated by summing all flows into and out of regions, respectively.

The number of layers of graph transformer in the topology and flow generation phases are 2 and 3 respectively, while the numbers of channels are both 64. The diffusion steps of two diffusion models in our method are 1000, and both use a cosine noise scheduler proposed by Nichol et al (Nichol and Dhariwal, 2021). The denoising networks in two diffusion models are trained with Adam optimizer (Kingma and Ba, 2014) with a learning rate of 3e-4.

The gravity model follows the way in (Barbosa et al., 2018) with 4 fitting parameters. The $n\_estimators$ of random forest are 100. The number of layers for DGM is 10. The number of layers in GNN-based models is 3 and the number of channels is 64. The DDPM baseline follow the parameter settings with the flow diffusion model in our method.