{'fraction_non_alphanumeric': 0.047252387193409476, 'fraction_numerical': 0.01802096985583224, 'mean_word_length': 4.152417701676113, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 19, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Rapid Detection of Hot-spots via Tensor Decomposition with applications to Crime Rate Data

We propose an efficient statistical method (denoted as SSR-Tensor) to robustly and quickly detect hot-spots that are sparse and temporal-consistent in a spatial-temporal dataset through the tensor decomposition. Our main idea is first to build an SSR model to decompose the tensor data into a Smooth global trend mean, Sparse local hot-spots, and Residuals. Next, tensor decomposition is utilized as follows: bases are introduced to describe within-dimension correlation, and tensor products are used for between-dimension interaction. Then, a combination of LASSO and fused LASSO is used to estimate the model parameters, where an efficient recursive estimation procedure is developed based on the large-scale convex optimization, where we first transform the general LASSO optimization into regular LASSO optimization and apply FISTA to solve it with the fastest convergence rate. Finally, a CUSUM procedure is applied to detect when and where the hot-spot event occurs. We compare the performance of the proposed method in a numerical simulation study and a real-world case study, which contains a dataset including a collection of three types of crime rates for U.S. mainland states during the year 1965-2014. In both cases, the proposed SSR-Tensor is able to achieve the fast detection and accurate localization of the hot-spots.


Introduction
The objective of our research is hot-spot detection when monitoring multiple data sources or streams across different spatial regions over time, where one is interested in quickly and accurately determining which data sources or streams change their patterns at which local regions and at which time. And in our paper, the definition of hot-spot is the structured outliers that are sparse over the spatial domain but persistent over time. A concrete motivating example is to monitor three types of annual crime rates from 1965 to 2014 for 48 mainland states in the United States, see Section 2 below CONTACT Author Email: yzhao471@gatech.edu arXiv:2004.11710v3 [stat.AP] 14 May 2020 for the detailed data description. There are two kinds of changes: one is at the global level, and the other is at the local level. For hot-spot detection, we are more interested in detecting those local changes with the following two properties: (1) spatial sparsity, i.e., the local changes are sparse in the spatial domain; (2) temporal persistence, i.e., the local changes last for a reasonably long time period unless one takes some actions.
Generally speaking, hot-spot detection in spatio-temporal data can be treated as a change-point detection problem. There are three major categories of methodologies and approaches in the literature. The first one is the LASSO-based control chart that integrates LASSO estimators for change point detection and declares non-zero components of the LASSO estimators as the hot-spot, see [22,[32][33][34]. Unfortunately, the LASSO-based control chart lacks the ability to separate the local hot-spots from the global trend mean in the spatio-temporal data. The second category of methods is the dimension-reduction-based control chart, where one monitors the features from Principal Component Analysis (PCA) or other dimension reduction methods. For example, [16] reduces the dimensionality in spatio-temporal data by constructing T2 and Q charts. [19] combines multivariate functional PCA with change-point models to detect the hot-spots. For other dimension reduction methods, please see [1,3,9,17,18] for more details. The drawback of PCA or other dimension reduction based methods are the restriction of the change detection problem and the failure to take full advantage of the spatial clustering property of hot-spot. The third category of hot-spot detection methods from spatio-temporal data is the decomposition-based method that decomposes the hot-spot from background events. For example, [27] and [28] proposed Smooth-Sparse Decomposition (SSD) model for hot-spot detection in the spatio-temporal data. SSD can separate hot-spot from the functional mean by utilizing the spatial structure of both the functional mean and hot-spot. More references can be found in [13,25,26,29,30]. However, these existing approaches investigate structured images or curves data and assume that the hot-spot events are independent over the time domain.
In this paper, we propose to develop a decomposition-based hot-spot detection method when the hot-spots are from autoregressive (AR) model, which is typical for time series data. It is worth noting that the spatio-temporal data can often be represented in 3-dimensional tensor format as "Spatial dimension × Temporal dimension × Attributes". Our main idea is to decompose this tensor into three components: smooth global trend mean, sparse local hot-spot, and residuals. We term our proposed decomposition model as SSR-Tensor. Furthermore, when fitting the raw data to the SSR-Tensor, we propose to add two penalty functions: one is the LASSO type penalty to guarantee the sparsity of hot-spots, and the other is the fused-LASSO type penalty to the autoregressive properties of a hot-spot or time-series data. Through our proposed SSR-Tensor model, we are able to (1) detect when the hot-spot happens (i.e., the change point detection problem); and (2) localize where and which type of the hot-spots occurs if the change happens (i.e., the spatial localization problem). We call the first capacity as hot-spot detection and the second capacity as hot-spot localization.
Considerable research has been done on modeling and prediction of the spatiotemporal data. Some popular time series models are AR, MA, ARMA model, etc. And the parameters can be estimated by the Yule-Walker method [7], maximum likelihood estimation or the least squares method [6]. In addition, spatial statistics have also been extensively investigated in their own right, see [5,12,14,21]. Spatio-temporal models are proposed by combining the time series models with spatial statistics [11,31]. Please see the textbook [4] for additional literature and detailed discussion. We emphasize that our proposed SSR-Tensor model is different from these existing spatio-temporal models in the sense that its primary objective is for hot-spot detection, not for estimation or prediction.
While our paper focuses only on a 3-dimensional tensor due to our motivating application in crime rates, our proposed SSR-Tensor model can easily be extended to any d-dimensional tensor ( d ≥ 3 ). If we have additional dimensions, such as the unemployment rate and economic performance, these can be modeled as the additional dimension in the tensor analysis. The reason is that our proposed model uses the basis to describe correlation within each dimension, and utilizes tensor products for interaction between different dimensions. Thus, as the dimension d increases, we just need to add the corresponding bases. The capability of extending to high-dimensional data is the main advantage of our proposed SSR-Tensor model.
The remainder of this paper is as follows: Section 2 introduces and visualizes the crime rate dataset, which is used as our motivating example. Section 3 presents our proposed SSR-Tensor model and discusses how to estimate model parameters from data. Section 4 describes how to use our proposed SSR-Tensor model to find hotspots, both for detection and localization. Our proposed methods are validated through extensive simulations in Section 5 and the case study using the crime rate dataset is shown in Section 6.

Motivating Example & Background
This section gives a detailed description of the crime rate dataset used in this paper. The dataset is available from the website of the U.S. Department of Justice Federal Bureau of Investigation, see https://www.ucrdatatool.gov/Search/Crime/State/ StateCrime.cfm. The crime rate dataset is recorded from 1965 to 2014 for 48 mainland states in the U.S. annually. In each year and for each state, three types of crime crates are recorded: a) Murder and non-negligent manslaughter ; b) Legacy rape; and c) Revised rape. These three annual crime rates are denoted by r 1 , r 2 and r 3 , respectively.
It is worth noting that our motivating dataset is of three dimensions, which includes the year (temporal dimension), state (spatial dimension) and different types of crime rates (attribute/category dimension). For the purpose of clarity and visual representation, we plot several figures which shows the characteristic of each dimension. To begin with, we first show the characteristic of the temporal domain (year), where we plot the time series of the annual crime rate for the entire U.S. in Figure 1(a). The x-axis is the year ranging from 1965 to 2014, and the y-axis is the annual crime rate of U.S.. It can be seen that the crime rates are increasing in the first ten years during 1965-1975, then become stationary during 1975-1995, and finally have a decreasing trend during 1995-2014. Furthermore, we also highlight the two peaks around 1980 and 1992 as well, because the reason to cause these two peaks fails to be determined easily. Probably it is caused by the global trend, or it is also possible to be caused by the local hot-spots.
Then we show the characteristic of our motivating data on the category domain (type of the crime rates) in Figure 1(b), where different bars represents different type of the crime rates, and the height of the bar represents the cumulative crime rate from 1965 to 2014. It can be seen that, these three crime rates overall happen with similar frequencies, which makes it different to detect the hot-spot if we compress the three-dimension data into one dimension only.
Finally, we show the characteristic of our motivating data on the spatial domain (state) in Figure 2. In Figure 2, each map shows the spatial information of the third crime rate (revised rape) in six different years. And the selected six years is starting from 1965, and ranging with a ten-year interval. For the sixth map, since the data in Year 2015 is not available yet, we use Year 2014 instead. We can see from the spatial plot in Figure 2 that the crime rate spatial patterns are different in different years. In addition, the crime rates are increasing in the first ten years and are decreasing in the last ten years, which is consistent with that in Figure 1(a).
From Figure 1 and Figure 2, there seems to have a brief increasing trend during 1984-1995, but it is difficult to visually conclude whether this is due to the global change or local hot-spots without refined analysis. If this is a local hot-spot event, we want to detect when the hot-spot appears, and then identify where and which type of crime rate accounts for this change. Since the data have three dimensions (states, rates and year), how to properly model the global trend mean and local hot-spot structures become very crucial, which is the motivation for our proposed model described in Section 3.1.   It is worth noting that we store our data as a three-dimension tensor, and we noted it as Y. Mathematically speaking, the element Y i,j,t represents the j-th crime rate of state i in year t, where i = 1, . . . , 48 for 48 mainland states, j = 1, . . . , 3 for three different type of crime rate and t = 1, . . . , 50 for 50 years from 1965 to 2014. In the next sections, we will model this multi-dimensional array through tensor decomposition. For the convenience of notation and ease of understanding, we first introduce some basic tensor algebra and notation, including basic notation, definitions, and operators in tensor (multi-linear) algebra that are useful in this paper. Throughout the paper, scalars are denoted by lowercase letters (e.g., θ), vectors are denoted by lowercase boldface letters (θ), matrices are denoted by uppercase boldface letter (Θ), and tensors by curlicue letter (ϑ). For example, an order-K tensor is represented by ϑ ∈ R I1×···×IK , where I k represent the mode-n dimension of ϑ for k = 1, . . . , K.
The mode-n product of a tensor ϑ ∈ R I1×...×IN by a matrix B ∈ R Jn×In is a tensor A ∈ R I1×...In−1×Jn×In+1×...IN , denoted as A = ϑ × n B, where each entry of A is defined as the sum of products of corresponding entries in A and B: A i1,...,in−1,jn,in+1,...,iN = in ϑ i1,...,iN B jn,in . Here we use the notation B jn,in to refer the (j n , i n )-th entry in matrix B, And the notation ϑ i1,...,iN is used to refer to the entry in tensor ϑ with index (i 1 , . . . , i N ). And the notation A i1,...,in−1,jn,in+1,...,iN to refer the entry in tensor A with index (i 1 , . . . , i n−1 , j n , i n+1 , . . . , i N ).
The mode-n unfold of the tensor ϑ ∈ R I1×...×IN is denoted by ϑ (n) ∈ R In×(I1×...In−1×In+1×IN ) , where the column vector of ϑ (n) are the mode-n vector of ϑ. And the mode-n vector of ϑ are defined as the I n dimensional vector obtained from ϑ by varying the index i n while keeping all the other indices fixed. For example, ϑ :,2,3 is a model-1 vector. A very useful technique in tensor algebra is Tucker decomposition, which decomposes a tensor into a core tensor multiplied by a matrix along each mode: is an orthogonal I k × I k matrix and is a principal component mode-k. Tensor product can be represented equivalently by a Kronecker product, i.e., vec(Y) = (B (K) ⊗ · · · ⊗ B (1) )vec(θ), where vec(·) is the vectorized operator defined as vec(Y) = Y (K+1) (an I 1 × I 2 × · · · × I K -dimension vector). The definition of Kronecker product is as follow: Suppose B 1 ∈ R m×n and B 2 ∈ R p×q are matrices, the Kronecker product of these matrices, denoted by B 1 ⊗ B 2 , is an mq × nq block matrix defined by

Our Proposed SSR-Tensor Model
This section presents our proposed methodology. Since the crime rate data is of three dimensions, namely states, rates and year, it will likely have complex withindimension and between-dimension relationships. A within-dimension relationship includes within-state correlation, within-crime-type correlation, and within-year correlation. Between-dimension correlations include between-state-and-crime-type interaction, between-state-and-year interaction, as well as between-year-and-crime-type interaction. In order to handle these complex "within-dimension" and "between-dimension" interaction structures, we use the tensor decomposition method, where basis is used to address "within" correlation, and the tensor product is used for "between" interaction.
The choice of basis is also important since different basis can represent different spatial or temporal patterns. The detailed discussion on the selection of basis is developed in Section 3.3. The structure of this section is as follows: Section 3.1 presents our proposed model that is able to characterize the complex correlation structures; Section 3.2 develops the optimization algorithm to solve the estimation problem in Section 3.1; Section 3.3 discusses the choice of basis in practice.

Our Proposed Model
We will use order-three tensor as an example to develop the methodology, as it applies to the crime rate dataset which has three components and can be represented as a three-dimension tensor Y n1×n2×T with n 1 = 48 mainland states, n 2 = 3 different types of crime rates, and T = 50 years.
Note that the i th ,j th , and k th slice of the 3-D tensor along the dimension of state, crime type, and year can be achieved as Y i:: , Y :j: , Y ::k correspondingly, where i = 1 . . . n 1 , j = 1 . . . n 2 and k = 1 . . . T . For simplicity, we denote Y k = Y ::k . We further denote y k as the vectorized form of Y k , and y as the vectorized form of Y.
The key idea of our proposed model is to separate the global trend mean from the local pattern by decomposing the tensor y into three parts, namely the smooth global trend mean µ, local hot-spots h, and residuals e, i.e., y = µ + h + e. For the first two of the components (e.g., the global time trend mean and local hot-spots), we introduce a basis decomposition framework to represent the structure of the within correlation in the global background and local hot-spots, please see [28] for a similar concept used for image defect detection.
To be more concrete, we assume that global trend mean and local hot-spots can be represented as µ = B m θ m and h = B h θ h , where B m and B h are two bases that will be discussed below. The vectors θ m and θ h are the model coefficients vector of length n 1 n 2 T and needed to be estimated, and we will discuss the estimation method later.
Here the subscript of m and h are abbreviations for the mean and hot-spots. Since the first parameter θ m is to estimate the global trend mean, and we refer it as global mean parameter. And the second parameter θ h is to estimate the local hot-spots, and we call it as local hot-spots parameter.
It is useful to discuss how to choose the bases B m and B h , so as to characterize the complex "within" and "between" correlation or interaction structures. For the "within" correlation structures, we propose to use pre-specified bases, B m,s and B h,s , for within-state correlation in the global trend mean and hot-spot, where the subscript of s is an abbreviation for states. Similarly, B m,r and B h,r are the pre-specifed bases for within-correlation of the same type of crime rates, whereas B m,t and B h,t are the bases for within-year correlation over time. As for the "between" interaction, we use a tensor product to describe this interaction, i.e, This Kronecker product has been proved to have better computational efficiency in the tensor response data [see 10]. With the well-structured "within" and "between" interaction, our proposed model can be written as: where e∼N (0, σ 2 I) is the random noise. Mathematically speaking, both B m,s and B h,s are n 1 × n 1 matrix, B m,r and B h,r are n 2 × n 2 matrix and B m,t and B h,t are n T × n T matrix, respectively. Besides, our proposed model in (1) can also be rewritten into a tensor format: where ϑ m and ϑ h is the tensor format of θ m and θ h with dimensional n 1 × n 2 × n T . Accordingly, the (i, j, t) th element in ϑ m , ϑ h ( or equivalently the ((k − 1)n 1 n 2 + (i − 1)n 1 +j) th entry of θ h , θ m ) can estimate the global mean and hot-spots in the i th state and j th crime rate in k th year respectively. The tensor representation in equation (2) allows us to develop computationally efficient methods for estimation and prediction. After developing the models above in equation (1), we now discuss how to estimate the parameters θs (the global mean parameter θ m and the local hot-spots parameter θ h ) in our model from the data via the penalized likelihood function. We propose to add two penalties in our estimation. First, because hot-spots rarely occur, we assume that the local hot-spots parameter θ h is sparse and the majority of entries in the hot-spot coefficient θ h are zeros. Thus, we propose to add the penalty term R 1 (θ h ) = λ θ h 1 to encourage the sparsity property of θ h . Second, we assume the hot-spots are temporal consistent, as the unusual phenomenon of last year is likely to affect the performance of hot-spots in the current year. Thus, we add the second penalty to ensure the temporal continuity of the hot-spot, where θ h,t is a sub-vector of length n 1 n 2 starting from the ((t−1)n 1 n 2 +1) th element to the (tn 1 n 2 ) th element in θ h , which represents the hot-spot parameter for the t th year. By combining both two penalties, we propose to estimate the parameters (θ m , θ h ) via the following optimization problem: arg min θm,θh where θ m = vec(θ m,1 , · · · θ m,T ) and θ h = vec(θ h,1 , · · · θ h,T ). Combining these two penalties in equation (3), 23] controlling both the sparsity and temporal consistency of the hot-spots. We will discuss how to efficiently optimize equation (3) for tensors in Section 3.2.

Optimization Algorithm for Estimation
In this section, we develop an efficient optimization algorithm for solving the optimization problem in equation (3). For notation convenience, we slightly adjust the notation above. Because θ m , θ h in equation (3) is solved under penalty we re-denote θ m , θ h as θ m,λ1,λ2 , θ h,λ1,λ2 to emphasis the penalty parameter λ 1 and λ 2 . Accordingly, θ h,0,λ2 refers to the estimator only under the second penalty λ 2 R 2 (θ h ), i.e, θ h,0,λ2 = arg min θm,θh The main idea of our proposed estimation algorithm can be summarized as follows.
First we reduce the number of unknown vectors in our model, i.e, find a closed-form correlation between θ m,λ1,λ2 given θ h,λ1,λ2 . After reducing the number of parameters, we focus on estimating θ h,0,λ2 , where FISTA [see 2] is the main tool in this stage. We use notation θ h,0,λ2 to describe the corresponding estimator of θ h,0,λ2 . Finally, the estimated θ h,λ1,λ2 , described by θ h,λ1,λ2 , is solved by the closed-form of θ h,λ1,λ2 given θ h,0,λ2 . Following the main idea above, we first reduce the number of unknown vectors. Although there are two sets of parameters, namely θ m,λ1,λ2 and θ h,,λ1,λ2 in the model, we note that given θ h,λ1,λ2 , the parameter θ m,λ1,λ2 is involved in the standard least squared estimation and thus can be solved in the closed-form solution, see equation (5) in the proposition below.
Proposition 3.1. In the optimization problem shown in equation (3), when given θ h,λ1,λ2 , the closed-form solution of θ m,λ1,λ2 is given by: To estimate the parameter θ h,λ1,λ2 , we plug equation (5) into equation (3). Then, the optimization problem for estimating θ h,λ1,λ2 in equation (3) becomes arg min θh,λ 1 ,λ 2 y * − Xθ h,λ1,λ2 Because solving the inverse of a matrix B m B m is computational expensive, therefore, in order to simply the calculation, we developed the tensor format of y * , which can be rewritten as y To develop an efficient optimization algorithm to solve the global optimum of equation (3), we first solve θ h,0,λ2 , i.e., with the sparsity penalty parameter λ 1 = 0, which is used to solve θ h,λ1,λ2 for general λ 1 later. The estimation of θ h,0,λ2 is to optimize arg min θh,0,λ 2 y * − Xθ h,0,λ2 where matrix D is of dimension n 1 n 2 (n 3 − 1) × n 1 n 2 n 3 , whose (i, i) th , (i, i + n 1 n 2 − 1) th elements are of value −1 and 1, respectively. Obviously, the optimization problem for equation (7) is a generalized LASSO problem, which can be transformed into a regular LASSO problem [see 24]. The details of the transformation procedure are shown in Proposition 3.2.

Proposition 3.2. Assume that
(1) matrix A is of dimension n 1 n 2 × n 1 n 2 n 3 and constructed through combining n 1 n 2 × n 1 n 2 identity matrix rowly n 3 − 1 times, (2) matrix D of dimension n 1 n 2 n 3 × n 1 n 2 n 3 and is defined as matrix X 1 is the first n 1 n 2 (T − 1) rows of matrix X D −1 , matrix X 2 is the remaining part of matrix X D −1 , (4) the projection onto the column space of X 2 is noted as P = X 2 (X 2 X 2 ) −1 X 2 , (5) the first n 1 n 2 (n 3 −1) entries of β (β = Dθ h,0,λ2 is noted as β 1 , and the remaining are noted as β 2 . The generalized LASSO problem in equation (7) can be solved by θ h,0,λ2 = D −1β , where β = ( β 1 , β 2 ) and Proposition 3.2 allows us to focus on efficiently solving the LASSO-type optimization problem in (8). To do so, we propose to use the FISTA algorithm in [2] due to its fast convergence rate. Indeed, [2] showed that the convergence rate of the FISTA algorithm is of order O(1/k 2 ) where k indicates the iterations. FISTA is a very efficient and standardized algorithm in optimization research to solve the optimization problem involving 1 penalty. The proof of Proposition 3.2 can be found in Appendix B.
Finally, with the solved θ h,0,λ2 for λ 1 = 0, we can easily compute θ h,λ1,λ2 for general λ 1 > 0 due to the their closed form relations in the following proposition: Proposition 3.3. The closed form relationship between θ h,λ1,λ2 and θ h,0,λ2 is where is an element-wise product operator.

Selection of Bases in Practice
This section discusses how to choose the proper bases B m,s , B m,r , B m,t , B h,s , B h,r , B h,t . Generally speaking, some reasonable choices of the bases can be: (1) identity matrix when one has little to no prior knowledge of the data structure; (2) Gaussian kernel if the data shares a very smooth background; (3) Other kernels, including Cosine, Silverman and etc., depending on the nature or characteristics of the data.
In the crime rate data, we begin with the basis for the global trend mean. Figure  2 shows that the global pattern is very smooth, where no distinctive abrupt changes between neighbor states. To model the smooth spatial correlation of the global state pattern, we propose to apply kernel matrix, defined as B m,s with exp{−d 2 /(2c 2 )} for the (i, j) th element, where d is the distance between the i th state and j th state, and c is the bandwidth chosen by cross-validation. For the correlation among the crime rates, we set B m,r to be an identity matrix, since we do not have any prior knowledge of it. Next, for the bases for the hot-spots, we assume there is no prior knowledge of the hotspots. Thus we set B h,s and B h,r to be an identity matrix. Moreover, for the temporal basis in both global trend mean and hot-spots, the identity matrix is used, while this reflects that we do not know when hot-spots will occur. Our optimization algorithm includes an autoregressive-type regularization term in the estimation procedure to guarantee that there is the temporal continuity of hot-spot.

Detection and Localization of Hot-spots
This section focuses on the detection and localization of the hot-spot, which includes detecting the year (when) and localizing the state (where) as well as finding the correct crime type (which) of the hot-spot. In our case study, we focus on the upward shift of crime rates, since the increasing crime rates are generally more harmful to the societies and communities. Of course, one can also detect the downward shift with a slight modification of our proposed algorithms by multiplying −1 to the raw data.
For the ease of presentation, we first discuss the detection of the hot-spot, i.e., detect when a hot-spot occurs in Section 4.1. Then, in Section 4.2, we consider the localization of the hot-spot, i.e., determine which states and which crime types are involved for the detected hot-spot.

Detect When the Hot Spot Occurs?
To determine when the hot-spot occurs, we consider the following hypothesis test and set up the control chart for the hot-spot detection in equation (10).
where r t is the expected residuals after removing the mean. The essence of this test is that, we want to detect whether r t has a mean shift in the direction of h t (estimated in Section 3.2). Please note that h t is the subvector of h, starting from the ((t − 1)n 1 n 2 + 1) th element to the (n 1 n 2 t) th element. And h is derived by h = B h θ h . To test this hypotheses, the likelihood ratio test is applied to the residual r t at each time t, i.e., r t = y t − µ t , where it assumes that the residuals r t is independent after removing the mean and its distribution before and after the hot-spot remains the same. Accordingly, the test statistics monitoring upward shift is designed as where h + t only takes the positive part of h t with other entries as zero. Here we put a superscript "+" to emphasis that it aims for upward shift.
Unfortunately, different choices of the penalty parameters λ 1 , λ 2 gives different test statistics P + t . In order to select the one with the most power, we propose to calculate a series of P + t under different combination of (λ 1 , λ 2 ) from the set )}. For better illustration, we denote the test statistics under the penalty parameters (λ 1 , λ 2 ) as P + t (λ 1 , λ 2 ). The test statistics [see 33] with the most power to detect the change, noted asP + t , can be computed by where E(P + t (λ 1 , λ 2 )), V ar(P + t (λ 1 , λ 2 )) respectively are the mean and variance of P t (λ 1 , λ 2 ) under H 0 (e.g., for phase-I in-control samples).
Note that the penalty parameters (λ 1 , λ 2 ) detect maximization in equation (11) is generally different under different times t. To emphasize such dependence of time t, we denote the parameter pair that attains the maximization in equation (11) at time t as (λ * 1,t , λ * 2,t ), i.e, Thus, the series of the test statistics for the hot-spot at time t is P + With the test statistic available, we design a control chart based on the CUSUM procedure for the following reasons: 1) we are interested in detecting the change with the temporal continuity, therefore, alignment with the objective of CUSUM. 2) In the view of social stability, we want to keep the crime rates at a target value without sudden changes, which makes the CUSUM chart is a naturally better fit.
Specifically, in the CUSUM procedure, we compute the CUSUM statistics recursively by and W + t=0 = 0, where d is a constant and can be chosen according to the degree of the shift that we want to detect. Next, we set the control limit L as the four times of the standard derivation of P + t (λ * 1,t , λ * 2,t )(t = 1 · · · T ). Finally, whenever W + t > L at some time t = t * , we determine that a hot-spot occurs at time t * .

Localize Where and Which the Hot Spot Occur?
After the hot-spot t * has been detected by the CUSUM control chart in the previous subsection, the next step is to localize where and which crime rate may account for this hot-spot. To do so, we propose to utilize the vector at the declared hot-spot time t * and the corresponding parameter λ * 1,t * , λ * 2,t * in equation (12). For the numerical computation purpose, it is often easier to directly work with the tensor format of the hot-spot h λ * 1,t * ,λ * 2,t * , denoted as H λ * 1,t * ,λ * 2,t * , which is a tenor of dimension n 1 × n 2 × T . If the (i, j, t * ) th entry in H λ * 1,t * ,λ * 2,t * is non-zero, then we declare that there is a hot-spot for the j th crime rate type in the i th state in t * th year.

Simulation Study
In this section, we conduct simulation studies to evaluate our proposed methodologies by comparison with several benchmark methods in the literature. The structure of this section is as follow. We first present the data generation mechanism for our simulations in Section 5.1, then discuss the performance of hot-spot detection and localization in Section 5.2, and finally investigate the fitness of the global trend mean in Section 5.3.

Data Generation
In our simulation, at each time index t(t = 1 · · · T ), we generate a vector y t of length n 1 n 2 by where y i,t denotes the i-th entry in vector y t , and (Bθ t ) i denotes the i-th entry in vector the Bθ t . Besides, parameter δ denotes the change magnitude. Here 1(A) is the indicator function, which has the value 1 for all elements of A and the value 0 for all elements not in A, and w i,t is the i-th entry in the white noise vector whose entries are independent and follow N (0, 0.1 2 ) distribution.
For the anomaly setup, 1{t ≥ τ } indicates that the hot-spots only occur after the hot-spot τ . This ensures that the simulated hot-spot is temporal consistent. The second indicator function 1 i {i ∈ S h } shows that only those entries whose location index belongs set S h are assigned as local hot-spots. This ensures that the simulated hot-spot is sparse. Here we assume the change happens at τ = 20 and the hot-spots index set S h = {3, 4, 5, 45, 46, 47, 57, 58, 59, 77, 78, 79, 119, 120, 121, 137, 138, 139}.
To match the dimension in the case study, we choose n 1 = 48, n 2 = 3, T = 50. For the three terms on the right side of equation (13), they serve for the global trend mean, local sparse anomaly, and white noise, respectively.
In our simulation, the matrix B is a fixed B-spline basis with the degree of three and ten knots. Note that the B-spline basis is only used in the generative model in simulation to generate data, but is not used in our proposed methodologies.
Moreover, the vector θ t is generated by a normal distribution, and we consider the following two scenarios: • Scenario 1: The global trend mean is stationary, in which θ t is generated by the normal distribution with mean 1 and standard deviation 0.1. • Scenario 2: The global trend mean is decreasing over time, in which θ t is generated by normal distribution with mean 0.95 t−1 and standard deviation 0.1.
Moreover, in each of these two scenarios, we further consider two subcases, depending on the value of change magnitude δ in (13): one is δ = 0.1 (small shift) and the other is δ = 0.5 (large shift).

Hot-spot Detection Performance
In this section, we compare the performance of our proposed method (denoted as 'SSRtensor') for the detection of hot-spots with some benchmark methods. Specifically, we compare our proposed method with Hotelling T 2 control chart [20] (denoted as 'T2'), LASSO-based control chart proposed by [33] (denoted as 'ZQ LASSO'), PCA-based control chart [3] and SSD proposed by [28](denoted as 'SSD').
For the basis choices of our proposed method, to model the spatial structure of the global trend mean, we choose B m,1 as the kernel matrix to describe the smoothness of the background, whose (i, j) entry is of value exp{−d 2 /(2c 2 )} where d is the distance between the i-th state and j-th state and c is the bandwidth chosen by cross-validation. In addition, we choose identity matrices for the temporal basis and crime type basis since we do not have any prior information. Moreover, we use the identity matrix for the spatial and temporal basis of the hot-spots. For SSD in [see 28], we will use the same spatial and temporal basis in order to have a fair comparison.
For evaluation, we will compute the following four criteria: (i) precision, defined as the proportion of detected hot-spots that are true hot-spots; (ii) recall, defined as the proportion of the hot-spots that are correctly identified; (iii) F-measure, a single criterion that combines the precision and recall by calculating their harmonic mean; and (iv) the corresponding average run length (ARL 1 ), a measure on the average detection delay in the special scenario when the change occurs at time t = 1. All simulation results below are based on 1000 Monte Carlo simulation replications. Table 1 and Table 2 report the simulation results of our proposed SSR-tensor method and four other baseline methods. Note that the two baseline methods, PCA and T2, cannot localize the hot-spots, and thus we do not report the corresponding values on the precision, recall, and F-measure. Moreover, in our simulation, if a method fails to detect any hot-spots within the entire T = 50 years (recall that the true hot-spot occurs at time τ = 20), we record its detection delay as 30. Thus, for those methods with a large standard deviation of ARL 1 , it is likely caused by the failure of reporting hot-spots.
From Tables 1 and 2, it is easy to see that our proposed SSR-tensor method achieves the smallest ARL 1 and largest F-measure due to the ability to capture both temporal consistency and spatial sparsity of the hot-spots. This implies that our proposed method not only provides a more immediate alarm when the hot-spots occur but also gives a more accurate estimation of the hot-spots location when they occur.
Meanwhile, the two baseline methods, SSD and ZQ LASSO, are worse than our proposed SSR-tensor method due to their inability to capture the temporal continuity of the hot-spots, particularly in Scenario 2 (decreasing global trend mean). The baseline methods, PCA and T2, perform the worst due to their inability to detect the sparse changes. In particular, T2 fails to detect the hot-spots even under the large shift scenario (i.e., ARL 1 is 30). The reason is that both T2 and PCA are designed based on a multivariate hypothesis test on the global mean change, which cannot consider the sparsity of the change as well as the non-stationary global mean trend.
In addition, we also visualize our hot-spot detection results in Figure 3. From Figure  3, we can see that our proposed SSR-Tensor method can accurately detect the hot-spot location with the smallest false positive (i.e., red). Figure 4 illustrates the trend of the detection delay, ARL 1 , of all methods, as δ   changes from 0.1 to 0.5 with the step size of 0.1. From the plot, our proposed SSR-Tensor method (the red curve) has the smaller detection delays than other baseline methods, particularly when there is a decreasing global trend mean and the magnitude of the hot-spot is small. Also, it is interesting to note that the detection delays of all methods are decreasing as the magnitude of the hot-spot is increasing, which is consistent with our intuition that it is easier to detect larger changes. In addition, note that the PCA method has the largest detection delays since it fails to consider the spatio-temporal correlation and the sparsity of the hot-spots.

Background Fitness
In this subsection, we will illustrate that our proposed SSR-tensor has a good estimation for the global trend mean. To do this, we compare the Squared-Root of Mean Square Error (SMSE) of the fitness of a global trend mean in Table 3. Here we only compare our proposed method with SSD, since other baseline methods (ZQ LASSO, PCA, and T2) can not model the global trend mean. It is clear from Tables 3 that our proposed SSR-tensor method does better in terms of the background fitness, especially in the scenario 2 (decreasing global trend mean).

Case Study
In this section, we apply our proposed SSR-tensor method to the crime rate dataset described in Section 2. Our proposed method is compared with other benchmarks (the same benchmarks we used in Section 5) from two aspects, one is the performance in the temporal detection of hot-spots (i.e., which year it occurs) and the other is the performance in the localization of the hot-spots (i.e., which state and which type of crime rates may trigger the alarm).
First, we compare the performance of the detection delay of the hot-spot. For our proposed SSR-Tensor method, we build a CUSUM control chart utilizing the test statistic in Section 4, which is shown in Figure 5. From this plot, we can see that the hot-spots are detected at 24 − 25-th, 32 − 35-th and 44-th years, i.e, 1989-1990, 1997-2000 and 2009.
For the benchmark methods for comparison, we also apply SSD [see 28], ZQ LASSO [see 33], PCA [see 3] and T2 [see 20] to the crime rate dataset and summarize the performance of the detection of a hot-spot in Table 4. Note that all the temporal changes reported in Table 4 is the year with the first alarm. Clearly, our proposed SSR-tensor method achieves the fastest detection of the hot-spot compared to all other benchmark methods. Numerically verified in simulation, we can say that it is of quite likely that 1989-1990, 1997-2000 and 2009 are indeed hot-spots, and our proposed SSR-tensor method works well in the crime rate dataset.  Table 4. Detection of hot-spot in crime rate dataset Next, after the detection of hot-spots, we need to further localize the hot-spots in Year Test Statistics Figure 5. Detection of hot-spot by SSR-Tensor the sense that we need to find out which state and which type of crime rate may lead to the occurrence of a hot-spot. Because the baseline methods, PCA and T2, can only detect when the changes happen and ZQ-LASSO fails to detect any changes, we only show the localization of hot-spot by SSR-Tensor and SSD, where the results are visualized in Figure 6. Compared to SSD, it can be seen that our proposed SSR-tensor method detects some localized sparse hot-spots, which might be useful to identify where the sudden increase of crime happens. As an example, let us consider Kansas, which is declared as a hot-spot in all the three temporal changes year 1989, 1997 and 2009. Figure 7 shows the time series plot of the first type of crime rate in Kansas. The plot shows that, during the decreasing trend of the last 30 years, Kansas experienced some sudden increase, which may cause it to be detected as a hot-spot by SSR-Tensor. Another example is Georgia, which is declared as a hot-spot in 1995. One possible explanation is that Gerogia experience huge floods caused by hurricane and tropical storm. The floods tended to damage people's houses and cut their food supply, which possibly leads to an increase in crime rates. So we can see from Figure 7(b) that, during the global decreasing trend, Georgia experienced an sudden increase around 1995. However, without our SSR-Tensor model, it may be difficult to tell whether it is a global trend or a hot-spots.

Conclusion
Hot-spot detection in spatio-temporal data is an important problem in real life. In this paper, we propose the 'SSR-Tensor' method for the detection of hot-spots in spatiotemporal data. Unlike existing methods for hot-spot detection, which is only workable for scalar or functional data, our proposed SSR-Tensor method is able to decompose the variate pattern of the multi-dimension data into the global trend mean, local hotspots and residuals. The estimation of hot-spot is solved by optimizing the sum of residuals with two penalty terms, which controls the sparsity of the hot-spots and the  temporal consistency of the hot-spots, respectively. To efficiently solve the above-mentioned high-dimensional optimization problems, we first reduce the unknown parameters and then simplify the generalized LASSO problem into the regular LASSO problem where many well-known LASSO algorithms can be used. We chose the FISTA algorithm in our paper because it currently has the fastest convergence rate up. We compare our proposed SSR-Tensor method with other benchmarks in terms of detection accuracy, computational time and background fitness. Based on Monte Carlo simulations and the case study of the crime rate dataset, we conclude that overall our proposed SSR-Tensor method outperforms other benchmarks. While the classical statistical process control (SPC) or sequential change-point detection problems have been studied for several decades, research on the hot-spot of tensor data is rather limited, mainly due to computational complexity and sparsity of hot-spots. Clearly, there are many opportunities to improve the algorithms and methodologies. For instance, it will be interesting to study the relationship between our proposed SSR-Tensor method and the autoregressive model. Moreover, in this paper, we specify the tensor basis, and then investigate its performance on the detection and localization of the hot-spot. It will be useful to investigate the robustness effects of different tensor bases. Finally, further research direction includes hot-spot detection problems with more complicated spatiotemporal dynamics.