arXiv:1802.04838v1 [stat.ML] 13 Feb 2018

Network Estimation from Point Process Data
Benjamin Mark1,4 , Garvesh Raskutti2,4 , and Rebecca Willett3,4
1
Department of Mathematics, 2 Department of Statistics,
3
Department of Electrical and Computer Engineering
4
Wisconsin Institute for Discovery
University of Wisconsin-Madison
February 15, 2018
Abstract
Consider observing a collection of discrete events within a network that
reflect how network nodes influence one another. Such data are common in
spike trains recorded from biological neural networks, interactions within a
social network, and a variety of other settings. Data of this form may be modeled as self-exciting point processes, in which the likelihood of future events
depends on the past events. This paper addresses the problem of estimating
self-excitation parameters and inferring the underlying functional network
structure from self-exciting point process data. Past work in this area was
limited by strong assumptions which are addressed by the novel approach here.
Specifically, in this paper we (1) incorporate saturation in a point process
model which both ensures stability and models non-linear thresholding effects;
(2) impose general low-dimensional structural assumptions that include sparsity, group sparsity and low-rankness that allows bounds to be developed in
the high-dimensional setting; and (3) incorporate long-range memory effects
through moving average and higher-order auto-regressive components. Using
our general framework, we provide a number of novel theoretical guarantees
for high-dimensional self-exciting point processes that reflect the role played
by the underlying network structure and long-term memory. We also provide
simulations and real data examples to support our methodology and main
results.

1

1

Introduction

In a variety of settings, our only glimpse of a network’s structure is through the
lens of discrete time series observations. For instance, in a social network, we may
observe a time series of members’ activities, such as posts on social media. In
electrical systems, cascading chains of power failures reveal critical information
about the underlying power distribution network. During epidemics, networks of
computers or of a population are reflected by the time at which each node becomes
infected. In biological neural networks, firing neurons can trigger or inhibit the firing
of their neighbors, so that information about the network structure is embedded
within spike train observations.
This paper focuses on estimating the influence network which models the extent
to which one node’s activity stimulates or inhibits activity in another node. For
instance, the network structure may indicate who is influencing whom within a social
network [50, 55, 57, 26, 8, 65], the connectivity of neurons [9, 27, 12, 56, 32, 14, 49,
38], interactions among financial instruments [10, 4, 39], how power failures may
propagate across the power grid [17], or patterns of criminal activity and military
engagements [57, 8, 16, 37, 39]. The interactions between nodes are thus critical
to a fundamental understanding of the underlying functional network structure and
accurate predictions of likely future events.
Learning the influence network presents a number of challenges both in terms
of formulating the model and developing suitable theory and methodology. First,
in the applications described above the number of network nodes is typically large
relative to the length of time they are observed, making the network parameter highdimensional. Furthermore, the most natural model in these settings are multivariate
self-exciting point processes (SEPPs). While empirical work has demonstrated the
efficacy of SEPP models in various applications (cf., [16, 10, 37, 39, 17]), little is
known about the statistical properties of these estimators. In this paper, we formulate
a model and provide a general framework for estimating network parameters in
discrete-time high-dimensional SEPP models.
Let M denote the number of nodes in the network and T the number of time
intervals over which we collect data. We observe Xt,m , the number of events at node
m during time period t, for m = 1, . . . , M and t = 1, . . . , T . We model these counts
as
Xt,m ∼ Poisson(λt,m )
where the logarithm of λt,m is a function of the previous counts of events in the
network and the interactions between nodes. For a simple example, we might have
P
log λt,m = m
m0 =1 Am,m0 Xt−1,m0 . However, a fundamental challenge associated
with SEPP models is that they can be highly unstable: due to the exponential link

2

function, the counts can diverge even when the interactions {Am,m0 } are small. In
[23] the authors give extensive justification for the interest in these models from a
neuroscience perspective, but also show how learned model parameters can result
in generative models that are highly inconsistent with physiological measurements.
Existing statistical learning bounds for SEPP models [28] guarantee stability by
assuming all network interactions are inhibitory.
A major contribution of this work is learning guarantees for SEPPs without
restrictive assumptions on the structure of the network or types of interactions
among nodes. We will address stability issues by introducing saturation effects
on the rate parameter λt,m . Saturated SEPP models were recently described in
application-driven work without theoretical guarantees [17]. In contrast, this work
aims to derive statistical learning guarantees for saturated point processes.
We study a fairly general class of saturated SEPPs whose parameters can be
estimated via regularized maximum likelihood estimation. We assume that the
number of possible interactions between nodes (i.e., graph edges) M 2 is large
relative to the number of time points T , but that the network has an underlying
low-dimensional structure that can be promoted via regularization. The question we
address in our theory is how many time points T are needed to guarantee a desired
level of statistical accuracy in terms of the number of nodes M , the underlying
network structure, the regularizer used, and the type of saturation effects introduced?

1.1

Relationship to Prior Work

A number of works have studied linear SEPPs (where λt,m is a linear function of
past events, in contrast to log-linear models, where log λt,m is a linear function of
past events) from a theoretical perspective. Examples include works on the Hawkes
process [11, 15, 54, 29, 6]. In a multivariate Hawkes process setting, one frequently
aims to learn the excitation matrix characterizing interactions within the network.
In [18] the authors establish that learning the excitation matrix is sufficient for
learning the directed information graph of the process. The linear Hawkes process is
frequently studied under an assumption there are no inhibitory interactions, although
recent work [11] was able to incorporate both inhibitory and stimulatory interactions.
Prior work on learning parameters in discrete high-dimensional time series models
requires linearity or Gaussianity assumptions (cf., [7]) which do not hold in our
model.
In contrast, we study log-linear SEPPs. Prior works have demonstrated the
empirical value of log-linear SEPPs [41, 60] and these models are frequently used in
the neuroscience community [23]. Moreover, log-linear point process models can
be advantageous from the perspective of optimization [47] and naturally allow for
inhibitory interactions. However, log-linear SEPPs can not easily model stimulatory
3

interactions while maintaining stability, and incorporating stimulatory interactions is
a major contribution of this paper.
There is limited work on learning rates for log-linear SEPPs, and much of it
is only applicable in the setting where M is small relative to T [20]. The most
related work is our recent work [28] which considers a special case of our SEPP
along with a sparsity assumption on the network and applies in the high-dimensional
setting. This prior work is limited since the model only considers recent memory,
sparsity regularization, and assumes only inhibitory influences to ensure stability
and learnability.

1.2

Main Contributions

Our paper makes the following major contributions.
• We provide a general upper bound (Theorem 3.1) for developing theoretical
guarantees for estimating SEPPs and build on the analysis in [28] in three
significant ways. First, we incorporate saturation effects in our model by using
a thresholding function in order to ensure stability, and account for these effects
in our theory. Second, we provide learning rates for a class of processes which
incorporate longer-range dependence effects in a variety of ways, improving
upon [28] which only considers first-order auto-regressive models. Finally,
we allow for several different regularization choices corresponding to various
prior beliefs about the structure of the network.
• We apply our general upper bound to a number of different processes and
regularization schemes. For processes with longer-range dependence, we
prove that a restricted eigenvalue condition holds for the ARMA(1, 1) and
AR(2) models in Lemmas 4.1 and 4.2 respectively.
• In terms of regularization schemes, we consider strict sparsity, group sparsity
and low-rank regularization and provide three novel guarantees stated in
Theorems 4.4, 4.6 and 4.8. All our mean-squared error bounds match the
optimal bounds in the independent case up to log factors.
• A thorough simulation study in Section VI provides support for our theoretical
mean-squared error bounds and also examines parameters associated with the
magnitude of the entries of A∗ and clipping thresholds.
• We further demonstrate the practical benefits of our regularized likelihood
framework on three real data examples. The first involves modeling the
interplay between crime events in different neighborhoods of Chicago, the
4

second modeling connections between different neurons in the brain within a
rat during sleep and wake states, and the third involving meme-tracker data
in social networks. The three examples illustrate the advantages of using
different regularizers.
• Finally, we show that our SEPP framework can be viewed as a discretization of the widely studied Hawkes process, and discuss some advantages of
considering point processes in discrete time.

1.3

Notation

For a matrix A, we let am. denote the mth row of A and a.m denote the mth column
of A. We then let kam. k1+ denote the sum of the positive entries of am. and kam. k1−
denote the absolute value of the sum of the negative entries of am. , so that
kam. k1 = kam. k1+ + kam. k1− .
Given a norm k · kR on a real vector space, we let k · kR∗ denote its dual norm defined
by
kvkR∗ = sup < u, v >
kukR ≤1

where < ·, · > denotes the standard inner product. Throughout the paper, we work
with mixed norms
X
1
kAkp,q = ( kam. kpq ) p ,
m

as well as the nuclear norm
kAk∗ =

M
X

σi (A)

i=1

where σi (A) denotes the ith singular value of A, and the operator norm
kAkop = sup kAxk2 .
kxk2 ≤1

The Frobenius norm, denoted by k · kF , is a special case of the k · kp,q norm with
p = q = 2.
Finally, we let kAk0 denote the number of nonzero elements of a matrix A.

2

Model Formulation

In this section, we present a class of SEPPs and discuss how saturation effects can be
included in order to ensure stability. Recall that Xt,m denotes the number of events
5

from node m during time period t. To start, consider the following model:
Xt+1,m ∼ Poisson(λt+1,m )
t X
M
X

log(λt+1,m ) = νm +

(2.1)

hm,m0 [t − s]Xs,m0 .

(2.2)

s=1 m0 =1

Here the logarithm of the rate for Xt+1,m is linear in all the previous observations.
For each node m0 in the network, that node’s count Xs,m0 at time s is scaled by an
influence function hm,m0 evaluated at t − s. The influence function hm,m0 describes
the relationship between nodes m and m0 . As in [11], we assume each influence
function can be written as the linear combination of K known basis functions
{φk }K
k=1 , i.e.,
hm,m0 [t] =

K
X

am,m0 ,k φk [t].

k=1

Hence estimating the network structure amounts to estimating the matrix A∗ ∈
RM ×M K where the mth row of A∗ is



a

rewrite (2.2) in matrix-vector form as

M
m,m0 ,k

K

m0 =1 k=1

log(λt+1 ) = ν + A∗ g(Xt ),

It will be convenient to

(2.3)

where Xt = [X1 , . . . , Xt ] denotes the history of the process up to time t and
g(XT ) ∈ RM K×1 is the vector defined as follows. For k = 1, . . . , K, let

gk (XT ) :=

 PT
s=1 Xs,1 φk [T
 PT
 s=1 Xs,2 φk [T

..


.
PT

− s]
− s] 






(2.4)

Xs,M φk [T − s]

s=1

and




g1 (XT )


 g2 (XT ) 

g(XT ) :=  .. 
.

. 

(2.5)

gK (XT )
A number of commonly studied discrete time models can be realized in this
manner. We briefly mention two which are discussed further in Section IV. As
a first example, K = 1 and φ[t] = αt corresponds to an autoregressive moving
average ARMA(1, 1) process. When K = p and φk [t] = I{k=t} where I{B} :=
6

(

1, B true
is the indicator function, we recover the AR(p) process. This second
0, otherwise

example shows the value in assuming that hm,m0 is in the span of a collection of
basis functions, rather than just one. Allowing for multiple basis functions allows
us to study processes which incorporate higher order effects in more sophisticated
ways than would be possible with only one basis function.
We let νmin and νmax be upper and lower bounds on the constant offset parameter
νm in (2.2) and we assume that A∗ lies within a set A which we define as follows.
Let amax be an upper bound on ka∗m. k1+ and similarly let amin be an upper bound on
ka∗m. k1− . We let A denote the set of M × M K matrices with kam. k1− ≤ amin and
kam. k1+ ≤ amax for all m. With the assumption that A∗ ∈ A we can search for an
estimate Ab of A∗ over the bounded set A.

2.1

Saturation

As discussed in the introduction, point process models along the lines of (2.3) are
widely used to describe count data in a variety of applications. However, due to
instability issues inherent to SEPPs of this form, these models can be highly unstable
and lead to unbounded counts. Hence, pure SEPPs make poor generative models (c.f.,
[23]) and are difficult to understand theoretically without making overly restrictive
assumptions about A∗ (c.f, [28]). We will address this problem by introducing
saturation effects to the vector g(Xt ) defined in Equation (2.5). The application
focused work [17] introduced saturated SEPPs, but to the best of our knowledge, this
is the first work to study the theoretical properties of saturated models. To address
stability issues we adjust the definition of gk (XT ) in (2.4) to the following:
 PT

s=1

gk (XT ) =

min(Xs,1 , Ũ )φk [T − s]

P
 T
 s=1 min(Xs,2 , Ũ )φk [T

..


.
PT
s=1




− s] 


.



(2.6)

min(Xs,M , Ũ )φk [T − s]

That is, each past count which exceeds some threshold Ũ ≥ 1 gets clipped to Ũ .
Further, we assume that
∞
X

φk [s] ≤ τ < ∞

s=1

for each basis function, so that each entry of g(Xt ) in (2.5) is bounded by
Ũ

T
X

φk [s] ≤ τ Ũ =: U.

s=1

7

In other words, with clipping we have kg(Xt )k∞ ≤ U , guaranteeing the stability of
our process.
In particular, this allows us to define the maximum and minimum Poisson rate
from which each observation can be drawn. We denote the maximum and minimum
rates by
Rmax = exp (νmax + amax U )
Rmin = exp (νmin − amin U ) .

(2.7a)
(2.7b)

Throughout this paper, we take min(., Ũ ) to be our saturation function for
simplicity. However, our theory extends to other saturation functions provided that
the function is bounded, which is crucial for our analysis. The details are provided
in Proposition 2 in the appendix.
While this framework has advantages, a central question we need to address
is how departing from the standard SEPP framework and incorporating non-linear
saturation effects change our estimation errors.

2.2

Regularized optimization formulation

In the high-dimensional setting, the number of potential pairwise interactions, M 2 , is
large relative to the number of time periods, T , making standard maximum likelihood
optimization techniques unsuitable. Instead, we assume some prior knowledge on
the parameter A∗ , which can be incorporated in estimation via a regularization term
k.kR . Specifically, we consider the estimator
Â = arg min
A∈A

> X
M
X





>
exp νm + a>
m. g(Xt ) − Xt+1,m (νm + am. g(Xt )) + λkAkR

t=0 m=1

(2.8)
where the first two terms of (2.8) are the negative log-likelihood of the observed data
given A. We discuss various choices of regularization k.kR in the next section. Note
that the optimization problem in (2.8) is convex. Further, it can easily be generalized
to unknown ν; we omit this discussion here for simplicity of presentation.

3
3.1

Statistical learning bounds
Decomposable Regularizers

Our learning bounds apply to general decomposable regularizers introduced in [46].
Given a subspace M ⊆ RM ×M K , we define its orthogonal complement as
⊥

M = {v ∈ RM ×M K | < u, v >= 0 for all u ∈ M}.
8

Given a normed vector space (RM ×M K , k · kR ) and subspaces M ⊆ M ⊆ RM ×M K ,
⊥
we say R is a decomposable regularizer with respect to (M, M ) if for A ∈ M
⊥
and B ∈ M we have
kA + BkR = kAkR + kBkR .
This definition encompasses widely-studied regularizers including the l1 norm,
nuclear norm, and the group sparsity inducing k · k1,2 norm. We refer the reader to
[46] for more details and intuition. While working in this general framework allows
us to incorporate a wide variety of prior beliefs about the structure of our network, a
fundamental question we need to address is how the specific choice of regularizer
effects our learning rates. Due to the temporal dependence and non-linearities in
SEPP models, deriving learning rates for various decomposable regularizers requires
us to leverage martingale concentration inequalities.

3.2

Assumptions

In Section II we presented a class of SEPPs which depends on a choice of basis
functions, and a general RMLE procedure which depends on a choice of regularization penalty. In this section, we introduce four assumptions which are needed
for our theoretical guarantees. We then give examples where we show that for
certain choices of basis functions {φ1 , . . . φK } and regularizers k.kR of interest the
assumptions hold with high probability.
Our first assumption depends on the basis functions but is not related to the
choice of regularizer.
Assumption 1 (Restricted Eigenvalue). There exists some ω > 0 and p ∈ N such
that smallest eigenvalue of E[g(Xt )g(Xt )> |Xt−p ] is lower bounded by ω for all t.
Assumption 1 is analogous to various restricted eigenvalue conditions in other
works. However, in much of the literature, one needs to lower bound the eigenvalues
of a sensing matrix whose columns are assumed to be independent. Dependence
introduced in our autoregressive model makes this a more complex condition to verify.
In past work on sparse autoregressive inference (e.g., [7]), restricted eigenvalue
conditions have been framed in terms of a stationary covariance matrix.
Informally, the value of ω measures the strength of the intertemporal dependence
of our process. If our network and basis functions are structured such that strong
long-range dependencies exist, then the smallest eigenvalue can be near zero, leading
to a poor bound on the error kAb − A∗ k2F .
The RE condition must also account for the level of clipping in our process: if
the network is so stimulatory that most observations are clipped, then the matrix
9

E[g(Xt )g(Xt )> |Xt−p ] will be nearly singular and ω will be close to zero. Thus, to
come up with an acceptable bound on ω, we need to establish that our network is
well-behaved enough that most observations will be unclipped. In other words, our
theory suggests that introducing non-linear saturation effects will not ruin our ability
to infer the structure of our network, provided that our network is not too stimulatory
and is usually stable without clipping. A further discussion of the intuition behind
the RE condition is provided in Example 1.
Next, we present assumptions which need to be verified in terms of the regularizer used. Recall that we assume the regularizer k · kR used in Equation (2.8) is
⊥
decomposable with respect to the pair of subspaces (M, M ).
Assumption 2 (Subspace Compatibility). There exists a constant Ψ(M) satisfying
kAkR
≤ Ψ(M).
A∈M kAkF
sup

Assumption 2 is a subspace compatibility condition as in [46], which controls
how large the Frobenius norm can be relative to the R norm on the subspace M.
Assumption 3 (Cone Row Sparsity). Let AM and AM⊥ denote the projections of a
⊥

matrix A onto the subspaces M and M respectively. Define
n

o

0
BR
= A ∈ RM ×M K : kAM⊥ kR ≤ 3kAM kR and kAkF = 1 .

Then there exists a constant µR such that
sup kBk22,1 ≤ µR .

0
B∈BR

Assumption 3 corresponds to assuming some notion of row sparsity on the error
matrix Ab − A∗ =: 4. It is needed to apply the empirical process techniques from
[28].
Assumption 4 (Deviation Bound). Let
t,m = Xt+1,m − exp(νm + a>
m. g(Xt )),
then there exists a constant λ < ∞ such that
T
1X
t g(Xt )>
T t=1

≤
R∗

λ
.
2

Assumption 4 is similar to deviation bound conditions found in the literature.
Due to the temporal dependence across observations, we must use martingale concentration inequalities under various norms in order to verify it.
10

3.3

General result

Provided our process and estimation procedure satisfy Assumptions 1-4 for reasonable constants, we can guarantee the learnability of our model.
Theorem 3.1. Assume (Xt )Tt=1 is generated by (2.3) and satisfies Assumptions 1-4
and assume A∗ is estimated according to (2.8) with a regularizer k · kR that is
⊥
decomposable with respect to the subspaces (M, M ). Then
kAb − A∗ k2F ≤
with probability at least 1 −

2
M2

36pΨ(M)λ2
2
Rmin
ω2

for

128p2 U 4 µ2R log M
T ≥
ω2
for constants C, c > 0 which are independent of M, T and Ψ(M).
Theorem 3.1 is a direct consequence of Theorem 1 in [46] combined with
Theorem 1 in [28]. Specifically, [46] gives Theorem 3.1 in a general decomposable
regularizer setup under a restricted strong convexity (RSC) assumption, which in
our language states that the error 4 =: Ab − A∗ satisfies
1 XX >
(4m. g(Xt ))2 ≥ kk4k2F
T t m
for some k > 0. Due to the fact that our process is neither linear nor Gaussian, many
techniques, e.g., [7] [44], used to establish an RSC condition directly are unworkable
in our setting. Instead, we use similar techniques to Theorem 1 in [28] which uses
empirical process results to turn the RSC assumption into the restricted eigenvalue
(RE) condition in Assumption 1.

4

Examples

In order to use Theorem 3.1, we need to prove that the four assumptions hold for
basis functions and regularizers of interest. First, we show that Assumption 1 is
satisfied for different point process models. Second we show Assumptions 2-4
are satisfied for a class of regularizers. Finally, we combine the results from this
section with Theorem 3.1 to give overall learning rates for ARMA(1, 1) and AR(2)
processes under different regularization schemes.
Recall that the constants Rmax and Rmin from (2.7) denote the maximum and
minimum Poisson rate for each observation.
11

4.1

Specific Point Process Models

Example 1: ARMA(1, 1) process First-order autoregressive moving average
(ARMA(1, 1)) point process models have been studied in a variety of settings
[37, 35]. Moreover, the corresponding continuous time model is one of the most
frequently studied point process models [4, 13]. Consider the following saturated
ARMA(1, 1) model with memory parameter α ∈ [0, 1):
Xt+1 ∼ Poisson(λt+1 )
log(λt+1 ) =ν + A∗ min(Xt , Ũ ) + α log(λt ).

(4.1)

Algebraic manipulation shows that (4.1) is a special case of (2.3) with K = 1 basis
function corresponding to φ[t] = αt . Here α is a memory parameter which captures
Ũ
the strength of the long-range dependence in the process, and kg(Xt )k∞ ≤ 1−α
=U
so that
Rmax

amax Ũ
= exp νmax +
1−α

Rmin

amin Ũ
= exp νmin −
1−α

!

!

.

An AR(1) process, corresponding to (4.1) where α = 0, was considered in [28].
However, due to the inherent instability of SEPPs without saturation, the authors
were forced to assume amax = 0.
Lemma 4.1. Suppose (Xt )Tt=1 is generated according to (4.1). Then Assumption 1
is satisfied with


1
ω = min Rmin , κ
2
where κ is a constant depending on Ũ , α and amax but independent of M .
The proof of Lemma 4.1 requires us to account for the effects of clipping. We
show that finding a lower bound on the eigenvalues of E[g(Xt )g(Xt )> |Xt−p ] can
be reduced to finding a lower bound on Var(min(Xt,m , Ũ )|Xt−1 ), which simplifies
the calculation since we rely on first-order dependence. Since we need to construct
a lower bound on Var(min(Xt,m , Ũ )|Xt−1 ) we consider the two cases when the
variance will be smallest.
In particular, if Xt,m ∼ Poisson(Rmin ), then its variance will be small because
the variance and mean of a Poisson random variable are equal. Specifically, when
Xt,m ∼ Poisson(Rmin ) we lower bound the variance of min(Xt,m , Ũ ) by 12 Rmin .
On the other hand, when Xt,m ∼ Poisson(Rmax ) and Rmax is large relative to Ũ ,
then min(Xt,m , Ũ ) is likely to be Ũ (clipped), so again the variance will be small.
12

We lower bound the variance by the constant κ that is the variance of a Bernoulli
random variable, where one outcome corresponds to a Poisson random variable
Z ∼ Poisson(Rmax ) exceeding Ũ (clipped) and the other outcome corresponds to
Z < Ũ .
One of these two worst case scenarios will give an absolute lower bound on the
variance. In both cases we construct a lower bound on the variance independent of
min Ũ
M . Note that ω increases with Rmin = exp(νmin − a1−α
) so ω grows inversely with
α. In other words, as the long range memory of the process increases, Lemma 4.1
suggests that network estimation becomes more difficult. This is consistent with
prior work [7].
The value κ−2 may be viewed as a proxy for the rate of clipping, and the
appearance of κ in Lemma 4.1 illustrates a tradeoff associated with clipped models.
On one hand, clipping ensures a stable process. However, as clipping increases, it
also increases the temporal dependencies among observations, leading to a smaller
κ and larger error bound.
In Figure 1, we fix α = 0 and get a sense of the value of κ for varying amax
and Ũ . (Recall that larger κ corresponds to a better-posed estimation problem.)
We see that for small Ũ , κ is not prohibitively small for a wide range of values of
amax . However, as Ũ increases the range of reasonable amax decreases, and as Ũ
approaches ∞, we approach the amax = 0 setting from [28]. To understand this
3
trend, we consider a special case of (4.1) with M = 1, ν = 0, α = 0, A∗ = 10
and
Ũ = 1000. In this case, the process follows:
Xt+1

3 min(Xt , 1000)
∼ Poisson exp
10

!!

.

Given a small X0 , this process will not be clipped for the first few observations,
but eventually the process will diverge and reach the clipping threshold of 1000.
This will happen within the first 100 observations with probability ≈ 1. Once an
observation reaches Xt ≥ 1000 and is clipped, the next observation will follow
Xt+1

3 min(Xt , 1000)
∼ Poisson exp
10
= Poisson(exp(300)).

!!

and so Xt+1 is virtually guaranteed to be clipped as well. In other words, once we
actually reach the clipping threshold, we enter a constant clipping regime which is
reflected in the small value of κ for amax = .3 and Ũ = 1000.
On the other hand, if Ũ = 6 our process follows
Xt+1

3 min(Xt , 6)
∼ Poisson exp
10
13

!!

.

Even when Xt ≥ 6, the Poisson rate exp( 18
) is approximately 6, and so the next
10
observation is reasonably likely to be unclipped and we do not enter the constant
clipping loop as when Ũ = 1000. In other words, over the long run we experience
less clipping for smaller Ũ , and thus κ is larger for smaller Ũ .

Figure 1: Colors indicate the value of κ for a given (amax , Ũ ) pair. Note that due to
our exponential link functions, elements of A∗ above one would be unreasonably
excitatory for many networks, and our process can be significantly stimulatory even
with coefficients well below one.

Example 2: AR(2) process As a second example, we consider an AR process
with two time lags:
Xt+1 ∼ Poisson(λt+1 )
log(λt+1 ) = ν +

A∗1

min(Xt , Ũ ) +

(4.2)
A∗2

min(Xt−1 , Ũ ).

(4.3)

This process fits within the framework of (2.3) with two basis functions corresponding to: φ1 [t] = I{t=1} and φ2 [t] = I{t=2} , |{z}
A∗ = [A∗1 , A∗2 ]. This example illustrates
M ×2M

the benefit of considering a basis with more than one element to describe the influence functions hm,m0 . A richer class of higher order models can be expressed with
multiple basis functions. Under this setup, the maximum and minimum possible
Poisson rates are
Rmax = exp(νmax + amax Ũ )
and
Rmin = exp(νmin − amin Ũ ).
14

Learning rates for high-dimensional linear AR(p) processes with Gaussian noise
were studied in [7]. However, the techniques used in that work to prove a restricted
eigenvalue condition relied heavily on the Gaussianity of the process. We prove that
the restricted eigenvalue condition in Assumption 1 holds for the AR(2) process in
Lemma 4.2.
In order to state Lemma 4.2 we first need several definitions. A node m is said to
be a parent of node m0 if it influences m0 through A∗1 , while m0 is said to be a child
of m. Furthermore, two nodes are said to be siblings if they share a parent node.
Lemma 4.2. Suppose (Xt )Tt=1 is generated according to (4.3). Let ρ(p)
m denote the
(s)
number of parents of m, let ρ(c)
denote
the
children
of
m
and
let
ρ
m
m denote the
number of siblings of m. Then
λmin (E[g(Xt )g(Xt )> |Xt−2 ]) ≥ rρ > 0
(c)
(s)
for a constant rρ depending on Rmax , Rmin , ρ(p)
m , ρm , ρm but independent of M .
(c)
(s)
The constant rρ scales inversely with ρ(p)
m , ρm , ρm and Rmax − Rmin . In the high
dimensional setting, this means a sparsity assumption on A∗ is necessary for our
bound to be useful.
We prove Lemma 4.2 by showing that the matrix Cov(g(Xt )|Xt−2 ) is strictly
diagonally dominant. A matrix B is said to be strictly diagonally dominant if there
P
exists a constant ω > 0 such that bi,i − j6=i |bi,j | ≥ ω for all i, and the eigenvalues
of a symmetric strictly diagonally dominant matrix are lower bounded by ω. With a
sparsity assumption on A∗ , almost all of the off diagonal terms in Cov(Xt ) will be
zero, and the remaining terms can be controlled with the techniques from Lemma
4.1 and appropriate assumptions on the size of Rmin and Rmax relative to the sparsity
(c) (s)
constants ρ(p)
m , ρm , ρm .

4.2

Regularization Examples

In this subsection, we verify Assumptions 2-4 under various regularization schemes.
Example 1: Element-wise Sparsity Regularization We first explore sparsity
regularization for these processes that accounts for the sparsity of A∗ natural to many
application domains. For the remainder of the section, we assume
kA∗ k0 = s  M 2 .
Sparse models of network structure encapsulate essential aspects of many common
statistical network models [24], and have connections to stochastic block models,
15

exponential random graph models, and various graphical models. We consider the
regularizer
XX
kAk1,1 =
|ai,j |
i

j

along with its dual
kAk∞,∞ = max max |ai,j |.
i

j

To see that k · k1,1 is decomposable we first define the set
S = {(i, j) ∈ {1, . . . , M } × {1, . . . , M K} : A∗i,j 6= 0},
and next define
S = {s ∈ RM ×M K : si,j = 0 for all (i, j) 6∈ S}.
Then k · k1,1 is decomposable with respect to the pair (S, S ⊥ ).
Note that the optimization problem corresponding to k · k1,1 regularization is
convex and can be solved with a variety of sparse regularization solvers. Furthermore,
it can trivially be parallelized across the rows of A.
Lemma 4.3. Suppose (Xt )Tt=1 is generated according to (2.3) with kA∗ k0 = s.
Further, assume A∗ is estimated according to (2.8) using k · k1,1 regularization. Then
(a) Assumption 2 is satisfied with
√
Ψ(S) = 4 s.
(b) Assumption 3 is satisfied with
√
µ(1,1) = 4 s.
(c) Assumption 4 is satisfied with
k

T
log3 (M T )
1X
t g(Xt )> k∞,∞ ≤ CRmax √
T t=1
T

with probability at least 1 −
dent of M, T and s.

1
(M T )c

for constants C, c > 0 which are indepen-

16

Proof Overview
• To verify Assumption 4, we rely on the fact that the Poisson rate can never
exceed Rmax . This allows us to bound the largest recorded observation by
C log(M T ) with high probability. From here, we are in a position to use martingale concentration inequalities developed in [33] to establish the deviation
bound.
• Assumptions 2 and 3 are straightforward consequences of the relation between
l1 and l2 norms.
Combining Lemma 4.3, Theorem 3.1 and the restricted eigenvalues results from
the previous subsection gives overall bounds for sparse SEPPs which are applicable
in the high-dimensional setting.
Theorem 4.4. (Learning rates for l1 regularization) Suppose (Xt )Tt=1 follows the
SEPP framework of (2.3) and A∗ is estimated using sparsity regularization.
(a) If (Xt )Tt=1 is generated according to the ARMA(1, 1) model in (4.1) then
kAb − A∗ k2

F

2
Rmax
s log6 (M T )
≤C 2
T
Rmin min( 21 Rmin , κ)2

with probability at least 1 −

1
(M T )c

for T, M satisfying

T ≥ 128U 4 s

log M
min( 12 Rmin , κ)2

for constants C, c > 0 which are independent of M, T and s.
(b) If (Xt )Tt=1 is generated according to the AR(2) model in (4.3) then
kAb − A∗ k2

F

with probability at least 1 −

2
Rmax
s log6 (M T )
≤C 2 2
Rmin rρ
T

1
(M T )c

for T, M satisfying

T ≥ 128U 4 s
6

log M
.
rρ2

The mean-squared error bound s log T(M T ) matches the minimax optimal rate in
the independent case [52] up to log factors. Theorem 4.4 extends results in Hall et
al. [28] to ARMA(1,1) and AR(2) processes.
17

Example 2: Group Sparsity Group lasso regularization is a popular tool used to
estimate a sparse parameter where one has prior knowledge on the structure of the
sparsity (see e.g., [62] for more details). We consider a special case of group lasso
regularization where the groups are the columns of the matrix. Let a.m denote the
mth column vector of a matrix A. Our structured sparsity assumption is that only
sG  M columns of A∗ contain nonzero entries.
In terms of network structure, this means that only a small number of hub nodes
have influence on other nodes in the network. To estimate networks of this form, we
consider l2 penalization on the columns vectors, followed by l1 penalization on the
resulting l2 norms. In other words, we have
kAkG = kA> k1,2 =

X

ka.m k2 .

m

The dual of this norm is
kAkG∗ = kA> k∞,2 = max
ka.m k2 .
m
Let
SG = {i : a∗.i 6= 0}.
Then k · kG is decomposable with respect to the subspaces
SG = {A : a.j = 0 for all j 6∈ SG }
and
SG⊥ = {A : a.j = 0 for all j ∈ SG }.
We show Assumptions 2-4 hold in Lemma 4.5 below.
Lemma 4.5. Suppose (Xt )Tt=1 is generated according to (2.3) where only sG
columns of A∗ contain nonzero entries. Further, assume A∗ is estimated according
to (2.8) using k · kG regularization. Then
(a) Assumption 2 is satisfied with
√
Ψ(SG ) = 4 sG .
(b) Assumption 3 is satisfied with
µG = 16sG .
(c) Assumption 4 is satisfied with
s

T
1X
M
k
t g(Xt )> kG∗ ≤ CRmax log2 (M T )
T t=1
T

with probability at least 1 −

1
.
(M T )c

18

Proof Overview
• For Assumption 4 we construct a high probability bound on the l2 norm of
each column of our noise matrix and take a union bound over all the columns
to get a final bound on the k · k∞,2 norm. To bound the norm of each individual
column, we rely on [51] which provides martingale concentration inequalities
for 2-smooth norms.
• For Assumption 3 we derive an error row sparsity constant which depends only
on sG rather than M . The k · k2,1 norm can be large relative to the Frobenius
norm in cases where the matrix
√ is row-dense. In this case, the l1 norm of each
row can be on the order of M larger than the l2 norm. However, we only
need to derive a compatibility constant on the cone
BG = {B ∈ RM ×M K : kBS

⊥
G

kG ≤ 3kBSG kG }.

Since elements of SG has at most sG nonzero entries in each row, we can think
of all matrices in the cone
√ BG as being “almost row sparse” and so the k · k2,1
norm should not be O( M ) larger than the Frobenius norm on the cone.
• Assumption 2 follows from the relationship between the l1 and l2 norms.
Combining Lemma 4.5, Theorem 3.1 and the restricted eigenvalue conditions from
the previous subsection gives the following result.
Theorem 4.6. (Learning rates for group lasso regularization)
Suppose (Xt )Tt=1 follows the SEPP framework of (2.3) and A∗ is estimated using
column group lasso regularization.
(a) If (Xt )Tt=1 is generated according to the ARMA(1, 1) model in (4.1) then
kAb − A∗ k2F ≤

C
sG M log4 (M T )
2
T
Rmin
min( 21 Rmin , κ)2

with probability at least 1 −

1
(M T )c

for T, M satisfying

T ≥ 128U 4 s2G

log M
.
min( 12 Rmin , κ)2

for constants C, c > 0 which are independent of M, T and sG .

19

(b) If (Xt )Tt=1 is generated according to the AR(2) model in (4.3) then
sG M log4 (M T )
2
Rmin
rρ2
T
C

kAb − A∗ k2F ≤
with probability at least 1 −

1
(M T )c

for T, M satisfying

T ≥ 128U 4 s2G

log M
.
rρ2

Example 3: Low-rank Regularization Estimation of high-dimensional but lowrank matrices is a widely studied problem with numerous applications [44, 53, 64, 36,
19]. Low-rank models can be seen as a generalization of sparse models, where the
matrix is sparse in an unknown basis. A standard technique to estimate a low-rank
matrix is to take a convex relaxation of an l0 penalty on the singular values [19]: the
nuclear norm penalty
kAk∗ =

M
X

σi (A),

i=1

where σi (A) denotes the ith singular value of A. The dual to the nuclear norm is the
operator norm
kAkop = sup kAxk2 .
kxk2 =1

As discussed in [46], the nuclear norm is decomposable with respect to the subspaces
W = {A ∈ RM ×M K : row(A) ⊆ row(A∗ ) and col(A) ⊆ col(A∗ )}
and
⊥

W = {A ∈ RM ×M K : row(A) ⊆ row(A∗ )⊥ and col(A) ⊆ col(A∗ )⊥ },
where row(A) and col(A) denote the row and column spaces of A respectively.
Unlike the previous two examples, here W =
6 W.
In this low-rank setup, there is no limitation on the number of nodes which can
influence a given node. This introduces challenges in establishing Assumption 3,
which guarantees near row sparsity of the error. In order to get around this, we
impose a technical assumption on kA∗ k2,1 in Lemma 4.7. An area of interest for
future work is to examine whether our estimation procedure is flawed when one
node can have many nodes influence it, or whether the need for this assumption is
an artifact of our analysis.

20

Lemma 4.7. Suppose
(Xt )Tt=1 is generated according to (2.3) with rank(A∗ ) = r
√
and kA∗ k22,1 = D M for a universal constant D. Further, assume A∗ is esti√
mated according to (2.8) over the ball {A : kAk22,1 ≤ D M } using nuclear norm
regularization. Then
(a) Assumption 2 is satisfied with
Ψ(W) =

√
2r.

(b) Assumption 3 is satisfied with
√
µ∗ = 2D M .
(c) Assumption 4 is satisfied with
s

T
1X
M
k
t g(Xt )> kop ≤ log4 (M T )
T t=1
T

with probability at least 1 −

1
.
(M T )c

Proof Overview
• The main challenge in Lemma 4.7 comes in the proof of the deviation bound
condition, which depends on the concentration properties of vector-valued
martingales. The concentration properties of 2-smooth norms was studied in a
number of works, including [21, 22]. We leverage recent work in [34] which
extends the concentration results for 2-smooth norms to operator norms.
b
• Assumption
3 follows from assuming kA∗ k2,1 and kAk
2,1 are on the order of
√
M . Without this assumption, we could potentially have µ∗ = O(M ). This
would give us a final bound in Corollary 4.8 which is only applicable when
T ≥ M 2 , so this technical assumption is crucial in constructing a meaningful
bound.

• The subspace compatibility constant in Assumption 2 was shown in [53]. In
the sparsity case, this condition is trivial because S = S and thus 4S is known
to lie in a subspace where every element is s-sparse. The condition is more
subtle in the nuclear norm regularization case because W = W if and only if
A∗ is symmetric. We do not assume symmetry of A∗ so 4W need not lie in
the subspace W where each element has rank at most r. However, [53] shows
that W only contains matrices of rank at most 2r.
21

Combining Lemma 4.7, Theorem 3.1 and the restricted eigenvalues results gives
the following Theorem.
Theorem 4.8. (Learning rates for nuclear norm regularization)
∗
Suppose (Xt )Tt=1 follows the SEPP framework of (2.3) and
√ A is estimated using
2
nuclear norm regularization over the ball {A : kAk2,1 ≤ D M }.
(a) If (Xt )Tt=1 is generated according to the ARMA(1, 1) model in (4.1) then
kAb − A∗ k2F ≤

C
rM log8 (M T )
2
T
Rmin
min( 21 Rmin , κ)2

with probability at least 1 −

1
(M T )c

for T, M satisfying

T ≥ 128U 4 M

log M
min( 21 Rmin , κ)2

for constants C, c > 0 which are independent of M, T and r.
(b) If (Xt )Tt=1 is generated according to the AR(2) model in (4.3) then
kAb − A∗ k2

F

with probability at least 1 −

rM log8 (M T )
≤ 2 2
Rmin rρ
T
C

1
(M T )c

for T, M satisfying

T ≥ 128U 4 M

log M
.
rρ2
8

Once again the mean-squared error bound rM logT (M T ) matches the minimax
optimal rate for independent design [45] up to log factors.

5

Numerical experiments

We validate our methodology and theory using a simulation study and real data
examples. The focus of the simulation study is to confirm that the rates in meansquared error in terms of s, r, T and κ scale as the theory predicts. We generate data
according to the ARMA(1, 1) model from 4.1.
Our focus with real data experiments is to demonstrate that the models we
analyze are sufficiently complex to capture real-world phenomena and enhance
prediction performance relative to naive models. Others have successfully used more
22

complex, difficult to analyze models (cf., [39, 17]) which are similar in spirit to those
analyzed here. Our claim is not that our approach leads to uniformly better empirical
performance than previous methods, but rather that our models capture essential
elements of all these approaches and hence our theoretical work provides insights
into a variety of approaches.
Our first real data example shows that our model and estimation procedure
determines interactions among shooting events across different communities of
Chicago that obeys sensible spatial structure (even though the algorithm does not use
any spatial information). Our second real data example looks at neuron firing patterns
in the brain of a rat and shows that our model can differentiate between the firing
patterns during a sleep period and the patterns during a wake period. Finally, we
examine a data set consisting of articles posted by different news websites and we try
to determine influences between the sites using a variety of different regularization
techniques. We implement the estimation procedure in (2.8) using the SpaRSA
algorithm from [61].

5.1

Simulation study

We generate data according to (4.1) with ν = 0, M = 50, Ũ = 6, α = .25 and
varying values for T and s. Recall that ν controls the background rate, M is the
number of nodes, Ũ is the clipping threshold, α is the memory of the process in
(4.1), T is the number of time steps observed, and s is the number of edges in the
network.
Each time we generate a matrix, we randomly select s entries to be nonzero, and
assign each value uniformly in [−.7, .3]. The sparsity ranges between 10 and 50.
With these parameters, our process is usually stable on its own, and only occasionally
relies on the clipping function to dampen the observations. Even at s = 50, only
around 1% of the observations exceed 6, and the clipping percentage is
√even lower
for smaller s. For each choice of s, T , we run 100 trials with λ = .1/ T . In the i
trial we form a ground truth matrix A∗i , compute Âi , and measure the mean squared
error (MSE) as kA∗i − Âi k2F .
In Figure 2(a) , we plot MSE vs T for several different values of s, and in
Figure 2(b) we plot MSE vs s for several values of T . The plots agree with our
theory, which suggests that the error scales linearly in s and T1 .
Next, instead of assuming that A∗ only has s non-zero entries, we assume that A∗
has rank r, and measure MSE as a function of r. We hold the remaining parameters
the same. To generate a rank r matrix A∗ , we randomly generate a M × r matrix
and multiply it by a randomly generated r × M matrix where the entries of both
matrices are uniformly drawn from [−.7, .3]. We then normalize the resulting matrix
so that amax is approximately .3. For all choices of rank considered, less than 5% of
23

(a) MSE vs T

(b) MSE vs s

Figure 2: (a) shows MSE vs T for varying values of s, while (b) shows MSE vs s for
varying values of T . Plots agree with theory which suggests that MSE scales linearly
in s and T1 . Median of 100 trials are shown, and error bars denote the standard
deviation.
q

observations are clipped. We set λ = .1 M
as guided by our theory, run 100 trials
T
for each (r, T ) pair, and plot the median in Figure 3. The simulations agree with
Lemma 4.7 which suggests that the MSE should scale linearly in r and T1 .

(a) MSE vs T

(b) MSE vs r

Figure 3: (a) shows MSE vs T for varying values of r, while (b) shows MSE vs r
for varying values of T . Plots agree with theory which suggests that MSE scales
linearly in r and T1 . Median of 100 trials are shown.
We now examine the relationship between amax , our theoretical MSE, and simulated MSE. For the remainder of the section fix Ũ = 6 and α = 0 and T = 400.
Recall from Figure 1 that there is a stark phase transition where κ goes from reasonably large to minuscule. For Ũ = 6, this transition occurs between amax = .3 and
amax = .4. Our MSE bound scales with κ−2 , so this phase transition controls where
our bound is reasonably small.
24

To get a sense of how tight the bound is, we consider two different methods to
generate a 50 × 50 matrix A∗ . The first is a block design, where A∗ is zero outside
of five 10 × 10 blocks on the diagonal. Within the blocks, each row has five nonzero
entries picked at random with values equal to amax
. Matrices with this structure
5
have strong feedback loops, where large observations from one node stimulate other
nodes which are likely to feedback to the original node. In other words, with this
block design method it is likely that many observations will actually be drawn close
to the maximum possible rate Rmax , so we expect the MSE to align closely with
our theoretical bound. We estimate A∗ using l1 regularization, for varying values of
amax .
As a second method, we consider a low-rank design and estimate A∗ using
nuclear norm regularization. We choose the first two rows of A∗ to be orthogonal,
both with row sums equal to amax . We then let each remaining row be a random
convex combination of the first two rows. In this case, feedback loops are less of
an issue; if one node has a large observation, the nodes it stimulates are less likely
to have strong connections feeding back to the original node. Since this design
method is not particularly stimulatory, we expect that most observations will not be
drawn close to the maximum rate Rmax . Our theoretical bound is potentially loose
in this case for the following reason. Our bound on κ which captures the amount of
clipping is worst case based on the size of the coefficients of A∗ , but does not take
into account the structure of A∗ . If the coefficients of A∗ are large but are structured
such that there aren’t a lot of feedback loops then our bound on κ will be loose.
We randomly generate 50 different A∗ over various amax for both design choices
and then evaluate their efficiency by plotting the fraction of trials for which the MSE
is above one. The results are shown in Figure 4. The simulated results also exhibit
strong phase transitions, with the fraction of accurate trials shifting from one to zero
with small changes in amax . This suggests that our theoretical results capture a real
phenomenon of our model. In the block case, the phase transition occurs almost
exactly where predicted by the MSE, whereas in the low-rank case there is a small
lag in the phase transition. In other words, while our theoretical results are fairly
tight for very stimulatory network structures, there appears to be some flexibility for
networks with weaker feedback loops.
Finally, we consider the block matrix design under a wide range of Ũ . We
consider Ũ between 3 and 30, and amax between 0 and .6 in increments of .02. For
each (amax , Ũ ) pair, we generate 20 matrices according to the block matrix design
strategy outlined above and estimate Ab via sparsity regularization. In Figure 5 we
plot a heat map displaying the fraction of trials for with the MSE is below one. The
red line shows the κ = .01 contour, so Lemma 4.3 suggests our model will be hard to
learn above the boundary line. Figure 5 generally resembles the heat map displaying
values of κ in Figure 1, suggesting that the role of κ in our theory reflects a true
25

Figure 4: In blue, κ−2 as a function of amax . In orange, the fraction of 50 trials
which have MSE above 1 for a block matrix design strategy and a low-rank matrix
design strategy. Our theoretical bound on the MSE scales with κ−2 , suggesting stark
phase transitions in the MSE which are confirmed in the simulated results. The
block matrices correspond to more stimulatory networks than the low-rank ones, and
therefore align more closely with our worst case bounds.
phenomenon that when amax is sufficiently large and clipping is frequent, then the
model becomes difficult to learn.

5.2

Real Data Example – Chicago Crime

A number of studies have used various self-exciting point processes to predict crime,
including [57, 43, 42]. We test our model on a data set [2] consisting of burglaries
in Chicago since 2004, broken down by the M = 77 community areas in the city.
In [3], the authors fit self-exciting processes to the Chicago homicide data broken
down by community area and performed clustering on the areas as we do below.
We estimate the network based on the data from January 2004 to August 2010 and
test it on the data from September 2010 to March 2017. To test results, we compare
the log-likelihood of events using our learned matrix on the test set data, with that
for a constant Poisson process. This gives√approximately 600 time periods for both
our training and test sets. We set λ = .1/ T using our theory as a guide. We show
results for a half-day time discretization period, with Ũ = 7 and α = .2.
The test set log-likelihood of our learned matrix shows an improvement over
the test set log-likelihood for a learned constant process, where λt,m = λm for all t
(−6.62 × 105 compared to −1.09 × 106 ).
26

Figure 5: Heat map displaying the fraction of trials for which the MSE is below
one for different (amax , Ũ ) pairs. Red line shows κ = .01 contour. Above this line
our theory predicts difficulty in learning, agreeing with the heat map which shows
inaccurate recovery of A∗ for (amax , Ũ ) pairs above the line. For each (amax , Ũ ) pair,
we run 20 trials with a block matrix design.
To examine the structure of our learned matrix, we treat the positive coefficients
of the matrix as edges of the adjacency matrix of a graph. We then perform spectral
clustering with four clusters. The results are shown in Figure 6, with colors indicating
cluster membership. We note that our data contains no information about the
geospatial location of the areas aside from index (not location) of the community
area. However, there are clear geographic patterns in the clusters, providing some
validation to the estimated influences between communities.
Finally, we test whether modeling these crime patterns as a multi-dimensional
point processes leads to stronger results than modeling the patterns as a collection of
independent univariate point processes, where
log λt+1,m = νm + a∗m min(Xt,m , Ũ ) + α log(λt,m ).
Specifically, we compare finding Ab as in (2.8) with choosing Ab to be the solution
to the optimization problem in (2.8) over the set of all diagonal matrices. We then
perform the log-likelihood analysis described above for α varying from 0 to .6 in
increments of .2 and for the time-discretization period varying from half a day to
three days in increments of half a day. The results are shown in Table 1. Note that
the multivariate model outperforms the univariate one whenever the discretization
period is at least a full day but does worse for the half day discretization period.

27

Figure 6: Clusters learned from crime data with a half-day time discretization. The
clusters are overlaid on a map of community areas in Chicago. The data contained
no geospatial information, but clusters show geographical patterns.

5.3

Real Data Example – Spike Train Data

SEPPs have been widely used in neuroscience to describe neuron spike train data
[38, 60, 41, 40]. In this section, we analyze a multi-neuron spike train dataset from
[58, 59]. The dataset consists of spike trains recorded from 51 neurons in the brain
of a rat. The recordings were divided into a wake period and a sleep period. Using
the data from the first half of the wake period, we learn a matrix Awake using equation
(4.1), a 100ms discretization period, α = .7 and Ũ = 5. We then follow the same
process to learn Asleep . We get a sense of the structure of the matrices Awake and
Asleep in Figure 7. We note that connections are much stronger during the wake
period, during which there is more frequent neural firing.
In previous work [35] the authors use a similar SEPP to analyze neural spikes
and discuss the significance of the time discretization period in more depth. In particular, they conclude that while models at this discretization length may have strong

28

Table 1: Difference between Log-Likelihood of multivariate process and univariate
process
Discretization
α=0
.5 days
−5.4 × 104
1 day
1.2 × 105
1.5 days
1.8 × 105
2 days
1.8 × 105
2.5 days
2.1 × 105
3 days
1.9 × 105

α = .2
−6.5 × 104
7.6 × 104
1.1 × 105
1.2 × 105
1.1 × 105
1.1 × 105

α = .4
−6.6 × 104
6.2 × 104
1.4 × 105
2.3 × 105
2.7 × 105
2.8 × 105

predictive power, the discretization period is sufficiently large that the connections
learned are not direct physical effects. In other words, if the connection between
neuron A and neuron B is negative, this suggests that neuron B is less likely to fire
in a 100ms interval after neuron A fires. However, there could be a complex chain
of interactions causing this effect, and it does not mean there is a direct physical
connection between neuron A and neuron B.
To validate Awake , we compute the log-likelihood of events for the second half
of the wake period using both Awake and Asleep as the ground truth matrix. We find
log p(Xwake |Âwake ) = −6.6 × 104 while log p(Xsleep |Âwake ) = −7.4 × 104 .
Following the same process for Asleep , we find that log p(Xsleep |Âsleep ) = −2.34×
105 while log p(Xwake |Âsleep ) = −3.04×105 . This suggests that our model is capable
of differentiating firing patterns in different sleep states. The log-likelihood of events
for a constant process was orders of magnitude smaller in both cases.

5.4

Real Data Example – Memetracker Data

As a final example, we consider a data set [1] which consists of metadata for a
collection news articles and blog posts. We only consider the time and website from
which each post occurs but omit all other data such as the content of the post and other
websites to which the post links. Further, we consider only articles posted by 198
popular news sites from http://www.memetracker.org/lag.html. Lowrank models have been applied in social network settings in a number of different
works [65, 66, 63]; in particular, the work [65] proposes low-rank regularization of
a point process model on this same data set.
To test the model, we collect all articles posted by 198 popular new websites
during October 2008. Using a one hour discretization period, we divide the month
into a training set and a test set, giving T = 500 training periods and 500 test periods.

29

1

1

5

5
0.5

10

0.5

10
0

0

15

15
-0.5

20
25

-0.5

20
25

-1

30

-1

30
-1.5

35

-1.5

35
-2

40
45

-2.5

50

-3

5

10

15

20

25

30

35

40

45

-2

40
45

-2.5

50

50

-3

5

10

15

(a) Awake

20

25

30

35

40

45

50

(b) Asleep

Figure 7: (a) shows Awake matrix charting estimated relationships between neurons in
a rat’s brain during a wake period, while (b) shows Asleep matrix charting estimated
relationships between neurons in a rat’s brain during a sleep period.
We train our model using the following regularization techniques. We perform the l1
regularization and nuclear norm regularization schemes described in Section 4.2, as
well as a low-rank plus sparse model where we use the regularizer
k · kR = k · k1,1 + k · k∗ .
This last model is optimized using alternating descent. Finally, we learn a multidimensional model with no regularization, where we simply use the negative loglikelihood as our loss function, and a one-dimensional model where all interactions
between different nodes are set to zero. The results are shown in Figure 8. The lowrank model performs best, followed by the low-rank plus sparse model, suggesting
that the interactions between websites exhibit some low-rank behavior.

6

Connections to Hawkes Process

In this section, we observe that the model in (2.3) can be seen as a discretized version
of the multivariate Hawkes process, in which there is much long-standing and recent
interest (e.g., [30, 31, 13, 11, 65, 48, 25]). By formulating our discrete time model in
this manner we aim to highlight the connections between the two classes of models.
There are advantages to analyzing point process models from both the continuous
and discrete perspective. Some advantages of the discrete perspective include:
1. Real world data comes inherently discretized. In some cases, e.g., social media
posts, one might record data accurately up to very fine time windows. However,
30

Figure 8: Log-likelihood of events on test set for matricies learned using memetracker
data set under a variety of structure assumptions. Data set consists of timestamps for
articles posted by 198 popular news websites and blogs during October 2008.
in other problems, the data collection process forces a coarse discretization.
For example in [37], which used Hawkes processes to model civilian deaths in
Iraq, reliable data was only obtained for the day on which attacks occurred.
2. In many works the authors discuss continuous Hawkes models, but their
algorithms work with discretized data for computational efficiency. Examples
include [4, 37, 41]. This provides additional motivation for our decision to
study the ARMA(1, 1) model because its continuous version is one of the
most widely studied types of Hawkes process.

6.1

Continuous Hawkes

In a multivariate Hawkes process, point process observations are drawn using an
intensity function λ(Xτ ), where Xτ is the collection of all events up to (continuous)
time τ . Each event i is associated with two components: (τi , mi ), where τi is the
time of the event and mi is the node or channel associated with the event. Nτ denotes
the number of events before time τ .
We model the log-linear Hawkes process intensity at node m as1
log λ(c)
m (Xτ ) =νm +

Nτ
X

hm,mi (τ − τi ),

(6.1)

i=1

where the (c) superscript denotes continuous time.
1 (c)
λm (Xτ )

(c)

would be more precisely written as λm (τ ; Xτ ); we let the dependence on τ be
understood for simplicity of presentation.

31

Here each function hm,m0 (τ ) measures the influence of node m0 on node m
after τ seconds since the event on m0 . This model is standard in the point process
literature. We write each of these functions as a linear combination of the basis
functions φ1 (τ ), . . . , φK (τ ):
K
X

hm,m0 (τ ) =

a∗m,m0 ,k φk (τ ),

(6.2)

k=1

yielding
log λ(c)
m (τ ) =νm +

Nτ X
K
X

a∗m,mi ,k φk (τ − τi )

i=1 k=1



=νm +

M
X
m0 =1

=νm +

K
 X X

a∗m,m0 


i<Nτ : k=1
mi =m0

M X
K
X
m0 =1




φk (τ − τi )


(c)

a∗m,m0 ,k gm,k (Xτ ),

(6.3)

k=1

where
(c)

gm,k (Xτ ) :=

X

φk (τ − τi ).

i<Nτ :
mi =m

Vectorizing across nodes and letting
(c)

(c)

λ(c) (Xτ ) := [λ1 (Xτ ), · · · , λM (Xτ )]> ∈ RM
+
ν := [ν1 , · · · , νM ]> ∈ RM
(c)

g (c) (Xτ ) := (gm,k (Xτ ))m∈{1,...,M },k∈{1,...,K} ∈ RM K
A∗ := (am,m0 ,k )m,m0 ∈{1,...,M },k∈{1,...,K} ∈ RM ×M K ,
we have
log λ(c) (Xτ ) = ν + A∗ g (c) (Xτ )
which exhibits the same general form as (2.3).
In order to formalize the connection between the multivariate Hawkes process in
(6.3) and the SEPP in (2.3) we first describe our sampling process and the Hawkes
and Poisson log likelihoods needed to prove Proposition 1: The Hawkes process can
be discretized by sampling λ(c) (Xτ ) at τ = t4 for some sampling period 4 > 0
and letting
Nt4

Xt,m =

X
i=N(t−1)4 +1

32

Im=mi

(6.4)

for t = 1, . . . , T .
Here IE is the indicator function which returns 1 if E is true and 0 if E is
false and Xt,m is the number of events on node m during the sampling interval
[(t − 1)4, t4). Overloading notation somewhat, let Xt = {Xs,m }s=1,...,tm=1,...,M be
the history of event counts up to time t. The log-likelihood of the original Hawkes
process observations up to time T 4 is
NT 4

`H (XT 4 |{λ(c)
m }m )

=

X

log λ(c)
mi (τi )

i=1

=

4m.
M NT
X
X

m=1

−

M Z T4
X
m=1 0

log λ(c)
mi (τi ) −

i=1

λ(c)
m (τ )dτ

Z T4
0

(c)
λm
(τ )dτ.

Further note that if Xt,m ∼ Poisson(λt,m ), then the Poisson log likelihood is proportional to
`p (XT )|{λm (Xt )}t,m ) =

M X
T
X

[Xt,m log(λm (Xt )) − λm (Xt )].

m=1 t=1

We consider in Proposition 1 a SEPP with the intensity
(c)
λt,m = 4λ(c)
m (X4t ) ≡ 4λm (4t; X4t ),

(6.5)

where the last equality makes the sampling time explicit. We now present a proposition which formalizes the connections between the SEPP and the log-linear Hawkes
process.
Proposition 1. The likelihood of the discretized multivariate Hawkes data in (6.4)
can be approximated by the likelihood of the Poisson autoregressive model (2.3)
with intensity (6.5), modulo terms independent of the unknown λ(c) , where the
approximation error depends on the sampling period 4.
This proposition suggests that the models and analysis we develop for SEPPs
also provides insight into related Hawkes process models provided that the sampling
period ∆ is sufficiently small.

7

Conclusion

The proposed saturated SEPP allows us to analyze statistical learning rates for a
large class of point processes, including discretized Hawkes processes, with longrange memory and saturation or clipping effects common in real-world systems.
33

The analysis presented in this paper addresses instability issues present in prior
works and incorporates a wide variety of structural assumptions on the ground truth
processes by allowing for arbitrary decomposable regularizers. The proposed bounds
provide novel insight not only into sample complexity bounds, but also into phase
transition boundaries dictated by stability and saturation effects that are supported by
simulation results. In addition, experiments on data from neuroscience, criminology,
and social media suggest that the models considered in this paper exhibit sufficient
complexity to model real-world phenomena.

References
[1] Memetracker. 2008. http://www.memetracker.org/data.html.
[2] City of chicago. 2017. https://data.cityofchicago.org/Public-Safety/Crimes2001-topresent/ijzp-q8t2.
[3] Ryan Adams and Scott. Linderman. Discovering latent network structure in
point process data. In Proc. International Conference on Machine Learning
(ICML), 2014.
[4] Y. Aït-Sahalia, J. Cacho-Diaz, and R. J. A. Laeven. Modeling financial contagion using mutually exciting jump processes. Technical report, National
Bureau of Economic Research, 2010.
[5] K. Azuma. Weighted sums of certain dependent random variables. Tohoku
Mathematical Journal, 19(3):357–367, 1967.
[6] E. Bacry and J.-F. Muzy. First- and second-order statistics characterization
of hawkes processes and non-parametric estimation. IEEE Transactions on
Information Theory, 62:2184–2202, 2016.
[7] S. Basu and G. Michailidis. Regularized estimation in sparse high-dimensional
time series models. Annals of Statistics, 43(4):1535–1567, 2015.
[8] C. Blundell, K. A. Heller, and J. M. Beck. Modelling reciprocating relationships
with Hawkes processes. In Proc. NIPS, 2012.
[9] E. N. Brown, R. E. Kass, and P. P. Mitra. Multiple neural spike train data
analysis: state-of-the-art and future challenges. Nature neuroscience, 7(5):456–
461, 2004.

34

[10] V. Chavez-Demoulin and J. A. McGill. High-frequency financial data modeling
using Hawkes processes. Journal of Banking & Finance, 36(12):3415–3426,
2012.
[11] S. Chen, A. Shojaie, E. Shea-Brown, and D. Witten. The multivariate
Hawkes process in high dimensions: Beyond mutual excitation. 2017.
https://arxiv.org/abs/1707.04928.
[12] T. P. Coleman and S. Sarma. Using convex optimization for nonparametric
statistical analysis of point processes. In Proc. ISIT, 2007.
[13] D. J. Daley and D. Vere-Jones. An introduction to the theory of point processes,
Vol. I: Probability and its Applications. Springer-Verlag, New York, second
edition, 2003.
[14] M. Ding, CE Schroeder, and X. Wen. Analyzing coherent brain networks with
Granger causality. In Conf. Proc. IEEE Eng. Med. Biol. Soc., pages 5916–8,
2011.
[15] J. Muzy E. Bacry, S. Gaïffas. A generalization error bound for sparse and
low-rank multivariate hawkes processes. arXiv:1501.00725, 2015.
[16] M. Egesdal, C. Fathauer, K. Louie, J. Neuman, G. Mohler, and E. Lewis.
Statistical and stochastic modeling of gang rivalries in Los Angeles. SIAM
Undergraduate Research Online, 3:72–94, 2010.
[17] Ş. Ertekin, C. Rudin, T. H. McCormick, et al. Reactive point processes: A new
approach to predicting power failures in underground electrical systems. The
Annals of Applied Statistics, 9(1):122–144, 2015.
[18] J. Etesami, N. Kiyavash, K. Zhang, and K. Singhal. Learning network of
multivariate Hawkes processes: A time series approach. arXiv preprint
arXiv:1603.04319, 2016.
[19] M Fazel.
Matrix rank minimization with applications.
http://faculty.washington.edu/mfazel/thesis-final.pdf.

2002.

[20] K. Fokianos and D. Tjøstheim. Log-linear poisson autoregression. Journal of
Multivariate Analysis, 102(3):563–578, 2011.
[21] D.H. Garling. Functional central limit theorems in banach spaces. Annals of
Probability, 4(4):600–611, 1976.

35

[22] D.H. Garling. Convexity, smoothness and martingale inequalities. Israel
Journal of Mathematics, 29:189–198, 1978.
[23] F. Gerhard, M. Deger, and W. Trucculo. On the stability and dynamics of
stochastic spiking neuron models: Nonlinear Hawkes process and point process
glms. PLos Computation Biology, 13(2):1069–1097, 2017.
[24] A. Goldenberg, A. X. Zheng, S. E. Fienberg, and E. M. Airoldi. A survey of
statistical network models. Foundations and Trends R in Machine Learning,
2(2):129–233, 2010.
[25] H. Zha H. Xu, M. Farajtabar. Learning granger causality for hawkes processes.
Proceedings of Machine Learning Research, 48:1717–1726, 2016.
[26] E. Hall and R. Willett. Dynamical models and tracking regret in online convex
programming. In Proc. International Conference on Machine Learning (ICML),
2013. http://arxiv.org/abs/1301.1254.
[27] E. Hall and R. Willett. Online learning of neural network structure from spike
trains. In Proceedings of the 7th International IEEE EMBS Neural Engineering
Conference (NER’15), 2015.
[28] E. C. Hall, G. Raskutti, and R. Willett. Inference of high-dimensional autoregressive generalized linear models. arXiv preprint arXiv:1605.02693, 2016.
[29] N. R. Hansen, P. Reynaud-Bouret, and V. Rivoirard. LASSO and probabilistic
inequalities for multivariate point processes. arXiv preprint arXiv:1208.0570,
2012.
[30] A. G. Hawkes. Point spectra of some self-exciting and mutually-exciting point
processes. Journal of the Royal Statistical Society. Series B (Methodological),
58:83–90, 1971.
[31] A. G. Hawkes. Point spectra of some mutually-exciting point processes. Journal of the Royal Statistical Society. Series B (Methodological), 33(3):438–443,
1971.
[32] M. Hinne, T. Heskes, and M. A. J. van Gerven. Bayesian inference of wholebrain networks. arXiv:1202.1696 [q-bio.NC], 2012.
[33] Christian Houdré and Patricia Reynaud-Bouret. Exponential inequalities, with
constants, for u-statistics of order two. In Stochastic inequalities and applications, pages 55–69. Springer, 2003.
36

[34] A. Juditsky and A. Nemirovski. Large deviations of vector-valued martingales
in 2-smooth normed spaces. 2008. https://arxiv.org/abs/0809.0813.
[35] R. Kelly, R. Kass, Ml Smith, and T. Lee. Accounting for network effects in neuronal responses using l1 regularized point process models. NIPS, 23(2):1099–
1107, 2010.
[36] K. Lee and Y Bresler. Guaranteed minimum rank approximation from linear
observations by nuclear norm minimization with an ellipsoidal constraint. 2009.
https://arxiv.org/abs/1602.07389.
[37] E. Lewis, E. Mohler, P. J. Brantingham, and A. L. Bertozzi. Self-exciting point
process models of civilian deaths in iraq. Security Journal, 25(3):244–264,
2012.
[38] S. Linderman, R. Adams, and J. Pillow. Bayesian latent structure discovery
from multi-neuron recordings. In Advances in neural information processing
systems, 2016.
[39] S. W. Linderman and R. P. Adams. Discovering latent network structure in
point process data. In ICML, pages 1413–1421, 2014. arXiv:1402.0914.
[40] S. Mensi. A new mathematical framework to understand single neuron computations. 2014. PhD Thesis.
[41] S. Mensi, R. Naud, and W. Gerstner. From stochastic nonlinear integrate-andfire to generalized linear models. NIPS, 2011(1):1377–1385, 2011.
[42] G. Mohler. Marked point process hotspot maps for homicide and gun crime
prediction in chicago. International Journal of Forecasting, 30(3):491–497,
2014.
[43] G. Mohler, M Short, P Brantingham, F Schoenber, and G. Tita. Self-exciting
point process modeling of crime. Journal of the American Statistical Association, 106(493):100–108, 2012.
[44] S. Negahban and M. Wainwright. Estimation of (near) low-rank matrices with
noise and high-dimensional scaling. The Annals of Statistics, 39(2):1069–1097,
2011.
[45] S. Negahban and M. J. Wainwright. Restricted strong convexity and weighted
matrix com- pletion. 13:1665–1697, 2012.

37

[46] Sahand Negahban, Pradeep Ravikumar, Martin Wainwright, and Bin Yu. A
unified framework for high-dimensional analysis of m-estimators with decomposable regularizers. In Statistical Science, volume 27, pages 538–557,
2010.
[47] A. K. Oh, Z. T. Harmany, and R. M. Willett. To e or not to e in Poisson image
reconstruction. In Proceedings of the IEEE International Conference on Image
Processing (ICIP), 2014.
[48] S. A. Pasha and V. Solo. Hawkes-Laguerre reduced rank model for point
processes. In ICASSP, 2013.
[49] J. W. Pillow, J. Shlens, L. Paninski, A. Sher, A. M. Litke, E. J. Chichilnisky,
and E. P. Simoncelli. Spatio-temporal correlations and visual signalling in a
complete neuronal population. Nature, 454:995–999, 2008.
[50] M. Raginsky, R. Willett, C. Horn, J. Silva, and R. Marcia. Sequential anomaly
detection in the presence of noise and limited feedback. IEEE Transactions on
Information Theory, 58(8):5544–5562, 2012.
[51] A. Rakhlin and K. Sridharan. On equivalence of martingale tail bounds and
deterministic regret inequalities. 2015. https://arxiv.org/abs/1510.03925.
[52] G. Raskutti, M. J. Wainwright, and B. Yu. Minimax rates of estimation for
high-dimensional linear regression over `q -balls. IEEE Trans. on Information
Theory, 57(10):6976–6994, 2011.
[53] B. Recht, M. Fazel, and P. Parrilo. Guaranteed minimum-rank solutions
of linear matrix equations via nuclear norm minimization. SIAM Review,
52(3):471–501, 2010.
[54] P. Reynaud-Bouret and S. Schbath. Adaptive estimation for Hawkes processes;
application to genome analysis. Annals of Statistics, 38(5):2781–2822, 2010.
[55] J. Silva and R. Willett. Hypergraph-based anomaly detection in very large
networks. IEEE Transactions on Pattern Analysis and Machine Intelligence,
31(3):563–569, 2009.
[56] A. C. Smith and E. N. Brown. Estimating a state-space model from point
process observations. Neural Computation, 15:965–991, 2003.
[57] A. Stomakhin, M. B. Short, and A. Bertozzi. Reconstruction of missing data in
social networks based on temporal patterns of interactions. Inverse Problems,
27(11), 2011.
38

[58] B. Watson, D. Levenstein, J. Greene, J. Gelinas, and G. Buzaki. Multiunit spiking activity recorded from rat frontal cortex (brain regions mpfc,
ofc, acc, and m2) during wake-sleep episode wherein at least 7 minutes
of wake are followed by 20 minutes of sleep. 2016. Dataset:CRCNS.org.
http://dx.doi.org/10.6080/K02N506Q.
[59] B. Watson, D. Levenstein, J. Greene, J. Gelinas, and G. Buzaki. Network
homeostasis and state dynamics of neocortical sleep. Neuron, 90(4):839–852,
2016.
[60] A. Weber and J. Pillow. Capturing the dynamical repertoire of single neurons
with generalized linear models. 2016. https://arxiv.org/abs/0903.4742.
[61] S. Wright, R. Nowak, and M. Figueiredo. Sparse reconstruction by separable
approximation. IEEE Transactions on Signal Processing, 57, 2009.
[62] Ming Yuan and Yi Lin. Model selection and estimation in regression with
grouped variables. 68:49–67, 02 2006.
[63] J. Zhang, J. Chen, S. Zhi, Y. Chang, P. Yu, and J. Han. Link prediction across
aligned networks with sparse and low rank matrix estimation. ICDE, 33, 2017.
[64] Liu. Zhang and L. Vandenberghe. Interior-point method for nuclear norm
approximation with application to system identification. SIAM J. Matrix Anal.
Appl., 31(3):1235–1256, 2009.
[65] K. Zhou, H. Zha, and L. Song. Learning social infectivity in sparse lowrank networks using multi-dimensional Hawkes processes. In Proceedings
of the 16th International Conference on Artificial Intelligence and Statistics
(AISTATS), 2013.
[66] J. Zhuang, T. Mei, S. Hoi, X. Hua, and Y. Zhang. Community discovery from
social media by low-rank matrix recovery. TIST, 5(4), 2015.

A
A.1

Proofs
Proof of Theorem 3.1

Theorem 3.1 is the combination of results from [46] and [28]. We give a proof for
the sake of completeness but claim no originality of techniques. For the first part of

39

the proof, we follow Theorem 1 from [46]. By the definition of Ab and properties of
Bregman divergence for strongly convex functions, we have
Rmin X X >
1 XX
∗
b
t,m 4>
(4m. g(Xt ))2 ≤
m. g(Xt ) + λ(kA kR − kAkR ).
2T m t
T m t
where Rmin is a strong convexity parameter for ex on the domain x ∈ [Rmin , Rmax ].
Next note that
XX
m

t,m 4>
m. g(Xt ) =

XX

4m,m0

m m0

t

= h4,

X

X

Xt,m0 t,m

t
>

t g(Xt ) i

t

≤ k4kR k

X

t g(Xt )> kR∗ .

t

Thus, assuming λ/2 > T1 k

>
t t g(Xt ) kR∗

P

we have

1 XX
λ
∗
b
b .
t,m 4>
k4kR + λkA∗ kR − λkAk
R
m. g(Xt ) + λ(kA kR − kAkR ) ≤
T m t
2
(A.1)
Then
∗
∗
b
kAk
R = kA + 4kR = kA + 4M⊥ + 4M kR .
⊥

Since R is decomposable with respect to the subspaces (M, M ), we have
∗
b
kAk
R ≥ kA kR + k4M ⊥ kR − k4M kR .

Thus
λ
λ
3λ
3λ
b
k4kR + λkA∗ kR − λkAk
k4M kR − k4M⊥ kR ≤
k4M kR . (A.2)
R ≤
2
2
2
2
Recalling that Ψ(M) is the subspace compatibility constant, we have
k4M kR ≤ Ψ(M)k4M kF ≤ Ψ(M)k4kF .
It follows that
λ
3λ
b
k4kR + λkA∗ kR − λkAk
Ψ(M)k4kF .
R ≤
2
2
Let k4k2T =

1
T

P

m

>
2
t (4m. g(Xt ))

P

k4k2T ≤

and therefore
3λ
Ψ(M)k4kF .
Rmin
40

From here, we reduce the lower bound into the restricted eigenvalue condition.
Denote the subsets
BR = {B ∈ RM ×M K : kBM⊥ kR ≤ 3kBM kR }
and
0
= {B ∈ BR : kBkF = 1}.
BR

Note that Equation (A.2) implies that 4 ∈ BR . Let T = {p, 2p, . . . , T } so
|T |/|T | = p1 (here we assume Tp is an integer for simplicity). By Assumption
1 we have
1 XX >
4m. E[g(Xt )g(Xt )> |Xt−p ]4m.
T t∈T m
1 XX >
2
− sup
(bm. g(Xt ))2 − E[(b>
m. g(Xt )) |Xt−p ]
B∈BR T t∈T m
ω
1 XX >
2
(bm. g(Xt ))2 − E[(b>
≥ k4k2F − sup
m. g(Xt )) |Xt−p ].
p
B∈BR T t∈T m

k4k2T ≥

We want to show
sup
B∈BR

ωk4k2F
1 XX >
2
(bm. g(Xt ))2 − E[(b>
g(X
))
|X
]
≤
t
t−p
m.
T t∈T m
2p

(A.3)

0
.
with high probability, and we note that it suffices to show this for all B ∈ BR
M ×M
Now we define the matrix G ∈ R
as follows:

G :=

1 X
(g(Xt )g(Xt )T − E[g(Xt )g(Xt )T |Xt−p ]).
T t∈T

Following the definition of G, and denoting each entry Gm,m0 ,
sup
0
B∈BR

M
h
i
X
1 XX >
2
(bm. g(Xt ))2 − E (b>
b>
m. g(Xt )) |Xt−p = sup
m. Gbm.
0
T t m
B∈BR
m=1

≤ sup kBk22,1 max0 |Gm,m0 |.
0
B∈BR

m,m

2
Recall that supB∈BR
0 kBk2,1 ≤ µR by Assumption 3. Hence

sup
0
B∈BR

i
h
1 XX >
2
(bm. g(Xt ))2 − E (b>
m. g(Xt )) |Xt−p ≤ µR max0 |Gm,m0 |.
m,m
T t m

41

Note that each entry Gm,m0 is a martingale and |Gm,m0 | ≤ 2U 2 . Therefore we
can apply the Azuma-Hoeffding inequality [5]. For completeness, we state the
Azuma-Hoeffding inequality as Theorem A.5 in Section A.9. If we let
Yn :=

n
1X
(g(Xt )g(Xt )T − E[g(Xt )g(Xt )T |Xt−p ]),
T t=0

where n = 0, 1, 2, .., |T |, and we set t =
Appendix A.9, we have
ω
P |Gm,m0 | ≥
2µR p

!

ω
2µR p

and cn =

2U 2
T

as in Theorem A.5 in
!

T ω2
≤ 2 exp −
.
32U 4 p2 µ2R

Applying a union bound,
ω
P max0 |Gm,m0 | ≥
m,m
2µR p

!

!

T ω2
≤ 2M exp −
.
32U 4 p2 µ2R
2

Hence if we set

128U 4 p2 µ2R log M
,
ω2
(A.3) holds with probability at least
2
1 − 2,
M
ω
2
2
guaranteeing that k4kT ≥ 2p k4kF with this same probability. Putting everything
together, we have
36pΨ(M)2 λ2
k4k2F ≤
2
Rmin
ω2
with probability at least 1 − M22 for
T >

128p2 U 4 µ2R log M
T ≥
.
ω2

A.2

Proof of Lemma 4.1

Before we prove Lemma 4.1, we first need the following supporting lemma.
Lemma A.1. Let Z = min(bλc , Ũ ). Define the random variables X ∼ Poisson(λ),
X = min(X, Ũ ) and

0 if X ≤ Z
Y =
.
1 if X > Z
Then
Var(Y ) ≤ Var(X).
42

Proof We write
X = (Y + Z) + (X − Y − Z).
Since Z is a constant we have
Var(X) = Var(Y ) + Var(X − Y ) + 2Cov(Y, X − Y )
and it suffices to show Cov(Y, X − Y ) ≥ 0. Conditioning on Y gives
h

i

E (Y − E[Y ])(X − Y − E[X − Y ]) =
h

p(Y = 1)(1 − E[Y ])EX,Y |Y =1 X − Y − E[X − Y ]
h

i

(A.4)

i

+p(Y = 0)(−E[Y ])EX,Y |Y =0 X − Y − E[X − Y ] .
Now observe that Y = 1 implies X − Y ≥ Z (where we rely on the fact that X and
Z are both integers) while Y = 0 implies X − Y ≤ Z; then
E[X − Y |Y = 1] ≥ E[X − Y ]

(A.5)

E[X − Y |Y = 0] ≤ E[X − Y ].

(A.6)

and

We argue that both terms in the sum in (A.4) are non-negative. For the first term, we
have
h

i

h

i

EX,Y |Y =1 X − Y − E[X − Y ] = EX,Y |Y =1 [X − Y ] − E[X − Y ] ≥ 0
by (A.5).
EX,Y |Y =0 X − Y − E[X − Y ] = EX,Y |Y =0 [X − Y ] − E[X − Y ] ≤ 0
by (A.6). Finally note that E[Y ] ∈ (0, 1) so that both terms in (A.4) are indeed
non-negative, and therefore Cov(Y, X − Y ) ≥ 0 as claimed.
We now prove Lemma 4.1. Note that
Cov(min(Xt,m , Ũ ), min(Xt,m0 , Ũ )|Xt−1 )) = 0
for m 6= m0 . We have
E[min(Xt , Ũ ) min(Xt , Ũ )> |Xt−1 ] =
E[min(Xt , Ũ )|Xt−1 ]E[min(Xt , Ũ )|Xt−1 ]> + Diag(Var(min(Xt , Ũ )|Xt−1 ))
43

where the first matrix is positive semi-definite because it is the outer product of a
vector with itself. Thus, to come up with a lower bound for our original matrix,
we just need to lower bound the smallest element of Var(min(Xt , Ũ )|Xt−1 ). This
amounts to lower bounding the variance of min(Xλ , Ũ ) where Xλ is a Poisson
random variable with mean λ ∈ [Rmin , Rmax ]. Define
Yλ =


0

if min(Xλ , Ũ ) ≤ min(bλc , Ũ )
1 if min(Xλ , Ũ ) > min(bλc , Ũ ).

By Lemma A.1, Var(min(Xλ , Ũ )) ≥ Var(Yλ ) so our problem reduces to lower
bounding the variance of Yλ . We argue that
Var(Yλ ) ≥ min(Var(YRmin ), Var(YRmax ))
by considering two cases. When analyzing these cases, we use the fact that Var(Yλ )
will be minimized when the probability of outcome (0) is either maximized or
minimized. We take Rmin ≤ 51 to make the exposition clearer. At the end of the
proof we discuss the Rmin > 15 scenario which is virtually identical.
Case 1: λ ∈ [Rmin , Ũ ) where Ũ may be either greater than or less than Rmax
In this scenario

0 if X ≤ bλc
λ
Yλ =
1 if Xλ > bλc .
We claim Var(YRmin ) ≤ Var(Yλ ) for λ ∈ [Rmin , Ũ ). To do this, we look at two
subcases.
First, if 1 ≤ λ ≤ Ũ , then basic properties of the median of the Poisson distribution imply that the probability of outcome (0) will be between 15 and 45 and so
4
Var(Yλ ) ≥ 25
. For the second case where Rmin ≤ λ < 1, outcome (0) corresponds
to
P (X = 0|X ∼ Poisson(λ)) = exp(−λ) ≤ exp(−Rmin ).
Since Rmin ≤ 15 we get that exp(−Rmin ) > 45 . Combining the two cases, we have
concluded that Var(Yλ ) is minimized on λ ∈ [Rmin , Ũ ] at λ = Rmin . Now, when
λ = Rmin ,
Var(Yλ ) = exp(−Rmin )(1 − exp(−Rmin )).
Since
f (x) =

e−x (1 − e−x )
x

44

decreases monotonically on the interval (0, 51 ], we have
min1 f (x) ≥ f

1
1
≥ .
5
2

 

x∈(0, 5 ]

Using the fact that Rmin ∈ (0, 15 ] we conclude
1
Var(YRmin ) = exp(−Rmin )(1 − exp(−Rmin )) ≥ Rmin .
2

(A.7)

Case 2: λ ∈ [Ũ , Rmax ]
Next we consider the variance when λ ≥ Ũ . By the same reasoning as in Case
1, outcome (0) can have probability no larger than 45 . It remains to consider when
outcome (1) can have probability larger than 45 . It is clear that for λ ∈ [Ũ , Rmax ],
outcome (1) will be maximized for λ = Rmax . When λ = Rmax , we directly
compute the variance as
Var(YRmax ) =

Ũ
−1
X
i=0

i
Rmax

−Rmax

exp
i!


1 −

Ũ
−1
X
i=0

i
Rmax

−Rmax

exp
i!




= κ.

(A.8)

Combining Case 1 and Case 2, we get that
Var(Yλ ) ≥ min (Var(Rmin ), Var(Rmax ))
and combining Equations (A.7) and (A.8) gives the final result.
If Rmin ≥ 51 an identical argument shows that for Case 1 where λ ∈ [Rmin , Ũ ),
4
Var(Yλ ) ≥ 25
and for Case 2 where λ ∈ [Ũ , Rmax ], Var(Yλ ) ≥ κ. Hence, a lower
4
bound on Var(Yλ ) covering all possible values of Rmin would be min( 12 Rmin , κ, 25
).
1
In main body of the paper we present the bound for the Rmin ≤ 5 scenario in order
to make the statement more interpretable.

A.3

Proof of Theorem 4.2

Recall that the AR(2) model is a special case of (2.3) with φ1 [t] = I{t=1} and
φ2 [t] = I{t=2} . With these choices of basis functions,
g(Xt ) = [min(X1 , Ũ ), min(X2 , Ũ )]> .
A computation shows that if we choose to condition on Xt−1 as in the proof of
Lemma 4.1 we get a singular matrix. However, Assumption 1 allows us to condition
45

on Xt−p for any p > 0 and so for this example it will be easiest to condition on Xt−2 .
We have
E[g(Xt )g(Xt )> |Xt−2 ] = E[g(Xt )|Xt−2 ]E[g(Xt )|Xt−2 ]> + Cov(g(Xt )|Xt−2 ).
The first matrix is an outer product of a vector with itself so it is positive semidefinite and it suffices to lower bound the eigenvalues of the covariance matrix
Cov(g(Xt )|Xt−2 ). Recall that a matrix B is said to be strictly diagonally dominant if
bi,i −

X

|bi,j | ≥ ω > 0

j6=i

for all i, and the eigenvalues of a symmetric strictly diagonally dominant matrix are
lower bounded by ω. To lower bound the eigenvalues of the covariance matrix, we
will show it is strictly diagonally dominant. We break the rows up into two cases.
The first case corresponds to rows where the diagonal depends on a lagged count
Xt−1,m while the second case corresponds to rows where diagonal depends on a
count Xt,m without a lag.
Case 1: Rows 1 through M The first M rows of Cov(g(Xt )|Xt−2 ) have their diagonal of the form Var(Xt−1,m |Xt−2 ) ≥ Rmin . We have Cov(Xt−1,m , Xt−1,m0 |Xt−2 ) =
0 for all m0 6= m. If node m is not a parent of m0 , then Xt−1,m and Xt,m0 are independent conditioned on Xt−2 , so Cov(Xt−1,m , Xt,m0 |Xt−2 ) = 0. All that remains is
to control Cov(Xt−1,m , Xt,m0 |Xt−2 ) for the ρ(c)
m children of m.
To do this, recall the decomposition Xt,m0 = exp(νm0 + a>
m.0 g(Xt−1 )) + t,m0 .
For the remainder of the proof we let fm0 (Xt ) = exp(νm0 + a>
g(X
t )) for notational
m.0
simplicity and note that the Poisson noise term t,m is zero mean conditioned on
Xt−1,m . Hence
Cov(Xt−1,m , Xt,m0 |Xt−2 ) = Cov(Xt−1,m , fm0 (Xt−1 )|Xt−2 ).
Since fm0 (Xt−1 ) takes values in the interval [Rmin , Rmax ], the variance of fm0 (Xt−1 )
is bounded by a scaled Bernoulli random variable which takes values 0 with proba2
min )
bility 12 and Rmax − Rmin with probability 12 . This variance is equal to (Rmax −R
4
and therefore
Cov(Xt−1,m , fm0 (Xt−1 )|Xt−2 ) ≤
q

√

Var(Xt−1,m |Xt−2 )Var(fm0 (Xt−1 )|Xt−2 ) ≤

Hence the off diagonal entries sum to at most
√
Rmax (Rmax − Rmin )
ρ(c)
m
2
46

Rmax (Rmax − Rmin )
.
2

so for these rows of the covariance matrix we have
Cov(g(Xt ))i,i −

X

Cov(g(Xt ))i,j

j6=i

√
ρ(c)
Rmax (Rmax − Rmin )
m
≥ Rmin −
.
2

(A.9)

Case 2: Rows M + 1 through 2M We next consider the final M rows of the
covariance matrix whose diagonal is of the form Var(Xt,m |Xt−2 ) ≥ Rmin . We know
Cov(Xt,m , Xt−1,m0 |Xt−2 ) will be zero whenever node m0 is not a parent of m, and
for the ρ(p)
m parents of m, the covariance is bounded below by
√
Rmax (Rmax − Rmin )
2
just as in the previous paragraph.
Finally, we need to consider Cov(Xt,m , Xt,m0 |Xt−2 ). When m and m0 are not
siblings this covariance will be zero. When they do share a parent, we again recall
the decomposition Xt,m = fm (Xt−1 ) + t,m and Xt,m0 = fm0 (Xt−1 ) + t,m0 and note
that the t,m and t,m0 are zero mean conditioned on Xt,m0 and Xt,m respectively.
Therefore
Cov(Xt,m , Xt,m0 |Xt−2 ) = Cov(fm (Xt−1 ), fm0 (Xt−1 )|Xt−2 )
and using the fact that each fi (Xt ) takes values in the interval [Rmin , Rmax ] it follows
that this covariance is bounded by
q

Var(fm (Xt−1 )|Xt−2 )Var(fm0 (Xt−1 )|Xt−2 ) ≤

(Rmax − Rmin )2
.
4

Recall that ρ(s)
m denotes the number of siblings of m. Overall we have concluded that
the sum of the off diagonal entries for the first m rows is at most
√
ρ(p)
(Rmax − Rmin )2
Rmax (Rmax − Rmin ) ρ(s)
m
+ m
2
4
so that for these rows we have
Cov(g(Xt ))i,i −

X

Cov(g(Xt ))i,j >

j6=i

√
ρ(p)
Rmax (Rmax − Rmin ) ρ(s)
(Rmax − Rmin )2
m
Rmin −
− m
. (A.10)
2
4
We conclude that the smallest eigenvalue of the covariance matrix is lower bounded
by the minimum of the two lower bounds in Equations (A.9) and (A.10), and we
define this minimum to be rρ .
47

A.4

Proof of Lemma 4.3

We know from Equation (A.2) that k4S ⊥ k1 ≤ 3k4S k1 and since
k4k1 = k4S ⊥ k1 + k4S k|1
√
it follows that k4k1 ≤ 4k4S k1 . Recall that kvk1 ≤ skvk2 for any s-sparse vector
v. Thus we have
√
√
k4k1 ≤ 4k4S k1 ≤ 4 sk4S kF ≤ 4 sk4kF .
For Assumption 4 we use a concentration result due to [33] in a similar manner as in
P
[28]. The result is restated as Theorem A.4 below. Define Yn = T1 nt=1 g(Xt )i t,j
and note the following values
Yn − Yn−1 =
and
n
X



g(Xn )i
n,j
T

g(Xt )i
Mnk =
E
t,j
T
t=1



!k

Xt−1  .

Thus Yn is a martingale. We have g(Xt )i ≤ U , and by Lemma 1 in [28],
t,j ≤ Xt,j ≤ C log(M T )
for all t, j with probability at least 1 − exp(−cM T ). Thus,
|Yn − Yn−1 | ≤

CU log(M T )
=: B
T

with this same probability. Next, note that
Mn2

n
X

"

g(Xt )2i 2
=
E
 Xt−1
T 2 t,j
t=1
n
1 X
≤ 2
U 2 Rmax
T t=1
n
d 2.
= 2 U 2 Rmax =: M
n
T

#

Here we use that E[2t,j |Xt−1 ] is the variance of a Poisson random variable with mean
bounded by Rmax , so it must also be bounded by Rmax . Next, we bound Mnk :
n
X



g(Xi )m
Mnk :=
E
(Xi,l − E[Xi,l |Xi−1 ])
T2
i=1
48

!k



Xi−1  ≤ B k−2 Mn2 .

In the language of Theorem A.4,
Dn :=

X γk

k!

k

Mnk

2

d X γkBk
M
n
c .
≤ 2
=: D
n
B k k!
d.
Let D̃n corresponds to the negative sequence of Dn , and so it is still bounded by D
n
Using Markov’s inequality, we get

P(|Yn | ≥ y) = P(Yn ≥ y) + P(−Yn ≤ y)
≤ E[exp(γYn )] exp(−γy) + E[exp(−γYn )] exp(−γy)
d − γy) + E[exp(−γY ) − D̃ ] exp(D
d − γy).
≤ E[exp(γYn − Dn )] exp(D
n
n
n
n

Using Theorem A.4 we conclude that
E[exp(γYn − Dn )] ≤ 1
and E[exp(γYn − D̃n )] ≤ 1 so
c − γy).
P(|Yn | ≥ y) ≤ 2 exp(D
n

We set
1
yB
γ = log 1 +
d2
B
M
n

!

and to simplify things note that (1 + x) log(1 + x) − x ≥
together gives


P(|YT ≥ y|) ≤ 2 exp 

Putting everything



−3y 2
2yB

3x2
.
2(x+3)

d2
+ 6M



n

!

3y 2 T
.
= 2 exp
2U C log(M T )y + 6Rmax
Now, recall that

λ
4CRmax U 2 log2 (M T )
√
=
2
T

49

and setting y =

λ
2

gives

48C 2 U 4 log4 (M T )
λ
√
P(YT ≥ ) ≤ 2 exp
2
2C 2 U 3 log3 (M T )/ T + 6U 2 Rmax


!



48U log(M T )

= 2 exp  √
max
2/ T + C 2 U6R
log3 (M T )
!

48U log(M T )
≤ 2 exp
.
8
Taking a union bound over all i, j gives us


P max
i,j

T
√ 
1 X
g(Xt )i t,j ≥ 4U 2 log2 (M T )/ T
T t=1

≤ exp(log(2M 2 ) − 6U log(M T ))
≤ exp(3 log(M T ) − 6U log(M T ))
= exp(−c log(M T ))
for c = 6U − 3 which is positive since U ≥ 1. In the final statement of the proof
we assume C log(M T ) ≥ U and replace U with C log(M T ) in order to limit the
number of constants and make the crucial dependencies clear. This assumption
should hold for reasonable choices of U in the settings we imagine in practice, but if
not, a factor of log2 (M T ) can be replaced by U 2 in the final bound.

A.5

Proof of Lemma 4.5

For Assumption 2, note that any A ∈ SG satisfies A.c = 0 for i 6∈ SG . Therefore
k4kG ≤ 4k4S G kG = 4

X

k4.i k2

i∈SG

√
≤ 4 sG k4S G kF
√
≤ 4 sG k4kF .
For Assumption 3, we have k4k2G ≤ 16sG k4k2F from the previous paragraph, and
we claim kAk2,1 ≤ kAkG for any matrix A. To see this, we compute
kAk2G

=

2

s
X X

a2 
r,c

c

=

XX
c

r

s

X

(

r

c0

X

a2r,c )(

a2r,c0 )

r

while
kAk22,1 =

X X

(

r

|ar,c |)2 =

c

XXX
c

50

c0

r

|ar,c kar,c0 |.

To complete the proof, we fix c, c0 and need to show
X

s

|ar,c kar,c0 | ≤

X

X

a2r,c )(

(

r

r

a2r,c0 ),

r

or equivalently that
X

(

|ar,c kar,c0 |)2 ≤ (

X

r

a2r,c )(

X

r

a2r,c0 ).

r

We have
X

(

|ar,c kar,c0 |)2 =

r

XX
r

|ar,c ar,c0 ar0 ,c ar0 ,c0 |.

r0

Let J denote all two element combinations of M and we can write
XX
r

|ar,c ar,c0 ar0 ,c ar0 ,c0 | =

X

(ar,c ar,c0 )2 +

r

r0

X

2|ai,c aj,c ai,c0 aj,c0 |.

(i,j)∈J

On the other hand,
X

(

r

a2r,c )(

X

a2r,c0 ) =

r

XX
r

=

X

a2r,c a2r0 ,c0

r0

(ar,c ar,c0 )2 +

r

X

(ai,c aj,c0 )2 + (aj,c ai,c0 )2 .

(i,j)∈J

The proof follows from noting that
(ai,c aj,c0 )2 + (aj,c ai,c0 )2 ≥ 2|ai,c aj,c ai,c0 aj,c0 |
for any real numbers ai,c , ai,c0 , aj,c , aj,c0 .
For Assumption 4, we rely on Theorem 1 in [51] which is restated as Theorem
P
A.3. In our setup, we need to bound the l2 norm of the mth column of T1 t t g(Xt )> .
Note that the l2 norm is 2-smooth, because for any x, y ∈ Rm we have
kx + yk2 + kx − yk2 = hx + y, x + yi + hx − y, x − yi
= 2hx, xi + 2hy, yi.
In the language of Theorem A.3, for a fixed m we form a martingale difference
sequence {Zt } with

1
Zt =
t g(Xt )>
.m
T
so that
s
s
1 X
g(Xt )m X 2
kZt k2 =
(g(Xt )m t,m0 )2 =
t,m0 .
T m0
T
m0
51

We know g(Xt )m ≤ U and by Lemma 1 from [28] t,m0 ≤ C log(M T ) with probability at least 1 − exp(−cM T ). We conclude
kZt k2 ≤
and thus

T
X

1 q
U M C log2 (M T )
T

kZt k22 ≤ CU 2 log2 (M T )

t=1

M
.
T

To compute the constant Qmax appearing in Theorem A.3 we let R(x) = x> x so
that ∇R(x) = x. Then for any x, y in the unit ball with respect to the k · k2 norm,
we have
DR (x, y) = kxk22 − kyk22 − hy, x − yi ≤ kxk22 + kyk2 kx − yk2 ≤ 3.
√
by Cauchy-Schwarz. Thus we can take Qmax = 3. To simplify, we note Wn ≤ Vn
and
q
(E[ Vn + Wn ])2 ≤ E[Vn + Wn ] ≤ 2Vn .
√
√
Further, 2.5Qmax ( Vn + 1) ≤ 5 Vn . With these simplifications, Theorem 1 from
[51] says that

1 X
P
t g(Xt )>
.m
T
t

!

> (5 + 2u)Vn ≤
2

√

!

u2
.
2 exp −
16

Setting u = log(T ) and plugging in our values for Vn we conclude that



1 X
P
t g(Xt )>
.m
T
t

s

2



M √
> CU log2 (M T )
≤ 2exp(− log2 (M T )).
T

Taking a union bound over all m, we get that
s


1 X
M
t g(Xt )>
≤ CU log2 (M T )
.m
T
T
t
2
√
for all m with probability at least 1 − 2 exp(− log(M T )) = 1 −

A.6

√
2
.
MT

Proof of Lemma 4.7

√
For Assumption 3, from the statement of Lemma 4.7 we have kA∗ k22,1 ≤ D M and
√
we search for Ab over the ball {A : kAk22,1 ≤ D M }. Thus
√
sup kBk22,1 ≤ 2D M = µR .
0
B∈BR

52

In contrast to the sparsity case, Assumption 2 is nontrivial to verify in the low-rank
case because W =
6 W. However, this condition was shown in Lemma 3.4 of [53].
For Assumption 4 we rely on the notion of a k-regular normed vector space
defined in Section A.9 as well as Theorem 2.1 from [34] which is stated in Theorem
A.2. Further, Example 3.1 in [34] establishes that (RM ×M K , k · k∗ ) is k-regular
for k = 3 log(min(M, N )). In the language of Theorem A.2 we form a martingale
difference sequence {ζt } with ζt = T1 t g(Xt )> and then
kζt kop =

M
>
1 X
t g(Xt )
=
t,m g(Xt )m .
T
T m=1

Consider the random variable
M
X

t,m g(Xt )m .

(A.11)

m=1

We have g(Xt )m ≤ U and t,m ≤ C log(M T ) for all t, m with probability at least
e−cM T by Lemma 1 from [28]. Further, conditioned on Xt , the t,m g(Xt )m are all
independent, so (A.11) is a sum of zero mean independent random variables bounded
by CU log(M T ). Hence, we have
M
X

E[t,m g(Xt )m |Xm−1 ] ≤ CM log2 (M T )U 2 ,

m=1

and applying Bernstein’s inequality gives


P |>
t g(Xt )| >

√



log4 (M T )/2



M log2 (M T ) ≤ 2 exp −

C log2 (M T )U 2 +

CU log(M
T)
√
3 M


.

Therefore


P |>
t g(Xt )| >

√

M log2 (M T ) for at least one t





≤ 2 exp log(T ) − log4 (M T ))/(2C log2 (M T )U 2 + 1)
log4 (M T )
≤ exp −
.
4C log2 (M T )U 2 + 2
!

We apply Theorem A.2 with k = 3 log(M ),
T
X
i=1

σi2

M log4 (M T )
=
T
53



and γ = log(T ). This gives
s


√
√
M
1 X
k
t g(Xt )> kop > (3 2 log(M )+ 2 log(T )) log4 (M T )
P
T t
T


− log2 (T )
≤ exp
.
2
!

A.7

Proof of Proposition 1

(c)
In this proof, we take λ(c)
m (τ ) to mean λm (τ ; Xτ ). Using the approximation

Z 4t
(t−1)4

(c)
(c)
λm
(τ )dτ ≈ 4λm
(4t)

we derive an approximate sampled Hawkes (SH) log-likelihood proportional to
`H (XT 4 |{λ(c)
m }m )

≈

M X
T
X

(c)
(c)
[Xt,m log λ(c)
m (4t)−4λm (4t)] =: `SH (XT |{λm }m ).

m=1 t=1

If Xt,m were generated according to (2.3) with intensity (6.5) for T = 1, . . . T , then,
ignoring terms independent of A∗ ,
`P (Xt |{4λ(c)
m (X4t )}t,m ) :=

M X
T
X

(c)
[Xt,m log 4λ(c)
m (4t) − 4λm (4t)].

m=1 t=1

Note that
(c)
`P (XT |{4λ(c)
m (X4t )}t,m = `SH (XT |{λm }m ) + C

where the constant C depends on 4 but is independent of λ(c)
m .

A.8

Extension to more general saturation functions

In the main body of the paper the only saturation function we consider is f (x) =
min(x, Ũ ) for purposes of simplicity, but our theory extends to a larger class of
saturation functions f . However, for our analysis it is crucial to assume that the
function f is bounded so that we can define the maximum and minimum rates, Rmin
and Rmax , from which each observation is drawn. The only place where we rely
on the structure of f beyond its boundedness is in proving the restricted eigenvalue
condition in Assumption 1. In the case of the ARMA(1, 1) model, we show our
results extend to monotonic differentiable functions in Proposition 2 below.
54

Proposition 2. Suppose (Xt )Tt=1 is generated according to the ARMA(1, 1) model
in (4.1) with a general saturation function f applied entrywise to the vector Xt .
Suppose f is bounded on R, monotonically increasing, and differentiable with
f 0 (x) ≥ c on [0, Rmax ]. Then
1
λmin [E[g(Xt )g(Xt )> |Xt−1 ] ≥ c2 min( Rmin , κ).
2
Proof We have
E[g(Xt )g(Xt )> |Xt−1 ] = E[g(Xt )|Xt−1 ]E[g(Xt )|Xt−1 ]> + Diag(Var(g(Xt )|Xt−1 ))
where the first matrix is positive semi-definite because it is the outer product of a
vector with itself. Thus, to come up with a lower bound for our original matrix, we
just need to lower bound the smallest element of Var(g(Xt )|Xt−1 ). This amounts to
lower bounding the variance of f (X) where X is a Poisson random variable with
mean λ ∈ [Rmin , Rmax ].
Let p = P(X ≤ bλc) so 1 − p = P(X ≥ dλe). Consider the random variable X 0
which takes the value bλc with probability p and dλe with probability 1 − p. Since f
is monotonic, the argument from Lemma A.1 shows that Var(f (X 0 )) ≤ Var(f (X))
so we reduce our problem to lower bounding the variance of f (X 0 ).
Note that this variance is equal to the shifted random variable X 00 defined by
f (X 0 ) = 0 with probability p and f (dλe) − f (bλc) with probability 1 − p which is
a scaled Bernoulli random variable with variance
(f (dλe) − f (bλc))2 p(1 − p).
Since f 0 (x) ≥ c on [0, Rmax ],
f (dλe) − f (bλc) ≥ c
and so the lower bound on our variance becomes c2 p(1 − p). Finally, by Lemma 4.1
we have p(1 − p) ≥ min( 21 Rmin , κ) which completes the proof.

A.9

Supplemental Theorems

Definitions Before introducing martingale concentration results, we give the following definitions.
Definition 1. A Banach space (E, k · k) is s-smooth if there exists C > 0 satisfying
kx + yks + kx − yks ≤ 2kxks + 2C s kyks
for all x, y ∈ E.
55

Note that (RM , k · k2 ) is 2-smooth with C = 1 because
kx + yk2 + kx − yk2 = hx + y, x + yi + hx − y, x − yi = 2hx, xi + 2hy, yi.
Definition 2. A Banach space (E, k · k) is k − regular if there exists k+ ∈ [1, k]
along with a norm k · k+ such that (E, k · k+ ) is k+ -smooth and
k
kxk2
k+

kxk2 ≤ kxk2+ ≤
for all x ∈ E.

By Example 3.3 from [34], the space (RM ×N , k · k∗ ) is k-regular for k =
3 log(min(M, N )).
Theorem A.2. (Theorem 2.1.iii in [34])
Let (E, k · k) be k-regular and let ζi be an E-valued martingale difference
P
sequence with kζi k ≤ σi . Let SN = N
i=1 ζi . Then
v





N
√
√ u
uX
γ2

P kSN k ≥ ( 2k + 2γ)t σi2  ≤ exp(− ).
2
i=1

Theorem A.3. (Theorem 1 in [51]) Let (E, k · k) be a 2-smooth Banach space. Let
R be a function which is 1-strongly convex on the unit ball in the dual norm of k · k.
DR : B∗ × B∗ → R
be the Bregman divergence with respect to R, and finally let Q2max =
supx,y∈B∗ DR (f, g). Let Z1 , . . . , Zn be a martingale difference sequence with
P
P
Vn = nt=1 kZt k2 and Wn = nt=1 Et−1 kZt k2 . Then


Pk

n
X



r

q

q

Zt k > 2.5Qmax ( Vn + 1)+u Vn + Wn + (E[ Vn + Wn ])2 

t=1

≤

√

2 exp(−

u2
).
16

Theorem A.4. (Lemma 3.3 in [33]) Let (Yn ) be a martingale and let
Mnk =

n
X

E[(Yi − Yi−1 )k |Yi−1 ].

i=1

Let γ be such that for all i ≤ n, we have
E[exp(|γ(Yi − Yi−1 )|)] ≤ ∞.
56

Then
n = exp(γYn −

γk K
Mn )
k≥2 k!
X

is a super-martingale. Moreover, if Y0 = 0 then E[n ] ≤ 1.
Theorem A.5. (Azuma-Hoeffding inequality) Let (Yn ) be a martingale and |Yn −
Yn−1 | < cn . Then
t2

!

P(|YN − Y0 | ≥ t) ≤ 2 exp − PN 2 .
2 n=1 cn

57

