# Hybrid Smoothing for Anomaly Detection in Time Series

Matthew Hofkes
Department of Statistics, Colorado School of Mines
and
Douglas Nychka
Department of Statistics, Colorado School of Mines  The authors gratefully
acknowledge MoWater, The National Science Foundation, and Denver Water.

###### Abstract

Many industrial and engineering processes monitored as times series have
smooth trends that indicate normal behavior and occasionally anomalous
patterns that can indicate a problem. This kind of behavior can be modeled by
a smooth trend such as a spline or Gaussian process and a disruption based on
a sparser representation. Our approach is to expand the process signal into
two sets of basis functions: one set uses $L_{2}$ penalties on the
coefficients and the other set uses $L_{1}$ penalties to control sparsity.
From a frequentist perspective, this results in a hybrid smoother that
combines cubic smoothing splines and the LASSO, and as a Bayesian hierarchical
model (BHM), this is equivalent to priors giving a Gaussian process and a
Laplace distribution for anomaly coefficients. For the hybrid smoother we
propose two new ways of determining the penalty parameters that use effective
degrees of freedom and contrast this with the BHM that uses loosely
informative inverse gamma priors. Several reformulations are used to make
sampling the BHM posterior more efficient including some novel features in
orthogonalizing and regularizing the model basis functions. This methodology
is motivated by a substantive application, monitoring the water treatment
process for the Denver Metropolitan area. We also test the methods with a
Monte Carlo study designed around the kind of anomalies expected in this
application. Both the hybrid smoother and the full BHM give comparable results
with small false positive and false negative rates. Besides being successful
in the water treatment application, this work can be easily extended to other
Gaussian process models and other features that represent process disruptions.

Keywords: Multiple change point detection, smoothing spline, process control,
LASSO, $L_{2}$/$L_{1}$ penalty, Bayesian, hybrid smoothing

##  1 Introduction

Many industrial and engineering processes that are monitored as times series
have smooth trends that indicate normal behavior and occasionally anomalous
patterns that can indicate a problem. In this work we combine a flexible model
for a smooth curve with the potential to include more discontinuous
interruptions. Although there are many possibilities for introducing anomalous
behavior, this work is motivated by a specific problem in water filtration for
the city of Denver. Here normal operation appears to be well described by a
Gaussian process (GP) and the anomalous intervals exhibit abrupt
discontinuities added to this GP baseline. From a frequentist point of view we
characterize a GP model as involving a basis function expansion with quadratic
penalties on the basis coefficients. The anomalous signal, however, tends to
have discontinuities that are limited in number and so suggests a basis
representation where many coefficients are small or zero. This sparse
representation can be captured by $L_{1}$ penalties on the basis coefficients.
Thus we are led to a hybrid process signal that is a sum of two components,
based on two separate sets of basis functions, where one set uses $L_{2}$
penalties to identify the trend while the other set uses $L_{1}$ penalties to
identify distrubances. The development of the LASSO and related models $L_{1}$
penalties have been very successful in building sparsity into a set of
parameters in an efficient and practical manner. However, the combination of
$L_{2}$ and $L_{1}$ models for the basis coefficients to identify a mixed
signal, to our knowledge, has not been exploited.

In the context of monitoring a process one often has a record of past system
performance and some prior knowledge of the normal and anomalous behavior.
Thus it is natural to approach this problem also from a Bayesian perspective
in order to take advantage of reference distributions and also to quantify the
uncertainty for the measurements over a given period. A Bayesian hierarchical
model (BHM) also builds in more symmetry between the trend and anomaly
signals. The trend is based on a multivariate normal distribution prior for
the basis coefficients and the anomaly uses a prior for the coefficients that
is double exponential (Laplace). Besides providing a useful solution to our
motivating application, the comparison between frequentist and Bayesian
approaches is interesting in that both are able to accurately detect anomalous
behavior and also have low rates of false positives. They accomplish this
identification, however, through different means.

The challenge for any urban area is to provide a safe and ample supply of
drinking water. The scale of this effort is often not appreciated. Denver
Water, providing clean drinking water to all of the Denver, Colorado
metropolitan area, must filter and treat an average of 180 million gallons a
day. The scale of the Denver Water system precludes an operator from tracking
in real time every filtration cycle on each filter, and therefore, it is
important to provide a statistical method to flag unusual filter cycles for
further examination. In this work we focus on the time series of the variable
headloss to track filter performance.

Large scale water filtration such as that used in the Denver facility is
achieved by gravity flow through filter media consisting of anthracite, sand,
and gravel. Over the course of a day, as sediment accumulates in the filter,
pressure at the filter outlet slowly increases, and the water level begins to
rise. Headloss is a measure of the filtration system water pressure associated
with rising water levels. When headloss in a filter reaches a certain level,
an operator shuts down the filter and a backwash is initiated in order to
clean the filter before a new cycle begins.

Figure 1: Five headloss cycles are seen here. The third is anomalous.

Headloss measurements in normal cycles increase steadily over time, however
cycle length, height, and concavity vary significantly between cycles. The
first four and the last four cycles in Figure 1 are normal. The third cycle is
anomalous. Anomalous cycles can take a variety of forms but generally consist
of at least one large jump in the measurements of headloss over time. We refer
to these jumps as change points, anomalies, or disturbances in the smooth
trend. Two other examples of anomalous cycles containing change points can be
found in Figure 2.

There are other ways in which a cycle can be outside of the norm including
being very short or not reaching the 8-8.5 units of headloss value before a
backwash is initiated, however those cycles can be easily found and flagged
using other methods. In contrast, our analysis looks for disturbances in the
trend and can often be subtle. Although the development of our approach is
focused on this type of time series, a general model that combines a GP with a
non Gaussian and sparse process may be useful in other applications. In
particular, this method can be adapted to other change point detection
problems in a time series. Moreover, it may be surprising to the reader that
this type of separation is possible along with the ability to sample from the
Bayesian posterior of either signal.

Figure 2: Two additional anomalous cycles with a disturbance in the mean.

In Section 2, we detail a general model, the relationship between cubic
smoothing splines and GPs, and LASSO regression. We then develop a model for
curve fitting that combines the $L_{1}$ and $L_{2}$ penalties, and outline an
optimization algorithm for estimating the coefficients. This frequentist
approach requires choosing penalty values, and we adapt a machine learning
idea for this task. In Section 3, we adopt a Bayesian hierarchical approach
and choose priors that are comparable to the $L_{1}$ and $L_{2}$ penalties in
Section 2. We then discuss adjustments that can be made which result in a full
Gibbs sampler and overcome the confounding that often occurs in Bayesian
models mixing fixed and random effects (Hodges and Reich (2010)). In Section
4, we compare the effectiveness of the frequentist and Bayesian approaches
based on a Monte Carlo study for the filter data. Finally, in Section 5, we
compare the effectiveness of the two approaches in identifying anomalous
cycles in the Denver Water headloss data and a present a general assessment of
the each approach.

##  2  A Frequentist Approach to Identifying Anomalies

In order to identify the change points in any anomalous cycle, we consider our
data to have come from the additive model

$y(t_{i})=f(t_{i})+g(t_{i})+\epsilon_{i}$ (1)

where $\mathbf{y}$ is a vector of observations, $f(t)$ is a smooth function,
$g(t)$ is a rough function, $\epsilon_{i}\sim N(0,\sigma^{2})$, and
$\\{t_{i}\\}$ are $n$ time points. The observation points are equally spaced
in our application but are not required to be. In our case we identify the
smooth function, $f(t)$, as the normal trend of the process over time, and the
rough function, $g(t)$ represents possible anomalies that occur having abrupt
change points. Figure 3 illustrates the general model and a typical anomaly
pattern that is encountered in the filter application.

Figure 3: An example of the components in the additive model in Equation 1.

The key to separating these components lies in our hybrid smoothing model for
$f$ and $g$. We will represent $f$ as cubic smoothing spline and $g$ as being
dominated by a sparse collection of step changes. Moreover, our BHM
formulation has a posterior mode that is comparable to the hybrid smoother.

###  2.1 Cubic Smoothing Splines and GPs

A cubic smoothing spline was originally proposed as a method of curve fitting
via a smoothness penalty on the second derivative of the model. The spline
curve is continuous in its first and second derivatives and is linear beyond
its outer knots. We begin with the additive model omitting $g$.

$y(t_{i})=f(t_{i})+\epsilon_{i}$

The cubic spline optimization can be is written

$\begin{array}[]{rclcl}\displaystyle\min_{f\in{\cal
H}}&\lx@intercol\sum_{i=1}^{n}(y(t_{i})-f(t_{i}))^{2}+\omega\int(f(t)^{\prime\prime})^{2}dt\hfil\lx@intercol\end{array}$
(2)

Formally ${\cal H}$ is a particular Hilbert space of functions but
operationally can taken to be all functions where the smoothness penalty is
finite. This solution can be characterized with a finite representation. There
is a set of $n$ basis functions, $\\{\phi_{k}(t)\\}$, and coefficients,
$\\{c_{k}\\}$, with $f(t)=\sum_{k=1}^{n}\phi_{k}(t)c_{k}$ and the coefficients
are found as

$\begin{array}[]{rclcl}\displaystyle\min_{\mathbf{c}}&\lx@intercol\|\mathbf{y}-\Phi\mathbf{c}\|^{2}_{2}+\omega\mathbf{c}^{T}R\mathbf{c}\hfil\lx@intercol\end{array}$
(3)

where $\Phi_{i,k}=\phi_{k}(t_{i})$ and $R$ is a Gram matrix of the smoothing
penalty applied to the basis functions

$R_{j,k}=\int\phi_{j}(t)^{\prime\prime}\phi_{k}(t)^{\prime\prime}dt$

There are several common ways to choose these basis functions that will give
equivalent results but may result in different degrees of conditioning in the
matrix $\Phi$. In all cases the spline is found as a form of ridge regression
and the coefficients are found using the usual ridge regression formulation.
Finally we note that $\Phi$ has full rank and so one can reparametrize this
minimization in terms of the spline function values at the data locations.

$\begin{array}[]{rclcl}\displaystyle\min_{\mathbf{f}}&\lx@intercol\|\mathbf{y}-\mathbf{f}\|^{2}_{2}+\omega\mathbf{f}^{T}\Gamma\mathbf{f}\hfil\lx@intercol\end{array}$
(4)

Here $\mathbf{f}_{i}=f(t_{i})$ and $\Gamma=\Phi^{-T}R\Phi^{-1}$. Some algebra
shows that the minimizer in this form is
$\hat{\mathbf{f}}=(I+\omega\Gamma)^{-1}\mathbf{y}$ and can be interpreted as a
linear smoother applied to the time series in $\mathbf{y}$. A subtlety
introduced in this version is that the penalty
$\mathbf{f}^{T}\Gamma\mathbf{f}$ will be zero when $\mathbf{f}$ is a linear
function and so $\Gamma$ will have a rank of $n-2$.

An important property of a cubic spline and many other splines is that they
also have an interpretation as a prediction from a Gaussian process model. For
clarity here we focus on the stochastic interpretation of a cubic spline and
prediction at the observation locations, although we note that these ideas
hold in general for the family of multidimensional thin plate splines and
prediction at arbitrary locations.

Assume that $f(t)=P(t)+h(t)$ with $P(t)$ is a linear function in $t$ and
$h(t)$ a mean zero, Gaussian process with covariance function

$COV(h(t),h(t^{*}))=k(t,t^{*})$

and define a $n\times 2$ regression matrix $\mathbf{X}$ with the $i^{th}$ row
given by $[1,t_{i}]$. That is, $\mathbf{X}$ is the usual regression matrix for
fitting a linear trend. In this setup we have the usual linear GP model

$\mathbf{y}=\mathbf{f}+\epsilon=\mathbf{X}\boldsymbol{\beta}+\mathbf{h}+\epsilon$

where $\mathbf{h}_{i}=h(t_{i})$. Now let $K$ be the covariance matrix of
$\mathbf{h}$ and it is a standard result that the best linear unbiased
estimator for f is given as

$\hat{\mathbf{f}}=\mathbf{X}\hat{\boldsymbol{\beta}}+K(K+\sigma^{2}I)^{-1}(\mathbf{y}-X\hat{\boldsymbol{\beta}})$

and where $\hat{\boldsymbol{\beta}}$ is the generalized least squares
estimator of $\boldsymbol{\beta}$. Under this model, the error vector,
$\mathbf{f}-\hat{\mathbf{f}}$ will be distributed $MN\sim(0,\Omega)$ where
$\Omega$ depends on $\sigma^{2}$, $X$, and $K$. The surprising result is that
one can also show the estimator, $\hat{\mathbf{f}}$ is also obtained as the
solution to a minimization problem.

$\begin{array}[]{rclcl}\displaystyle\min_{\boldsymbol{\beta},\mathbf{h}}&\lx@intercol\|\mathbf{y}-\mathbf{X}\boldsymbol{\beta}-\mathbf{h}\|^{2}_{2}+\omega\mathbf{h}^{T}K^{-1}\mathbf{h}\hfil\lx@intercol\end{array}$
(5)

Using a specific choice of covariance function, $k$, this inference for a
Gaussian process model will be the same as a cubic smoothing spline and was an
early landmark result connecting numerical splines to statistics. The
covariance function that gives a cubic spline is identified on page 39 of
Smoothing Splines for Observational Data (Wahba (1990)) and implemented with
the Tps.cov() function in the R package fields (Nychka et al. (2021)). Thus we
have connected the initial smoothing problem for functions to an estimator
that assumes a stochastic process model for the unknown function. This
equivalence is useful when interpreting the BHM in Section 3.

The last elaboration needed for our work is the extension of this
correspondence to a complete Bayesian model where all unknown quantities have
stochastic, prior distributions. Suppose one places a prior on
$\boldsymbol{\beta}$ as $MN(0,bI)$ and conditions on all other statistical
parameters, e.g. $\sigma^{2}$ and $K$. The conditional distribution of f given
$\mathbf{y}$ is then, again, multivariate normal. Moreover in the limit as
$b\rightarrow\infty$ this conditional distribution will converge to
$MN(\hat{\mathbf{f}},\Omega)$ with the same mean and covariance as detailed
above.

###  2.2 LASSO and a Non-Gaussian Model.

Tibshirani (1996) was the first to introduce the LASSO as a method of
constraining the estimates of the linear regression coefficients via an
$L_{1}$ penalty.

$\begin{array}[]{rclcl}\displaystyle\min_{\boldsymbol{\gamma}}&\lx@intercol\|\mathbf{y}-\mathbf{Q}\boldsymbol{\gamma}\|^{2}+\lambda\|\boldsymbol{\gamma}\|_{1}\hfil\lx@intercol\end{array}$
(6)

Here $\mathbf{Q}$ denotes the usual form for covariates in a linear model
where $E(\mathbf{y})=\mathbf{Q}\boldsymbol{\gamma}$.

The main frequentist interest in LASSO is due to its property that the
parameter estimates can be identically zero where the number of zeroes are
controlled by the size of the penalty parameter $\lambda$. In our application
however, covariate selection is not as important as the effect of down
weighting many of the components of $\boldsymbol{\gamma}$ and leaving a sparse
number with larger values.

In parallel with the cubic spline development, let $\boldsymbol{\psi}_{k}(t)$
be a basis and expand $g$ as

$g(t)=\sum_{k=1}^{n}\boldsymbol{\psi}_{k}(t)\boldsymbol{\gamma}_{k}$

To model steps in $g$, we take the basis to be

$\boldsymbol{\psi}_{k}(t)=\left\\{\begin{array}[]{ll}0&\mbox{ for }t<t_{k}\\\
1&\mbox{ for }t\geq t_{k}\\\ \end{array}\right.$ (7)

Let $\Psi_{i,k}=\boldsymbol{\psi}_{k}(t_{i})$ and substitute this matrix for
$\mathbf{Q}$ in the generic LASSO formulation above. From a frequentist
perspective we use the $L_{1}$ penalty to seek a model with a limited number
of steps in $g$. Although this simple device works well for the water
filtering application one could also extend this to more complicated anomalies
by modifying the basis functions. More about this flexibility will be
discussed in Section 5.

Tibshirani at the time also suggested that the LASSO could be interpreted as a
double-exponential (Laplace) prior on $\boldsymbol{\beta}$ in the Bayesian
framework. Figueiredo (2003) expanded on this work and was one of the first to
publish a Bayesian hierarchical model using Laplace priors. Park and Casella
(2008) followed proposing a Bayesian LASSO with a full Gibbs sampler. To our
knowledge however, this model has not been used in a curve fitting context
where it is combined with an additional smooth component.

###  2.3 Combined $L_{2}$ and $L_{1}$ Hybrid Curve Estimator

Based on the preceding development we consider the combined model as

$\mathbf{y}=\Phi\mathbf{c}+\Psi\boldsymbol{\gamma}+\epsilon$

and the combined optimization for these coefficients as

$\begin{array}[]{rclcl}\displaystyle\min_{\boldsymbol{\gamma},\mathbf{c}}&\lx@intercol\|y-(\Phi\mathbf{c}+\Psi\boldsymbol{\gamma})\|^{2}+\lambda\|\boldsymbol{\gamma}\|_{1}+\omega\mathbf{c}^{T}R\mathbf{c}\hfil\lx@intercol\end{array}$
(8)

where $R$ is the cubic smoothing spline penalty matrix discussed in section
2.1. Previous work by Shen et al. (2014) proposed a similar use of LASSO,
modeling the trend with a linear function. Here however, we allow the unknown
trend to be modeled by any smooth function. This new expanded model is not a
particularly easy optimization problem to solve as stated. However, one can
take advantage of the fact that minimization over $\mathbf{c}$ in the
quadratic terms has an analytic solution. Substituting this result back into
(8) reduces the program to a straight $L_{1}$ optimization over only
$\boldsymbol{\gamma}$. Accordingly, the minimization with respect to
$\mathbf{c}$ is

$\hat{\mathbf{c}}=(\Phi^{T}\Phi+\omega
R)^{-1}\Phi^{T}(\mathbf{y}-\Psi\boldsymbol{\gamma})$

Setting $S(\omega)=(\Phi^{T}\Phi+\omega R)^{-1}$ and
$\mathbf{u}=(\mathbf{y}-\Psi\boldsymbol{\gamma})$ gives
$\hat{c}=S\Phi^{T}\mathbf{u}$ and we have the reduced problem in
$\boldsymbol{\gamma}$ as

$\min_{\boldsymbol{\gamma}}{l}{\|\mathbf{u}-\Phi
S(\omega)\Phi^{T}\mathbf{u}\|^{2}_{2}+\lambda\|\boldsymbol{\gamma}\|_{1}+\omega(\mathbf{u})^{T}\Phi
S(\omega)^{T}RS(\omega)\Phi^{T}\mathbf{u}}$

Some straight forward algebra gives a surprisingly simple reduction

$\min_{\boldsymbol{\gamma}}{(\mathbf{u})^{T}(I-\Phi
S(\omega)\Phi^{T})\mathbf{u}+\lambda\|\boldsymbol{\gamma}\|_{1}}$

Finally, if we let $\mathbf{W}(\omega)$ be the symmetric square root matrix of
$I-\Phi S(\omega)\Phi^{T}$, we arrive at a conventional LASSO problem.

$\min_{\boldsymbol{\gamma}}{\|\mathbf{W}(\omega)\mathbf{y}-\mathbf{W}(\omega)\Psi\boldsymbol{\gamma}\|_{2}^{2}+\lambda\|\boldsymbol{\gamma}\|_{1}}$
(9)

Given the solution to (9) one then uses the relationship with
$\hat{\mathbf{c}}$ given above to find the spline component. An important
feature in this derivation is of course the presence of the penalty parameters
$\omega$ and $\lambda$ throughout, and for this reason, we have carried the
dependence on the spline parameter, $\omega$, forward to this last expression.
The main statistical challenge is to estimate these penalty parameters based
on the data and the computational challenge is to recompute the solution to
the optimization problem for each new combination of $(\omega,\lambda)$. These
aspects are addressed in the following sections.

###  2.4 Solving the LASSO

The optimization problem in (9) can be solved using a standard LASSO
implementation after transforming $\mathbf{y}$ and $\Psi$ by multiplication by
$\mathbf{W}(\omega)$. We adopt the FISTA algorithm proposed by Beck and
Teboulle (2009).

Setting $\mathbf{y}^{*}=\mathbf{W}(\omega)\mathbf{y}$ and
$\mathbf{X}=\mathbf{W}(\omega)\Psi$. The steps of the FISTA algorithm are as
follows:

  * •

Step 0: Choose an initial value for $\boldsymbol{\gamma}$. $s_{0}=1$ and
$q_{0}=\boldsymbol{\gamma}_{0}$.

  * •

Step 1: Gradient decent on $\boldsymbol{\gamma}$

$\mathbf{z}_{k}=\boldsymbol{\gamma}_{k}-\tau\mathbf{X}^{T}(\mathbf{X}\boldsymbol{\gamma}_{k}-\mathbf{y}^{*})$

where $\tau$ is the reciprocal of the largest eigenvalue of
$\mathbf{X}^{T}\mathbf{X}$. See Appendix C.

  * •

Step 2: Soft thresholding on $\mathbf{z}_{k}$

$q_{k}=(|\mathbf{z}_{k}|-\alpha)_{+}\text{sign}(\mathbf{z}_{k})$

where $\alpha=\frac{\lambda\tau}{2}$

  * •

Step 3: Update $s_{k}$.

$s_{k+1}=\frac{1+\sqrt{1+4s_{k}^{2}}}{2}$

  * •

Step 4: Update $\boldsymbol{\gamma}$.

$\boldsymbol{\gamma}_{k+1}=q_{k}+\frac{s_{k}-1}{s_{k+1}}(q_{k}-q_{k-1})$

  * •

Step 5: Repeat steps 1 through 4 until convergence.

###  2.5 Choosing Values of $\lambda$ and $\omega$

It is necessary to choose values of $\lambda$ and $\omega$ that balance the
signals and noise and also distinguish between $f$ and $g$. The goal is to
determine values that simultaneously result in a small RMSE, a sparse
$\boldsymbol{\gamma}$, and a smooth trend. In the case of a smoothing spline
or GP, this is typically done by cross validation or maximum likelihood and in
the case of the LASSO penalty one uses an information criterion or cross
validation. The combination of the LASSO and smoothing spline also suggests
using cross validation (CV). We are hesitant however about this approach. In
cross validation an estimate for the removed point is made at the omitted
locations. For smooth functions this is straight forward, however in the
presence of an anomaly, a removed point at or adjacent to a jump is
problematic. For this reason, we considered two methods in choosing $\lambda$
and $\omega$ distinct from CV. The first is a modified machine learning
approach and the second is a version of the Akaike’s Information Criterion
(AIC) (Bozdogan (1987)). The machine learning algorithm does an excellent job
of recovering the smooth trend, but as the reader will see in Section 4, it
lacks the sensitivity to identify small disturbances. The AIC approach
however, while very sensitive to anomalies, tends to choose less parsimonious
models resulting in a poorer estimate of the smooth trend.

####  2.5.1 The Elbow - A Machine Learning Algorithm

We modify an approach based on a machine learning method by Satopaa et al.
(2011) that provides a good choice of $\lambda$ and $\omega$ with the aim of
simultaneously achieving a small RMSE, a parsimonious signal and a small
number of non-zero $\boldsymbol{\gamma}$. For any specific $\lambda$ and
$\omega$ define the root mean squared error (RMSE)

$E(\lambda,\omega)=\sqrt{(1/n)\sum_{i=1}^{n}(y_{i}-\hat{f}(t_{i})-\hat{g}(t_{i}))^{2}}$

of the fit to the data using the hybrid model and also let $E_{0}$ denote the
the RMSE for the base model with just a linear trend and fit by least squares.
Let $H_{2}(\lambda)$ be the number of nonzero elements in
$\hat{\boldsymbol{\gamma}}$ and $N(\lambda,\omega)$ denote the effective
degrees of freedom for a given pair of penalty parameters defined as

$N(\lambda,\omega)=H_{1}(\omega)+H_{2}(\lambda)$

with $H_{1}(\omega)=tr(\Phi(\Phi^{T}\Phi+\omega R)^{-1}\Phi)$. For
justification, recall

$\displaystyle\hat{\mathbf{y}}$
$\displaystyle=\Phi\hat{\mathbf{c}}+\Psi\hat{\boldsymbol{\gamma}}$
$\displaystyle=\Phi(\Phi^{T}\Phi+\omega
R)^{-1}\Phi^{T}(\mathbf{y}-\Psi\hat{\boldsymbol{\gamma}})+\Psi\hat{\boldsymbol{\gamma}}$
$\displaystyle=\Phi(\Phi^{T}\Phi+\omega
R)^{-1}\Phi^{T}\mathbf{y}+(I-\Phi(\Phi^{T}\Phi+\omega
R)^{-1}\Phi^{T})\Psi\hat{\boldsymbol{\gamma}}$ (10)
$\displaystyle=H_{1}\mathbf{y}+H_{2}\hat{\boldsymbol{\gamma}}$

where $H_{2}(\lambda)=(I-\Phi(\Phi^{T}\Phi+\omega R)^{-1}\Phi^{T})\Psi$.

The effective degrees of the freedom in a model with a smoothing penalty is
the rank of the hat matrix which can be estimated to be the sum of the its
diagonal (Tibshirani and Hastie (1987)). In this model however, we have the
hat matrix, $H_{1}$, and a second matrix, $H_{2}$, we need to consider. It is
clear from Equation 2.5.1 that $H_{1}$ and $H_{2}$ are orthogonal to one
another, and therefore it is reasonable to estimate the effective degrees in
the model as the sum of effective degrees of freedom in $H_{1}$ and the
degrees of freedom associated with $H_{2}$. Zou et al. (2007) proposed using
the number of active coefficients in a LASSO model as an unbiased estimate of
the degrees of freedom in $H_{2}$.

Because $N(\lambda,\omega)$ represents the degrees of freedom in the model, we
omit any choices of $\lambda$ and $\omega$ that result in
$N(\lambda,\omega)>n$. Now consider the three dimensional coordinates

$\mathbf{v}(\lambda,\omega)=(N(\lambda,\omega)/n,S(\lambda)/n,E(\lambda,\omega)/E_{0}).$

By this normalization, $0\leq N/n\leq 1$ and $0\leq S/n\leq 1$. the base
spline model also includes a linear trend. Thus the RMSE will be decreasing as
a function of $N$ and $S$, and $E(\lambda,\omega)/E_{0}$ will also be bounded
between 0 and 1. Taken together $\mathbf{v}$ will be under the standard
$2-simplex$ in the positive octant of $\mathcal{R}^{3}$. Estimates for
$\lambda$ and $\omega$ are now determined by searching
$\mathbf{v}(\lambda,\omega)$ over a grid of values and finding the point that
is farthest from the simplex plane. The intuition behind this choice is that
it makes the RMSE small while keeping both the overall model complexity ($N$)
and the sparsity of $\boldsymbol{\gamma}$ ($S$) small. The distance to the
simplex plane for a given $(\lambda,\omega)$ is

$\text{distance}(\lambda,\omega)\propto
1-\sum_{k=1}^{3}\mathbf{v}(\lambda,\omega)_{k}$

Typically the “elbow” point in a function for RMSE is identified as a good
choice for a penalty parameter and is found quantitatively as the maximum in
the second derivative of the RMSE curve. However, the sparsity of
$\boldsymbol{\gamma}$ results in $S(\lambda)$ not being twice differentiable,
so instead the elbow point is inferred by the farthest distance from the
simplex plane.

####  2.5.2 Akaike Information Criterion

Here we adapt a version of AIC that arises from the normality assumption. We
give each combination of $\lambda$ and $\omega$ a score using

$AICc(p)=log(\hat{\sigma}^{2}_{p})+\frac{n+p}{n-p}$

where $\hat{\sigma}^{2}_{k}=\frac{SSE(k)}{n}$ and $p=N(\lambda,\omega)$. We
implement this criterion using the same estimate for effective degrees of
freedom used with the elbow approach.

###  2.6 Identifying Anomalous Intervals

Figure 4 shows the results of estimating $\lambda$ and $\omega$ and the
corresponding signals for the third cycle in Figure 1 using the elbow method.
The ordinary path of the cycle is well modeled by the cubic smoothing spline
and the anomalous portion of the cycle is seen in $g(t)$. Change points are
found for this frequentist approach as those $\boldsymbol{\gamma}$ with values
that are non-zero.

Figure 4: The decomposition of the third cycle in figure 1. The model found
the time indices $73,74$ and $102,103$ to be change points based on the elbow
method for determining $\lambda$ and $\omega$.

These frequentist (hybrid) approaches to identifying change points have proven
very effective. However due to the difficulty of quantifying uncertainty, we
develop an alternative, fully Bayesian approach were $\lambda$ and $\omega$
can be identified as parameters in a hierarchical model. This switch to a
Bayesian framework will also provide error estimates on the model coefficients
and thus quantify the uncertainty in the detection process.

##  3 A Bayesian Approach to Identifying Anomalies

We now propose a Bayesian hierchical model (BHM) for identifying disturbances
where the penalties in the hybrid approach are interpreted as the log of prior
densities while the mode of the log posterior distribution is related to the
solution to the optimization problem in Section 2.3. In this way the two can
be compared.

###  3.1 A Hierarchical Model

We propose the top level to be the same observational model but now assume a
Laplace prior for $\boldsymbol{\gamma}$ and a multivariate normal prior on
$\mathbf{c}$.

$\displaystyle y|\mathbf{c},\boldsymbol{\gamma},\sigma^{2}$ $\displaystyle\sim
N_{n}(\Phi\mathbf{c}+\Psi\boldsymbol{\gamma},\sigma^{2}I_{n})$
$\displaystyle[\boldsymbol{\gamma}|\sigma^{2},\lambda^{2}]$
$\displaystyle=\left(\frac{\lambda^{2}}{2\sqrt{\sigma^{2}}}\right)^{n-1}e^{-\frac{\lambda^{2}}{\sqrt{\sigma^{2}}}||\boldsymbol{\gamma}||_{1}}=\prod_{j=1}^{n-1}\frac{\lambda^{2}}{2\sqrt{\sigma^{2}}}e^{-\frac{\lambda^{2}}{\sqrt{\sigma^{2}}}|\boldsymbol{\gamma}_{j}|}$
$\displaystyle[\mathbf{c}|\omega]$
$\displaystyle\propto\omega^{1/2}exp(-\frac{\omega}{2}\mathbf{c}^{T}R\mathbf{c})$
(11) $\displaystyle[\sigma^{2}]$ $\displaystyle\propto 1$
$\displaystyle\lambda^{2}$ $\displaystyle\sim
gamma(\alpha_{\lambda^{2}},\beta_{\lambda^{2}})$ $\displaystyle\omega$
$\displaystyle\sim gamma(\alpha_{\omega},\beta_{\omega})$

Notice that the prior on $\boldsymbol{\gamma}$ is an uncorrelated multivariate
$Laplace(0,\frac{\sqrt{\sigma^{2}}}{\lambda^{2}})$. Additionally, the prior on
$\mathbf{c}$ is a multivariate Gaussian distribution, however $R$ is not full
rank. The rank deficiency of $R$ poses a problem that needs to be considered.
One way to address this problem is to utilize the elegant GP representation of
the smoothing cubic spline discussed in Section 2.1. The smooth trend
$\Phi\mathbf{c}$ becomes $X\boldsymbol{\beta}+\mathbf{g}$ with priors
$[\boldsymbol{\beta}]\propto 1$ and $\mathbf{g}\sim N(0,K)$.

These priors work to achieve sparseness in $\boldsymbol{\gamma}$ and
smoothness for the trend. Finally, we use an improper prior on
$\boldsymbol{\beta}$ to be consistent with the BLUE frequentist estimator.
Following the terminology in spatial statistics we refer to
$\boldsymbol{\beta}$ as a fixed effect even though it is a random quantity in
the BHM.

####  3.1.1 Reformulating the $\boldsymbol{\gamma}$ parameter

Unfortunately, the prior on $\boldsymbol{\gamma}$ is not conjugate, and if
this model is implemented in its current form, sampling $\boldsymbol{\gamma}$
would require a Metropolis-Hastings step. Our goal is to formulate this model
so that all sampling involves Gibbs updates. By adding an additional level to
this hierarchy with some extra parameters, it is possible to avoid a
Metropolis-Hastings update on $\boldsymbol{\gamma}$, thereby achieving a Gibbs
sampler but functionally maintaining the Laplace prior. Here we use the
insightful result of Andrews and Mallows (1974) and utilized by Park and
Casella (2008) where a Laplace distribution can be written as an exponential
scale mixture of Gaussian distributions.

$\frac{a}{2}e^{-a|z|}=\int_{0}^{\infty}\frac{1}{\sqrt{2\pi
s}}e^{-z^{2}/(2s)}\frac{a^{2}}{2}e^{-a^{2}s/2}ds,\;\;\;a>0$ (12)

Using this representation and the GP version of the cubic smoothing spline
results in an equivalent hierarchical model

$\displaystyle y|\boldsymbol{\beta},\mathbf{g},\boldsymbol{\gamma},\sigma^{2}$
$\displaystyle\sim
N_{n}(X\boldsymbol{\beta}+\Psi\boldsymbol{\gamma}+\mathbf{g},\sigma^{2}I_{n})$
$\displaystyle\boldsymbol{\gamma}|\sigma^{2},\tau_{1}^{2},...,\tau_{n-1}^{2}$
$\displaystyle\sim N_{n-1}(0_{n-1},\sigma^{2}F)\text{ where
}F=diag(\tau_{1}^{2},...,\tau_{n-1}^{2})$
$\displaystyle\tau_{i}^{2}|\lambda^{2}$ $\displaystyle\sim
exp(\frac{\lambda^{2}}{2})\text{ for $i=1,...,n-1$}$
$\displaystyle\mathbf{g}|\omega$ $\displaystyle\sim
N_{n}(0_{n},\frac{1}{\omega}K)$ $\displaystyle[\sigma^{2}]$
$\displaystyle\propto 1$ $\displaystyle\lambda^{2}$ $\displaystyle\sim
gamma(\alpha_{\lambda^{2}},\beta_{\lambda^{2}})$ $\displaystyle\omega$
$\displaystyle\sim gamma(\alpha_{\omega},\beta_{\omega})$
$\displaystyle[{\boldsymbol{\beta}}]$ $\displaystyle\propto 1$

Here the second and third levels give the same conditional distribution for
$\boldsymbol{\gamma}$ as the second line of the first hierarchy, and a proof
of this fact is given in Appendix D.

####  3.1.2 Reducing Confounding of Fixed and Random effects

The BHM in $3.1.1$, however, suffers from an additional problem. It does not
mix well due to confounding between the fixed and random effects (Hodges and
Reich (2010)). This problem is a well known feature in these types of BHMs and
is distinct from the frequentist/optimization where modfications are not
needed. The introduction of the Gaussian process into this model is known to
result in spatial confounding with the fixed effects in $X$. Not accounting
for this confounding results in poor mixing of $\boldsymbol{\beta}$,
$\mathbf{g}$, and $\boldsymbol{\gamma}$ and a general failure for the MCMC
chains to converge. The process of orthogonalizing the fixed effects with the
random effects and then (nearly) orthogonalizing the Gaussian process,
$\mathbf{g}$ with the discontinuous stochastic component
$\Psi\boldsymbol{\gamma}$, is accomplished with projection matrices based on
$\mathbf{X}$ and $\Psi$. In the process, the key elements
$\boldsymbol{\beta}$, $\mathbf{g}$, and $\boldsymbol{\gamma}$ are
reparametrized in a way that improves mixing and convergence. A novel element
in this implementation is the use of a regularized version of the projection
for the anomaly basis function to make the computations numerically stable.

Let $\mathbf{P}_{X}=X(X^{T}X)^{-1}X^{T}$ . Then

$\displaystyle\mathbf{y}$
$\displaystyle=X\boldsymbol{\beta}+\Psi\boldsymbol{\gamma}+\mathbf{g}+\epsilon$
$\displaystyle=X\boldsymbol{\beta}+P_{\mathbf{X}}(\Psi\boldsymbol{\gamma}+\mathbf{g})+(I-P_{\mathbf{X}})(\Psi\boldsymbol{\gamma}+\mathbf{g})+\epsilon$
$\displaystyle=X\boldsymbol{\beta}^{*}+\Psi^{*}\boldsymbol{\gamma}+(I-P_{\mathbf{X}})\mathbf{g}+\epsilon$

with the reparamterization
$\boldsymbol{\beta}^{*}=\boldsymbol{\beta}+(X^{T}X)^{-1}X^{T}(\Psi\boldsymbol{\gamma}+\mathbf{g})$
and $\Psi^{*}=(I-P_{X})\Psi$, a projection of the anomaly basis functions onto
the space orthogonal to the linear trend.

We now use the same operation to further improve sampling, letting

$P_{\Psi^{*}}=\Psi^{*}({\Psi^{*}}^{T}\Psi^{*}+I\delta)^{-1}{\Psi^{*}}^{T}$

where $\delta I$ functions as a small regularization term, on the order of
single precision zero (e.g. $\delta=10^{-8}$) because $\Psi^{*}$ is
computationally singular.

$\displaystyle\mathbf{y}$
$\displaystyle=X\boldsymbol{\beta}^{*}+\Psi^{*}\boldsymbol{\gamma}+(I-P_{X})\mathbf{g}+\epsilon$
$\displaystyle=X\boldsymbol{\beta}^{*}+\left(\Psi^{*}\boldsymbol{\gamma}+P_{\Psi^{*}}(I-P_{\mathbf{X}})\mathbf{g}\right)+(I-P_{\Psi^{*}})(I-P_{\mathbf{X}})\mathbf{g}$
$\displaystyle=X\boldsymbol{\beta}^{*}+\Psi^{*}\boldsymbol{\gamma}^{*}+H\mathbf{g}$
(13)

with the final reparametrizations

$\boldsymbol{\gamma}^{*}=\boldsymbol{\gamma}+J\mathbf{g}$,
$J=({\Psi^{*}}^{T}\Psi^{*}+I\delta)^{-1}{\Psi^{*}}^{T}(I-P_{\mathbf{X}})$, and
$H=(I-P_{\Psi^{*}})(I-P_{X})$

Summarizing these modifications we arrive at an “orthogonalized” BHM (OBHM)
that is equivalent to the first model proposed, up to the approximation for
$P_{\Psi^{*}}$, but it is much more efficient to apply Gibbs sampling and a
Markov Chain Monte Carlo (MCMC) algorithm.

$\displaystyle
y|\boldsymbol{\beta}^{*},\mathbf{g},\boldsymbol{\gamma}^{*},\sigma^{2}$
$\displaystyle\sim
N_{n}(X\boldsymbol{\beta}^{*}+\Psi^{*}\boldsymbol{\gamma}^{*}+H\mathbf{g},\sigma^{2}I_{n})$
$\displaystyle\boldsymbol{\gamma}^{*}|\sigma^{2},\tau_{1}^{2},...,\tau_{n-1}^{2},\mathbf{g}$
$\displaystyle\sim N_{n-1}(J\mathbf{g},\sigma^{2}F)\text{ where
}F=diag(\tau_{1}^{2},...,\tau_{n-1}^{2})$
$\displaystyle\tau_{i}^{2}|\lambda^{2}$ $\displaystyle\sim
exp(\frac{\lambda^{2}}{2})\text{ for $i=1,...,n-1$}$
$\displaystyle\mathbf{g}|\omega$ $\displaystyle\sim N(0,\frac{1}{\omega}K)$
$\displaystyle\omega$
$\displaystyle\sim\text{gamma}(\alpha_{\omega},\beta_{\omega})$
$\displaystyle\lambda^{2}$
$\displaystyle\sim\text{gamma}(\alpha_{\lambda^{2}},\beta_{\lambda^{2}})$
$\displaystyle[\boldsymbol{\beta}^{*}]$ $\displaystyle\propto 1$
$\displaystyle[\sigma^{2}]$ $\displaystyle\propto 1$

###  3.2 Gibbs Updates and Posterior Distributions

All of the full conditional distributions for the OBHM are straightforward to
work out with the possible exception of $\tau_{j}^{2}$ which is derived out in
Appendix E. Below is a list of the full conditional posterior distributions of
each level. Note that in the posteriors of $\lambda^{2}$ and $\omega$, the
second parameter is the rate.

$\displaystyle\boldsymbol{\beta}^{*}|\cdot$ $\displaystyle\sim
N_{n}(A^{-1}b,A^{-1})\text{ where }A=\frac{X^{T}X}{\sigma^{2}}\text{ and
}b=\frac{X^{T}(Y-\Psi^{*}\boldsymbol{\gamma}^{*}-H\mathbf{g})}{\sigma^{2}}$
(14) $\displaystyle\boldsymbol{\gamma}^{*}|\cdot$ $\displaystyle\sim
N_{n}(A^{-1}b,A^{-1})\text{ where
}A=\frac{\Psi^{*T}\Psi^{*}+F^{-1}}{\sigma^{2}}\text{ and
}b=\frac{\Psi^{*T}(y-X\boldsymbol{\beta}^{*}-H\mathbf{g})+F^{-1}J\mathbf{g}}{\sigma^{2}}$
(15) $\displaystyle\sigma^{2}|\cdot$ $\displaystyle\sim\text{inv-
gamma}\left(\frac{2n-1}{2},r\right)$ (16) $\displaystyle\text{where }r$
$\displaystyle=\frac{(y-X\boldsymbol{\beta}^{*}-\Psi^{*}\boldsymbol{\gamma}^{*}-H\mathbf{g})^{T}(y-X\boldsymbol{\beta}^{*}-\Psi^{*}\boldsymbol{\gamma}^{*}-H\mathbf{g})}{2}+\frac{(\boldsymbol{\gamma}^{*}-J\mathbf{g})^{T}F^{-1}(\boldsymbol{\gamma}^{*}-J\mathbf{g})}{2}$
$\displaystyle\lambda^{2}|\cdot$ $\displaystyle\sim\text{gamma
}(\alpha_{\lambda^{2}}+n-1,\beta_{\lambda^{2}}+\sum_{j=1}^{n}\tau_{j}^{2}/2)$
(17) $\displaystyle\frac{1}{\tau_{j}^{2}}|\cdot$ $\displaystyle\sim\text{inv-
Gaussian}\left(\sqrt{\frac{\sigma^{2}\lambda^{2}}{(\boldsymbol{\gamma}^{*}-J\mathbf{g})_{j}^{2}}},\lambda^{2}\right)$
(18) $\displaystyle\mathbf{g}|\cdot$ $\displaystyle\sim
N_{n}(A^{-1}b,A^{-1})\text{ where }A=\omega
K^{-1}+\frac{H^{T}H+J^{T}F^{-1}J}{\sigma^{2}}$ (19) $\displaystyle\text{and
}b=\frac{H^{T}(Y-X\boldsymbol{\beta}^{*}-\Psi^{*}\boldsymbol{\gamma}^{*})+J^{T}F^{-1}\boldsymbol{\gamma}^{*}}{\sigma^{2}}$
$\displaystyle\omega$
$\displaystyle\sim\text{gamma}\left(\frac{n}{2}+\alpha_{\omega},\frac{\mathbf{g}^{T}K^{-1}\mathbf{g}}{2}+\beta_{\omega}\right)$
(20)

Upon completion of the MCMC it is necessary to parametrize back to
$\boldsymbol{\beta}$ and $\boldsymbol{\gamma}$. This is done by finding
$\boldsymbol{\gamma}=\boldsymbol{\gamma}^{*}-Jf$ and
$\boldsymbol{\beta}=\boldsymbol{\beta}^{*}-(X^{T}X)^{-1}X^{T}(\Psi\boldsymbol{\gamma}+\mathbf{g})$.
Although the draws were made of $\boldsymbol{\gamma}^{*}$ and
$\boldsymbol{\beta}^{*}$ the inferences can still be made on
$\boldsymbol{\gamma}$ and $\boldsymbol{\beta}$.

###  3.3 Convergence, Priors, and Initial Values

We implemented a uniform prior on $\boldsymbol{\beta}$ so that the combination
of the linear trend identified by $X\boldsymbol{\beta}$ and the GP,
$\mathbf{g}$, is precisely a cubic smoothing spline. This equality was
discussed in Section 2.1.

A uniform prior on $\sigma^{2}$ was chosen for the sake of simplicity. We note
however that although we achieved convergence of the MCMC with this prior, we
did notice a higher effective sampling rate when using a loosely informative
inverse gamma prior on $\sigma^{2}$. We decided on loosely informative priors
on $\lambda^{2}$ and $\omega$ making use of the analysis from the hybrid
approach. These loosely informative priors are relatively flat yet place
greater likelihood near the mean of each parameter gathered from the
frequentist approach. Improper, uniform priors were tested for $\lambda^{2}$
and $\omega$ and rejected because they greatly reduced the effective sample
sizes in each chain to the point that the MCMC struggled to converge.

The initial values of $\boldsymbol{\gamma}$ were chosen by finding
$z_{i}=|y_{i+1}-y_{i}|$ for $i=2,..,n$ and the $95^{th}$ percentile of these
values, $z_{.95}$. We then defined $\boldsymbol{\gamma}$ as

$\boldsymbol{\gamma}_{i}=\left\\{\begin{array}[]{ll}0&\mbox{ for
}z_{i}<z_{.95}\\\ z_{i}&\mbox{ for }z_{i}\geq z_{.95}\\\ \end{array}\right.$
(21)

Given the initial guess of the disturbances, $\boldsymbol{\gamma}$, natural
starting values for $\boldsymbol{\beta}$ are the OLS coefficient estimates
regressing $y-\Psi\boldsymbol{\gamma}$ on $X$. We finally set the initial
values of $\lambda^{2}$, $\omega$, and all $\tau^{2}_{i}$ to the mean values
of their distributions and sampled $\sigma^{2}$ first in our MCMC. Convergence
and mixing happen quickly in this formulation requiring a burn-in of less than
$200$ iterations and chains of only $1000$ dependent samples.

###  3.4 Identifying Disturbance Locations

Figure 5 shows $50$ posterior predictions from the MCMC modelling the third
cycle in Figure 1, the smooth portion of the model on the left and the rough
portion on the right. Again, we see a clear separation between the two. An
alternative form of $\Psi$ discussed in Appendix B was used to reduce the size
of the confidence intervals on the edges. Figure 6 shows the $95\%$ confidence
intervals of each $\boldsymbol{\gamma}$ coefficient. These posterior results
were then used in combination with a threshold of $.15$ units headloss for
identifying components of $\boldsymbol{\gamma}$ as disturbance locations. This
threshold was originally chosen for judging the methods because smaller
anomalies were not discernible by visual inspection.

Figure 5: The two plots are the decomposition of the third cycle in Figure 1
using the Bayesian GP approach. The first plot contains posterior predictions
of $X\boldsymbol{\beta}+\mathbf{g},$ and the second contains posterior
predictions of $\Phi\boldsymbol{\gamma}$. Figure 6: Displayed are the $95\%$
credible intervals of each $\boldsymbol{\gamma}$ coefficient. The dotted lines
are placed at $\pm.15$. Any locations at which the credible interval does not
include zero and the mean is outside of the dotted lines is flagged as a
relevant change point. This approach, similar to the hybrid approach, found
the time indices $73,74$ and $102,103$ to be change points.

##  4 Simulation Study

In order to test the effectiveness of the proposed methods, we ran a
simulation study using the trend found in Figure 3 and adding a disturbance at
time index $150$. We ran the study over a grid of $\sigma$ and disturbance
sizes, producing $100$ datasets for each combination of noise and disturbance.
We allowed the disturbance direction to vary randomly. For each dataset, we
recorded whether the disturbance was detected using the criteria that a non-
zero $\boldsymbol{\gamma}$ value is detected in the hybrid approach or a
central $95\%$ credible interval for $\boldsymbol{\gamma}_{k}$ does not
include $0$ in the Bayesian approach. The estimated probability of detection
surface can be seen in Figure 7 along with contour lines at $p=.1$, $.5$, and
$.9$.

The results of the Bayesian and hybrid AICc approaches are nearly identical
while the Hybrid approach utilizing the elbow method for the selection of
$\lambda$ and $\omega$ lacks the sensitivity needed to identify smaller
disturbances as can be seen in the lower left of the plots in Figure 7. For
disturbances greater than $.5$ all three methods are essentially the same.
Based on the results, there appears to be a lower limit of detection for the
hybrid elbow approach that does not exist for other two approaches.

Figure 7: The results of the simulation study displaying the probability of
detection under given disturbance and noise sizes with contours at $p=.1$
(lowest line), $.5$ (middle line), and $.9$ (highest line).

##  5 Analysis of Denver Water Dataset

Denver Water provided one year of headloss data on sixteen filters, measured
at five minute intervals, for analysis. The lengths of the filtration cycles
ranged from 30 minutes to 48 hours. A number of factors appear to influence
cycle lengths including: water quality, weather, season, operator preferences,
and maintenance events. Operators typically ended a cycle and initiated a
cleaning (a backwash) when headloss reached 8-8.5 units, but on numerous
occasions cycles were ended after only a short time, long before reaching the
8-8.5 range. The net result is that these data are quite heterogenous,
differing in cycle length and across the different filters, and so makes the
detection of anomalous cycles more difficult.

The two methods were applied to $3225$ cycles of which $64$ anomalous cycles
were identified by the facility environmental engineers. The Bayesian approach
used four MCMC chains of $1000$ with a burn-in of $200$ and required that
disturbances be at least .15 units headloss. The engineers had not deemed
smaller disturbances significant. We attached a stopping limit on the FISTA
algorithm in the hybrid approach of $1e-4$ and required that disturbances be
greater than .15 units headloss, although this second requirement was
unnecessary when using the elbow method of selection $\lambda$ and $\omega$
due to the less sensitive nature of this method.

The Bayesian, the hybrid AICc, and the hybrid elbow approaches were successful
in identifying nearly all of the anomalous cycles. Given that we selected to
only flag disturbances greater than $.15$, it is not surprising that the elbow
and AICc hybrid approaches yielded nearly identical results. Two of the three
missed cycles by the Bayesian and hybrid elbow approaches can be seen in
Figure 8. All three cycles share a common feature that the anomaly is not an
abrupt jump in the mean, so it is understandable why these cycles were missed
as they are not what our model expects.

Figure 8: Examples of two false positives.

Before setting a minimum size on disturbance identification, both the hybrid
AIC and BHM approaches yielded a large number of false positives, roughly
$400$ of the $3141$ normal cycles. Upon closer examination however, the models
had identified entire months of anomalous data which had gone previously
unrecognized by Denver Water staff. Many of these newly identified cycles have
a subtle, sawtooth nature to their measurements that is only barely visible.
Although small, this structure signifies less than ideal operation conditions
and should be addressed (Figure 9). This level of sensitivity speaks to the
effectiveness of these methods, and also suggests the need to properly
identify what level of disturbance signifies problematic behavior.

Figure 9: Cycle 3009, one of many cycles identified only by the Bayesian
method as anomalous and originally thought to be a false positive. Upon
zooming in however, one can see the irregular, sawtooth nature of the
measurements.

When running on four cores on a MacBook Pro (13-inch, M1, 2020), the
computation time required for analyzing cycles is seen in Figure 10. Overall,
the hybrid approach is roughly twenty times faster than the BHM. Beside length
of cycle, the main factor in computation time for the hybrid approach is the
grid size of $\lambda$ and $\omega$; for the Bayesian GP approach the main
factor is the length of the MCMC chains.

Hybrid Approach (Elbow/AICc)
IDENTIFIED ACTUAL Anomaly Normal Disturbance in Mean 61/63 3/1 Normal 63/49
3098/3112

Bayesian Approach
IDENTIFIED ACTUAL Anomaly Normal Disturbance in Mean 61 3 Normal 23 3138

Table 1: Confusion matrices for each approach (Section 5) in identifying
anomalous cycles identified by Denver Water facility engineers. Figure 10:
Computation times for anomaly detection. The Hybrid approach is approximately
twenty times faster than the Bayesian approach.

##  6 Conclusions

Recovering a smooth trend from a time series using quadratic roughness
penalties, an $L_{2}$ type criteria, or selecting a subset of variables using
the LASSO, an $L_{1}$ type criteria, are both well established and successful
data analysis tools. Here we have combined these models to detect anomalous
signals in process monitoring. At the outset, it may seem that this is not
possible due to the ambiguity of separating a signal into smooth and anomalous
components, however we have found that, for our motivating application in
water treatment, the method works well showing false negative rates (anomalous
behavior not detected) of a few percent. The false positive rate is also low
on the order of a couple percent when tuned to remove very small disturbances.
Additionally, both the BHM and the hybrid AICc approaches provide an
opportunity to detect more subtle discontinuities in the filtering cycles that
are not easily spotted by a visual inspection of a plot. Besides having
success for a substantial application, this work also provides two
methodological contributions. One is testing in a Monte Carlo setting the
simplex strategy for choosing the two smoothing/sparsity parameters in the
frequentist method. This parameter choice drawn from machine learning worked
surprisingly well and is also fast. Its success motivates some additional
investigation into deeper statistical principles as to why this works. The
other contribution is in transforming the hybrid model into equivalent forms
that make it efficient to estimate and perform an MCMC computation. This
includes reducing the minimization for the hybrid estimator to a linearly
transformed problem that just has the $L_{1}$ minimization. For computing the
BHM we were able to rephrase the model to avoid more expensive MH updates and
also to reparametrize the linear models to accelerate mixing. The net result
is an efficient MCMC sampler for a complex hierarchical model.

The places where our models misidentify filter cycles suggest two areas for
further work. The model is sensitive to step changes in the signal but we
noticed that the three false negative cases are better represented by
discontinuities in the derivative, e.g. a change in slope. This feature could
be included by basis functions and penalties that focus on this additional
kind of anomaly. More generally our approach is flexible in what basis
functions are used for the anomaly signal and we expect that the BHM will
still be efficient if one switches to other covariance functions for the GP
smooth trend and different anomaly signatures. Future work could be applying
this model for the detection of multiple change points in observational
climate records.

One current topic in spatial statistics is models for non-Gaussian fields. We
note that the BHM is not tied to a particular GP and one could add extra
parameters for the covariances typically used in geostatistics. The
$L_{2}$/$L_{1}$ model would make the estimated fields more interpretable given
the separation into a GP and non-Gaussian process components. One strategy for
handling the analysis of large spatial data sets is to expand the spatial
fields in basis functions with a multivariate normal prior on the
coefficients. This is similar to the $L_{2}$ part of our hybrid smoother and
it is reasonable to extend this to additional basis functions to capture
discontinuities. This additional part of the model could capture transitions
between land/water boundaries for climate data or demographic and municipal
regions when one considers social data. In general, this hybrid model has
potential to go beyond representations of data that are not supported by just
a smooth surface and random noise.

##  7 Acknowledgements

The authors gratefully acknowledge Craig McGonagill from Denver Water for
making the filter data available, identifying the anomalous cycles, and and
more generally providing insight into the water water filtration process at
his facility. The authors also acknowledge Mike Wakin and Mevin Hooten for
useful conversations.

##  8 Funding

This research is supported in part by the National Science Foundation under
NSF Award 1924146 HDR DSC: Collaborative Research: Modernizing Water and
Wastewater Treatment through Data Science Education & Research (MoWaTER)

## References

  * Andrews and Mallows [1974] David F Andrews and Colin L Mallows.  Scale mixtures of normal distributions.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , 36(1):99–102, 1974.
  * Beck and Teboulle [2009] Amir Beck and Marc Teboulle.  A fast iterative shrinkage-thresholding algorithm for linear inverse problems.  _SIAM journal on imaging sciences_ , 2(1):183–202, 2009.
  * Bozdogan [1987] Hamparsum Bozdogan.  Model selection and akaike’s information criterion (aic): The general theory and its analytical extensions.  _Psychometrika_ , 52(3):345–370, 1987.
  * Chhikara and Folks [1989] RS Chhikara and JL Folks.  The inverse gaussian distribution: theory, 1989.
  * Figueiredo [2003] Mário AT Figueiredo.  Adaptive sparseness for supervised learning.  _IEEE transactions on pattern analysis and machine intelligence_ , 25(9):1150–1159, 2003.
  * Hodges and Reich [2010] James S Hodges and Brian J Reich.  Adding spatially-correlated errors can mess up the fixed effect you love.  _The American Statistician_ , 64(4):325–334, 2010\.
  * Nychka et al. [2021] Nychka, Reinhard Furrer, John Paige, and Stephan Sain.  fields: Tools for spatial data, 2021.  URL https://github.com/dnychka/fieldsRPackage.  R package version 15.2.
  * Park and Casella [2008] Trevor Park and George Casella.  The bayesian lasso.  _Journal of the American Statistical Association_ , 103(482):681–686, 2008.
  * Satopaa et al. [2011] Ville Satopaa, Jeannie Albrecht, David Irwin, and Barath Raghavan.  Finding a ”kneedle” in a haystack: Detecting knee points in system behavior.  In _2011 31st International Conference on Distributed Computing Systems Workshops_ , pages 166–171, 2011.  doi: 10.1109/ICDCSW.2011.20.
  * Shen et al. [2014] Jie Shen, Colin M Gallagher, and QiQi Lu.  Detection of multiple undocumented change-points using adaptive lasso.  _Journal of Applied Statistics_ , 41(6):1161–1173, 2014.
  * Tibshirani [1996] Robert Tibshirani.  Regression shrinkage and selection via the lasso.  _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , 58(1):267–288, 1996.
  * Tibshirani and Hastie [1987] Robert Tibshirani and Trevor Hastie.  Local likelihood estimation.  _Journal of the American Statistical Association_ , 82(398):559–567, 1987.
  * Wahba [1990] Grace Wahba.  _Spline models for observational data_.  SIAM, 1990.
  * Zou et al. [2007] Hui Zou, Trevor Hastie, and Robert Tibshirani.  On the “degrees of freedom” of the lasso.  2007\.

##  Appendix A Cycle Separation

The headloss data from Denver Water were not reported as individual cycles, so
the first step in the analysis process was separating the cycles. There are a
couple approaches to separating the headloss data into individual cycles that
can be used. Ideally, headloss data is accompanied by a measurement of whether
the filter is on or off which makes the separation process simple. If however
the on/off information is not present, then cycles can be separated by passing
a mathematical filter over the data and analyzing the inner products of the
data and the filter. The one used in this paper can be seen in Figure 11. Once
the inner products of the filter and date are computed at each point in the
cycle, these values can be combined in various ways to identify the start and
end points of the cycles as seen in the figure. The start and end points
identified using this filter approach are a near perfect match of the on/off
data reported by Denver Water, however for a number of cycles the on/off data
was found to be inaccurate and the filter approach did a better job of finding
the correct start and end times.

Figure 11: The three filters passed over the headloss data to identify the
start and end points of each cycle. The dotted lines show the filter scores
used to determine the start and end of each cycle.

Once the cycles were separated, a decision was made to eliminate the warm up
and cool down phases of each cycle and preserve only the running phase because
all of the anomalous events occur in this phase. This further trimming of the
data was achieved by finding the first differences in the data and removing
the sections on the ends of each cycle associated with the largest absolute
differences since the most extreme changes in headloss occur in the warm up
and cool down phases.

##  Appendix B Other Constructions of $\Psi$

The piecewise stochastic portion of the model, $g(t)=\Psi\boldsymbol{\gamma}$,
designed to detect disturbances in the mean holds some intriguing
possibilities if we allow for other formulations of $\Psi$.

In the recovery of $g(t)$, we found that we can dramatically reduce the error
bounds by reformulating $\Psi$ to reconstruct $g(t)$ from the center of the
time series instead of the beginning. $\Psi$ is an upper triangular matrix of
ones on the left half of the matrix and a lower triangular matrix of ones on
the right side of the matrix, and the middle column is removed form the matrix
because the mean is already identified elsewhere in the model. This reduced
the size of the largest credible intervals around any location in $g(t)$ by a
factor of two. This is the construction of $\Psi$ that we used in the
implementation of all three approaches. An important note, however, is that
the choice of $\Psi$ does not change the detection of significant
$\boldsymbol{\gamma}$ coefficients or their distributions.

##  Appendix C The Lipschitz Constant in FISTA

Convergence of the FISTA algorithm requires that
$h:\mathcal{R}^{n}\rightarrow\mathcal{R}$ is convex and differentiable and
that $\nabla h$ is Lipschitz continuous with constant $L\geq 0$, that is that

$\|\nabla h(x)-\nabla h(y)\|\leq L\|x-y\|$

for all $x$ and $y$. $L$ can be taken to its limit, and the Lipschitz constant
becomes the largest eigenvalue of the Hessian matrix of $h$. In both programs
$f=\|\mathbf{W}(\omega)\mathbf{y}-\mathbf{W}(\omega)\boldsymbol{\gamma}\|^{2}_{2}$,
so
$\nabla^{2}h=\Psi^{T}\mathbf{W}(\omega)^{T}\mathbf{W}(\omega)\Psi=2\mathbf{X}^{T}\mathbf{X}$.

##  Appendix D Proof Model Has Functionally Laplace Prior

Demonstrating that this proposed hierarchical model is the functional
equivalent of the double exponential requires integrating out each of the
$\tau_{j}$’s.

$\displaystyle[\boldsymbol{\gamma}|\sigma^{2},\lambda^{2}]$
$\displaystyle\propto\prod_{j=1}^{n}\int_{0}^{\infty}\frac{1}{\sqrt{\sigma^{2}\tau_{j}^{2}}}exp(-\frac{1}{\sigma^{2}\tau_{j}^{2}}\boldsymbol{\gamma}_{j}^{2})\frac{\lambda^{2}}{2}exp(-\frac{\lambda^{2}\tau_{j}^{2}}{2})d\tau_{j}^{2}$
(22)
$\displaystyle\propto\prod_{j=1}^{n}\int_{0}^{\infty}\frac{1}{\sqrt{\sigma^{2}\tau_{j}^{2}}}exp(-\frac{1}{\sigma^{2}\tau_{j}^{2}}\boldsymbol{\gamma}_{j}^{2})\frac{\lambda^{2}}{2\sigma^{2}}exp(-\frac{\lambda^{2}\tau_{j}^{2}\sigma^{2}}{2\sigma^{2}})\sigma^{2}d\tau_{j}^{2}$
(23)

using the change of variables $s=\sigma^{2}\tau_{j}^{2}$,
$a=\sqrt{\frac{\lambda^{2}}{\sigma^{2}}}$, we get

$\displaystyle[\boldsymbol{\gamma}|\sigma^{2},\lambda^{2}]$
$\displaystyle=\prod_{j=1}^{n}\int_{0}^{\infty}\frac{1}{\sqrt{2\pi
s}}e^{-\boldsymbol{\gamma}^{2}/(2s)}\frac{a^{2}}{2}e^{-a^{2}s/2}ds$ (24)

And by Identity (12),

$\displaystyle[\boldsymbol{\gamma}|\sigma^{2},\lambda^{2}]$
$\displaystyle=\frac{a}{2}e^{-a|\boldsymbol{\gamma}|}$ (25)

##  Appendix E Conditional Posterior of $\tau_{j}^{2}$

The posterior of $\frac{1}{\tau_{j}^{2}}$ is the inverse Gaussian distribution
discussed by Chhikara and Folks [1989]. Here is an explanation of the
formulation found in Park and Casella [2008].

$\displaystyle[\tau_{j}^{2}]$
$\displaystyle\propto\frac{1}{|\sigma^{2}F|^{1/2}}exp(\frac{-1}{2\sigma^{2}}\boldsymbol{\gamma}^{T}F^{-1}\boldsymbol{\gamma})\prod_{j=1}^{n}\frac{\lambda^{2}}{2}e^{-\lambda^{2}\tau_{j}^{2}/2}$
(27)
$\displaystyle\propto\prod_{j=1}^{n}\frac{\lambda^{2}}{2(\sigma^{2}\tau_{j}^{2})^{1/2}}exp(-\frac{\lambda^{2}\tau_{j}^{2}}{2}-\frac{1}{2\sigma^{2}}\frac{\boldsymbol{\gamma}_{j}^{2}}{\tau^{2}_{j}})$
(28)

Applying the change of variables $x=\frac{1}{\tau^{2}_{j}}$ results in

$\displaystyle[x]$
$\displaystyle\propto\prod_{j=1}^{n}x^{1/2}exp(-\frac{\boldsymbol{\gamma}_{j}^{2}}{2\sigma^{2}}(x+\frac{\lambda^{2}\sigma^{2}}{\boldsymbol{\gamma}_{j}^{2}x}))x^{-2}$
(29)
$\displaystyle\propto\prod_{j=1}^{n}x^{-3/2}exp(-\frac{\boldsymbol{\gamma}_{j}^{2}\lambda^{2}}{2\sigma^{2}\lambda^{2}x}(x^{2}+\frac{\lambda^{2}\sigma^{2}}{\boldsymbol{\gamma}_{j}^{2}}))$
(30)
$\displaystyle\propto\prod_{j=1}^{n}x^{-3/2}exp(-\frac{\boldsymbol{\gamma}_{j}^{2}\lambda^{2}}{2\sigma^{2}\lambda^{2}x}(x-\sqrt{\frac{\lambda^{2}\sigma^{2}}{\boldsymbol{\gamma}_{j}^{2}}})^{2})$
(31)

Letting
$\mu^{\prime}=\sqrt{\frac{\lambda^{2}\sigma^{2}}{\boldsymbol{\gamma}_{j}^{2}}}$
and $\lambda^{\prime}=\lambda^{2}$, we get

$\displaystyle[x]$
$\displaystyle\propto\prod_{j=1}^{n}x^{-3/2}exp(-\frac{\lambda^{\prime}}{2\mu^{\prime}x}(x-\mu^{\prime})^{2})$
(32) $\displaystyle\frac{1}{\tau_{j}^{2}}$ $\displaystyle\sim\text{Inv-
Gaussian}(\mu^{\prime},\lambda^{\prime})$ (33)