Emilien Joly
email: emilien.joly@u-paris10.fr

Gábor Lugosi
email: gabor.lugosi@upf.edu

Roberto Imbuzeiro Oliveira

On the estimation of the mean of a random vector

We study the problem of estimating the mean of a multivariate distribution based on independent samples. The main result is the proof of existence of an estimator with a non-asymptotic sub-Gaussian performance for all distributions satisfying some mild moment assumptions.

Introduction

Let X be a random vector taking values in R d . We assume throughout the paper that the mean vector µ = EX and covariance matrix Σ = (Xµ)(Xµ) T exist. Given n independent, identically distributed samples X 1 , . . . , X n drawn from the distribution of X, one wishes to estimate the mean vector.

A natural and popular choice is the sample mean (1/n) n i=1 X i that is known to have a near-optimal behavior whenever the distribution is sufficiently light tailed. However, whenever heavy tails are a concern, the sample mean is to be avoided as it may have a suboptimal performance. While the one-dimensional case (i.e., d = 1) is quite well understood (see [START_REF] Catoni | Challenging the empirical mean and empirical variance: a deviation study[END_REF], [START_REF] Devroye | Sub-Gausssian mean estimators[END_REF]), various aspects of the multidimensional problem are still to be revealed. This paper aims at contributing to the understanding of the multi-dimensional case.

Before stating the main results, we briefly survey properties of some mean estimators of real-valued random variables. Some of these techniques serve as basic building blocks for the estimators we propose for the vector-valued case.

Estimating the mean of a real-valued random variable

When d = 1, the simplest and most popular mean estimator is the sample mean µ n = (1/n) n i=1 X i . The sample mean is unbiased and the central limit theorem guarantees an asymptotically Gaussian distribution. However, unless the distribution of X has a light (e.g., sub-Gaussian) tail, there are no non-asymptotic sub-Gaussian performance guarantees for µ n . We refer the reader to Catoni [START_REF] Catoni | Challenging the empirical mean and empirical variance: a deviation study[END_REF] for details. However, perhaps surprisingly, there exist estimators of µ with much better concentration properties, see Catoni [START_REF] Catoni | Challenging the empirical mean and empirical variance: a deviation study[END_REF] and Devroye, Lerasle, Lugosi, and Oliveira [START_REF] Devroye | Sub-Gausssian mean estimators[END_REF].

A conceptually simple and quite powerful estimator is the so-called median-of-means estimator that has been proposed, in different forms, in various papers, see Nemirovsky and Yudin [START_REF] Nemirovsky | Problem complexity and method efficiency in optimization[END_REF], Hsu [START_REF] Hsu | Robust statistics[END_REF], Jerrum, Valiant, and Vazirani [START_REF] Jerrum | Random generation of combinatorial structures from a uniform distribution[END_REF], Alon, Matias, and Szegedy [START_REF] Alon | The space complexity of approximating the frequency moments[END_REF]. The median-of-means estimator is defined as follows. Given a positive integer b and x 1 , . . . , x b ∈ R, let q 1/2 denote the median of these numbers, that is,

q 1/2 (x 1 , . . . , x b ) = x i , where #{k ∈ [b] : x k ≤ x i } ≥ b 2 and #{k ∈ [b] : x k ≥ x i } ≥ b 2 .
(If several i fit the above description, we take the smallest one.) For any fixed δ ∈ [e 1-n/2 , 1), first choose b = ⌈ln(1/δ)⌉ and note that b ≤ n/2 holds. Next, partition [n] = {1, . . . , n} into b blocks B 1 , . . . , B b , each of size |B i | ≥ ⌊n/b⌋ ≥ 2. Given X 1 , . . . , X n , we compute the sample mean in each block

Y i = 1 |B i | j∈B i X j
and define the median-of-means estimator by µ

(δ) n = q 1/2 (Y 1 , . . . , Y B ).
One can show (see, e.g., Hsu [START_REF] Hsu | Robust statistics[END_REF]) that for any n ≥ 4,

P | µ (δ) n -µ| > 2e 2Var(X) (1 + ln(1/δ)) n ≤ δ , (1)
where Var(X) denotes the variance of X.

Note that the median-of-means estimator µ (δ)

n does not require any knowledge of the variance of X. However, it depends on the desired confidence level δ and the partition B 1 , . . . , B b . Any partition satisfying ∀i, |B i | ≥ ⌊n/b⌋ is valid in order to get [START_REF] Alon | The space complexity of approximating the frequency moments[END_REF]. Hence, 2 we do not keep the dependence on the partition B 1 , . . . , B b in the notation µ (δ) n . Devroye, Lerasle, Lugosi, and Oliveira [START_REF] Devroye | Sub-Gausssian mean estimators[END_REF] introduce estimators that work for a large range of confidence levels under some mild assumptions. Catoni [START_REF] Catoni | Challenging the empirical mean and empirical variance: a deviation study[END_REF] introduces estimators of quite different flavor and gets a non-asymptotic result of the same form as [START_REF] Alon | The space complexity of approximating the frequency moments[END_REF]. Bubeck, Cesa-Bianchi and Lugosi [START_REF] Bubeck | Bandits with heavy tail[END_REF] apply these estimators in the context of bandit problems.

Estimating the mean of random vectors

Consider now the multi-dimensional case when d > 1. The sample mean µ n = (1/n) n i=1 X i is still an obvious choice for estimating the mean vector µ.

If X has a multivariate normal distribution with mean vector µ and covariance matrix Σ, then µ n is also multivariate normal with mean µ and covariance matrix (1/n)Σ and therefore, for δ ∈ (0, 1), with probability at least 1δ,

µ n -µ ≤ Tr(Σ) n + 2λ max log(1/δ) n , (2)
where Tr(Σ) and λ max denote the trace and largest eigenvalue of the covariance matrix, respectively (Hanson and Wright [START_REF] Hanson | A bound on tail probabilities for quadratic forms in independent random variables[END_REF]). For non-Gaussian and possibly heavy-tailed distributions, one cannot expect such a sub-Gaussian behavior of the sample mean. The main goal of this paper is to investigate under what conditions it is possible to define mean estimators that reproduce a (non-asymptotic) sub-Gaussian performance similar to (2). Lerasle and Oliveira [START_REF] Lerasle | Robust empirical mean estimators[END_REF], Hsu and Sabato [START_REF] Hsu | Loss minimization and parameter estimation with heavy tails[END_REF], and Minsker [START_REF] Minsker | Geometric median and robust estimation in Banach spaces[END_REF] extend the median-ofmeans estimator to more general spaces. In particular, Minsker's results imply that for each δ ∈ (0, 1) there exists a mean estimator µ (δ) n and a universal constant C such that, with probability at least 1δ,

µ (δ) n -µ ≤ C Tr(Σ) log(1/δ) n . (3)
While this bound is quite remarkable-note that no assumption other than the existence of the covariance matrix is made-, it does not quite achieve a sub-Gaussian performance bound that resembles (2). An instructive example is when all eigenvalues are identical and equal to λ max . If the dimension d is large, ( 2) is of the order of (λ max /n)(d + log(δ -1 )) while (3) gives the order (λ max /n)(d log(δ -1 )). The main result of this paper is the construction of a mean estimator that, under some mild moment assumptions, achieves a sub-Gaussian performance bound in the sense of (2). More precisely, we prove the following.

Theorem 1 For all δ ∈ (0, 1) there exists a mean estimator µ

(δ)
n and a universal constant C such that if X 1 , . . . , X n are i.i.d. random vectors in R d with mean µ ∈ R d and covariance matrix Σ such that there exists a constant

K > 0 such that, for all v ∈ R d with v = 1, E (X -µ) T v 4 ≤ K(v T Σv) 2 , 3 then for all n ≥ CK log d (d + log(1/δ)), µ (δ) n -µ ≤ C Tr(Σ) n + λ max log(log d/δ) n .
The theorem guarantees the existence of a mean estimator whose performance matches the sub-Gaussian bound (2), up to the additional term of the order of (1/n)λ max log log d for all distributions satisfying the fourth-moment assumption given above. The additional term is clearly of minor importance. (For example, it is dominated by the first term whenever Tr(Σ) > λ max log log d.) With the estimator we construct, this term is inevitable. On the other hand, the inequality of the theorem only holds for sample sizes that are at least a constant times d log d. This feature is not desirable for truly high-dimensional problems, especially taking into account that Minsker's bound is "dimension-free".

The fourth-moment assumption can be interpreted as a boundedness assumption of the kurtosis of (Xµ) T v. The same assumption has be used in Catoni [START_REF] Catoni | Pac-bayesian bounds for the gram matrix and least squares regression with a random design[END_REF] and Giulini [START_REF] Giulini | Robust dimension-free gram operator estimates[END_REF] for the robust estimation of the Gram matrix. The fourth-moment assumption may be weakened to an analogous "(2 + ε)-th moment assumption" that we do not detail for the clarity of the exposition.

We prove the theorem by constructing an estimator in several steps. First we construct an estimator that performs well for "spherical" distributions (i.e., for distributions whose covariance matrix has a trace comparable to dλ max ). This estimator is described in Section 2. In the second step, we decompose the space-in a data-dependent way-into the orthogonal sum of O(log d) subspaces such that all but one subspaces are such that the projection of X to the subspace has a spherical distribution. The last subspace is such that the projection has a covariance matrix with a small trace. In each subspace we apply the first estimator and combine them to obtain the final estimator µ (δ) n . The proof below provides an explicit value of the constant C, though no attempt has been made to optimize its value.

The constructed estimator is computationally so demanding that even for moderate values of d it is hopeless to compute it in reasonable time. In this sense, Theorem 1 should be regarded as an existence result. It is an interesting an important challenge to construct estimators with similar statistical performance that can be computed in polynomial time (as a function of n and d). Note that the estimator of Minsker cited above may be computed by solving a convex optimization problem, making it computationally feasible, see also Hsu and Sabato [START_REF] Hsu | Loss minimization and parameter estimation with heavy tails[END_REF] for further computational considerations.

An estimator for spherical distributions

In this section we construct an estimator that works well whenever the distribution of X is sufficiently spherical in the sense that a positive fraction of the eigenvalues of the covariance matrix is of the same order as λ max . More precisely, for c ≥ 1, we call a distribution c-spherical if dλ max ≤ cTr(Σ).

For each δ ∈ (0, 1) and unit vector w ∈ S d-1 (where S d-1 = {x ∈ R d : x = 1}), we may define m (δ) n (w) as the median-of-means estimate (as defined in Section 1.1) of w T µ = Ew T X based on the i.i.d. sample w T X 1 , . . . , w T X n .

Let N 1/2 ⊂ S d-1 be a minimal 1/2-cover, that is, a set of smallest cardinality that has the property that for all u ∈ S d-1 there exists w ∈ N 1/2 with ux ≤ 1/2. It is well known (see, e.g., [START_REF] Matoušek | Lectures on discrete geometry[END_REF]Lemma 13

.1.1]) that |N 1/2 | ≤ 8 d .
Noting that Var(w T X) ≤ λ max , by ( 1) and the union bound, we have that, with probability at least 1δ,

sup w∈N 1/2 m (δ/8 d ) n (w) -w T µ ≤ 2e 2λ max ln(e8 d /δ) n .
In other words, if, for λ > 0, we define the empirical polytope

P δ,λ = x ∈ R d : sup w∈N 1/2 m (δ/8 d ) n (w) -w T x ≤ 2e 2λ ln(e8 d /δ) n ,
then with probability at least 1δ, µ ∈ P δ,λmax . In particular, on this event, P δ,λmax is nonempty. Suppose that an upper bound of the largest eigenvalue of the covariance matrix λ ≥ λ max is available. Then we may define the mean estimator

µ (δ) n,λ = any element y ∈ P δ,λ if P δ,λ = ∅ 0 otherwise .
Now suppose that µ ∈ P δ,λ and let y ∈ P δ,λ be arbitrary. Define u = (y -µ)/ y -µ ∈ S d-1 , and let w ∈ N 1/2 be such that wu ≤ 1/2. (Such a w exists by definition of N 1/2 .) Then

y -µ = u T (y -µ) = (u -w) T (y -µ) + w T (y -µ) ≤ (1/2) y -µ + 4e 2λ ln(e8 d /δ) n ,
where we used Cauchy-Schwarz and the fact that y, µ ∈ P δ,λ . Rearranging, we obtain that, on the event that µ ∈ P δ,λ ,

µ (δ) n,λ -µ ≤ 8e 2λ d ln 8 + ln(e/δ) n ,
provided that λ ≥ λ max . Summarizing, we have proved the following.

Proposition 1 Let λ > 0 and δ ∈ (0, 1). For any distribution with mean µ and covariance matrix Σ such that λ max = Σ ≤ λ, the estimator µ In particular, if the distribution is c-spherical and λ ≤ 2λ max , then

µ (δ) n,λ -µ ≤ 16e
cTr(Σ) ln 8 + λ max ln(e/δ) n .

The bound we obtained has the same sub-Gaussian form as (2), up to a multiplicative constant, whenever the distribution is c-spherical. To make the estimator fully datadependent, we need to find an estimate λ that falls in the interval [λ max , 2λ max ], with high probability. This may be achieved by splitting the sample in two parts of equal size (assuming n is even), estimating λ max using samples from one part and computing the mean estimate defined above using the other part. In the next section we describe such a method as a part of a more general procedure.

Empirical eigendecomposition

In the previous section we presented a mean estimate that works well for "spherical" distributions. We will use this estimator as a building block in the construction of an estimator that has the desirable performance guarantee for distributions with any covariance matrix. In addition to finite covariances, we assume that there exists a constant

K > 0 such that, for all v ∈ R d with v = 1, E (X -µ) T v 4 ≤ K(v T Σv) 2 . (4)
In this section we assume that n ≥ 2(400e

) 2 K log 3/2 d d log 25 + log(2 log 3/2 d) + log(1/δ) .
The basic idea is the following. We split the data into two equal halves. We use the first half in order to decompose the space into the sum of orthogonal subspaces such that the projection of X into each subspace is 4-spherical. Then we may estimate the projected means by the estimator of the previous section.

Next we describe how we obtain an orthogonal decomposition of the space based on n i.i.d. observations X 1 , . . . , X n .

Let s = ⌈log 3/2 d 2 ⌉ and m = ⌊n/s⌋. Divide the sample into s blocks, each of size at least m. In what follows, we describe a way of sequentially decomposing R d into the orthogonal sum of s + 1 subspaces

R d = V 1 ⊕ • • • ⊕ V s+1 .
First we construct V 1 using the first block X 1 , . . . , X m of observations. Then we use the second block to build V 2 , and so on, for s blocks. The key properties we need are that (a) the random vector X, projected to any of these subspaces has a 4-spherical distribution; (b) the largest eigenvalue of the covariance matrix of X, projected on V i is at most λ max (2/3) i-1 .

To this end, just like in the previous section, let N γ ⊂ S d-1 be a minimal γ-cover of the unit sphere S d-1 for a sufficiently small constant γ ∈ (0, 1). The value γ = 1/100 is sufficient for our purposes and in the sequel we assume this value. Note that |N γ | ≤ (4/γ) d (see [START_REF] Matoušek | Lectures on discrete geometry[END_REF]Lemma 13.1.1] for a proof of this fact).

Initially, we use the first block X 1 , . . . , X m . We may assume that m is even. Using these observations, for each u ∈ N γ , we compute an estimate 2 , where X ′ is an i.i.d. copy of X. We may construct the estimate by forming m/2 i.i.d. random variables (1/2)(u T (X 1 -X m/2+1 )) 2 , . . . , (1/2)(u T (X m/2 -X m )) 2 and estimate their mean by the median-of-means estimate V (δ) m (u) with parameter δ/(s(4/γ) d ). Then (1), together with assumption (4) implies that, with probability at least 1δ/s,

V (δ) m (u) of u T Σu = E(u T (X -µ)) 2 = (1/2)E(u T (X -X ′ ))
sup u∈Nγ u T Σu -V (δ) m (u) u T Σu ≤ 4e K log(s(4/γ) d /δ) m def.
= ε m .

Our assumptions on the sample size guarantee that ε m < 1/100. The event that the inequality above holds is denoted by

E 1 so that P{E 1 } ≥ 1 -δ/s. Let M δ,m be the set of all symmetric positive semidefinite d × d matrices M satisfying sup u∈Nγ u T Mu -V (δ) m (u) u T Σu ≤ ε m .
By the argument above, Σ ∈ M δ,m on the event E 1 . In particular, on E 1 , M δ,m in nonempty. Define the estimated covariance matrix

Σ (δ) m = any element of M δ,m if M δ,m = ∅ 0 otherwise Since on E 1 both Σ (δ)
m and Σ are in M δ,m , on this event, we have

u T Σu 1 -ε m 1 + ε m ≤ u T Σ (δ) m u ≤ u T Σu 1 + ε m 1 -ε m for all u ∈ N γ . (5)
Now compute the spectral decomposition

Σ (δ) m = d i=1 λ i v i v T i ,
where λ 1 ≥ • • • ≥ λ d ≥ 0 are the eigenvalues and v 1 , . . . , v d the corresponding orthogonal eigenvectors.

Let u ∈ S d-1 be arbitrary and let v be a point in N γ with smallest distance to u. Then

u T Σ (δ) m u = v T Σ (δ) m v + 2(u -v) T Σ (δ) m v + (u -v) T Σ (δ) m (u -v) ≤ v T Σ (δ) m v + λ 1 (2γ + γ 2 ) ( 6
) (by Cauchy-Schwarz and using the fact that uv ≤ γ)

≤ (v T Σv) 1 + ε m 1 -ε m + 3γ λ 1 (by (5)) ≤ 1 + ε m 1 -ε m λ max + 3γ λ 1 .
In particular, on E 1 we have λ 1 ≤ βλ max where β = 1+εm 1-εm /(1 -3γ) < 1.1.

constructed on the sample X 1 , . . . , X n . Then, by Proposition 1, with probability 1δ/(s + 1), In the last subspace V s+1 , we may use Minsker's estimator, based on Π s+1 (X 1 ), . . . , Π s+1 (X n ) to compute an estimator µ s+1 = µ (δ/(s+1)) n of Π s+1 (µ). Since the largest eigenvalue of the covariance matrix of Π s+1 (X) is at most λ max /d 2 , using (3), we obtain that, with probability 1δ/(s + 1), µ s+1 -Π s+1 (µ) 2 ≤ C λ max log((2 log 3/2 d + 1)/δ) n .

µ i -Π i (µ) 2 ≤ (8e) 2
Our final estimator is µ (2/3) i-1 ≤ max .

On the other hand, since

Tr(Σ) = E X -µ 2 = s+1 i=1 E Π i (X) -Π i (µ) 2
and for i ≤ s each Π i (X) has a 4-spherical distribution, we have that This concludes the proof of Theorem 1.

  defined above satisfies, with probability at least 1δ, µ (δ) n,λµ ≤ 8e 2λ d ln 8 + ln(e/δ) n .

1 d

 1 i ln 8 + ln(e(2 log 3/2 d + 1)/δ) n .

2 ≤ (8e) 2 ( 8 1 +C

 2281 i=1 µ s+1 . By the union bound, we have that, with probability at least 1δ,µ (δ) nµ 2 = s+1 i=1 µ i -Π i (µ) λ max log((2 log 3/2 d + 1)/δ) nFirst notice that, by properties (iii) and (iv) at the end of the previous section,

* Supported by the French Agence Nationale de la Recherche (ANR), under grant ANR-13-BS01-0005 (project SPADRO).

† Supported by the Spanish Ministry of Economy and Competitiveness, Grant MTM2015-67304-P and FEDER, EU.

‡ Support from CNPq, Brazil via Ciência sem Fronteiras grant # 401572/2014-5. Supported by a Bolsa de Produtividade em Pesquisa from CNPq, Brazil. Supported by FAPESP Center for Neuromathematics (grant# 2013/ 07699-0 , FAPESP -S. Paulo Research Foundation).

By a similar argument, we have that for any u ∈ S d-1 , if v is the point in N γ with smallest distance to u, then on E 1 ,

In particular, λ max ≤ β λ 1 ≤ (4/3) λ 1 . Similarly,

Let d 1 be number of eigenvalues λ i that are at least λ 1 /2 and let V 1 be the subspace of R d spanned by v 1 , . . . , v d 1 . Denote by Π 1 (X) the orthogonal projection of the random variable X (independent of the

and therefore

In particular, the ratio of the largest and smallest eigenvalues of the covariance matrix of Π 1 (X) is at most 4 and therefore the distribution of Π 1 (X) is 4-spherical.

On the other hand, on the event E 1 , for any unit vector

(by ( 5), ( 6), and a similar argument for the last term)

In other words, the largest eigenvalue of the covariance matrix of Π ⊥ 1 (X) (the projection of X to the subspace V ⊥ 1 ) is at most (2/3)λ max .

In the next step we construct the subspace V 2 ⊂ V ⊥ 1 . To this end, we proceed exactly as in the first step but now we replace R d by V ⊥ 1 and the sample X 1 , . . . , X m on the first block by the variables Π ⊥ 1 (X m+1 ), . . . ,

. Just like in the first step, with probability at least 1δ/s we obtain a (possibly empty) subspace V 2 , orthogonal to V 1 such that Π 2 (X), the projection of X on V 2 , has a 4-spherical distribution and largest eigenvalue of the covariance matrix of Π ⊥ 2 (X) (the projection of X to the subspace

We repeat the procedure s times and use a union bound the s events. We obtain, with probability at least 1δ, a sequence of subspaces V 1 , . . . , V s , with the following properties:

(i) V 1 , . . . , V s are orthogonal subspaces.

(ii) For each i = 1, . . . , s, Π i (X), the projection of X on V i , has a 4-spherical distribution.

(iii) The largest eigenvalue of the covariance matrix of Π i (X) is at most λ

Note that it may happen for some T < s, we have

Putting it all together

In this section we construct our final multivariate mean estimator and prove Theorem 1. To simplify notation, we assume that the sample size is 2n. This only effects the value of the universal constant C in the statement of the theorem.

The data is split into two equal halves (X 1 , . . . , X n ) and (X n+1 , . . . , X 2n ). The second half is used to construct the orthogonal spaces V 1 , . . . , V s as described in the previous section. Let d 1 , . . . , d s denote the dimension of these subspaces. Recall that, with probability at least 1δ, the construction is successful in the sense that the subspaces satisfy properties (i)-(iv) described at the end of the previous section. Denote this event by E. In the rest of the argument we condition on (X n+1 , . . . , X 2n ) and assume that E occurs. All probabilities below are conditional.

If

. . , Π s+1 denote the projection operators on the subspaces V 1 , . . . , V s+1 , respectively. For each i = 1, . . . , s + 1, we use the vectors Π i (X 1 ), . . . , Π i (X n ) to compute an estimator of the mean E [Π i (X)|(X n+1 , . . . , X 2n )] = Π i (µ).

For i = 1, . . . , s, we use the estimator defined in Section 2. In particular, within the d i -dimensional space V i , we compute µ i = µ (δ/(s+1)) n,(4/3) λ i . Note that since λ i comes from an empirical estimation of Σ restricted to an empirical subspace V i , µ i is an estimator


Dmitry Ostrovsky
email: dmitrii.ostrovskii@imag.fr

Zaid Harchaoui
email: zaid.harchaoui@nyu.edu

Anatoli Juditsky
email: juditsky@imag.fr

Arkadi Nemirovski
email: nemirovs@isye.gatech.edu

Structure-Blind Signal Recovery

. However, when the set is unspecified, designing an estimator that is blind to the hidden structure of the signal remains a challenging problem. We propose a new family of estimators to recover signals observed in Gaussian noise. Instead of specifying the set where the signal lives, we assume the existence of a well-performing linear estimator. Proposed estimators enjoy exact oracle inequalities and can be efficiently computed through convex optimization. We present several numerical illustrations that show the potential of the approach.

Introduction

We consider the problem of recovering a complex-valued signal (x t ) t∈Z from the noisy observations

y τ = x τ + σζ τ , -n ≤ τ ≤ n. (1)
Here n ∈ Z + , and ζ τ ∼ CN (0, 1) are i.i.d. standard complex-valued Gaussian random variables, meaning that ζ 0 = ξ 1 0 + ıξ 2 0 with i.i.d. ξ 1 0 , ξ 2 0 ∼ N (0, 1). Our goal is to recover x t given the sequence of observations y -n , ..., y t up to instant t, a task usually referred to as (pointwise) filtering in machine learning, statistics, and signal processing [START_REF] Haykin | Adaptive filter theory[END_REF].

The traditional approach to this problem considers linear estimators, or linear filters, which write as

x t = t+n+1 τ =0
φ τ y t-τ , -n ≤ t ≤ n.

Linear estimators have been thoroughly studied in various forms, they are both theoretically attractive [START_REF] Tsybakov | Introduction to Nonparametric Estimation[END_REF][START_REF] Ibragimov | minskii. Estimation of linear functionals in gaussian noise[END_REF][START_REF] Donoho | Renormalization exponents and optimal pointwise rates of convergence[END_REF][START_REF] Donoho | Statistical estimation and optimal recovery[END_REF] and easy to use in practice. If the set X of signals is well-specified, one can usually compute a (nearly) minimax on X linear estimator in a closed form. In particular, if X is a class of smooth signals, such as a Hölder or a Sobolev ball, then the corresponding estimator is given by the kernel estimator with the properly set bandwidth parameter [START_REF] Tsybakov | Introduction to Nonparametric Estimation[END_REF] and is minimax among all possible estimators. Moreover, as shown by [START_REF] Ibragimov | Nonparametric estimation of the value of a linear functional in gaussian white noise[END_REF][START_REF] Donoho | Statistical estimation and optimal recovery[END_REF], if only X is convex, compact, and centrally symmetric, the risk of the best linear estimator of x t is within a small constant factor of the minimax risk over X . Besides, if the set X can be specified in a computationally tractable way, which clearly is still a weaker assumption than classical smoothness assumptions, the best linear estimator can be efficiently computed by solving a convex optimization problem on X . In other words, given a computationally tractable set X on the input, one can compute a nearly-minimax linear estimator and the corresponding (nearly-minimax) risk over X . The strength of this approach, however, comes at a price: the set X still must be specified beforehand. Therefore, when one faces a recovery problem without any prior knowledge of X , this approach cannot be implemented.

We adopt here a novel approach to filtering, which we refer to as structure-blind recovery. While we do not require X to be specified beforehand, we assume that there exists a linear oracle -a well-performing linear estimator of x t . Previous works [START_REF] Juditsky | Nonparametric denoising of signals with unknown local structure, I: Oracle inequalities[END_REF][START_REF] Harchaoui | Adaptive recovery of signals by convex optimization[END_REF], following a similar philosophy, proved that one can efficiently adapt to the linear oracle filter of length m = O(n) if the corresponding filter φ is time-invariant, i.e. it recovers the target signal uniformly well in the O(n)-sized neighbourhood of t, and if its 2 -norm is smallbounded by ρ/ √ m for a moderate ρ ≥ 1. The adaptive estimator is computed by minimizing the ∞ -norm of the filter discrepancy, in the Fourier domain, under the constraint on the 1 -norm of the filter in the Fourier domain. Put in contrast to the oracle linear filter, the price for adaptation is proved to be of order O(ρ 3 √ ln n), with the lower bound of O(ρ √ ln n) [START_REF] Juditsky | Nonparametric denoising of signals with unknown local structure, I: Oracle inequalities[END_REF][START_REF] Harchaoui | Adaptive recovery of signals by convex optimization[END_REF]. We make the following contributions:

• we propose a new family of recovery methods, obtained by solving a least-squares problem constrained or penalized by the 1 -norm of the filter in the Fourier domain;

• we prove exact oracle inequalities for the 2 -risk of these methods;

• we show that the price for adaptation improves upon previous works [START_REF] Juditsky | Nonparametric denoising of signals with unknown local structure, I: Oracle inequalities[END_REF][START_REF] Harchaoui | Adaptive recovery of signals by convex optimization[END_REF] to O(ρ 2 √ ln n) for the point-wise risk and to O(ρ √ ln n) for the 2 -risk.

• we present numerical experiments that show the potential of the approach on synthetic and real-world images and signals.

Before presenting the theoretical results, let us introduce the notation we use throughout the paper.

Filters. Let C(Z) be the linear space of all two-sided complex-valued sequences x = {x t ∈ C} t∈Z . For k, k ∈ Z we consider finite-dimensional subspaces

C(Z k k ) = {x ∈ C(Z) : x t = 0, t / ∈ [k, k ]} .
It is convenient to identify m-dimensional complex vectors, m = k -k + 1, with elements of C(Z k k ) by means of the notation:

x k k := [x k ; ...; x k ] ∈ C k -k+1 .
We associate to linear mappings

C(Z k k ) → C(Z j j ) (j -j + 1) × (k -k + 1) matrices with complex entries. The convolution u * v of two sequences u, v ∈ C(Z) is a sequence with elements [u * v] t = τ ∈Z u τ v t-τ , t ∈ Z.
Given observations (1) and ϕ ∈ C(Z m 0 ) consider the (left) linear estimation of x associated with filter ϕ: Discrete Fourier transform. We define the unitary Discrete Fourier transform (DFT) operator F n :

x t = [ϕ * y] t ( x t is
C n+1 → C n+1 by z → F n z, [F n z] k = (n + 1) -1/2 n t=0 z t e 2πıkt n+1 , 0 ≤ k ≤ n.
The inverse Discrete Fourier transform (iDFT) operator F -1 n is given by F -1 n := F H n (here A H stands for Hermitian adjoint of A). By the Fourier inversion theorem, F -1 n (F n z) = z. We denote • p usual p -norms on C(Z):

x p = ( t∈Z |x t | p ) 1/p , p ∈ [1, ∞].
Usually, the argument will be finite-dimensional -an element of C(Z k k ); we reserve the special notation x n,p := x n 0 p . Furthermore, DFT allows to equip C(Z n 0 ) with the norms associated with p -norms in the spectral domain:

x * n,p := x n 0 * p := F n x n 0 p , p ∈ [1, ∞]
; note that unitarity of the DFT implies the Parseval identity: x n,2 = x * n,2 . Finally, c, C, and C stand for generic absolute constants.

2 Oracle inequalities for adaptive recovery procedures

Constrained recovery

Given observations [START_REF] Bhaskar | Atomic norm denoising with applications to line spectral estimation[END_REF] and > 0, we first consider the constrained recovery x con given by

[ x con ] t = [ ϕ * y] t , t = 0, ..., n,
where ϕ is an optimal solution of the constrained optimization problem min

ϕ∈C(Z n 0 ) y -ϕ * y n,2 : ϕ * n,1 ≤ / √ n + 1 . (2)
The constrained recovery estimator minimizes a least-squares fit criterion under a constraint on ϕ * n,1 = F n ϕ n 0 1 , that is an 1 constraint on the discrete Fourier transform of the filter. While the least-squares objective naturally follows from the Gaussian noise assumption, the constraint can be motivated as follows.

Small-error linear filters. Linear filter ϕ o with a small 1 norm in the spectral domain and small recovery error exists, essentially, whenever there exists a linear filter with small recovery error [START_REF] Juditsky | Nonparametric denoising of signals with unknown local structure, I: Oracle inequalities[END_REF][START_REF] Harchaoui | Adaptive recovery of signals by convex optimization[END_REF]. Indeed, let us say that x ∈ C(Z n 0 ) is simple [START_REF] Harchaoui | Adaptive recovery of signals by convex optimization[END_REF] with parameters m ∈ Z + and ρ ≥ 1 if there exists

φ o ∈ C(Z m 0 ) such that for all -m ≤ τ ≤ 2m, E |x τ -[φ o * y] τ | 2 1/2 ≤ σρ √ m + 1 . (3)
In other words, x is (m, ρ)-simple if there exists a hypothetical filter φ o of the length at most m + 1 which recovers x τ with squared risk uniformly bounded by

σ 2 ρ 2 m+1 in the interval -m ≤ τ ≤ 2m. Note that (3) clearly implies that φ o 2 ≤ ρ/ √ m + 1, and that |[x -φ o * x] τ | ≤ σρ/ √ m + 1 ∀τ, -m ≤ τ ≤ 2m. Now, let n = 2m, and let ϕ o = φ o * φ o ∈ C n+1 .
As proved in Appendix C, we have

ϕ o * n,1 ≤ 2ρ 2 / √ n + 1, (4)
and, for a moderate absolute constant c,

x -ϕ o * y n,2 ≤ cσρ 2 1 + ln[1/α] (5)
with probability 1 -α. To summarize, if x is (m, ρ)-simple, i.e., when there exists a filter φ o of length ≤ m + 1 which recovers x with small risk on the interval [-m, 2m], then the filter ϕ o = φ o * φ o of the length at most n + 1, with n = 2m, has small norm ϕ o * n,1 and recovers the signal x with (essentially the same) small risk on the interval [0, n].

Hidden structure. The constrained recovery estimator is completely blind to a possible hidden structure of the signal, yet can seamlessly adapt to it when such a structure exists, in a way that we can rigorously establish. We formalize the hidden structure as an unknown shift-invariant linear subspace of C(Z), ∆S = S, of a small dimension s. We do not assume that x belongs to that subspace. Instead, we make a more general assumption that x is close to this subspace, that is, it may be decomposed into a sum of a component that lies in the subspace and a component whose norm we can control.

Assumption A We suppose that x admits the decomposition

x = x S + ε, x S ∈ S,
where S is (an unknown) shift-invariant, ∆S = S, subspace of C(Z) of dimension s, 1 ≤ s ≤ n + 1, and ε is "small", namely,

∆ τ ε n,2 ≤ σκ, 0 ≤ τ ≤ n.
Examples. Shift-invariant subspaces of sequences C(Z) are actually sets of solutions of difference equations

s i=0 p i x τ -i = 0, τ ∈ Z. (6)
For instance, discrete-time polynomials x τ = s i=0 a i τ i , τ ∈ Z of degree s form a linear space of dimension s + 1 of solutions of the equation [START_REF] Donoho | Renormalization exponents and optimal pointwise rates of convergence[END_REF] with binomial coefficients p i = (-1) i n i . Harmonic oscillations are another example: x τ = Ce ıωτ for some frequency ω ∈ [0, 2π) is the set of solutions of the equation x τ -e ıω x τ -1 = 0.

We can now state an oracle inequality for the constrained recovery estimator; see Appendix B for the proof.

Theorem 2.1. Let ≥ 1, and let ϕ o ∈ C(Z n 0 ) be such that

ϕ o * n,1 ≤ / √ n + 1.
Suppose that Assumption A holds for some s ∈ Z + and κ < ∞. Then for any α, 0 < α ≤ 1, it holds with probability at least 1 -α:

x -

x con n,2 ≤ x -ϕ o * y n,2 + Cσ s + κ ln [1/α] + ln [n/α] . (7)
When considering simple signals, Theorem 2.1 gives the following.

Corollary 2.1. Assume that signal x is (m, ρ)-simple, ρ ≥ 1 and m ∈ Z + . Let n = 2m, ≥ 2ρ 2 , and let Assumption A hold for some s ∈ Z + and κ < ∞. Then for any α, 0 < α ≤ 1, it holds with probability at least 1 -α:

x -

x con n,2 ≤ Cσρ 2 ln[1/α] + C σ s + κ ln [1/α] + ln [n/α] .
Adaptation and price. The price for adaptation in Theorem 2.1 and Corollary 2.1 is determined by three parameters: the bound on the filter norm , the deterministic error κ, and the subspace dimension s. Assuming that the signal to recover is simple, and that = 2ρ 2 , let us compare the magnitude of the oracle error to the term of the risk which reflects "price of adaptation". Typically (in fact, in all known to us cases of recovery of signals from a shift-invariant subspace), the parameter ρ is at least √ s. Therefore, the bound (5) implies the "typical bound" O(σ √ γρ 2 ) = σs √ γ for the term x -ϕ o * y n,2 (we denote γ = ln(1/α)). As a result, for instance, in the "parametric situation", when the signal belongs or is very close to the subspace, that is when κ = O(ln(n)), the price of adaptation O σ[s + ρ 2 (γ + √ γ ln n)] 1/2 is much smaller than the bound on the oracle error. In the "nonparametric situation", when κ = O(ρ 2 ), the price of adaptation has same order of magnitude as the oracle error. Finally, note that under the premise of Corollary 2.1 we can also bound the pointwise error; see Appendix B for the proof. We state the result for = 2ρ 2 for simplicity. Theorem 2.2. Assume that signal x is (m, ρ)-simple, ρ ≥ 1 and m ∈ Z + . Let n = 2m, = 2ρ 2 , and let Assumption A hold for some s ∈ Z + and κ < ∞. Then for any α, 0 < α ≤ 1, the constrained recovery x con satisfies

|x n -[ x con ] n | ≤ C σρ √ m + 1 ρ 2 ln[n/α] + ρ κ ln [1/α] + √ s .

Penalized recovery

The constrained recovery estimator assumes that the parameter is known. If the noise variance is known (or can be estimated from data), then we can build a more practical estimator that still enjoys an oracle inequality.

The penalized recovery estimator [ x pen ] t = [ ϕ * y] t is an optimal solution to a regularized least-squares minimization problem, where the regularization penalizes the 1 -norm of the filter in the Fourier domain:

ϕ ∈ Argmin ϕ∈C(Z n 0 ) y -ϕ * y 2 n,2 + λ √ n + 1 ϕ * n,1 . (8)
We establish an oracle inequality for the penalized recovery estimator, similarly to Theorem 2.1:

Theorem 2.3. Let Assumption A hold for some s ∈ Z + and κ < ∞, and let

ϕ o ∈ C(Z n 0 ) satisfy ϕ o * n,1 ≤ / √ n + 1 for some ≥ 1.
1 o . Suppose that the regularization parameter of penalized recovery x pen satisfies λ ≥ λ,

λ := 60σ 2 ln[63n/α].
Then, for 0 < α ≤ 1, it holds with probability at least 1 -α:

x -

x pen n,2 ≤ x -ϕ o * y n,2 + C λ + C σ s + ( + 1)κ ln[1/α],
where

:= √ n + 1 ϕ * n,1 . 2 o . Moreover, if κ ≤ κ, κ := 10 ln[42n/α] ln [16/α] ,
and λ ≥ 2λ, one has

x -x pen n,2 ≤ x -ϕ o * y n,2 + C λ + C σ √ s.
The proof of Theorem 2.3 closely follows that of Theorem 2.1 and is omitted.

Discussion

There is some redundancy between "simplicity" of a signal, as defined by (3), and Assumption A. Usually a simple signal or image x is also close to a low-dimensional subspace of C(Z) (see, e.g., [13, section 4]), so that Assumption A holds "automatically". Likewise, x is "almost" simple when it is close to a low-dimensional time-invariant subspace. Indeed, if x ∈ C(Z) belongs to S, i.e. Assumption A holds with κ = 0, one can easily verify that for n ≥ s there exists a filter

φ o ∈ C(Z n -n ) such that φ o 2 ≤ s/(n + 1), and x τ = [φ o * x] τ , τ ∈ Z . (9)
See Appendix C for the proof. This implies that x can be recovered efficiently from observations (1):

E |x τ -[φ o * y] τ | 2 1/2 ≤ σ s n + 1 .
In other words, if instead of the filtering problem we were interested in the interpolation problem of recovering x t given 2n + 1 observations y t-n , ..., y t+n on the left and on the right of t, Assumption A would imply a kind of simplicity of x. On the other hand, it is clear that Assumption A is not sufficient to imply the simplicity of x "with respect to the filtering", in the sense of the definition we use in this paper, when we are allowed to use only observations on the left of t to compute the estimation of x t . Indeed, one can see, for instance, that already signals from the parametric family X α = {x ∈ C(Z) : x τ = cα τ , c ∈ C}, with a given |α| > 1, which form a one-dimensional space of solutions of the equation x τ = αx τ -1 , cannot be estimated with small risk at t using only observations on the left of t (unless c = 0), and thus are not simple in the sense of (3). Of course, in the above example, the "difficulty" of the family X α is due to instability of solutions of the difference equation which explode when τ → +∞. Note that signals x ∈ X α with |α| ≤ 1 (linear functions, oscillations, or damped oscillations) are simple. More generally, suppose that x satisfies a difference equation of degree s:

0 = p(∆)x τ = s i=0 p i x τ -i , (10)
where p(z) = s i=0 p i z i is the corresponding characteristic polynomial and ∆ is the right shift operator. When p(z) is unstable -has roots inside the unit circle -(depending on "initial conditions") the set of solutions to the equation ( 10) contains difficult to filter signals. Observe that stability of solutions is related to the direction of the time axis; when the characteristic polynomial p(z) has roots outside the unit circle, the corresponding solutions may be "left unstable" -increase exponentially when τ → -∞. In this case "right filtering" -estimating x τ using observations on the right of τ -will be difficult. A special situation where interpolation and filtering is always simple arises when the characteristic polynomial of the difference equation has all its roots on the unit circle. In this case, solutions to [START_REF] Ibragimov | minskii. Estimation of linear functionals in gaussian noise[END_REF] are "generalized harmonic oscillations" (harmonic oscillations modulated by polynomials), and such signals are known to be simple. Theorem 2.4 summarizes the properties of the solutions of (10) in this particular case; see Appendix C.1 for the proof.

Theorem 2.4. Let s be a positive integer, and let p = [p 0 ; ...; p s ] ∈ C s+1 be such that the polynomial p(z) = s i=0 p i z i has all its roots on the unit circle. Then for every integer m satisfying m ≥ m(s) := Cs 2 ln(s + 1), one can point out q ∈ C m+1 such that any solution to (10) satisfies

x τ = [q * x] τ , ∀τ ∈ Z, and q 2 ≤ ρ(s, m)/ √ m where ρ(s, m) = C min s 3/2 √ ln s, s ln[ms] . (11)

Numerical experiments

We present a preliminary simulation study of the proposed recovery algorithms in several application scenarios. We illustrate the performance of constrained 2 -recovery of Sec. 2.1 and the penalized 2 -recovery of Sec. 2.2 and compare it to that of the constrained ∞ -recovery of [START_REF] Harchaoui | Adaptive recovery of signals by convex optimization[END_REF] and the Lasso recovery of [START_REF] Bhaskar | Atomic norm denoising with applications to line spectral estimation[END_REF] for signal and image denoising. Some implementation details are discussed in Appendix D. Details and discussion of the discretization approach underlying the competing Lasso method can be found in Sec. 3.6 of [START_REF] Bhaskar | Atomic norm denoising with applications to line spectral estimation[END_REF]. First, we investigate the behavior of the different methods on a simple one-dimensional signal denoising example.

Signal denoising

We consider two different scenarios. In Scenario 1, the signal is a sum of k = 4 sine waves of equal amplitudes and random frequencies. In Scenario 2, k/2 = 2 pairs of close frequencies are sampled, with the frequencies in each pair separated by only 0.1 of the DFT bin 2π/n, making recovery harder due to high signal coherency. The signal-to-noise ratio (SNR) is defined as the ratio x2 n,2 /(σ 2 n). In both scenarios, we perform 100 Monte-Carlo trials for a signal-to-noise ratio in the range 10 -1 , . . . , 10 1 . The parameter of the constrained 2 -and constrained ∞ -recoveries was set = 4. We use the same regularization parameter of the penalized 2 -recovery as in the dimension reduction example above. For the Lasso, we use the theoretically recommended value [START_REF] Bhaskar | Atomic norm denoising with applications to line spectral estimation[END_REF]. As Fig. 1 suggests, by comparing results with = 4 and = 16, the constrained 2 -recovery method is in fact not too sensitive to the choice of (k), where k is the number of sines.

Image Denoising

We now consider the recovery of an unknown regression function f on the regular grid T on [0, 1] 2 given noisy observations:

y τ = f (τ ) + σζ τ , τ ∈ T := 1 n , 2 n , ..., 1 2 . ( 12
)
The signal-to-noise ratio (SNR) is defined here as x 2 n,2 /(σ 2 n).

Denoising textures. In this experiment, we apply the proposed recovery methods to denoise two images from the original Brodatz texture database, observed in white Gaussian noise.

The signal-to-noise ratio (SNR) is set to SNR = 1. We use the blockwise implementation of the constrained 2 -recovery algorithm, as described in Appendix D. We set the constraint parameter to = 4, and we use the adaptation procedure of [START_REF] Harchaoui | Adaptive recovery of signals by convex optimization[END_REF]Sec. 3.2] to define the estimation bandwidth. We use the Lasso method of [START_REF] Bhaskar | Atomic norm denoising with applications to line spectral estimation[END_REF] with λ = 0.5σ √ 2 log N , halving the theoretically recommended value to prevent edge oversmoothing. The resulting images are presented in Fig. 2. Despite comparable quality in the mean square sense, the two methods significantly differ in their local behaviour. In particular, constrained 2 -recovery better restores the local signal features, whereas the Lasso tends to oversmooth, even with the chosen value of the regularization penalty.

Dimension reduction. The purpose of this experiment is to illustrate the performance of the penalized

f (t) = g(θ T t), g(•) ∈ S 1 β (L), θ 1 = 1. ( 13
)
True signal Observations Constrained recovery Lasso Here, S 1 β (L) = {g : R → R, g (β) (•) 2 ≤ L} is the Sobolev ball of smooth periodic functions on [0, 1], and θ, called the index, is some unknown direction.

Note that if it is known a priori that the regression function possesses the structure ( 13), but the index is unknown, one can adapt to it, obtaining the one-dimensional rates of recovery; see e.g. [START_REF] Lepski | Adaptive estimation under single-index constraint in a regression model[END_REF] In our experiments, a random direction θ is sampled uniformly from the unit sphere and then renormalized. We use the blockwise strategy, with one block corresponding to the entire image. For penalized Denoising harmonic oscillations. We consider the recovery of noisy images which are sums of k = 4 harmonic oscillations in R2 with random frequencies. We compare constrained 2 -recovery (2), with a single block and the constraint parameter set to (k) = k 2 , to the Lasso, with the regularization penalty parameter set to the theoretically recommended value [START_REF] Bhaskar | Atomic norm denoising with applications to line spectral estimation[END_REF]. The signal-to-noise ratio (SNR) is set to SNR = 0.25. We present the results in Fig. 4.

True signal

Observations

Penalized recovery Lasso

A Preliminaries

We start with introducing several objects used in the sequel. We denote Convolution matrices. We will extensively use the matrix-vector representation of the discrete convolution.

• Given y ∈ C(Z), we associate to it an (n + 1) × (m + 1) Toeplitz matrix

T (y) =     y 0 y -1 ... y -m y 1 y 0 ... y 1-m ... ... ... ... y n y n-1 ... y n-m     . ( 14
)
such that [ϕ * y] n 0 = T (y)ϕ m 0 for ϕ ∈ C(Z m 0 ). Its squared Frobenius norm satisfies

T (y) 2 F = m τ =0 ∆ τ y 2 n,2 . (15)
• Given ϕ ∈ C(Z m 0 ), consider an (n + 1) × (m + n + 1) matrix

M (ϕ) =     ϕ m ϕ m-1 ... ϕ 0 0 0 ... 0 0 ϕ m ... ...
For y ∈ C(Z) we have [ϕ * y] n 0 = M (ϕ)y n -m and

M (ϕ) 2 F = (n + 1) ϕ 2 m,2 . (17)
• Given ϕ ∈ C(Z m 0 ), consider the following circulant matrix of size m + n + 1:

C(ϕ) =           ϕ m ... ...
          . ( 18
)
One has C(ϕ)

2 F = (m + n + 1) ϕ 2 m,2
. This matrix is useful since C(ϕ)y n -m encodes the circular convolution of y n -m and the zero-padded filter ϕ m+n 0 (recall that ϕ ∈ C(Z m 0 )) which is diagonalized by the DFT. Specifically,

C(ϕ) = F H m+n D(ϕ)F m+n , where D(ϕ) = √ m + n + 1 diag(F m+n ϕ m+n 0 ). ( 19
)
Deviation bounds. We use the following simple facts about Gaussian random vectors.

• Let ζ ∼ CN (0, I n ) be a standard complex Gaussian vector meaning that ζ = ξ 1 + ıξ 2 where ξ 1,2 are two independent draws from N (0, I n ). We will use a simple bound

Prob ζ ∞ ≤ √ 2 ln n + 2u ≥ 1 -e -u (20)
which may be checked by explicitly evaluating the distribution since |ζ 1 | 2 2 ∼ Exp(1/2).

• The following deviation bounds for ζ 2 2 ∼ χ 2 2n are due to [15, Lemma 1]:

Prob ζ 2 2 2 ≤ n + √ 2nu + u ≥ 1 -e -u , Prob ζ 2 2 2 ≥ n - √ 2nu ≥ 1 -e -u . (21)
By simple algebra we obtain an upper bound for the norm:

Prob ζ 2 ≤ √ 2n + √ 2u ≥ 1 -e -u . (22)
• Further, let K be an n × n Hermitian matrix with the vector of eigenvalues λ = [λ 1 ; ...

Prob ζ H Kζ 2 ≤ Tr(K) + (u + √ 2u) K F ≥ 1 -e -u . (23)
Further, when K is positive semidefinite, we have

K F ≤ Tr(K), whence Prob ζ H Kζ 2 ≤ Tr(K)(1 + √ u) 2 ≥ 1 -e -u . (24)
B Proof of Theorems 2.1 and 2.2

B.1 Proof idea

Despite the striking similarity with the Lasso [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF], [START_REF] Candes | The Dantzig selector: statistical estimation when p is much larger than n[END_REF], [START_REF] Bickel | Simultaneous analysis of lasso and dantzig selector[END_REF], the recoveries of section 2 are of quite different nature. First of all, the 1 -minimization in these methods is aimed to recover a filter but not the signal itself, and this filter is not sparse. 1 The equivalent of "regression matrices" involved in these methods cannot be assumed to satisfy Restricted Eigenvalue or Restricted Isometry conditions, usually imposed to prove statistical properties of "classical" 1 -recoveries (see e.g. [START_REF] Candes | Stable signal recovery from incomplete and inaccurate measurements[END_REF], [START_REF] Van De Geer | On the conditions used to prove oracle results for the lasso[END_REF], and references therein). Moreover, being constructed from the noisy signal itself, these matrices depend on the noise, what introduces an extra degree of complexity in the analysis of the properties of these estimators. Yet, proofs of Theorem 2.1 and 2.3 rely on some simple ideas and it may be useful to expose these ideas stripped from the technicalities of the complete proof. Given y ∈ C(Z n -n ) let T (y) be the (n + 1) × (n + 1) "convolution matrix," as defined in [START_REF] Juditsky | On detecting harmonic oscillations[END_REF] such that for ϕ ∈ C(Z n 0 ) [ϕ * y] n 0 = T (y)ϕ n 0 . When denoting f = F n ϕ, the optimization problem in (2) can be recast as a "standard" 1 -constrained least-squares problem with respect to f :

min f ∈C n+1 y -A n f 2 2 : f 1 ≤ / √ n + 1 ( 25
)
1 Unless we consider recovery of signal composed of harmonic oscillations with frequencies on the DFT grid.

where

A n = T (y)F -1 n . Observe that f o = F n ϕ o is feasible in (25) so that y -A n f 2 n,2 ≤ y -A n f o 2 n,2 ,
where f = F n ϕ, so that

x -A n f 2 n,2 -x -A n f o 2 n,2 ≤ 2σ ζ, x -A n f o n -ζ, x -A n f n ≤ 2σ ζ, A n (f o -f ) n ≤ 2σ A H n ζ n 0 ∞ f o -f 1 ≤ 4σ A H n ζ n 0 ∞ √ n+1
In the "classical" situation, where ζ n 0 is independent of A n (see, e.g., [START_REF] Juditsky | Functional aggregation for nonparametric regression[END_REF]), the norm c α is a logarithmic in α -1 factor. This would rapidly lead to the bound (7) of the theorem. In the case we are interested in, where A n incorporates observations y n -n and thus depends on ζ n 0 , curbing the cross term is more involved and requires extra assumptions, e.g. Assumption A.

A H n ζ n 0 ∞ is bounded by c α √ ln n max j [A n ] j 2 ≤ c α √ n ln n A n ∞ where A ∞ = max i,j |A ij | and

B.2 Extended formulation

We will prove a simple generalization of the theorem in the case where the length n of "validation sample" may be different from the length m of the adjusted filter. We consider the "skewed" sample

y τ = x τ + σζ τ -m ≤ τ ≤ n, (26)
with κ m,n := n+1 m+1 . Accordingly, we assume that the regular recovery ϕ * y uses the filter

ϕ ∈ Argmin ϕ∈C(Z m 0 ) y -ϕ * y n,2 : ϕ * m,1 ≤ / √ m + 1 . (27)
The corresponding modification of Assumption A is as follows:

Assumption A Let S be a (unknown) shift-invariant linear subspace of C(Z), ∆S = S, of dimension s, 1 ≤ s ≤ n + 1. We suppose that x admits the decomposition:

x = x S + ε,
where x S ∈ S, and ε is "small", namely,

∆ τ ε n,2 ≤ σκ, 0 ≤ τ ≤ m.
In what follows we use the following convenient reformulation of Assumption A (reformulation of Assumption A when m = n):

There exists an s-dimensional (complex) subspace S n ⊂ C n+1 and an idempotent Hermitian (n + 1) × (n + 1) matrix Π Sn of rank s -the projector on S n -such that

(I n+1 -Π Sn ) [∆ τ x] n 0 2 = ∆ τ ε n,2 ≤ σκ, τ = 0, ..., m (28)
where I n+1 is the (n + 1) × (n + 1) identity matrix.

Theorem B.1. Let m, n ∈ Z + , ≥ 1, and let ϕ o ∈ C(Z m 0 ) be such that ϕ o * 1 ≤ / √ m + 1.
Suppose that Assumption A holds for some s ∈ Z + and κ < ∞. Then for any α, 0 < α ≤ 1, there is a set

Ξ ⊂ C 2n+1 , Prob{ζ n -m ∈ Ξ} ≥ 1 -α, of "good realisations" of ζ such that whenever ζ n -m ∈ Ξ, x -ϕ * y n,2 ≤ x -ϕ o * y n,2 + 2σ V 2 α + ( + 1)c α κ + √ 2s + c α , (29)
where c α := 2 ln[16/α], and

V 2 α = 2 (1 + 4κ m,n ) 2 ln [55(m + n + 1)/α] .
B.3 Proof of Theorem B.1

1 o . The oracle filter ϕ o is feasible in (27), hence, x -ϕ * y 2 n,2 = (1 -ϕ o ) * y 2 n,2 -σ 2 ζ 2 n,2 -2 σ ζ, x -ϕ * y n = x -ϕ o * y 2 n,2 -2 σ ζ, x -ϕ * y n δ (1) +2 σ ζ, x -ϕ o * y n δ (2) . ( 30
)
Let us bound δ (1) . Denote for brevity I := I n+1 , and recall that Π Sn is the projector on S n from (28). We have the following decomposition:

δ (1) = σ ζ n 0 , Π Sn [x -ϕ * y] n 0 δ (1) 1 + σ ζ n 0 , (I -Π Sn )[x -ϕ * x] n 0 δ (1) 2 -σ 2 ζ n 0 , (I -Π Sn )[ ϕ * ζ] n 0 δ (1) 3 (31)
One can easily bound δ

1 under the premise of the theorem:

δ (1) 1 ≤ σ Π Sn ζ n 0 2 Π Sn [x -ϕ * y] n 0 2 ≤ σ Π Sn ζ n 0 2 x -ϕ * y n,2 .
Note that Π Sn ζ n 0 ∼ CN (0, I s ), and by ( 22) we have

Prob Π Sn ζ n 0 2 ≥ √ 2s + √ 2u ≤ e -u , obtaining the bound Prob δ (1) 1
≤ σ x -ϕ * y n,2 √ 2s + 2 ln [1/α 1 ] ≥ 1 -α 1 . (32)
2 o . We are to bound the second term of (31). To this end, note first that

δ (1) 2 = σ ζ n 0 , (I -Π Sn )x n 0 -σ ζ n 0 , (I -Π Sn )[ ϕ * x] n 0 .
By (28), (I -Π Sn )x n 0 2 ≤ σκ, thus with probability 1 -α,

| ζ n 0 , (I -Π Sn )x n 0 | ≤ σκ 2 ln[1/α]. (33)
On the other hand, using the notation defined in ( 14), we have [ ϕ * x] n 0 = T (x) ϕ m 0 , so that

ζ n 0 , (I -Π Sn )[ ϕ * x] n 0 = ζ n 0 , (I -Π Sn )T (x) ϕ m 0 . Note that [T (x)] τ = x -τ +n -τ
for the columns of T (x), 0 ≤ τ ≤ m. By (28), (I -Π Sn )T (x) = T (ε), and by ( 15),

(I -Π Sn )T (x) 2 F = T (ε) 2 F = m τ =0 ε -τ +n -τ 2 2 ≤ (m + 1)σ 2 κ 2 .
Due to (24) we conclude that

T (x) H (I -Π Sn )ζ n 0 2 2 ≤ 2(m + 1)σ 2 κ 2 1 + ln[1/α] 2 with probability at least 1 -α. Since | ζ n 0 , (I -Π Sn )T (x) ϕ m 0 | ≤ √ m + 1 T (x) H (I -Π Sn )ζ n 0 2 ,
we arrive at the bound with probability 1 -α:

| ζ n 0 , (I -Π Sn )T (x) ϕ m 0 | ≤ √ 2σκ 1 + ln[1/α] .
Along with (33) this results in the following bound:

Prob δ (1) 2 ≤ √ 2σ 2 κ( + 1) 1 + ln [2/α 2 ] ≥ 1 -α 2 . ( 34
)
3 o . Let us rewrite δ (1)
3 as follows:

δ (1) 3 = σ 2 ζ n 0 , (I -Π Sn )M ( ϕ)ζ n -m = σ 2 ζ n -m , QM ( ϕ)ζ n -m ,
where M ( ϕ) ∈ C (n+1)×(m+n+1) is defined by ( 16), and Q ∈ C (m+n+1)×(n+1) is given by

Q = 0 m,n+1 ; I -Π Sn ;
hereinafter we denote 0 m,n the m × n zero matrix. Now, by the definition of ϕ and since the mapping

ϕ → M (ϕ) is linear, δ (1) 3
= 1 2 (ζ n -m ) H (QM ( ϕ) + M ( ϕ) H Q H K1( ϕ) )z n -m ≤ σ 2 2 √ m + 1 max u ∈ C(Z m 0 ), u * m,1 ≤ 1 (ζ n -m ) H K 1 (u)ζ n -m = σ 2 √ m + 1 max 1≤j≤m+1 max θ∈[0,2π] 1 2 (ζ n -m ) H K 1 (e ıθ u j )ζ n -m , (35)
where u j ∈ C(Z m 0 ), and [u j ] m 0 = F -1 m e j , e j being the j-th canonical orth of R m+1 . Indeed, M (ϕ) attains its maximum over the convex set

B * m,1 = {u ∈ C(Z m 0 ), u * m,1 ≤ 1}. (36)
at an extremal point e ıθ u j , θ ∈ [0, 2π]. It is easy to verify that

K 1 (e ıθ u) = K 1 (u) cos θ + K 2 (u) sin θ for the Hermitian matrix K 2 (u) = ı QM (u) -M (u) H Q H . Denoting q j i (ζ) = 1 2 (ζ n -m ) H K i (u j )ζ n -m for i = 1, 2, we have max θ∈[0,2π] 1 2 (ζ n -m ) H K 1 (e ıθ u j )ζ n -m = max θ∈[0,2π] q j 1 (ζ) cos θ + q j 2 (ζ) sin θ = q j 1 (ζ) 2 + q j 2 (ζ) 2 ≤ √ 2 max q j 1 (ζ) , q j 2 (ζ) . (37)
By simple algebra and using (17), we get

Tr K i (u j ) 2 ≤ 4 Tr[M (u j )M (u j ) H ] = 4(n + 1) u j 2 m,2 ≤ 4(n + 1), i = 1, 2.
Now let us bound Tr[K i (u)], i = 1, 2, on the set (36). One may check that for the circulant matrix C(u), cf. [START_REF] Nesterov | On first-order algorithms for 1 /nuclear norm minimization[END_REF], it holds:

QM (u) = QQ H R C(u),
where R = QQ H is an (m + n + 1) × (m + n + 1) projection matrix of rank s defined by

R = 0 m,m 0 m,n+1 0 n+1,m I -Π Sn .
Hence, denoting • * the nuclear norm, we can bound Tr[K i (u)], i = 1, 2, as follows:

Tr[K i (u)] ≤ 2 Tr[RC(u)] ≤ 2 R C(u) * ≤ 2 C(u) * = 2 √ m + n + 1 u * m+n,1 ,
where in the last transition we used [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF]. The following technical lemma gives an upper bound on the norm of a zero padded filter (see Appendix B.4 for the proof):

Lemma B.1. For any u ∈ B * m,1 , see (36), and n ≥ 1, we have

u * m+n,1 ≤ 1 + κ 2 m,n (ln[m + n + 1] + 3).
Thus we arrive at

Tr[K i (u j )] ≤ 2 √ m + 1(κ 2 m,n + 1)(ln[m + n + 1] + 3), i = 1, 2.
By (23) we conclude that for any fixed pair (i, j) ∈ {1, 2} × {1, ..., m + 1}, with probability 1 -α,

q j i (ζ) ≤ Tr[K i (u j )] + K i (u j ) F 1 + ln[2/α] 2 .
With α 0 = 2(m + 1)α, by the union bound together with ( 35) and (37) we get Prob δ

(1)

3 ≤ 2 √ 2σ 2 (κ 2 m,n + 1)(ln[m + n + 1] + 3) + κ m,n 1 + ln [4(m + 1)/α 0 ] 2 ≥ 1 -α 0 . ( 38
)
4 o . Bounding δ (2) is a relatively simple task since ϕ o does not depend on the noise. We decompose

δ (2) = σ ζ, x -ϕ o * x n -σ 2 ζ, ϕ o * ζ n . Note that ζ, x -ϕ o * x n ∼ CN (0, x -ϕ o * x 2 n,2 ), therefore, with probability 1 -α, ζ, x -ϕ o * x n ≤ 2 ln[1/α] x -ϕ o * x n,2 . (39)
On the other hand,

x -ϕ o * x n,2 ≤ x -ϕ o * y n,2 + σ ϕ o * ζ n,2 ≤ x -ϕ o * y n,2 + √ 2σ κ m,n 1 + ln[1/α] (40)
with probability 1 -α. Indeed, one has

ϕ o * ζ 2 n,2 = M (ϕ o )ζ n -m 2 2
,

where for M (ϕ o ) by ( 17) we have

M (ϕ o ) 2 F = (n + 1) ϕ o 2 m,2 ≤ κ 2 m,n 2 . ( 41
)
Using (24) we conclude that, with probability at least 1 -α,

ϕ o * ζ 2 n,2 ≤ 2κ 2 m,n 2 1 + ln[1/α] 2 ,
which implies (40). Using (39) and (40), we get that with probability at least 1 -α 3 ,

ζ, x -ϕ o * x n ≤ 2 ln [2/α 3 ] x -ϕ o * y n,2 + √ 2σ κ m,n 1 + ln [2/α 3 ] ≤ x -ϕ o * y n,2 2 ln [2/α 3 ] + 2σ κ m,n 1 + ln [2/α 3 ] 2 . ( 42
)
The indefinite quadratic form

ζ, ϕ o * ζ n = (ζ n -m ) H K 0 (ϕ o )ζ n -m 2 ,
where

K 0 (ϕ o ) = [0 m,m+n+1 ; M (ϕ o )] + [0 m,m+n+1 ; M (ϕ o )]
H , can be bounded similarly to 3 o . We get

| Tr[K 0 (ϕ o )]| ≤ 2(n + 1) |ϕ o m | ≤ 2κ 2 m,n
. Indeed, for e m+1 = [0; ...; 0; 1] ∈ R m+1 one has

|ϕ o m | = | [ϕ o ] m 0 , e m+1 | ≤ ϕ o * m,1 F m e m+1 ∞ ≤ m + 1 since F m e m+1 ∞ = 1/ √ m + 1. By (41), K 0 (ϕ o ) 2 F ≤ 4 M (ϕ o ) 2 F ≤ 4κ 2 m,n 2 .
Hence by (23),

Prob -ζ, ϕ o * ζ n ≤ 2κ 2 m,n + 2κ m,n 1 + 2 ln [1/α 4 ] 2 ≥ 1 -α 4 . ( 43
)
5 o . It remains to combine the bounds obtained in 1 o -4 o . For any α ∈ (0, 1], putting α 0 = α 1 = α 4 = α/4, α 2 = α 3 = α/8, and using the the union bound, we get from (30) with probability ≥ 1 -α:

x -ϕ * y 2 n,2 ≤ x -ϕ o * y 2 n,2 + 2δ (2) -2δ (1) [by ( 42 +υ (1) + υ (2) .

The latter implies that

x -ϕ * y n,2 ≤ x -ϕ o * y n,2 + 2 √ 2σ √ s + ln[16/α] + υ (1)
+ υ (2) .

Finally, we arrive at (29) using the bounds

υ (1) ≤ 4 √ 2σ 2 ( +

B.4 Proof of Lemma B.1

The function u * m+n,1 is convex on (36), so its maximum over this set is attained at one the extreme points F m [u j ] m 0 = e ıθ e j where e j is the j-th canonical orth of R m+1 and θ ∈ [0, 2π]. Since u j τ = 1 √ m+1 exp ıθ -2πıτ j m+1 , we obtain

u j * m+n,1 = 1 γ m+n k=0 m τ =0 exp ı 2π k m + n + 1 - j m + 1 ω jk τ = 1 γ m+n k=0 D m ω jk 2 ,
where γ = (m + n + 1)(m + 1), and the Dirichlet kernel D m (•) is defined as

D m (x) := sin((m+1)x) sin(x) , x = πl, m + 1, x = πl.
Hence, γ u j * m+n,1 ≤ max ∈[0,π] S m+n ( ), where

S m,n ( ) = m+n k=0 D m πk m + n + 1 - , ∈ [0, π]. (44)
Note that |D m (x)| is upper bounded by the following (positive) function on the circle R/πZ: where in the last transition we used the simple bound H n ≤ ln n + 1 for harmonic numbers.

B m (x) = π 2 min(x,π-x) , x ∈ (0, π), m + 1, x = 0.

B.5 Proof of Theorem 2.2

We decompose

|[x -ϕ * y] n | = |[(φ o + (1 -φ o )) * (x -ϕ * y)] n | ≤ |[φ o * (x -ϕ * y)] n | + |[(1 -ϕ) * (1 -φ o ) * x] n | + σ|[ ϕ * ζ] n | + σ|[ ϕ * φ o * ζ] n | := δ (1) + δ (2) + δ (3) + δ (4) . (45)
We have δ

(1) ≤ φ o 2 x -ϕ * y) m,2 ≤ ρ √ m + 1 x -ϕ * y m,2 .
When using the bound of Corollary 2.1 with = 2ρ 2 , we conclude that, with probability ≥ 1 -α/3,

δ (1) ≤ c σρ √ n ρ 2 ln[1/α] + ρ κ ln 1/α + ln n/α + √ s .
Next we get

δ (2) ≤ (1 + ϕ 1 ) (1 -φ o ) * x n,∞ ≤ (1 + 2ρ 2 ) σρ √ m + 1
.

By the Parseval identity,

δ (3) = σ| F n [ ϕ * ] n 0 , F n [ζ] 0,n | ≤ σ ϕ * n,1 ζ * n,∞ ≤ 2σρ 2 √ m + 1 σ 2 ln [3(n + 1)/α],
where the last inequality, holding with probability ≥ 1 -α/3, is due to [START_REF] Tsybakov | Introduction to Nonparametric Estimation[END_REF]. Finally, observe that, with probability ≥ 1 -α/3 (cf. ( 46)),

φ o * ζ n,2 ≤ √ 2ρ 1 + ln[3/α] .
Therefore, we have for δ (4) :

δ (4) ≤ σ ϕ n,2 φ o * ζ n,2 ≤ σ 2ρ 2 √ m + 1 √ 2ρ 1 + ln[3/α] = σ 2 √ 2ρ 3 √ m + 1 1 + ln[3/α]
with probability 1 -α/3. When substituting the bound for δ (k) , k = 1, ..., 4, into (45) we arrive at the result of the theorem.

C Miscellaneous proofs

Proof of relations (4) and (5

). Let n = 2m, φ ∈ C(Z m 0 ), and let ϕ ∈ C(Z n 0 ) satisfy ϕ = φ * φ. Then ϕ * n,1 = (2m + 1) -1/2 2m k=0 |[F 2m ϕ 2m 0 ] k | = √ 2m + 1 2m k=0 |[F 2m φ 2m 0 ] k | √ 2m + 1 2 = √ 2m + 1 φ * 2 2m,2 = √ 2m + 1 φ 2 2m,2 = √ 2m + 1 φ 2 m,2 ≤ √ 2m + 1 ρ 2 m + 1 , implying (4)
. Moreover, since 1 -φ * φ = (1 + φ) * (1 -φ), for all x ∈ C(Z) one has for all τ ∈ Z:

|x τ -[ϕ o * x] τ | = |[(1 + φ o ) * (1 -φ o ) * x] τ | = m j=0 [1 + φ o ] j [x -φ o * x] τ -j ≤ 1 + φ o 1 max 0≤j≤m |[x -φ o * x] τ -j | ≤ σ(1 + ρ)ρ √ m + 1
(we have used (3) to obtain the last inequality), and

x -[ϕ o * x] n,2 ≤ σ(1 + ρ)ρ. Next note that ϕ o * ζ 2 n,2 = ζ, M (ϕ o )ζ n
, where M (ϕ) is defined as in [START_REF] Lepski | Adaptive estimation under single-index constraint in a regression model[END_REF]. When taking into account that

ϕ 2 ≤ ϕ * n,1 ≤ 2ρ 2 √ n + 1 n + 2 , we get (cf. ( 17
)) M (ϕ o ) 2 F = (n + 1) ϕ 2 2 ≤ 4ρ 4
, so that the concentration inequality (24) now implies that, given 0 < α ≤ 1, with probability at least 1 -α,

ϕ o * ζ n,2 ≤ 2 √ 2ρ 2 1 + ln[1/α] , (46)
and we arrive at [START_REF] Donoho | Statistical estimation and optimal recovery[END_REF].

Proof of [START_REF] Ibragimov | Nonparametric estimation of the value of a linear functional in gaussian white noise[END_REF]. Assume that Assumption A holds true for some n ≥ s and ε ≡ 0. Let Π Sn be the Euclidean projector on the space S n of elements of S restricted on C(Z n 0 ). Since dim(S n ) ≤ s, Π Sn C.1 Proof of Theorem 2.4

Note that to prove the lemma we have to exhibit a vector q ∈ C m+1 of small 2 -norm and such that the polynomial 1 -q(z) = 1 -m i=0 q i z i is divisible by p(z), i.e., that there is a polynomial r(z) of degree m -s such that 1 -q(z) = r(z)p(z).

Indeed, this would imply that

x t -[q * x] t = [1 -q(∆)]x t = r(∆)p(∆)x t = 0 due to p(∆)x t = 0,

The bound q 2 ≤ C s 3/2 ln s m of ( 11) is proved in [14, Lemma 6.1]. Our objective is to prove the "remaining" inequality

q 2 ≤ C s ln[ms] m .
So, let θ 1 , ..., θ s be complex numbers of modulus 1 -the roots of the polynomial p(z). Given δ = 1-∈ (0, 1), let us set δ = 2δ/(1 + δ), so that δ

δ -1 = 1 -δ > 0. ( 47
)
Consider the function q(z) = s i=1 z -θ i δz -θ i .

Note that q(•) has no singularities in the circle B = {z : |z| ≤ 1/ δ}; besides this, we have q(0) = 1.Let |z| = 1/ δ, so that z = δ-1 w with |w| = 1. We have

|z -θ i | |δz -θ i | = 1 δ |w -δθ i | |w -δ δ θ i | .
We claim that when |w| = 1, |w -δθ i | ≤ |w -δ δ θ i |. Indeed, assuming w.l.o.g. that w is not proportional to θ i , consider triangle ∆ with the vertices A = w, B = δθ i and C = δ δ θ i . Let also D = θ i . By (47), the segment AD is a median in ∆, and ∠CDA is ≥ π 2 (since D is the closest to C point in the unit circle, and the latter contains A), so that |w -δθ i | ≤ |w -δ δ θ i |.

Alternatively, one may split the signal into blocks, and process the points of each block using the same filter (cf. e.g. Theorem 2.1). For instance, a one-dimensional signal can be divided into blocks of length, say, 2m + 1, and to recover x ∈ C(Z m -m ) in each block one may fit one filter of length m + 1 recovering the right "half-block" x m 0 and another filter recovering the left "half-block" x -1 -m . When recovering a signal or an image, the "optimal" filter bandwidth -the length of the filter for which the balance of the stochastic and the approximation error is attained -depends on the local signal parameters (e.g. smoothness) and vary from one point to another. In the experiments of Sec. 3, we implement pointwise and blockwise bandwidth adaptation using a procedure similar in spirit to the celebrated Lepski adaptation procedure [START_REF] Lepskii | On a problem of adaptive estimation in gaussian white noise[END_REF]. Details are given in e.g. [START_REF] Juditsky | Nonparametric denoising signals of unknown local structure, II: Nonparametric function recovery[END_REF]Sec. 3.3] and [START_REF] Harchaoui | Adaptive recovery of signals by convex optimization[END_REF]Sec. 3.2].

Note that the optimization problems ( 2) and ( 8) underlying recovery algorithms are well structured Second-Order Conic Programs (SOCP) and can be solved using Interior-point methods (IPM). However, the computational complexity of IPM applied to SOCP with dense matrices grows rapidly with problem dimension, so that large problems of this type arising in signal and image processing are well beyond the reach of these techniques. On the other hand, these problems possess nice geometry associated with complex

  recovery, rho=4 Constrained l2-recovery, rho=16 = 4 vs. = 16.

Figure 1 :

 1 Figure 1: Denoising harmonic oscilations.

Figure 2 :

 2 Figure 2: Recovery of two instances of the Original Brodatz database, cut by half to 320 × 320 and observed with SNR = 1. 2 -error: 1.35e4 for the constrained 2 -recovery 1.25e4 for the Lasso in the first row (inst. D73); 1.97e4 for the constrained 2 -recovery, 2.02e4 for the Lasso method in the second row (inst. D75).

Figure 3 :Figure 4 :

 34 Figure 3: Recovery of the single-index signal (13), observed with SNR = 1, for β = 2 (first row) and β = 1 (second row). Results are in Tab. 1.

  •, • the Hermitian scalar product: for a, b ∈ C n , a, b = a H b. For a, b ∈ C(Z) we reserve the shorthand notation a, b n := a n 0 , b n 0 = [a n 0 ] H b n 0 .

2 F = 2 K 2 F

 222 ; λ n ]. Then the real-valued quadratic form ζ H Kζ has the same distribution as ξ T Bξ, where ξ = [ξ 1 ; ξ 2 ] ∼ N (0, I 2n ), and B is a real 2n × 2n symmetric matrix with the vector of eigenvalues [λ; λ]. Hence, we have Tr(B) = 2Tr(K), B and B = K ≤ K F , where • and • F denote the spectral and the Frobenius norm of a matrix. Invoking [15, Lemma 1] again (a close inspection of the proof shows that the assumption of positive semidefiniteness can be relaxed), we have

For

  any ∈ [0, π], the summation in (44) is over a regular (m + n + 1)-grid on R/πZ. The contribution to the sum of each of two closest to zero points of the grid is at most max x∈[0,π] D m (x) = m + 1. For the remaining points, we can upper bound |D m (x)| ≤ B m (x) noting that B m (x) decreases over [ π m+n+1 , π 2 ] as long as n ≥ 1. These considerations result in S m,n ( ) ≤ 2(m + 1) +

2 2 =

 2 Tr(Π Sn ) ≤ s, there is ι ∈ {0, ..., n} such that the ι-th column [Π Sn ] ι of Π Sn satisfies [Π Sn ] ι 2 ≤ s/(n + 1). Note that one has x ι -[Π Sn ] ι , x n 0 = 0, and S n is time-invariant, implying thatx τ -[Π Sn ] ι , x τ -ι+n τ -ι = 0, ∀τ ∈ Z.We conclude that there isφ o ∈ C(Z n -n ),φ o = [0; ...; 0; [Π Sn ] ι ; 0; ...; 0] (i.e., vector [Π Sn ] ι completed with zeros in such a way that ι-th element of [Π Sn ] ι becomes the central (n + 1-th) entry of φ o ) such that φ o 2 ≤ s/(n + 1), and x τ -[φ o * x] τ = 0, ∀τ ∈ Z.

Table 1 :

 1 and reference therein. 2 -error in the dimension reduction experiment.

	β Observations Penalized recovery Lasso
	0.5	45.30	54.22	54.94
	1.0	64.71	17.05	16.94
	2.0	65.66	4.02	7.49
	3.0	65.63	2.72	6.95

  )] ≤ x -ϕ o * y 2 n,2 + 2σ x -ϕ o * y n,2 2 ln[16/α] -ϕ * y n,2 + x -ϕ o * y n,2

	[by (42), (43)]		+4σ 2 κ 2 m,n + 2κ m,n 1 + 2 ln[16/α]	2
	[by (32)] [by (34)]		+2σ x -ϕ * y n,2 +2 √ 2σ 2 ( + 1) 1 + ln[16/α] κ √ 2s + 2 ln[16/α]
	[by (38)]		+4	√	2σ 2 (κ 2 m,n + 1)(ln[m + n + 1] + 3) + κ m,n 1 + ln [16(m + 1)/α]	2	.
	Hence, denoting υ (1) = 2	√	2σ 2 ( + 1) 1 + ln[16/α] κ,
	υ (2) = 4	√	m,n + 1)(ln[m + n + 1] + 3) + κ m,n 1 + ln [16(m + 1)/α] 2σ 2 (κ 2	2
	+4σ 2 κ 2 m,n + 2κ m,n 1 + 2 ln[16/α]	2
	we obtain x -ϕ * y 2 n,2 ≤	x -ϕ o * y 2 n,2 + 2σ	√	2s + 2 ln[16/α] x

-recovery in the problem of estimating a function with a hidden structure. We consider the single-index model of the regression function f :

-recovery[START_REF] Haykin | Adaptive filter theory[END_REF], the regularization penalty parameter is set to λ = 2σ 2 log(21/α) with a confidence level 1 -α = 0.9. For Lasso, we use the same parameter setting as previously. The signal-to-noise ratio (SNR) is set to SNR = 1. The corresponding results are presented in Tab. 1 and Fig.

-norm. Moreover, their first-order information -the value of objective and its gradient at a given ϕcan be computed using Fast Fourier Transform is time which is almost linear in problem size. Therefore, we used first-order optimization algorithms, such as Mirror-Prox (MP) and Nesterov's accelerated gradient algorithms (see[START_REF] Nesterov | On first-order algorithms for 1 /nuclear norm minimization[END_REF] and reference therein), in our recovery implementation. A complete description of the application of these optimization algorithms to our problems is beyond the scope of the paper; we shall present it elsewhere.

The first, second and third authors were supported by the LabEx PERSYVAL-Lab (ANR-11-LABX-0025) and the project Titan (CNRS-Mastodons). The second author was also supported by the project Macaron (ANR-14-CE23-0003-01) and the program "Learning in Machines and Brains" (CIFAR). Research of the fourth author was supported by NSF grants CMMI-1262063, CCF-1523768.

As a consequence, we get (δz

on the boundary of B. It follows that the coefficients r j of the Taylor series of r satisfy

for |z| ≤ 1, utilizing the trivial upper bound |p(z)| ≤ 2 s , we get

Note that q (0) = p(0)r (0) = p(0)r(0) = 1, that q is a polynomial of degree + s, and that q is divisible by p(z). Besides this, on the unit circumference we have, by (51),

We can upper-bound R:

Now, given positive integer and positive α such that

Furthermore,

When invoking (52) and utilizing ( 56) and ( 55) we get

On the other hand, denoting by q 0 , q 1 ,...,q +s the coefficients of the polynomial q and taking into account that q0 = q (0) = 1, we have

We are done: when denoting m = + s, and q(z) = m i=1 q j z j , we have the vector of coefficients q = [0; q 1 ; ...; q m ] ∈ C m+1 of q(z) such that, by (57),

and such that the polynomial q (z) = 1 + q(z) is divisible by p(z) due to (50).

D Comments on algorithm implementation

In Sec. 2 we considered only recoveries which estimated the value x t of the signal via the observations at n + 1 points t -n, ..., t "on the left" (filtering problem). To recover the whole signal, one may consider a more flexible alternative -interpolating recovery -which estimates x t using observations on the left and on the right of t. In particular, if the objective is to recover a signal on the interval {-n, ..., n}, one can apply interpolating recoveries using the same observations y -n , ..., y n to estimate x τ for each τ ∈ {-n, ..., n}. Ideally, when using pointwise recovery, a specific filter is constructed for each time instant t. This may pose a tremendous amount of computation, for instance, when recovering a high-resolution image.