Seasonal Disorder in Urban Traffic Patterns: A Low Rank Analysis

This article proposes several advances to sparse nonnegative matrix factorization (SNMF) as a way to identify large-scale patterns in urban traffic data. The input to our model is traffic counts organized by time and location. Nonnegative matrix factorization additively decomposes this information, organized as a matrix, into a linear sum of temporal signatures. Penalty terms encourage this factorization to concentrate on only a few temporal signatures, with weights which are not too large. Our interest here is to quantify and compare the regularity of traffic behavior, particularly across different broad temporal windows. In addition to the rank and error, we adapt a measure introduced by Hoyer to quantify sparsity in the representation. Combining these, we construct several curves which quantify error as a function of rank (the number of possible signatures) and sparsity; as rank goes up and sparsity goes down, the approximation can be better and the error should decreases. Plots of several such curves corresponding to different time windows leads to a way to compare disorder/order at different time scalewindows. In this paper, we apply our algorithms and procedures to study a taxi traffic dataset from New York City. In this dataset, we find weekly periodicity in the signatures, which allows us an extra framework for identifying outliers as significant deviations from weekly medians. We then apply our seasonal disorder analysis to the New York City traffic data and seasonal (spring, summer, winter, fall) time windows. We do find seasonal differences in traffic order.


Motivation
Traffic management is one of the persistent challenges of the modern industrialized world. It simultaneously reflects both a critical infrastructural necessity and a problem involving a wide range of scales and interactions. Over the past 2 decades, floating-car data (for example, from GPS-equipped taxis and other vehicles) has become an important source of real-time and large-scale city traffic information Deri and Moura (2015), Donovan and Work (2015). These data streams can be used to manage traffic signals, influence equilibrium traffic states Ban et al. (2011), Zheng and Liu (2017), Herman and Prigogine (1979), Mahmassani et al. (1984), Geroliminis and Daganzo (2008), Krichene et al. (2016), Deri and Moura (2015), Zhu et al. (2016), Zhan et al. (2014), Guan et al. (2016), Alonso-Mora et al. (2017), Ferreira et al. (2013), and optimally route traffic.
Our interest here is in part to quantify and compare regularity of traffic behavior, particularly across different broad temporal windows. Traffic behavior changes as seasons change, and we would like to build a framework for comparing the regularity of these behaviors. With reliable methods which allow us to identify broad trends in traffic behavior, we can understand responses to urban planning decisions and singular events, and hopefully help cities be more responsive.

Problem Statement and Contribution
We would like to find broad citywide patterns which can be used to describe macroscopic behavior of traffic counts; that is, the number of cars on a given link (directed edge) in a road network at a certain time. Two-way roads are, of course, represented as two links. The challenge is to find low-dimensional descriptions of this potentially large dataset, and to quantify disorder as the inability to fully describe traffic behavior in this low-dimensional way.
Our goal here is to decompose traffic counts on different links into different behavioral signatures. The weighted summation of such signatures will represent a low-rank approximation of the traffic behavior. In this additive decomposition method, a given behavioral signature can, for example, represent heavy traffic volumes during the morning rush hour, medium traffic for the rest of the workday, and light traffic in the evening. However, as a signature spans the entire year, this example suggests a segment of the overall pattern which regularly repeats (every workday). Both the behavioral signatures and the weights of links on the signatures should be non-negative; each signature describes a pattern of a non-negative number of vehicles on a link (e.g., activity near a public gathering place) and there are a non-negative number of people behaving according to that pattern; the mathematics of non-negative matrix factorization will play a key role.
We are interested in the error of this approximation; a large error, roughly, corresponds to disorder in traffic behavior. Our low-rank approximation depends on several parameters. By looking at this error across several seasonal windows and constructing a curve which captures the error of this approximation as it depends on these parameters, we can compare disorder in traffic patterns.

Dataset
Our method uses traffic counts on links, broken up into time increments. In this paper, we will apply our calculations to New York City traffic data given in Donovan et al. (2016). See Section 4 for more information about the dataset (and its limitations). An added component of our analysis is that it fills in some missing values via matrix completion; see [ZWFM], Xu et al. (2012).
Another somewhat novel focus of our work is using low rank factorization to quantify disorder (and in particular compare seasonal disorder); see Sect. 5. We carry out a range of low rank factorizations, finding that there are robust seasonal biases in the error.
Our work follows the ideas of Lee and Seung (2001) and Kim and Park (2007). Matrix factorization can be viewed as an approach to reduce the dimensionality or compress Asif et al. (2013) traffic data. Traffic data dimensionality techniques span more classical techniques building on principal component analysis Li et al. (2007), Yang and Qian (2019), Li et al. (2015), to more recent developments using variational autoencoders Boquet et al. (2020). Exploiting spatio-temporal patterns can both support data reduction and also prediction of future traffic states on the network, Yang and Qian (2019). For a recent review of approaches to exploit structure in traffic data for prediction, we refer to the reviews Ermagun and Levinson (2018), Nagy and Simon (2018), Pavlyuk (2019).

Outline
In Sect. 2, we set up the theory and background for the algorithm we use. In Sect. 3, we describe metrics and tools we use for analyzing error and sparsity associated with the approximation. In Sect. 4, we apply our analysis to the dataset of New York City traffic described in Sect. 1.3. We further analyze the error of our algorithm in relation to the New York City dataset in Sect. 5. We give some conclusions in Sect. 6. The appendix reviews the mathematics behind the iterative algorithm we develop in Sect. 2.3.
The supporting source code for this work is published at [KYA + ].

Traffic Counts
We start by organizing traffic counts into a matrix D, where D t, is the traffic count on link during time period t. The time index ranges through T def = {1, 2 … T} and the link index ranges through a finite set L def = {1, 2 … L} of links. Our interest is seasonal fluctuations, so we assume that T represents a large time horizon. If the traffic count on link at time t is unknown (i.e., missing), we set D t, def = NaN and define as the set of indices for which we have traffic estimates.

Non-Negative Matrix Factorization
We want to decompose D as where N is the rank of the factorization, with N ≤ min{T, L} , and where the entries of W and H are nonnegative, and where WH denotes the standard matrix product.
To quantify the error in (2), let us first define with I as in (1). In other words, we replace NaNs with zeroes. Then, define with ‖ ⋅ ‖ F being the Frobenius norm.

Sparse Nonnegative Matrix Factorization
The rank N denotes the size of the "universe" of possible signatures. We can think of the columns of W as temporal signatures and the rows of H as weights. Within this collection of temporal signatures, we want to represent each column of D (i.e., the traffic count on each link) using as few of these signatures as possible.
For any matrix A ∈ ℝ M×N , let A ⋅,n denote the n-th column of A and let A m,⋅ denote the m-th row of A. Sparse Nonnegative Matrix Factorization (SNMF), as seen in Park (2008) (cf. Hoyer (2002)), fixes positive parameters and and minimizes Sparsity in the columns of H is encouraged by the penalty term involving (which involves an L 1 LASSO-type penalty), while the term involving is used to prevent the values in W from being too large. The choice of these parameters is problem-specific.

Normalization
To provide a common reference for comparing behavior in different columns (i.e., links), we finally require that the columns of W sum to 1 (i.e., have L 1 -norm of 1). This gives us a relative traffic count; allowing us to better understand how the traffic count in the signature is temporally broken down (i.e., by time, week, and season). It also explicitly resolves ambiguity in the multiplicative decomposition D ≈ WH.
Due to this normalization, our H contains the scale-factor of D . In other words, scaling D by a constant factor leaves W unchanged but scales up H by that same factor. The choice for N and therefore are independent of the scale of entries in D.

Algorithm
Our iterative algorithm is based on Lee and Seung (2001) with sparsity modifications as laid out in Kim and Park (2008). Details of the calculation are in the appendix.
For positive integers R (number of rows) and C (number of columns), let R×C be the (R × C)-matrix whose entries are all 1. For matrices A and B in ℝ R×C , let denote Hadamard (i.e., elementwise) multiplication and division.
Starting with any full-rank (W, H) ∈ ℝ T×N + × ℝ N×L + , our iterative update is given by a sequence of recursive update rules: followed by

Algorithm Termination
There are three natural indicators that the algorithm has converged: We use a combination of the first two, that is, we stop our algorithm when for a fixed positive threshold .

Sparsity
As we increase the parameter of (5), which increases the size of the LASSO-type penalty on the columns of H, we expect H to become more sparse; i.e. to have more small entries. We can measure this via the calculations of Hoyer (2004), which are based on the relationship between the L 1 and L 2 norms. For nonzero x ∈ ℝ N , define Then, 0 ≤ Sparsity (x) ≤ 1 , with Sparsity (x) = 1 if and only if all except one of the entries of x is nonzero, and Sparsity (x) = 0 if and only if all entries of x take on a common (nonzero) value. We also note that Sparsity ( x) = Sparsity (x) for nonzero ∈ ℝ ; i.e., Sparsity is scale-invariant.
For H ∈ ℝ N×L having no nonzero columns, let us similarly define Sparsity (H) is the average of the Sparsity (H ⋅, ) over all . We note that this measure of sparsity is unaffected by the rescaling of (7).

Tradeoff
Once we have found a reasonable value for , a range of values of the sparsity penalty and the rank N (see Sect. 4.1), we can carry out a perturbative analysis of error. For a given low-rank decomposition (W, H), we can construct the error matrix We are interested in the seasonal fluctuations of the error (10) as an indication of how "ordered" traffic behavior is at different times. A small error would mean that traffic behavior can, in fact, be additively decomposed into different signatures.
We first fix a reference (see Sect. 4.1). Let's also fix a set N of possible values of the rank parameter N, and a set B of possible values of the penalty parameter . For each (N, ) ∈ N × B (i.e., a grid of values for N and ) let (W * (N, ), H * (N, )) ∈ ℝ T×N + × ℝ N×L + be the result of the iterative scheme of Sect. 2.5, with this and N, terminated according to the criterion laid out in Sect. 2.6. If in some cases, the algorithm fails to converge, then we say that (W * (N, ), H * (N, )) is undefined.
Fix T ′ ⊂ T , and define We want to compare E * for different subsets T ′ of T (e.g. different T ′ 's corresponding to different seasons). Informally, we want to compare E * (N, , T � ) and E * (N, , T �� ) for the same values of N and . A larger value of E * corresponds to more disorder on the original data. Note that our matrix factorization does not depend on the subset T ′ ; we rather are evaluating the error over different subsets of time. This allows us to use as much data as possible in the factorization, and resolve some of the challenges in dealing with sparse data. Informally, we are using a common collection of H weights, while selecting seasonal parts of the W signatures.
Plot 1 (Monotonicity of sparsity). For each N, we can plot For a fixed N ∈ N , we expect that Sparsity * should be increasing in (see Fig. 10 on page 13). If so, we can uniquely reparameterize effects of as effects of Sparsity * of (9).
Next, we can understand how Sparsity * and E * vary for each given N ∈ N .
Plot 2 (Dependence of E * on rank and ) For each N ∈ N , we can parametrically plot We expect this to be nondecreasing; more sparsity reflects more restrictions on the factors in the matrix product, leading to more error (see Figs. 11-14 on pages 13-15).
We finally can approximately plot error as a function of N for fixed sparsity values. Fix N ∈ N and let S N ∶ B → [0, 1] be the piecewise linear function with knots at (12).
Plot 3 (Dependence of E * on rank and sparsity). If S N is nondecreasing for each N ∈ N , and s ∈ (0, 1) is in the range of all of the S n 's (for n ∈ N ), we can plot (See Figs. 15,16,17, 18 on pages 15-16).
We expect that as sparsity increases, so does the error.

Analysis of New York City Data: Matrix Factorization
Let's apply our analysis to the data of Sect. 1.3 (i.e., Donovan et al. (2016)), which reverse-engineers estimates of traffic counts from origin-destination pairs for taxi trips. Our taxi dataset is an illustrative (and perhaps scaled) proxy for true traffic counts, but we recognize that it is potentially biased. Our methodology could readily be modified when (12) vs. Sparsity * (N, )) | ∈ B(N) . have empty roads at any time of the day, we reset these zero entries to NaN . To restrict our calculations to a subset for which we have somewhat reliable data, we will consider only those links for which there are at most a total of 720 missing hours (30 days worth) of data. We get L with L def = 2, 302 . Selecting more links comes with the benefit of incorporating more traffic data into our analysis while also introducing the cost of having a higher number of missing data entries. The distribution of entries is shown in Fig. 1. Multiplying these links by the time horizon of 8, 760 hours gives over 2 × 10 7 data points. Our goal is to understand how to appropriately identify and understand complex patterns in this data. Figure 1 shows that 2, 302 of the most-traveled links captures 19% of the entries in D. Adding more (sparsely traveled) links would correspond to adding sparse columns to D, perhaps with diminishing returns for data interpretation.

Constants, Initialization, and Computational Considerations
For our year-long taxi traffic dataset, on the links of L , we use: These values stem from a grid search of the results of the algorithm of Sect. 2.5. The initial entries of W and H (i.e., the initial condition of W and H) are taken to be i.i.d. Uniform(0, 1) random variables (thus ensuring that W and H start with full rank). We disregard runs which lead to zero columns of W; such zero columns are invariant under our multiplicative update rule and indicate sub-optimal use of rank. Note that, by (5) and (6), an iteration leading to a zero column of W will produce a zero row of H and vice versa. The algorithm is relatively insensitive to changes in , and the above value of leads to values of W which are not too large. Fixing the value of (to the above value), we carry out a refined grid search on the (N, ) parameters. Table 2 gives specific numerical results for two pairs of (N, ) parameters.
By thresholding H, we can forcibly express our approximation of each D as a linear combination of as few signatures as possible. For each column of H, we set all entries in each column below the 40th percentile of that column to zero. In practice, this reduces each link to a linear combination of at most eight signatures. Table 3 updates Table 2 once we have thresholded.
The results of the grid search in (N, ) are summarized in Figs. 2 and 3. These figures confirm a natural tradeoff: higher N allows lower error and higher leads to better sparsity, and these are competing objectives, however. We generally consider a sparsity value between 0.8 and 1.0 to be sufficient as it allows for only a handful of dominant signatures per link. Informally, higher values of are likely to mean sparser H, and higher values of N mean a larger universe of signatures. Hence, in cases where the algorithm produces a zero signature, we take this to mean that the matrix factorization needs a smaller rank.
To focus on the interpretation of the output of the algorithm, let us now fix N = 50 (the smaller value of N) and = 5000 for the rest of this section. The output consists of two matrices W and H of sizes 8, 760 × 50 and 50 × 2, 302 , respectively. The columns of W are L 1 -normalized and the columns of H are sparse.
Recall that columns of W represent traffic signatures over time. Each signature is a time-series for the entire year and hence need not be periodic. For example, the signatures can capture traffic anomalies during holidays, hurricanes, and blackouts.
Furthermore, the entries of a column of H are coefficients for the linear decomposition of a link into distinct signatures. For example, if the 4-th column of H is (0, 7, 2, 0, … , 0) T , then the traffic in link 4 of L can be written as seven times the second signature plus two times the third signature. This decomposition allows us to identify spatial patterns in traffic across the city. These matrices and the patterns derived from them can then aid in making specific observations about the large-scale behavior of traffic (as detailed in Sect. 4.4).

Independence of Signatures
To quantify the linear independence of the signatures obtained from the algorithm, we can compute the condition number of WHagen et al. (2000). A high condition number (in the thousands) would indicate that the rank should be reduced. By performing several runs of our algorithm for rank N = 50 with different (W, H)-initializations, we determine This number is low enough to give us confidence in the W returned by the algorithm., and further validates our choice of N = 50.

Robustness of Algorithm
The factors produced by our algorithm are not unique and can differ by permutations. If W 1 and W 2 are produced by two runs of the algorithm (starting with different initial conditions), we can calculate the Pearson product-moment correlation coefficient between the columns of W 1 and the columns of W 2 . A coefficient value close to one implies that the signatures follow the same pattern up to a scale factor. We can construct a greedy algorithm for searching for a permutation which will maximize correlations by sequencing through the columns of W 1 , finding the column of W 2 which maximizes the correlation (with respect to the selected column of W 1 ), and then removing that column of W 2 from future computations. We can thus permute the columns of W 2 to better match those of W 1 . After doing this, we can construct a heatmap which shows the correlation between the columns of W 1 and the columns of (the permuted) W 2 . Figure 4 shows this heatmap for two runs of the Algorithm of Sect. 2.5 with different (random) initial conditions and with the prevailing values of = 4171 , N = 50 and = 5000 (see Sect. 4.1). 4 gives the correlation between the i-th column of W 1 and the j-th column of the permuted W 2 . We see high correlation on the diagonal, meaning that the original matrix W 2 is very close to a permutation of W 1 .
We make observations about the low-rank decomposition are in Sect. 4.4.

Periodicity and Anomalous Observations
We note that the columns of D (and hence the signatures) are roughly periodic, with a period of 7 days. A power spectrum periodogram of D (averaged over all columns) is shown in Fig. 5.
In light of this periodicity, we can look at one day of the week across the entire year (e.g., all Mondays) and compute the hour-wise median traffic for that day and then identify anomalous behavior. See Figs. 6, 7a, 8a, and 8b (in dark red). In gray, we plot the relative taxi counts i.e., entries of the normalized signatures. We then determine which dates have relative taxi counts that differ significantly from the hourwise median traffic in the sense that those (signature, day) pairs have the highest sum of absolute deviations from the median weekly traffic. In the subsections below, we identify some possible origins of anomalous behavior. Figure 6 shows Signature 0, which captures a near-shutdown of taxi traffic on August 27, 2011. This may have Fig. 4 Heatmap of correlation coefficients of W-columns from two runs of the algorithm been caused by Hurricane Irene hitting NYC. There was an early warning and all subways and buses were shut down at noon on Saturday, August 27. A zoned taxi system was implemented at 9 am and taxis were thereafter running flat fares instead of meters [wny]. All other signatures also show similar behavior on and around August 28. Figure 7a shows the behavior of Signature 21 on February 26, 2011. The traffic deviates from the median Saturday trend. This may have been caused by a Labor Rally that took place near the New York City Town Hall [wis]. Figure  7b shows that Signature 21 is used by links near the Town Hall, which can be seen as further evidence connecting the rally to this traffic deviation.

Christmas Day
Anomalous behavior was also observed on Christmas Day. This can be seen in Fig. 8a for Signature 0 and Fig. 8b for Signature 4. Note how the Christmas Day traffic is reduced by about half throughout the day.

Endemic Signatures
After thresholding as discussed in Sect. 4.1, we can show the support { ∈ L ∶ H j, > 0} of signature n on a map, for any given n ∈ {1, 2 … N} . See Fig. 9. We note that of the 50 signatures, some tend to be geographically restricted (called endemic), while others are spread out over larger areas (called dispersed). For example, Signature 0, as seen in Fig. 9a is dispersed. Endemic signatures might sometimes explain traffic densities only on a single but long stretch of road. For example, Fig. 9b shows that Signature 10 is largely used by the the northbound 3rd Avenue and streets like Bowery, Lafayette St. and the southernmost part of Broadway that feed into 3rd Avenue. Similarly, Fig. 9d shows Signature 40 being used exclusively by a small section of the south-bound Broadway traffic near Central Park.
In some other cases, signatures can be seen as having a lateral sphere of influence in that they affect not only one street but also others feeding into or out of the street transversally-Signature 24, for example, as seen in Fig. 9c.

Analysis of New York City Data: Error Analysis
We second carry out the error analysis of Sect. 3. We will use Fig. 5 A power-spectrum periodogram showing 7-day periodicity in taxi count data averaged over all links. It is standard practice to measure y-axis values in periodograms in terms of Volt-squared. The peak at 7 shows weekly periodicity. The peaks at 3.5 days (i.e. half-timestep), 2.33 days (i.e. a third-time-step), and 1.75 days (i.e. a quartertime-step) can be seen as overtones of the 7-day periodicity Fig. 6 Signature 0 during Hurricane Irene (August 27). The other low traffic volume day marked by the dashed line is July 09 which saw a Yankees vs. Tampa Bay Rays baseball game in Bronx and was attended by more than than 48 thousand people as the sets of (11). Figure 10 gives us Plot 1. We see that Sparsity increases with an increase in for every value of N ∈ N . This confirms that we can reparametrize as Sparsity * from Sect. 3.2 for fixed values of N. Figures 11, 12, 13 14 gives Plot 2. Namely, for each season T ′ , it gives ( Sparsity * (N, ), E * (N, , T � ) . We compare each season with the yearly average and note that in spring and winter, the seasonal (relative) error of our approximation is lesser then the yearly (relative) error. This indicates that spring and winter traffic are comparatively low-rank and thus more orderly. The fall and summer traffic are more disordered and have higher approximation error than the yearly average. We also note that the gap between the season and annual approximation errors is more pronounced in fall and winter compared to spring and summer. This indicates that spring and summer are more representative of the annual traffic trends while fall and winter capture more of the anomalous behavior in our traffic dataset. Figures 15, 16, 17, 18 gives Plot 3. As with Figures 11,  12, 1314, we see a bit more order in the spring and winter, and a bit more disorder in the summer and fall. (Since we restrict the analysis to links missing at most 30 days (1/4 of a season), these conclusions we are likely to be statistically meaningful).

Variables Other Than Taxi-Counts
We looked at variables other than taxi counts and observed similar low-rank (periodic) behavior in the data. Figure 19 shows 7-day periodicity in taxi travel times, speeds, pace (reciprocal of speed) and taxi density (which was computed as a ratio of the taxi counts to link lengths).
A matrix factorization similar to the one shown here for taxi counts might be mathematically carried out for each of

Conclusions
In this paper, we studied a New York City taxi traffic dataset using SNMF techniques. This gave us some insight into underlying behavior with a small number of signatures. It also enabled identification of anomalous traffic patterns which we visualized in Sect. 4.4. These visualizations captured several widespread anomalous events.
We also developed an analysis which might suggest comparisons between traffic patterns in different circumstances (viz. different time intervals). Various city attributes can be connected to different traffic conditions and signatures. This might be relevant for adaptive urban planning and scheduling.
One result of any coarse-graining analysis, such as ours, is to help focus exploration of large datasets around patterns and events of interest. If one already understands the data and is designing, e.g., an incident detection algorithm, then other (perhaps model-based) approaches might give more refined and/or complementary insights.
This work did not directly address generalizeability of this approach to other cities. Comparable datasets from other cities would enable us to search for city-specific idiosyncratic behavior as opposed to behavior which might be somewhat invariant to the urban environment.
Starting with our analysis of seasonal error in Sect. 5, another interesting future direction would be a closer look at how different signatures (which last over the interval (13) of a year) might preferentially represent behavior of different seasons. Informally, the rank of signatures representing a particular season's behavior might be another indication of the complexity of traffic in that season. We seek to minimize this by alternating between minimization problems in W and H. Namely, if we start with a fixed (W, H) ∈ ℝ T×N + × ℝ N×L + , we can construct a descent step for the function E , (W, ⋅) and then, letting H ′ be the result, we can construct a descent step for E , (⋅, H � ) . This should decrease the value of E , , and we can then proceed iteratively.
The gradients of E , in the directions of W and H are given by Fig. 19 A power-spectrum periodogram of taxi-traffic data for travel times, pace (i.e. inverse of speed), speed and density (i.e. taxi count divided by link lengths). The data were averaged over all links. These periodograms are very similar to the one shown in Fig. 5 for taxi counts and As in Kim and Park (2008), we want to iteratively find the critical points of E , , i.e. the solutions of The above formulae suggest a multiplicative descent rule (which need not be gradient descent; see Lee and Seung (2001)