Submitted to Statistical Science

Estimating linear functionals of a
sparse family of Poisson means
Olivier Collier, Arnak S. Dalalyan

arXiv:1712.01775v2 [math.ST] 17 Jan 2018

Modal’X, Université Paris-Nanterre and CREST, ENSAE

Abstract. Assume that we observe a sample of size n composed of pdimensional signals, each signal having independent entries drawn from a
scaled Poisson distribution with an unknown intensity. We are interested
in estimating the sum of the n unknown intensity vectors, under the
assumption that most of them coincide with a given “background” signal.
The number s of p-dimensional signals different from the background
signal plays the role of sparsity and the goal is to leverage this sparsity
assumption in order to improve the quality of estimation as compared
to the naive estimator that computes the sum of the observed signals.
We first introduce the group hard thresholding estimator and analyze
its mean squared error measured by the squared Euclidean norm. We
establish a nonasymptotic upper bound showing that the risk is at
√
most of the order of σ 2 (sp + s2 p log3/2 (np)). We then establish lower
bounds on the minimax risk over a properly defined class of collections
of s-sparse signals. These lower bounds match with the upper bound,
up to logarithmic terms, when the dimension p is fixed or of larger order
than s2 . In the case where the dimension p increases but remains of
smaller order than s2 , our results show a gap between the lower and
√
the upper bounds, which can be up to order p.
MSC 2010 subject classifications: Primary 62J05; secondary 62G05.
Key words and phrases: Nonasymptotic minimax estimation, linear
functional, group-sparsity, thresholding, Poisson processes.

1. INTRODUCTION AND PROBLEM FORMULATION
Let µ = (µ1 , . . . , µp )> be any vector from Rp with positive entries. In what follows, we denote
by Pp (µ) the distribution of a random vector X with independent entries Xj drawn from a
Poisson distribution P(µj ), for j = 1, . . . , p. We consider that the available data is composed
Modal’X, UPL, Univ Paris Nanterre, F92000 Nanterre France. 3 avenue Pierre Larousse,
92245 Malakoff, France.

1

2

O. COLLIER AND A. DALALYAN

of n independent vectors X 1 , . . . , X n randomly drawn from Rp such that
iid

X i ∼ σ 2 Pp (σ −2 µi ),

i = 1, . . . , n,

where σ > 0 is a known parameter that can be interpreted as the noise level. The use of this

term is justified by the fact that when σ goes to zero, we have σ 2 Pp (σ −2 µ) ≈ Np µ, σ 2 diag(µ) .
We will work under the assumption that most intensities are known and equal to a given
vector µ0 , which can be thought of as a background signal. We denote by X = [X 1 , . . . , X n ]
the p × n matrix obtained by concatenating the vectors X i for i = 1, . . . , n.
Assumption (S): for a given vector µ0 ∈ (0, ∞)p , the set S = {i : µi 6= µ0 } is a very small
subset of [n] := {1, . . . , n}.
The unknown parameter in this problem is the p × n matrix M = [µ1 , . . . , µn ] obtained by
concatenating the signal vectors µi . However, we will not aim at estimating the matrix M.
Instead, we aim at estimating a linear functional of the intensities:
X
X
L(M) =
(µi − µ0 ) =
(µi − µ0 ) = M1n − nµ0 .
i∈S

i∈[n]

This functional has a clear meaning: it is the superposition of all the signals contained in X i ’s
beyond the background signal µ0 . Notice that if we knew the set S, it would be possible to
use the oracle
X
bS =
L
(X i − µ0 ),
(1)
i∈S

which would lead to the risk
X

 XX
2
b S − L(M)k2 =
E kL
Var(Xij ) = σ 2
µ>
2
i 1p ≤ µ∞ σ sp,
i∈S j∈[p]

i∈S

where µ∞ is an upper bound on maxi kµi k∞ and s is the cardinality of S. On the other hand,
one can use the naive estimator
X
b [n] =
L
(X i − µ0 ) = L(X)
i∈[n]

which does not use at all the sparsity assumption (S). It has a quadratic risk given by
X

 X X
2
b [n] − L(M)k2 =
E kL
Var(Xij ) = σ 2
µ>
2
i 1p ≤ µ∞ σ np.
i∈[n] j∈[p]

i∈[n]

One can remark that the rate σ 2 sp we have got for the oracle is independent of n, which
obviously can not continue to be true when S is unknown and should somehow be inferred
from the data. An important statistical question is thus the following.
Question (Q1): What kind of quantitative improvement can we get upon the rate σ 2 np of
the naive estimator by leveraging the sparsity assumption?

ESTIMATING LINEAR FUNCTIONALS OF POISSON MEANS

3

In the framework considered in the present work, for answering Question (Q1), we are allowed
to use any estimator of L(M) which requires only the knowledge of σ and µ0 . However, in
practice, it is common to give advantage to estimators having small computational complexity.
Therefore, another question to which we will give a partial answer in this work is:
Question (Q2): What kind of quantitative improvement can we get upon the rate σ 2 np of
the naive estimator by leveraging the sparsity assumption and by constraining ourselves to
estimators that are computable in polynomial time?
In the statement of the last question, the term “polynomial time” refers to the running time
as a function of s, n and p.
1.1 Relation to previous work
Various aspects of the problem of estimation of a linear functional of an unknown highdimensional or even infinite-dimensional parameter were studied in the literature, mostly
focusing on the case of a functional taking real values (as opposed to the vector valued
functional considered in the present work). Early results for smooth functionals were obtained
by Koshevnik and Levit (1977). Minimax estimation of linear functionals over various classes
and models were thoroughly analyzed by Butucea and Comte (2009); Cai and Low (2004, 2005);
Donoho and Liu (1987); Efromovich and Low (1994); Golubev and Levit (2004); Juditsky
and Nemirovski (2009); Klemela and Tsybakov (2001); Laurent et al. (2008). There is also a
vast literature on studying the problem of estimating quadratic functionals (Bickel and Ritov,
1988; Cai and Low, 2006; Donoho and Nussbaum, 1990; Laurent and Massart, 2000). Since
the estimators of (quadratic) functionals can be often used as test statistics, the problem of
estimating functionals has close relations with the problem of testing that were successfully
exploited in (Collier and Dalalyan, 2015; Comminges and Dalalyan, 2012, 2013; Lepski et al.,
1999). The problem of estimation of nonsmooth functionals was also tackled in the literature,
see (Cai and Low, 2011).
Most papers cited above deal with the Gaussian white noise model, Gaussian sequence
model, or the model in which the observations are iid with a common (unknown) density
function. Surprisingly, the results on estimation of functionals of the intensity of a Poisson
process are very scarce. (We are only aware of the paper (Kutoyants and Liese, 1998), which
in an asymptotic set-up established the asymptotic normality and efficiency of a natural
plug-in estimator of a real valued smooth linear functional.) This is surprising, given the
relevance of the Poisson distribution for modeling many real data sets arising, for instance, in
image processing and biology. Furthermore, there are by now several comprehensive studies
of nonparametric estimation and testing for Poisson processes, see Birgé (2007); Ingster
and Kutoyants (2007); Kutoyants (1998); Reynaud-Bouret and Rivoirard (2010) and the
references therein. In addition, several statistical experiments related to Poisson processes
were proved to be asymptotically equivalent to the Gaussian white noise or other important
statistical experiments (Brown et al., 2004; Grama and Nussbaum, 1998). However, the models
related to the Poisson distribution demonstrate some notable differences as compared to the
Gaussian models. From a statistical point of view, the most appealing differences are the
heteroscedasticity of the Poisson distribution1 and the fact that it has heavier tails than the
1

By heteroscedasticity we understand here the fact that if the mean is unknown, then the variance is

4

O. COLLIER AND A. DALALYAN

Gaussian distribution. The impact of these differences on the statistical problems is amplified
in high dimensional settings, as it can be seen in (Ivanoff et al., 2016; Jia et al., 2013).
To complete this section, let us note that the investigation of the statistical problems related
to functionals of high-dimensional parameters under various types of sparsity constraints was
recently carried out in several papers. The case of real valued linear and quadratic functionals
was studied by Collier et al. (2017) and Collier et al. (2016), focusing on the Gaussian sequence
model. Verzelen and Gassiat (2016) analyzed the problem of the signal-to-noise ratio estimation
in the linear regression model under various assumptions on the design. In a companion paper,
Collier and Dalalyan (2017) consider the problem of a vector valued linear functional estimation
in the Gaussian sequence model. The relations with this work will be discussed in more details
in Section 4 below.
1.2 Agenda
The rest of the paper is organized as follows. In Section 2, we introduce an estimator of the
functional L(M) based on thresholding the columns of X with small Euclidean norm. The
same section contains the statement and the proof of a nonasymptotic upper bound on the
expected risk of the aforementioned estimator. Lower bounds on the risk are presented in
Section 3. Some avenues for further research and a summary of the contributions of the present
work are given in Section 4. The proofs of the technical lemmas are gathered in Section 5.
1.3 Notation
We use uppercase boldface letters for matrices and uppercase italic boldface letters for vectors.
For anyP
vector v and for every q ∈ (0, ∞), we use the standard definitions of the `q -norms
q
kvkq = j |vj |q , along with kvk∞ = maxj |vj |. For every n ∈ N, 1n (resp. 0n ) is the vector in
Rn all the entries of which are equal to 1 (resp. 0). For a matrix M, we denote by kMk the
spectral norm of M (equal to its largest singular value), and by kMkF its Frobenius norm
(equal to the Euclidean norm of the vector composed of its singular values). For any matrix
M, we denote by L(M) the linear functional defined as the sum of the columns of M.
2. GROUP HARD THRESHOLDING ESTIMATOR
This section describes the group hard thresholding estimator, that is computationally tractable
and achieves the optimal convergence rate in a large variety of cases. The underlying idea
is rather simple: the Euclidean distance between each signal X i and the background signal
µ0 is compared to a threshold; the sum of all the signals for which the distance exceeds the
threshold is the retained estimator.
Thus, the group hard thresholding method analysed in this section can be seen as a two-step
estimator, that first estimates the set S and then substitutes the estimator of S in the oracle
b S considered in (1). To perform the first step, we choose a threshold λ > 0 and define
L
n
o
b
S(λ)
= i ∈ [n] : kX i − µ0 k2 ≥ σ(kX i k1 + λ)1/2 .
(2)
unknown as well and these two parameters are interrelated

5

ESTIMATING LINEAR FUNCTIONALS OF POISSON MEANS

b b as the group hard thresholding estimator, and
We will refer to the plug-in estimator L
S(λ)
will denote it by
X
b GHT =
L
(X i − µ0 ).
(3)
b
i∈S(λ)

It is clear that the choice of the threshold λ has a strong impact on the statistical accuracy of
b GHT . One can already note that the threshold to which the distance kX i −µ0 k2
the estimator L
is compared depends on the signal X i ; this is due to the heteroscedasticity of the models
based on observations drawn from the Poisson distribution.
Prior to describing in a quantitative way the impact of the threshold on the expected risk, let
b GHT requires O(np) operations; so it can be done in
us emphasize that the computation of L
polynomial time. Indeed, the computation of each of the n norms kX i − µ0 k2 requires O(p)
b
operations, so S(λ)
can be computed using O(np) operations. The same is true for the sum in
the right hand side of (3).
The following theorem establishes the performance of the previously defined thresholding
estimator, for a suitably chosen threshold λ.
−3/2

Theorem 1. Define µ∞ = maxi,j |Mij | and suppose σ 3 µ∞
The group hard thresholding estimator (2)-(3) with

3/2

≤ 40 log3/2 (2np) ≤ 0.9σ −3 µ∞ .

√
λ = 40kµ0 k∞ p log3/2 (2np),

(4)

has a mean squared error bounded as follows:


 GHT

√
b
E kL
− L(M)k22 ≤ µ∞ σ 2 170s2 p log3/2 (2np) + 6sp .
Proof. Let us introduce the auxiliary notation ξ i = σ −1 (X i −µi ) and write the decomposition
X
X
X
b GHT − L(M) =
L
(X i − µi ) −
(X i − µ0 ) +
(X i − µ0 )
i∈S

X

=

σξ i −

b
i∈S∩S(λ)

b
i∈S\S(λ)

b
i∈S(λ)\S

X

X

(µi − µ0 ) +

b
i∈S\S(λ)

σξ i .

b
i∈S(λ)\S

We have
b GHT − L(M)
L

2

X

≤σ

ξi

2

X

+

b
i∈S∩S(λ)

|

{z
T1

(µi − µ0 )

2

X

+σ

b
i∈S\S(λ)

}

|

ξi

2

.

(5)

b
i∈S(λ)\S

{z
T2

}

|

{z
T3

}

Let us use the notation ΞS = [ξ i : i ∈ S] for the matrix obtained by concatenating the vectors
ξ i with i ∈ S, and denote by ξ jS its j-th row. We can write
2
2
T12 = kΞS 1S(λ)
b k2 ≤ skΞS k = s

X
j∈[p]

(ξ jS )> ξ jS ,

6

O. COLLIER AND A. DALALYAN

b
where 1S(λ)
is the vector in {0, 1}S having zeros at positions i ∈ S(λ)
and zeros elsewhere.
b
One can observe that the s × s matrices (ξ jS )> ξ jS are independent and each of them is in
expectation equal to the diagonal matrix diag(µjS ). Since the operator norm of a diagonal
matrix is equal to its largest (in absolute value) entry, we have
 X

 2
 j > j
j
E T1 ≤ µ∞ sp + sE
(ξ S ) ξ S − diag(µS )
j∈[p]

  X
 j > j
≤ µ∞ sp + s E
(ξ S ) ξ S − diag(µjS )

2

1/2
.

j∈[p]

Upper bounding the spectral norm by the Frobenius norm, we arrive at

 X
 X
2
 j > j
 j > j
j
E
(ξ S ) ξ S − diag(µS )
≤E
(ξ S ) ξ S − diag(µjS )
j∈[p]

j∈[p]

=

X


E (ξ jS )> ξ jS − diag(µjS )

2



F

2
.
F

j∈[p]

Using the fact that the fourth order central moment of the Poisson distribution P(a) is a + 3a2
and that ξij ∼ σ(P(σ −2 µ0j ) − σ −2 µ0j ), one can check that
4
E[ξij
] = σ 2 µij + 3µ2ij ≤ (2µ∞ )2 ,

(6)

where the last inequality holds true since σ 2 ≤ (0.9)1/3 µ∞ .
In view of (6) and the Cauchy-Schwarz inequality, it holds that E[k(ξ jS )> ξ jS − diag(µjS )k2F ] ≤
2.5sµ2∞ + (s2 − s)µ2∞ ≤ (s + 1)2 µ2∞ . Putting these estimates together, we get
 
√
√
E T12 ≤ µ∞ sp + µ∞ s(s + 1) p ≤ µ∞ sp + 2µ∞ s2 p.
(7)
We would like to stress here that this upper bound is quite rough, because of the step in which
we upper bounded the spectral norm by the Frobenius one. Using more elaborated arguments
on the expected spectral norm of a random matrix, such as the Rudelson inequality Rudelson
(1999) and (Boucheron et al., 2013, Cor. 13.9) combined with the symmetrisation argument
√
√
√
(Koltchinskii, 2011, Theorem 2.1), it is possible to replace the rate s2 p by s( s + p).
√
However, we do not give more details on this computations since the rate s2 p appears
anyway as an upper bound on E[T2 ].
To upper bound the term T2 in (5), we define τi = σ 2 (kX i k1 + λ) and u = 40 log3/2 (2np), so
√
b
that τi = σ 2 (kX i k1 + kµ0 k∞ p u). Note that i ∈ S \ S(λ)
implies
kµi − µ0 k22 < τi − σ 2 kξ i k2 − 2σ(µi − µ0 )> ξ i
< τi − σ 2 kξ i k2 + 2σ|(µi − µ0 )> ξ i |
< τi − σ 2 kξ i k2 + 0.5kµi − µ0 k22 + 2σ 2

|(µi − µ0 )> ξ i |2
.
kµi − µ0 k22

ESTIMATING LINEAR FUNCTIONALS OF POISSON MEANS

7

In other terms, setting θ i = (µi − µ0 )/kµi − µ0 k2 , we have


2
2
2
2
>
2
b
S \ S(λ)
⊂ kµi − µ0 k2 < 2τi − 2σ kξ i k + 4σ |θ i ξ i | .
As a consequence,
T22 ≤ s

X

kµi − µ0 k22

b
i∈S\S(λ)

X
X
√
2
≤ 2µ∞ s2 σ 2 p u + 2σ 2 s
(kX i k1 − kξ i k22 )+ + 4σ 2 s
(θ >
i ξi ) .
i∈S

i∈S

Taking the expectation of both sides, in view of the independence of the components of ξ i and
those of X i and the condition σ 2 ≤ µ∞ , we arrive at
p
 
√
E T22 ≤ 2µ∞ σ 2 s2 p u + 2µ∞ σ 2 s2 6p + 4µ∞ σ 2 s2
√
≤ 89µ∞ σ 2 s2 p log3/2 (2np).
(8)
The third term in (5) corresponds to the second kind error in the problem of support estimation
and its expectation can be upper bounded using the Cauchy-Schwarz inequality as follows
X
E[T32 ] ≤ n
E[kξ i k22 1kξi k22 ≥τi /σ2 ]
i6∈S

=n

XX

≤n

XX

2
E[ξij
1kξi k22 ≥τi /σ2 ]

i6∈S j∈[p]
4
E[ξij
]P(kξ i k22 ≥ τi /σ 2 )

1/2

.

i6∈S j∈[p]

In view of (6), this yields
E[T32 ] ≤ 2µ∞ np

X

P(kξ i k22 ≥ τi /σ 2 )1/2 .

(9)

i6∈S

For the last probability, we have
√
P(kξ i k22 ≥ τi /σ 2 ) = P(kξ i k22 − kX i k1 ≥ kµ0 k∞ p u).
To upper bound this probability, we apply Lemma 1 below2 to η = σ −2 X j and ν = σ −2 µj .
Lemma 1. Let η1 , . . . , ηp be independent random variables drawn from the Poisson distributions P(ν1 ), . . . , P(νp ), respectively. We set η = (η1 , . . . , ηp ) and ν = (ν1 , . . . , νp ). Then, for
−3/2
3/2
every u satisfying ν∞ ≤ u ≤ 0.9ν∞ , we have

√ 
2/3
P kη − νk22 − kηk1 ≥ ν∞ p u ≤ (2p + 1)e−cu .
where c ≥ 0.38 is a universal constant.
2

The proof of Lemma 1 is postponed to Section 5

8

O. COLLIER AND A. DALALYAN
−3/2

This yields that for σ 3 µ∞

3/2

≤ u ≤ 0.9σ −3 µ∞ , we have
2/3

P(kξ i k22 ≥ τi /σ 2 ) ≤ (2p + 1)e−cu

,

where c ≥ 0.38 is an absolute constant. In view of (9), this implies the following upper bound
on the expectation of T3 :
2/3

E[T32 ] ≤ 2µ∞ n2 p(2p + 1)1/2 e−0.5cu
2/3

≤ 4µ∞ n2 p2 e−0.5cu

.

The choice u = 40 log3/2 (2np) ensures that
E[T32 ] ≤ 4µ∞ n2 p2 e−2 log(2np) ≤ µ∞ ,

(10)

which is negligible with respect to the upper bound on E[T12 ] obtained above.
The result follows from upper bounds (5), (7), (8) and (10), as well as the choice of u. Indeed,
by the Minkowskii inequality, we get

1/2
1/4
b GHT − L(M)k2 ]
E[kL
≤ 11.85σµ1/2
log3/4 (2np) + σ(µ∞ sp)1/2 .
2
∞ sp
Taking the square of both sides and using the inequality (a + b)2 ≤ (6/5)a2 + 6b2 , we get the
claim of the theorem.
3. LOWER BOUNDS ON THE MINIMAX RISK OVER THE SPARSITY SET
Let us denote by M(n, p, s, µ0 , µ∞ ) the set of all p × n matrices with positive entries bounded
by µ∞ and having at most s columns different from µ0 . This implicitly means, of course, that
µ∞ ≥ kµ0 k∞ , otherwise the set is empty. The parameter set M(n, p, s, µ0 , µ∞ ) is the sparsity
class over which the minimax risk will be considered. The result of the previous section implies
√
that the minimax risk over M(n, p, s, µ0 , µ∞ ) is at most of order µ∞ σ 2 s(p + s p), up to a
logarithmic factor. In the present section, we will establish two lower bounds on the minimax
risk over M(n, p, s, µ0 , µ∞ ).
The first lower bound shows that the minimax risk is at least of order s2 , up to log factor. This
implies that the group hard thresholding estimator defined in previous section is rate-optimal
when p is fixed and s goes to infinity. The proof of this lower bound relies on the comparison
of the Kullback-Leibler divergence between the probability distribution of X corresponding
to M0 = [µ0 , . . . , µ0 ], and the one corresponding to a matrix M, the n columns of which are
randomly drawn from the set {µ0 , µ1 } with proportion (n − s) : s. When µ1 is well chosen,
the two aforementioned distributions are close, whereas the corresponding linear functional
L(M0 ) is sufficiently away from L(M).
Theorem 1.

Assume that s ≥ 128, then
inf

sup

b M∈M(n,p,s,µ0 ,µ∞ )
L

b − L(M)
EM L

2

≥


kµ0 k∞ 2 2
n 
σ s log 1 + 2 ,
513
2s

where the infimum is taken over all estimators of the linear functional.

ESTIMATING LINEAR FUNCTIONALS OF POISSON MEANS

9

Proof. We will make use of some lemmas, the statement of which requires additional notation.
As usual in lower bounds, we will make use of the Bayes risk: for every measure π on the set
of p × n real matrices, we define the probability measure
Z
Pπ ( · ) =
PM ( · ) π(dM).
Rp×n

Let F : Rp×n
→ Rm be any functional, not necessarily linear, defined on the set of matrices.
+
Lemma 2. Assume that π0 and π1 are two probability measures on the set of p × n matrices,
1
Rp×n
+ , and set M = maxi Varπi (F ). Furthermore, suppose that KL(Pπ1 , Pπ0 ) ≤ /8 and, for
some v > 0,
√ 
Eπ0 (F ) − Eπ1 (F ) 2 > 2 v + 2 M .
b of F (M), we have
Then, for any estimator F
b − F (M)
sup EM F
M∈M

2
2

≥

1/8


− max πi (M{ ) v 2 .
i

Lemma 3. Let µ0 > 0 and  > 0 be two real numbers. Define π0 = δµ⊗n
as the Dirac mass
0
concentrated in (µ0 , . . . , µ0 ) ∈ Rn+ and π1 = λ⊗n as the product measure of the mixture

 s 
s 
δµ0 +
δ
λ= 1−
2 .
2n
2n µ0 +σ ε
Then
KL(Pπ1 kPπ0 ) ≤


s2 σ2 ε2 /µ0
e
−1 .
4n

We are now in a position to establish the lower bound of Theorem 1. For simplicity and without
loss of generality, let us assume that the first entry of µ0 , µ01 , is the largest one. We choose two
probability measures π0p and π1p on the set of p × n (intensity) matrices with positive entries,
p
Rp×n
+ , in the following way. The measure π0 is simply the Dirac mass at M0 = (µ0 , . . . , µ0 ).
As for π1p , its marginal distribution of the first row equals the measure π1 defined in Lemma 3,
whereas all the other rows are deterministically equal to the corresponding rows of the matrix
M0 . So, for two matrices M and M0 randomly drawn from π0p and π1p , all the rows except the
first are equal (and equal to the corresponding row of M0 ). This implies that choosing the
real ε by the formula
ε2 = (µ01 /σ 2 ) log(1 + n/2s2 ),
we get (in view of Lemma 3) KL(Pπ1p kPπ0p ) ≤ 1/8.
Furthermore, we have
Eπ0 [L(M)] = 0p ,

Eπ1 [L(M)] =

and
Varπ0 [L(M)] = 0,

Varπ1 [L(M)] =



1 2
ε
sσ
,
0p−1
2
1
s  4 2
1−
sσ ε .
2
2n

10

O. COLLIER AND A. DALALYAN

This means that, after upper bounding s by s2 /128, the constant M of Lemma 2 can be chosen
as M = 2−8 s2 σ 4 ε2 . Hence, for v = (1/8)σ 2 sε, we have
Eπ0 [L(M)] − Eπ1 [L(M)]

2

=

√
sσ 2 ε
= 2v + 4 M .
2

This implies that the two conditions of Lemma 2 are fulfilled with
√
σs µ01
σ 2 sε
=
log1/2 (1 + n/2s2 ).
v=
8
8
Therefore, taking into account the fact that π0p (M(n, p, s, µ0 , µ∞ )) = 1, we have
sup

b − L(M)
EM L

M∈M(n,p,s,µ0 ,µ∞ )

2
2

≥

1/8


− π1p (M(n, p, s, µ0 , µ∞ ){ ) v 2 .

To evaluate the probability π1p (M(n, p, s, µ0 , µ∞ ){ ), we remark that it is equal to P(ζ > s),
where ζ is a random variable drawn from the binomial distribution B(n, (s/2n)). Using the
Bernstein inequality (Shorack and Wellner, 1986, p. 440), one can check that P(ζ > s) ≤
e−3s/5 ≤ e−76 , provided that s ≥ 128. This completes the proof.
It follows from Theorem 1 that the factor s2 present in the upper bound of the minimax risk is
unavoidable. The next result establishes another lower bound showing the the term sp present
in the upper bound is unavoidable as well.
Theorem 2.

For all 1 ≤ s ≤ n, if p ≥ 16, then
inf

sup

b − L(M)
EM L

2

≥ 2−14 µ∞ σ 2 sp.

b M∈M(n,p,s,µ ,µ∞ )
L
0

Proof. We set ε2 = 2−7 (σ2/sµ∞ ) and, for any T ⊂ [p], define the matrix MT by


if (i, j) ∈ T × [s]
µ∞ (1 − ε),

MT i,j = µ∞ ,
if (i, j) ∈ T { × [s],


µ0,i ,
j 6∈ [s].
We also set M0 = M∅ . It is clear that all these matrices MT belong to the sparsity class
M(n, p, s, µ0 , µ∞ ). In addition, for any T , we have
i spε2 µ
 sµ∞ X h
∞
KL PMT kPM0 = 2
(1 − ε) log(1 − ε) + ε ≤
.
σ
σ2
j∈T

According to the Varshamov-Gilbert bound (Tsybakov, 2009, Lemma 2.9), there exist m =
[ep/8 ] + 1 sets T0 , . . . , Tm ⊂ [p] such that T0 = ∅ and3
i 6= j
3

=⇒

p
|Ti 4 Tj | ≥ .
8

In what follows, 4 stands for the symmetric difference of two sets.

ESTIMATING LINEAR FUNCTIONALS OF POISSON MEANS

11

Fig 1. The rate optimality diagram in the minimax sense for the group hard thresholding estimator.

With this choice of sets, we have for all i 6= j,


L MTi − L MTj

2

≥ s2 |Ti 4 Tj |ε2 µ2∞ ≥

s2 pε2 µ2∞
.
23

Finally, using Theorem 2.5 in (Tsybakov, 2009) with α = 1/16, we get
inf

b − L(M)
EM L

sup

b M∈M(n,p,s,µ0 ,µ∞ )
L

2


1
1  s2 pε2 µ2∞
≥ (1 − e−p/16 ) × 1 − − √
8
p
25
≥ 2−14 σ 2 spµ∞ .

and the result follows.

4. CONCLUDING REMARKS AND OUTLOOK
Combining the results of the previous sections one can draw the optimality diagram, see Fig. 1.
In particular, Theorem 1 and Theorem 1 imply that the group hard thresholding estimator
(2), (3) with the adaptively chosen threshold (4) is rate optimal when the term sp dominates
√
s2 p log3/2 (2np). This corresponds to the “dimension dominates squared sparsity” regime
p ≥ s2 log3 (2np). According to Theorem 1 and Theorem 2, the estimator defined by (2), (3),
(4) is rate optimal in the low dimensional setting p = O(1) as well. The range of values of p in
which the determination of the minimax rate remains open corresponds to the regime p → ∞
so that p/s2 = O(1).
In a companion paper (Collier and Dalalyan, 2017), the problem of estimating the linear
functional L(M) is considered, when the observations are Gaussian vectors with unknown
means µi . It turns out that in the Gaussian setting the optimal scaling of the minimax
risk is ps + s2 , up to a logarithmic factor. However, this rate is achieved by an estimator
that has a computational complexity that scales exponentially with n. Furthermore, in the
Gaussian model, the upper bound relies heavily on the fact that the Gaussian distribution
is sub-Gaussian. The fact that the Poisson distribution is not sub-Gaussian (even though it

12

O. COLLIER AND A. DALALYAN

is sub-exponential), does not allow us to apply the methodology developed in the Gaussian
model.
It is worth mentioning here that although the risk bounds (both lower and upper) depend
on the sparsity level s and the largest intensity µ∞ , the group hard thresholding estimator is
completely adaptive with respect to these parameters.
The problem we have considered in this work assumes that the noise variance σ 2 and the
background signal µ0 are known. It is not very difficult to check that the result of Theorem 1
can be extended to the case where µ0 is unknown, but an additional data set X̄ 1 , . . . , X̄ m is
available, such that X̄ i ’s are iid vectors drawn from the distribution σ 2 Pp (µ0 /σ 2 ). One can
then estimate the background signal µ0 by the empirical mean of X̄ i ’s. As for σ, one can
easily estimate this parameter using the fact that the entries of the data matrix X take their
values in the set σ 2 N. Next, these estimates can be plugged in the group hard thresholding
estimator. If m is large enough, the error of estimating µ0 and σ 2 is negligible with respect to
the error of estimating L(M). Determining more precise risk bounds clarifying the impact of
m on the risk is an interesting avenue for future research.
5. PROOFS OF THE LEMMAS
Proof of Lemma 1. Let us introduce an auxiliary real number γ > 1.52 and denote
Zj = (ηj − νj − 0.5)2 ,

Zj,γ = Zj ∧ γ,

Z̄j,γ = (Zj − γ)+

as well as
νj,γ = E[Zj,γ ],

ν̄j,γ = E[Z̄j,γ ].

It is clear that νj,γ + ν̄j,γ = νj + 0.25 and we have
X
X
kη − νk22 − kηk1 =
(Zj,γ − νj,γ ) +
(Z̄j,γ − ν̄j,γ )
j∈[p]

j∈[p]

= kZ γ k1 − kν γ k1 + kZ̄ γ k1 − kν̄ γ k1 .

(11)

Since the random variables Zj,γ are independent and bounded, with P(Zj,γ ∈ [0, γ]) = 1, the
Hoeffding inequality yields


2
2
P kZ γ k1 − kν γ k1 > t ≤ e−2t /(pγ ) .
(12)
On the other hand, it is clear that
X


P
(Z̄j,γ − ν̄j,γ ) > 0 ≤ P maxj∈[p] Zj > γ
j∈[p]

√ 
γ
j∈[p]

√
≤ p max P |ηj − νj | > γ − 0.5
j∈[p]
√ 
≤ p max P |ηj − νj | > (2/3) γ .
≤ p max P |ηj − νj − 0.5| >

j∈[p]

13

ESTIMATING LINEAR FUNCTIONALS OF POISSON MEANS

Using the Bennet inequality for the Poisson random variables—see, for instance, (Pollard,
2015, Problem 15)— we get4
√ 
√ 
√ 
P |ηj − νj | > (2/3) γ ≤ P ηj > νj + (2/3) γ + P ηj < νj − (2/3) γ
≤ e−1.54γ/(9νj ) + e−2γ/(9νj ) ≤ 2e−γ/(6νj ) ,
Combining (11), (12) and (13), we get


P kη − νk22 − kηk1 ≥ t ≤

inf

2 ]
γ∈[(3/2)2 ,(9/4)ν∞

2
∀γ ≤ (9/4)ν∞
.

(13)

n
o
2
2
e−2t /(pγ ) + 2pe−γ/(6ν∞ ) .

√
3 . For u such that ν −3/2 ≤ u ≤ 0.9ν 3/2 , the
Let us set t = ν∞ p u and γ 3 = 12 u2 ν∞
∞
∞
aforementioned value of γ belongs to the range over which the infimum of the last display is
computed. This yields

√ 
1/3 2/3
P kη − νk22 − kηk1 ≥ ν∞ p u ≤ (2p + 1)e−(12) u /6 ,
and the claim of the lemma follows.
Proof of Lemma 2. Using Markov’s inequality, we get
b − F (M)
sup v −2 EM F
M∈M

2
2

b − F (M)k2 ≥ v
≥ sup PM kF



M∈M


b − F (M)k2 ≥ v − max πi (M{ ).
≥ max Pπi kF
i=0,1

i=0,1

b − Eπ (F ) and the notation M =
Now, using the test statistic ψ = arg mini∈{0,1} F
i
2
maxi Varπi (F ) = maxi Eπi kF − Eπi [F ]k22 , we have


b − Eπ (F ) ≤ F
b − Eπ (F )
Pπi (ψ 6= i) = Pπi F
1−i
i
2
2


p
b − Eπ (F ) ≥ v + M/α .
≤ Pπi F
i
2
Finally, using the triangle inequality again,



b − F (M)k2 ≥ v + Pπ F (M) − Eπ (F )
Pπi (ψ 6= i) ≤ Pπi kF
i
i


b − F (M)k2 ≥ v ,
≤ α + Pπi kF

√
2

>

α−1 M



so that

b − F (M)k2 ≥ v ≥ max Pπ (ψ 6= i) − α − max πi (M{ ).
sup PM kF
i

M∈M

i=0,1

i=0,1

The result is now a direct application of Theorem 2.2 in Tsybakov (2009) using the assumption
on the Kullback-Leibler divergence.
4

We use here the fact that the Bennet ψ ∗ satisfies ψ ∗ (x) ≥ 0.772 when x ∈ [0, 1].

14

O. COLLIER AND A. DALALYAN

Proof of Lemma 3. Let usR denote by Qµ , for every vector µ ∈ (0, ∞), the Poisson distribution P(σ −2 µ), and by Qλ = (0,∞) Qµ λ(dµ), the λ-mixture of Poisson distributions. We have,
for every X ∈ N,

s  d Qµ0 +σ2 
d Qλ
(X) = 1 +
(X) − 1
d Qµ0
2n
d Qµ0

n
s
σ 2 ε oX −ε
e −1 .
=1+
1+
2n
µ0
Using the inequality log(1 + t) ≤ t, we get
Z n
s
s
σ 2 ε oX −ε
KL(Qλ kQµ0 ) ≤
e Qλ (dX) −
1+
.
2n N
µ0
2n
On the other hand, for every µ ∈ (0, ∞)p , we have
Z n
Xn

σ 2 ε ok µk −µ/σ2
σ 2 ε oX
Qµ (dX) =
1+
1+
e
= exp µε/µ0 ,
2k
µ0
µ0
σ k!
k∈N

so that
s   s2  σ2 2 /µ0
s
s 
1−
+
e
−
2
2n
2n
4n
2n

s2 σ2 ε2 /µ0
= 2 e
−1 .
4n

KL(Qλ kQµ0 ) ≤

The result follows since, by independence, KL(Pπ1 kPπ0 ) = nKL(Qλ kQµ0 ).
ACKNOWLEDGMENTS
The authors thank an anonymous referee for his/her careful reading of the paper and fruitful
suggestions. O. Collier’s research has been conducted as part of the project Labex MMEDII (ANR11-LBX-0023-01). The work of A. Dalalyan was partially supported by the grant
Investissements d’Avenir (ANR-11IDEX-0003/Labex Ecodec/ANR-11-LABX-0047).
REFERENCES
Bickel, P. J. and Ritov, Y. (1988). Estimating integrated squared density derivatives: Sharp best order of
convergence estimates. Sankhy: The Indian Journal of Statistics, Series A (1961-2002), 50(3):381–393.
Birgé, L. (2007). Model selection for Poisson processes, volume Volume 55 of Lecture Notes–Monograph Series,
pages 32–64. Institute of Mathematical Statistics, Beachwood, Ohio, USA.
Boucheron, S., Lugosi, G., and Massart, P. (2013). Concentration Inequalities: A Nonasymptotic Theory of
Independence. OUP Oxford.
Brown, L. D., Carter, A. V., Low, M. G., and Zhang, C.-H. (2004). Equivalence theory for density estimation,
poisson processes and gaussian white noise with drift. Ann. Statist., 32(5):2074–2097.
Butucea, C. and Comte, F. (2009). Adaptive estimation of linear functionals in the convolution model and
applications. Bernoulli, 15(1):69–98.
Cai, T. T. and Low, M. G. (2004). Minimax estimation of linear functionals over nonconvex parameter spaces.
Ann. Statist., 32(2):552–576.
Cai, T. T. and Low, M. G. (2005). On adaptive estimation of linear functionals. Ann. Statist., 33(5):2311–2343.

ESTIMATING LINEAR FUNCTIONALS OF POISSON MEANS

15

Cai, T. T. and Low, M. G. (2006). Optimal adaptive estimation of a quadratic functional. Ann. Statist.,
34(5):2298–2325.
Cai, T. T. and Low, M. G. (2011). Testing composite hypotheses, hermite polynomials and optimal estimation
of a nonsmooth functional. Ann. Statist., 39(2):1012–1041.
Collier, O., Comminges, L., and Tsybakov, A. B. (2017). Minimax estimation of linear and quadratic functionals
on sparsity classes. Ann. Statist., 45(3):923–958.
Collier, O., Comminges, L., Tsybakov, A. B., and Verzélen, N. (2016). Optimal adaptive estimation of linear
functionals under sparsity. ArXiv e-prints, ArXiv:1611.09744.
Collier, O. and Dalalyan, A. S. (2015). Curve registration by nonparametric goodness-of-fit testing. J. Statist.
Plann. Inference, 162:20–42.
Collier, O. and Dalalyan, A. S. (2017). Rate-optimal estimation of p-dimensional linear functionals in a sparse
Gaussian model. ArXiv e-prints.
Comminges, L. and Dalalyan, A. S. (2012). Tight conditions for consistency of variable selection in the context
of high dimensionality. Ann. Statist., 40(5):2667–2696.
Comminges, L. and Dalalyan, A. S. (2013). Minimax testing of a composite null hypothesis defined via a
quadratic functional in the model of regression. Electron. J. Statist., 7:146–190.
Donoho, D. and Liu, R. (1987). On Minimax Estimation of Linear Functionals. Technical report (University of
California, Berkeley. Department of Statistics). Department of Statistics, University of California.
Donoho, D. L. and Nussbaum, M. (1990). Minimax quadratic estimation of a quadratic functional. Journal of
Complexity, 6(3):290 – 323.
Efromovich, S. and Low, M. G. (1994). Adaptive estimates of linear functionals. Probability Theory and Related
Fields, 98(2):261–275.
Golubev, Y. and Levit, B. (2004). An oracle approach to adaptive estimation of linear functionals in a gaussian
model. Mathematical Methods of Statistics, 13(01):392–408.
Grama, I. and Nussbaum, M. (1998). Asymptotic equivalence for nonparametric generalized linear models.
Probability Theory and Related Fields, 111(2):167–214.
Ingster, Y. I. and Kutoyants, Y. A. (2007). Nonparametric hypothesis testing for intensity of the poisson
process. Mathematical Methods of Statistics, 16(3):217–245.
Ivanoff, S., Picard, F., and Rivoirard, V. (2016). Adaptive lasso and group-lasso for functional poisson regression.
Journal of Machine Learning Research, 17(55):1–46.
Jia, J., Rohe, K., and Yu, B. (2013). The lasso under poisson-like heteroscedasticity. Statistica Sinica,
23(1):99–118.
Juditsky, A. B. and Nemirovski, A. S. (2009). Nonparametric estimation by convex programming. Ann. Statist.,
37(5A):2278–2300.
Klemela, J. and Tsybakov, A. B. (2001). Sharp adaptive estimation of linear functionals. Ann. Statist.,
29(6):1567–1600.
Koltchinskii, V. (2011). Oracle inequalities in empirical risk minimization and sparse recovery problems, volume
2033 of Lecture Notes in Mathematics. Springer, Heidelberg. Lectures from the 38th Probability Summer
School held in Saint-Flour, 2008,.
Koshevnik, Y. A. and Levit, B. Y. (1977). On a non-parametric analogue of the information matrix. Theory of
Probability & Its Applications, 21(4):738–753.
Kutoyants, Y. (1998). Statistical Inference for Spatial Poisson Processes, volume 134 of Lecture Notes in
Statistics. Springer New York.
Kutoyants, Y. and Liese, F. (1998). Estimation of linear functionals of poisson processes. Statistics & Probability
Letters, 40(1):43 – 55.
Laurent, B., Ludena, C., and Prieur, C. (2008). Adaptive estimation of linear functionals by model selection.
Electron. J. Statist., 2:993–1020.
Laurent, B. and Massart, P. (2000). Adaptive estimation of a quadratic functional by model selection. Ann.
Statist., 28(5):1302–1338.
Lepski, O., Nemirovski, A., and Spokoiny, V. (1999). On estimation of the Lr norm of a regression function.
Probability Theory and Related Fields, 113(2):221–253.
Pollard, D. (2015). A few good inequalities. Technical report, arXiv:1306.0187.
Reynaud-Bouret, P. and Rivoirard, V. (2010). Near optimal thresholding estimation of a poisson intensity on
the real line. Electron. J. Statist., 4:172–238.
Rudelson, M. (1999). Random vectors in the isotropic position. Journal of Functional Analysis, 164(1):60 – 72.
Shorack, G. R. and Wellner, J. A. (1986). Empirical Processes with Applications to Statistics. Wiley, New York.
Tsybakov, A. B. (2009). Introduction to Nonparametric Estimation. Springer Publishing Company, Incorporated,

16

O. COLLIER AND A. DALALYAN

1st edition.
Verzelen, N. and Gassiat, E. (2016). Adaptive estimation of High-Dimensional Signal-to-Noise Ratios. ArXiv
e-prints.

