Proceedings

2014 Workshop on Optimization

Methods for Anomaly Detection

OMAD 2014

Held as part of the 2014 SIAM International Conference on Data Mining

in Philadelphia, PA
24-26 April 2014

Edited By

Sanjay Chawla, Kamalika Das, Aris Gionis

  i  
  

   

  ii  
 Foreward
Anomaly Detection techniques are playing an increasingly important role in the analysis of large
data sets across many application domains. In particular, the use of anomaly detection
techniques for monitoring digital and physical infrastructure is growing rapidly, Other important
application domains include health and climate informatics. The aim of the 2014 workshop on
Optimization Methods for Anomaly Detection (OMAD) is to bring together researchers to study
the detection of anomalies in large data sets in a systematic optimization framework. The
workshop consists of four papers and two keynote addresses.

In ParitoSVR: Parallel Iterated Optimizer for Support Vector Regression in the Primal, the
authors will present a distributed algorithm for support vector regression using the ADMM
framework. The SVM model is then used to identify anomalies as those data points, which have
a large residual value vis-à-vis the model. The authors apply paritoSVR on a real (and large) data
set consisting of fuel consumption patterns in airline flights.

In Anomaly Detection Using Tripoint Arbitration Similarity Method, the author, proposes a
tripoint similarity function to identify outliers. The similarity function is used in a MinMaxCut
optimization framework. The proposed method is validated on an application related to
monitoring computing infrastructure.

A new measure of anomalousness (called q-value) is proposed in Measuring Anomalousness in

Statistical Models. The measure, which is related to p-value, provides a natural way to find
anomalies in clustered data.

In Identifying Precursors to Anomalies Using Inverse Reinforcement Learning, the authors

propose a method for determining pre-cursor signals just before the advent of an anomalous
event. The application domain includes the monitoring of airplane flight data.

The workshop will also host two keynote talks. The first by Professor Vipin Kumar from the
University of Minnesota will highlight the role of anomaly detection in getting a better
understanding of climate data. The second, by Dr. Dragos Margineantu, from Boeing, will focus
on the application of anomaly detection techniques in the airline industry.

The organizing committee would like to thank (i) the authors for participating in the workshop,
(ii) the organizers of the SDM 2014, especially Professor Tina Elliasi-Rad, the SDM workshops
chair, (iii) colleagues who served in the program committee

Finally we would like to thank the sponsors of the workshop: NASA Ames and National ICT
Australia (NICTA) for their generous support.

Thanks!

Sanjay Chawla, Kamalika Das and Aris Gionis

  iii  
 Workshop Committee

Program Chairs:

Sanjay Chawla, University of Syndey

Kamalika Das, UARC, NASA Ames Research Center
Aris Gionis, Aalto University

Program Committee:

Leman Akoglu SUNY, Stonybrook

Arindam Banerjee, University of Minnesota
Kanishka Bhaduri, Netflix Inc.
Varun Chandola, SUNY, Buffalo
Tina Eliassi-Rad , Rutgers University
Jing Gao, SUNY, Buffalo
Manish Gupta, Microsoft India
Aditya Menon, NICTA
Emmanuel Müller , Karlsruhe Institute of Technology
Khoa Nguyen, HCMUT
Nikunj Oza, NASA Ames Research Center
Aditya Prakash, Virginia Tech
Eric Schubert, Ludwig-Maximilians-Universität München
Nikolaj Tatti, Aalto University
Matthijs van Leeuwen, KU Leuven
Hamed Valizadegan , UARC, NASA Ames
Jilles Vreeken, University of Antwerp
Arthur Zimek, Ludwig-Maximilians-Universität München

  iv  
2014 Workshop on Optimization Methods for
Anomaly Detection

OMAD 2012
Table of Contents
Foreword..………………………………………………………………………………..……...iii

Workshop Committee…………………………………………………………………………..iv

Invited Talk Abstracts

Understanding Global Change: Opportunities and Challenges for Data Driven Research……….1
Vipin Kumar

Data Mining Research Questions for Maintenance Tasks………………………………………...2

Dragos Margineantu

Accepted Abstracts

1. ParitoSVR: Parallel Iterated Optimizer for Support Vector Regression in the Primal.……....3

2. Anomaly Detection Using Tripoint Arbitration Similarity Method……………………...…..6

3. Measuring Anomalousness in Statistical Models…………………………………...………..9

4. Identifying Precursors to Anomalies Using Inverse Reinforcement Learning……...………13

  v  
 Understanding Global Change: Opportunities and Challenges for
Data Driven Research
 
Vipin Kumar

The world's population is growing steadily and many countries are simultaneously
industrializing, developments that have been ongoing at varying rates for two
centuries but have accelerated over the past several decades. These processes are
increasingly straining already scarce natural and food resources, which must scale up
to keep pace with growing demand. The consequences of such large-scale changes
include tremendous stresses on the environment that would be calamitous at the
current rate of change if they are not managed sustainably. As a result, scientists are
tasked with providing answers to challenging questions such as: What is the effect of
urbanization on regional land use and ecology? What is the impact of climate change
on global water resources? How does deforestation affect the net carbon balance?
How does increased biofuel production impact crop patterns and food availability?
Addressing these interconnected, societally-relevant questions requires development
of new computational methods that enable monitoring, analysis and understanding of
changes in the Earth system, interactions between different processes, and their
impacts on factors such as the carbon cycle, hydrology, air quality, and biodiversity.

This talk will present an overview of research being done in a large interdisciplinary
project on the development of novel data driven approaches that take advantage of the
wealth of climate and ecosystem data now available from satellite and ground-based
sensors, the observational record for atmospheric, oceanic, and terrestrial processes,
and physics-based climate model simulations. These information-rich datasets offer
huge potential for monitoring, understanding, and predicting the behavior of the
Earth's ecosystem and for advancing the science of global change. This talk will
discuss some of the challenges in analyzing such data sets and our early research
results.

1
 User-in-the-loop Learning and Optimization for Anomalous Action
Detection
Dragos Margineantu

An increasing number of users collect transaction data and need scalable tools that
assist them in identifying abnormalities. This talk will present an interactive user-in-
the-loop approach based on inverse reinforcement learning and linear optimization
methods for detecting anomalies and intent in data. We implemented and tested our
algorithms on real-world GMTI and AIS sensor data.

2
ParitoSVR: Parallel Iterated Optimizer for Support Vector Regression
in the Primal
Kamalika Das∗ Kanishka Bhaduri† Nikunj Oza‡

Abstract on these input parameters. One such popular regression

Regression problems on massive data sets are ubiquitous method is Support vector machines (SVM) [1] which is a
in many application domains including the Internet, earth class of maximum margin classifiers, that demonstrates
and space sciences, and aviation. Support vector regression good generalization performance. SVM’s can also ex-
(SVR) is a popular technique for modeling the input- ploit the kernel trick, thereby making them suitable for
output relations of a set of variables under the added non-linear model learning as well. SVMs however are
constraint of maximizing the margin, thereby leading to computationally expensive for large datasets.
a very generalizable and regularized model. However, for In this paper we propose Parallel Iterated Optimizer
a dataset with m training points, it is challenging to for Support Vector Regression in the Primal (Pari-
build SVR models due to the O(m3 ) cost involved in toSVR), a new support vector regression algorithm that
building them. In this paper we propose ParitoSVR — a can be deployed over a network of machines, where each
parallel iterated optimizer for Support Vector Regression machine solves a small (sub-)problem based only on the
in the primal that can be deployed over a network of data observed locally and these solutions are then com-
machines, where each machine iteratively solves a small bined to form the solution to the global problem. Our
(sub-)problem based only on the data observed locally and proposed method is based on the Alternating Direction
these solutions are then combined to form the solution to the Method of Multipliers (ADMM) optimization technique
global problem. Experiments on real datasets demonstrate [2][3], which is parallelizable for separable convex prob-
the accuracy and scalability of our algorithm. As a real lems, and converges to the exact solution as the central-
application, we use ParitoSVR to detect flights having ized version with theoretical guarantees.
abnormal fuel consumption from a fleet-wide commercial
aviation database. 2 Background
Our ParitoSVR algorithm uses as a building block
1 Introduction two components: (1) Alternating Direction Method of
In many application domains, it is important to predict Multipliers (ADMM), and (2) SVR. In this section, we
the value of one feature based on certain other mea- discuss these two topics.
sured features. For example, in commercial aviation, it ADMM: ADMM [3] is a decomposition algorithm for
solving separable convex optimization problems of the
is very important to model the fuel consumption based
form:
on input parameters such as aircraft speed, wind speed,
control surfaces, engine power, pitch, roll, yaw etc. This (2.1) min G1 (x) + G2 (y)
x,y
is because according to the Air Transportation Associ-
ation (ATA), fuel is an airline’s largest expense at a subject to Ax − y = 0, x ∈ Rn , y ∈ Rm
staggering 17.5 billion gallons per year1 . Identifying where A ∈ Rm×n and G1 and G2 are convex functions.
flights with abnormal fuel consumption may help the ADMM is an iterative technique and the update equa-
airlines to do proper maintenance of these aircrafts and tions are: n 2 o
save operating costs. For such problems, a regression xt+1 = min G1 (x) + ρ/2 Ax − yt + pt 2
x
n 2 o
model can be learned that predicts the fuel flow based yt+1 = min G2 (y) + ρ/2 Axt+1 − y + pt
y 2

pt+1 = pt + Axt+1 − yt+1

∗ UARC, NASA Ames Research Center. kama- where p = (1/ρ)z. ADMM effectively decouples the
[Link]@[Link] x and y updates such that parallel execution becomes
† Netflix Inc. [Link]@[Link]
‡ NASA Ames Research Center. [Link]@[Link] possible. In a distributed computing framework, this
1 [Link] becomes even more interesting since each computing
[Link] node can now solve a (smaller) subproblem in x inde-

3
pendently, and then, these solutions can be efficiently 100
ADMM Iteration 1
ADMM final

gathered to compute the consensus variable y and the 50

y = x β + noise
dual variable p. ADMM converges within a few itera- 0
tions when moderate precision is required. This can be −50
Centralized

particularly useful for many large scale problems, simi- −100

lar to what we consider here. −150
−8 −6 −4 −2 0 2 4 6 8
SVR: Give m data tuples (training set) D = (xi , yi )m
i=1 , x
where xi ∈ R is the input and yi ∈ R is the cor-
n

responding output or target, SVR solves the following Figure 1: Models formed by node 1 on synthetic dataset
optimization problem: as the algorithm progresses.
" m
#
X
2
(2.2) min λ||w|| + `ε (w · xi + b − yi )
w,b
i=1 data and scatter-gather operations on z until they
converge to the same result.
where λ is a constant and `ε is the ε-insensitive loss
Theorem 3.1. The ADMM update rules for the linear
function defined as, `ε (r) = max(|r| − ε, 0). This is a support vector regression primal optimization are:
convex optimization problem which can be solved using Xmj   ρ 2

(j) (j)
convex optimization solvers such as CVX2 . wjt+1 = min `ε wj · xi − yi + wj − zt − utj 2
wj
i=1
2
In the next section we show how to build SVR mod-  

els for very large datasets using distributed computing zt+1 = min λ kzk22 + N ρ z − wt+1 − ut 22
z 2
via the ADMM technique. t+1 t t+1 t+1
uj = uj + wj −z

3 ParitoSVR formulation where u ∈ Rn is the (scaled) dual variable and wt+1

and ut+1 are the averages of the variables over all the
For the linear ParitoSVR algorithm setup, we assume
nodes.
that the training data is distributed among N client
processors (nodes) P1 , . . . , PN with a central machine Proof. We omit the proof here due to shortage of space.
P0 acting as the server or collector. The dataset at
machine Pj , denoted byoDj , consists of mj data points The w update can be executed in parallel for each
n mj machine. It involves solving a convex optimization
(j) (j)
i.e. Dj = xi , yi . It is assumed that the
i=1
T SN problem in n + 1 variables at each node. This solution
datasets are disjoint: Di Dj = ∅ and j=1 Dj = D, depends only on the data available at that partition.
where D is the total (global) data set. The goal is The z update step involves computing the average of
to learn a linear support vector regression model on D the w and u vectors in order to combine the results
without exchanging all of the data among all the nodes. from the different partitions. Critical to the working
Given Eqn. 2.2," the optimization problem # is now: of ADMM is the convergence criteria. The primal and
Xm 2
min `ε (w · xi − yi ) + λ||w|| 2 dual
residuals can be written as: rpt = kwt − zt k2 , rdt =
w
i=1 ρ(zt − zt−1 ) Also, given the thresholds pri and dual ,
 
XN Xm j 
(j) (j)
 the primal√and dual thresholds can be written as,
⇔ min  ` ε w · x i − yi + λ kwk2  = abs m + rel max(kwk , kzk) and and dual =
pri √
w
j=1 i=1
abs m+ρrel kuk . The iterations terminate when rpt <
The inner sum can be computed by each node indepen- t
dently (assuming that w is known). We next write it in pri and rd < dual .
a form such that  it is decoupled across the nodes: 
mj
XN X 
(j) (j)

2
4 Experiments
(3.3) min  `ε wj · xi − yi + λ kzk 
w1 ,...,wN ,z
j=1 i=1
In this section we demonstrate the performance of the
subject to wj = z ParitoSVR algorithm.
In the ADMM decomposition, each node can solve ParitoSVR has been implemented in MATLAB
its local problem using its own data and optimization 2011b. The experiments have been executed in NASA
variable and then coordinate the results across the nodes Pleiades supercomputer facility3 . For solving the con-
to drive them into consensus. The nodes update the vex problems at each iteration, we have4 used the convex
consensus variable z iteratively, based on their local optimization toolbox CVX for Matlab .
3 [Link]
2 [Link] 4 [Link]

4
 2 4000 4000

Average instantaneous

Average instantaneous
3500 Observed fuel consumption Observed fuel consumption
1.5
Normalized mean

Predicted fuel consumption Predicted fuel consumption

fuel consumption

fuel consumption
squared error

3000
1 3000

0.5 2500
2000
0 2000

−0.5 1500 1000

0 500 1000 1500 0 500 1000 1500 0 500 1000 1500 2000
Number of test flights Time (in seconds) in cruise Time (in seconds) in cruise

(a) Squared error for all flights (b) Outlier flight fuel consumption (c) Normal flight fuel consumption

Figure 2: Fuel flow study on CarrierX dataset. Fig. (a) shows squared error for all test flights, the 3-σ bound and
flights which cross the threshold. Fig. (b) shows the observed and predicted fuel flow of top ranked anomalous
flight. Fig. (c) shows the same for a normal flight.

Fig. 1 shows the sample dataset generated from gression in the primal. Our formulation is paralleliz-
a linear model following y = w × x + noise, where able among a number of computing nodes connected to
w is the weight of the regression model. We have a central computing node. Empirical study show that
used 2 nodes in this experiment and, for each node, our algorithm is accurate and scalable, ideal for large
chosen a different w vector so that each node sees scale deployment. As future work, we plan to develop
a different data distribution. The data of the two asynchronous version of this problem for peer-to-peer
nodes are shown in two different colors (circle and plus architectures.
markers). Also shown in the figure are the models
(straight lines) formed by node 1 at different iterations Acknowledgements
of linear ParitoSVR algorithm. The material is based upon work supported by the
ARMD seedling fund from the National Aeronautics
4.1 Anomaly detection on CarrierX dataset We and Space Administration under Prime Contract Num-
use the linear ParitoSVR algorithm to detect anomalous ber NAS2-03144 awarded to the University of Califor-
fuel consumption in a commercial aircraft. We model nia, Santa Cruz, University Affiliated Research Center.
the average fuel flow as a function of 29 different
parameters that measure system parameters such as
lateral and longitudinal acceleration, roll and pitch References
angle, air pressure, and velocity, as well as external
parameters such as wind speed and direction. We have [1] V. V. N. Vapnik, The Nature of Statistical Learning
used all 1500 flights (≈ 4.5 million training instances) for Theory. Springer-Verlag New York, Inc., 1995.
a specific tail number for a particular year for training, [2] D. Bertsekas and J. Tsitsiklis, Parallel and distributed
and tested subsequent years’ flights for predicting fuel computation: numerical methods. Upper Saddle River,
consumption. Flights for which the mean squared errors NJ, USA: Prentice-Hall, Inc., 1989.
of the predicted instantaneous fuel consumption fall [3] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eck-
outside the 3-σ boundary of the average mean squared stein, “Distributed Optimization and Statistical Learn-
ing via the Alternating Direction Method of Multipli-
error, are tagged anomalous (σ is the standard deviation
ers,” Found. and Trends in Mach. Learn., vol. 3, no. 1,
of the mean predictions). Out of approximately 1800 pp. 1–122, 2011.
flights for a test year, 14 flights were determined to
be anomalous. Figure 2(a) shows the mean squared
errors for each of the flights in blue and the 3-σ bounds
in green. The instantaneous fuel flow for the top
ranked anomalous flight among these 14 flights is shown
in Figure 2(b). The red graph depicting observed
fuel flow is significantly higher than the predicted fuel
consumption, shown in blue.

5 Conclusion
In this paper we have proposed ParitoSVR — a par-
allel iterated optimizer which solves support vector re-

5
 Anomaly Detection Using Tripoint Arbitration Similarity Method

Aleksey Urmanov∗

Abstract can be used with a specified missed detection rate (type

The tripoint arbitration similarity method uses every point II error).
in a sample as an observer to evaluate the similarity of a Distributional and possibly other data-generating
pair of points of the sample. The similarity of the pair is assumptions and tuning of various critical parameters
aggregated over all observers in the sample. The resulting are required to use existing anomaly detection methods.
pairwise similarity matrix captures information about clus- For example, when using the Mahalanobis distance,
ters of similar points. An anomalous point is defined as an a mutivariate Gaussian assumption is made for the
observer point for which all points in the sample are pair- data generating mechanism. When using clustering,
wise similar. The method is independent of the underlying a number of clusters must be specified and a specific
joint distribution of the sample points and does not require cluster formation mechanism must be assumed.
the user to tune any parameters other than selecting the ap- The analysis is becoming more laborious when ob-
propriate distance function and setting the admissible false servations are represented by heterogeneous data. For
detection rate. The proposed method handles heterogeneous instance, a health monitoring system of a computing
data by computing a combined similarity score which is of infrastructure that provides cloud services must contin-
interest for many industrial, social, scientific, web, retail, fi- uously monitor diverse types of data about thousands
nance and health sciences applications. The work in progress of targets. The monitored data may include physical
on anomaly detection using the tripoint arbitration similar- sensors, soft error rates of communication links, data
ity method is reported. paths, memory modules, network traffic patterns, inter-
nal software state variables, performance indicators, log
1 Introduction files, workloads, user activities etc, all combined into a
heterogeneous observation describing a target within a
Anomaly/outlier detection is one of the practical prob-
time interval. An anomaly detection system must con-
lems of data analysis. Applications range from cleansing
sume all this data and alert the system administrator
of data in statistical hypothesis testing and modeling,
about anomalously behaving targets. In such environ-
performance degradation detection in systems prognos-
ments it is unpractical to expect that the system ad-
tics, workload characterization and performance opti-
ministrator will possess sufficient skills to set and tune
mization for computing infrastructures, intrusion detec-
various anomaly detection parameters.
tion in network security applications, medical diagnosis
The tripoint arbitration similarity-based anomaly
and clinical trials, social network analysis and market-
detection system is developed to address these new
ing, optimization of investment strategies, filtering fi-
challenges. It has the following features designed-in:
nancial market data, fraud detection in insurance and e-
commerce applications. Methods for anomaly detection • Makes no distributional or other assumptions about
utilize statistical approaches such as hypothesis testing the data-generating mechanism.
[1] and machine learning approaches such as one-class
classification and clustering [2]. See [3] for a review. • Operates without tuning of any parameters by the
An anomaly is defined qualitatively as an observa- user. The appropriate distance function is selected
tion that significantly deviates from the rest of the sam- based on the type of the data.
ple. To quantify significant deviation a model is created
that represents nominal observations and allows to com- • Detects anomalies with a desired false detection
pute deviation from it with a given false detection rate rate. The user may specify the admissible false
(type I error). In rare cases when instances of actual detection rate, otherwise < 1% is used by default.
outliers are available in quantities sufficient to create
a model describing the outlier observations, likelihood • Handles seamlessly observations composed of het-
ratio-based statistical tests and two-class classification erogeneous components (numeric, text, categorical,
time series, other) as long as an appropriate dis-
∗ Oracle Labs, San Diego. Email: [Link]@[Link] tance function is available for each data type.

6
2 Tripoint Similarity and Clustering Using this definition of the pairwise similarity ma-
Tripoint arbitration similarity method is based on a trix for a data set, the clustering problem can be
novel definition of similarity of data points. Consider formulated as follows. Given a set of points D =
a collection of samples x1 , x2 , . . . , xn in Rm with the {x1 , x2 , . . . , xn }, xi ∈ Rm , the problem is to par-
Euclidian distance dij = d(i, j) = d(xi , xj ) as the tition the set into an unknown number of clusters
closeness measure for two points. Given a pair of points, C1 , C2 , . . . , CL so that points in the same cluster are
(xi , xj ), and an arbiter point a ∈ Rm , the tripoint similar and points in different clusters are dissimilar.
arbitration similarity is defined as This clustering problem can be casted into an optimiza-
tion problem that can be efficiently solved using matrix
min(d(i, a), d(j, a)) − dij spectral analysis methods
(2.1) Sa (xi , xj ) =
max(min(d(i, a), d(j, a)), dij ))
(2.6) min J(C1 , C2 , . . . , CL )
Sa takes values between −1 and +1 with the
(2.7) SD (Cp , Cp ) ≥ 0, 1 ≤ p ≤ L
following interpretation. Sa (xi , xj ) = −1 means that
points xi and xj are completely dissimilar for the (2.8) SD (Cp , Cq ) ≤ 0, 1 ≤ p < q ≤ L
arbiter point a. Sa (xi , xj ) = +1 means the points
SD (Cp , Cq ) is the average of all pairwise similarities of
are completely similar. Sa (xi , xj ) = 0 means the
points from clusters Cp and Cq
arbiter point cannot decide whether the points similar or
(2.9)
dissimilar. All other non-zero values reflect the degree 1 X X
of similarity (positive values) or dissimilarity (negative SD (Cp , Cq ) = SD (xi , xj )
| Cp | | Cq |
values) of the pair for the arbiter. i:xi ∈Cp j:xj ∈Cq
For a set of arbiter points A = {a1 , a2 , . . . , al } the
similarity is aggregated over all arbiters and the objective function J is constructed to simulta-
neously satisfy min SD (Cp , Cq ) for 1 ≤ p < q ≤ L and
l max SD (Cp , Cp ) for 1 ≤ p ≤ L. One such objective
1X
(2.2) SA (xi , xj ) = Sak (xi , xj ) function is [5]
l
k=1
X SD (Cp , Cq ) SD (Cp , Cq )
Let D = {x1 , x2 , . . . , xn } be a random sample from (2.10) J= +
SD (Cp , Cp ) SD (Cq , Cq )
an unknown population. The empirical pairwise tripoint 1≤p<q≤L
arbitration similarity matrix of sample D is defined as
Dropping the constraints in (2.6) leads to the prob-
(2.3) SD = [SD\{xi ,xj } (xi , xj )] lem similar to the MinMaxCut formulation in [5] with
pairwise associations given by the tripoint similarity.
Since similarity (2.1) ranges from -1 to +1 for any To deal with the constraints in (2.6) an iterative
type of data, it is possible to combine similarities of dif- approach was adopted. At the initial iteration the orig-
ferent modalities of multimodal data into a single overall inal problem is solved by partitioning the data set into
similarity. One rule for combining modal similarities is two clusters using an appropriate objective function, for
to follow common sense (analogously to [4]). For modal example, the MinMaxCut objective function in (2.10).
similarities with the same sign, the overall similarity be- At the next iteration each of the two clusters is parti-
comes bigger than either of the modal similarities but tioned in two and so forth. At each iteration the con-
still remains ≤ 1. straints are checked. Violation of the constrains serves
(1) (1) (2) (2) as a stopping criterion for the iterations. The process of
Sa (xi , xj ) = Sa(1) (xi , xj ) + Sa(2) (xi , xj )
(2.4) splitting clusters is stopped when no more clusters can
(1) (1) (2) (2)
− Sa(1) (xi , xj ) · Sa(2) (xi , xj ) be split without violating the inter-cluster dissimilarity
constraint. This iterative procedure automatically pro-
When modal similarities have different signs, the duces the appropriate number of clusters. The tripoint
overall similarity is determined by the maximum abso- arbitration clustering method uses matrix spectral anal-
lute value but the degree of similarity or dissimilarity ysis results to iteratively find the appropriate number of
weakens. clusters by solving the problem (2.6).
(2.5)
Sa (xi , xj ) = 3 Anomaly Detection System
(1) (1) (2) (2) The proposed anomaly detection system relies on tri-
Sa(1) (xi , xj ) + Sa(2) (xi , xj ) point clusters to determine a possible global structure
(1) (1) (2) (2)
1 − min(| Sa(1) (xi , xj ) |, | Sa(2) (xi , xj ) |) in the data. Tripoint clustering finds automatically the

7
Figure 1: Outlier detection on artificial data with FAR Figure 3: Estimated sampling distribution of Sz .
< 1%. All observations that lie in the yellow region
will be detected as outliers given the nominal data that from which are considered as outliers with < 1% FAR
comprise two clusters. for the data set compose of two clusters shown with
blue crosses and green triangles. Tripoint clustering
finds automatically two clusters in the data set and
assigns cluster labels to the corresponding observations.
Any new observation z for which Sz (C1 ) > 0.5 or
Sz (C2 ) > 0.5 is a detected outlier or anomaly.
Figure 2 shows the results of clustering and anomaly
detection in time series which represent certain at-
tributes of monitored targets in a computing infrastruc-
ture. Tripoint clustering found 2 clusters in a pool of
17 targets shown on the left and right sides of the plot.
A new target (#18) when presented to the anomaly de-
tection system was correctly detected as anomalous.
And finally Figure 3 shows the estimated sampling
distribution of Sz for a multivariate Gaussian data set
Figure 2: Anomalous target detection using time series with n = 500 points. By varying n, the values of
data representing operating targets in a computing tα can be tabulated for different α for use in (3.11).
infrastructure. Target #18 is correctly detected as For more rigorous calculation of tα , the exact sampling
anomalously behaving compared to nominally behaving distribution of Sz can be determined through Monte-
targets #1-17. Carlo simulations or asymptotic distribution theory.

References
appropriate number of clusters and labels the observa-
tions with a cluster label l = 1, 2, . . . , L. The resulting [1] P.J. Rousseeuw and A. Leroy, Robust regression and
clusters C1 , C2 , . . . , CL constitutes the nominal model outlier detection. New York: Wiley (1987).
based on the sample. [2] B. Schlkopf, J.C. Platt, J. Shawe-Taylor, A.J. Smola,
An anomaly is defined as an observation z for which and R.C. Williamson, Estimating the support of a high-
all cluster-average similarities are positive or all points dimensional distribution, Neural computation, 13(7),
from clusters C1 , C2 , . . . , CL are pairwise similar on pp. 1443-1471 (2001).
average, i.e. [3] V.J. Hodge and J. Austin, A survey of outlier detec-
(3.11) tion methodologies, Artificial Intelligence Review 22.2,
1 X pp. 85-126 (2004).
Sz (Cl ) = Sz (xi , xj ) > tα , i, j : xi , xj ∈ Cl , [4] I. B. Sirodza, Quantum models and methods of artificial
| Cl | i,j
intelligence for decision-making and control, Nauch-
naya Mysl (2002), pp. 92.
where α is the desired false detection rate. The thresh- [5] C. Ding, et al., A minmaxcut spectral method for data
old tα is determined from Prob(Sz > tα ) < α . clustering and graph partitioning, Lawrence Berkeley
Figure 1 illustrates the area (in yellow color) points National Laboratory, Tech. Rep 54111 (2003).

8
 Measuring Anomalousness in Statistical Models

Thomas Veasey∗ Stephen Dodson∗

1 Introduction. sponding to unusual system or user behaviour given the

Broadly speaking, unsupervised anomaly detection ordering defined by the q-value.
techniques fall into two categories: distance based ap- Various characteristics are ubiquitous in the data
proaches, which look at the distance between points or sets we work with, and we highlight those which we
local density, for example [4, 7, 8], and statistical ap- have found are particularly important to capture in
proaches, which are usually based on robust estimators the statistical model in order to get accurate anomaly
and hypothesis testing, see for example [2]. In this pa- detection. Specifically, non-Gaussian distribution tails,
per, we study a statistical technique that can be used proper handling of integer data, brakes and/or highly
for unsupervised anomaly detection, based on a varia- variable data rates and seasonality. Furthermore, the
tion of the concept of the p-value of the value of a test data sets are typically very large: monitoring data for
statistic, see for reference [12]. large computer networks can comprise tens or even
We show how this quantity, which we will refer to hundreds of thousands of performance metrics; common
henceforth as the q-value of an event, is well defined and data related to network security, such as proxy logs
generates an intuitive measure of the events’ anomalous- have transaction rates in the thousands of events per
ness in the presence of distribution modes, which would second. Any time series model must be highly compact,
correspond to anomaly detection on clustered data. and fast enough to compute on these data volumes.
Furthermore, unlike many distance based approaches We have found that summarising the time series by
for anomaly detection, such as kNN proposed in [10], small numbers of statistics, such as the mean, minimum
it naturally captures the relationship between the num- and maximum of n metric values, allows us to scale
ber of items in a cluster and their anomalousness. We to enormous data volumes with little loss in detection
discuss how to compute the q-value for some specific performance. In fact, varying the resolution, varying
distributions and also numerical approaches that can n for our example, effectively provides different insight
be used to compute it for arbitrary distributions. into the data: different types of anomaly emerge for
Finally, we study a class of high dimensional different choices.
anomaly detection problems where the events of pri- We give results for a set of performance metrics
mary interest are statistically significant deviations in generated monitoring an internet banking system over a
one, or a small number of, the dimensions. For such three day period. This is around 12 GB and comprises
problems, many of the difficulties associated with high around 33,500 distinct metric time series.
dimensional anomaly detection, data sparsity, choice of
distance metric [1], runtime [6], and model size, can be 2 A Definition of Anomalousness.
circumvented for the proposed measure of anomalous- The q-value is defined for any statistical model for which
ness, without using dimension reduction techniques. In a distribution function exists. Specifically, the model
particular, we show how applying the q-value to order must be some random variable from a probability space
statistics on the individual dimension values leads to a (Ω, F , P ) to some measure space (X, A ) and there must
natural measure for solving exactly this problem. exist a measurable function f : X → R+ , where R+
The authors’ primary interest is in anomaly detec- denotes the non-negative real line with Borel algebra,
tion for application performance monitoring and net- which recovers the probabilities of the measurable sets
work security. The data sets are nearly always col- of X. We define the q-value of an event x ∈ X as:
lections of time series, and a common requirement for
anomaly detection in this context is to provide alerts q(x) = P ({y : f (y) ≤ f (x)})
about anomalous system behaviour. We discuss a deci- This is clearly well defined, since the closed interval
sion criterion, which we use to identify anomalies corre- [0, f (x)] is Borel measurable and so its preimage is
A measurable. Since q(x) is a probability it takes
∗ Prelert Ltd., 156 Blackfriar’s Road, London, UK. E-mail: values in the interval [0, 1]. Any subset of [0, 1] has
{tveasey,steve}@[Link] the usual strict total ordering of the reals, and so we

9
can define a strict weak ordering of events by their Proof. By definition we have that
anomalousness, i.e. x >a y if and only if q(x) < Z
q(y). In particular, anomalousness can be defined as q(x) = 1{f (y) ≤ f (x)}f (y)dy
some monotonic decreasing function of the q-value, for
example − log(q(x)). X  fi fix

fi
= 1 ≤ V (Bi )
It is interesting to compare this definition with V (B i ) V (B ix ) V (B i)
i
the p-value, which is always defined in terms of some X
test statistic, say T (x), of the event. In particular, = fi
the p-value is the probability of the test statistic ex- {i:di ≤d(x)}

ceeding its observed value given the null hypothesis:

where 1{·} denotes the indicator function, and in the
P (T (y) > T (x)|H0 ). In the case that the null hy-
last line we have defined di = fi/V (Bi ). We can interpret
pothesis is the data set is Gaussian distributed with
this summation as the fraction of items for which the
known mean m and covariance V , and the test statistic
density is less than or equal to the density at the item
is T (x) = kx − mk2 , the so called two-tailed test, then
x.
the definitions coincide, since the probability density
function is less than f (x) in exactly the region where We note, also, that any random variable can be ap-
the test statistic is greater than T (x). Note, however, proximated in distribution by a mixture of uniforms,
that this case, a single mode symmetric model, is one of since we may approximate the cumulative density func-
the few cases where the values coincide. Furthermore, tion by a sequence of piecewise constant functions. So
the q-value makes no specific appeal to a null hypothe- this result can be used, in conjunction with binary space
sis. The idea is, given a statistical model of a data set, partitioning, to estimate q-values for arbitrary small-
to define a quantity that naturally relates to the anoma- ish dimensional multivariate models. Also, the density
lousness of observed events, in much the same way as function for any reasonable mixture model describing
say then mean distance to k-nearest neighbours does for clustered data will be proportional to the fraction of
distance based anomaly detection. items in a mode, in the vicinity of that mode, so items
If the data set to be analysed for outliers has been from clusters with fewer items will naturally have lower
clustered then the corresponding statistical model will q-values.
be multimodal. In particular, each cluster will typically
generate a mode of the distribution, or local maximum 3 Numerical Schemes for Calculating q-values.
in the density function. The q-value for an item from
For many univariate distributions the q-value can be
such a data set is therefore the (Lebesgue) integral of
evaluated in closed form. Otherwise, efficient numerical
the density function over some region that (usually)
methods often exist for finding the density function level
contains multiple holes, corresponding to modes of the
sets. However, for general multivariate distribution and
distribution.
mixture models numerical methods must be used. We
The exact value for an item depends on the choice
review a couple of approaches. Perhaps the simplest
of statistical model. However, for reasonable models the
scheme for computing the q-value, if the statistical
value should be close to the fraction of items in lower
model can be sampled is the following: generate n
density regions. In particular, we show P that for any mix- independent samples of the distribution, say Yn = {y},
ture of uniform random variable, i fi × U (Bi ) where
and define:
fi denotes the fraction of items in an m-dimensional
cuboid Bi = [a1,i , b1,i ] × [a2,i , b2,i ] × ... × [am,i , bm,i ] and |{y ∈ Yn : f (y) ≤ f (x)}|
qn (x) =
U (Bi ) is a uniform random variable on that cuboid, then n
the q-value of an item x is exactly the fraction of items a.s.
for which the density f (y) = fiy /V (Biy ) is less than or We show that qn (x) −−→ q(x) as n → ∞. In fact, we
equal to f (x), where iz denotes the index of the box can compute the asymptotic error distribution.
which contains item z and Q V denotes the volume func- Proof. As before, define our system model to be a ran-
tion defined as V (Bi ) = j |bj,i − aj,i |.
dom variable Y with probability density function f .
Let, A(x) denote the A measureable set f −1 [{z : z ≥
0, z ≤ f (x)}], and 1{A(x)} denote the indicator func-
tion of A(x). Given a random sample y of Y then, by
definition, 1{A(x)}(y) = 1 with probability q(x) and
0 otherwise. Therefore, 1{A(x)}(Y ), which we under-
stand as 1{A(x)} ◦ Y , is a Bernoulli random variable

10
with success probability p = q(x). By definition, performance metrics. Similarly, many types of network
attack amount to a small set of users doing highly un-
n
|{y ∈ Yn : f (y) ≤ f (x)}| 1X usual things at a given instant. For example, a port scan
∼ 1{A(x)}(Y )
n n i=1 attack would correspond to a client sending requests to
an unusually large range of server port addresses on
Pn a host in a relatively short period of time. For these
Furthermore, i=1 1{A(x)}(Y ) ∼ B(n, p), i.e. it is
a binomial random variable with number of trials n problems, it is not necessary to try and model the dis-
a.s. tribution on the full space. Instead, we can accurately
and probability of success p. Noting that B(n, p) −−→
N (np, np(1 − p)) as n → ∞ it follows that model its marginals, for example the individual perfor-
mance metrics, or the distribution of a population for
n  
1X a.s. q(x)(1 − q(x)) particular attributes, such as unique server port address
1{A(x)}(Y ) −−→ N q(x), requests in a fixed time interval. Then compute q-values
n i=1 n
on these, and aggregate the individual q-values to get
an effective measure for overall anomalous at any given
In particular, it is normally distributed with mean q(x)
instance.
and variance q(x)(1 − q(x))/n. The variance is maximised
a.s. To understand why we must account for the number
when q(x) = 1/2, and so qn (x) −−→ q(x) as n → ∞ for
of dimensions in this aggregation process, consider the
all x and we are done.
simple case that each marginal is a Gaussian. If there
This is just a particular Monte Carlo scheme for is no anomaly, we expect N independent samples from
R
evaluating the integral 1{f (y) ≤ f (x)}f (y)dy. It has a Gaussian, where N is the number of dimensions. In
one important advantage over the other schemes we dis- the case that N = 10, 000 the most extreme sample we
cuss: the same set of samples can be used for evaluating expect to see is about 4 standard deviations, where as
q(x) for any x. This means it is particularly well suited in the case that N = 10 the most extreme sample we
to Sequential Monte Carlo methods for estimating the expect to see is around 1.5 standard deviations. These
√
statistical model. Otherwise, the following methods will correspond to q-values of 1 − erf (4/ 2) = 6.3 × 10−5 and
√
have lower error for a given running time. 1 − erf ( / 2) = 0.13, respectively.
1.5

A recursive stratified sampling, such as the MISER Given a collection of N q-values, {qi }, a heuristic
algorithm [9], will yield lower variance for the same we have found to be useful for computing an aggregate
sample size. Further speedup can be obtained by q-value is to compute the q-value on the order statistic
importance sampling. In particular, the integrand is q (N ) (x) = P ({y : fX (N ) (y) ≤ fX (N ) (x)})
identically zero where f (y) > f (x). Therefore, if
N −1
we are using
P a mixture model, with density function Here, fX (N ) (y) = 1!(NN−1)!
!
(F (y) − F (−y)) f (y) de-
f (y) = i πi fi (x) and can compute the regions {Ri }, notes the distribution of the most extreme sample as a
bounded by zi : zi = fi−1S (f (x)/πi ) , we should only function of y, from a collection of N independent identi-
sample outside the region i Ri . Finally, we note for cally distributed samples from a symmetric single mode
a mixture of uniforms approximation, then storing the distribution, and x is the value of the most extreme sam-
region densities in a red-black tree augmented with the ple. In our case, we are not interested when the smallest
fraction of items in each subtree means it is possible to q-value is too large, given the sample size. Therefore,
compute the q-value for a new item in O(log(N )) and we evaluate P ({y : |y| ≥ |x|}). This is equal to
update the data structure in O(log(N )) where N is the  N
2N 1
number of uniforms, see [5] for details. [(2t − 1)]FX (x) = 1 − 1 − min{qi }
2N i

4 Statistically Significant Anomalies in Low In particular, we set our aggregate q-value to be 1 −

Dimensional Subspaces. (1 − mini {qi })N . Note, 1 − FX (x) = mini {qi }/2 follows
For many problems in application performance moni- from the assumptions about sample distributions and
toring and network security, anomalies of particular in- the definition of x. This result can be generalized
terest are statistically significant deviations in a small to compute the aggregate q-value from the M most
number of the raw measurement dimensions. For exam- extreme samples for N dimensions under the same
ple, if a system has ten thousand performance metrics, assumptions.
which might comprise average response times of differ-
ent database queries, responses per interval, errors per 5 Test Data and Methodology.
interval and so on, a system problem is likely to man- The data set we analysed was gathered by the CA
ifest itself as highly unusual values in a subset of the APM product monitoring three servers of an internet

11
banking site. Every performance metric is reported at noise ratio of around 65dB; however, this was not sig-
60s intervals, although some record transactions and are nificant enough to result in user noticeable system per-
not necessarily available at this granularity. It contains formance degradation.
33,159,939 distinct records and 33,456 distinct time We generated two alerts, corresponding to these two
series. The data cover a period of 72 hours and so the incidents for this data set. Our algorithm to generate
total data rate is around 500,000 values per hour. There alerts from the raw aggregate q-values is based on both
are 38 categories of metric; these include responses per the signal-to-noise and the historical quantiles of the
interval, average response time, errors per interval, stall aggregate q-value, for which we use the data structure
counts and average result processing. Note that various proposed in [11].
categories are split out by SQL command, host and so
on, which accounts for the total number of distinct time References
series.
For this data set, we found it was sufficient to as-
sume that the series were stationary. For other prob- [1] C. C. Aggarwal, A. Hinneburg, and D. A. Keim,
lems, capturing diurnal and weekly variation is impor- On the Surprising Behavior of Distance Metrics in
tant, for which we use radial basis function interpolation High Dimensional Space, Proc. Int. Conf. on Database
to fit the periodic temporal patterns. We chose to fit Theory, (2001), pp. 420–434.
[2] V. Barnett, and T. Lewis, Outliers in Statistical Data,
either a Gaussian distribution with unknown mean and
Third Edition, John Wiley & Sons, Chichester, 2006.
precision, a gamma distribution with unknown shape
[3] C. M. Bishop, Pattern Recognition and Machine Learn-
and rate, or a log-normal distribution with unknown ing, Springer, 2006.
location and scale to the (assumed) stationary distribu- [4] M. M. Breunig, H. P. Kriegel, R. T. Ng, and J. Sander,
tion of each time series. We use standard Bayesian tech- LOF: Identifying Density-Based Local Outliers, SIG-
niques to estimate the parameters, and Bayesian model MOD00: Proc. ACM SIGMOD Int. Conf. on Manage-
selection to choose among the models, see [3] for details ment of data, (2000), pp. 93–104.
on Bayesian model selection. Finally, on this data set [5] T. H. Cormen, C. E. Leiserson, R. Rivest, and C. Stein,
we found it was very important to accurately account Introduction to Algorithms, Third Edition, The MIT
for time series comprising integer data with low vari- Press, 2009.
ation, in particular, series for which particular integer [6] M. Datar, N. Immorlica, P. Indyk, and V. S. Mirrokni,
Locality-Sensitive Hashing Scheme Based on p-Stable
values have significant probability. Such data are gen-
Distributions, Proc. ACM Symp. on Computational
erally badly modelled by continuous distributions. We Geometry, (2004), pp. 253–262.
automatically detect this case, and model these data us- [7] H. Fan, O. Zaı̈ane, A. Foss, and J. Wu, A Nonpara-
ing a latent variable. In particular, we assume that the metric Outlier Detection for Efficiently Discovering
observed values are described by X + U ([0, 1]), where Top-N Outliers from Engineering Data, Proc. Pacific-
we estimate X, and U ([0, 1]) denotes a uniform random Asia Conf. on Knowledge Discovery and Data Mining
variable on the interval [0, 1]. (PAKDD), (2006), pp. 557–566.
A large anomaly manifested itself as system perfor- [8] E. M. Knorr, and R. T. Ng, Algorithms for Mining
mance degradation during the interval 32 to 35 hours af- Distance-Based Outliers in Large Datasets, Proc. of
ter the start of the data set. In terms of the raw anomaly the 24th International Conference on Very Large Data
scores, which were obtained by aggregating individual Bases, (1998), pp. 392–403.
[9] W. H. Press, and G. R. Farrar, Recursive stratied sam-
time series q-values, this corresponded to a signal-to-
pling for multidimensional Monte Carlo integration,
noise ratio of around 330dB, where the noise level was Computers in Physics, vol. 4, (1990), pp. 190–195.
taken as the median aggregate q-value. If the time se- [10] S. Ramaswamy, R. Rastogi, and K. Shim. Efficient
ries are ordered by their q-values at that time, then 560 Algorithms for Mining Outliers from Large Data Sets,
of the 33,456 time series are significantly anomalous. Proc. ACM SIGMOD Int. Conf. on Management of
These results indicated there was an operational issue data, (2000), pp. 427–438.
with a specific component of the backend, which re- [11] N. Shrivastava, C. Buragohain, D. Agrawal, and
sulted in the response time of a subsection of the website S. Suri. Medians and Beyond: New Aggregation Tech-
(6 JSPs) having dramatically increased response times. niques for Sensor Networks, Proc. of the 2nd Int.
In addition, there was a precursor to the main anomaly, Conf. on Embedded Network Sensor Systems, (2004),
at 27 hours after the start of the data set, which pro- pp. 239–249
[12] R. E. Walpole, R. H. Myers, S. L. Myers, and K. Ye,
vided the system administrators with early warning of
Probability and Statistics for Engineers and Scientists,
the specific problem before the main failure. This was Ninth Edition, Prentice Hall, 2011.
detected in the performance metrics with a signal-to-

12
 Identifying Precursors to Anomalies Using Inverse Reinforcement Learning∗

Vijay Manikandan Janakiraman† Santanu Das‡ Bryan Matthews§ Nikunj Oza¶

Abstract introducing some background in inverse reinforcement

In this paper, we consider the problem of discovering can- learning, using its solution to perform value function es-
didate precursors to anomalies in a set of time sequenced timation and using the optimal value function, discover
data. Typical scenarios involving time sequential data in- precursors.
clude dynamical systems and general monitoring systems.
In such scenarios, a precursor could be any event that fre- 2.1 Inverse Reinforcement Learning The goal of
quently precedes a given event of interest. Anomalies are inverse reinforcement learning (IRL) is to determine
rare but significant events in time series data and identify- the underlying reward function using observed behavior
ing precursors to anomalies is vital in proactive management. of the agent making decisions in a Markov Decision
In this work, an inverse reinforcement learning (IRL) based Process (MDP). A finite MDP is a tuple (S, A, Ps,a ,
method is formulated to succinctly represent the nominal γ and R(s)) where S is a state space with n states, A
behavior and identify sequences that preceded the anoma- is an action space with k actions, {Ps,a } are the state
lous events. A preliminary evaluation is performed on flight transition probabilities corresponding to an action a at
recorded data identifying challenges and future directions for state s, γ ∈ [0, 1) is the discount factor, R(s) ∈ R is the
application. underlying reward function. A policy π can be defined
as any map π : S 7→ A and the corresponding value
1 Introduction function at any state s1 can be given by
In many applications including finance, study of natu- (2.1) V π (s1 ) = E[R(s1 ) + γR(s2 ) + γ 2 R(s3 ) + ...|π]
ral calamities and extreme weather, network security [1]
etc., finding precursors to an event of interest (a phe- where the expectation is over the distribution of state
nomenon) is a task of high importance. The knowledge sequences (s1 , s2 , s3 , ...) following the policy π starting
about precursors to these phenomena can be vital to from s1 .
proactive management of risk. If precursor events could Given the setting above, the goal of standard re-
be identied, appropriate alarming mechanisms can be inforcement learning is to determine a policy π ∗ that
designed to either prevent or at least minimize the dele- maximizes V π (s) among all policies for all s ∈ S. When
terious consequences of the phenomenon. Anomalous the agent’s reward function is known, this task can be
events are rare but significant events which in many achieved using existing techniques for value function es-
cases, lead to an abnormal behavior or a risky situa- timation [2]. However, in several situations, the agent’s
tion. In such cases, it is important to analyze and iden- behavior is not completely known, i.e., the reward func-
tify precursors that lead to anomalies for proactive risk tion cannot be defined easily. In such situations, the
management. This paper considers anomalies in time expert’s observed behavior can be used to either recon-
sequenced data and attempts to discover candidate pre- struct the underlying reward function as in the case of
cursors to the anomalous events. inverse reinforcement learning [3] or construct optimal
policies directly as in the case of apprenticeship learning
2 Discovering Precursors to Anomalies [4].
In this section, an algorithm using inverse reinforcement Assuming availability of sampled trajectories (rele-
learning is proposed to identify candidate precursors to vant to the problem involving time series in this paper),
anomalies in time series data. The section proceeds by the IRL problem can be posed as in [3]. The sampled
trajectories can be considered as demonstrations of both
the expert and non-expert acting in the MDP. Using the
∗ Supported by the NASA System-wide Safety and Assurance
trajectories, the value functions of the expert and non-
Technologies (SSAT) Project.
† UARC, Nasa Ames Research Center, Moffett Field, CA expert policies can be determined as follows. Let the
‡ Verizon, Palo Alto, CA unknown reward function be parameterized as
§ SGT Inc., Nasa Ames Research Center, Moffett Field, CA
¶ Nasa Ames Research Center, Moffett Field, CA (2.2) R(s) = α1 φ1 (s) + α2 φ2 (s) + .. + αd φd (s)

13
where the φi represent the features of the reward etc. In this paper, considering sample time series from
function. The expert value function following policy a policy as monte carlo samples, the value function is
πE at state s1 can be given by approximated as follows. For each policy πj including
the expert policy πE , a sample trajectory is used to
(2.3) V πE (s1 ) = E[R(s1 ) + γR(s2 ) + ...|π] identify the state sequences and using the reward func-
= E[α1 φ1 (s1 ) + α2 φ2 (s1 ) + .. + αd φd (s1 ) tion obtained above, the values of every state in S is
updated. This is repeated for several trajectories from
+ γα1 φ1 (s2 ) + γα2 φ2 (s2 ) + .. + γαd φd (s2 ) + ..|πE ]
the selected policy and the average returns are stored
as state values.
= E[α1 (φ1 (s1 )+γφ1 (s2 )+..)+α2 (φ2 (s1 )+γφ2 (s2 )+..)
+ ... + αd (φd (s1 ) + γφd (s2 ) + ..)|πE ] 2.3 Precursor identification The value function of
the expert policy πE obtained above can be used to
compare a non-expert behavior to identify a possible
= α1 E[(φ1 (s1 ) + γφ1 (s2 ) + ..)|πE ]+ precursor sequence. V πE can be thought of as the
+ α2 E[(φ2 (s1 ) + γφ2 (s2 ) + ..)|πE ] expert’s value function and any action that is greedy
+ ... + αd E[(φd (s1 ) + γφd (s2 ) + ..)|πE ] with respect to the expert’s value function gives the
optimal policy π∗ [2]. Let the greedy action at s be
a∗ (s) and the corresponding value be V π∗ (s). In our
= α1 λ1 + α2 λ2 + .. + αd λd problem involving time series, the time sequence and
the physics of the problem can be used to restrict the
where λi represent the feature expectations, i.e., the state space for searching optimal actions in some cases.
value function if the reward function is composed of A given test time series can be analyzed as follows.
φi (s) only. After calculating the feature expecta- Using the state sequences of the test data and the
tions knowing the state sequences, the value function obtained reward function, the state values V πtest can
can be defined as a function of the unknown α = be estimated. By comparing the V πtest with V π∗ , we
[α1 , α2 , .., αd ]T as follows. can indirectly evaluate the actions taken by the agent
d in the test trajectory. Let
X
πE
(2.4) V (α) = αi λi (2.7) ∆V = V πtest (s) − V π∗ (s)
i=1
and if ∆V ≤ 0, then it would mean that a sub-optimal
Similarly, by knowing the sequence of states (trajecto-
action has been taken by the agent executing the test
ries) for J sub-expert/non-optimal policies, the V πj (α)
policy and by comparing over the state sequence, we
can be calculated. The objective of IRL is to deter-
can identify a sequence of bad actions by the agent.
mine the coefficients αi so that V πE (α) ≥ V πj (α) for
As defined earlier, an optimal action is one that cor-
j = 1, 2, ..J. A linear programming problem can be
responds to a nominal time series while a non-optimal
solved for αi as follows
action would correspond to an anomalous sequence as
XJ defined in the IRL problem. It should however be noted
(2.5) min ζj that the test policy is evaluated just based on one time
α
j=1 series and hence not an expectation. However, the goal
is to identify the level of sub-optimality in the state se-
 quences specifically executed by the test trajectory to
 πj πE
V (α) − V (α) − ζj ≤ 0 identify the precursor and not for the policy in general.
(2.6) subject to ζj ≥ 0, j = 1, 2, .., J This assumption needs to be analyzed more in detail


|αi | ≤ 1, i = 1, 1, .., d and will be considered in the future. Further, if the
action space is well defined, instead of comparing the
2.2 Value Function Estimation The IRL problem value functions as above, the actions of the test agent
gives an optimal α which gives a model of the underly- can be directly compared against the optimal actions of
ing expert’s reward function. The reward function can the expert and precursors can be identified by noting
then be used to determine the expert’s value function their difference.
using a regular reinforcement learning algorithm. Any
of the methods described in [2] such as dynamic pro- 3 Application to Flight Anomalies
gramming, monte carlo or temporal difference depend- In this section, the IRL based precursor discovery al-
ing on availability of the system model, ability to sample gorithm is evaluated on flight time series data sets ob-

14
tained from a FOQA (Flight Operations Quality As- on a hold-out data set. Following section 2.1, the ∆V
surance) archive. Typical FOQA parameters consist for a given test flight is calculated. A negative value for
of both continuous and discrete (categorical) data from ∆V indicates that the given test flight performs inferior
the avionics, propulsion system, control surfaces, land- to the optimal policy π∗ and a negative rate of ∆V indi-
ing gear, the cockpit switch positions, and other critical cates a sequential inferior behavior. These two features
systems. Each flight record can have up to 500 param- are used in defining precursor candidates for the given
eters in the form of time sequences and are sampled at test flight. It should be noted that the problem in hand
1 Hz. uses FOQA data that only records the state of the flight
Flight anomalies are of significant interest within and no explicit information about the intentions/actions
the NASA System-wide Safety and Assurance Technolo- of the agent (a pilot) is available and hence we were re-
gies (SSAT) project to assess the health of large com- stricted to comparing the value functions as mentioned
mercial fleets of aircraft. In this paper, flights that in section 2.3
violated exceedance thresholds on computed air-speed Using the identified precursor sequences of a given
are considered as operational anomalies. A specific ex- test flight in terms of the states s, the FOQA historical
ceedance defined as computed air-speed above a cer- data can be used to identify the flight parameters that
tain threshold (in knots) at an altitude of 1000 feet is are abnormal. The identified precursor sequence points
considered an operationally significant high-energy ap- to a section of the flight prior to the adverse event where
proach. The goal of this study is to discover precursors interesting precursor events can be discovered. By mod-
to such high-energy approach flights [5] for use in proac- eling a nominal distribution of the FOQA parameters,
tive flight management. The data set consists of about any abnormality can be detected by comparison against
20000 nominal flights (flights that did not violate the the nominal. The identified abnormal parameters may
exceedance and considered optimal with respect to the contain information about possible factors that lead to
exceedance) and about 250 anomalous flights. the adverse event. This is algorithmically analyzed and
validated by a domain expert.
3.1 Discovery of candidate precursor sequences
The FOQA raw data consists of more than 400 parame- 4 Results and Discussion
ters recorded as time sequences during the flight. How- In this section, a high-energy approach flight is analyzed
ever, to overcome the curse of dimensionality in solving for precursors from 35 nautical miles until touchdown.
the Markov decision process in the IRL, the FOQA data Figure 1 shows the evolution of the fight in terms of
is abstracted to represent the various events happening the parameters reported by the algorithm as candidate
in a flight using a high level parameter such as the air- precursors. The blue shaded region represents the nom-
craft energy. With the given definition of an anomaly, inal distribution of that parameter (99 percentile of the
the flight data is considered as a sample from an expert non-exceedance flights) The green shaded regions of the
policy (πE ) if it doesn’t flag the exceedance or a sample figure represent the sequence of precursors (a precursor
from a non-expert policy (πj ) if it flags the exceedance. window) as identified by comparing the flight’s state
A reward function R(s) can be defined as a linear com- values to V π∗ .
bination of several Gaussian functions defined with re- It can be observed from Figure 1 that the computed
spect to the states s. It has to be noted that the state air-speed of the flight is high compared to the nominal
definition is given by s = [E D]T where E represents distribution of the non-exceedance flights indicating
the kinetic energy of the aircraft while D represents the that the test flight is indeed an example of a high-energy
distance in nautical miles to touchdown. The reward approach. Further, out of the 56 chosen parameters
function R(s) can be represented as in equation (2.2) from the FOQA list, only 11 were listed as possible
where φi could represent Gaussian functions with mean precursor parameters as these parameters were out
µi and spread σi and d represents the total number of of the nominal distribution in the precursor window.
Gaussian functions in the state space. Using the re- The algorithm also reported ground speed which is
ward function with unknown coefficients αi the value correlated with the computed air-speed, vertical speed,
function of each trajectory is calculated and the IRL stabilizer position, engine speed, flight director specified
problem is solved as in section 2.1. The optimal value speed etc. However, a close look at the discrete
of α gives a model of the underlying reward function parameters reported by the algorithm gives a clear
that when used to solve the associated MDP, results in picture of the actions responsible for the anomalies,
maximum possibility of avoiding the given exceedance. i.e., the landing gear has been deployed a little earlier
The model hyper-parameters including d, µi , σi are de- compared to nominal flights and the flaps were deployed
termined based on cross-validating the learned model very late causing the aircraft to slow down late leading

15
Figure 1: Figure showing the test flight trajectory (black curve) along with the precursor window (green region)
as identified by the IRL algorithm. The nominal distribution of the continuous parameters such as computed
air-speed, engine RPM, stabilizer position, vertical speed and ground speed are shown in blue - light blue region
represents 0 - 99 percentile while dark blue region represents 25 - 75 percentiles. The nominal distribution of
discrete variables including landing gear, flaps, auto speed control are shown by blue curve with markers indicating
the probability of the variable having a value of 1.

to the exceedance (In the discrete plots, the marked be identified for evaluation of this algorithm in future.
blue curve represents the probability that a nominal Finally, some of the underlying hypotheses/assumptions
discrete event takes a value of 1). Both these factors of the algorithm will be tested in future.
were validated by domain experts as probable precursors
to the high speed exceedance. The initial high computed References
air-speed followed by a lack of optimal action (which is
to deploy landing gear and flaps on time) in this case [1] J. B. D. Cabrera, L. Lewis, X. Qin, W. Lee, R. Pras-
can be concluded as a valid precursor for flights violating anth, B. Ravichandran, and R. Mehra, “Proactive de-
the high-speed exceedance at 1000 feet altitude. tection of distributed denial of service attacks using
mib traffic variables-a feasibility study,” in Integrated
5 Conclusions Network Management Proceedings, 2001 IEEE/IFIP
International Symposium on, 2001, pp. 609–622.
In this paper, a novel method to discover precursors to [2] R. S. Sutton and A. G. Barto, Introduction to Rein-
anomalies has been formulated using inverse reinforce- forcement Learning, 1st ed. Cambridge, MA, USA:
ment learning. It is argued that a value function of a MIT Press, 1998.
non-expert, if compared against the optimal value func- [3] A. Y. Ng and S. Russell, “Algorithms for inverse
tion of an expert, can be used to identify instances of reinforcement learning,” in in Proc. 17th International
a “bad” or sub-optimal actions/situations in time se- Conf. on Machine Learning. Morgan Kaufmann, 2000,
ries data. A high dimensional FOQA time series data pp. 663–670.
has been abstracted and used for preliminary evalua- [4] P. Abbeel and A. Y. Ng, “Apprenticeship learning
via inverse reinforcement learning,” in Proceedings of
tion of the algorithm. The results indicate that the al-
the Twenty-first International Conference on Machine
gorithm indeed finds the precursors that were validated
Learning, ser. ICML ’04. New York, NY, USA: ACM,
by a domain expert. Although the analysis on a cou- 2004, pp. 1–.
ple of flights gave us promising results, the algorithm is [5] S. Das, L. Li, A. Srivastava, and R. J. Hansman, “Com-
at infancy and requires extensive validation for which parison of algorithms for anomaly detection in flight
data sets with ground truth information about the pre- recorder data of airline operations,” in 12th AIAA Avi-
cursors and anomalies are required. Also, for precursor ation Technology, Integration, and Operations (ATIO)
identification, an appropriate performance metric will Conference. American Institute of Aeronautics and
Astronautics, 2012.

16