1

Recursive Sparse Point Process Regression with
Application to Spectrotemporal Receptive Field
Plasticity Analysis

arXiv:1507.04727v1 [] 16 Jul 2015

Alireza Sheikhattar, Jonathan B. Fritz, Shihab A. Shamma, and Behtash Babadi

Abstract—We consider the problem of estimating the sparse
time-varying parameter vectors of a point process model in an
online fashion, where the observations and inputs respectively
consist of binary and continuous time series. We construct a novel
objective function by incorporating a forgetting factor mechanism
into the point process log-likelihood to enforce adaptivity and
employ ℓ1 -regularization to capture the sparsity. We provide a
rigorous analysis of the maximizers of the objective function,
which extends the guarantees of compressed sensing to our
setting. We construct two recursive filters for online estimation
of the parameter vectors based on proximal optimization techniques, as well as a novel filter for recursive computation of
statistical confidence regions. Simulation studies reveal that our
algorithms outperform several existing point process filters in
terms of trackability, goodness-of-fit and mean square error. We
finally apply our filtering algorithms to experimentally recorded
spiking data from the ferret primary auditory cortex during
attentive behavior in a click rate discrimination task. Our analysis
provides new insights into the time-course of the spectrotemporal
receptive field plasticity of the auditory neurons.
Index Terms—Adaptive filtering; point process models; compressed sensing; neural signal processing; receptive field plasticity.

I. I NTRODUCTION
Analyses of spiking activity recorded from sensory neurons
have revealed three main features: first, neuronal activity is
stochastic and exhibits significant variability across trials;
second, the spiking statistics often undergo rapid changes
referred to as neuronal plasticity, in order to adapt to changing
stimulus salience and behavioral context; and third, the tuning characteristics of sensory neurons to the stimuli exhibit
a degree of sparsity. Examples include place cells in the
hippocampus [1] and spectrotemporally tuned cells in the
primary auditory cortex [2]. Hence, in order to gain insight
into the functional mechanism of the underlying neural system,
it is crucial to have a mathematical theory to simultaneously
capture the stochasticity, dynamicity and sparsity of neuronal
activity.
On one hand, the theory of point processes [3] has been
recently adopted as a mathematical framework to model the
stochasticity of neuronal data. Traditionally, these models have
A. Sheikhattar is with the Department of Electrical and Computer Engineering (ECE), University of Maryland, College Park, MD; J. B. Fritz is with the
Institute for Systems Research (ISR), University of Maryland, College Park,
MD; S. A. Shamma and B. Babadi are with ECE and ISR, University of
Maryland, College Park, MD (e-mails: arsha89@umd.edu; ripple@umd.edu;
sas@umd.edu; behtash@umd.edu).
Corresponding author: B. Babadi (e-mail: behtash@umd.edu).

been used to predict the likelihood of self-exciting processes
such as earthquake occurrences [4], [5], but have recently
found significant applications in the analysis of neuronal data
[6]–[12].
On the other hand, classic results in signal processing such
as the Least Mean Squares (LMS) and Recursive Least Squares
(RLS) algorithms [13] have created a framework to efficiently
capture the dynamics of the parameters in linear observation
models. Existing solutions in computational neuroscience have
adopted this framework to estimate the dynamics of neuronal
activity. For instance, in [7] an LMS-type point process filter
was introduced to study plasticity in hippocampal neurons.
In [14], more general adaptive filtering solutions based on
approximations to the Chapman-Kolmogorov equation were
introduced. Although quite powerful in analyzing neuronal
data, these solutions do not account for the sparsity of the
underlying parameters.
Finally, the theory of compressed sensing (CS) has provided
a novel methodology for measuring and estimating statistical
models governed by sparse underlying parameters [15]–[17].
In particular, for static linear and generalized linear models
(GLM) with random covariates and sparsity of the parameters,
the CS theory characterizes sharp trade-offs between the number of measurement, sparsity, and estimation accuracy [16],
[18]. The sparse solutions of CS are typically achieved using
batch-mode convex programs and greedy techniques. In online
settings, sparse adaptive filters have only been introduced in
the context of linear systems governed by sparse parameters
such as communication channels [19]–[21].
Despite significant progress in all these research fronts, a
unified framework to simultaneously capture the stochasticity,
dynamicity and sparsity of neuronal data is lacking. In this
paper, we close this gap by integrating techniques from point
process theory, adaptive filtering, and compressed sensing. To
this end, we consider the problem of estimating time-varying
stimulus modulation coefficients (e.g., receptive fields) from
a sequence of binary observations in an online fashion. We
model the spiking activity by a conditional Bernoulli point
process, where the conditional intensity is a logistic function of
the stimulus and its time lags. We will then design a novel objective function by incorporating the forgetting factor mechanism of RLS-type algorithms into the ℓ1 -regularized maximum
likelihood estimation of the point process parameters. We will
present theoretical guarantees that extend those of CS theory
and characterize fundamental trade-offs between the number
of measurements, forgetting factor, model compressibility, and

2

estimation error of the underlying point processes in the nonasymptotic regime. We will next develop two adaptive filters
for recursive estimation of the objective function based on
proximal gradient techniques, as well as a filter for recursive
computation of statistical confidence regions.
In order to validate our algorithms, we provide simulation
studies which reveal that the proposed adaptive filtering algorithms significantly outperform existing point process filters
in terms of goodness-of-fit, mean square error and trackability. We finally apply our proposed filters to multi-unit
spiking data from ferret primary auditory cortex (A1) during
passive stimulus presentation and during performance of a
click rate discrimination task [22] in order to characterize
the spectrotemporal receptive field (STRF) plasticity of A1
neurons. Application of our algorithm to these data provides
new insights into the time course of attention-driven STRF
plasticity, with over 3 orders of magnitude increase in temporal
resolution from minutes to centiseconds, while capturing the
underlying sparsity in a robust fashion. Aside from their theoretical significance, our results are particularly important in
light of the recent technological advances in neural prostheses,
which require real-time robust neuronal system identification
from limited data.
The rest of the paper is organized as follows: In Section
II, we present our notational conventions, preliminaries and
problem formulation. In Section III, we introduce the main
theoretical results of this paper, including the construction
and stability analysis of the objective function, recursive filter
development, and computation of confidence regions. Section
IV provides numerical simulations as well as application to
real data, followed by our concluding remarks in Section V.
Technical details of Section III are presented in Appendices
A–C.
II. P RELIMINARIES

AND

P ROBLEM D EFINITION

We first give a brief introduction to point process models
(see [3] for a detailed treatment). We will use the following notation throughout the paper. Parameter vectors are denoted by
bold-face greek letters. For example, ω = [ω1 , ω2 , · · · , ωM ]′
denotes an M -dimensional parameter vector, with [·]′ denoting
the transpose operator.
Consider a stochastic process defined by a sequence of
discrete events at random points in time, noted by tJ1 =
[t1 , t2 , · · · , tJ ]′ , and a counting measure given by
Z t
J
X
dN (t) =
δ(t − tk ), and N (t) =
dN (u),
k=1

0

where δ(.) is the Dirac’s measure. The Conditional Intensity
Function (CIF) for this process, denoted by λt|Ht , is defined
as
P (N (t + ε) − N (t) = 1|Ht )
,
(1)
λ(t|Ht ) := lim
ε→0
ε
where Ht denotes the history of the process as well as the
covariates up to time t. The CIF can be interpreted as the
instantaneous rate given the history of the process and the
covariates. A point process with a CIF given by λ(t|Ht ) is
defined as:

1) N (0) = 0
2) Given 0 = t0 < t1 < t2 < · · · , the random
variables N (tk ) − N (tk−1 ) are conditionally mutually
independent.
3) For any 0 ≤ t1 < t2 , N (t2 )−N (t1 ) is a Poisson random
variable with probability distribution
R
k R t2
t2


λ(t|H
)dt
e− t1 λ(t|Ht )dt
t
t1
.
P N (t2 ) − N (t1 ) = k =
k!
A point process model is fully characterized by its CIF.
For instance, λ(t|Ht ) = λ corresponds to the homogenous
Poission process with rate λ. A discretized version of this
process can be obtained by binning N (t) within an observation
interval of [0, T ] by bins of length ∆, that is
nt := N (t∆) − N ((t − 1)∆),

t = 1, 2, · · · , T

(2)

where T := ⌈T /∆⌉. Throughout this paper, {nt }Tt=1 will be
considered as the observed spiking sequence, which will be
used for estimation purposes. Also, by approximating Eq. (1)
for small ∆ ≪ 1, and defining λt := λ(t∆|Ht∆ ), we have:
P(nt = 0) = 1 − λt ∆ + o(∆),
P(nt = 1) = λt ∆ + o(∆),
P(nt ≥ 2) = o(∆).

(3)

In discrete time, the orderliness of the process is equivalent
to the requirement that with high probability not more than
one event fall into any given bin. In practice, this can always
be achieved by choosing ∆ small enough. An immediate
consequence of Eq. (3) is that {nt }Tt=1 can be approximated
by a sequence of Bernoulli random variables with success
probabilities {λt ∆}Tt=1 .
A popular class of models for the CIF is given by Generalized Linear Models (GLM). In its general form, a GLM consists of two main components: an observation model (which
is given by (3) in this paper) and an equation expressing some
(possibly nonlinear) function of the observation mean as a
linear combination of the covariates. In neuronal systems, the
covariates consist of extrinsic covariates (e.g., neural stimuli)
as well as intrinsic covariates (e.g., the history of the process).
In this paper, we only consider GLM models with purely
extrinsic covariates, although most of our results can be
generalized to incorporate intrinsic covariates as well.
Let st denote the stimulus at time bin t, [θ0 , θ1 , · · · , θM−2 ]′
denote the vector of stimulus modulation parameters, and µ
denote the baseline firing rate. We adopt a logistic regression
model for the CIF as follows:


M−2
X
λt ∆
θi st−i (4)
=µ+
logit(λt ∆) := log
1 − λt ∆
i=0
By defining ω := [µ, θ0 , θ1 , · · · , θM−2 ]′ and xt
[1, st , · · · , st−M+2 ], we can equivalently write:
λt ∆ = logit−1 (ω ′ xt ) :=

exp(ω ′ xt )
1 + exp(ω ′ xt )

:=

(5)

The model above is also known as the logistic-link CIF
model. Another popular model in the computational neuroscience literature is the log-link model where λt ∆ =
exp(ω ′ xt ). The significance of the logistic-link model is that

3

logit−1 (.) maps the real line (−∞, +∞) to the unit probability
interval (0, 1), making it a feasible model for describing
statistics of binary events independent of the scaling of the
covariates and modulation parameters.
Despite capturing the stimulus dependence in quite a general
form, the GLM model in (5) represents a static model. We
therefore generalize this model to the dynamic setting by
allowing temporal variability of the modulation parameters:
λt ∆ = logit−1 (ω ′t xt ) =

exp(ω ′t xt )
1 + exp(ω ′t xt )

(6)

′

where ω t := [µt , θt,0 , θt,1 , . . . , θt,M−2 ] represents the timevarying parameter vector at time t. Throughout the rest of the
paper, we refer to xt and ω t as the covariate vector and the
modulation parameter vector at time t, respectively.
In our applications of interest, the modulation parameter
vector ω exhibits a degree of sparsity [23], [24]. That is, only
certain components in the stimulus modulation have significant
contribution in determining the statistics of the process. These
components can be thought of as the preferred or intrinsic
tuning features of the underlying neuron. To be more precise,
for a sparsity level L < M , we denote by S ⊂ {1, 2, · · · , M }
the support of the L highest elements of ω in absolute value,
and by ω L the best L-term approximation to ω. We also define
σL (ω) := kω − ω L k1

(7)

to capture the compressibility of the parameter vector ω.
M
Recall
PM that for x ∈ R , the ℓ1 -norm is defined as kxk1 :=
i=1 |xi |. When σL (ω) = 0, the parameter ω is called Lsparse, and when σL (ω) is small compared to kωL k1 , the
parameter is called L-compressible [25].
Finally, the main estimation problem of this paper can
be stated as follows: given binary observations {nt }Tt=1 and
covariates {xt }Tt=−M+1 from a point process with a CIF given
by Eq. (6), the goal is to estimate the M -dimensional Lcompressible parameter vectors {ωt }Tt=1 in an online and
stable fashion.
III. M AIN R ESULTS
In this section, we will first describe the construction of
an appropriate objective function for addressing our main
estimation problem. We will then present a rigorous analysis
of the maximizers of the objective function, which extends
the results of CS to our setting. Next, we will introduce two
adaptive filters to recursively maximize the objective function
based on proximal gradient techniques. Finally, we will outline
how statistical confidence bounds can also be constructed in
a recursive fashion for our estimates.
A. ℓ1 -regularized Exponentially Weighted Maximum Likelihood (ML)
Before proceeding with the construction of the objective
function, we need to introduce more notational conventions. In
order to have a framework allowing multi-timescale dynamics,
we consider piece-wise constant dynamics for the parameter
ω t . That is, we assume that ω t remains constant over windows
of arbitrary length W ≥ 1 samples, for some integer W . By
segmenting the corresponding spiking data {nt }Tt=1 into K :=

T
W

windows of length W samples each, we assume that the
CIF for each time point (k − 1)W + 1 ≤ t ≤ kW is governed
by ωt = ω k , for k = 1, 2, · · · , K. Note that number of spiking
samples K is assumed to be an integer multiple of window
size W , without loss of generality.
Invoking the Bernoulli approximation to the Poisson statistics for ∆ ≪ 1, the log-likelihood of the observation nt at
time t can be expressed as:
log p(nt ) ≈ nt log(λt ∆) + (1 − nt ) log(1 − λt ∆)
= nt (x′t ω t ) − log (1 + exp (x′t ω t )) .

(8)

Assuming conditional independence of the spiking events, the
joint log-likelihood of the observations within window k is
given by:
Lk (ω k ) :=

W
X
j=1

n(k−1)W +j x′(k−1)W +j ω k

− log 1 + exp(x′(k−1)W +j ω k )

(9)

In order to explicitly enforce adaptivity in the log-likelihood
function, we adopt the forgetting factor mechanism of the
RLS algorithm, where the log-likelihood of each window is
exponentially weighted regressively in time, with a forgetting
factor 0 < β ≤ 1. That is, the effective data log-likelihood up
to and including window k is taken to be:
Lβ (ω k ) :=

k
X
i=1

β k−i Li (ω k )

(10)

for some 0 < β ≤ 1. Note what for β = 1, L1 (ω k ) coincides
with the natural data log-likelihood. Moreover, if we replace
the Bernoulli log-likelihood with the Gaussian log-likelihood,
then Lβ (ω k ) coincides with the conventional RLS objective
function.
Next, in order to explicitly enforce sparsity, we adopt the
ℓ1 -regularization mechanism of CS. That is, at window k, we
seek an estimate of the form:
b k = argmax
ω
ωk



Lβ (ω k ) − γkωk k1

(11)

where γ is a regularization parameter controlling the trade off
between the log-likelihood fit and the sparsity of estimated
parameters. Our theoretical analysis in the next subsection
reveals appropriate choices for γ, β and the trade-offs therein.
B. Stability Analysis of the Objective Function
In order to quantify the trade-offs involving our choice of
the objective function in Eq. (11), we proceed in the tradition
of performance analysis result of the RLS algorithm [13] by
characterizing the geometric properties of the estimates ωk in a
stationary environment where ωk = ω for all k. Our analysis,
however, is quite general and avoids ad hoc assumptions
such as direct averaging or covariate independence which are
usually invoked in the analysis of least squares problems.
We need to make the following technical assumptions for
our analysis:
1) The stimulus sequence {st }Tt=−M+1 consists of independent (but not necessarily identically distributed) random variables with a variance of σ 2 which are uniformly

4

bounded by a constant B > 0 in absolute value. Note
that with this assumption, two successive covariate vectors, say at times t and t + 1, given respectively by
xt = [1, st−M+2 , st−M+3 , st−M+4 , · · · , st ] and xt+1 =
[1, st−M+3 , st−M+4 , · · · , st , st+1 ] are highly dependent, as
they have M − 3 random variables in common. Hence,
the independent assumption used in studying least squares
problem is violated.
2) We further assume that for all times t, 0 < pmin ≤
λt ∆ ≤ pmax < 1, for some constants pmin and pmax , i.e., the
probability of spiking does not reach its extremal values of 0
and 1, but can get arbitrarily close. This assumption can be
realized due to the boundedness of the covariates and appropriate normalization of the stimulus modulation coefficients,
and does not result in any practical loss of generality.
We have the following theoretical result regarding the
stability of the maximizers of the objective function:
Theorem 1: Suppose that binary observations from a point
process with a CIF given by Eq. (6) are given over K windows
of length W each. Consider a stationary environment with
ω k = ω for all k and suppose that ω is L-compressible. Then,
under assumptions (1) and (2), for a fixed positive constant
d > 0, there exist constants C, C ′ , and C ′′ such that for
log 2
C′
1 − L2 log
M ≤ β < 1, K ≥ log( β1 ) and a choice of γ ≥
q
M
b
C ′′ log
1−β , any solution ω to (11) satisfies the bound
p
p
p
b −ωk2 ≤ C (1−β)L log M + CσL (ω) 4 (1−β)L log M ,
kω

with probability at least 1 − M5d . The constants C, C ′ , and
C ′′ are only functions of d, pmin , pmax , σ 2 , B, and W , and
are explicitly given in Appendix A.
Proof: The proof is given in Appendix A.

Remarks. The result of Theorem 1 has
√ four major implications. First, the error bound scales with L, the sparsity level,
as opposed to M , the ambient dimension of the parameter
vector, which is consistent with results from CS, and results in
the robustness of the estimate when the underlying parameter
is sparse. Note that the bounds holds for general non-sparse ω,
but is sharpest when σL (ω) is negligible, i.e., the parameter
vector is nearly L-sparse.
Second, the theorem prescribes a lower bound on the
forgetting factor akin to the bounds obtained in CS theory
for the total number of observations. For instance, the result of
[26] for CS under Toeplitz sensing measurements for the linear
model requires T = O(L2 log M ) number of measurements to
achieve a similar scaling of the error bound. In our case, the
role of the number of measurements is transferred to forgetting
1
as the effective length of the measurefactor by taking 1−β
ments. In the absence of the forgetting factor (β = 1), by a
careful limiting process, our results require T = O(L2 log M )
measurements. The latter case can be compared to the result of
[18] for point process models with independent and identically
distributed covariate vectors, which requires O(L log M ) for
stability. The loss of O(L) is incurred due to the shift structure
and hence high dependence of the covariate vectors in our
case, as exemplified in assumption (1).

Third, the theorem reveals the scaling of the regularization
parameter in terms of M and β. In particular, this scaling is
significant as it reveals another role for the forgetting factor
mechanism: not only the forgetting factor mechanism allows
for adaptivity of the estimates, it also controls the scaling of
the ℓ1 -regularization term with respect to the log-likelihood
term. Fourth, unlike conventional results in the analysis of
adaptive filters which concern the expectation of the error in
the asymptotic regime, our result holds for a single realization
with probability polynomially approaching 1, in the nonasymptotic regime.
Note that the objective function is clearly concave, and
assuming that the matrix of the covariate vectors is full-rank,
will be strictly concave with a unique maximizer. However,
the result of Theorem 1 does not require the uniqueness of
the maximizer and holds for any maximizer of the objective
function. In the next section, we will proceed with the development of recursive filters to track the maximizer of the
objective function in the more general time-varying setting.
C. Algorithm Development
Several standard optimization techniques, such as interior
point methods, can be used to find the maximizer of (11).
However, most of these techniques operate in batch mode and
do not meet the real-time requirements of the adaptive filtering
setting where the observations arrive in a streaming fashion. In
order to avoid the increasing runtime complexity and memory
requirements of the batch-mode computation, we seek a recursive approach which can perform low-complexity updates in
an online fashion upon the arrival of new data in order to form
the estimates. To this end, we adopt the proximal gradient
approach. A version of the proximal gradient algorithm is
given in Appendix B. Each iteration of the algorithm moves
the previous iterate along the gradient of the log-likelihood
function, which will then pass through a shrinkage operator.
Before describing further details, we introduce
a more compact notation for convenience. Let
nk := [n(k−1)W +1 , n(k−1)W +1 , . . . , nkW ]′ denote the vector
of observed spikes within
′
 window k, for k = 1, 2, . . . , K.
Similarly, let λk := λ(k−1)W +1 , λ(k−1)W +1 , . . . , λkW
denote the vector of CIFs within window k. By extending
the domain of the logit−1 (·) to vectors in a component-wise
fashion, we can express λk ∆ as:


(12)
λk ∆ = logit−1 Xk ω k

′
where Xk := x(k−1)W +1 , x(k−1)W +2 , . . . , xkW is the data
matrix of size W ×M with rows corresponding to the covariate
vectors in window k. Suppose that at window k, we have an
(ℓ)
b k , for ℓ = 0, 1, · · · , R, with R being an
iterate denoted by ω
integer denoting the total number of iterations. The gradient
(ℓ)
b k can be written as:
of Lβ (·) evaluated at ω
k





 X
(ℓ)
(ℓ)
(ℓ)
bk
b k =: gk ω
bk =
(13)
β k−i X′i εi ω
∇ω Lβ ω
i=1

where εi (·) := ni − λi (·)∆ represents the innovation vector
of the point process at window i. The innovation vector εi can
be thought of as the counterpart of the conventional innovation

5

vector in adaptive filtering of linear models. The proximal
gradient iteration for the ℓ1 -regularization can be written in
the compact form as:



(ℓ)
(ℓ)
(ℓ+1)
bk
b k + αgk ω
bk
(14)
ω
= Sγα ω

where Sτ (·) is the element-wise soft thresholding operator at a
level of τ given in Appendix B. The final estimate at window
k is obtained following the Rth iteration, and is denoted by
(R)
b k := ω
b k . In order to achieve a recursive updating rule for
ω
gk , we can rewrite Eq. (13) as:






(ℓ)
(ℓ)
(ℓ)
bk
b k + X′k εk ω
b k = β gk−1 ω
(15)
gk ω

However, in an adaptive setting, we only have access to values
(1:L)
b k−1 (·)! In order to turn Eq. (15) into
of gk−1 evaluated at ω
a fully recursive updating rule, all the previous CIF vectors
k−1
{λi (·)}i=1
should be recalculated at the most recent set of
(1:L)
b
iterates ω k (·). In order to overcome this computational
burden, we exploit the smoothness of the logistic function and
employ the Taylor series expansion of the CIF to approximate
the required recursive update. In what follows, we consider
the zeroth order and first order expansions, which result in
two distinct, yet fully recursive, updating rules for Eq. (15).
Zeroth Order Expansion: By retaining only
 thefirst term in
(ℓ)
bk
b i,
around ω
the Taylor series expansion of the CIF λi ω


we get:
(ℓ)
b k ∆ ≈ λi (ω
b i) ∆
(16)
λi ω

)
b i ). Substituting this approximation
where λi ∆ = logit−1 (Xi ω
in Eq. (13), we can express the zeroth order approximation to
the gradient at window k, denoted by gk0 (·), as:
k

 X
(ℓ)
bk =
β k−i X′i εi (b
ωi )
(17)
gk0 ω
i=1

It is then straightforward to obtain a recursive form as:






(ℓ)
(ℓ)
(ℓ)
0
bk
b k = β gk−1
b k + X′k εk ω
(18)
gk0 ω
ω

The shrinkage step will be then given by:



(ℓ)
(ℓ)
(ℓ+1)
bk
b k + αgk0 ω
bk
ω
= Sγα ω

(19)

We refer to the resulting filter as the ℓ1 -regularized Point
Process Filter of the Zeroth Order (ℓ1 -PPF0 ). A pseudo-code
is given in Algorithm 1.

First Order Expansion: If instead, we retain the first two
terms in the Taylor expansion, Eq. (16) will be replaced by:




(ℓ)
(ℓ)
bk − ω
bi
b k ∆ ≈ λi (b
(20)
ω i ) ∆ + Λi (b
ωi ) Xi ω
λi ω

where Λi (b
ω i ) is a diagonal W × W matrix with the (m, m)th
diagonal element given by λ(i−1)W +m ∆(1 − λ(i−1)W +m ∆).
Using the first order approximation above, we can improve the
resulting approximation to the gradient, denoted by gk1 , as:
k

 X


(ℓ)
(ℓ)
1
bk − ω
bi
bk =
β k−i X′i εi (b
ωi ) − Λi (b
ω i )Xi ω
gk ω
i=1

By defining:
(ℓ)

uk :=
(ℓ)

k
X

i=1
k
X



bi
β k−i X′i εi (b
ω i ) + Λi (b
ω i )Xi ω

β k−i X′i Λi (b
ω i )Xi ,

i=1 
(ℓ)
bk
we can express gk1 ω
as:


(ℓ)
(ℓ)
(ℓ) (ℓ)
b k = uk − B k ω
bk .
gk1 ω
Bk :=

(0)

b k , and R .
Inputs: nk , Xk , gk−1 , ω
1: for ℓ = 0, . . . , R −1 do

(ℓ)
bk
2:
λk ∆ = logit−1 Xk ω
3:
ε k = nk − λ k ∆
′
4:
gk = β gk−1 +
i
h Xk εk
(ℓ)
(ℓ+1)
b k + αgk
bk
5:
ω
= Sγα ω
6: end for
(R)
b k := ω
bk .
Output: ω

(22)

(23)

It is then straightforward to check that both uk and Bk can
be updated recursively [19] as:




 
(ℓ)
(ℓ)
(ℓ)
(ℓ)
(R)
bk
b k Xk ω
b k + Λk ω
uk = β uk−1 + X′k εk ω


(ℓ)
(ℓ)
(R)
b k Xk
Bk = β Bk−1 + X′k Λk ω
(ℓ)

(ℓ)

Note that the update rules for both Bk and uk involve
simple rank-W operations. The shrinkage step is then given



by:
(ℓ)
(ℓ)
(ℓ)
(ℓ+1)
b k + αuk
bk
(24)
ω
= Sγα I − αBk ω

We refer to the resulting filter as the ℓ1 -regularized Point
Process Filter of the First Order (ℓ1 -PPF1 ). A pseudo-code
is given in Algorithm 2.
Remark. The computational complexity of ℓ1 -PPF0 and ℓ1 PPF1 algorithms can be shown to be linear and quadratic
in M per iteration, respectively. Our results in Section IV
will reveal that both filters outperform existing filters of the
same complexity, respectively. Furthermore, ℓ1 -PPF1 exhibits
superior performance over ℓ1 -PPF0 as expected, although with
a cost of O(M ) in computational complexity per iteration.
Algorithm 2 ℓ1 -regularized Point Process Filter of the First Order
(ℓ1 -PPF1 )
(0)

Algorithm 1 ℓ1 -regularized Point Process Filter of the Zeroth Order
(ℓ1 -PPF0 )

(21)

b k , and R .
Inputs: nk , Xk , uk−1 , Bk−1 , ω
1: for ℓ = 0, . . . , R −1 do

(ℓ)
bk
2:
λk ∆ = logit−1 Xk ω
3:
εk = nk − λk ∆
4:
(Λk )m,m = (λk )m ∆ (1 − (λk )m ∆), 
m = 1, · · · , W
(ℓ)
′
bk
5:
uk = β uk−1 + Xk εk + Λk Xk ω
6:
Bk = β Bk−1 + X′k Λk Xk
(ℓ)
bk
7:
gk = uk − Bkhω
i
(ℓ)
(ℓ+1)
b k + αgk
bk
8:
ω
= Sγα ω
9: end for
(R)
b k := ω
bk .
Output: ω

6

D. Constructing Confidence Intervals
Characterizing the statistical confidence bounds for the
estimates is of utmost importance in neural data analysis, as
it allows to test the validity of various hypotheses. Although
construction of confidence bounds for linear models in the
absence of regularization is well understood and widely applied, regularized ML estimates are usually deemed as point
estimates for which the construction of statistical confidence
regions is not straightforward. A series of remarkable results in
high-dimensional statistics [27]–[29] have recently addressed
this issue by providing techniques to construct confidence
intervals for ℓ1 -regularized ML estimates of linear models as
well as GLMs. These approaches are based on a careful inspection of the Karush-Kuhn-Tucker (KKT) conditions for the
regularized estimates, which admits a procedure to decompose
the estimates into a bias term plus an asymptotically Gaussian
term (referred to as ’de-sparsifying’ in [28]), which can be
computed using a nodewise regression [30] of the covariates.
In what follows, we give a brief description of how the
methods of [28] apply to our setting, and leave the details to
Appendix C. Using the result of [28], the estimate ω k as the
maximizer of (11) can be decomposed as:
b k gk (b
bk = Θ
bk
ω
ωk ) + w

(25)

b k is an approximate inverse to the Hessian of Lβ (ω)
where Θ
b k , gk is the gradient of Lβ (ω) previously
evaluated at ω
b k is an unbiased and asymptotdefined in Eq. (13), and w
ically Gaussian random vector with a covariance matrix of
b k Gk (b
b ′ , with
b k) = Θ
cov(w
ω k )Θ
k
Gk (b
ω k ) :=

k
X

β 2(k−i) X′i εi (b
ω k )εi (b
ω k )′ Xi .

(26)

i=1

The first term in Eq. (25) is a bias term which can be directly
b k . Given cov(w
b k ), statistical confidence
computed given Θ
bounds for the second term at desired levels can be constructed
Algorithm 3 Recursive Construction of the Confidence Reb k .a
gions for the mth Component of w

b (0) .
b k , Gk−1 , m, γm , and ψ
Inputs: nk , Xk , uk , Bk , and ω
m
bk
1: gk = uk − Bk ω
2: Gk = β 2 Gk−1 + X′k εk ε′k Xk
3: for ℓ = 0, . . . , R −
i
h 1(ℓ)do 
b −α (Bk )m,\m−(Bk )\m,\m ψ
b (ℓ)
b (ℓ+1) = Sγ α ψ
4:
ψ
m
m
m
m
5: end for
2
b (R) (Bk )′
6: τm
= (Bk )m,m − ψ
m
m,\m
(R)
b
7: (c)m = 1, and (c)\m = −ψ m
b k )m = 12 c
8: (Θ
τm

2
b k )m Gk (Θ
b k )′m
σ
bk,m
:= (Θ
b k )m gk
b k )m = (b
10: (w
ω k )m − (Θ
b k )m ± Φ−1 (1 − α/2) σ
Output: CRk,m := [(w
bk,m ]

9:

a For a matrix A ∈ RM ×M , we denote by (A)
m,\m the mth row with
the mth element removed, and by (A)\m,\m the submatrix of A with both
the mth row and column removed.

in a standard way. The main technical issue in the aforemenb k in a
tioned procedure in our setting is the computation of Θ
b
recursive fashion. Since the rows of Θk are computed using
ℓ1 -regularized least squares, we use the SPARLS algorithm
[19] as an efficient method to carry out the computation in
a recursive fashion. Algorithm 3 summarized the recursive
computation of confidence intervals for the mth component
b k.
of w
IV. A PPLICATIONS
In this section, we will apply the proposed algorithms to
simulated data as well as real spiking data from the ferret
primary auditory cortex. In our simulation studies, we compare
the performance of our proposed filters with two of the stateof-the-art point process filters, namely the steepest descent
point process filter (SDPPF) [7] and the stochastic state point
process filter (SSPPF) [14]. These adaptive filters are based
on approximate solutions to the Chapman-Kolmogorov forward equation obtained by a steepest descent and a Gaussian
approximation procedure, respectively.
A. Simulation Study 1: MSE and Sparse Recovery Learning
Curves
First, we consider a stationary environment where ω is constant over time. We use a bin size of ∆ = 1 ms and window
size of W = 1 sample, for a total observation window of
T = 30 sec (K = 30000). The length of the parameter vector
ω = [µ, θ] is chosen as M = 101. For each realization, we
draw a sparse parameter vector θ of fixed length M −1 = 100
and sparsity L = 3. The support S and values of the nonzero
components of θ are chosen randomly and the values are
normalized so that kθk2 = 10. The stimulus input sequence
{sk }K
k=−M+1 is drawn from an i.i.d. Gaussian distribution
N (0, σ 2 ). The binary spike train {nk }K
k=1 is generated as a
single realization of conditionally independent Bernoulli trials
with success rate λk ∆. The stimulus variance σ 2 is chosen
as σ 2 = 0.01 small enough so that the average spiking rate
λavg ∆ = 0.13 ≪ 1 to ensure that the Bernoulli approximation
is valid. All the simulations are done with R = 1 iteration per
time step. The step size α is chosen as α ≃ 9 × 10−4 (See
Appendix B for details).
For a given pair of (α, β) parameters, we select an optimal
value for the regularization parameter γ by performing a twofold even-odd cross validation procedure: first, the data are
split into two sets of even and odd samples in an interleaved
manner. Then, one set is used as the training set for estimation
of the parameter vector ω k and the other is used to assess the
b k with respect to the loggoodness-of-fit of the estimates ω
likelihood of the observations. We repeat the process switching
the role of the two sets and take the average as the overall
measure of fit.
b denote the averaging operator with respect to realizaLet E
tions. We consider two performance metrics: the normalized
b ω
bk −
mean squared error (MSE) defined as MSEk := Ek
2 b
2
ωk k /Ekωk k to evaluate MSE performance at time step k;
bk −
b θ
and the out-of-support energy defined as SPMk := Ek
2
2
bk )S k /Ekθ
bk k to represent a sparsity metric (SPM). Ideally,
(θ
SPMk must be equal to zero at all times. The averaging is
carried out over a sufficiently large number of runs.

7

MSE (dB)

A 0
-4
-8
-12

SPM

B 1

0
0

20

10

30

Time (s)
Fig. 1. Learning curves of the adaptive filtering algorithms in a stationary
environment. A) MSE vs. time, B) SPM vs. time.

Figure 1 shows the corresponding learning curves for
the four algorithms. According to Figure 1–A, the ℓ1 -PPF1
achieves the lowest stationary MSE measure of −10.4 dB,
followed ℓ1 -PPF0 which achieves an MSE of −9.2 dB. The
SSPPF and SDPPF algorithms respectively achieve an MSE
of −3.35 dB and −1.9 dB, which reveals a gap of ≈ 7 dB
with respect to our proposed filters. Note that this result is
consistent with the prediction of Theorem 1 and highlights
the MSE gain achieved by ℓ1 -regularization, as opposed to
ML, when the underlying parameter is sparse.
B. Simulation Study 2: Tracking and Goodness-of-fit Performance
In the second simulation scenario, we consider a more
realistic setting where ωk evolves in time. Furthermore, as in
the case of real data applications, we assume that the support
of ω k is not available as a performance benchmark and resort
to statistical goodness-of-fit test. These tests for point process
models have been developed as an application of the timerescaling theorem [31], [32] and consist of the KolmogorovSmirnov (KS) test for assessing the conditional intensity
estimation accuracy, and the Autocovariance Function (ACF)
test to assess the conditional independence assumption. We
skip the details in the interest of space, and refer the readers
to the aforementioned references for a detailed treatment.
As in the previous case, we consider a bin size of ∆ = 1ms,
window size of W = 1, and a total observation window of
T = 60sec (K = 60000 bins). The stimulus is generated as
in the previous case. For the parameter vector ω k , we choose
a fixed baseline rate of µk = −2.51 to set the baseline spiking
rate to λavg ∆ ≈ 0.1, and select a sparse modulation vector
θk of length M = 100 with a support S = {1, 10, 20} of size
L = 3, and respective values of (θk ){1,10,20} = {10, −5, 5}
for k ≤ K/2. Halfway through the test, at k = K/2 + 1,
the largest component (θ k )1 , drops rapidly and linearly to 0,
within a window of length 1 sec and remains zero for the rest
of the run.
Figure 2 shows the performance of all four algorithms in
the aforementioned setting. Each row (A through D) shows the
true time-varying parameter vector (dashed traces) as well as
the filtered estimates (solid traces) in the left panel. In particular, the gray solid traces show the out-of-support components
which must ideally be equal to zero. The colored hulls around
bk )1 show the 95% confidence intervals (note that confidence
(θ

intervals for SDPPF cannot be directly obtained and require
averaging over multiple realizations). The middle and right
panels show the KS and ACF test results at a 95% confidence,
respectively. For the quadratic algorithms ℓ1 -PPF1 and SSPPF,
a forgetting factor of β = 0.9995 is chosen. The regularization
parameter for ℓ1 -PPF1 is chosen as γ = 0.5, obtained by the
aforementioned two-fold even-odd cross validation. For the
first order algorithm ℓ1 -PPF0 , a smaller forgetting factor of
β = 0.995 is chosen to ensure stability, and a value of γ = 0.1
is used based on cross validation. These settings ensure that
all the algorithms are tuned in their optimal operating point
for fairness of comparison.
Figure 2–A and 2–B reveal three striking performance
gaps between the two second-order algorithms (with the same
computational complexity, quadratic in M ): first, the out-ofsupport components (gray traces) of ℓ1 -PPF1 are significantly
smaller than those of SSPPF; second, the confidence regions
of ℓ1 -PPF1 are narrower than those of SSPPF; and third, ℓ1 PPF1 fully passes the KS test, while SSPPF marginally does
so. Similarly, comparing the two first order algorithms (with
the same computational complexity, linear in M ) Figure 2–C
and 2–D reveal that the ℓ1 -PPF0 significantly suppresses the
out-of-support components as compared to SDPPF. Moreover,
ℓ1 -PPF0 provides confidence bounds, which cannot be directly
obtained for SDPPF. Finally, ℓ1 -PPF0 marginally fails the KS
test, whereas SDPPF does so significantly. Both algorithms fail
the ACF test, which shows that the second-order corrections
embedded in ℓ1 -PPF1 and SSPPF is necessary to achieve a
better goodness-of-fit, which a price of higher computational
complexity.
We also inspect the estimated firing probability λk (b
ω k )∆
for the four algorithms in Figure 3. In addition, we include
the probability estimated by the normalized reverse correlation
(NRC) method, which is commonly used in neural data
analysis, and fits the modulation parameters using a linear
model. Figure 3 shows the true spiking probability (blue
solid trace) and the resulting spikes (black vertical lines). In
the subsequent rows (B through F), the true and estimated
probabilities are shown by dashed blue and solid red traces,
respectively. A comparison of all the rows reveals that ℓ1 -PPF1
and ℓ1 -PPF0 outperform SSPPF and SDPPF, respectively, in
terms of estimating the true probability. The NRC method
is inferior to the preceding four algorithms, and results in
negative estimates of the probability due to its use of a linear
model (as opposed to logistic).
C. Application to Real Data: Dynamic Analysis of Spectrotemporal Receptive Field Plasticity
The responses of neurons in the primary auditory cortex
(A1) can be characterized by their spectrotemporal receptive
fields (STRFs), where each neuron is tuned to a specific region
in the time-frequency plane, and only significantly spikes when
the acoustic stimulus contains spectrotemporal contents matching its tuning region [2] (See Figure 4, top row, leftmost panel).
Several experimental studies have revealed that receptive fields
undergo rapid changes in their characteristics during attentive
behavior in order to capture salient stimulus modulations [22],
[33], [34]. In [22], it is suggested that this rapid plasticity

8

1

1

10
5
0
-5

ACF

Emp. Quantiles

A

0
0
1
1

10
5
0
-5

ACF

Emp. Quantiles

B

0
0
1

1

10
5
0
-5

ACF

Emp. Quantiles

C

0
0
1
1

10
5
0
-5

ACF

Emp. Quantiles

D

0
0

0

30
Time (s)

60

0

1
True Quantiles

0

50
Lag (ms)

Fig. 2. Performance comparison of the adaptive filtering algorithms: A) ℓ1 -PPF1 , B) SSPPF, C) ℓ1 -PPF0 , and D) SDPPF. In each row, the left panel shows
the trie parameter vector with dashed traces and the estimates with solid traces. Colored hulls show the 95% confidence intervals for one of the components.
The middle and right panels show the corresponding KS and ACF test plots, respectively. Red traces show confidence regions at a level of 95% for both tests.

A

0.4

B

0
0.4

C

0
0.4

D

0
0.4

E

0
0.4

F

0
0.3

-0.1
34.2

has a significant role in the functional processes underlying
active listening. However, most of the widely-used estimation
techniques (e.g., normalized reverse correlation) provide static
estimates of the receptive field with a a temporal resolution
of the order of minutes. Moreover, they do not systematically
capture the inherent sparsity manifested in the receptive field
characteristics.

34.3

34.4

Time (s)

Fig. 3.
Firing rate estimates for adaptive filtering algorithms within an
interval of 34.2 s ≤ t ≤ 34.4 s: A) true rate (blue solid trace) and spikes
(black vertical lines), B) ℓ1 -PPF1 , C) SSPPF, D) ℓ1 -PPF0 , E) SDPPF, and
F) normalized reverse correlation (NRC). In rows B through F, the dashed
blue traces and solid red traces show the true rate and the estimated rate,
respectively.

In the context of our model, the STRF can be modeled as an
(I ×J)-dimensional matrix, where I and J denote the number
of time lags and frequency bands, respectively. By vectorizing
this matrix, we obtain an (M − 1)-dimensional vector θ k at
window k, where M = I ×J +1. Augmenting the baseline rate
parameter µk , we can model the activity of the A1 neurons
using the logistic CIF with a parameter ω k := [µk , θk ]′ .
The stimulus vector at time t, st is given by the vectorized
version of the spectrogram of the acoustic stimulus with J
frequency bands and I lags. In order to capture the sparsity of
the STRF in the time-frequency plane, we further represent θk
over a Gabor time-frequency dictionary consisting of Gaussian
windows centered around a regular subset of the I × J timefrequency plane. That is, for θk = Fξk , where F is the
dictionary matrix and ξ k is the sparse representation of the
STRF. The estimation procedures of this paper can be applied
to ξ k , by absorbing the dictionary matrix into the data matrix
Xk at window k.

9

Frequency (Hz)

16k

0.02

B

8k

E
A

4k

C
0.01

D

0

2k

-0.01

1k

-0.02

0.5k

10

20

30

40

-0.03

50

STRF Components

Time lag (ms)
0.02

A
B
C

0

D

-0.02

E

-0.04

0

100

200

300

400

500

600

700

800

900

Time (s)
Fig. 4. The time-course of task-dependent STRF plasticity of a ferret A1 neuron. The top row shows snapshots of the STRF at five selected points in time,
marked by the dashed vertical lines in the bottom graph. The bottom graph shows the time-course of five selected points (A through E) in the STRF marked
on the leftmost panel of the top row.

We apply our proposed adaptive filter ℓ1 -PPF1 to multi-unit
spike recordings from the ferrets A1 during a series of passive
listening conditions and active auditory task conditions. During each active task, ferrets attended to the temporal dynamics
of the sounds, and discriminated the rate of acoustic clicks
[33]. The STRFs were estimated from the passive condition,
where the quiescent animal listened to a series of broadband
noise-like acoustic stimuli known as Temporally Orthogonal
Ripple Combinations (TORC). The experiment consisted of 2
active and 11 passive blocks. Within each passive block, 30
TORCs were randomly repeated a total of 4-5 times each. In
our analysis, we pool the spiking data corresponding to the
same repeated TORC within each block. Therefore, the time
axis corresponds to the experiment time modulo repetitions
within each block. We discretize the resulting duration of
T = 990s to time bins of size ∆ = 1 ms, and segment
data to windows of size W = 10 samples (10 ms). The
STRF dimensions are 50 × 50, regularly spanning lags of
1 to 50 ms and frequency bands of 0.5 kHz to 16 kHz
(in logarithmic scale). The dictionary F consists of 13 × 13
Gabor atoms, evenly spaced within the STRF domain. Each
atom is a two-dimensional Gaussian kernel with a variance of
D2 /4 per dimension, where D denotes the spacing between
the atoms. We selected a forgetting factor of β = 0.9998, a
4(1−β)
2
step size of α = MW
σ̄2 , where σ̄ is the average variance of
the spectrogram components, R = 1 iterations per sample, and
a regularization parameter of γ = 40 via two-fold even-odd
cross validation.
Figure 4, top row, depicts five snapshots taken at
{180, 360, 540, 630, 990} sec corresponding to the end-points
of the {2, 4, 6, 7, 11}th passive tasks. The bottom row shows
the time-course of five selected points (marked as A through
D in the leftmost panel of the top row) of the STRF during
the experiment. The STRF snapshots at times 180 and 540 sec

correspond to 90 secs after the two active tasks, respectively,
and verify the sharpening effect of the excitatory region
(∼ 30 msec, 8 kHz) due to the animal’s attentive behavior
following the active task reported in [22]. Moreover, the STRF
snapshots at times 360 and 630 sec reveal the weakening of
the excitatory region long after the active task and returning to
the pre-active state, highlighting the plasticity of A1 neurons.
Previous studies have revealed the STRF dynamics with a
resolution of the order of minutes [34]. Our result in Figure
4 provides a temporal resolution of the order of centiseconds
(3 orders of magnitude increase), while capturing the STRF
sparsity in a robust fashion.
V. C ONCLUDING R EMARKS
In this paper, we considered recursive estimation of the
time-varying parameter vectors in a logistic regression model
for binary time series driven by continuous input. To this end,
we integrated several techniques from compressed sensing,
adaptive filtering, optimization and statistics. We constructed
an objective function which enjoys from the trackability
features of the RLS-type algorithms, sparsifying features
of ℓ1 -minimization, and unlike the rate-based linear models
commonly used to analyze spiking data, takes into account
the binary statistics of the observations. We analyzed the
maximizers of the objective function in a rigorous fashion,
revealing novel trade-offs between various model parameters.
We constructed two adaptive filters, with respective linear and
quadratic complexity requirements, for recursive maximization
of the objective function in an online setting. Moreover,
we characterized the statistical confidence regions for our
estimates, and devised a recursive procedure to compute them
efficiently.
Although we specialized our treatment to logistic statistics
and ℓ1 -regularization, our approach to algorithm development
has a plug-and-play feature: other GLM link functions (e.g.,

10

log-link) with possibly history dependent covariates and other
regularization schemes (e.g., re-weighted ℓ1 , or group-sparse
regularization) can be used instead and result in a large class of
adaptive filters for sparse point process regression. We tested
the performance of our algorithms on simulated as well as
experimentally recorded spiking data. Our simulation studies
revealed that the proposed filters outperform several existing
point process filters. Application of our filters to real data from
the ferret primary auditory cortex provided a high-resolution
characterization of the time-course of spectrotemporal receptive field plasticity, with 3 orders of magnitudes increase
in temporal resolution. Although we focused on auditory
neurons, we expect a similar superior performance of our
filters when applied to other sensory or motor neurons (e.g.,
cells in primary or supplementary motor cortex [35]).
A PPENDIX A
P ROOF OF T HEOREM 1
The proof is mainly based on the beautiful treatment of
Negahban et al. [18]. The major difficulty in our case lies
in the high inter-dependence of the covariates, which form
a Toeplitz structure due to the setup of adaptive filtering.
We address the latter issue by adopting techniques from
another remarkable paper by Haupt et al. [26] to deal with the
underlying interdependence. In the process, we also employ
concentration inequalities for dependent random variables due
to van de Geer [36].
In order to proceed, we adopt the notion of Strong Restricted
Convexity (RSC) introduced in [18]. For a twice differentiable
log-likelihood with respect to ω, the RSC property or order
L implies the existence of a lower quadratic bound on the
negative log-likelihood:
DL (∆, ω) := −Lβ (ω+∆)+Lβ (ω)+∆′ ∇Lβ (ω) ≥ κk∆k22 ,
for a positive constant κ > 0 and all ∆ ∈ RM satisfying:
k∆S c k1 ≤ 3k∆S k1 + 4σS (ω).

(27)

for any index set S ⊂ {1, 2, · · · , M } of cardinality L. The
following key lemma establishes the RSC for Lβ (ω):

Lemma 1: Let {xt }KW
t=1 denote a sequence of covariates
and let ω denote the corresponding logistic parameters. Then,
for a fixed positive constant d > 0, there exist constants C ′
log 2
C′
and κ > 0 such that for β ≥ 1 − L2 log
M and K ≥ log( 1 )
β
the negative log-likelihood −Lβ (ω) satisfies the RSC of order
κ
with probability greater than 1 − M3d .
L with constant 1−β
′
The constants C and κ are only functions of d, pmin , pmax ,
σ 2 , B, W , and are explicitly given in the proof.
Proof: The proof is inspired by the elegant treatment of
Negahban et al. [18]. The major difficulty in our setting is
the high interdependence of successive covariates due to the
shift structure induced by the adaptive setting, whereas in [18],
the matrix of covariates is composed of i.i.d. rows. Using the
Taylor’s theorem, DL (∆, ω) can be written as:

K X
W
X
exp x′(i−1)W +j ω ⋆ |∆′ x(i−1)W +j |2
K−i
, (28)
β

2
′
⋆
i=1 j=1
ω
1 + exp x(i−1)W +j

with ω ⋆ = ω + τ ∆ for some τ ∈ (0, 1). Since by hypothesis
0 < pmin ≤ λi ∆ ≤ pmax < 1, we have:

exp x′(i−1)W +j ω ⋆
2 ≥ pmin(1 − pmax ).

(29)
′
1+exp x(i−1)W +j ω ⋆

We can therefore further lower bound DL (∆, ω) by:

DL (∆, ω) ≥ pmin (1 − pmax )σ 2 Nβ ∆′ Cβ ∆ ,

(30)

K+1

where Nβ := W 1−β
1−β , and
Cβ :=

K W
1 X X K−i
β
x(i−1)W +j x′(i−1)W +j .
σ 2 Nβ i=1 j=1

(31)

Note that the matrix Cβ has highly inter-dependent elements
due to the Toeplitz structure in the adaptive design. In order to
establish the RSC condition, we show the stronger Restricted
Eigenvalue (RE) property, which in turn implies RSC [37]. Let
δ ∈ (0, 1) be fixed. To do so, we need to bound the eigenvalues
of (Cβ )S , the restriction of Cβ to a subset of columns and
rows corresponding to indices in S ⊂ {1, 2, · · · , M } with
|S| = rL, for some integer r > 1 + 9(1+δ)
1−δ .
Without loss of generality, we assume that the first entry
of the covariate vectors xt is replaced by σ instead of 1,
for presentational simplicity of the following treatment. For
m, m′ 6= 1, we have:
(Cβ )m,m′ =

K W −1
1 X X K−i
β
s(i−1)W +j+m−M
σ 2 Nβ i=1 j=0

× s(i−1)W +j+m′ −M .

For m = m′ = 1,
(Cβ )1,1 =

K
1 X
1
1 − β K+1
W β K−i σ 2 =
W
= 1,
2
σ Nβ i=1
Nβ
1−β

and for m 6= 1,

(Cβ )m,1 = (Cβ )1,m =

K W −1
1 X X K−i
β
s(i−1)W +j+m−M .
σNβ i=1 j=0

We also have E{(Cβ )m,m′ } = δmm′ . Using Hoeffding’s
inequality [38] we get:
!
2t2 Nβ2
P (|(Cβ )m,m − 1| > t) ≤ 2 exp −
P
2(K−i)
B4 K
i=1 W β
!
2Nβ2 t2 σ 4
= 2 exp − 4
B Nβ 2


2Nβ t2 σ 4
≤ 2 exp −
,
(32)
B4
K+1

since Nβ 2 = Nβ 1+β
1+β

≤ Nβ , for β ∈ [0, 1]. Similarly,


2Nβ t2 σ 2
P (|(Cβ )1,m | > t) ≤ 2 exp −
,
B2

By adopting the partitioning technique of Theorem 4 in [26],
we also get for m 6= m′ :


N β t2 σ 4
P (|(Cβ )m,m′ | > t) ≤ 4 exp −
.
8B 4

11

4

Let B0 := max{B 2 , B
Proposition 1: Consider a sequence of σ-fields F0 ⊂ F1 ⊂
σ2 }. Now, using the union bound we
have:
· · · . Suppose that Xi is Fi -measurable with |Xi | ≤ Bi for
!

 some constant Bi , i = 1, 2, · · · and E{Xi |Fi−1 } = 0. Then
M 
[
Nβ δ 2 σ 2
δ
≤ 2M 2 exp −
, for all t > 0,
P
(Cβ )m,m′ >
2rL
32B0 r2 L2
!


n
m,m′ =1
X
2t2
m6=m′
.
Xi ≥ t ≤ exp − Pn
P
2
 M2
i=1 Bi
M
i=1
where we have used 2 < 2 and,
Proof: This result is a special case of Theorem 2.5 of
!


M 
[
δ
Nβ δ 2 σ 2
[36] for bounded and possibly dependent random variables,
|(Cβ )m,m − 1| >
P
≤ 2M exp −
.
which generalizes Hoeffding’s inequality.
2rL
4B0 r2 L2
m=1
In
our
case,
we
can
similarly
show
that
Now, by invoking the Gershgorin’s disc theorem, the eigneval- E [s {n − λ ∆} |F ] = E [s E [{n − λ ∆} |F , F ]] =
t
t
t
t−1
t
t
t
t−1
t
ues of any sub-matrix of Cβ restricted to an index set S with 0. Moreover, each summand is bounded by 2β K−i B. Hence,
|S| = rL, lie in the interval [1 − δ, 1 + δ] with probability at using the result of Proposition 1, by taking n = T W and
least:
Xi = si [ni − λi ∆], we get:




2 2
2 2
!
N
δ
σ
N
δ
σ
β
β
1 − 2M 2 exp −
− 2M exp −

t2 Nβ2

2
2
2
2
β
32B0 r L
4B0 r L
P ∇L (ω) m > tNβ ≤ 2 exp − PK


2 i=1 W β 2(K−i)
2 2
N
δ
σ
β
2


.
(33)
≥ 1 − 3M exp −
N β t2
32B0 r2 L2
≤ 2 exp −
.
2
32B0 σ2 r 2 (d+2) 2
Hence, by choosing Nβ ≥
L log M , the probδ2
Using the union bound, we have:
ability above is greater than 1 − M3d .



N β t2
Next, by invoking Lemma 4.1 (ii) of [37], we have that Cβ
.
(37)
P ∇Lβ (ω) ∞ > tNβ ≤ 2M exp −
satisfies the RSC condition over the set given by Eq. (27) with
2
q
a constant given by:
2(d+1) log M
By
choosing
t
=
, we have that
q
2

Nβ
1+δ
(34)
κ0 = (1 − δ) 1 − 3 (r−1)(1−δ) .
q
(38)
∇Lβ (ω) ∞ < 2(d + 1)Nβ log M ,
Hence, the negative log-likelihood satisfies the RSC with a
constant given by pmin (1 − pmax )σ 2 Nβ κ0 . Finally, by taking with probability at least 1 − M2d .
log 2
W
, which makes κ
K ≥ log
, we have that Nβ ≥ 2(1−β)
Hence, for a fixed δ < 1, d > 0, and r > 1 + 9(1+δ)
1−δ ,
( β1 )
W δ2
C′
′
independent of K and β, given by:
by taking β ≥ 1 − L2 log M with C := 64B0 σ2 r2 (d+2) and
q
p
2
M
′′
pmin (1 − pmax )σ κ0 W
b
4(d + 1)W , any maximizer ω
γ > C ′′ log
1−β with C :=
,
(35)
κ :=
2
satisfies:
p
p
p
C′
W δ2
′
and β ≥ 1 − L2 log
b −ωk2 ≤ C (1−β)L log M + CσL (ω) 4 (1−β)L log M ,
kω
M with C := 64B0 σ2 r 2 (d+2) .
Next, the result of Theorem 1 of [18] implies:
r
√
2γ L
2γσL (ω)
b − ωk2 ≤
+
,
kω
κ
κ

for γ > 2k∇Lβ (ω)k∞ . We have, for m 6= 1,

K X
W
X

β K−i s(i−1)W +j+m−M+1
∇Lβ (ω) m =



i=1 j=1

(36)

with probability at least 1 − M3d − M2d , where C is given by
p
128(d + 1)
C := √
2 .

q
1+δ
W pmin (1 − pmax )σ 2 (1 − δ) 1 − 3 (r−1)(1−δ)
A PPENDIX B
T HE P ROXIMAL G RADIENT A LGORITHM

In this appendix, we give an overview of the proximal
n(i−1)W +j+m−M+1 − λ(i−1)W +j+m−M+1 ∆ . gradient algorithm for minimization of convex function. The
Now, let Ft be the σ-field generated by s−M+1 , · · · , st , i.e., corresponding algorithm for maximization of concave functions can be obtained by negating the objective functions.
σ(s−M+1 , · · · , st ). We have that
Consider the general optimization problem
E [{nt − λt ∆} st ] = E [E [{nt − λt ∆} st |Ft ]]
min f (x) + g(x),
(39)
x
= E [st E [{nt − λt ∆} |Ft ]]
i
h 

where functions f (x) : RM → R and g(x) : RM → R ∪ {∞}
= E st E {λt ∆ − λt ∆} |Ft = 0.
{z
}
|
are assumed to be closed proper convex functions. Suppose
=0
that f is differentiable with a Lipschitz continuous gradient


Hence for all m, E ∇Lβ (ω) m = 0. We next invoke ∇f with constant L(∇f ). The function g can be possibly
the following result for concentration of dependent random non-smooth. A wide range of practical optimization problems
variables:
can be cast in this form, particularly in the context of machine

12

learning [39], where the objective function can be decomposed
into a loss function and a regularization term.
The proximal gradient algorithm provides an iterative procedure for solving (39) in the following form:
h
i
(40)
x(ℓ+1) = Pα(ℓ) g x(ℓ) − α(ℓ) ∇f (x(ℓ) ) ,

where the parameter α(ℓ) is an appropriately chosen step size
1
at iteration ℓ so that α(ℓ) < L(∇f
) , and the proximal operator
Pαg (.) of function g with parameter α is defined as
o
n
1
(41)
ku − xk22 .
Pαg (x) := argmin g(u) +
2α
u

Among the several interpretations available for the proximal
gradient method, we have adopt a quadratic approximationbased model to derive the main iterative scheme in (40).
This interpretation [40], [41], is based on the MajorizationMinimization algorithm (see [42] for a detailed discussion).
In the approximation-based derivation, the ℓ-th iteration for
solving the general problem (39) can be written in the following form:
x(ℓ+1) = argmin fbα (x, x(ℓ) ) + g(x),

(42)

u

where the original objective function f is replaced with
a quadratically-regularized linear approximation around the
previous iterate x(ℓ) , given by
1
kx − yk22 , (43)
fbα (x, y) := f (y) + ∇f (y) (x − y) +
2α
where the quadratic term is referred to as the trust region
penalty. Modulo constants, the objective function in (42) can
be rearranged to get the proximal gradient form
n
o
1
x(ℓ+1) = argmin g(x) +
kx − x(ℓ) + α∇f (x(ℓ) )k22
2α
x
h
i
(ℓ)
(44)
= Pαg x − α∇f (x(ℓ) ) .
′

The proximal operator often admits closed form expressions. As for ℓ1 -regularization, the proximal operator takes
the simple form of the soft thresholding shrinkage operator
Pαk.k1 =: Sα whose ith component is given by
(Sα (x))i := sgn(xi )(xi − α)+ ,

(45)

with sgn denoting the standard signum function, and (a)+ :=
max{a, 0}. In this case, the proximal algorithm leads to a family of algorithms referred to as iterative shrinkage algorithms
[43]–[45], where each iteration involves a simple gradient
descent step followed by a shrinkage operation.
Finally, in our setting, the function f is taken to be
the exponentially weighted log-likelihood Lβ (·). Due to the
smoothness of the logistic function, the Lipschitz constant
for ∇Lβ (ω k ) can be upper bounded by the trace of the
Hessian Bk (ω k ) given in Eq. (21). Noting that the elements
ofP
Λi are
equal to 1/4, we get L(∇Lβ (ω k )) ≤
PWat most
k
1
k−i 2
x(i−1)W +j . Using assumption (1) of Seci=1
j=1 β
4
tion III-B and an application of Hoeffding’s inequality, shows
σ2
that the sum is concentrated around its mean given by MW
4(1−β) ,
for large enough k. Therefore, we choose the step size
(1−β)
α = cMW
σ2 , for some constant c ≥ 1/4.

C OMPUTATION

A PPENDIX C
C ONFIDENCE I NTERVALS

OF

The ℓ1 -regularized ML estimate of Eq. (11) can be written
in the following form
b k = argmax {Pβ L(ω k ) − γkωk k1 } ,
ω

(46)

ωk


where L(ω) := log p n|X, ω denotes the log-likelihood
function over a generic window with spiking vector n,
data matrix X and parameter vector ω, and the operator
Pβ f (n, X, ω) is defined for a function f : {0, 1}W ×
RW ×M × RM → R as the empirical expectation exponentially
weighted with a forgetting factor β:
Pβ f (ω k ) =

k
X
i=1


β k−i f ni , Xi ; ω k ,

(47)

where we have suppressed the dependence of f on n and X
on the left hand side for notational simplicity. Following the
treatment of Theorem 3.1 of [36], the corresponding empirical
gradient vector and Hessian are respectively given by:
gk (ω k ) := Pβ ∇L(ω k ) =

k
X

β k−i X′i εi (ω k ),

(48)

i=1

Bk (ω k ) := Pβ ∇2 L(ω k ) = −

k
X

β k−i X′i Λi (ω k )Xi . (49)

i=1

b k can be then written
The KKT conditions for the estimator ω
as:
gk (b
ω k ) − γbsk = 0, kb
sk k∞ ≤ 1.
(50)
b k k1 is a subgradient vector from the subdifferwhere b
sk ∈ ∂kω
ential of the ℓ1 norm, with components (b
sk )m = sgn ((b
ω k )m )
for (b
ω k )m 6= 0 and |(b
sk )m | ≤ 1 otherwise, for m =
1, 2, . . . , M . Substituting Pβ L(ω k ) by its quadratic approximation around the true parameter vector ω k , and inverting the
corresponding KKT conditions, the ’de-sparsified’ estimator
b k can be obtained as:
w
b k gk (b
b k := ω
bk − Θ
w
ω k ),
(51)
b k is the approximate inverse of Hessian
where the matrix Θ
matrix Bk (b
ω k ), and can be computed using the following
node wise regression procedure [36]. To compute the mth row
b k , first the solution to the following LASSO problem is
of Θ
obtained:



b := argmin − 2(Bk )m,\m ψ + ψ ′ (Bk )\m,\m ψ + 2γm kψk1 ,
ψ
m
ψ∈RM −1

b k is suppressed for notawhere the dependence of Bk on ω
tional convenience, and the subscript notations are the same
as those described in the footnote of Algorithm 3. Then, we
define the vector c ∈ RM as:
b ,
(c)m = 1, (c)\m = −ψ
(52)
m

and the scaling constant

2
τm

as

2
b (L) (Bk )′
τm
:= (Bk )m,m − ψ
m
m,\m .

(53)

b k is given by (Θ
b k )m := 12 c. The
Finally, the mth row of Θ
τm
variance and the confidence interval at a level of α for the

13

b k can then be computed as given in lines
mth component of ω
9 and the output of Algorithm 3 [36], where
k
X
′
2
β 2(k−i) X′i εi (ω)εi (ω)′ Xi .
Gk (ω) := Pβ ∇L(ω)∇L (ω) =
i=1

Using Taylor expansion similar to that in the development
of ℓ1 -PPF1 , the matrix Gk (b
ω k ) can be recursively updated
as given in line 2 of Algorithm 3. Finally, the node wise
regression can be recursively computed using the SPARLS
algorithm [19], which is given in lines 3–5 of Algorithm 3.
The parameter γm can be chosen to be in the same order of
γ in (11).
R EFERENCES
[1] L. M. Frank, G. B. Stanley, and E. N. Brown, “Hippocampal plasticity
across multiple days of exposure to novel environments,” The Journal
of neuroscience, vol. 24, no. 35, pp. 7681–7689, 2004.
[2] D. A. Depireux, J. Z. Simon, D. J. Klein, and S. A. Shamma, “Spectrotemporal response field characterization with dynamic ripples in ferret
primary auditory cortex,” Journal of neurophysiology, vol. 85, no. 3, pp.
1220–1234, 2001.
[3] D. J. Daley and D. Vere-Jones, An introduction to the theory of point
processes: volume II: general theory and structure. Springer Science
& Business Media, 2007, vol. 2.
[4] Y. Ogata, “Statistical models for earthquake occurrences and residual
analysis for point processes,” Journal of the American Statistical Association, vol. 83, no. 401, pp. 9–27, 1988.
[5] D. Vere-Jones, “Stochastic models for earthquake occurrence,” Journal
of the Royal Statistical Society. Series B, pp. 1–62, 1970.
[6] 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, vol. 7, no. 5, pp. 456–461, 2004.
[7] E. N. Brown, D. P. Nguyen, L. M. Frank, M. A. Wilson, and V. Solo,
“An analysis of neural receptive field plasticity by point process adaptive
filtering,” Proceedings of the National Academy of Sciences, vol. 98,
no. 21, pp. 12 261–12 266, 2001.
[8] A. Smith and E. N. Brown, “Estimating a state-space model from point
process observations,” Neural Comp., vol. 15, no. 5, pp. 965–991, 2003.
[9] L. Paninski, “Maximum likelihood estimation of cascade point-process
neural encoding models,” Network: Comp. in Neural Systems, vol. 15,
no. 4, pp. 243–262, 2004.
[10] L. Paninski, J. Pillow, and J. Lewi, “Statistical models for neural
encoding, decoding, and optimal stimulus design,” Progress in brain
research, vol. 165, pp. 493–507, 2007.
[11] J. W. Pillow, Y. Ahmadian, and L. Paninski, “Model-based decoding,
information estimation, and change-point detection techniques for multineuron spike trains,” Neural Comp., vol. 23, no. 1, pp. 1–45, 2011.
[12] W. Truccolo, U. T. Eden, M. R. Fellows, J. P. Donoghue, and E. N.
Brown, “A point process framework for relating neural spiking activity to
spiking history, neural ensemble, and extrinsic covariate effects,” Journal
of neurophysiology, vol. 93, no. 2, pp. 1074–1089, 2005.
[13] S. S. Haykin, Adaptive filter theory. Pearson Education India, 2008.
[14] U. T. Eden, L. M. Frank, R. Barbieri, V. Solo, and E. N. Brown, “Dynamic analysis of neural encoding by point process adaptive filtering,”
Neural comp., vol. 16, no. 5, pp. 971–998, 2004.
[15] D. L. Donoho, “Compressed sensing,” Information Theory, IEEE Trans.
on, vol. 52, no. 4, pp. 1289–1306, 2006.
[16] E. J. Candès et al., “Compressive sampling,” in Proceedings of the
International Congress of Mathematicians, vol. 3. Madrid, Spain, 2006,
pp. 1433–1452.
[17] E. J. Candès and M. B. Wakin, “An introduction to compressive
sampling,” Signal Processing Magazine, IEEE, vol. 25, no. 2, pp. 21–30,
2008.
[18] S. N. Negahban, P. Ravikumar, M. J. Wainwright, and B. Yu, “A
unified framework for high-dimensional analysis of M-estimators with
decomposable regularizers,” Statistical Science, vol. 27, no. 4, pp. 538–
557, 2012.
[19] B. Babadi, N. Kalouptsidis, and V. Tarokh, “SPARLS: The sparse RLS
algorithm,” Signal Processing, IEEE Trans. on, vol. 58, no. 8, pp. 4013–
4025, 2010.
[20] N. Kalouptsidis, G. Mileounis, B. Babadi, and V. Tarokh, “Adaptive
algorithms for sparse system identification,” Signal Processing, vol. 91,
no. 8, pp. 1910–1919, 2011.

[21] B. Dumitrescu, A. Onose, P. Helin, and I. Tăbuş, “Greedy sparse RLS,”
Signal Processing, IEEE Trans. on, vol. 60, no. 5, pp. 2194–2207, 2012.
[22] J. Fritz, S. Shamma, M. Elhilali, and D. Klein, “Rapid task-related
plasticity of spectrotemporal receptive fields in primary auditory cortex,”
Nature neuroscience, vol. 6, no. 11, pp. 1216–1223, 2003.
[23] W. Truccolo, U. T. Eden, M. R. Fellows, J. P. Donoghue, and E. N.
Brown, “A point process framework for relating neural spiking activity to
spiking history, neural ensemble, and extrinsic covariate effects,” Journal
of neurophysiology, vol. 93, no. 2, pp. 1074–1089, 2005.
[24] Z. Chen, D. F. Putrino, S. Ghosh, R. Barbieri, and E. N. Brown,
“Statistical inference for assessing functional connectivity of neuronal
ensembles with sparse spiking data,” Neural Systems and Rehabilitation
Engineering, IEEE Trans. on, vol. 19, no. 2, pp. 121–135, 2011.
[25] D. Needell and J. A. Tropp, “CoSaMP: Iterative signal recovery from
incomplete and inaccurate samples,” Applied and Computational Harmonic Analysis, vol. 26, no. 3, pp. 301–321, 2009.
[26] J. Haupt, W. U. Bajwa, G. Raz, and R. Nowak, “Toeplitz compressed
sensing matrices with applications to sparse channel estimation,” Information Theory, IEEE Trans. on, vol. 56, no. 11, pp. 5862–5875, 2010.
[27] A. Javanmard and A. Montanari, “Confidence intervals and hypothesis testing for high-dimensional regression,” The Journal of Machine
Learning Research, vol. 15, no. 1, pp. 2869–2909, 2014.
[28] S. Van de Geer, P. Bühlmann, Y. Ritov, and R. Dezeure, “On asymptotically optimal confidence regions and tests for high-dimensional models,”
The Annals of Stat., vol. 42, no. 3, pp. 1166–1202, 2014.
[29] C.-H. Zhang and S. S. Zhang, “Confidence intervals for low dimensional
parameters in high dimensional linear models,” Journal of the Royal
Statistical Society: Series B (Statistical Methodology), vol. 76, no. 1,
pp. 217–242, 2014.
[30] N. Meinshausen and P. Bühlmann, “High-dimensional graphs and variable selection with the lasso,” The Annals of Stat., pp. 1436–1462, 2006.
[31] E. N. Brown, R. Barbieri, V. Ventura, R. E. Kass, and L. M. Frank,
“The time-rescaling theorem and its application to neural spike train
data analysis,” Neural comp., vol. 14, no. 2, pp. 325–346, 2002.
[32] R. Haslinger, G. Pipa, and E. Brown, “Discrete time rescaling theorem:
determining goodness of fit for discrete time statistical models of neural
spiking,” Neural comp., vol. 22, no. 10, pp. 2477–2506, 2010.
[33] J. Fritz, M. Elhilali, and S. Shamma, “Active listening: task-dependent
plasticity of spectrotemporal receptive fields in primary auditory cortex,”
Hearing research, vol. 206, no. 1, pp. 159–176, 2005.
[34] N. Mesgarani, J. Fritz, and S. Shamma, “A computational model of rapid
task-related plasticity of auditory cortical receptive fields,” Journal of
computational neuroscience, vol. 28, no. 1, pp. 19–27, 2010.
[35] F. Gandolfo, C.-S. Li, B. Benda, C. P. Schioppa, and E. Bizzi, “Cortical
correlates of learning in monkeys adapting to a new dynamical environment,” Proceedings of the National Academy of Sciences, vol. 97, no. 5,
pp. 2259–2263, 2000.
[36] S. A. van de Geer, “On Hoeffding’s inequality for dependent random variables,” in Empirical Process Techniques for Dependent Data,
H. Dehling and W. Philipp, Eds. Springer, 2001.
[37] P. J. Bickel, Y. Ritov, and A. B. Tsybakov, “Simultaneous analysis of
Lasso and Dantzig selector,” The Annals of Stat., pp. 1705–1732, 2009.
[38] W. Hoeffding, “Probability inequalities for sums of bounded random
variables,” Journal of the American statistical association, vol. 58, no.
301, pp. 13–30, 1963.
[39] B. Scholkopf and A. J. Smola, Learning with kernels: Support vector
machines, regularization, optimization, and beyond. MIT press, 2001.
[40] M. A. Figueiredo and R. D. Nowak, “A bound optimization approach to
wavelet-based image deconvolution,” in Image Processing, 2005. ICIP
2005. IEEE International Conf. on, vol. 2. IEEE, 2005, pp. II–782.
[41] M. A. Figueiredo, J. M. Bioucas-Dias, and R. D. Nowak, “Majorization–
minimization algorithms for wavelet-based image restoration,” Image
Processing, IEEE Trans. on, vol. 16, no. 12, pp. 2980–2991, 2007.
[42] D. R. Hunter and K. Lange, “A tutorial on MM algorithms,” The
American Statistician, vol. 58, no. 1, pp. 30–37, 2004.
[43] M. A. Figueiredo and R. D. Nowak, “An EM algorithm for waveletbased image restoration,” Image Processing, IEEE Trans. on, vol. 12,
no. 8, pp. 906–916, 2003.
[44] A. M. Bruckstein, D. L. Donoho, and M. Elad, “From sparse solutions of
systems of equations to sparse modeling of signals and images,” SIAM
review, vol. 51, no. 1, pp. 34–81, 2009.
[45] I. Daubechies, M. Defrise, and C. De Mol, “An iterative thresholding
algorithm for linear inverse problems with a sparsity constraint,” Communications on Pure and Applied Mathematics, vol. 57, no. 11, pp.
1413–1457, 2004.

