Deep-Ensemble-Based Uncertainty Quantification in
Spatiotemporal Graph Neural Networks for Traffic Forecasting

The first step of our approach is the adaptation of SQR in DCRNN to predict the quantiles of the forecasting traffic distribution.

Consider a typical regression setting consisting of input features $X$ and output variable $Y$ taking real values $y$. Let $F(y)=P(Y\leq y)$ be the cumulative distribution function of $Y$; the $q^{th}$ quantile of $Y$ can be calculated by the quantile distribution function $F^{-1}(q)=inf(y:F(y)\geq q)$, where $q\in(0,1)$. The quantile regression approach seeks to model $\hat{y}=\hat{f}_{q}(x)$ that approximates the conditional quantile distribution function $y=F^{-1}(q\mid X=x)$. Specifically, the quantile regression estimates the value of $Y$ at a given quantile $q$ conditioned on the given $x$ values of input features $X$. This estimation can be achieved by using the pinball loss function [15, 16]:

Building on this, the recently proposed SQR approach estimates all the quantiles simultaneously by solving the following optimization problem:

This optimization problem is solved by using minibatch stochastic gradient descent, where random quantile levels $q\sim U[0,1]$ (uniform distribution) are sampled afresh for each training point and for each minibatch.

In DCRNN, for each sensor location and for each forecast time step, we seek to estimate the traffic metric at all quantiles conditioned on the past $K^{\prime}$ traffic metrics defined over $N$ sensor locations. This is achieved by replacing the MAE loss minimization with the SQR minimization. Applying Eq. 5 in Eq. 3, we get

The key advantages of SQR include well-calibrated prediction intervals; joint estimation of the quantiles that can eliminate the unacceptable crossing quantiles [17]; and the ability to model non-Gaussian, skewed, asymmetric, multimodal, and heteroskedastic noise in the data.

The second step of our approach is the use of a scalable hyperparameter search for tuning the hyperparameters of DCRNN.

DCRNN has several training hyperparameters, including batch size, a threshold max_grad_norm to clip the gradient norm to avoid exploring gradient problem of RNN [18], initial learning rate and learning rate decay, maximum diffusion steps, and architecture configurations (network topology parameters) such as number of encoder and decoder layers, number of RNN units per layers, and filter type (i.e., random walk, Laplacian). The hyperparameters are crucial in determining DCRNN’s forecasting accuracy.

To tune the hyperparameters of the DCRNN model with SQR, we use the DeepHyper [19] open source software that has a scalable Bayesian optimization (BO) method. Given the training and validation data, a neural network model, and the feasible set of values of the training hyperparameters, DeepHyper uses an asynchronous BO method seeks to minimize the validation error by repeatedly training the neural network model with various hyperparameter values. The DeepHyper BO uses an incrementally updated random forest model to learn the relationship between the hyperparameter values and their corresponding validation errors. The BO method uses the random forest model to predict the validation errors of the evaluated hyperparameter configurations along with the uncertainty associated with the predicted value. The BO search involves both evaluation of hyperparameter configurations, where the random forest model is most uncertain, and exploitation, focusing on hyperparameter configurations that are closer to the previously found high-performing configurations. The balance between exploration and exploitation is achieved by selecting an unevaluated hyperparameter configuration $\boldsymbol{x}^{*}$ that minimizes the lower confidence bound ($\mathrm{LCB}$) acquisition function:

where $\hat{\mu}(\boldsymbol{x})$ and $\hat{\sigma}(\boldsymbol{x})$ are respectively the predicted estimate and the standard deviation of that estimate from the random forest model of the given unevaluated hyperparameter configuration $\boldsymbol{x}$ sampled from the hyperparameter space $\mathcal{D}$ and where $\lambda$ is a user-defined parameter that controls the trade-off between the exploration and exploitation.

A key reason for choosing DeepHyper for the DCRNN hyperparameter tuning is scalability. The DeepHyper BO follows a single-manager multiple-workers parallization scheme, where the manager runs the search, generates hyperparameter configurations, and sends them to workers for evaluation; each worker is responsible for evaluating a hyperparameter configuration by training the neural network model with the given configuration, computing the validation error, and returning that value to the manager. The manager generates hyperparameter configurations for evaluation in an asynchronous manner: as soon as one or more workers returns the validation error, the manager uses the returned values to retrain the surrogate model within the BO and generates the next hyperparameter configuration for the worker. To enable simultaneous multiworker evaluations, the BO uses a constant liar scheme, where configurations are generated sequentially by refitting the random forest model with a lie for the validation error, given by a predicted value from the model itself, for the previously selected unevaluated hyperparameter configuration.

We use DeepHyper to tune the hyperparameters of DCRNN not only to find the best hyperparameter configuration with minimal validation error but also to extract a set of hyperparameter values with minimal validation errors for generative modeling. Moreover, we scale the hyperparameter search on an increasing number of compute nodes to see how scaling impacts the overall uncertainty estimation in DCRNN.

The third step of our approach is generative modeling. We model the joint probability distribution of the high-performing hyperparameter configurations from the previous step and use the distribution for generating more high-performing hyperparameter configurations to train an ensemble of DCRNN models.

Copula is a statistical approach used to model the joint probability distribution of a given set of random variables by analyzing the dependencies between their marginal distributions. We use a multivariate Gaussian copula model [20] to learn the joint probability distribution between high-performing hyperparameter configurations. The copula model $C$ takes multiple marginal distributions of the hyperparameters $u_{1},u_{2},...,u_{hp}$ and returns a multivariate distribution as follows:

where $\Phi_{\sum}$ represents the cumulative distribution function (CDF) of a multivariate normal, with covariance $\sum$ and mean 0, and $\Phi^{-1}$ is the inverse CDF for the standard normal. This captures the joint probability distribution of the high-performing hyperparameter values such as learning rate, batch size, number of layers, and objective function.

The rationale for generative modeling through multivariate Gaussian copula is twofold. First, in hyperparameter search problems characterized by integer, real, and categorical hyperparameters, there are multiple high-performing regions; while BO cannot sample these regions uniformly, it can sample a few configurations from different high-performing regions. By selecting these high-performing configurations and fitting a Gaussian copula, we can sample from the sparsely sampled high-performing regions as well. Second, if we have to find more high-performing configurations from different regions directly with BO, then we have to run the search multiple times with different initialization; this process will be both computationally expensive and resource intensive.

Once the copula model $C$ is trained, we sample a number of hyperparameter configurations from $C$ and run multiple DCRNN training sessions with the sampled hyperparameter configurations. From these trained DCRNN models, we select a set of high-performing models for ensemble construction.

The fourth step of our approach is to use the ensemble of DCRNN models to estimate the traffic metrics and the associated aleatoric and epistemic uncertainty using variance decomposition.

Let $\mathcal{E}={\theta_{i},i=1,2,...,M}$ be the ensemble of $M$ models. Previous studies have showed that if each model is configured to predict the mean and variance of an output distribution that is assumed to be Gaussian and if the model is trained with negative log likelihood loss [21], the empirical estimation of the mean and the total variance can be written as follows.

The total uncertainty quantified by $\sigma_{\mathcal{E}}^{2}$ is a combination of aleatoric and epistemic uncertainty, which is given by the mean of the predictive variance of each model in the ensemble and the predictive variance of the mean of each model in the ensemble.

Our proposed SQR-based DCRNN seeks to predict quantiles of the distribution as opposed to the mean and variance of the distribution. Thus it is more generic and free from Gaussian assumptions. We derive the empirical estimate of the median ($q_{50\%}$) and the total variance of the ensemble as

where $q_{33\%}$ and $q_{66\%}$ are the 33% and 66% quantiles of the output distribution, respectively.

Note that Eq. 10 becomes Eq. 9 when the output distribution is Gaussian, where the mean and median will be the same, (${q_{50\%}}_{\theta}=\mu_{\theta}$), and $({q_{66\%}}_{\theta}-{q_{33\%}}_{\theta})/2$ will be equal to $\sigma_{\theta}$.

Here we show that scaling the HPS results in an increased number of high-performing hyperparameter configurations that can improve the quality of Gaussian copula model and the overall ensemble accuracy.

For the scaling experiments, we ran DeepHyper HPS runs with 10, 25, 50, and 100 workers, each with 12 hours of wall time as the budget. We note that each compute node in the cluster can run up to 2 workers simultaneously. These runs resulted in 49, 108, 203, and 453 model evaluations, respectively. The default hyperparameter configuration achieved a validation SQR loss of $2.46$. Therefore, we consider the hyperparameter configurations as high performing when they have a validation loss less than $2.46$.

Figure 2 shows the cumulative number of high-performing configurations obtained by DeepHyper as a function of search time. We observe that the number of high-performing configurations increases almost linearly with the increase in the number of workers. In particular, with 100 workers, the search achieves 50 high-performing configurations within 6 hours, whereas the search with 50 workers takes 10 hours to reach that number. In 6 hours, the 100-worker search finds more high-performing configurations than the 25- and 10-worker searches find. The number of high-performing configurations obtained by 10-, 25-, 50-, and 100-worker searches is 23, 45, 98, and 172, respectively.

Next we selected a number of high-performing configurations from the DeepHyper HPS runs for generative modeling. For a given worker-count search run, we computed a 10% quantile of the validation loss values of all hyperparameter configurations. We then selected the configurations that were less than or equal to that cut-off value. These configurations were used to fit a Gaussian copula (GC) model. We applied this procedure for each worker-count search run and denoted the resulting five GC models as GC-10, GC-25, GC-50, and GC-100. Table 1 lists the time taken by each GC model to fit the configurations and sample new hyperparameter configurations from the distribution. We observe that time required for fitting the GC model and sampling from the model is negligible (less than 3 seconds for fitting and 0.04 second for sampling on a single CPU node). Once they were fitted, we sampled 100 new hyperparameter configurations from each GC model and trained all them simultaneously on 50 K80 nodes (two model training per node).

Figure 3 shows a box plot of the validation loss distribution of 100 trained models from GC-10, GC-25, GC-50, and GC-100. The results clearly show that the validation loss significantly decreases as we move to GC models trained with a larger number of workers. The median values of the validation loss distributions are 3.24, 3.0, 2.77, and 2.47 for GC models trained with 10-, 25-, 50-, and 100-workers HPS results, respectively.

For model selection, we selected the top 25 models from the generated 100 based on the validation loss. The number 25 is motivated by the fact that the state-of-the-art SG-MCMC method used 25 posterior sample models from the trained model.

We also compared the validation loss distributions of (1) all the hyperparameter configurations in the HPS run, (2) 10% quantile subset of the HPS run used for training the GC-100 model, (3) 100 models from GC-100, and (4) the top 25 of 100 models from GC-100. The results are shown in Figure 4. The results show that the validation loss distribution achieved by the top 25 models is superior to the hyperparameter values found by the search alone. Moreover, using all 100 models from the GC model results in a validation loss distribution that is better than that of HPS run. The strategy of generating models from GC and selecting a subset from them seems a good strategy. We expect that by generating a large number of high-performing configurations by scaling HPS further and/or by using faster GPUs, we can train better GC models, which will result in overall improvement in the ensemble validation loss distribution.

Using the selected 25 models, we performed inference on the test data and computed the test loss. The results are shown in Figure 5. We observe that the 25 models generated with a larger number of workers for the HPS significantly reduce the test loss. The trend is similar to the validation loss distribution with 100 models. The medians of the distribution are 3.96, 3.89, 3.45, and 2.8 for 10, 25, 50, and 100 workers, respectively. The findings lead us to conclude that when conducted on a large number of workers, HPS yields better configurations, which result in improved overall predictive accuracy of the ensembles.

In this section we compare the uncertainty estimates from our proposed DCRNN-SQRUQ with several other DCRNN UQ techniques that were used in a recent study [25] and show that our method outperforms all other methods. For the experimental comparison, we adopt the same setup in terms of error metrics, forecasting horizons, and confidence interval as described in [25].

The uncertainty estimation techniques reported in [25] are as follows. (1) Bootstrap [28]: This method randomly samples 50% of training data points to train DCRNN model. Prediction is performed on the entire test dataset, and 25 samples are collected in order to calculate mean predictions and confidence intervals. (2) Quantile [29]: This method uses pinball loss for DCRNN, and the prediction is performed for three different quantiles 2.5%, 50%, and 97.5% (95% confidence interval). (3) Spline quantile regression (SQ) [30]: This method models the quantile functions using linear splines and optimizes the continuous ranked probability score as a loss function. (4) Mean Interval Score [25]: MIS is directly used in the loss function along with MAE, and a multiheaded model is used to jointly output the upper bound, lower bound, and prediction for a given input. (5) Monte Carlo Dropout [31]: This method uses the algorithm develops by Zhu et al. [32]. A 5% random drop rate is used during the testing time over 50 iterations to get a stable prediction. (6) Stochastic Gradient Markov Chain Monte Carlo (SG-MCMC) [33]: This is the current state of the art Bayesian neural network for DCRNN and uses a subsampling technique to minimize the cost of Markov chain Monte Carlo (MCMC) per iteration. A Gaussian prior $\mathcal{N}$(0, 0.4) with randomly initialization of $\mathcal{N}$(0, 0.2) is used for the model parameters. The mean and confidence intervals are obtained by averaging from 25 posterior model samples.

The comparison results are shown in the Table 2. The results of the six methods (Bootstrap, Quantile, SQ, MIS, MC Dropout, SG-MCMC) are taken directly from [25]. The performance of each method is evaluated by using the error metrics MAE and MIS (95% confidence interval). The comparison is done on three different forecasting horizons: 15 minutes, 30 minutes, and 60 minutes.

The results show that our DCRNN-SQRUQ method outperforms the state-of-the-art SG-MCMC method and a number of other UQ methods. The mean and the median MAE of our DCRNN-SQRUQ shows a better prediction accuracy in MAE than does SG-MCMC across all three forecasting horizon. Moreover, the MIS score and interval (defined in Eqs. 12 and 13) show that our method achieves a low MIS score and interval with a narrow interval compared with other techniques. An exception is for the MIS method that directly minimizes the MIS score in the loss function; as expected, it achieves the best MIS performance. However, the prediction accuracy and interval score of DCRNN-SQRUQ are better than those found by MIS. Overall, our proposed DCRNN-SQRUQ is more robust in prediction accuracy and achieves MIS and interval scores that are better than the state of the art SG-MCMC method.

Here we estimate aleatoric and epistemic uncertainty using DCRNN-SQRUQ and investigate the relationships between the traffic dynamics and the estimated uncertainty measures. According to this analysis, increased traffic dynamics produced by fast changes in traffic behavior result in high aleatoric and epistemic uncertainty. Currently, none of the existing UQ methods developed for DCRNN provide separate aleatoric and epistemic uncertainties.

Traffic dynamics describe the change in traffic behavior observed by the variability of the traffic speed. High dynamics indicate instability in the speed caused by rapid change of traffic behavior, whereas low dynamics indicate stable behavior of the traffic. The traffic dynamics at a given sensor location can be measured by the coefficient of variation (COV) [23] of the time series data (speed). The COV becomes large (small) when the variations/change in the speed values are high (low). To measure COV, we calculated the mean ($\mu$) and standard deviation ($\sigma$) of the speed of each sensor location over the entire timeline. The COV is estimated by $\sigma$/$\mu$. Afterwards, the COV on METR-LA dataset are binned in four ranges: (1) $0.2<$ COV $<0.3$; (2) $0.3<$ COV $<0.4$; (3) $0.4<$ COV $<0.5$; and (4) $0.5<$ COV. The binning helps us identify the sensor location with low to high traffic dynamics. Of a total of 207 sensor locations, 34 locations have COV between 0.2 and 0.3, 120 locations have COV between 0.3 and 0.4, 42 locations have COV between 0.4 and 0.5, and 11 locations have COV more than 0.5.

We estimate the aleatoric and epistemic uncertainty using Eq. 10. The aleatoric uncertainty is the mean of the prediction variance of the 2.5% and 97.5% quantile over the top 25 models selected from the ensemble. For each sensor location we now have the COV value and the aleatoric uncertainty measure. We separate the sensor locations based on the COV values into the four bins mentioned above and plot the distribution of aleatoric uncertainty across three forecasting horizons in box-and-whisker plot in Figure 6. We have three forecasting horizons and four COV bins, for a total of 12 box plots per figure. Across all the COV bins we can observe that 25%, median, and 75% quantiles of the aleatoric uncertainty increase with the increasing COV values. The box plots for the 15-, 30-, and 60-minute forecasting horizons show a similar trend. This result indicates that the aleatoric uncertainty is high on the sensor locations with high traffic dynamics or COV values. Since the aleatoric uncertainty is irreducible, we cannot mitigate it. In a proactive traffic management system, however, aleatoric uncertainty can be used to estimate worst-case traffic management and planning.

Next, we estimate the epistemic uncertainty by using Eq. 10. Epistemic uncertainty measures the variance of the mean across top 25 models selected from the ensemble. Similar to the aleatoric uncertainty, we separate the sensor locations based on COV values into four bins and plot the distribution of epistemic uncertainty across three forecasting horizons in box-and-whisker plots in Figure 7. A similar trend can be observed here: the epistemic or model uncertainty increases for the sensor locations having high COV value or traffic dynamics. The epistemic uncertainty for 60-minute forecasting is high for COV $>0.5$ (9th box in the Figure 7). The predictive uncertainty of the model is high due to the high traffic dynamics. Therefore, we recommend performing 15-minute forecasts for those sensor location where COV $>0.5$ (11 sensors belong to this bin). The 15-minute forecast (7th box in Figure 7) has a median of the distribution of the epistemic uncertainty below 10. Similarly, uncertainty is high for 60-minute forecasting in the COV bin between 0.4 to 0.5. Hence, we recommend performing 30-minute forecasting for those sensor locations (42 sensors belong to this bin). The ability to provide these uncertainty measures and consequent recommendations is crucial to putting these models into practice. Limiting the forecasting horizon for highly dynamic sensors will stabilize the DCRNN prediction and minimize the risk of using the forecasting results for crucial decision making.

Further, we plot 4 sets of time series data taken from 4 different COV bins in Figure 8. Each plot contains the observed and predicted timeseries data, and the epistemic and aleatoric uncertainty intervals estimated by our proposed method DCRNN-SQRUQ methods. The results shown in the Figure 8 are for the 60-minute forecasting horizon. The epistemic uncertainty interval is calculated by using the 2.5% and 97.5% quantiles of the median prediction of the top 25 models - shown in blue. The aleatoric uncertainty intervals are the medians of the 2.5% and 97.5% quantile of the top 25 models- shown in yellow. Figure 8(a) shows the sensor location with the lowest COV value (between 0.2 and 0.3). Traffic behavior does not change rapidly in this COV bin. Hence, the interval band is narrow for both of the uncertainty measures. Figure 8(b) shows the sensor location with COV value between 0.3 and 0.4. Traffic dynamics are higher compared with the previous bin. Both the aleatoric and epistemic uncertainty measures are high around 9 am and 3 pm. The traffic dynamic increases more in Figure 8(c), and the uncertainty band becomes even wider here. Epistemic uncertainty estimates are high during the drop/rise in traffic speed. In contrast, aleatoric uncertainty is high around 6 am, 9 am, 12 pm, and 11 pm. Specifically, during 9 am the uncertainty is high but the traffic is stable. This uncertainty is coming from the data. Therefore, to validate it, we plot the timeseries data of the same sensor (sensor ID: 716960) over multiple days in box-and-whisker plot shown in the Figure 9. The results show huge variation in the speed profile during 3 am, 9 am, 12 pm, and 11 pm. We also plot the aleatoric uncertainty on top of the daily variation of the data. We see that our uncertainty measure covers the variation in the data. This explains the high uncertainty at 9 am in Figure 8(c). Likewise, the traffic dynamics increase more in Figure 8(d) where the COV is greater than 0.5. The aleatoric uncertainty is high from 3 am to 7 am because of the variety of patterns in the data. The epistemic uncertainty is high during the drop/rise of traffic speed, which is difficult for the model to capture.