{'fraction_non_alphanumeric': 0.02287045370321351, 'fraction_numerical': 0.01328815712848748, 'mean_word_length': 4.5937222386733785, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Sparse Density Representations for Simultaneous Inference on Large Spatial Datasets

Large spatial datasets often represent a number of spatial point processes generated by distinct entities or classes of events. When crossed with covariates, such as discrete time buckets, this can quickly result in a data set with millions of individual density estimates. Applications that require simultaneous access to a substantial subset of these estimates become resource constrained when densities are stored in complex and incompatible formats. We present a method for representing spatial densities along the nodes of sparsely populated trees. Fast algorithms are provided for performing set operations and queries on the resulting compact tree structures. The speed and simplicity of the approach is demonstrated on both real and simulated spatial data.


Introduction and motivation
Many large streaming datasets are generated by systems that continually record the location of spatiallycentric events. Examples include the address where ambulances are requested through an emergency dispatch center, location names detected in a query through a web search engine, and the GPS coordinates recorded by a taxicab's routing software. When spatial measurements come from a relatively small and discrete set, they can often be handled just like any other covariate. Spatial data recorded from a very large set, such as addresses, or as a continuous measurement of latitude and longitude, may require more specialized treatment. Predicting the location of the next observation, for example, is often accomplished by modeling historical data as a spatial point process and estimating a density function with techniques such as kernel density estimators or Gaussian mixtures.
With enough data, it may be both possible and desirable to generate alternative density estimates for different classes of events. When crossing several variables, densities can quickly become quite granular; consider, for example,the density of ambulance requests 'related to traffic accidents from 7am-8am on Tuesdays'. An increase in granularity comes with a corresponding increase in the number of density estimates that must be stored. On its own, this increase is not a particular point of concern. The estimates can be stored in an off-the-shelf database solution, a small subset relevant to a given application can be queried, and the results used in a standard fashion.
Computational issues do arise when there is a need to simultaneously work with a large set of densities. Inverse problems, where an event is detected in a location and the type of event must be predicted, present a common examples of such an application. Performing set operations, such as intersections and unions, to construct even more complex density estimates is another example. A large set of non-parametric density estimates, represented by predictions over a finegrained lattice, can lead to a substantial amount of memory consumption. Parametric estimates, perhaps compactly represented by Gaussian mixtures, may be safely loaded into memory but eventually cause problems as set operations between densities yield increasingly complex results. These make inference and visualization intractable as the scale of the data increases.

Prior work
The estimation of general density estimation has a long history, with many spatial density routines being a simple application of generic techniques to the two-dimensional case. Mixture models were studied as early as 1894 by Karl Pearson [14], with nonparametric techniques, such as kernel density estimation, being well known by at least the 1950s [13], [16]. Further techniques have also been developed in the intervening years, such as smoothed histograms [10], splines [2,5], Poisson regression [5], and hierarchical Bayes [19]. A variant on mixture models that allow for densities to depend on additional covariates was proposed by Schellhase and Kauermann [18], though their work focuses only on one-dimensional models.
Computational issues regarding density estimation have also been studied in the recent literature. Some attention has focused on calculating the kernel density estimates given by a large underlying training dataset [15,22]. As an active area of computer vision, there as been a particular interest on two-dimensional problems [4,20,23]. These often utilize some variation of the fast Gauss transform of Greengard and Strain [9]. To the best of our knowledge, no prior work has focused on the computational strains of working simultaneously with a large set of density estimates.
A significant amount of literature and software also exists from the perspective of manipulating generic spatial data within a database. These store spatial data in a format optimized for some set of spatial queries, such as k-nearest neighbors or spatial joins [21], and can be made to handle and query fairly large sets of data.
Recent developments have even made spatial queries possible over data distributed across horizontally scalable networks [12], and integrated for fast real-time visualizations [11]. The elements in such systems are typically points, lines, or polygons. Density functions can be encoded into such a system by storing the centroid of the density modes within mixture models, or saving predictions over a fine grid for non-parametric models. The problem in using this solution for our goal is that operations on the parametric models do not translate into fast queries on the database, and the size of the grid required to store non-parametric models quickly grows prohibitively large. Our solution builds a density estimate that can be queried within a database system, but stored in a significantly smaller space.

Approach
For the remainder of this article, we focus on density estimation over a rectangular grid of points. We assume that there is some observed sample density y i over the grid, which generally corresponds to assigning observations to the nearest grid point, and an initial estimate of the density given by z i . The latter typically comes from a kernel density estimate, but may be generated by any appropriate mechanism. Our goal is to calculate predictions x i such that the following are all true: 1. The new estimates are expected to be nearly equally as predictive for a new set of observed data y new as the estimates z i , The estimates x i can be represented as a sparse vector over some fixed dictionary D, which has a total dimension of size O(n), the number of data points in the grid.
3. Queries of the form i∈A x i can typically be calculated faster than O(|A|) for A corresponding to any subset of points over the original grid.
To satisfy these requirements, we construct an overcomplete, hierarchical dictionary and calculate an estimate of z using an 1 -penalized regression model.

Quadtree dictionary
In order to describe our target dictionary, it is best to consider the case where the grid of points is a square with 2 k points on each side. We otherwise embed the observed grid in the smallest such square and proceed as usual knowing that we can throw away any empty elements in the final result.
For any integer m 0 between 0 and k, we define the following sets: For a fixed m 0 , these sets produce a disjoint, equally sized, partition of the grid of points, with m 1 and m 2 giving the horizontal and vertical coordinates of the superimposed grid. We refer to each T m 0 m1,m2 as a tile, as it represents a square subset of the original lattice. For a visualization of the partitioning scheme over a small grid, see Figure 1. It will be helpful to have a way of referencing the set of all tiles that share a given m 0 parameter, often referred to as a zoom level: Combining all of the zoom levels, we can construct the set of all tiles from which we will define the final dictionary: We will assume that there is some fixed order of this set, so that we may refer unambiguously to the j-th tile τ j in the set T . The elements m 0 j , m 1 j , and m 2 j refer to the corresponding indices of that j-th tile. The set T consists of all tiles corresponding to a fixed depth quadtree over the grid of all points. They mirror the structure of tiles used in tilemap servers, such as slippy map generated by OpenStreetMap [17]. When the grid is defined over the entire globe, our tiles directly coincide with the slippy tiles.
Finally, we define the dictionary D α as a sparse matrix in R 4 k ×d , with d = |T |. It is defined such that: For α equal to 1, element D (1) i,j encodes the proportion of the tile τ j that is covered by the point i,j is instead a simple indicator for whether a the point x i is in tile τ j . Other values of α provide a continuous scale of weightings to the tiles that moves between these two extremes, and will be useful in the penalized estimation routine. Notice that the dictionary is an over-representation of the space of all possible estimators, and for α equal to 0, contains a copy of the identity matrix as a permuted subset of its columns.

Estimation algorithm
In order to learn a sparse representation of z given our dictionary D α , we first use an 1 -penalized estimator. This well-known technique produces a parsimonious estimator via convex optimization. The parameter α is set to a number between 0 and 1 to control the degree to which the penalty should be proportional to tile size. A large parameter will yield a model with many small tiles, whereas a small parameter gives a smaller model with more large tiles. We find that a value near 0.5 typically works well.
In order to solve Equation 6, we calculate the regularization path for a sequence of λ values using software that provides a customized application of coordinate decent [8]. Choosing the final tuning parameter can be done by a number of methods; we have found the one standard deviation rule, originally suggested by Breiman, with 5-fold cross validation provides good predictability without overfitting the model [3].
To increase interpretability and reduce prediction bias, we calculate the non-negative least squares estimator over the support of the 1 -penalized estimator. With a large training dataset, this should be refit on a holdout set. Empirically, we observe good performance even when refitting on the same training data. Without the 1 -penalty the choice of α does not effect the predicted values, so here we use D 1 to facilitate the interpretability and normalization of the results.
We then hard threshold the non-negative least squares solution by a value of δ.
And finally, the density estimator is normalized to have a sum of 1: The predicted values, x, can be calculated by projecting b by the dictionary D (1) .
Clearly we do not want to save the predicted values explicitly. Otherwise, we would have simply saved the raw predicted values z. We instead show, in the next subsection, that predictions x can be generated quickly from the sparse vector b.

Fast query techniques
Given the quadtree nature of dictionary D, any particular grid point x i is represented by at most k + 1 = log 4 (n) + 1 elements. This implies that we can calculate the predicted value at that point in O(log(n)) time. Due to the nature of the 1 -penalty, the number of non-zero terms of β should be small; we denote this by || β|| 0 = s. Given the hardthresholding in our example, s is at most δ −1 but will typically be much smaller. From this it is possible to calculate i∈A x i in O(s log(n)) time.
Assume now that we have two estimates density estimates, b A and b B , which represent the densities P(x i |A) and P(x i |B) for disjoint events A and B. In order to estimate P(x i |A∪B), we can take the weighted sum of their sparse representations: The final representation may be further hard thresholded by δ to guarantee that b A∪B has no more than δ −1 elements. We can analogously use the same method to take the union of an large set of densities.
Finally, consider observing two independent samples in which we first observe event A and then observe the event B. Conditioned on the fact that both were observed in the same location, we wish to calculate the density over space given the estimators P(x i |A) and P(x i |A). We call the event of observing A followed by B in the same location' C'. Applications of this type of event arise when a single person or device is observed displaying two types of events in close temporal proximity to one another (in other words, too fast to have moved given the granularity of the grid).
We cannot directly calculate this density of the nonzero elements of the sparse estimators, however clearly we have the following relationship between the sparse representations: Calculating this directly would require computing estimates of the original densities at all n grid points. However, notice that because of the hierarchical structure of the dictionary, D (1) b A has at most || b A || 0 + 1 unique values. This is due to the fact that for any two tiles either one contains the other or they are disjoint. Likewise, the product

Simulation from a Gaussian mixture model
We sampled a set of two hundred thousand points from a Gaussian mixture model with 6 modes. Using δ = 0.001, we calculated the sparse representation of the density over the quadtree based dictionary. A scatterplot of the sampled points superimposed over the predicted densities are shown in Figure 2 for a particular grid size and value of α. Table 1 displays the total variation distance between the estimated density and the true Gaussian distribution over a grid of tile sizes and values of α. These are compared to the total variation distances gained by directly using the observed densities at each point. We see that, with the exception of the smallest grid size, the total variation of both estimators performs more poorly with larger grid sizes. This seems reasonable as larger grid sizes provide a harder estimation problem with many more degrees of freedom. The total variation as a function of α shows that values too close to 1 perform poorly regardless of the grid size. Values near 0 function okay, but the optimal value tends to be from 0.25 to 0.75.
In Table 2, we show the number of non-zero elements in the prediction vector by grid size and the α parameter. The number of non-zero elements generally   increases with the grid size. Given the exponentially growing dictionary size this is not surprising, however it is interesting that the sizes do not grow exponentially. This shows the ability of our method to represent complex, granular densities in a relatively compact way regardless of the grid size. When α is equal to 1 (or very close to it) the predicted model size is very small. This is the result of a poor model fit due to overfitting on small tile sizes; the cross-validation procedure eliminates the overfitting and leaves a rel-atively constant model. We recommend setting α to anything between 0 and the value that maximized the model size. In this simulation notice that maximizing the model size with respect to α would generate the lowest total variation of the resulting model. This is a useful heuristic because while Table 1 cannot be generated without knowing the true generating distribution, Table 2 can be created without an external information.  peace  vehicle deceptive  theft narcotics battery violation  theft  practice  theft  143  226  210  195  154  171  narcotics  239  241  254  219  257  battery  156  226  195  244  peace violation  191  197  228  Union   vehicle theft  107  187  deceptive practice  186  theft  143  227  190  199  157  188  narcotics  231  266  287  246  151  battery  152  218  157  168  peace violation  174  197  166  vehicle theft  106 168 Intersection deceptive practice 173

Chicago crime data
The city of Chicago releases incident level data for all reported crimes that occur within the city [1]. We fit sparse density estimates to six classes of crimes; these estimated densities are shown in Figure 3. We have picked classes that exhibit very different spatial densities over the city. For example, deceptive practice crimes occur predominantly in the center of the city whereas narcotics violations are concentrated in the western and southern edges of Chicago. Table 3 shows the number of non-zero terms in the predicted density vectors of pairwise unions and intersections. Of particular importance, note that the complexity of the unions and intersections are not significantly larger than the complexity of the original estimates (given by diagonal terms of the table of unions). This property would not hold for most other density estimation algorithms.

Uber pickup locations in NYC
In response to a freedom of information request, the New York City government released a dataset showing the requested pickup locations from 5 million rides commissioned by the transit company Uber [7]. We used this data to construct two-hour density buckets, three of which are shown in Figure 4. Notice that, unlike the Chicago crime data, the non-zero tiles are relatively consistent from image to image, roughly following the population density of the city. Temporal differences do exist: for example the heavily neighborhood of the Upper West Side has a particularly high density only during the morning commute. The ability to detect localized spikes at the airports exhibit the adaptive nature of the sparse learning algorithm. Table 4 shows the model sizes when computing the union of densities from any continuous time interval during the day. Due to the truncating of small densities by δ and the fact that the non-zero tiles generally line up across time periods, the overall size of the unions never grows much larger than the original estimates. Table 5 shows the same information over arbitrary unions. These intersections would be useful, for example, when trying to determine where taxicab waiting spots should be constructed as they indicate areas of high density throughout periods of the day.
Overall, we are able to quickly calculate 144 complex densities by only estimating and storing 12 of them (or 4095, when considering non-contiguous time periods).

Conclusions and future extensions
We have presented an algorithm for calculating sparse representations of spatial densities. This method has been shown to be able to compute fast density estimations over arbitrary regions and to support union and intersections over a large set of independent density estimates. These claims have been illustrated theoretically, by controlled simulations, and over two real datasets. We are now looking to generalize the quadtree approach to work with alternative hierarchical partitions of space. In particular, this could create a finer grained dictionary near places of high density (i.e., roads and city centers) allowing for a smaller dictionary and further decreasing the estimation error.