Fabio Cuzzolin
email: fabiocuzzolin@gmail.com

Learning Pullback Metrics for Linear Models

In this paper we present an unsupervised differential-geometric approach for learning Riemannian metrics for dynamical models. Given a training set of models the optimal metric is selected among a family of pullback metrics induced by the Fisher information tensor through a parameterized diffeomorphism. The problem of classifying motions, encoded as dynamical models of a certain class, can then be posed on the learnt manifold. Experimental results concerning action and identity recognition based on simple scalar features are shown, proving how learning a metric actually improves classification rates when compared with Fisher geodesic distance and other classical distance functions.

Introduction

Manifold learning [START_REF] Bar-Hillel | Learning distance functions using equivalence relations[END_REF][START_REF] Bilenko | Integrating constraints and metric learning in semi-supervised clustering[END_REF][START_REF] Schultz | Learning a distance metric from relative comparisons[END_REF][START_REF] Tsang | Distance metric learning with kernels[END_REF][START_REF] Zhang | Learning metrics via discriminant kernels and multidimensional scaling: Toward expected euclidean representation[END_REF][START_REF] Eick | Using clustering to learn distance functions for supervised similarity assessment[END_REF] has become a popular topic in machine learning and computer vision in the last few years, as many objects of interests (like natural images, or sequences representing walking persons), in spite of their apparent high dimensionality, live in a non-linear space of reduced dimension. Many unsupervised algorithms (e.g. locally linear embedding [START_REF] Roweis | Nonlinear dimensionality reduction by locally linear embedding[END_REF]) take an input dataset and embed it into some other space, implicitly learning a metric. However, they fail to learn a full metric for the whole input space (even though extensions which take this issue into account have been recently formulated [START_REF] Bengio | Out-of-sample extensions for LLE, Isomap, MDS, eigenmaps, and spectral clustering[END_REF]). On the other side, approaches to successfully reduce metric learning to constrained least square optimization in the linear case have been proposed [START_REF] Xing | Distance metric learning with applications to clustering with side information[END_REF][START_REF] Shental | Adjustment learning and relevant component analysis[END_REF].

In particular, it makes a lot of sense to represent videos or image sequences as realizations of some sort of dynamical model, either stochastic (e.g. HMMs) or deterministic (e.g. ARMA). This approach already proved to be effective in contexts like video coding (e.g. dynamic textures [START_REF] Doretto | Editable dynamic textures[END_REF]). A number of metrics or distance functions for linear systems has been already introduced, in particular in the context of system identification: cepstrum distances [START_REF] Martin | A metric for ARMA processes[END_REF], subspace angles [START_REF] Cock | Subspace angles and distances between arma models[END_REF], gap metric [START_REF] Zames | Unstable systems and feedback: The gap metric[END_REF] and its variants nu-gap [START_REF] Vinnicombe | A ν-gap distance for uncertain and nonlinear systems[END_REF] and graph metric [START_REF] Vidyasagar | The graph metric for unstable plants and robustness estimates for feedback stability[END_REF], kernel methods [START_REF] Smola | Hilbert space embeddings in dynamical systems[END_REF]. Besides, a vast literature exists about dissimilarity measures between Markov models, variants of the Kullback-Leibler divergence [START_REF] Kullback | On information and sufficiency[END_REF].

Consider though the problem of classifying a dynamical model (as representative of an input image sequence). A simple mental experiment is enough to understand that no single distance function can possibly outperform all the others in each and every classification problem, since data-points can be endowed with different labeling while maintaining the same geometrical structure. In gait analysis, for instance, each image sequence possesses several different labels: the identity of the moving person, the action performed, the viewpoint (when several cameras are presented [START_REF] Gross | The CMU motion of body (Mobo) database[END_REF]), the emotional status, etc.

The most reasonable thing to do when possessing some a-priori information is then try and learn in a supervised fashion the "best" distance function for a specific classification problem. As even linear dynamical models live in a nonlinear space, the need for a principled way of learning Riemannian metrics from the data naturally arises. An interesting tool is provided by the formalism of pullback metrics. If the models belong to a Riemannian manifold M , any diffeomorphism of M onto itself induces such a metric on M . By designing a suitable family of diffeomorphisms depending on a parameter p we then obtain a family of pullback metrics on M we can optimize on.

In this paper we propose to apply this formalism to the case of linear dynamical models endowed with the classical Fisher metric. Analytical expressions of Fisher metric and related geodesics for several classes of linear dynamical systems have been provided by Hanzon [START_REF] Hanzon | Identifiability, recursive identification and spaces of linear dynamical systems[END_REF], Ravishanker [START_REF] Ravishanker | Differential geometry of ARMA models[END_REF], and Rijkeboer [START_REF] Rijkeboer | Fisher optimal approximation of an AR(n)-process by an AR(n-1)process[END_REF]. Given a dataset of dynamical models a "natural" Riemannian metric is learned by inverse volume element minimization on a family of pullback metrics induced by a class of diffeomorphisms. The learnt metric can later be used in any classification scheme. We argue that this improves the classification performance by posing the problem in the region where the systems actually live. We apply this learning scheme to action recognition and identity recognition from gait. We use the image sequences collected in the Mobo database [START_REF] Gross | The CMU motion of body (Mobo) database[END_REF] to show experiments in which (as a reference) simple nearest-neighbor classifiers based on pullback Fisher metric between stochastic models outperform analogous classifiers based on classical a-priori distances.

Learning non-linear manifolds of dynamical models

Metric and manifold learning. Learning distance functions or metrics to improve, for instance, classification performance is a sensible approach when some a-priori information is available [START_REF] Bar-Hillel | Learning distance functions using equivalence relations[END_REF][START_REF] Bilenko | Integrating constraints and metric learning in semi-supervised clustering[END_REF][START_REF] Schultz | Learning a distance metric from relative comparisons[END_REF][START_REF] Tsang | Distance metric learning with kernels[END_REF][START_REF] Zhang | Learning metrics via discriminant kernels and multidimensional scaling: Toward expected euclidean representation[END_REF][START_REF] Eick | Using clustering to learn distance functions for supervised similarity assessment[END_REF]. A natural optimization criterion consists on maximizing the classification performance achieved by using the learnt metric in a simple nearest-neighbor classifier. Xing et al. [START_REF] Xing | Distance metric learning with applications to clustering with side information[END_REF] have recently proposed a way to solve this optimization problem for linear maps y = A 1/2 x, when some pairs of points are known to be "similar". This leads to a parameterized family of Mahalanobis distances x -y A . Given the linearity of the map the problem of minimizing the squared distance between similar points turns out to be convex, hence solvable with efficient numerical algorithms. Shental et al. [START_REF] Shental | Adjustment learning and relevant component analysis[END_REF] have successfully faced the same problem in a linear framework, by posing an optimization problem based on information theory (relevant component analysis). However, in several real-world applications the data are inherently nonlinear. This is exactly the case of dynamical models (and human motions, in particular). The need for a way of learning Riemannian metrics from such data arises.

Pullback metrics. Some notions of differential geometry provide us indeed with a tool for building a structured family of metrics on which to define an optimization problem, the basic ingredient of a metric learning algorithm. The idea is the following. Suppose your dataset already lives on a Riemannian manifold M of some sort: a Riemannian metric is defined in any point of the manifold. Any diffeomorphism (a differentiable map) from M to itself induces a new metric, called "pullback metric" (see below). If we manage to define an entire class of diffeormorphisms depending on some parameter λ we then get the corresponding family of pullback metrics, also depending on λ. We can then define an optimization problem over this family of metrics in order to select an "optimal" metric, which in turn determines the desired manifold. Of course, the nature of this manifold depends on the objective function we choose to optimize. Following Lebanon [START_REF] Lebanon | Metric learning for text documents[END_REF] we propose to maximize the inverse volume of the space around the dataset to select the manifold on which the latter lives.

Proposed methodology. In this paper the dataset is composed by dynamical models representing motions. Many classes of dynamical systems do live in a Riemannian manifold M associated with the Fisher information metric [START_REF] Amari | Differential geometric methods in statistics[END_REF]. We can then apply the formalism of pullback metrics, and learn a manifold associated with a training set D of movements. When new motions are acquired they can then be classified on this reduced space within a three-step procedure:

1. each sequence of measurements representing a dynamical data-set is mapped into the parameters of a dynamical model describing the sequence by parameter identification; 2. a parametric metric is learned for the resulting set of dynamical systems; 3. standard classifiers (e.g. k-nearest neighbor, SVM, etc.) are used to classify the systems according to the new metric.

Volume minimization for pullback metric learning

Consider a family of diffeomorphisms between the Riemannian manifold M in which the dataset D = {m 1 , ..., m N } ⊂ M resides and itself:

F p : M → M , m → F p (m), m ∈ M . Let us denote by T m M the tangent space to M in m.
Any such diffeomorphism F is associated with a push-forward map

F * : T m M → T F (m) M v ∈ T m M → F * v ∈ T F (m) M (1)
defined as

F * v(f ) = v(f • F
) for all the smooth functions f on M (see Figure 1). Consider now a Riemannian metric g : T M × T M → R on M . Roughly speaking, a Riemannian metric determines how to compute scalar products of tangent vectors v ∈ T m M . The map F induces a pullback metric on M as

g * m (u, v) . = g F (m) (F * u, F * v). (2)
Scalar products of two vectors u, v according to g * m are computed as the scalar product with respect to the original metric g of the pair of vectors F * u, F * v which are mapped onto u, v by F . Any parametric family of differentiable maps naturally generates an entire parametric family of metrics on M . The pullback geodesic between two points is then the geodesic connecting their images with respect to the original metric.

Volume element minimization. From this parameterized family of metrics, we wish to determine the Riemannian manifold which "best" represents the dataset, according to some criterion. We propose here to minimize the volume element associated with a metric, as suggested by G. Lebanon [START_REF] Lebanon | Metric learning for text documents[END_REF] in the context of document retrieval (see also [START_REF] Murray | Differential Geometry and Statistics[END_REF]). Equivalently, we seek to maximize

O(D) = N k=1 (det g(m k )) -1 2 M (det g(m)) -1 2 dm (3)
where g(m k ) denotes the Riemannian metric in the point m k of the data-set D living on a Riemannian manifold M . This amounts to finding a lower dimensional representation of the dataset, in a similar fashion to locally linear embedding [START_REF] Roweis | Nonlinear dimensionality reduction by locally linear embedding[END_REF] or laplacian eigenmaps [START_REF] Belkin | Semi-supervised learning on riemannian manifolds[END_REF], where dimensionality reduction is often considered a factor in improving classification.

Computing the Gramian. The computation of (3) requires that of the Gramian det g. To find the expression of the Gramian associated with the pullback metric (2) we first need to choose a base of the tangent space T m M to M . Let us then denote with {∂ i }, i = 1, ..., dim M the base of T m M . We can compute the push-forward of the vectors of this base, getting in turn a base for T F (m) M . By definition, the push-forward of a vector v ∈ T m M is [START_REF] Lebanon | Metric learning for text documents[END_REF]

F p (v) . = d dt F p (m + t • v) t=0 , v ∈ T m M. (4)
The diffeomorphism F p induces a base for the space of vector fields on M , w i . = {F p (∂ i )}, for i = 1, ..., dim M . We can rearrange these vectors as rows of a matrix

J = [w 1 ; • • • ; w dim M ]. ( 5
)
The volume element of the pullback metric g * in a point m ∈ M is then the determinant of the Gramian [START_REF] Lebanon | Metric learning for text documents[END_REF]:

det g * (m) . = det[g(F * p (∂ i ), F * p (∂ j ))] ij = det(J T gJ).
If J is a square matrix (as in the following) we get simply

det g * (m) = det(J) 2 • det g(m). (6)
Plugging ( 6) in ( 3) we obtain the function to minimize. Algorithm. We have then a method to select an optimal metric for a dataset D = {m 1 , ..., m N } of points of a Riemannian manifold M with (basis) metric g.

1. First, a family {F p , p ∈ P } of diffeomorphisms (P the parameter space) from M to itself is designed to provide a large search space of metrics (the variable in our optimization scheme) from which to select the optimal one; 2. F p induces a family g * p of pullback metrics (2) on M ; 3. we can then pose an optimization problem max p∈P O(p), where

O(p) = N k=1 (det g * p (m k )) -1 2 M (det g * p (m)) -1 2 dm (7)
is the objective function, the inverse volume of the manifold in the neighborhood of the N points {m k , k = 1, ..., N } of the dataset D; 4. this yields an optimal pullback metric ĝ * ; 5. knowing the geodesics of M we can compute the distance between two points according to ĝ * ; 6. this distance function can then be used to cluster or classify the dataset.

In the following Section we will show how to apply this optimization scheme to the case of linear dynamical models, endowed with the Fisher metric. In particular we will consider the simple case of autoregressive models of order 2, for which both Fisher metric and geodesics are analytically known, and recover closed-form expressions for the objective function [START_REF] Roweis | Nonlinear dimensionality reduction by locally linear embedding[END_REF].

Learning pullback metrics for linear models

As linear dynamical models live in a Riemannian space, we can apply to them the pullback metric formalism and obtain a family of metrics on which to optimize.

Fisher metric. The study of the geometrical structure of the space formed by a family of probability distribution was first due to Rao, and was developed by Nagaoka and Amari [START_REF] Amari | Differential geometric methods in statistics[END_REF]. A family S of probability distributions p(x, ξ) depending on a n-dimensional parameter ξ can be regarded in fact as an n-dimensional manifold. If the Fisher information matrix

g ij . = E ∂logp(x, ξ) ∂ξ i ∂logp(x, ξ) ∂ξ j (8)
is nondegenerate, G = [g ij ] is a Riemannian tensor, and S is a Riemannian manifold. In particular, the analytic expressions of the entries of the Fisher information matrix for several manifolds of linear MIMO systems have been obtained by Hanzon et al. [START_REF] Hanzon | Aspects of Fisher geometry for stochastic linear systems, problem 25[END_REF].

An interesting class of linear dynamical models. Let us consider in particular the class of stable autoregressive discrete-time processes of order 2, AR(2), in a stochastic setting in which the input signal is a Gaussian white noise with zero mean and unit variance. This set can be given a Riemannian manifold structure under Fisher metric. A natural parametrization uses the nonunit coefficients (a 1 , a 2 ) of the denominator of the transfer function, h(z) = z 2 /(z 2 + a 1 z + a 2 ) (which corresponds to the AR difference y(k) = -a 1 y(k -1) -a 2 y(k -2)). To impose stability the necessary conditions are 1 + a 1 + a 2 > 0, 1 -a 1 + a 2 > 0, and 1 -a 2 > 0. The manifold is then composed by a single connected component (see Figure 2-left). The Riemannian Fisher tensor can be expressed, in alternative local coordinates given by the Schur parameters γ 1 = a 1 /(1 + a 2 ), γ 2 = a 2 , as [START_REF] Peeters | On the Riemannian manifold structure of classes of linear systems[END_REF]

[0,-1] [2,1] [-2,1]
g(γ 1 , γ 2 ) = 1 (1 -γ 2 2 ) (1+γ2) 2 (1-γ 2 1 ) 0 0 1 . (9)
Geodesics. To compute the distance between two points of a Riemannian manifold (and in particular two dynamical models) it is not sufficient to know the metric: it is necessary to compute (analytically or numerically) the shortest path connecting them on the manifold (geodesic). All the geodesics of stable AR(2) systems endowed with the Fisher metric (9) as a function of the Schur parameters have been analytically computed by Rijkeboer [START_REF] Rijkeboer | Fisher optimal approximation of an AR(n)-process by an AR(n-1)process[END_REF].

A family of diffeomorphisms for AR(2). To build a parameterized family of Riemannian metrics for AR(2) we can apply the optimal pullback metric scheme of Section 2 to the case in which the dataset D is a collection of linear systems which leaves in a Riemannian manifold M = AR(2), the triangle of Figure 2. It is then necessary to choose a family of diffeomorphisms of M onto itself. Clearly the choice of a class of diffeomorphisms depends on the class of dynamical systems we adopt, but it is not univocal. The more sophisticated the set of diffeomorphisms, the larger is the search space to optimize the metric on. Some hints can be provided by the geometry itself of the space M of dynamical models. As a matter of fact one possible choice for a diffeomorphism of AR(2) onto itself is suggested by the triangular form of the manifold (see Figure 2),

F p (m) = F p ([m 1 , m 2 , m 3 ]) = [λ 1 m 1 , λ 2 m 2 , λ 3 m 3 ] λ • m (10)
where p = λ = [λ 1 , λ 2 , λ 3 ] with λ 1 +λ 2 +λ 3 = 1, while m = [m 1 , m 2 , m 3 ] collects the simplicial coordinates of a system a ∈ AR(2) in the manifold:

a = [a 1 , a 2 ] ′ = m 1 [0, -1] ′ + m 2 [2, 1] ′ + m 3 [-2, 1] ′ (11)
and λ • m denotes the scalar product of the two vectors. The application [START_REF] Shental | Adjustment learning and relevant component analysis[END_REF] stretches the triangle towards the vertex with the largest λ i (Figure 2-right).

Volume element for AR [START_REF] Bilenko | Integrating constraints and metric learning in semi-supervised clustering[END_REF]. It is then possible to apply the method of Section 3 to find the analytical expression of the determinant of the Gramian detg * λ (λ, m) (6), as a function of the parameter λ of the family of diffeomorphisms [START_REF] Shental | Adjustment learning and relevant component analysis[END_REF]. By plugging it in the expression for the inverse volume (3) we finally obtain the objective function O(λ) to optimize to find an optimal metric. Theorem 1. The volume element of the pullback metric on AR(2) induced by the diffeomorphism ( 10) is given by

detg * λ (λ, m) ∝ (λ 1 λ 2 λ 3 ) 2 (λ • m) 6 • 1 m 2 1 m 2 m 3 . ( 12
)
Proof. Let us choose as base of the tangent space in AR(2) the set

∂ 1 = [1/2, 1/2] ′ , ∂ 2 = [-1/2, 1/2] ′ .
To compute the Gramian we need to express the diffeomorphism with respect to a. From Equation ( 11)

a 1 = 2(m 2 -m 3 ), a 2 = m 2 + m 3 -m 1 , while the inverse relation is m 2 = 1+a1+a2 4 , m 3 = 1-a1+a2 4 , m 1 = 1-a2
2 . Hence F λ can be expressed in terms of a 1 , a 2 as

F λ (a) = 1 ∆ 2λ 2 (1 + a 1 + a 2 ) -2λ 3 (1 -a 1 + a 2 ), λ 2 (1 + a 1 + a 2 ) + λ 3 (1 -a 1 + a 2 ) -2λ 1 (1 -a 2 ) ′ ( 13
)
where

∆ = 2λ 1 (1 -a 2 ) + λ 2 (1 + a 1 + a 2 ) + λ 3 (1 -a 1 + a 2 ), so that F λ (a + tv) = [2λ 2 (1 + a 1 + tv 1 + a 2 + tv 2 ) -2λ 3 (1 -a 1 -tv 1 + a 2 + tv 2 ), λ 2 (1 + a 1 + tv 1 + a 2 + tv 2 ) + λ 3 (1 -a 1 -tv 1 + a 2 + tv 2 ) -2λ 1 (1 -a 2 -tv 2 )] ′ . We can then compute 4 d dt F λ (m + t • v) t=0 ,

and in particular

w 1 = d dt F λ (m + t∂ 1 ) t=0 = [ 2λ 1 λ 2 (3 -a 2 + a 1 ) + 2λ 3 (2λ 2 -λ 1 )(1 -a 1 + a 2 ), 2λ 1 λ 2 (3 -a 2 + a 1 ) + 2λ 1 λ 3 (1 -a 1 + a 2 )] (14) while w 2 = d dt F λ (m + t∂ 2 ) t=0 = [ -2λ 1 λ 3 (3 -a 2 + a 1 ) + 2λ 2 (λ 1 -2λ 3 )(1 + a 1 + a 2 ), 2λ 1 λ 3 (3 -a 2 + a 1 ) + 2λ 1 λ 2 (1 + a 1 + a 2 )]. ( 15
)
The determinant of the matrix J = [w 1 ; w 2 ] ( 5) is then (after a few passages)

det J = 32 λ 1 λ 2 λ 3 ∆ 3 = 1 2 λ 1 λ 2 λ 3 (λ • m) 3 . ( 16
)
Finally, the function (3) to maximize assumes the form

O(λ) = N k=1 (λ • m k ) 3 AR(2) (λ • m) 3 m 1 √ m 2 m 3 dm (17)
where the normalization factor

I(λ) = AR(2) (λ • m) 3 m 1 √ m 2 m 3 dm forbids triv-
ial solutions in which the volume element is minimized by shrinking the whole space. It can be computed by decomposing (λ • m) 3 using Tartaglia's formula:

I(λ) = c1+c2+c3=3 3! c 1 !c 2 !c 3 ! • 3 j=1 λ cj j AR(2) m 1+c1 1 m 1/2+c2 2 m 1/2+a3 3 dm. ( 18
)
The objective function ( 17) can be maximized by means of any numerical optimization scheme, simple (like gradient descent) or more sophisticated.

Extension to multi-dimensional systems.

It is important to stress that the proposed scheme is by no means limited to scalar dynamical models. It is true that the analytic expressions for the Fisher metric and the associated geodesics are known for scalar outputs. This allows to derive elegant analytic expressions for the inverse volume [START_REF] Smola | Hilbert space embeddings in dynamical systems[END_REF] to maximize. However, in the case of multi-dimensional linear systems both Fisher metric and its geodesics can still be computed by means of an iterative numerical scheme [START_REF] Peeters | On the Riemannian manifold structure of classes of linear systems[END_REF][START_REF] Hanzon | Aspects of Fisher geometry for stochastic linear systems, problem 25[END_REF]. The extension of the pullback manifold learning scheme to multi-dimensional systems is then straightforward. In any case, using the Fisher information metric as basis Riemannian metric is not mandatory. We can as easily adopt the standard Euclidean metric as initial distance, and build a family of pullback Euclidean metrics to optimize upon.

Effect on classification performances

It can be argued that maximizing classification performance and minimizing volume elements are indeed correlated optimization problems. In other words, that classification is easier on the pullback manifold which better interpolates the dataset of dynamical systems.

Experiments on human motion classification. To test this conjecture we considered the problem of recognizing actions and identities from image sequences. We used the Mobo database [START_REF] Gross | The CMU motion of body (Mobo) database[END_REF], a collection of 600 image sequences of 25 people walking on a treadmill in four different variants (slow walk, fast walk, walk on a slope, walk carrying a ball), seen from 6 different fixed viewpoints (see Figure 3-middle). Each sequence of the database possesses three labels: action, view, and identity. As the objective function to optimize has been computed in Section 4 for single input-single output linear systems, we extracted from each image a single scalar feature, the width of the minimum box containing the silhouette (Figure 3-left). Each scalar sequence was passed as input to an identification algorithm which generated a dataset of AR(2) linear systems, one for each labeled motion sequence. However, as we just pointed out, this is not a limitation of the proposed approach.

We empirically assessed the conjecture on the effect of pullback manifold learning on classification by measuring the performance of a nearest-neighbor classifier based on the optimal pullback metrics induced by the diffeomorphism [START_REF] Shental | Adjustment learning and relevant component analysis[END_REF] and comparing the results with those of NN classifiers based on several apriori distances between models (including the Fisher geodesic distance itself). Identity recognition. In a first experiment we selected a training set of models, and used the pullback geodesic distance to classify the identity of the person appearing in a different set of randomly selected sequences. This is a very difficult problem, as there are 25 different people, and the one-dimensional signal we chose to represent sequences with (the series of bounding box widths along the image sequence) clearly provides insufficient information. However, measuring the comparative performance of the metrics can be useful to see how learning appositively a metric improves classification. We implemented a naive Frobenius norm of the system matrices in canonical form, gap and ν-gap metrics [START_REF] Zames | Unstable systems and feedback: The gap metric[END_REF][START_REF] Vinnicombe | A ν-gap distance for uncertain and nonlinear systems[END_REF], subspace angles [START_REF] Cock | Subspace angles and distances between arma models[END_REF][START_REF] Martin | A metric for ARMA processes[END_REF], and of course the Fisher geodesic distance together with the associated optimal pullback metric. Figure 3-right shows the average percentage of correctly classified testing systems over several dozens of runs in which we randomly selected an increasing number of systems in both training and testing set. The number of training sequences ranged from 10 to 50, while the size of the testing set would vary from 50 to 300. The purpose of the random selection was avoiding biased results due to a particular choice of training and testing sets. Standard variations were very small and they are not reported. As a comparison we also tested a second diffeomorphism for AR(2) systems, namely

F λ (m) = [λm 1 + (1 -λ)m 2 , λm 2 + (1 -λ)m 3 , λm 3 + (1 -λ)m 1 ], (19)
with 0 < λ < 1. Figure 3-right shows that the larger the parameter space (as ( 19) has just one parameter, while [START_REF] Shental | Adjustment learning and relevant component analysis[END_REF] has two) the greater the chance of selecting a better metric. All other distances have rather similar performances as they obviously have not been designed to solve classification tasks but for identification purposes. In the following only the performance of the second best distance is reported for comparison. Action recognition. In a second series of experiments we used the pullback geodesic distance to classify the action performed by the person. The correct classification rate was measured after several dozens of repeated trials in which training and testing set of increasing size (as in the ID experiment) were randomly selected from the database in order to avoid any bias in the results. Figure 4 compares again the average recognition rate over such repetitions of the optimal pullback metric for autoregressive systems (in particular under deformation map [START_REF] Gross | The CMU motion of body (Mobo) database[END_REF]) and its best competitor among a-priori metrics (typically Frobenius norm or cepstrum distance), as a function of testing and training size, 10 repetitions for each pair of such sizes. Rates are not very high, but remember that actions in the Mobo database are just slightly different variations of the walking gait, and we were using a single scalar feature. We want to stress relative performances here. To test the approach more thoroughly we also conducted six separate experiments by selecting the portion of the dataset associated with a single view, for each possible view. Figure 5 illustrates the average performance of the classifiers associated with the optimal pullback metric for AR(2) models with (this time) diffeomorphism [START_REF] Shental | Adjustment learning and relevant component analysis[END_REF] and the best competing non-learnt, a-priori distance (gap, ν-gap, etc.) for four view-dependent experiments, as a function of the testing set size. As usual the average (for each size) is computed over 10 repeated random selections. The optimal pullback metric performs far better than all the others.

Conclusions

In this paper we proposed a differential-geometric framework for manifold learning given a data-set of linear dynamical models. We pose an optimization problem in which the pullback metric induced by a diffeomorphism which minimizes the volume element around the available data is learned. We adopt as basis metric tensor the classical Fisher information matrix. This yields a global embedding, while usual spectral methods only provide images of training points. We showed experimental results concerning identity and action recognition, which support the claim that learning a metric in such a way actually improves classification performance with respect to the Fisher metric and all other known a-priori distances between dynamical models.

Fig. 1 .

 1 Fig. 1. The push-forward map associated with a diffeomorphism on a Riemannian manifold M .

Fig. 2 .

 2 Fig. 2. Left: Graphical representation of the manifold of stable autoregressive systems of order 2, AR(2), with the non-unit coefficients of the denominator of h(z) as parameters. It forms a simplex with vertices [-2, 1], [2, 1], [0, -1]. Right: Effect of a diffeomorphism of the form (10) on the AR(2) simplex.

Fig. 3 .

 3 Fig. 3. Left: from each image the width of the bounding box containing the silhouette is extracted. Middle: location and orientation of the six cameras in the Mobo experiment (the origin of the frame is roughly in the position of the walking person on the treadmill). Right: Performance of the NN classifier in the ID recognition experiment, based on several a-priori distances between dynamical models and two pullback metrics induced by the Fisher geometry. Classification rate for the used metrics averaged on the size of the training set. 1 -Frobenius norm of the matrices in canonical form; 2gap metric; 3ν-gap matric; 4 -subspace angles; 5 -basis Fisher metric; 6 -pullback metric with diffeomorphism (19); 7 -pullback metric with diffeomorphism (10).

Fig. 4 .

 4 Fig. 4. Average recognition performance of second-best a-priori distance (left, in this case the Frobenius distance between the matrices representing two systems) and optimal pullback metric for AR(2) systems and diffeomorphism (19) (middle) in the action recognition experiment. Here the whole dataset is considered, regardless of the viewpoint. Rates are plotted against the size of testing (x axis) and training (y axis) set. Right: recognition rates for view 5 only. The learnt pullback metric always outperforms its best competitor.

Fig. 5 .

 5 Fig.5. Correct action classification rates for sequences coming from three different viewpoints[START_REF] Bar-Hillel | Learning distance functions using equivalence relations[END_REF][START_REF] Zhang | Learning metrics via discriminant kernels and multidimensional scaling: Toward expected euclidean representation[END_REF][START_REF] Eick | Using clustering to learn distance functions for supervised similarity assessment[END_REF]. Sequences are represented as autoregressive order 2 systems. The performance (in red) of the pullback metric associated with the map (10) is shown outperforming the best competitor (in blue, typically Frobenius norm or Martin's distance) for increasing sizes of the testing set (from 20 to 50).


Michalis Raptis
email: mraptis@cs.ucla.edu

Kamil Wnuk
email: kwnuk@cs.ucla.edu

Stefano Soatto
email: soatto@cs.ucla.edu

Flexible Dictionaries for Action Classification

We present a simple approach to action classification which constructs a vector quantization of primitive motions from time series data corresponding to relative limb position estimates. The temporal scale, mean, and shape of primitive motion trajectories are independently modeled, thus creating a flexible dictionary of action primitives. We then explore two inference techniques that leverage our action dictionary representation, and evaluate their performance on both motion capture and video benchmark data. Our results indicate that even simplistic algorithms can outperform significantly more sophisticated ones in existing benchmark datasets.

Introduction

Classification of human activity has become a core interest domain in recent years, with a wide range of applications including surveillance, human-machine interfaces and video indexing. Human actions can be particularly difficult to classify due to potentially large variability in the time dimension in addition to the usual difficulties associated with processing of images. However, the temporal evolution contains a significant amount of information, as illustrated by [START_REF] Johansson | Visual perception of biological motion and a model for its analysis[END_REF]. In this work, we focus on classification of events with distinct temporal signatures.

Our core hypothesis is that the most complex human actions can be decomposed into short primitive actions. A simple, yet representative, example is the motion of the legs during walking. This motion can be decomposed into two components: one when the foot is in contact with the ground and the other when it is swinging above the ground [START_REF] Mcgeer | Passive Dynamic Walking[END_REF]. The decomposition of actions into elementary motions has been explored previously in the literature. Marr and Vaina [START_REF] Marr | Representation and Recognition of the Movements of Shapes[END_REF] proposed a method to segment movements into pieces, [START_REF] Bissacco | Recognition of human gaits[END_REF] defined primitive actions in terms of linear dynamical systems, and [START_REF] Bregler | Learning and recognizing human dynamics in video sequences[END_REF] proposed a probabilistic action decomposition framework.

In contrast to some recent work that reasons in the spatio-temporal domain, we do not deal with pixel or feature-based appearance during motions deemed to be interesting [START_REF] Niebles | Unsupervised learning of human action categories using spatial-temporal words[END_REF][START_REF] Nowozin | Discriminative subsequence mining for action classification[END_REF][START_REF] Dollar | Behavior recognition via sparse spatio-temporal features. Visual Surveillance and Performance Evaluation of Tracking and Surveillance[END_REF]. Instead we focus on the source of such appearance variations: the deformation of the human body. However, we do not seek to explicitly identify dynamical systems that drive a particular human skeleton [START_REF] Bissacco | Recognition of human gaits[END_REF]. Instead we seek to learn primitives of the motion of human limbs through examples, and perform classification decisions using a dictionary-based representation.

In an attempt to capture the available temporal information, we represent actions as multi-dimensional time series. We capitalize on this representation to easily extract trajectories of simple body movements. By clustering such primitive trajectories, we seek to construct a nonparametric quantization of the physically possible primitive motions, as they appear in the actions that we are interested in classifying. Because primitive limb motions contain a large amount of variability in temporal scale (speed of the action), as well as their relative location with respect to the center of the body, we explicitly and independently quantize these informative components to more effectively represent and match primitive trajectories.

We pursue the above goals by constructing dictionaries of action primitives. The dictionaries are obtained by clustering windows extracted at interesting points in our time series representation of actions. We demonstrate the performance of our dictionaries by applying two types of inference techniques. Our results show that:

1. Both of our inference techniques discriminate actions with well above 90% accuracy on commonly used benchmark datasets, using only a simple lowdimensional parameterization of human motion that we can track over time in a fully automatic fashion. 2. When exact pose is unavailable, our approach is robust to using a very rough approximation (bounding box), allowing us to succeed where other pose-based classifiers require supervision in the pose tracking procedure. 3. When exact pose information is available, such as in motion capture data, we show that we outperform recent dynamics-based approaches by up to 8.33% mean classification accuracy.

Related Work

Human action recognition algorithms can be separated into two main classes. The first one is the holistic approach, where a collection of statistics is extracted from a video sequence and used for classification [START_REF] Zelnik-Manor | Statistical analysis of dynamic actions[END_REF][START_REF] Blank | Actions as Space-Time Shapes[END_REF][START_REF] Niebles | A hierarchical model of shape and appearance for human action classification[END_REF][START_REF] Dollar | Behavior recognition via sparse spatio-temporal features. Visual Surveillance and Performance Evaluation of Tracking and Surveillance[END_REF][START_REF] Efros | Recognizing action at a distance[END_REF][START_REF] Jhuang | A biologically inspired system for action recognition[END_REF]. The majority of such works represent a video sequence as spatio-temporal words [START_REF] Niebles | Unsupervised learning of human action categories using spatial-temporal words[END_REF][START_REF] Ikizler | Human action recognition using distribution of oriented rectangular patches[END_REF] without taking into account causal constraints. Consequently, they lack the ability to discriminate between actions that share common basic movements in different order. In an attempt to incorporate the order of movements, Nowozin et al. [START_REF] Nowozin | Discriminative subsequence mining for action classification[END_REF] encoded the relative temporal sequence of such spatio-temporal features by finding a discriminative subsequence, again without regard to how the features were generated.

The second general approach is model-based, where the motion is characterized by the parameters of a model that is fit to the data. Hidden Markov Models [START_REF] Wilson | Parametric hidden markov models for gesture recognition[END_REF][START_REF] Song | Towards detection of human motion[END_REF] and finite-state models [START_REF] Ikizler | Searching video for complex activities with finite state models[END_REF][START_REF] Hongeng | Video-based event recognition: activity representation and probabilistic recognition methods[END_REF] have been used to model the temporal variability of human motion. However, these methods lose the details necessary to infer the dynamics as a result of the coarseness of their representation. Others model the dynamics of human gaits using hybrid linear models [START_REF] Pavlovic | Learning switching linear models of human motion[END_REF][START_REF] Vidal | An algebraic geometric approach to the identification of a class of linear hybrid systems[END_REF][START_REF] Agarwal | Tracking articulated motion using a mixture of autoregressive models[END_REF][START_REF] Bissacco | Classifying human dynamics without contact forces[END_REF][START_REF] Li | Simultaneous learning of nonlinear manifold and dynamical models for high-dimensional time series[END_REF], which exhibit great generative power; however, their discriminative power decreases with increased model complexity. This tradeoff makes it difficult to classify actions such as dancing.

Others in the second class of approaches propose to design invariants of the motion sequence [START_REF] Ali | Chaotic invariants for human action recognition[END_REF]. However, these methods are difficult to generalize and require many samples of each action.

Our work attempts to capitalize on the best of both techniques. With our flexible dictionaries of action primitives we create a nonparametric model of possible human motions. The unique element of our approach is that such a model allows us to use statistical inference techniques typically reserved for holistic action approaches.

Actions as Time Series

As noted in our introductory statements, in this work we take a less common view of actions. Instead of representing an action as a collection of spatio-temporal patches [START_REF] Niebles | Unsupervised learning of human action categories using spatial-temporal words[END_REF][START_REF] Nowozin | Discriminative subsequence mining for action classification[END_REF][START_REF] Dollar | Behavior recognition via sparse spatio-temporal features. Visual Surveillance and Performance Evaluation of Tracking and Surveillance[END_REF], we interpret human motion as a multi-dimensional time series that captures the deformations of the body during activity. In particular, we track the positions of limb endpoints relative to the center of the body. As in [START_REF] Ali | Chaotic invariants for human action recognition[END_REF], these include the arms, legs, and head (only 5 points).

In the case of motion capture data, where positions are expressed in three dimensions, our chosen representation produces a 15 dimensional time series (x, y, z for each limb endpoint). For video we do not attempt to estimate 3D pose. Instead we roughly capture body deformations in the image plane. Because extraction of time series from video can be dataset dependent, details are provided below in the experimental section.

It is important to note that our methods are not bound to the above representation and can be used with time series generated from video or motion capture data using any number of techniques.

Dictionary of Action Primitives

Throughout this paper we define an action primitive as a time series subsequence that encodes a single dimension of a commonly occuring deformation. Due to our chosen method for obtaining time series, in our case such a primitive actually corresponds to a projection of a simple limb trajectory onto a single axis. For example, for a 3D trajectory of a single limb the time series resulting from the x, y and z coordinates each have their own action primitives. Independent action primitive dictionaries are constructed for all dimensions. This, for example, means 15 dictionaries in the case of motion capture data.

In the following sections we explore three incrementally more complex dictionary building procedures and their associated strengths and weaknesses for particular tasks. Each successive technique incorporates an additional degree of invariance by explicitly modeling a particular informative characteristic of action primitives.

Construction of a Simple Dictionary

The most basic approach to dictionary construction is to simply extract and cluster time series sub-sequences of a fixed size. Such clustering must be done with care, as recent literature on time series classification has shown that failing to tailor the sub-sequence sampling procedure to the signal can introduce artifacts that overpower the input, producing sinusoids independent of the underlying data [START_REF] Lin | Clustering of streaming time series is meaningless[END_REF][START_REF] Fujimaki | Theoretical analysis of subsequence time-series clustering from a frequency-analysis viewpoint[END_REF]. For this reason we sample the signal by extracting sub-windows only around extrema points, which achieves local translation invariance as illustrated in [START_REF] Soatto | On the distance between non-stationary time series[END_REF]. To counter the effects of noise when detecting extrema, a small initial smoothing is applied. Sampling at the extrema of all the time series captures all significant changes of motion direction within each dimension of an action.

Once sub-windows are extracted for all dimensions in all action samples, we quantize them independently for each dimension using K-means. The L 1 distance was used to compare sub-windows. As shown in Fig. 1(a), the representative power of this simple approach suffers because several dictionary elements are wasted on encoding the same primitive shape occurring at different amplitudes. Thus, in our next approach we seek to increase representative power by independently quantizing the mean of the sub-windows.

Increasing Representative Power by Mean Quantization

To achieve clustering of trajectories invariant to their average amplitude, we subtract the mean from all extracted sub-windows before quantizing them as above. The average amplitude of a primitive trajectory is likely to be informative for discriminating actions, as it captures the relative position of the limb during motion. Thus we quantize the subtracted means. Hence, in addition to one dictionary of primitive trajectories per dimension of the time-series, we also add quantization of the mean. Each window extracted from the time series composing an action can now be represented with a best matching canonized trajectory and an assignment to a particular bin of the quantized mean.

Fig. 1(b) shows that this approach generates a dictionary spanning a greater variability of trajectories, thus improving its representational and discriminative power. This dictionary construction approach was found to perform well in classification tasks. However, modeling the temporal scale of dictionary elements still proved essential for achieving best results.

Modeling Time Scale for Temporal Flexibility

Human actions can have large temporal variability. Even within a single action instance the same primitives may be performed at greatly varying speeds. Picking up an object, for example, may be performed very rapidly when the object is solid or very slowly when the object is fragile and caution is required. Therefore, to increase the accurracy with which the dictionary can capture our data it is Fig. 1. The action primitive dictionaries constructed without, 1(a), and with, 1(b), mean quantization for the x-coordinate of the left hand in the FutureLight motion capture dataset. The dictionaries were constructed with K-means, using K = 20, and a sub-window size of 21. Notice that in 1(a) several clusters (elements 2, 5, 6, and 7 for example) capture the same trajectory, but at different mean locations. In addition to being inefficient, this means that when we represent a signal with this dictionary, much of the effort will go to matching the mean, diminishing the significance of the trajectory (shape). In 1(b), the mean was quantized (lower right) into 5 bins independent of the trajectory quantization. One can see, by comparing the two figures, that once the mean was modeled out from the trajectory quantization much more variability was captured in the 20 trajectory elements.

important to be able to represent action primitives at different time scales. For this reason, we extend our sub-window extraction procedure.

As before, the centers of time series sub-windows are selected at extrema points. However, instead of using a preset size we automatically select the size of each sub-window. This is done by symmetrically growing window borders until the nearest extremum is encountered.

Once the length of all sub-windows is known, we quantize this temporal scale using K-means. Each extracted sub-window is then matched to the closest quantized window size and resampled to this canonical size. Interpolation with B-splines was used to minimize loss during sub-window resampling. By this quantization of the time scale, we avoid comparing large scale phenomena with small scale trajectories, while maintaining invariance within local scale ranges. This procedure of rescaling should be understood in the context of achieving invariance with respect to reparameterization of the time axis, as an alternative to performing dynamic time warping when comparing two time series. However, both the extraction of invariants and comparing modulo reparameterization of the temporal axis should be performed in a way that respects the dynamics of the underlying signal [START_REF] Soatto | On the distance between non-stationary time series[END_REF].

Following the rescaling, the mean quantization is performed as before, across all the extracted sub-windows. Finally, the clustering of the shape of the trajectories is performed independently within each canonical window size, thus creating a dictionary of action primitives spanning multiple time scales.

Adaptive window size selection does not only offer us a way of modeling the temporal variability, but also narrows down the importance of the choice of the window size, which needed to be manually specified in the previous approaches.

Moreover, as hypothesized, the power of the multi-scale dictionary becomes evident when applied to the reconstruction task, e.g. Fig. 2. The dictionary is able to capture the major components of variation of the time series. In order to generate the new canonical signal, a B-spline is fitted to the concatenated best matched dictionary elements. Fig. 2. Reconstruction of a time series representing a coordinate of the relative position of one of the limbs using the best matching dictionary elements. The original signal is shown in blue, and the approximation by the dictionary is shown in red. Fig. 2(a) shows a case where the quantization of the time scale, shape, and mean contains the correct elements needed for reconstruction. The example in Fig. 2(b) shows a case where several very rough approximations are made due to the coarseness of our quantization (mean K = 5, time scale K = 5 and shape at each scale K = 4).

Classification

We show that our dictionary of action primitives representation naturally lends itself to the action classification task, enabling a range of inference approaches. In particular, we demonstrate two inference techniques: A bag-of-words approach, which ignores temporal constraints, and a string classification technique that explicitly takes into account the sequence in which action primitives occur.

Characterizing an Action

Once dictionaries are constructed for all dimensions, any action can be expressed by projecting each component time series onto the appropriate dictionary. The projection procedure follows a similar pipeline to dictionary construction. First, sub-windows are extracted in the same manner as during dictionary construction. Each sub-window is then matched to the closest quantized scale, mean, and trajectory in the given dictionary and assigned unique labels indicating the correspondences. For the simplest dictionaries only the trajectories are matched, and for the mean-quantized dictionaries only mean and trajectory are considered. In the case of multi-scale dictionaries the trajectory is resampled to the closest canonical size before comparing mean and shape. The result of projecting a time series onto a dictionary is thus a sequence of labels annotating all interest points.

Bagging Action Primitives

Our first classification approach discards all causal constraints. It is thus meant to serve as a baseline for more complex techniques. In this approach an action is characterized by comparing the distribution of dictionary elements within the whole action. This is known as the bag-of-words approach. It has shown to be very effective in the domains of document classification, category recognition in images, as well in behavior analysis [START_REF] Dollar | Behavior recognition via sparse spatio-temporal features. Visual Surveillance and Performance Evaluation of Tracking and Surveillance[END_REF][START_REF] Ikizler | Searching video for complex activities with finite state models[END_REF]. To classify actions with this technique, we build a histogram of labels for each dimension of the time series. The dimensionality of the histogram for each time series is equal to the number of elements in the dictionary used. Each histogram bin records the number of times a particular dictionary element occurred during the action. For mean quantized dictionaries, we consider distributions of mean and shape in independent histograms. We found empirically that this independence assumption shows almost no difference from results given when performing inference with multi-dimensional histograms to couple the mean and trajectory labels. Each histogram is turned into a distribution independent of action duration via normalization. Finally, all histograms from each dimension are stacked to create a single action descriptor.

Once the stack of histograms characterizing each action is computed, we use a standard supervised SVM classification approach [START_REF] Vapnik | The nature of statistical learning theory[END_REF]. Given a partition of the action samples into labeled testing and training data, we train a multi-class χ 2 kernel SVM on the histograms. Any action or action segment can then be classified with the SVM, once it is converted to the histogram of action primitives representation.

Exploiting Temporal Constraints

To incorporate temporal constraints into our inference, we map the problem of classifying time series to a problem of classifying strings [START_REF] Teo | Fast and space efficient string kernels using suffix arrays[END_REF][START_REF] Yang | Learning actions using robust string kernels[END_REF]. By projecting the interest points of a time series to our dictionary, as described, we obtain a sequence of labels which can be viewed as a "string" representation of the time series. The alphabet elements are the elements of our dictionary. Our approach is motivated by algorithms from the domain of protein sequence similarity detection [START_REF] Liao | Combining pairwise sequence similarity and support vector machines for remote protein homology detection[END_REF]. Our classification procedure can be described in three steps:

First, we compute the pairwise similarity between all the string sequences that describe each sample of an action. There are several techniques for string matching in the literature [START_REF] Smith | Identification of common molecular subsequences[END_REF][START_REF] Needleman | An efficient method applicable to the search for similarities in the amino acid sequences of two proteins[END_REF], that are commonly used to align proteins or nucleotide sequences. In our implementation, we used the Smith-Waterman (SW) algorithm; a dynamic programming algorithm which finds the local alignment of two string sequences while allowing gaps. In order for the similarity score to be independent of the length of the two sequences that we compare, we normalize the score by the length of the aligned sequences. Similarities computed for each dimension of time series of an action are equally weighted so as to obtain an overall similarity score for each action sample. We also experimented with jointly aligning all channels with respect to a single temporal warping, as opposed to independently aligning each channel. While this approach is physically more plausible (there is, after all, one temporal dimension), the resulting model achieves classification results that are slightly worse than those using independent alignment on each channel. This empirical finding remains to be fully understood.

The second step is to represent each action via its similarity to the actions in the training set. More specifically, the action F is represented by a vector of scores:

X F = [x f1 , x f2 , . . . , x fn ] (1)
where f i is the ith training set sequence, and x fi is the value of the SW score between sequence F and f i . This procedure can be viewed as the vectorization step of a kernel method approach, such as a SVM. The third step is the definition of the positive definite kernel function:

K(X, Y ) = X • Y (X • X)(Y • Y ) . (2)
This kernel K(•, •) is then transformed into a radial basis kernel using the induced distance:

D(X, Y ) = K(X, X) + K(Y, Y ) -2K(X, Y ) (3)
Once we have the radial basis kernel, we train a multi-class SVM, and then classify actions.

Experiments

We applied both of our inference approaches to human motion classification in both motion capture and video data. Both fixed and adaptive window size dictionaries were used. We show that when using motion capture data, the ideal case for our approach, we achieve state of the art performance, demonstrating 5 to 8% improvement in classification over recent dynamics-based approaches. We also show that our methods perform competitively in video, despite the availability of only a rough estimate of relative limb locations.

Datasets

As presented in [START_REF] Ali | Chaotic invariants for human action recognition[END_REF], the FutureLight motion capture dataset [34] contains a total of 158 action samples which represent 5 actions, with a variable number of samples of each. The actions include dance, jump, run, sit and walk. Each action contains considerable intra-class variability, including several fairly ambiguous samples, such as a ballet dance that appears visually similar to a walk. As mentioned earlier, we choose to represent each recorded action with a 15 dimensional time series.

To demonstrate the applicability of our approach to video-based action classification, we evaluate performance on the Weizmann Human Action dataset provided by [START_REF] Blank | Actions as Space-Time Shapes[END_REF]. This dataset contains video of 9 different people, performing 9 unique actions, including: run, walk, sideways run, jump, jump in place, jumping jack, one-handed wave, two-handed wave, and bend. The dataset is annotated with silhouettes for each sample action, obtained using a background subtraction algorithm. Based on the silhouette properties we represent each action as a 5 dimensional time series. Specifically we split the silhouette into 4 quadrants at its center of mass and track the max width of the silhouette in each quadrant. We also track the absolute height of the silhouette over time. These features provide a very rough estimate of the positions of the body extrema over time and it is computationally more efficient than the spatio-temporal "cubes" that [START_REF] Blank | Actions as Space-Time Shapes[END_REF] used or the orientations of rectangular patches of the silhouette used by [START_REF] Ikizler | Searching video for complex activities with finite state models[END_REF].

Results

Performance was evaluated on the datasets above using a leave-one-out cross validation procedure. Below we report the best results achieved for each combination of dictionary type and inference technique after exploring a range of various parameters: Surprisingly, when ignoring global temporal information and using the bagof-words classification approach, we were able to achieve excellent classification performance in both motion capture and video data. Using dictionaries with fixed window size (BAG+SVM) we obtained 95.3% mean classification accuracy on the FutureLight dataset (Table 1, window size: 21 samples) and 91.36% mean accuracy on the Weizmann action dataset (Table 5, window size: 17 samples). Using multi-scale dictionaries (BAG+Temporal Scale+SVM) we obtained 97.40% mean classification accuracy on the FutureLight dataset (Table 3), and 95.06% mean classification on the Weizmann dataset (Table 7).

Incorporating temporal constraints (SW) with our second inference method, we obtained improved performance in both datasets for fixed window and adaptive window size dictionaries. This is shown in Tables 2 and4 for the Future-Light dataset and Tables 6 and8 for the Weizmann dataset. For the FutureLight dataset performance increased by a little over a half percent, but for classification in the Weizmann dataset, we obtained a significant improvement of roughly 5 to 6%. The small increase in improvement in FutureLight is due to the already excellent performance of the bag-of-words approach. It is also worthy to note that these datasets do not contain actions where the global temporal order of the action is crucial. Confusion Matrices for FutureLight (Tables 1234) and Weizmann (Tables 5678) results.

Conclusion

This work proposed an approach for constructing dictionaries of primitive actions suitable for action classification from time series data. We explored the representative and discriminative power of these dictionaries using them in reconstruction and classification tasks, respectively. We achieved state of the art classification results in the FutureLight motion capture dataset using both bag-of-words and string matching inference tech-niques. Also, we demonstrated competitive results in video data, while using significantly simpler features with only 5 dimensions. Together, these results make a compelling argument for viewing actions from the perspective of multidimensional time series.

1 :Run 1 :Frame

 11 Reconstruction based on the Best Matching Dictionary Elements ( Left Foot Z-Coordinate) Reconstruction based on the best Matching Dictionary Elements (Right Hand Y-Coordinate)

Table 1 .

 1 BAG+SVM

		Dance Jump Run Sit Walk	Dance Jump Run Sit Walk	Dance Jump Run Sit Walk
	Dance 28	1	2	Dance 29	2		Dance 28	1	2
	Jump		13	1	Jump	14		Jump	14
	Run	1	28	1	Run	4	26	Run	29	1
	Sit		35		Sit		35	Sit	35
	Walk			48	Walk		48	Walk		48
					Table 2.
					SW+SVM

Table 3 .

 3 BAG+Temporal Scale+SVM

	Dance Jump Run Sit Walk
	Dance 28	2
	Jump	14
	Run	29	1
	Sit	35
	Walk		48
	Table 4. SW+Temporal
	Scale+SVM

Table 5

 5

		A1 A2 A3 A4 A5 A6 A7 A8 A9
	A1 8 1
	A2	9
	A3	9
	A4	9
	A5	9
	A6	9
	A7	9
	A8	9
	A9	1	8
	.
	BAG+SVM

Table 6

 6

	A1 A2 A3 A4 A5 A6 A7 A8 A9
	A1 9
	A2	9
	A3	9
	A4	8 1
	A5	1 8
	A6	9
	A7	1	7 1
	A8		9
	A9		9
	.
	SW+SVM

Table 7

 7

	A1 A2 A3 A4 A5 A6 A7 A8 A9
	A1 9
	A2	9
	A3	9
	A4	9
	A5	9
	A6	9
	A7	9
	A8	9
	A9	9
	. BAG+Temporal
	Scale+SVM

Table 8 .

 8 SW+Temporal Scale+SVM

A1: Run, A2: Walk, A3: Side, A4: Jump, A5: PJump, A6: Jack, A7: Wave 1, A8: Wave 2, A9: Bend

Acknowledgments. We would like to acknowledge the support of AFOSR FA9550-06-1-0138, ONR 67F-1080868 and NSF ECS-0622245.


Chen Change Loy
email: ccloy@dcs.qmul.ac.uk

Tao Xiang
email: txiang@dcs.qmul.ac.uk

Shaogang Gong

From Local Temporal Correlation to Global Anomaly Detection

In this paper, we propose a novel framework tailored towards global video behaviour anomaly detection in complex outdoor scenes involving multiple temporal processes caused by correlated behaviours of multiple objects. Specifically, given a complex wide-area scene that has been segmented automatically into semantic regions where behaviour patterns are represented as discrete local atomic events, we formulate a novel cascade of Hidden Markov Models to model behaviours with complex temporal correlations by utilising combinatory evidences collected from local atomic events. Using a cascade configuration not only allows for accurate detection of video behaviour anomalies, more importantly, it also improves the robustness of the model in dealing with the inevitable presence of errors and noise in the behaviour representation resulting less false alarms. We evaluate the effectiveness of the proposed framework on a real world traffic scene. The results demonstrate that the framework is able to detect not only anomalies that are visually obvious, but also those that are ambiguous or supported only by very weak visual evidence, e.g. those that can be easily missed by a human observer.

Introduction

Automated anomaly detection has drawn increasing attention in computer vision research [START_REF] Hu | A survey on visual surveillance of object motion and behaviors[END_REF][START_REF] Dee | How close are we to solving the problem of automated visual surveillance[END_REF][START_REF] Moeslund | A survey of advances in vision-based human motion capture and analysis[END_REF] due to its potential applications in security and visual surveillance, video highlight extraction and video summarisation. These applications have traditionally relied on human operators to monitor the videos continuously and determine whether behaviours being monitored are normal or abnormal. Anomaly detection systems are anticipated to ease the burden of human operators by automatically sending alerts on abnormal behaviours that should be given further examination.

Abnormal behaviours may involve either single object or multiple objects. Anomalies in the former case are mostly caused by a large spatial deviation and/or an obvious change of visual appearance or motion characteristic triggered by objects in isolation, e.g. a person running whilst others walking, or a car doing U-turn on a one way street. They are visually more distinctive and can be easily detected either by human operators or most of the existing techniques such as template/trajectory matching [START_REF] Bobick | The recognition of human movement using temporal templates[END_REF][START_REF] Johnson | Learning the distribution of object trajectories for event recognition[END_REF]. Detection of abnormal behaviours involving multiple objects is more challenging, especially when the anomalies are caused by abnormal correlations among behaviours of multiple objects, e.g. a vehicle changes lane at a wrong time forcing other vehicles to slow down. These anomalies are more subtle to detect because they often take place in a larger spatial and temporal context. Conventional methods that consider only behaviour of individual object in isolation are inadequate in detecting this type of anomalies [START_REF] Brand | Discovery and segmentation of activities in video[END_REF][START_REF] Nguyen | Recognising behaviours of multiple people with hierarchical probabilistic model and statistical data association[END_REF] because each individual object's behaviour can be perfectly normal without being considered in a global context, e.g. in the aforementioned example, neither a vehicle changing lane or a vehicle slowing down is abnormal when viewed in isolation.

Moreover, the characteristics of different types of abnormal behaviour involving multiple objects are visually more difficult to define and quantify. Until now, there is little if any work done to address this problem systematically. To that end, we consider that multi-object anomalies can be grouped into three categories based on their visual distinctiveness and their frequency of occurrence in the training set. The first category (Category-A) of abnormal correlations often causes an apparent visual changes which are very different from what are observed from the training set (e.g. an accident in a traffic junction causing interruptions to the overall traffic flow). Category-B corresponds to anomalies that are visually ambiguous due to their rare occurrence in the training set. The third category (Category-C) of anomalies are supported only by very weak visual evidence, i.e., featured with very subtle deviation from the normal temporal order/durations of different correlated temporal processes (e.g. a train arrives slightly ahead of schedule on a platform). Anomalies in both Category-A and Category-C refer to those that have never occurred in the training set, whilst anomalies in Category-B are those that appear in the training set but are statistically under-represented. From the perspective of a human observer, anomalies in Category-A are visually obvious thus easy to detect. In contrast, both anomalies in Category-B and Category-C are more likely to be missed due to the lack of visual evidence.

A few approaches have been proposed for modelling activities involving multiple objects [START_REF] Brand | Discovery and segmentation of activities in video[END_REF][START_REF] Nguyen | Recognising behaviours of multiple people with hierarchical probabilistic model and statistical data association[END_REF][START_REF] Oliver | A bayesian computer vision system for modeling human interactions[END_REF][START_REF] Gong | Recognition of group activities using dynamic probabilistic networks[END_REF], little effort has been taken in detecting abnormal interactions that are either ambiguous or weakly supported by visual evidence. Previous work mainly focused on classifying interactive behaviours [START_REF] Nguyen | Recognising behaviours of multiple people with hierarchical probabilistic model and statistical data association[END_REF][START_REF] Oliver | A bayesian computer vision system for modeling human interactions[END_REF][START_REF] Gong | Recognition of group activities using dynamic probabilistic networks[END_REF] or detecting anomalies with apparent visual evidence (Category-A) [START_REF] Brand | Discovery and segmentation of activities in video[END_REF]. In contrast, we focus on the detection of abnormal interactions involving multiple objects (we refer to as global anomalies), particularly those belonging to Category-B and Category-C, for which modelling the correlations of multiple temporal processes in an effective and efficient manner becomes crucial. In particular, considering that in a busy public space different objects' behaviours are correlated either implicitly or explicitly with each other, we propose to exploit the combinatory evidence obtained from each individual temporal process to support the decision in global anomaly detection.

Related Work

A number of approaches have been proposed to model behaviours of multiple objects without considering the temporal correlations among multiple objects [START_REF] Brand | Discovery and segmentation of activities in video[END_REF][START_REF] Nguyen | Recognising behaviours of multiple people with hierarchical probabilistic model and statistical data association[END_REF]. The assumption of different temporal processes being independent is clearly invalid when dealing with behaviours involving multiple interacting objects. One of the early work in modelling correlations among objects was undertaken by Oliver et al. [START_REF] Oliver | A bayesian computer vision system for modeling human interactions[END_REF][START_REF] Oliver | Statistical modeling of human interactions[END_REF], in which a Coupled Hidden Markov Model (CHMM) [START_REF] Brand | Coupled Hidden Markov Models for complex action recognition[END_REF] is used to model people interactions based on their proximity and trajectories in an outdoor scene. Causal relationships among multiple temporal streams are taken into account by fully coupling pairs of Hidden Markov Models (HMMs) so that each hidden state is conditionally dependent on all past states of all processes in the previous time instance. An alternative solution based on Dynamically Multi-Linked HMM (DML-HMM) is introduced by Gong and Xiang [START_REF] Gong | Recognition of group activities using dynamic probabilistic networks[END_REF]. In contrast to CHMM, causal temporal relationships are established only among a subset of hidden states with high relevance across multiple temporal processes. Both CHMM and DML-HMM have been shown to be effective in modelling correlated behaviours. However, the behaviours being modelled in [START_REF] Oliver | A bayesian computer vision system for modeling human interactions[END_REF][START_REF] Gong | Recognition of group activities using dynamic probabilistic networks[END_REF] are limited to small local regions. Temporal correlations in a larger spatial context are not investigated. In addition, there are no attempt to detect ambiguous anomalies and subtle anomalies supported with weak evidence. In contrast, our cascade HMMs are designed to learn efficiently from combinatory evidence extracted from local temporal processes, thus providing a mechanism to detect abnormal correlations in a complex scene with a large spatial context. Moreover, thanks to the cascade model structure, our model is less susceptible to noise and errors in behaviour representation compared with an alternative single stage model.

More recently, a static model based on a hierarchical probabilistic latent semantic analysis (pLSA) is introduced by Jian et al. [START_REF] Li | Global behaviour inference using probabilistic latent semantic analysis[END_REF] for global phase inference and anomaly detection in traffic scenes. A scene segmentation algorithm [START_REF] Li | Scene segmentation for behaviour correlation[END_REF] is used to decompose a scene into semantic regions, where each of them encloses an area that contains a set of correlated atomic events. The bottom layer of the hierarchical pLSA is employed to model the local atomic events, whilst the top layer is used to model the global correlations among them. Despite the proposed method has been shown to be capable in detecting abnormal behaviours in a global context, it is limited to modelling static causal relationships without taking the temporal ordering and duration aspects of behaviours into account. The model is thus unable to detect anomalies embedded in the temporal structure of correlation.

In this work, we propose a novel framework based on a cascade of HMMs to model the multiple temporal processes in busy outdoor scenes. The use of cascade structure offers a number of unique advantages over the existing Dynamic Bayesian Network (DBN) based models: (1) it provides a mechanism that exploits evidences from local temporal processes to detect video behaviour anomalies supported only by weak visual evidence; (2) it provides effective temporal correlation modelling, enabling it to detect anomalies embedded in the temporal structure of correlation; (3) the cascade structure is less susceptible to noise compared with a single model that learns directly from the observational space. Our approach overcomes the problems associated with a static model such as Jian et al. [START_REF] Li | Global behaviour inference using probabilistic latent semantic analysis[END_REF], as this is a dynamic temporal model. Experiments are carried out on a busy traffic scene. The results show that the framework is able to detect not only anomalies that are visually obvious, but also those that are ambiguous or supported only by very weak visual evidence, i.e., those that can be easily missed even by a human observer.

Methodology

Behaviour Representation

Since reliable tracking is difficult and discontinuity in object trajectories is prevalent in a crowded outdoor scene, trajectory-based representation is not employed in this study.

Instead, we adopted a discrete event based approach introduced by Jian et al. [START_REF] Li | Scene segmentation for behaviour correlation[END_REF] to represent global behaviours in a wide-area outdoor scene. We provide a concise description of the approach here to facilitate explanations in other sections appeared later. Interested readers may refer to the original paper [START_REF] Li | Scene segmentation for behaviour correlation[END_REF] for more details.

A continuous video sequence V is first segmented uniformly into T non-overlapping video clips V = {v 1 , . . . , v t , . . . , v T }, with each video clip v t containing N f frames. Foreground blobs with similar spatio-temporal features are then clustered using kmeans into a set of atomic events (e.g. vehicles stop at the middle of intersection waiting for right-turning) across all the frames in a clip v t . The spatio-temporal features employed in the clustering include the centroid position and the width and height of the object bounding box, the ratio of width and height, the occupancy ratio, the optical flow vectors for the bounding box and the scaled optical flow vectors according to the bounding box size. Upon obtaining atomic events for all clips, global clustering based on Gaussian Mixture Model (GMM) is performed to group the atomic events into coarse global event classes. Based on the distributions of the clusters, spatial scene segmentation is carried out to decompose a scene S into R semantic regions, where S = {s 1 , . . . , s r , . . . , s R }. Specifically, two pixels are considered to be similar and grouped together if similar classes of events occurred there. Consequently, object behaviours within each segmented region are similar to each other whilst being different from those in other regions. To detect the atomic events more accurately, the aforementioned clustering method is repeated within each region with automatic feature selection to group foreground blobs into finer regional event classes. As a result, a video clip v t is spatially represented by segmented regions (see Fig. 1(b)), where each region contains a set of correlated regional atomic events. To construct the input features for the proposed cascade model, we represent the behaviours captured in a video clip v t using R binary vectors y 1 t , . . . , y r t , . . . , y R t , which correspond to the occurrence of regional events in each region. Specifically, the binary vector y r t is given as

y r t = [O 1 , . . . , O m , . . . , O Mr ],
where

O m = 1 if atomic event e r m occurs, 1 ≤ m ≤ M r 0 otherwise . ( 1
)
where e r m is the atomic event belonging to the mth regional event class in region s r and M r is the total number of regional event classes in region s r .

Cascade Model Structure

The proposed framework combines two stages of Multi-Observation HMMs (MOHMMs) [START_REF] Gong | Recognition of group activities using dynamic probabilistic networks[END_REF] in a cascade, and the model is represented as a DBN [START_REF] Murphy | Dynamic Bayesian Networks: Representation, Inference and Learning[END_REF]. The first stage of the model is composed of multiple MOHMMs, each of which is used to model the temporal evolvement of regional atomic events within each region. The second stage consists of a single MOHMM for modelling the state sequence inferred from the first stage, and is responsible for learning the global correlations among temporal processes across regions. In contrast to conventional HMM that emits a deterministic symbol in a given state, a MOHMM allows a state to produce multiple symbols in every time step. It is thus ideal for modelling multiple atomic events temporally in each region in the firststage cascade, and modelling the temporal correlations of behaviours across multiple regions in the second-stage cascade.

The MOHMMs in both stages are ergodic models and characterised by a tuple λ = (Q, O, π, A, B) where Q is the state variable, O is the set of observation symbols, π is the initial state distribution, A the state transition distribution (a ij = P (Q t+1 = q j |Q t = q i )), and

B is observational distribution (b j (m) = P (O t = o m |Q t = q j )).
In the following text, the MOHMMs at the first stage are known as local models and denoted as Λ l , where Λ l = {λ 1 , . . . , λ r , . . . , λ R }, and R is the number of regions. The second stage MOHMM is known as global model denoted as Λ g . We denote the state space of hidden variable Q l r of a local model as q l r i where 1 ≤ r ≤ R, and the state space of hidden variable Q g of a global model is denoted as q g i , where i is the state index. Figure 2 the number of clusters k equals to the number of hidden states N s . However, k-means clustering itself converges to numerous local minima depending on the initial cluster centroids. To overcome this problem, we repeat the clustering procedure for 100 times with the initial cluster centroids being selected randomly. We then select the cluster formation that minimises the sum of squared error (SSE) distances between the data and the centroid of their clusters. Specifically, the SSE is defined as:

SSE = k i=1 y∈ci (y i -c i ) 2 . ( 2
)
where c i is the centroid of cluster c i , where

1 ≤ i ≤ k.
Subsequently, we initialise the model λ r based on the result obtained from k-means clustering. Specifically, the observation probability distribution B = {b i (m)} is estimated by calculating the frequency of occurrence of the symbol in the corresponding cluster, which can be written as:

b i (m) = P (O m at t|Q l r t = q l r i ) = # O m in c i T . (3)
where 1 ≤ i ≤ N s and 1 ≤ m ≤ M r . We use uniform estimates for the initial state distribution π = 1/N s and state transition distribution A. Upon determining the initial values of the parameters, learning proceeds to estimate the parameters through EM algorithm.

After the training of local models Λ l , we proceed to the training of the global model Λ g . First, the most probable explanation (MPE) through a local model λ r is obtained by using Viterbi decoding, given as:

Q l r 1:T = argmax Q l r 1:T P (Q l r 1:T |y r 1:T ). (4)
The MPE obtained across all first-stage models represents the understanding of the learned models regarding the phases of individual temporal processes. With the help of Viterbi algorithm, the temporal and causal relationships learned during the model training can be employed to explain away the errors in the observational space. This would reduce the influence of noise during second-stage model learning. The MPE obtained across all local models are then augmented, forming an intermediate observational vector, given as:

z t = [Q l 1 t , . . . , Q l r t , . . . , Q l R t ]. (5)
The input is then transmitted to the global model Λ g for training. Parameters initialisation and estimation of the global model are carried out with the similar steps applied to the local models.

Anomaly Detection

On-line abnormality detection can be achieved by examining the log-likelihood of a given sequence with respect to the model over time. A low log-likelihood value may signify an occurrence of anomaly. Given an unseen sequence y r 1:t , where 1 ≤ t ≤ T , the normalised log-likelihood LL l r t at clip t with respect to the local model λ r and the normalised log-likelihood LL g t at clip t with respect to the global model Λ g can be computed as follows:

LL l r t = 1 t log P (y r 1:t |λ r ). ( 6
)
LL g t = 1 t log P (z 1:t |Λ g ). (7)
Both log-likelihood values can be used for abnormality detection. However, since we are interested in global behaviour anomalies that are defined in the correlations of multiple temporal processes, the use of global log-likelihood LL g t for anomaly detection would be more appropriate.

To obtain LL g t in an on-line manner so that anomalies can be detected on the fly, we need to estimate the local hidden state Q l r t instantly at every clip t. By on-line estimation we mean the past evidence before the current clip t need not be presented to the model for the estimation of most likely state. This results in a constant computational time in computing LL g t over time. To estimate Q l r t , the marginal probabilities P (Q l r t = q l r i |y r 1:t ) are first computed, where q l r i represents the ith hidden state and y r 1:t is the unseen sequence between times 1 and t. The most likely hidden state is then determined by choosing the hidden state that yields the highest marginal probability:

Q l r t = argmax q l r i P (Q l r t = q l r i |y r 1:t ). (8)
The computation of marginal probability takes constant space and time at every time instance by using recursive estimation to compute P (Q l r t |y r 1:t ) as a function of current input y r t and prior belief state P (Q l r t-1 |y r 1:t-1 ). Specifically,

P (Q l r t |y r 1:t ) ∝ P (y r t |Q l r t , y r 1:t-1 )P (Q l r t |y r 1:t-1 ) ∝ P (y r t |Q l r t )[ Q l r t-1 P (Q l r t |Q l r t-1 )P (Q l r t-1 |y r 1:t-1 )]. (9)
where the constant of proportionality is 1/P (y r t |y r 1:t-1 ). The procedure described here is traditionally known as "filtering", because we are filtering out the noise from the observations [START_REF] Murphy | Dynamic Bayesian Networks: Representation, Inference and Learning[END_REF].

The most likely hidden state sequences obtained from each local model at clip t are augmented according to [START_REF] Johnson | Learning the distribution of object trajectories for event recognition[END_REF] to form an observational input for global model. Subsequently, we can compute the LL g t following [START_REF] Nguyen | Recognising behaviours of multiple people with hierarchical probabilistic model and statistical data association[END_REF]. Global anomaly detection can be achieved by comparing LL g t with a pre-defined threshold T H. Specifically, if LL g t < T H, the unseen clip is detected as abnormal.

Experiments

Dataset

Experiments were conducted on a busy road traffic junction scene. The video footage was recorded at 25Hz and resized to a resolution of 360 × 288 pixels. The total duration of the recording is approximately 22 minutes, showing a busy road junction regulated by traffic lights, dominated by four types of traffic flows as illustrated in Fig. 3. Specifically, Flow A corresponds to traffic in vertical directions. Flow B, C and D are regarded as traffic flows in horizontal directions. In particular, Flow B represents left-turning and right-turning traffics by vehicles from vertical directions. Flow C corresponds to rightward traffic and flow D corresponds to leftward traffic. The order of the traffic flow depends upon how busy the vertical traffics are. During most of the recording, the scene is extremely crowded. Consequently, Flow B can only take place after Flow A finishes and is followed by C and D, (i.e. the typical temporal order is A, B, C, D). However, it is noted that very occasionally, there is gap in Flow A which is big enough for Flow B to take place until the gap closes. In other words, Flow A and B can occur alternatively during the vertical traffic phase. This makes global behaviour modelling and anomaly detection challenging in this scene as vehicles do not interact simply following the traffic light cycles, but also driven by the traffic volume as well as the driving habits and reactions of drivers. A total of 112 non-overlapping clips were segmented from the video. In particular, 73 clips (22000 frames) were used for training, whilst 39 clips (12000 frames) were reserved for testing. Scene segmentation was applied on the dataset and a total of 30 local atomic events were discovered in six regions shown in Fig. 1(b).

Prior to the evaluation of the proposed method, ground truth was first obtained by manually performing exhaustive frame-wise examination on both the training and testing set. Consequently, 5 out of 39 testing clips were found to be abnormal, they are clips 3, 4, 10, 25 and 35. Among them, clips 3 and 4 were categorised in Category-A since they contain abnormal behaviours that were visually recognisable. In particular, the abnormality was caused by a fire engine approaching from the left entrance towards the right exit, causing interruptions to the vertical traffic at both directions (see Fig. 4(a Clip 10 was labelled as abnormal behaviours belonging to Category-C, where the abnormality is almost undetectable by human without comprehensive examination of the traffic cycle duration over time. Clip 10 (see Fig. 4(e)) shows a white van on its way from the left entrance to right exit. After careful examination on the traffic cycle duration and the distance between the van and the vehicles behind, it is discovered that the van ran the red light before traffic flow B finished its cycle. The ground truth is summarised in Table 1.

Category Description

Clip A Abnormal behaviours that are visually obvious 3, 4 B Rare and ambiguous behaviours 25, 35 C Abnormal behaviours supported by weak evidence 10 Table 1. Ground truth was defined in accordance to the anomaly's visual distinctiveness and its frequency of occurrence in the training set.

Anomaly Detection

We employed a single MOHMM as a baseline method in this study. The number of hidden states were set to 2 with each of them represent the horizontal and vertical traffic flows, respectively. Unlike the proposed cascade model, the single MOHMM learns directly from the observational space and ignores the regional segmentation. Thus, each hidden state of the model consists of 30 observation nodes, where each node corresponds to an atomic event detected. Similar to the setting of cascade model, we used uniform initialisation for the prior distribution π and state transition distribution A for the single MOHMM. The observational distribution B was initialised according to [START_REF] Moeslund | A survey of advances in vision-based human motion capture and analysis[END_REF].

For the cascade model, the number of hidden states in global model was set to 2 since the global traffic flows have two phases. To obtain the optimal number of hidden states in each local model, we varied the number of local states from 2 to 10, and observed the matching accuracies of the global traffic phases inferred using the training data in each setting with the ground truth global phases. It was found that local model with five hidden states yielded the best accuracies. Figure 5(a) shows the plot of the matching accuracies averaged over 10 runs along with the standard deviations, against the number of hidden states. We trained the single MOHMM and the cascade model using the training set with 20 iterations of EM learning. Both models were then evaluated using the test set to obtain the log-likelihood values (LL g in the case of cascade model). The threshold variable T H were varied to obtain the Receiver Operating Characteristic (ROC) curves of the two models, as shown in Fig. 6(a). As can be seen from the ROC curves, the detection result of the proposed cascade model is significantly better than those obtained using the single MOHMM. In particular, the area under ROC (AUROC) achieved by using the single MOHMM was 0.5794, compared with 0.9765 obtained by using cascade MOHMMs.

To further investigate the miss detections and false alarms, the global log-likelihoods (LL g in the case of cascade model) obtained from the both models are plotted as shown in Fig. 6(b) and 6(c), respectively. From the log-likelihood plots, we observed a high false alarm rate in anomaly detection with the single MOHMM. In particular, video clips 13, 15, 14, 18, 3, and 1 are those clips that returned the lowest log-likelihoods, despite the fact that most of them (except clip 3) are normal. The high false alarm rate suggests that the model is susceptible to noise since it learns directly from the observational space. For instance, in clip 13, vehicles in regions 1 and 4 were mistakenly grouped as a large blob across regions (see Fig. 7). The same error was also observed in clip 14. As a result of imperfect blob detection, the atomic events were grouped into the wrong clusters and consequently led to a false alarm using the single MOHMM. In contrast, clips 13 and 14 yielded much higher log-likelihood values using the cascade model, indicating the model is able to cope with the erroneous input effectively. Apart from the noise problem, the large observational space per state of the single model implies that a single state has to represent implicitly multiple sources of atomic event variations at any given time instance. Unless we have sufficient amount of clean training instances, the poor structure would leads to ineffective learning of temporal correlation among atomic events. In contrast, as can be seen from Fig. 6(c), the log-likelihood plot of the cascade model is more selective and the model correctly picked up clips 3, 4, 10, 25 and 35 as abnormal behaviours. As expected, the cascade model to identified clips 3 and 4 as abnormal behaviour since these behaviours are supported with strong visual evidences. In addition, the proposed model was able to detect clip 25 and 35 as being abnormal, despite the fact that they are rare, ambiguous and harder to be identified without careful examination on the traffic cycles. More importantly, the proposed model successfully discover clip 10 as being abnormal, despite the behaviour is visually much less obvious. This demonstrates the advantage of our cascade model, i.e., although the evidence observed from a single atomic event is inadequate, the model is still able to detect global anomalies by exploiting combinatory temporal evidences from individual atomic events.

Conclusions

We have presented a novel framework for anomaly detection in outdoor crowded scene by modelling temporal correlations of multiple objects. Experiments on complex traffic scene have demonstrated the promising potential of the proposed model in detecting abnormal behaviours that are ambiguous or supported only by weak visual evidence. This is achieved by modelling temporal correlations among multiple atomic events and integrating lower-level evidences to increase the certainty in detecting anomalies. Our experimental results have also shown that cascade model is more accurate in detecting abnormal behaviours compared with single MOHMM model. The better performance is mainly due to the factorisation of observational space via semantic scene segmentation and the robustness to noise and errors in behaviour representation.

We should point out that different types of HMMs can be employed at both stages depending on detection requirement and scene complexity. The choice of models could range from a simple first-order HMM to more sophisticated models such as Hierarchical Hidden Markov Model (HHMM) or Hidden semi-Markov Model (HSMM). Our ongoing work is focused on investigating the incorporation of duration modelling capability into the framework to achieve a more accurate modelling of behaviour correlations with long-term temporal dependency.

  (a) A traffic scene (b) Six segmented regions

Fig. 1 .

 1 Fig. 1. Semantic scene segmentation.

  illustrates the structure of the proposed framework. Observation nodes are shown as shaded squares and hidden nodes are shown as clear squares.

Fig. 2 .

 2 Fig. 2. Dynamic Bayesian Network (DBN) representation of the proposed cascade model with one time slice unrolled.

3. 3

 3 Model Training Model training is accomplished through two steps, i.e., the local models training and global model training. We use training sequence y r 1:T = {y r 1 , . . . , y r t , . . . , y r T } to estimate parameters of the model λ r through Expectation Maximisation (EM) algorithm. However, it is well known that EM algorithm may converge to a poor local optimum solution if the parameters are not initialised properly. One possible solution for this is to perform k-means clustering on the training data for model parameter initialisation, with

Fig. 3 .

 3 Fig. 3. Traffic flows observed in the dataset.

  ) and 4(b)). The vehicles from vertical directions gave way to the fire engine which has priority, and the normal traffic flow resumed after the fire engine exited. Clips 25 and 35 correspond to abnormal traffic flow where vehicles did not follow the typical temporal order of A-B-C-D. In particular, clip 25 shows a motorbike making a turn to the left exit from the upper direction during vertical traffic flow A. In clip 35, vehicles were using a gap in the middle of traffic flow A to make right turn and left turn at the same time interval. Both clips were grouped into Category-B because they belong to rare/unusual behaviour with low frequency of occurrence in the training set (out of 73 training clips, only 3 clips correspond to left-turn and 2 clips correspond to turning both ways).

Fig. 4 .

 4 Fig. 4. Clips 3 and 4 are abnormal behaviours that are visually detectable (Category-A), whilst clips 25 and 35 are abnormal behaviours that are ambiguous (Category-B). Clip 10 contains subtle anomaly that is almost undetectable by human (Category-C). The corresponding objects that caused the anomalies are manually highlighted.

Figure 5 (

 5 b) shows an example of the global phases inferred from the training data and the ground truth. As can be seen, the cascade model managed to learn the implicit temporal transition of the vertical and horizontal traffic phases accurately, achieving accuracy rate of 97.26% when the number of local states was set to 5.

Fig. 5 .

 5 Fig. 5. (a) The matching accuracies of global traffic phases inferred from training data with the ground truth global phases averaged over ten runs, along with the standard deviations being shown as error bars. (b) Global phases inferred from training data and the ground truth global phases.

Fig. 6 .

 6 Fig. 6. (a) ROC curves for single MOHMM and the proposed cascade model. Subfigures (b) and (c) show the log-likelihood plots for single MOHMM and cascade MOHMMs, respectively.

Fig. 7 .

 7 Fig. 7. Imperfect blob detections result in local atomic events being grouped into wrong clusters. The corresponding bounding box is marked in black colour.