Learning Patterns for Detection with Multiscale Scan Statistics

arXiv:1802.06054v2 [math.ST] 1 Mar 2018

James Sharpnack
jsharpna@ucdavis.edu
Department of Statistics, UC Davis

Abstract
This paper addresses detecting anomalous patterns in images, time-series, and tensor data
when the location and scale of the pattern is unknown a priori. The multiscale scan statistic
convolves the proposed pattern with the image at various scales and returns the maximum of
the resulting tensor. Scale corrected multiscale scan statistics apply different standardizations
at each scale, and the limiting distribution under the null hypothesis—that the data is only
noise—is known for smooth patterns. We consider the problem of simultaneously learning and
detecting the anomalous pattern from a dictionary of smooth patterns and a database of many
tensors. To this end, we show that the multiscale scan statistic is a subexponential random
variable, and prove a chaining lemma for standardized suprema, which may be of independent
interest. Then by averaging the statistics over the database of tensors we can learn the pattern
and obtain Bernstein-type error bounds. We will also provide a construction of an -net of
the location and scale parameters, providing a computationally tractable approximation with
similar error bounds.

1

Introduction

Detection is the statistical task of determining if there is some structured signal within noisy data.
If classification answers the question, “what am I seeing?”, detection answers the question, “do I
see anything at all?”. In a sensor network (see [6]), one is often interested in the dual problems
of noticing an anomaly (detection) and then determining its location and extent (classification).
Sensor networks are deployed in natural environments for contaminant detection, [32, 31], realtime surveillance ([4]), radiation monitoring ([3]), and fire detection ([22]). In medical imaging, the
critical task is often to test the existence of an anomaly ([17, 14]). Quick detection of outbreaks
of pathogens, [13, 24], can lead to early intervention. Yet, given that this is such a fundamental
task in many applications, the development of detection methodology lags behind the sophisticated
tools for classification.
In images, time series, and tensors, it is natural to assume that there is a structured signal, such
as blob-like objects, but we do not know its location or scale. In Figure 1, we can see a chemical
plume from a multispectral image where each pixel value is lighter if a certain spectral signature
is present (see [16]). We did not know a-priori where this chemical would appear, nor how large it
would be. A natural approach to detecting signals of this type would be to center a regular shape
such as a square or an ellipse around every pixel and detect anomalously large quantities of the
spectral signature.
The aptly named scan statistic is a test statistic that is based on scanning the image or time
series for a pattern that may be centered anywhere in the domain. For each location, one can
form a likelihood ratio, and test if this statistic exceeds some predetermined threshold (thus, it is a
generalized likelihood ratio test). Scan statistics are widely used in spatial detection applications;
1

Figure 1: An image with an anomalous region of contaminant (left) and a simulated time series
with an embedded sinusoidal signal with values on the y-axis (right).
see [10] for a thorough introduction to the topic. They are commonly used to detect patterns in point
clouds, [18, 20], a closely related problem to our own. Because we often do not know the size as well
as the location of the anomaly, we will scan all locations and scales simultaneously—this is called
the multiscale scan statistic. Hence, the multiscale scan will translate various scaled versions of a
pattern, such as rectangles with many options of side lengths or circles of varying radii. In NeymanPearson testing, we attempt to control the probability of false rejection under a null hypothesis,
when our data only consists of noise. To this end, we either must approximate the distribution of
the scan statistic asymptotically, or set the detection threshold by simulation or resampling. [26]
provided a weak limit for the scan consisting of all intervals in 1 dimension (other approximations
can be found in [19] and [21]). In 2 dimensions, [11, 12, 30, 15], provided approximations of the null
distribution for the multiscale scan statistic. [1, 2] analyzed scan statistics for blob-like patterns
and determined thresholds for detectability in this context.
It was observed that if one naively tests all scales in the multiscale scan at the same threshold,
then the rejection events will be dominated by the finest scale. It was shown in [8] that by separately
standardizing the 1-dimensional scan at each scale one can detect with a signal-to-noise ratio that
adapts to the scale of the anomalous signal. This scale correction for multiscale scan statistics with
rectangular shapes was further studied in [9, 25]. [23] showed that similar results can be achieved
for a large class of patterns that satisfy an average Hölder continuity assumption.
In traditional detection applications, a tensor is scanned for a known function or for specific
blob-like patterns. In a database of many tensors with a repeated anomalous pattern in a database
then it may be possible to simultaneously test for all of the patterns in a class F. The catch is that
each time-series or image—more generally, tensor—will have the anomalous pattern in a different
location and at a different scale. For example, we have many time series, all of which have some
embedded smooth signal (such as the sinusoid in Figure 1) where the sinusoid begins at various
time points, and has different periodicities. Without knowing a priori that we are looking for a
sinusoid this problem seems intractable, but given that we know enough about the signal (such
as it comes from a finite dictionary of functions), then with enough data we may both learn and
detect the pattern within very noisy data.

1.1

Contributions

We will introduce the multiscale scan statistic, with scale correction, in Section 2, and describe an
-net construction using repeated dilation and convolution operations. We begin our theoretical
analysis, in Section 3, by proving a chaining result which shows that the standardized supremum
of certain subGaussian random fields is a subexponential random variable. This is then used, in
Section 4 to provide a finite sample bound on the scale corrected multiscale scan statistic (until
2

now it was only known that it was A.S. bounded). We conclude by demonstrating that our -net
construction is indeed correct and a type 2 error control in Section 5.

2
2.1

Method and Model
Continuous scan statistic

Let’s begin with the basic scan statistic over an image. For an image, which can be represented
by a matrix, {Yk,l : k, l = −L, . . . , L} we can convolve a pattern {Pk,l : k, l = −H, . . . , H} with the
image,
H
X
(P ? Y )k,l =
Yk−k0 ,l−l0 Pk0 ,l0 , k, l = −L + H, . . . , L − H.
k0 ,l0 =−H

Then the simple scan statistic is maxk,l (P ?Y )k,l (one can take the absolute value by also considering
P ← −P ). In the case that Pk,l = 1/(2H) then this scans a square activation pattern over the image.
A common assumption in this problem, is that under the null hypothesis, Yk,l are independent
standard subGaussian random variables. For an arbitrary, pattern matrix P we would like to scale
both dimensions, so that H ← Hj0 in dimension j for some 1 ≤ Hj0 ≤ L (e.g. stretching the
square to form a rectangle). For general patterns, this can be difficult, and one approach is to scale
the dimensions for a continuous function, f , and then rasterize it. This analysis can become very
cumbersome and not terribly enlightening, so we will approximate the pixels with a continuous
domain and the image with a random field.
In order to implement our scan, we begin by proposing a given pattern, which is a function
f ∈ F over the d-dimensional domain. For images, d = 2, and for time-series, d = 1, but we will only
assume that the d is fixed in our asymptotics. We assume that for every f ∈ F, kf kL2 = 1 and is
supported over Ω := [−1, 1]d . We will further assume that f has continuous gradient (F ⊂ C 1 (Ω)).
For a given field, X i , over ΩL := [−L, L]d , we propose a scale parameter,
h ∈ H := ×j [1, L),
such that hj is the scale parameter for dimension j. (Throughout j = 1, . . . d and j will always
indicate dimension.) Given an h ∈ H, we select
t ∈ Th := ×j [−(L − hj ), L − hj ]
respectively, and test if the pattern f centered at t and scaled by h is hidden within image X i .
−1/2
This is accomplished by convolving the field, X i , with the scaled function fh := h• f (./h),
Z
Z
1
i
i
√ f (τ )dX i (t − hτ ),
(fh ? dX )(t) =
fh (τ )dX (t − τ ) =
h•
ΩL
Ω
Q
where h• = j hj and vector operations are such that t/h are performed elementwise. (Throughout,
fh (t) = 0 if t is outside of the domain of fh .) The multiscale scan statistic takes the form,


i
i
s(X ; f ) := max vh max(fh ? dX )(t) − vh .
(1)
h∈H

t∈Th

q P
where in this work we will take vh = 2 j log(n/hj ). Hence, when the scale is coarse (hj is large)
then vh is smaller. At the finest scale, vh is large, the maximum at this scale concentrates about vh
3

with a rate parameter vh (meaning that this concentrates more tightly than the pixel noise variance
of 1).
Throughout, we will denote the set of all valid scale and location parameters as
D := {(t, h) : h ∈ H, t ∈ Th }.
i
i
We will assume that the fields,
R X ,i are independent and have additive white noise terms, W .
Our noise model is such that f dW is a zero mean subGaussian random field indexed by f with
variance kf k2L2 (see Section 3 for a definition). We will test if the field X i is just noise (null
hypothesis) or the noise is added to the function, f , which is translated by ti and scaled by hi ,

H0 : dX i (τ ) = dW i (τ ), i = 1, . . . , n
H1 : dX i (τ ) = µfhi (ti − τ )dτ + dW i (τ )
for some f ∈ F, and (ti , hi ) ∈ D, i = 1, . . . , n.
Critically, the location and scale parameters ti , hi , differ for each image. Notice that kfh kL2 = 1, so
under the null hypothesis, the convolution (fh ? dX i )(t) has mean 0 and variance 1. Because each
statistic, s(X i ; f ), is independent and standardized by vh , then we can average these in order to
increase the power of our final test statistic. To this end, let’s define the pattern adapted multiscale
scan statistic (PAMSS),
n
1 X
Sn (X; F) := max √
sn (X i ; f ).
(2)
f ∈F
n
i=1

We will show that sn (X i ; f ) is subexponential, and so we can obtain probabilistic bounds on Sn
with the subexponential Bernstein-type inequality.

2.2

Multiscale -net construction

In practice, data is discrete, and the scan statistic must be computed over a finite set of scales
and locations. Suppose that for sample, X i , we have a draw from the null hypothesis, so that
EdX i (τ ) = µfhi (ti − τ ). Instead of scanning over all (t, h) ∈ D, we use a finite subset, Dnet ⊂ D.
Let t0 , h0 be values in Dnet that are close to ti , hi in some sense. Then the expectation of the scan
at this approximating location and scale is
Z
i
0
E(fh0 ? dX )(t ) = µ fh0 (τ )fhi (τ + ti − t0 )dτ.
If we define the shift operator (St f )(τ ) := f (τ +t) then we see that this expectation is µhSt0 fh0 , Sti fhi iL2 .
(This metric used for -net construction in [2].) Incidentally, this metric appears when we consider
the variation of our scan statistic under the null hypothesis,
νf ((t, h), (t0 , h0 )) := kSt fh − St0 fh0 kL2 = V (fh ? dW )(t) − (fh0 ? dW )(t0 )

 1
2

,

where dW is the multivariate Wiener process. This fact will be useful for when we provide type
1 error control for our scan statistic. With this metric we will say that a finite subset Dnet ⊂ D
is an -net if for any f ∈ F and any point (t, h) ∈ D there exist points (t0 , h0 ) ∈ Dnet such that
νf ((t, h), (t0 , h0 )) ≤ .
The sensitivity of νf to small changes in t, h will depend on the smoothness of f . If f has large
gradients then a small shift, St , can misalign the function with the unshifted version. To this end,
4

we will consider two different notions of smoothness for the functions in our dictionary, F. Define
the isotropic total variation (recall that the functions have continuous gradients),
Z
k∇f (u)k2 du.
kf kTV :=
Ω

Then we may assume that all of the function are of bounded variation,
∃γ1 > 0 s.t. ∀f ∈ F,

kf kTV ≤ γ1 .

(TVC)

The bounded variation assumption is consistent with the assumption in [8], although we generalize this work to the multidimensional case. Another notion of smoothness is the average Hölder
condition of [23]. Define the Hölder functional,
Z
At,s (f ) :=
|f (t − z) − f (s − z)|2 dz.
ΩL

Then an alternative to the bounded variation assumption is that
∃0 < γ2 ≤ 1 s.t. ∀f ∈ F,

2
At,s (f ) ≤ cA kt − sk2γ
2

(AHC)

where cA is some constant. Throughout the paper we will refer to both conditions, and denote γ
as either γ1 or γ2 depending on context.
Given that our functions F satisfy some smoothness assumptions, we can specify an -net
construction, which was developed first in [29]. For a given β > 1 define `max = blogβ Lc. We will
begin with a subset of scales,
Hβ = {β ` : ` = {0, . . . , `max }d }.
where β > 1 is some parameter. At a given scale, we will consider a grid of evenly spaced locations,
where the distances between grid points increases with the scale, h,
Tα,h := (×j (αhj · Z)) ∩ Th .
The spacing, α > 0, is a tuning parameter as well in this construction. With β, α specified, we can
consider the -net to be
Dβ,α = {(t, h) : h ∈ Hβ , t ∈ Tα,h }.
(Notationally, we exchange Dnet ← Dβ,α .) This is similar to the construction in [25], which has
a fast implementation on GPUs using a hierarchy of convolution and downsampling layers. The
main idea is that instead of expanding the function f → fh can instead downsample the tensor (but
the details on a finite image are somewhat arduous). We expect that a similar implementation is
possible in this context, when the functions are adaptively rasterized, but such developments are
outside of the scope of this paper.
Given an -net, we can compute the approximate scan by restricting the evaluations to the finite
set of location and scales,

eβ,α (X i ; f ) := max vh (fh ? dX i )(t) − vh .
(t,h)∈Dβ,α

Then we similarly can define the -net pattern adapted multiscale scan statistic (-PAMSS),
n

1 X
En (X; F) = max √
eβ,α (X i ; f ).
f ∈F
n
i=1

It is clear that En ≤ Sn since we are maximizing over a strict subset of the continuous scan. Hence,
any type 1 error bound on Sn is also conferred to En .
5

3

A chaining bound for standardized suprema

The scale correction in (1) is based on a precise characterization of the rate and location of a
supremum of the random field resulting from the convolution. This standardization complicates
the analysis significantly, and until now, it was only known that s(X i ; f ) was almost surely bounded.
We initially select the location and scale parameters (t, h) ∈ D and take the supremum over this
selection in forming the scan of f over X i . Because under H0 , the random variables s(X i ; f )
form iid copies of the same random variable, we will seek to exponential concentration inequalities
bounding s. With this in hand, we can hope to obtain PAC-style bounds on their average and
the resulting selection of pattern from a finite function class F. We will restrict this work to
|F| < ∞, but a continuous class of location and scale parameters, D. While this bar seems to be
set pretty low compared to the rich developments in classification, such as Vapnik-Chernovenkis
theory, as we will see, controlling our statistic under finite function classes is a challenging first
step to a more complete understanding of learning patterns in detection. Let us begin with a
formal Rdefinition of a subGaussian random field, and recall that this is what we assume for the
field, { fh dW i (t − .) : t, h ∈ D}, which satisfied when W i is the d-dimensional Weiner process.
Definition 1. We say that a random field, {Z(ι)}ι∈I , is a (zero mean) standard subGaussian
process if there exists a constant u0 > 0 such that


u2
P {|Z(ι0 ) − Z(ι1 )| ≥ u} ≤ 2 exp −
,
(3)
2dZ (ι0 , ι1 )
 2
u
,
(4)
P {Z(ι0 ) ≥ u} ≤ exp −
2
p
for any ι0 , ι1 ∈ I, u > u0 , and dZ (ι0 , ι1 ) = E(Z(ι0 ) − Z(ι1 ))2 , is the canonical distance.
Our noise model assumption can be formally stated as
{(fh ? dW i )(t) : (t, h) ∈ D} is a subGaussian random field with canonical distance νf .

(5)

The generic chaining is a tool for bounding the suprema of random fields with subGaussian
tails, [27] (it also has generalizations to other types of concentration). For general subGaussian
random fields, the expectation of the supremum is bounded by a quantity, E sup Z ≤ TZ , and
one can also show that (sup Z − TZ )/TZ is a standard subGaussian random variable (using the
standard chain construction). Hence, as the expectation bound, TZ , increases the scale can also
grow with TZ (meaning that the supremum becomes more dispersed). This is in contrast to what
we know from extreme value theory, where the maximum of independent random variables tends to
concentrate more tightly (with a lower asymptotic variance), not less. For example, let {zi }N
i=1 be
iid standard normal random
variables.
Then
we
know
that
b
(max
z
−
b
)
approaches
a
Gumbel
i i
N
N
√
distribution where bN = 2 log n + o(1) is a specific sequence (see [7]). Notice that the scale of the
max decreases like 1/bN , which is in contrast to what we obtain from the standard construction
in the generic chaining, which has increasing scale. Multiscale scan statistics also have a Gumbel
limiting distribution, [25], and so we know that the scale of the statistic should be decreasing,
and not increasing. The following theorem significantly modifies the construction in the generic
chaining in order to provide an exponential inequality under parametric-style conditions on the
entropy of the random field, Z.
Theorem 2. Let Z(ι) be a standard subGaussian process over an index set I. Suppose that the
metric space (I, dZ ) has the following bound on the -covering number (ρ > 0),
N (I, dZ , ) ≤ Γ−ρ .
6

(6)

Then there exists an Γ0 > 0 such that for any Γ ≥ Γ0 , the following supremum is bounded in
probability,




p
p
P
c0 log Γ sup Z(ι) − 2 log Γ − a0 log log Γ > u ≤ e−u ,
(7)
ι∈I

for u > u0 where u0 , c0 , a0 are constant depending on ρ (but not on Γ). In words, the supremum of
such a subGaussian process is subexponential with location and rate parameter, (2 log Γ)1/2 .
Proof Sketch. The generic chaining consists of a clever use of the union bound, subGaussian
concentration, and a detailed chain construction. In order to illustrate how we can obtain an
exponential inequality for the max of subGaussian random variables, let us consider the max of
iid standard Gaussian random variables, {zi }N
i=1 . Notice that by union bound and subGaussian
concentration,


p
u2
u2
2
P max zi > 2 log N + u ≤ N e− log N − 2 = e− 2 .
i

p
√
2 ≤ u + 2 log N /(2u) as is done in the generic chaining,
2 log N + u√
Now, instead of bounding,
p
√
we will use the bound 2 log N + u2 ≤ 2 log N + u2 /(2 2 log N ) (we Taylor expand the square
root around 2 log N instead of around u2 ). Hence, we obtain,




p
p
P 2 2 log N max zi − 2 log N > u ≤ e−u .
i

This is the main technique that we use for obtaining subexponential bounds from the max of
independent subGaussian random variables.
The chain construction refers to a sequence of partitionings of the space I. Given a partition
Ak of I, we let Ak (ι) ∈ Ak be the element that contains the point ι ∈ I and ∆(Ak (ι)) be its
k
radius. In the standard generic chaining, a partition is called admissible if |Ak | ≤ 22 . Then the
supremum is controlled by uniformly bounding the centers of the single element A0 ∈ A0 , and then
the differences between the centers of Ak ∈ Ak and their closest centers in at level k − 1. The
standard result is a bound on the supremum based on a functional of the radii ∆(Ak (ι)). For our
modifications, the subexponential bound above requires a growing number of independent points
to work, so we begin our chain at a deeper level than at k = 0. Furthermore, we have to modify
k
the definition of an admissible partition to be |Ak | ≤ aa where a → 1 as Γ → ∞. Technical details
regarding the chain introduce the log log Γ term. See the appendix for a complete proof. Although,
we assume a specific form for the covering numbers, it may be possible to generalize this technique
to other entropy bounds.

4

Type 1 error guarantees for learning patterns

A type 1 error—detecting an anomaly under the null hypothesis—is typically the first error to
be controlled in the Neyman-Pearson testing framework. Our statistic Sn will be compared to a
threshold, which is determined through calibration, simulation, or theoretical guarantees. In this
section, we will provide a finite-sample probabilistic bound for the multiscale scan statistic with
exponential tail probability. We will then use this result to obtain a finite sample bound on the
PAMSS, Sn , which increases logarithmically with the number of functions, log |F|.
Lemma 3. Suppose that f ∈ F satisfies either (TVC) or (AHC), and that W satisfies the noise
assumption, (5). Let ` ∈ {0, . . . , blog2 Lc}d , and H2 (`) = ×j [2`j , 2`j +1 ]. Then when L ≥ L0 ,


P c1 · max vh ((fh ? dW )(t) − vh ) − a1 log log L > u ≤ e−u
(8)
h∈H` ,t∈Th

7

for constants L0 , a1 , c1 > 0 depending on γ, d only.
Proof. Let D0 denote the index set in the above display. Throughout, let c1 and a1 denote arbitrary
constants depending on d, γ alone. By assumption, {(fh ? dW )(t) : (t, h) ∈ D0 } is a subGaussian
random field with canonical distance νf . An -net of D, by definition, will be an -covering of
0
D0 ⊂ D, so we just need to bound the size of the -net, Dβ,α
:= D0 ∩ Dβ,α . By construction,
0
|Dβ,α
|=

X
h∈H2 (`)∩Hβ

X

|Tα,h | ≤ c1

Y L
(L logβ 2)d
P
≤ c1
.
αhj
α d 2 j `j

h∈H2 (`)∩Hβ j

and we can take α, β as specified in the proof of Theorem 6. Furthermore, notice that 2`j is within
0 | ≤ C̃Ld /(h c̃ ). We
a factor of 2 of any hj in H2 (`). Then, there are constants C̃, c̃ such that |Dβ,α
•
can see this because αd (log β)d ≤ αd (β − 1)d ≤ c̃ for some constant depending on d. By Theorem
2, we have the subexponential bound for the random variable,
s X
sup yh ((fh ? dW )(t) − yh ) where yh = 2
log(L/hj ) + 2 log C̃
(t,h)∈D0

j

Notice that yh = vh + O(1/vh ) which gives us the result, along with log vh ≤ C log log L.
From here, we simply apply the union bound with the bound in Lemma 3, for every element
in the partitioning H = ∪` H2 (`). There are on the order of log2 (d log2 L) = log2 log L + O(1) such
elements. This gives us


sn (W i , f )
P c2 ·
− a2 > u ≤ e−u ,
log log L
for some constants c2 , a2 depending on d, γ. Hence, zi (f ) := sn (X i , f )/ log log L is subexponential
with K := kzi (f )kψ1 only depending on d, γ (k.kψ1 is the Orlitz 1-norm). Thus, by the subexponential Bernstein inequality (Prop. 5.16 in [28]),
( n
)


 2
X
t
t
,
,
P
zi (f ) ≥ t ≤ exp −c3 min
nK 2 K
i=1

where c3 is an absolute constant. This gives us our main result (we absorb c3 into K below).
Theorem 4. Let F be finite and assume that either all functions in F satisfy either (TVC) or
(AHC). Let
 r
 

 K log |F | , log |F| ≤ n + log δ
δ
K
Fn (δ) :=
(9)
 

|F
|
K
n
 √ log
, log |F| > K + log δ
δ
n
then for some constants K, L0 depending on d, γ, and L > L0 ,
P {Sn (X, F) > Fn (δ) · log log L} ≤ δ,

(10)

when X is drawn according to the null hypothesis, H0 .
Theorem 4 proves this paper’s main hypothesis, that we can learn patterns from a finite dictionary where the type 1 error bound increases logarithmically with |F|. In fact if log |F| = o(n) then
we have subGaussian concentration of the final test statistic Sn .
8

5

-net approximation and type 2 error

We have provided a construction of the -net with parameters α > 0 and β > 1, but we did not
specify the selection of either or prove the veracity of our claim that this indeed produces an -net.
We used the construction of our -net in the proof of Lemma 3, so we will be careful in this section
to prove the correctness of our construction from first principles. The following technical lemma is
the main driver of these results.
Lemma 5. There is a constant C depending on d alone such that
1. Suppose that (TVC) holds for the class F, then

t − t0 2
h − h0
νf ((t, h), (t0 , h0 ))2 ≤ Cγ1 
+
h 2
h

2

+
2

2. [[23]] Suppose that (AHC) holds for the class F, then

d
d
X
tj − t0j 2γ2 X
hj − h0j
0 0 2

q
νf ((t, h), (t , h )) ≤ C
+
hj
hj h0
j=1

j=1

r

2

j

h0•
h•

!2 
− 1 .

d
X
hj − h0j
q
+
hj h0j
j=1

2γ2 

.

Proof. (2) is proven in [23] (pg. 32) so we will focus on (1). We will partially follow the arguments
in [8] (pg. 145). Let W denote a standard Wiener process in d dimensions. Notice that we can
define νf as,
Z

Z
νf ((t, h), (t0 , h0 ))2 := V
fh (t − z)dW (z) −
fh0 (t0 − z)dW (z) .
ΩL

Recall that

ΩL

1
(fh ? dW )(t) = √
h•

Z
f (u)dW (t − hu).
Ω

By integration by parts,
1
√
h•

Z

Z
f (u)dW (t − hu) =

Ω

Ω

W (t − hu)
√
· ∇f (u)du.
h•

Therefore,
(fh ? dW )(t) − (fh0 ? dW )(t0 ) =

W (t − hu) W (t0 − h0 u)
√
p
−
h•
h0•

Z
Ω

Z
≤
Ω

W (t − hu) W (t0 − h0 u)
√
p
−
h•
h0•

W (t − hu) W (t0 − h0 u)
√
p
≤ sup
−
h•
h0•
u∈Ω

!
· ∇f (u)du

· k∇f (u)k2 du
2

Z
·

k∇f (u)k2 du.
Ω

2

Hence,
0

0

2

νf ((t, h), (t , h )) ≤

kf k2TV

W (t − hu) W (t0 − h0 u)
√
p
· E sup
−
h•
h0•
u∈Ω
9

2

.
2

Decompose the supremum term,
W (t − hu) W (t0 − h0 u)
√
p
−
h•
h0•

E sup
u∈Ω

W (t0

W (t − hu) −
√
h•

≤ E sup
u∈Ω

−

2
h0 u)

+ kW (t0 − h0 u)k2
2

1
1
√ −p
h•
h0•

!

0


√
0
d
By Brownian scaling, W (t − hu) − W (t0 − h0 u))/ h• = W ht − u − W th − hh u . Hence, the
jth coordinate
inthe

 above LHS is a 1-dimensional Brownian motion equal in distribution to,
h0
t−t0
B h − 1 − h u . By the reflection principle, we have that for some constant c0 ,

E sup B
u∈Ω

 2

t − t0
h0
u ≤ c0
− 1−
h
h

tj − t0j
hj

2

hj − h0j
+
hj

2

!
.

Hence,
E sup
u∈Ω

W (t0

W (t − hu) −
√
h•

−

h0 u) 2
2



X hj − h0j
≤ c0 
hj
j

2


tj − t0j 2
.
+
hj

Without loss of generality, we can select t = 0 by translation invariance of νf . Furthermore by a
similar argument, E supu kW (h0 u)k22 ≤ c0 h0• , for another universal constant c0 . Hence,
0

0

E sup kW (t − h u)k2
u∈Ω

1
1
√ −p
≤ c0
0
h•
h•

r

h0•
−1 .
h•

Combining these bounds completes the proof.
Theorem 6. Suppose that either one of (TVC) or (AHC) holds. Let  > 0, then there exists a
β > 1, α > 0 depending on d,  such that for any (t, h) ∈ D, there exists (t0 , h0 ) ∈ Dβ,α (the -net)
with νf ((t, h), (t0 , h0 )) ≤ .
Proof. Throughout, let γ mean either γ1 , γ2 depending on context. By the triangle inequality let
us bound,
νf ((t, h), (t0 , h0 )) ≤ νf ((t, h0 ), (t0 , h0 )) + νf ((t, h), (t, h0 )).
In our construction, notice that for any t there is a grid point in Tα,h0 that is within αh0j from it in
dimension j. Hence,
2γ2
2
d
d
X
X
tj − t0j
tj − t0j
2γ2
≤ dα ;
≤ dα2 .
h0j
h0j
j=1

j=1

Furthermore, by the -net construction, there exists an h0j such that | log hj /h0j | ≤ log β for every
j. Thus,
r
hj − h0j
d
h0•
q
≤ β − 1;
− 1 ≤ β 2 − 1.
0
h
•
hj hj
Thus there is a constant, C, depending on γ, such that for  small enough α = C1/γ2 and β =
C(1 + )2/d is sufficient.

10

Suppose that we are under the alternative hypothesis, H1 , so that there is some embedded
signal f in each image. Consider evaluating the scan at the true location and scale, (ti , hi ) for a
given field X i , then (fhi ? dX i )(ti ) is normally distributed with mean µ and variance 1. In the event
that we scan over an -net, then by arguments in Section 2, there is an element in the approximate
scan with mean µ(1 − 2 /2). Hence, we have the following type 1 error bound, by summing the
resulting n normal random variables.
i n
Proposition 7. Suppose that {X
are drawn from
P}i=1
P H1 (with possibly different location and scale
2
parameters). Then define Vn = i vhi and Wn = i vhi ,
(
r )


 
1
2
Vn
Wn − Vn < u
≤ Φ(u)
P En (X; F) − √
µ· 1−
2
n
n

where Φ is the standard Normal CDF. The above display is also true if we let  = 0 and substitute
En,β ← Sn .
In order to have diminishing type 1 error probability, we set a threshold for Sn at Fn (δ)·log log L
for δ → 0. Assume that  → 0 (however slowly) then to have the type 2 error probability to decrease
as well, we require that
r !
µWn − Vn
Vn
√
− Fn (δ) log log L = ω
n
n
Suppose that n = 1 and |F| = 1 (the standard multiscale scan setting), then this would require
µ − vh1 −

K
1
log · log log L = ω(1),
vh1
δ

which is consistent up to some constants with previously known rates, [23]. Critically, the dependence on δ is logarithmic, and this provides the first such finite sample bound.
Let us conclude with a remark about the restrictiveness of the assumption that we have a finite
function class F. It is known that functions of bounded variation have Haar wavelet coefficients
that are bounded in a weak `1 norm, [5]. It is reasonable to discretize the allowed coefficient values
and then restrict our function class to functions with k-sparse wavelet coefficients of m then the
size of the class scales like k log m which is very managable. One advantage with this approach is
that the sparse Haar wavelets will naturally satisfy condition (TVC). It is outside of the scope of
this work to extend the result to infinite function classes, but this would present a very interesting
and important extension.

6

Conclusions

We have addressed learning and detecting patterns from a function class, F, using multiscale
scan statistics. We have introduced the multiscale scan statistic and proved a subexponential
concentration bound for it, which relied on a novel chaining result for standardized suprema of
subGaussian random fields (a result that may be of independent interest). We introduced the
pattern adapted multiscale scan statistic, that can learn patterns in a database of tensors (when
the locations and scales vary). This result allowed us to prove Bernstein-type concentration for
the PAMSS, meaning that we can learn finite function classes that grow exponentially with the
sample size, n. With evidence that representation learning and detection are not incompatible, we
anticipate that efficient methods for learning functions in this setting will emerge, by using modern
tools from deep learning and multiscale methods.
11

Acknowledgements: JS is supported in part by DMS-1712996.

References
[1] E. Arias-Castro, D.L. Donoho, and X. Huo. Near-optimal detection of geometric objects by
fast multiscale methods. IEEE Trans. Inform. Theory, 51(7):2402–2425, 2005.
[2] Ery Arias-Castro, Emmanuel J. Candès, and Arnaud Durand. Detection of an anomalous
cluster in a network. Ann. Statist., 39(1):278–304, 2011.
[3] Sean M. Brennan, Angela M. Mielke, David C. Torney, and Arthur B. Maccabe. Radiation
detection with distributed sensor networks. Computer, 37(8):57–59, 2004.
[4] Y. Caron, P. Makris, and N. Vincent. A method for detecting artificial objects in natural
environments. In Proceedings 16th International Conference on Pattern Recognition, volume 1,
pages 600–603. IEEE Comput. Soc., 2002.
[5] Albert Cohen, Wolfgang Dahmen, Ingrid Daubechies, Ronald DeVore, et al. Harmonic analysis
of the space bv. Revista Matematica Iberoamericana, 19(1):235–263, 2003.
[6] D. Culler, D. Estrin, and M. Srivastava. Overview of sensor networks. IEEE Computer,
37(8):41–49, 2004.
[7] Laurens De Haan and Ana Ferreira. Extreme value theory: an introduction. Springer Science
& Business Media, 2007.
[8] Lutz Dumbgen and Vladimir G Spokoiny. Multiscale testing of qualitative hypotheses. Annals
of Statistics, pages 124–152, 2001.
[9] Lutz Dümbgen and Günther Walther. Multiscale inference about a density. The Annals of
Statistics, pages 1758–1785, 2008.
[10] Joseph Glaz, Joseph I Naus, and Sylvan Wallenstein. Scan statistics. Springer, 2001.
[11] Joseph Glaz and Zhenkui Zhang. Multiple window discrete scan statistics. Journal of Applied
Statistics, 31(8):967–980, 2004.
[12] G Haiman and C Preda. Estimation for the distribution of two-dimensional discrete scan
statistics. Methodology and Computing in Applied Probability, 8(3):373–382, 2006.
[13] R. Heffernan, F. Mostashari, D. Das, A. Karpati, M. Kulldorff, and D. Weiss. Syndromic
surveillance in public health practice, New York City. Emerging Infectious Diseases, 10(5):858–
864, 2004.
[14] D. James, B. D. Clymer, and P. Schmalbrock. Texture detection of simulated microcalcification
susceptibility effects in magnetic resonance imaging of breasts. Journal of Magnetic Resonance
Imaging, 13(6):876–881, 2001.
[15] Zakhar Kabluchko. Extremes of the standardized gaussian noise. Stochastic Processes and
their Applications, 121(3):515–533, 2011.

12

[16] Dimitris Manolakis and Gary Shaw. Detection algorithms for hyperspectral imaging applications. IEEE signal processing magazine, 19(1):29–43, 2002.
[17] Nathan Moon, Elizabeth Bullitt, Koen van Leemput, and Guido Gerig. Automatic brain
and tumor segmentation. In MICCAI ’02: Proceedings of the 5th International Conference on
Medical Image Computing and Computer-Assisted Intervention-Part I, pages 372–379, London,
UK, 2002. Springer-Verlag.
[18] Joseph I. Naus. The distribution of the size of the maximum cluster of points on a line. J.
Amer. Statist. Assoc., 60:532–538, 1965.
[19] Joseph I Naus and Sylvan Wallenstein. Multiple window and cluster size scan procedures.
Methodology and Computing in Applied Probability, 6(4):389–400, 2004.
[20] Daniel B Neill. Fast subset scan for spatial pattern detection. Journal of the Royal Statistical
Society: Series B (Statistical Methodology), 74(2):337–360, 2012.
[21] Vladimir Pozdnyakov, Joseph Glaz, Martin Kulldorff, and J Michael Steele. A martingale
approach to scan statistics. Annals of the Institute of Statistical Mathematics, 57(1):21–37,
2005.
[22] D. Pozo, F.J. Olmo, and L. Alados-Arboledas. Fire detection and growth monitoring using a
multitemporal technique on AVHRR mid-infrared and thermal channels. Remote Sensing of
Environment, 60(2):111–120, 1997.
[23] Katharina Proksch, Frank Werner, and Axel Munk. Multiscale scanning in inverse problems.
arXiv preprint arXiv:1611.04537, 2016.
[24] L.D. Rotz and J.M. Hughes. Advances in detecting and responding to threats from bioterrorism
and emerging infectious disease. Nature Medicine, pages S130–S136, 2004.
[25] James Sharpnack and Ery Arias-Castro. Exact asymptotics for the scan statistic and fast
alternatives. Electronic Journal of Statistics, 10(2):2641–2684, 2016.
[26] David O Siegmund and Keith J Worsley. Testing for a signal with unknown location and scale
in a stationary gaussian random field. The Annals of Statistics, pages 608–639, 1995.
[27] Michel Talagrand. The generic chaining: upper and lower bounds of stochastic processes.
Springer Science & Business Media, 2006.
[28] Roman Vershynin. Introduction to the non-asymptotic analysis of random matrices. arXiv
preprint arXiv:1011.3027, 2010.
[29] Guenther Walther. Optimal and fast detection of spatial clusters with scan statistics. Ann.
Statist., 38(2):1010–1033, 2010.
[30] Xiao Wang and Joseph Glaz. Variable window scan statistics for normal data. Communications
in Statistics-Theory and Methods, 43(10-12):2489–2504, 2014.
[31] Brian A White, Antonios Tsourdos, Immanuel Ashokaraj, S Subchan, and Rafal Zbikowski.
Contaminant cloud boundary monitoring using network of uav sensors. IEEE Sensors Journal,
8(10):1681–1692, 2008.

13

[32] Y Jeffrey Yang, Roy C Haught, and James A Goodrich. Real-time contaminant detection and
classification in a drinking water pipe using conventional water quality sensors: Techniques
and experimental results. Journal of environmental management, 90(8):2494–2506, 2009.

A

Proof of Theorem 2

Proof. We will follow the construction of generic chains as in [27], but will be significantly more
careful about the details of the construction. This proof is similar in spirit to [8], but we will get
probabilistic bounds and prefer this proof because it uses only first principles. Throughout this
proof we will call variables that dependent only on ρ, constants, and some like C may change from
line to line. Variables from the main body of the paper, other than those defined in Theorem 2,
may appear in this proof and mean something different (we suppose that they are in a different
scope).
Let’s begin by defining G = Γ1/ρ and
a = (1 − log−1 G)−1 > 1.
k

Let an admissible partition, Ak , be any partition of I of size at most aa . Let Ak (ι) be the element
of the partition containing ι and let the center of this element be
τn (ι) :=

sup dZ (τ, ι0 ),

inf

τ ∈An (ι) ι0 ∈An (t)

and ∆(Ak (t)) = supι∈Ak (t) dZ (τk (ι), ι) be its radius. Let k0 = 1, . . . , ∞ then Z(ι) ≤ Z(τk0 (ι)) +
P
∞
k=k0 |Z(τk+1 (ι)) − Z(τk (ι))|. By the union bound and (4),

n
o

p
1
ak0
k0
2
k
2
0
P ∃ι ∈ I : Z(τk0 (ι)) > 2 log a · a + u ≤ a exp − 2 log a · a + u
2
= e−
Let

u2
2

.

k

− aρ

k = G · a

,

(11)

then there exists an admissible partition where the radius of the balls are k satisfy N (I, dZ , k ) ≤
k
aa . Hence,
n
o
p
P ∃ι ∈ I : |Z(τk+1 (ι)) − Z(τk (ι))| > k 2 log a · ak+1 + dk u2


1
ak+1
k+1
2
≤ 2a
+ dk u
exp − 2 log a · a
2
= 2e−
where

dk u2
2

,
1

dk = e 2 (1+log a)·(k−k0 ) .
Define the quantities,
Ak0 =

∞
p
p
X
k 2 log a · ak+1 + dk u2 ,
2 log a · ak0 + u2 +
k=k0

Bk0

 2 X


∞
u
u2
= exp −
+
exp −dk
,
2
2
k=k0

14

(12)

so that
P{sup Z(ι) > Ak0 } ≤ 2Bk0 .

(13)

ι∈I

Because

√

b1 + b2 ≤

√

Ak0

√
b1 + b2 /(2 b1 ) for b1 , b2 > 0 we have
∞
a−k0 /2 2 X
≤ a0 ak0 /2 +
u +
2a0

k+1

a0 k a

a− 2
+
k dk u2
2a0

k+1
2

k=k0

for a0 =

√

!
(14)

2 log a. Let k0 satisfy
ak0 ≤ ρ(loga G + k0 ) ≤ ak0 +1 ,

(15)

such that
k0 log a ≤ C 0 log log G
for some constant C 0 possibly depending on ρ, which is guaranteed by Lemma 8 for Γ large enough.
Then we have that
a0 a

k0
2

≤

p
p
ρk0 log a p
C 0 log log Γ
2 log Γ + 2ρk0 log a ≤ 2 log Γ + √
≤ 2 log Γ + √
2 log Γ
2 log Γ

(16)

Furthermore, by (11) and (15)
∞
X

a

k+1
2

k ≤ G

k=k0

∞
X

k

− aρ + k+1
2

= Ga

a

k0 +1
2

j

a

−ak0 aρ + 2j

j=0

k=k0

≤ Ga

∞
X

k0 +1
2

∞
X

j−1 + j
2

a−(loga G+k0 )a

=a

1−k0
2

j=0

=a

1−k0
2

∞
X

j−1 )+ j
2

a(loga G+k0 )(1−a

j=0



∞

1−a−1 √ X
j
j
 Gak0
+ a
a(loga G+k0 )(1−a )+ 2  .
j=0

The first term is
a

1−k0
2



Gak0

1−a−1

=e·a

1−k0
2

k0

a log G .

Furthermore, for any C < 1 there exists a G0 > 0 such that log a ≥ C/ log G for all G ≥ G0 , thus
k0
 1
1
a log G ≤ (ρ(loga G + k0 )) log G ≤ ρ(C log2 G + k0 ) log G → 1,

as G → ∞. Considering the second term,
∞
X

j )+ j
2

a(loga G+k0 )(1−a

≤1+

j=0

∞
X

j )+ j
2

a(loga G)(1−a

j=1

Using the first order approximation,
j

a −1=



1
1−
log G

15

−j
−1≥

j
.
log G

.

So we may bound the series,

  

−1
∞
∞
X
X
log a
log a
1
(loga G)(1−aj )+ 2j
≤
exp
− 1 j ≤ 1 − exp
−1
,
a
→
2
2
1 − e−1
j=0

j=0

as G → ∞. Hence, there are constants C, G0 > 0 depending on ρ such that for G > G0 ,
s
√
∞
X
k0
k+1
2a log a
C 2a log a
−
a0 k a 2 ≤ Ca0 a 2 ≤ C
.
≤ √
ρ(loga G + k0 )
2ρ
log
G
k=k

(17)

0

Also,
k0

C
a− 2
≤√
.
a0
2ρ log G

(18)

Consider the fourth term in Ak0 ,
∞
X

a

− k+1
2

∞
X

k dk = G

k=k0

a−

k+1
2

k

a

− aρ

1

e 2 (1+log a)·(k−k0 )

k=k0
1

≤ Ga− 2 (k0 +1)
3

1

≤ a− 2 k0 − 2

∞
X

j

j−1

a− 2 a−(loga G+k0 )a

j=0
∞
X j
− 2 (loga G+k0 )·(1−aj−1 )

a

a

1

e 2 (1+log a)·j
1

e 2 (1+log a)·j

j=0

Isolating the summation,
 

∞
∞
X
√ X
1 log G + k0 log a
(loga G+k0 )·(1−aj−1 ) 2j
−
·j
exp
a
e ≤ e
2
log G
j=0
j=1
√
e
√ ≤ 4.2.
≤
1 − 1/ e
Furthermore,
a

(loga G+k0 )·(1−a−1 )

≤e



k
− 1+ log 0 G
a

≤

1
< 0.4.
e

Hence, there exists a G0 > 0 such that
∞
3
1
u2 X − k+1
u2
4.6
a 2 k dk ≤ u2 √
· a− 2 k0 − 2 ≤ √
,
a0
2
log
a
ρ
log
G
k=k

(19)

0

for G ≥ G0 . Hence, by (14) combined with (16), (17), (18), and (19) yields
p
a0 log log Γ + u2 /2
√
Ak0 ≤ 2 log Γ +
,
c0 log Γ
for constants C0 , C1 , G0 with G ≥ G0 . Consider Bk0 ,

 X


 2 
 2 X
∞
∞
u
u2
1
u
exp
·
exp −dk
=
exp
1 − exp
(1 + log a) · j
·
2
2
2
2
j=0
k=k0



∞
∞
√ X
√ X
1
1
j
≤ e
exp − exp
j
= e
b−b ≤ 2.4,
2
2
j=0

(20)

(21)

j=0

(22)
16

where b =

√

e. Hence,
1

2

Bk 0 ≤ C · e − 2 u .
Recalling (13) gives us our result.
Lemma 8. Let ρ > 0, and define the function

a(G) = 1 −

1
log G

−1
.

(23)

There exists a constant G0 depending on ρ alone such that if G > G0 then there exists an integer
k0 satisfying
ak0 ≤ ρ(loga G + k0 ) ≤ ak0 +1 ,
(24)
and there is a constant C depending on ρ such that k0 log a ≤ C log log G.
Proof. Set G > 0 and a = a(G). Let κ1 be the root of the following function,
f (κ) := aκ − ρ(loga G + κ).
Notice that f is strictly convex and so let κ0 be its unique minimizer.
κ1 = loga (ρ(loga G + κ1 )) ≥ loga (ρ loga G) ,
Notice that log a ≤ 1/ log G, hence,
loga G =
and
κ1 ≥

log G
≥ log2 G,
log a

log log G
log ρ
+2
≥ log G · (log ρ + 2 log log G) =: x1 .
log a
log a

Moreover, κ0 satisfies κ0 aκ0 = ρ,
(κ1 − 1)aκ1 −1 =

ρ
ρ
(κ1 − 1)(loga G + κ1 ) ≥ (x1 − 1)(log2 G + x1 ).
a
a

The limit of the right hand side approaches +∞ as G increases. Hence, there exists an G0 > 0
depending on ρ only such that for any G ≥ G0 ,
(κ1 − 1)aκ1 −1 ≥ ρ = κ0 aκ0 .
Because κaκ is a monotone increasing function of κ > 0, this implies that κ1 − 1 ≥ κ0 .
Let k0 = bκ1 c. We know that f is increasing and convex over [κ0 , ∞). Let G ≥ G0 . Because
κ0 ≤ k0 ≤ κ1 ≤ k0 + 1 we have that f (k0 ) ≤ 0 ≤ f (k0 + 1). Hence,
ak0 − ρ(loga G + k0 ) ≤ 0 ≤ ak0 +1 − ρ(loga G + k0 + 1) ≤ ak0 +1 − ρ(loga G + k0 ),
which proves the first result. Finally,


log(a2k0 ) ≤ log ak0 − ρk0 ≤ log(ρ loga G)
for some G (k0 ) large enough. Notice that loga G ≤ C(log G)2 by the limit of log a → 1/ log G,
demonstrating the second result.
17

A note on the derivation of Theorem 4. Notice that t2 /nK < t/K when t < n, and that
let
s
|F|
nK
log
t=
.
c3
δ
Then


t2
≤δ
|F| exp −c3
nK 2
and t < n when
log |F| ≤

n
+ log δ.
c3 K

In the other case, we define t so that the probabilistic bound holds, and absorb c3 into K.

18

