Detection of Anomalous Crowd Behavior Using Spatio-Temporal Multiresolution Model and Kronecker Sum Decompositions

by   Kristjan Greenewald, et al.

In this work we consider the problem of detecting anomalous spatio-temporal behavior in videos. Our approach is to learn the normative multiframe pixel joint distribution and detect deviations from it using a likelihood based approach. Due to the extreme lack of available training samples relative to the dimension of the distribution, we use a mean and covariance approach and consider methods of learning the spatio-temporal covariance in the low-sample regime. Our approach is to estimate the covariance using parameter reduction and sparse models. The first method considered is the representation of the covariance as a sum of Kronecker products as in (Greenewald et al 2013), which is found to be an accurate approximation in this setting. We propose learning algorithms relevant to our problem. We then consider the sparse multiresolution model of (Choi et al 2010) and apply the Kronecker product methods to it for further parameter reduction, as well as introducing modifications for enhanced efficiency and greater applicability to spatio-temporal covariance matrices. We apply our methods to the detection of crowd behavior anomalies in the University of Minnesota crowd anomaly dataset, and achieve competitive results.



page 8

page 9


Kronecker PCA Based Spatio-Temporal Modeling of Video for Dismount Classification

We consider the application of KronPCA spatio-temporal modeling techniqu...

Kronecker Sum Decompositions of Space-Time Data

In this paper we consider the use of the space vs. time Kronecker produc...

Conformal Anomaly Detection on Spatio-Temporal Observations with Missing Data

We develop a distribution-free, unsupervised anomaly detection method ca...

Identifying Coherent Anomalies in Multi-Scale Spatio-Temporal Data using Markov Random Fields

Many physical processes involve spatio-temporal observations, which can ...

Integration of spatio-temporal contrast sensitivity with a multi-slice channelized Hotelling observer

Barten's model of spatio-temporal contrast sensitivity function of human...

A spatio-temporal model to understand forest fires causality in Europe

Forest fires are the outcome of a complex interaction between environmen...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

I Introduction

The detection of changes and anomalies in imagery and video is a fundamental problem in machine vision. Traditional change detection has focused on finding the differences between static images of a scene, usually observed in a pair of images separated in time (frequently hours or days apart). As video surveillance sensors have proliferated, however, it has become possible to analyze the spatio-temporal characteristics of the scene. The analysis of the behavior of crowds (particularly of humans and/or vehicles) has become essential, with the detection of both abnormal crowd motion patterns and individual behaviors being important surveillance applications. In this paper, we focus on the problem of modeling and detecting changes and anomalies in the spatio-temporal characteristics of videos. As an application, we consider the detection of anomalous crowd behaviors in several video datasets. Our methods provide pixel-based spatio-temporal models that enable the detection of anomalies on any scale from individuals to the crowd as a whole.

I-a Approach

The two main approaches to crowd video anomaly detection are those based on first extracting the tracks and finding anomalous track configurations (microscopic methods) and those that are based directly on the video without extracting the tracks (macroscopic)

[4]. An example of the microscopic approach is the commonly used social force model of [5]. The macroscopic methods tend to be the most attractive in dense crowds because track extraction can be computationally intensive in a crowd setting, and in dense crowds track association becomes extremely difficult and error-prone [4]. In addition, it is frequently the case that long-term individual tracks are irrelevant because the characteristics of the crowd itself are of greater interest. We thus follow the macroscopic approach.

In this work, we focus on learning generative joint models of the pixels of the video themselves as opposed to the common approach of feature extraction. Our reasons for this include improving the expressivity of the model, improving the general applicability of our methods, minimizing the assumptions that must be made about the data, and eliminating preprocessing and the associated loss of information.

We follow a statistical approach to anomaly detection [6], that is, we learn the joint distribution of the nonanomalous data and declare data that is not well explained by it in some sense to be anomalies. A block diagram of the process is shown in Figure 1.

We thus need to learn the joint distribution of the video pixels. In order to learn the temporal as well as the spatial characteristics of the data, the joint distribution of the pixels across multiple adjacent video frames must be found. For learning, we make the usual assumption that the distribution is stationary over the learning interval. To make the assumption valid, it may be necessary to limit the length of the learning interval and hence the number of samples. Our approach is to learn the distribution for a finite frame chunk size , that is, the frame spatio-temporal pixel joint distribution . Once this distribution is learned, it can be efficiently extended to larger frame chunk sizes according to either an AR (Markov), MA, or ARMA process model [7]. Limiting

reduces the number of learned parameters and thus the order of the process, hence reducing the learning variance.

In order to reduce the number of samples required for learning, we use the parametric approach of learning only the mean and covariance of . The number of samples required by standard covariance methods to achieve low estimation variance grows as . Hence, the spatio-temporal “patch size” of the learned distribution that these methods can handle is still severely limited. Making the “patch size” as large as possible is highly desirable as it allows for the modeling of interactions (such as between different regions of a crowd) across much larger spatial and temporal intervals, thus allowing for larger-scale relational type anomalies to be detected.

The number of samples required for covariance learning can be vastly reduced if prior assumptions are made, usually in the form of structure imposition and/or sparsity. In this paper, we use the sums of Kronecker products covariance representation, as well as a sparse tree-based multiresolution model which incorporates both structural and learned sparsity and has a degree of sparsity that is highly tunable. The multiresolution approach is particularly attractive because it defines hidden variables that allow analysis of the video at many different scales.

Once we have the estimate of the spatio-temporal pixel distribution, the typical statistical approach [6] is to determine whether or not a video clip is anomalous by evaluating its likelihood under the learned distribution. This approach is based on the fact that the anomalous distribution is unknown, so likelihood ratio tests are inappropriate. As a measure of how well the model explains the data, the likelihood is an attractive feature. Thresholding the likelihood is supported by theoretical considerations [6].

In this work, we use the Mahanolobis distance


which is the negative loglikelihood under the multivariate Gaussian assumption. Although the data is not exactly multivariate Gaussian, the metric is convenient and robust and its use in non-Gaussian problems is well supported. In addition, it allows for evaluation of the likelihood of local subregions by extracting submatrices from (conditional distribution) or (marginal distribution), and hence allowing for relatively efficient localization of the anomalous activity. We discuss the details of likelihood thresholding in Section IV.

Figure 1: Anomaly detection block diagram

I-B Previous Work on Detection of Crowd Anomalies

In this section, we briefly review select previous work on the detection of crowd anomalies using non track based (macroscopic) approaches.

Many macroscopic techniques are based on computing the optical flow in the video, which attempts to estimate the direction and magnitude of the flow in the video at each time instant and point in the scene. In [8] optical flow is used to cluster the movement in the scene into groups (hopefully of people), and models inter group interactions using a force model. Anomalies are declared when the observed “force” is anomalously large or unexpected. Particle advection, which is based on optical flow, has also been used. The authors of [9] compute the optical flow of the video and use it to advect sets of particles. Social force modeling is then performed on the particles, and anomalies are declared based on a bag of words model of the social force fields. [4]

computes chaotic invariants on the particle advection trajectories. The chaotic invariants are then modeled using a Gaussian mixture model for anomaly declaration.

Various spatio-temporal features have also been used. In [10], the authors model the video by learning mixtures of dynamic textures, that is, patchwise multivariate state space models for the pixels. This is slightly similar to our approach in that it models the pixels directly using a Gaussian model. The patch size they use, however, is severely limited ([10] uses blocks) due to the sample paucity issues which are our main focus in this work. A spatio-temporal feature-based blockwise approach using K nearest neighbors for anomaly detection is given in [11], and [12] uses a cooccurence model. A gradient based approach is used in [13]. They divide the video into cuboids, compute the spatiotemporal gradients, and model them using sparse KDE to get likelihoods and declare anomalies accordingly.

I-C Outline

The outline of the remainder of the paper is as follows. In Section II we discuss the sums of Kronecker products covariance representation and its application to video, as well as introduce a new estimation algorithm. In Section III we review the sparse multiresolution model of Choi et al. In Section III-D we present our modifications of and application to spatio-temporal data. Our approach to anomaly detection using the learned model is presented in Section IV. Section V presents video anomaly detection results, and Section VI concludes the paper.

Ii Kronecker Product Representation of Multiframe Video Covariance

In this section, we consider the estimation of using sums of Kronecker products. Additional details of this method are found in our paper [1] and in [14].

Ii-a Basic Method

As the size of can be very large, even for moderately large and

the number of degrees of freedom (

) in the covariance matrix can greatly exceed the number of i.i.d. samples available to estimate the covariance matrix. One way to handle this problem is to introduce structure and/or sparsity into the covariance matrix, thus reducing the number of parameters to be estimated. In many spatio-temporal applications it is expected (and confirmed by experiment) that significant sparsity exists in the inverse pixel correlation matrix due to Markovian relations between neighboring pixels and frames. Sparsity alone, however, is not sufficient, and applying standard sparse methods such as GLasso directly to the spatio-temporal covariance matrix is computationally prohibitive [15].

A natural non-sparse alternative is to introduce structure is by modeling the covariance matrix as the Kronecker product of two smaller matrices, i.e.


thus reducing the number of parameters from to where (). The equivalent graphical model decomposition is shown in Figure 2

. When the measurements are Gaussian with covariance of this form they are said to follow a matrix-normal distribution

[15]. This model lends itself to coordinate decompositions [16, 14, 1]. For spatio-temporal data, we consider the natural decomposition of space (pixels) vs. time (frames) as done in [1]. In this setting, the matrix is the “spatial covariance” and is the “time covariance.”

Previous applications of the model of Equation (2) include MIMO wireless channel modeling as a transmit vs. receive decomposition [17], geostatistics [18], genomics [19], multi-task learning [20], collaborative filtering [21]

, face recognition

[22], mine detection [22], recommendation systems [16], wind speed prediction [23] and prediction of video features [1].

An extension to the representation (2) introduced in [14] approximates the covariance matrix using a sum of Kronecker product factors


where is the separation rank.

This allows for more accurate approximation of the covariance when it is not in Kronecker product form but most of its energy is in the first few Kronecker components. An algorithm for fitting the model (3) to a measured sample covariance matrix was introduced in [14] called Permuted Rank Least Squares (PRLS) and was shown to have strong high dimensional guarantees in MSE performance. The LS algorithm is an extension of the Kronecker product estimation method of [24] and is based on the SVD. In [14], some regularization is also introduced. In this work, we use different regularization methods. A pictorial representation of the basic LS algorithm is shown in Figure 3.

Figure 2: Gaussian graphical model representation of Kronecker product decomposition
Figure 3: Basic LS sum of Kronecker products approximation algorithm.

Ii-B Diagonally Corrected Method

In [1], we presented a method of performing sum of Kronecker products estimation with a modified LS objective function that ignores the errors on the diagonal. The diagonal elements are then chosen to fit the sample variances, with care to choose them so as to guarantee positive semidefiniteness of the the overall estimate. The main motivation for this method is that uncorrelated variable noise occurs in most real systems but damages the Kronecker structure of the spatio-temporal covariance. This also allows for gain in expressivity when doing diagonal regression [1]. The weighted LS solution is given by the alternating projections method (iterative).

Ii-C Temporal Stationarity

Since we are modeling a temporal process with a length much longer than , the spatio-temporal covariance that we learn should be stationary in time, that is, should be block Toeplitz [7]. If all the matrices are Toeplitz, then is block Toeplitz. Furthermore, if is block Toeplitz and has separation rank , a Kronecker expansion exists where every is Toeplitz. We thus only need to estimate the value of the diagonal and each superdiagonal of the ( parameters). To estimate these parameters, we use the method of [25], which uses the same LS objective function as in [24, 14, 1] with the additional Toeplitz constraint. The equivalent (rearranged) optimization problem is given by


The Toeplitz requirement is thus equivalent to


for some vector



Clearly . Let


for reasons that will become apparent.

We now formulate this problem as a LS low rank approximation problem. We use the notation to denote the th row of . Let be given by






and the low rank constraint is clearly manifested as


This is now in the desired low-rank form and thus solvable using the SVD or by one of the other weighted LS methods in this section.

The matrices are found from the (left singular vectors) by unweighting according to (7) and expanding them into the using (5) and then rearranging.

The use of this method in the context of the PRLS regularization is clear. An additional benefit is that the size of the low rank problem has been decreased by a factor of .

Using this method, it is clear that any block Toeplitz matrix can be expressed without error using a sum of Kronecker products. Due to symmetry, however, the matrix can be completely determined using Kronecker products.

Due to the structure of the result (9) being the same as the original optimization problem (i.e. low rank approximation, where the right singular vectors are still ), Toeplitz structure can be enforced in both the and by starting with the result of (9) and repeating the derivations in (9) except with the constrained.

When it is desirable to learn an infinite-length AR, MA, or ARMA model by finding the banded block Toeplitz covariance as in [7], this method is particularly attractive because in order to do the LS estimate, it is only necessary to find the frame sample covariance and do Toeplitz approximation with the weights removed (i.e. ).

Ii-D Sum of Kronecker Products for Nonrectangular Grids

While with enough terms any covariance can be represented as a sum of Kronecker products, the separation rank is significantly lower for those matrices that are similar to a single Kronecker product. When “flow” is occuring through the variables in time, or equivalently the best tree defined on the spatiotemporal data is nonrectangular across frames, variables (pixels) of the same index do not correspond across frames. This results in a flow of correlations through the variables as the time interval increases. This situation is produces a covariance matrix that has a very non-Kronecker structure. If, however, we shift the indexes of the variables in each frame so that corresponding (highly correlated, or adjacent in a tree graph) pixels have the same index, the approximate Kronecker structure usually returns or at least improves. The mapping is usually not one-to-one, however, so some variables in a frame will not have corresponding variables in the next frame.

To handle this, we set up a larger () space-time rectangular grid of variables that contains within it the nonrectangular (

) grid of pixel variables defined by the required index shifting, and has the remaining variables be dummy variables. This unified grid of variables is then indexed according to the rectangular grid. See the two leftmost images in Figure

4. The covariance matrix of the complete set of variables then has valid regions corresponding to the real variables, and dummy regions corresponding to the dummy variables. As an example, see the last pane in Figure 4, where the dark regions correspond to the dummy regions. If the dummy regions are allowed to take on any values we choose, it is clear that there is indeed much better Kronecker structure using the Kronecker dimensions and . To handle this for Kronecker approximation, we merely remove (don’t penalize) the terms of the standard LS objective (approximation error) function corresponding to the dummy regions of the covariance, thus allowing them in a sense to take on the values most conducive for low separation rank. After the rearrangement operator in the LS algorithm, the problem is a low-rank approximation problem with portions of the matrix to be approximated having zero error weight. This is the standard weighted LS problem, which can be solved iteratively as mentioned above. Finally, after the approximating covariance is found, the valid regions are extracted and reindexed to obtain the covariance of the nonrectangular grid.

Figure 4:

Kronecker approximation method for nonrectangular variable grids. The nonrectangular grid is embedded in a rectangular grid with dummy variable padding. The Kronecker product representation of the new rectangular grid is then found using weighted least squares.

Ii-E Examples of Low Separation Rank Processes

An example of multiple Kronecker structure is the traveling wave field


Since the Kronecker representation is exact for separable processes, two Kronecker components are required to perfectly capture the covariance of a spatio-temporal wave, although one will capture necessary information such as wavelength, nature of , speed, and amplitude.

Iii Sparse Multiresolution Model

While the sum of Kronecker products representation reduces the number of parameters considerably, most images are too large to be able to estimate or even form (due to memory issues) the matrices directly. In addition, since video characteristics do vary over space, the Kronecker decomposition can break down as the spatial patch size increases. Hence further parameter reduction is needed. Simple approaches include considering block diagonal covariance estimation, where the video is divided into spatial blocks for estimation, and/or by enforcing spatial stationarity as done in the temporal dimension. Additional reduction can be achieved by using sums of triple Kronecker products which forces slowly changing characteristics over sets of blocks using windowing. An issue, however, with doing these blockwise decompositions is that correlations between neighboring pixels in different blocks are ignored.

We thus consider the use of tree based multiresolution models. As a starting point, we consider the multiresolution model of Choi et al [2] and modify it for our problem. Choi et al’s sparse covariance model, which they refer to as a sparse in-scale conditional covariance multiresolution (SIM) model, starts with a Gaussian tree with the observed variables on the bottom row and adds sparse in-scale covariances (conditioned on the other levels) to each level (see Figure 5). The added in-scale covariances are introduced because Gaussian trees are not expressive enough and introduce artifacts such as blocks.

Figure 5: Multiresolution models: Tree, Tree with sparse in-scale covariance (conjugate graph), Equivalent graphical model

Iii-a Trees

In order to clarify the next section, we briefly review the basic Gaussian tree model. The model is a Gaussian graphical model with the variables and connections arranged in a tree shape, that is, every variable is either the root node or has a single parent

, and can have multiple children. The edge and node parameters are usually expressed implicitly by viewing the tree as a Markov chain beginning at the root node. That is, each child variable is given by



The extension to multivariate nodes is simple [26].

Iii-B Inference

In order to use this model (especially for evaluating likelihoods), it is necessary to be able to infer the hidden variables given observed variables. For our application, we observe the bottom level variables and infer the upper levels.

The general formulation is that we observe a linear combination of a set of the variables plus noise with covariance , and infer the variables via maximum likelihood estimation under a multivariate Gaussian model. In what follows, we refer to the information matrix associated with the tree with the diagonal elements removed as and the added in-scale conditional covariance matrices arranged blockdiagonally as . is a blockwise positive semidefinite matrix with nonzero values only between variables in the same level. As a result, the overall information matrix of the multiresolution model is since the inversion of a blockdiagonal matrix is blockdiagonal with the blocks being the inverses of the original blocks.

The MLE solution is [2] is to solve for in




The approach of [2] is to exploit the sparsity of the tree model and the in-scale covariance corrections via matrix splitting. In particular, the solution is found iteratively by alternating between solving (Figure 6)


for (between scale inference) using an appropriate iterative algorithm and computing the sparse matrix vector multiplication (in-scale inference)


The term can be computed by solving the sparse system of equations for . Hence each iteration is performed in approximately linear time relative to the nonzero elements in the sparse matrices.

In the case where we observe a portion of the variables () without noise () (as in our application), it is straightforward to modify the above equations using the standard conditional MLE approach as in [27], resulting in greatly enhanced computational complexity and numerical precision.

Figure 6: Representation of multiresolution inference model of [2]. Alternate between between scale and in scale inference.

Iii-C Learning

The learning algorithm proposed in [2] is based on learning the tree first and then correcting it with the in-scale covariances. To learn the tree, first specify (or learn) a tree structure where the bottom layer is the observed variables. Once the tree structure is specified, the parameters (edge weights and node noise variances) are learned from the spatio-temporal training samples using the standard EM algorithm [26]. We use the tree to represent the covariance only, thus the training data has the mean subtracted to make it zero mean.

The in-scale covariances are then learned to eliminate the artifacts that arise in tree models. The first step of the approach is to pick a target covariance of the observed variables (the sample covariance in [2]) and determine the target in-scale conditional information matrices that would result in the target covariance being achieved. The target in-scale conditional information matrices are found using a recursive bottom up approach


where and are determined by [2].

The target information matrix can be computed to be [2]


by setting the marginal covariance equal to the target covariance.

Regarding computational complexity, it is important to note that the method requires the inversion of the target covariances at each level, thus making the learning complexity at least . This can be a severe bottleneck for our application. We propose a method of dramatically reducing this cost in the next section.

Secondly, the target in-scale conditional covariances are sparsified. It is well known that applying GLASSO style logdet optimization (regularization) [2] to a matrix sparsifies its inverse while maintaining positive semidefiniteness. Hence, applying the method to the target information matrices results in sparse target covariances [2]. This gives the sparse conditional in-scale covariances as required for the model. Methods of determining all the in-scale covariances jointly are also presented in [2], but we do not consider them here due to their computational complexity.

Iii-D Modifications for Space-Time Data

In this section, we develop appropriate modifications of the multiresolution model of the previous section in order to apply it to spatio-temporal data.

Iii-D1 Structure: Spatial Tree Only

Complete stationarity in the tree model itself can be achieved by decoupling the frames (with the interframe connections filled in later using the in-scale covariance corrections). This is equivalent to having the tree model the spatial covariance only. Hence there are substantial computational and parameter reduction benefits to this approach. Another option which we do not employ here is to use space-based priors when learning the tree.

Iii-D2 Subtree Based Learning

In many of the applications we consider, it is desirable to estimate the multiframe covariance more than once (e.g. to compare different portions of frame sequences). Hence, the complexity of the inversion required to compute the target information matrix at the lower scales is prohibitive for video data. Hence, we propose to only consider local in-scale connections.

Our approach is to force the inscale conditional information matrices to be blockdiagonal, at least for the lower levels. As a result, Equation indicates that only the corresponding block diagonal elements of the right hand side are required, giving substantial computational savings immediately. Additional savings are achieved by using a local estimate for the blocks of , i.e. estimating a block of containing but somewhat larger than the block of interest, inverting, and extracting the relevant portion. This is based on the notion that the interactions relevant to local conditional dependencies should also be local, and makes the algorithm much more scalable. This is a particularly good approximation when the dominant interactions are local.

To get the upper level target covariances, it is not necessary to form the bottom level sample covariance, as, following the recursion of (18),


where is the number of samples, is the matrix of samples, is the current level, and . If necessary, regularization etc. is applied to the first term of (20). This allows for interblock connections at as low a level as possible to minimize blockwise artifacts.

Iii-D3 Kronecker

In the original multiresolution model, parameter reduction for the in-scale covariances is achieved using sparsity. Naturally, we wish to use the Kronecker PCA representation for the covariance to reduce the number of parameters. We use DC-KronPCA on the first term of (20). This is possible because the tree is in the spatial dimension only, hence the multiplication with matrices to move through the levels does not affect the temporal basis. This allows for the direct use of the Kronecker product representation without needing to invert the target information matrix first.

Iii-D4 Regularization

It should be noted that thus far we have imposed no notion of spatial stationarity or slowly varying characteristics in the model. In behavior learning, it is frequently desirable due to the paucity of samples to incorporate information from adjacent areas when learning the covariance. To achieve this type of gain using the multiresolution model, we obtain additional samples by using slightly shifted copies of the original samples.

Iv Anomaly Detection Approach

Iv-a Model: AR process

Given the mean and covariance, the standard Mahanalobis distance (Gaussian loglikelihood) is given by Equation (1).

Video is a process, not a single instance. Hence, it is frequently desirable to evaluate the Mahanalobis distance for clips longer than the learned covariance. In order to do this, the larger covariance matrix needs to be inferred from the learned one. A common approach is to assume the process is a multivariate AR process [7]. Time stationary (block Toeplitz where each block corresponds to the correlations between complete frames) covariances define a length multivariate (in space) AR (in time) process [7]. Using the AR process model, the inverse covariance is achieved by block Toeplitz extension [7] of the learned with zero padded blocks on the super and sub block diagonals. The result is then substituted into Equation (1). A memory efficient implementation is achieved by using


where and .

Iv-B Anomaly Detection

Once the likelihood of a video clip has been determined, the result is used to decide whether or not the clip is anomalous. It is common in anomaly detection to merely threshold the loglikelihood and declare low likelihoods to correspond to anomalies. In the high dimensional regime, however, the distribution of the loglikelihood of an instance given that the instance is generated by the model under which the likelihood is evaluated becomes strongly concentrated about its (nonzero) mean due to concentration of measure. For example, the loglikelihood of a

dimensional Gaussian distribution follows a chi square distribution with

parameters. As a result, high likelihoods are highly unlikely, and are thus probably anomalous. These types of anomalies frequently occur due to excessive reversion to the mean, for example, when everyone in a video leaves the scene. This situation is clearly anomalous (change has occurred) but has a likelihood close to the maximum. Hence, we threshold the likelihood both above and below.

Combination of regions with abnormally high and abnormally low likelihoods can cancel each other out in some cases, resulting in a declaration as normal using the overall likelihood alone. To address this problem, if an instance is determined to be nonanomalous using the overall likelihood, we propose to divide the video into equal sized spatial patches, extract the marginal distributions of each, and compute the loglikelihoods. If the sample variance of these loglikelihoods is abnormally large, then the instance is declared anomalous.

V Results

V-a Detection of Crowd Escape Patterns

To evaluate our methods, we apply them to the University of Minnesota crowd anomaly dataset [3]. This widely used dataset consists of surveillance style videos of normal crowds moving around, suddenly followed by universal escape behavior. There are 11 videos in the dataset, collected in three different environments. Example normal and abnormal frames for each environment are shown in Figure 7. Our goal is to learn the model of normal behavior, and then identify the anomalous behavior (escape). Since our focus is on anomaly detection rather than identifying the exact time of the onset of anomalous behavior we consider short video clips independently.

Our experimental approach is divide each video into short clips (20-30 frames) and for each clip, use the rest of the video (with the exclusion of a buffer region surrounding the test clip estimate) to estimate the normal space-time pixel covariance. Since the learning the model of normality is unsupervised, the training set always includes both normal and abnormal examples. In essence, then, we are taking each video by itself and trying to find which parts of it are least like the “average.” For simplicity, we convert the imagery to grayscale. The original videos have a label that appears when the data becomes anomalous. We remove this by cutting a bar off the top of the video (see Figure 13). The likelihood of the test clip is then evaluated using the Mahanalobis distance based on the learned spatio-temporal covariance extended into an AR process as in (21).

Since anomalous regions are included in the training data, the learning of normal behavior is dependent on the preponderance of normal training data, which is the case to some degree. Anomaly detection ROC curves are obtained by optimizing the above and below thresholds following a Neyman-Pearson approach.

In our first experiment, we use an 8 frame covariance and compare anomaly detection results for the 3 term regularized Toeplitz sum of Kronecker products and the regularized sample covariance with Toeplitz constraint (in other words an 8 term sum of Kronecker products). For mean and covariance estimation, we divide the video into 64 spatial blocks and learn the covariance for each. The test samples are obtained by extracting 30 frame sequences using a sliding window incremented by one frame. The covariance is forced to be the same over sets of 4 blocks in order to obtain more learning examples. The negative loglikelihood profile of the first video as a function of time (frame) is shown in Figure 8 using the Kronecker approach. Note the significant jump in negative loglikilihood when the anomalous behavior begins. Figure 9 shows the ROC curves for the entire dataset for the case that the thresholds are allowed to be different on each video. Figure 10 shows the results for the case that the thresholds are forced to be constant over the videos in the same environment. This reduces the performance as expected due to less overfitting. Notice that in both cases the use of the Kronecker product representation significantly improves the performance, and that the false alarm rates are quite low.

We then examined the variation of performance with the frame length of the covariance. In this experiment, 20 frame test clips were used and the covariances was held constant over all 64 blocks. Results for 1, 4, and 8 frames are shown in Figures 11 and 12 for individual video and environment thresholding respectively. Note the major gains achieved by incorporating multiframe information.

We also considered localization of the anomalies. This was accomplished by dividing the video into pixel blocks and evaluating the spatio-temporal likelihood of each. Then simple thresholding is used to determine whether or not the patch is anomalous. The results are shown in Figure 13. Note the successful detection of the individuals and only those individuals who have begun running.

Figure 7: Example video frames from each environment in the CMU dataset.
Figure 8: Example negative loglikelihood profile (first video). Anomalous behavior begins at the large jump in likelihood. The subsequent decrease is due to people leaving the scene.
Figure 9: ROC curve for anomaly detection on the entire dataset for 8 frame covariance. Thresholds are set for each video individually. The blue curve corresponds to using 3 term sum of Kroneckers (AUC .9995), and the red to the sample covariance with regularization and Toeplitz (AUC .9989). Note the superiority of the Kronecker methods.
Figure 10: ROC curves for anomaly detection on the entire dataset for 8 frame covariance. Thresholds are set for each environment (set of videos) individually. The blue curve corresponds to using 3 term sum of Kroneckers (AUC .995), and the red to the sample covariance with regularization and Toeplitz (AUC .988). Note the superiority of the Kronecker methods.
Figure 11: ROC curves for 1, 4, and 8 frame covariances. Thresholds are set for each video individually. Note the superiority of multiframe covariance to single frame covariance due to the use of temporal information.
Figure 12: ROC curves for 1, 4, and 8 frame covariances. Thresholds are set for each environment (set of videos) individually. Note the superiority of multiframe covariance to single frame covariance due to the use of temporal information.
Figure 13: Example individual detection results. Blocks declared anomalous are indicated by red boxes. The anomalous behavior is just beginning. Notice the marking of running individuals as anomalous while avoiding the walking individuals.

V-B Detection of Anomalous Patterns in Marathon Videos

As an example of crowd videos with locally steady optical flow, we consider a video of the start of a marathon, and apply our multiresolution model to learning its covariance. Since steady flow is present, we use nonrectangular tree grids. It was found this was necessary for low separation rank structure to emerge. The model was trained using the same leave out and buffer approach in the previous section. Considering only the portion of the video after the start, we are able to easily determine that clips from the original video are not anomalous whereas the same clips played backwards are anomalous (Figure 14).

Figure 14: Marathon video results using multiresolution model. Upper left: Frame before marathon starts. Upper right: Frame after steady flow has been established. Lower left: Locations of nonzero entries of multiresolution information matrix. Lower right: Negative loglikelihoods as a function of time for test clips from the video and from the video played backwards.

Vi Conclusion

We considered the use of spatio-temporal mean and covariance learning to reliable statistical behavior anomaly detection in video. A major issue with spatio-temporal pixel covariance learning is the large number of variables, which makes sample paucity a severe issue. We found that the approximate pixel covariance can be learned using relatively few training samples using several prior covariance models.

It was found that the space-time pixel covariance for crowd videos can be effectively represented as a sum of Kronecker products using only a few factors, when adjustment is made for steady flow if present. This reduces the number of samples required for learning significantly.

We also used a modified multiresolution model based on [2] and incorporating Kronecker decompositions and regularization to decrease the number of required samples to a level that made it possible to estimate the spatio-temporal covariance of the entire image. The learning algorithm in [2] was modified to enable significantly more efficient learning.

Using the blockwise Kronecker covariance for the University of Minnesota crowd anomaly dataset, it was found that state of the art anomaly detection performance was possible, and the use of temporal modeling and the sums of Kroneckers representation enabled significantly improved performance. In addition, good anomaly localization ability was observed.

Vii Acknowledgements

This research was partially supported by ARO under grant W911NF-11-1-0391 and by AFRL under grant FA8650-07-D-1220-0006.


  • [1] K. Greenewald, T. Tsiligkaridis, and A. Hero, “Kronecker sum decompositions of space-time data,” in Proc. of IEEE CAMSAP (to appear), available as arXiv 1307.7306, Dec 2013.
  • [2] M. J. Choi, V. Chandrasekaran, and A. S. Willsky, “Gaussian multiresolution models: Exploiting sparse markov and covariance structure,” Signal Processing, IEEE Transactions on, vol. 58, no. 3, pp. 1012–1024, 2010.
  • [3] “Unusual crowd activity dataset made available by the university of minnesota at:
  • [4] S. Wu, B. E. Moore, and M. Shah, “Chaotic invariants of lagrangian particle trajectories for anomaly detection in crowded scenes,” in Computer Vision and Pattern Recognition (CVPR), 2010 IEEE Conference on.   IEEE, 2010, pp. 2054–2060.
  • [5] D. Helbing and P. Molnar, “Social force model for pedestrian dynamics,” Physical review E, vol. 51, no. 5, p. 4282, 1995.
  • [6] V. Chandola, A. Banerjee, and V. Kumar, “Anomaly detection: A survey,” ACM Computing Surveys (CSUR), vol. 41, no. 3, p. 15, 2009.
  • [7] A. Wiesel, O. Bibi, and A. Globerson, “Time varying autoregressive moving average models for covariance estimation,” 2013.
  • [8] D.-Y. Chen and P.-C. Huang, “Motion-based unusual event detection in human crowds,” Journal of Visual Communication and Image Representation, vol. 22, no. 2, pp. 178–186, 2011.
  • [9] R. Mehran, A. Oyama, and M. Shah, “Abnormal crowd behavior detection using social force model,” in Computer Vision and Pattern Recognition, 2009. CVPR 2009. IEEE Conference on.   IEEE, 2009, pp. 935–942.
  • [10] V. Mahadevan, W. Li, V. Bhalodia, and N. Vasconcelos, “Anomaly detection in crowded scenes,” in Computer Vision and Pattern Recognition (CVPR), 2010 IEEE Conference on.   IEEE, 2010, pp. 1975–1981.
  • [11] V. Saligrama and Z. Chen, “Video anomaly detection based on local statistical aggregates,” in Computer Vision and Pattern Recognition (CVPR), 2012 IEEE Conference on.   IEEE, 2012, pp. 2112–2119.
  • [12] Y. Benezeth, P.-M. Jodoin, V. Saligrama, and C. Rosenberger, “Abnormal events detection based on spatio-temporal co-occurences,” in Computer Vision and Pattern Recognition, 2009. CVPR 2009. IEEE Conference on.   IEEE, 2009, pp. 2458–2465.
  • [13] B. Luvison, T. Chateau, J.-T. Lapreste, P. Sayd, and Q. C. Pham, “Automatic detection of unexpected events in dense areas for video surveillance applications,” 2011.
  • [14] T. Tsiligkaridis and A. Hero, “Covariance estimation in high dimensions via kronecker product expansions,” to appear IEEE Trans on SP in 2013, Also see arXiv 1302.2686, Feb 2013.
  • [15] T. Tsiligkaridis, A. Hero, and S. Zhou, “On convergence of kronecker graphical lasso algorithms,” IEEE Trans. Signal Proc., vol. 61, no. 7, pp. 1743–1755, 2013.
  • [16]

    G. I. Allen and R. Tibshirani, “Transposable regularized covariance models with an application to missing data imputation,”

    Annals of Applied Statistics, vol. 4, no. 2, pp. 764–790, 2010.
  • [17] K. Werner and M. Jansson, “Estimation of kronecker structured channel covariances using training data,” in Proceedings of EUSIPCO, 2007.
  • [18] N. Cressie, Statistics for Spatial Data.   Wiley, New York, 1993.
  • [19] J. Yin and H. Li, “Model selection and estimation in the matrix normal graphical model,”

    Journal of Multivariate Analysis

    , vol. 107, 2012.
  • [20] E. Bonilla, K. M. Chai, and C. Williams, “Multi-task gaussian process prediction,” in NIPS, 2007.
  • [21] K. Yu, J. Lafferty, S. Zhu, and Y. Gong, “Large-scale collaborative prediction using a nonparametric random effects model,” in ICML, 2009, pp. 1185–1192.
  • [22] Y. Zhang and J. Schneider, “Learning multiple tasks with a sparse matrix-normal penalty,” Advances in Neural Information Processing Systems, vol. 23, pp. 2550–2558, 2010.
  • [23] M. G. Genton, “Separable approximations of space-time covariance matrices,” Environmetrics, vol. 18, no. 7, pp. 681–695, 2007.
  • [24] K. Werner, M. Jansson, and P. Stoica, “On estimation of covariance matrices with kronecker product structure,” Signal Processing, IEEE Transactions on, vol. 56, no. 2, pp. 478–491, 2008.
  • [25] J. Kamm and J. G. Nagy, “Optimal kronecker product approximation of block toeplitz matrices,” SIAM Journal on Matrix Analysis and Applications, vol. 22, no. 1, pp. 155–172, 2000.
  • [26] A. Kannan, M. Ostendorf, W. C. Karl, D. A. Castañon, and R. K. Fish, “Ml parameter estimation of a multiscale stochastic process using the em algorithm,” Signal Processing, IEEE Transactions on, vol. 48, no. 6, pp. 1836–1840, 2000.
  • [27] H. Yu, J. Dauwels, X. Zhang, S. Xu, and W. I. T. Uy, “Copula gaussian multiscale graphical models with application to geophysical modeling,” in Information Fusion (FUSION), 2012 15th International Conference on.   IEEE, 2012, pp. 1741–1748.