Greedy Variance Estimation for the LASSO

arXiv:1803.10878v1 [math.ST] 28 Mar 2018

Christopher Kennedy∗ and Rachel Ward†
Department of Mathematics, University of Texas at Austin
March 30, 2018

Abstract
Recent results have proven the minimax optimality of LASSO and related algorithms for noisy linear
regression. However, these results tend to rely on variance estimators that are inefficient or optimizations
that are slower than LASSO itself. We propose an efficient estimator for the noise variance in high
dimensional linear regression that is significantly faster than LASSO, only requiring p matrix-vector
multiplications. We prove this estimator is consistent with a good rate of convergence, under the condition
that the design matrix satisfies the Restricted Isometry Property (RIP). In practice, our estimator vastly
outperforms state of the art methods in terms of speed while incurring only a modest bias.

1

Introduction

The LASSO [22] is a classical algorithm for doing noisy linear regression in the case when the number
of regression coefficients is larger than the number of response variables p > n. The analysis of LASSO
has recently surged with much work on establishing oracle inequalities for ℓ2 estimation over sparse
vectors, and corresponding minimax rates. Typically, such results rely on knowledge of the variance of
the noise, which is unknown in practice. The full extent of the literature on LASSO is immense and
beyond the scope of this paper, but we point to a few important references on the oracle inequalities and
corresponding minimax error rates (see [3], [15], [30], [23], [17], [31], [27], [13], [29], [26], [4]).
A good review of variance estimators for LASSO is given in [19], where variance estimation using
cross-validated LASSO is highlighted as particularly strong in many sparsity regimes. This method
typically uses 5 or 10-fold cross-validation to train the hyperparameters in LASSO and analysis relies
on the restricted eigenvalue condition on the design matrix. The above work was later complemented
by a theoretical analysis of a slightly modified variant of cross-validated LASSO in [6] (see also [8] [9],
e.g.). The main drawback of existing methods such as cross-validation for variance estimation is that
they take multiple LASSO estimators to compute. This is costly when the dimensionality of the data
is large. The method of moments (see [7]) is a fast alternative to cross-validated LASSO. It relies on
the assumption that the design matrix is Gaussian and exploits statistical properties to formulate an
estimator. It is consistent with a good rate of convergence [7], but the design matrix has to be Gaussian
which is restrictive. We should also mention a variant of the LASSO - the square-root LASSO (see [2])
- whose penalty level doesn’t depend on the variance of the noise. However, the resulting estimator is
formulated as a conic programming problem which can be inefficient in practice and is beyond the scope
of this work.

1.1

Our Contribution

The main contributions of our a paper are the following:
• We provide an efficient variance estimator. In fact, our variance estimator only requires p matrixvector multiplications. This is a significant improvement over more standard estimators that iterate
the LASSO multiple times with various penalty levels. Our method is faster than a single iteration
of LASSO.
∗ Email:
† Email:

ckennedy@math.utexas.edu.
rward@math.utexas.edu.

1

• Our estimator is consistent in the sense that it converges in probability to the true variance. We
have a quantitative bound on the rate of convergence.
• We only require a deterministic assumption on the design matrix (the Restricted Isometry Property)
which holds with high probability over many standard matrix ensembles over appropriate parameter
regimes. In particular, the condition holds for any orthonormal design matrix.
Remark 1. Precisely, the term consistent means that for some ρ, c ≥ 0,
sup kb
σ 2 − σ 2 k = 0.

lim

n/p→ρ 0≤σ 2 ≤c
n,p→∞

Moreover, our estimator admits a surprisingly simple theoretical argument for convergence using
standard compressed sensing-type results and concentration of measure. We note that we assume the
design matrix satisfies the restricted isometry property, which is stronger than the restricted eigenvalue
condition typically considered in the literature (see [24] e.g.) but still admits a wide range of random
matrix ensembles.
In practice, although our estimator exhibits a higher bias than more standard estimators like crossvalidated LASSO, it does well in the high-dimensional regime where the cross-validated LASSO and even
the method of moments estimator become prohibitive to compute.

1.2

Notation

For a matrix X ∈ Rn×p , and a subset Ω ⊂ {1, .., p}, XΩ ∈ Rn×|Ω| will denote the restriction of X to its
columns indexed by Ω.
For a vector v ∈ Rp , Ωv is defined to be the support of v.
For each j = 1, ..., p/L, we use Ωj := {(j − 1)L + 1, ..., jL} to denote the jth “window” of the signal.

2

Problem Statement

Suppose β ∈ Rp is s-sparse, and that we are given a noisy, transformed version of this signal:
y = Xβ + η,
where η ∈ Rn has i.i.d. Gaussian entries ηj ∼ N (0, σ 2 ) and X ∈ Rn×p is a known design matrix. For
the purpose of analysis, we define a notion of a well-behaved design matrix X. We will assume that the
matrix X satisfies the Restricted Isometry Property (RIP), which was introduced in [5] and is a common
property used in Compressed Sensing. It guarantees that a matrix is a near-isometry on sparse vectors.
Specifically, we say X satisfies the RIP of order s and level δ > 0 if for all z such that kzk0 ≤ s,
(1 − δ)kzk22 ≤ kXzk22 ≤ (1 + δ)kzk22 .
With high probability, RIP is obtained with
s = s(m, N ) := O(n/ log(p/n))

(1)

on an n × p matrix X whose entries Xi,j are independent realizations of a subgaussian random variable,
such as a Gaussian or Bernoulli random variable [1]. The Restricted Isometry Property is obtained
with high probability on many classes of structured random matrices, such as random partial Fourier
matrices [18, 20], but with a slightly smaller (by factors of log(p)) constant s.
It is in general nontrivial to recover the true signal β. We consider here the standard LASSO algorithm
[22] to return a denoised version of β:
βb = arg min kXβ − yk22 + 2λkβk1 .
β

(2)

The magnitude of the parameter λ in the objective (2) controls the balance between the ℓ1 term which
b and a mean squared error term kXβ − yk22 which promotes
promotes sparsity in the recovered signal β,
consistency with the observed measurements. It is important to balance the two terms appropriately so
that one doesn’t overfit to the transformed signal y but also doesn’t over-enforce sparsity. The standard
analysis of the LASSO is conditioned on the event {λ : λ/4 ≥ kX T ηk∞ /n} (see [3]). In particular,
for the case that η is Gaussian with variance σ 2 and X is orthogonal, with high probability we have

2

kX T k∞ /n = Θ(σ 2 log(n)/n). Thus, with the choice λ = 4σ 2 log(n)/n, the LASSO will provably produce
a good estimate β.
However, in applications, the variance σ, and hence a proper choice of λ, is not known a priori. We
consider the case where σ is not known in advance, and needs to be estimated from the signal y. It should
be clear from the above observations that precision in estimating the parameter σ improves recovery of
the true signal.

3

Greedy Variance Estimation – The Orthonormal Case

For the moment we focus on the case where X ∈ Rp×p is an orthonormal matrix (p = n) and the problem
reduces to recovering the noisy signal y = β +η (by rotational invariance of the Gaussian). In this regime,
the LASSO has the closed form solution
βbi = sign(yi )(|yi | − λ)+ ,

where βbi = βbi (λ) implicitly depends on λ. A standard approach is to minimize the cross-validation error:
b
min ky − β(λ)k
2,
λ

which has nice practical and theoretical properties (see [11] e.g.). Moreover, given the optimal λ one can
infer a good estimate of the variance as kβb − yk2 /p. However, this approach still requires one to compute
the LASSO minimizer over a range of λ values, whereas one would like to perform a single computation
to estimate the variance (and thus optimal λ). We formulate a method to estimate the variance which
only needs a single pass over the input y.

Algorithm 1 Greedy Variance Estimator – Orthonormal Design Matrix
1 P
2
j ∈ {1, 2, . . . , p/L}.
1: Compute the window estimators Sj = L
i∈Ωj |yi | ,
P
p/(2L)
1
2L
c
2
2: Let σ = (1 + log(p) ) p
S(j) , where {S(j) }j is a non-decreasing arrangement of {Sj }j .
j=1
The basic idea behind the above algorithm is that we want to capture a noise estimator that avoids the
entries of y affected by signal (hence in the second step we take the average of the smaller 50% of the
1
window estimates). We multiply the resulting estimator by 1 + log(p)
to correct the downward bias that
results from averaging only over the smallest windows.
Remark 2. (Total variation denoising) Suppose we receive image-type data and instead of taking the
LASSO minimizer we want to instead want to regularize by the total variation seminorm:
βb = arg min kβ − yk22 + 2λTV(β),
β

(3)

P
where TV(β) := n kβn − βn−1 k. The typical assumption in this model is that the discrete derivative of
true signal is sparse, which is promoted by the above objective. In this case, we can apply our estimator
to the discrete derivative (which as observed is essentially a sparse signal plus noise) to get a reasonable
estimate of the variance of the noise in this setting. This approach originally appeared in [21] and
statistical guarantees on the resulting estimator βb have been developed in [14], [16], [28], culminating
most recently in [10]. These papers give a framework that allows one to generalize the estimator (3) to
when the signal is 2-D image data. We note that our estimators can also be easily adapted to 2-D image
data by replacing window estimates with box estimates.
c2 .
We have the following result which guarantees accuracy of the estimator σ

Theorem 3. Suppose y = Xβ + η where X ∈ Rp×p is orthonormal, ηj ∼ N (0, σ 2 ) are independent, and
p
β is s-sparse. Consider window size L ≥ log3 (p), and suppose that s ≤ 2L
. Then the Greedy Variance
Estimator produced by Algorithm 3 satisfies

with probability 1 − p2 .

c2 − σ 2 | ≤
|σ

3

6
σ2,
log p

4

Greedy Variance Estimation – RIP Design Matrix

We now turn to the more general case where the design matrix X ∈ Rn×p is possibly underdetermined
n ≤ p, but satisfies the Restricted Isometry Property with the appropriate constants (indeed this is a
more general case, as an orthonormal matrix satisfies the RIP with constant δ = 0). We define the
regularized design matrix as Z := [ZΩ1 , ..., ZΩp/L ] where each ZΩi ∈ Rn×L ,
ZΩi := Ui In×L Vi

such that

(4)

XΩi = Ui Σi Vi is the SVD of XΩi .
Then, we run a conditioning step based on the (block orthonormal) matrix Z and then run the algorithm
similar to the orthonormal case: In practice, we use the matrix X instead of Z, however using Z allows

Algorithm 2 Greedy Variance Estimator
1: Compute ỹ = Z T y.
1 P
2
j ∈ {1, 2, . . . , n/L}.
2: Compute the window estimators Sj = L
i∈Ωj |ỹi | ,
P
p/(2L)
c2 = (1 + 1 ) 2L
3: Let σ
S(j) , where {S(j) }j is a non-increasing arrangement of the window
j=1
log(p) p
estimators {Sj }j .
us to do a more streamlined theoretical analysis. To see why this should work intuitively, assume that we
precondition just on X that satisfies RIP for a large enough sparsity level s0 . Note that X T y = X T Xβ +
X T η, so the obstruction to estimating the noise is the X T X term. Then, kXβk2 = kXΩβ βk2 ≈ kβk2 ,
and if we assume our window set Ωj is disjoint from Ωβ , RIP implies the restricted matrices XΩTj , XΩβ
satisfy kXΩTj XΩβ k ≤ δ for δ > 0 small. Thus, for a “good” window estimator, we only see the noise X T η.

Theorem 4. Suppose y = Xβ + η where X ∈ Rn×p has the RIP of order s and level δ, ηj ∼ N (0, σ 2 )
n
are i.i.d., and β is s-sparse. Assume that L ≥ log3 (p), n ≥ L, s ≤ 2L
, and that X satisfies (RIP) with
order s = 2 max{L, s} and level δ > 0. Then, the variance estimator from the above algorithm satisfies




√
p
6σ 2
1
kβk2
1
+
+ max 4σ 2 log(p), 8 δσkβk2 log(p)
2δ
σ
b2 − σ 2 ≤ 1 +
log(p)
L
log(p)
L
with probability 1 − p4 .

Remark 5. The constants in Theorem 4 are chosen for neatness of presentation and are in no way
optimized.
Remark 6. Although the right hand side of Theorem 4 contains factors involving kβk2 (as opposed to
kβk1 which one finds in typical LASSO results), we do not expect this to be a problem in practice. In
particular, one can assume cσ ≤ |βj | ≤ Cσ for all j and some absolute constants C, c > 0. If the |βj |
are below this threshold, they are essentially noise and difficult to detect in general (this is called the
beta-min assumption). On the other hand, one can naturally expect the√entries of β to have a uniform
upper bound even as the problem size goes to infinity. Since kβk2 ≤ s Cσ, we just need that δ < 1s
which will hold for our sparsity regime and standard matrix models (i.i.d. normalized Gaussian entries,
for example) with high probability.

5

Experiments

Our experimental methodology is based off of the results in [19]. In particular, we generate a design
matrix X ∈ Rn×p with i.i.d. entries Xij ∼ N (0, n−1/2 ) so that X satisfies RIP with sufficiently small
constants with high probability. The sparsity level s = ⌈nα ⌉, with α < 1, and the non-zero entries of
β (chosen uniformly at random) are distributed according to a Laplace(1) distribution. The resulting β
is scaled to have the specified norm. The experiments are over the following grid of parameter values,
where n = 100 in all experiments.
• p = 100, 200, 500, 1000,
• kβk2 = 0.1, 1, 2, 5, 10,

• α = 0.1, 0.3, 0.5, 0.7, 0.9.

We use the following estimators in our analysis:

4

√
• oracle: the oracle estimator β̂ = kηk2 / n.

• window: the standard window estimator with the transformation ỹ = X T y.

• window-svd: the theoretical window estimator with the transformation ỹ = Z T y where Z is given
by (4).
• cv-lasso: 10-fold cross-validated LASSO.

• moment: method of moments estimator (see [7]).

We include the cross-validated LASSO because it was shown to be the most robust to changes in sparsity/dimension by [19] and the method of moments estimator because it is a fast alternative. The window
size is chosen based on an inflection point in the values of the estimator for a specific set of parameters
as the window size varies.
As we can see in Figure 1, the window and window-svd estimators have reasonable performance compared to the cv-LASSO with slightly larger biases. In particular, we do quite well for α = 0.1, β = 1,
performing similarly to cv-Lasso, and with a much smaller variance than the method of moments.
Remark 7. We only include results for α = 0.1 because the algorithm performs similarly for α ≤ 0.5.
Moreover our theory only covers up to roughly α = 0.5 for reasonable choices of window size. The
performance for dense signal α = 0.9 is covered in its own section below.
It is important to note that our estimators are significantly faster to compute than cv-LASSO. The
method of moments estimator is generally accurate but has unreasonably high variance for large p. We
now report computation times for each estimator averaged over 100 trials with n = 100, α = 0.9, β = 1,
p = 100, 1000, 10000, 100000. These experiments were run on an Intel i5 3.80 GHz, 4 core CPU with
16 GB of ram. As expected, the window and window-svd estimators scale incredibly well, whereas the

p
window
window-svd
moment
cv-lasso

100
0.0049
0.1559
0.0387
556.5134

1000
0.0048
1.3158
3.2771
186.1420

10000
0.0433
11.8787
1.1632e+03
N/A

100000
0.4227
135.1963
N/A
N/A

Table 1: Computation times for each estimator averaged over 100 trials.
method of moments estimator is fast for small scale problems but does not scale beyond p = 10000 and
the cv-lasso does not scale beyond p = 1000.

5.1

Optimal Window Size

It is notable to see how well our method can perform when the window size is optimized. Here, we give
some representative plots (Figure 2) to show what happens to performance when replacing the window
size with the optimal window size using prior knowledge of the variance. In all experiments, n=100 and
p=1000. For the low SNR regimes, we see a similar downward bias to the oblivious choice of window
size, although with a smaller bias. Similarly, for high SNR, the upward bias is also smaller than when
choosing an oblivious window size. In table 2 we report the optimal window size for various values of
α and kβk2 . The optimal window size was found by a grid search over all possible window sizes using
knowledge of the true variance.

5.2

High Dimension

In this section we highlight the regime in which our estimator is most useful - when p ≫ n is large. In
particular, we chose n = 100, p = 100000 in all experiments. In this setting all other estimators cannot
be computed efficiently. In this regime, it is inefficient to even compute an optimal box size based on an
inflection point in the value of the estimator, so instead the choice L = 25 was fixed for all experiments.
The results are shown in Figure 3. Although the bias remains, the estimator performs well, especially in
low SNR regimes. This is likely due to the strength of the compressed sensing properties for the design
matrix as the dimension grows.

5

1.6

1.6

1.4

1.4

1.2

1.2

1

1

0.8

0.8

0.6

0.6

0.4

0.4
oracle

window

window-svd

moment

cv-lasso

oracle

window

window-svd

moment

cv-lasso

oracle

window

window-svd

moment

cv-lasso

2.5
1.6
2
1.4
1.5

1.2

1

1

0.8

0.5

0.6
0
0.4
-0.5
oracle

window

window-svd

moment

cv-lasso

Figure 1: Window size chosen based on inflection point, p = 1000. Signal-less (kβk = 0), low SNR (α = 0.1,
kβk = 1), medium SNR (α = 0.1, kβk = 5), high SNR (α = 0.1, kβk = 10) respectively, left to right, top to
bottom.

6

1.5
1.4

2

1.3
1.2

1.5

1.1
1
1
0.9
0.8
0.5
0.7
0.6
0
0.5
oracle

window

window-svd

moment

cv-lasso

1.6

1.6

1.4

1.4

1.2

1.2

1

1

0.8

0.8

0.6

0.6

0.4

0.4
oracle

window

window-svd

moment

cv-lasso

oracle

window

window-svd

moment

cv-lasso

oracle

window

window-svd

moment

cv-lasso

Figure 2: Optimal window size, p = 1000. Left to right, top to bottom: Signal-less (kβk = 0), low SNR
(α = 0.1, kβk = 1), medium SNR (α = 0.1, kβk = 5), high SNR (α = 0.1, kβk = 10) respectively.

7

1.6

1.6

1.4

1.4

1.2

1.2

1

1

0.8

0.8

0.6

0.6

0.4

0.4
oracle

window

window-svd

oracle

window

window-svd

oracle

window

window-svd

2
1.6

1.8
1.6

1.4

1.4
1.2
1.2
1

1
0.8

0.8
0.6
0.6

0.4
0.2

0.4

0
oracle

window

window-svd

Figure 3: High-dimensional regime, p = 100000 L = 25. Left to right, top to bottom: Signal-less (kβk = 0),
low SNR (α = 0.1, kβk = 1), medium SNR (α = 0.1, kβk = 5), high SNR (α = 0.1, kβk = 10) respectively.

8

α

0.1
0.3
0.5
0.7
0.9

0.1
100
100
100
100
100

1
100
100
100
100
100

kβk2
2
100
100
100
100
100

5
20
22
18
18
14

10
4
4
4
4
3

Table 2: Optimal window sizes as a function of α and kβk2 for p = 200. We note that the optimal window
size is generally decreasing as a function of both the signal to noise ratio and the sparsity. Moreover, choosing
the maximal window size is optimal in modest regimes.

5.3

Orthogonal Design Matrix

We find our estimator performs quite well in the case where the design matrix is orthogonal, as shown in
Figure 4. In all experiments, p = n = 200 and the window size is chosen via inflection point in the value of
the estimator. The method of moments still performs reasonably well, but suffers a strong upwards bias
for large SNR. We note that in all regimes, our estimator performs better than cross-validated LASSO.
Moreover, it is more robust to changes in SNR than when the design matrix is RIP (but not necessarily
orthogonal).

5.4

Dense Signal

Our theory does not cover high sparsity levels (α ≥ 0.9), but nonetheless our estimator performs well.
Although more prone to high levels of SNR, we are still competitive with cv-LASSO in low SNR regimes
as seen in Figure 5.

6

Future Work

Our estimator has been shown to be a useful tool to use in high dimensional variance estimation, and
comes with nice theoretical properties that leverage results from the compressed sensing literature. Moreover, it is extremely fast and is competitive with cv-LASSO in most parameter regimes. Based on our
experimental/theoretical results there are some obvious directions to go in the future:
• Develop an efficient estimator that has theoretical guarantees for a more general design matrix, in
particular that satisfies the restricted eigenvalue condition.
• Find a choice of box size that is more robust to sparsity and SNR, which is still efficient to compute.

• Find an efficient estimator that has similar robustness to cross-validated LASSO but also scales
past p > 10000.

Although this estimator is by no means a replacement for existing estimators in typical regimes, it scales
extremely well into high dimensions and performs as well if not better when p ≫ n. This regime seems
the most interesting for developing more robust estimators.

Acknowledgments
We thank Abhinav Nellore for discussions on parameter selection in high dimensional problems which
motivated this work. R. Ward and C. Kennedy were partially supported during this work by NSF
CAREER grant #1255631.

References
[1] Richard Baraniuk, Mark Davenport, Ronald DeVore, and Michael Wakin. A simple proof of the
restricted isometry property for random matrices. Constructive Approximation, 28(3):253–263, 2008.
[2] Alexandre Belloni, Victor Chernozhukov, and Lie Wang. Square-root lasso: pivotal recovery of
sparse signals via conic programming. Biometrika, 98(4):791–806, 2011.

9

1.5

1.5

1

1

0.5

0.5

0

0
oracle

window

window-svd

moment

cv-lasso

1.5

1.5

1

1

0.5

0.5

0

oracle

window

window-svd

moment

cv-lasso

oracle

window

window-svd

moment

cv-lasso

0
oracle

window

window-svd

moment

cv-lasso

Figure 4: Orthogonal design matrix, p = 200. Left to right, top to bottom: Signal-less (kβk = 0), low
SNR (α = 0.1, kβk = 1), medium SNR (α = 0.1, kβk = 5), high SNR (α = 0.1, kβk = 10) respectively.
respectively.

10

1.6

1.4

1.2

1

0.8

0.6

0.4
oracle

window

window-svd

moment

cv-lasso

2.5
1.6
2
1.4
1.5

1.2

1

1

0.8

0.5

0.6
0
0.4
-0.5
oracle

window

window-svd

moment

cv-lasso

oracle

window

window-svd

moment

cv-lasso

Figure 5: Dense signal, p = 200. Left to right, top to bottom: Low SNR (α = 0.9, kβk = 1), medium SNR
(α = 0.9, kβk = 5), high SNR (α = 0.9, kβk = 10), respectively.

11

[3] Peter J Bickel, Ya’acov Ritov, and Alexandre B Tsybakov. Simultaneous analysis of lasso and
dantzig selector. The Annals of Statistics, pages 1705–1732, 2009.
[4] Emmanuel J Candes and Mark A Davenport. How well can we estimate a sparse vector? Applied
and Computational Harmonic Analysis, 34(2):317–323, 2013.
[5] Emmanuel J Candes and Terence Tao. Decoding by linear programming. IEEE transactions on
information theory, 51(12):4203–4215, 2005.
[6] Sourav Chatterjee and Jafar Jafarov. Prediction error of cross-validated lasso. arXiv preprint
arXiv:1502.06291, 2015.
[7] Lee H Dicker. Variance estimation in high-dimensional linear models. Biometrika, 101(2):269–284,
2014.
[8] Jianqing Fan, Shaojun Guo, and Ning Hao. Variance estimation using refitted cross-validation
in ultrahigh dimensional regression. Journal of the Royal Statistical Society: Series B (Statistical
Methodology), 74(1):37–65, 2012.
[9] Darren Homrighausen and Daniel McDonald. The lasso, persistence, and cross-validation. In International Conference on Machine Learning, pages 1031–1039, 2013.
[10] Jan-Christian Hütter and Philippe Rigollet. Optimal rates for total variation denoising. In Conference on Learning Theory, pages 1115–1146, 2016.
[11] Ron Kohavi. A study of cross-validation and bootstrap for accuracy estimation and model selection.
In Ijcai, volume 14, pages 1137–1145, 1995.
[12] Beatrice Laurent and Pascal Massart. Adaptive estimation of a quadratic functional by model
selection. Annals of Statistics, pages 1302–1338, 2000.
[13] Karim Lounici, Massimiliano Pontil, Sara Van De Geer, Alexandre B Tsybakov, et al. Oracle
inequalities and optimal inference under group sparsity. The Annals of Statistics, 39(4):2164–2204,
2011.
[14] Enno Mammen, Sara van de Geer, et al. Locally adaptive regression splines. The Annals of Statistics,
25(1):387–413, 1997.
[15] Nicolai Meinshausen and Bin Yu. Lasso-type recovery of sparse representations for high-dimensional
data. The Annals of Statistics, pages 246–270, 2009.
[16] Deanna Needell and Rachel Ward. Near-optimal compressed sensing guarantees for total variation
minimization. IEEE transactions on image processing, 22(10):3941–3949, 2013.
[17] Garvesh Raskutti, Martin J Wainwright, and Bin Yu. Minimax rates of estimation for highdimensional linear regression over ℓq -balls. IEEE transactions on information theory, 57(10):6976–
6994, 2011.
[18] Holger Rauhut. Compressive sensing and structured random matrices. Theoretical foundations and
numerical methods for sparse recovery, 9:1–92, 2010.
[19] Stephen Reid, Robert Tibshirani, and Jerome Friedman. A study of error variance estimation in
lasso regression. Statistica Sinica, pages 35–67, 2016.
[20] Mark Rudelson and Roman Vershynin. On sparse reconstruction from fourier and gaussian measurements. Communications on Pure and Applied Mathematics, 61(8):1025–1045, 2008.
[21] Leonid I Rudin, Stanley Osher, and Emad Fatemi. Nonlinear total variation based noise removal
algorithms. Physica D: nonlinear phenomena, 60(1-4):259–268, 1992.
[22] Robert Tibshirani. Regression shrinkage and selection via the lasso. Journal of the Royal Statistical
Society. Series B (Methodological), pages 267–288, 1996.
[23] Sara A Van de Geer. High-dimensional generalized linear models and the lasso. The Annals of
Statistics, pages 614–645, 2008.
[24] Sara A Van De Geer, Peter Bühlmann, et al. On the conditions used to prove oracle results for the
lasso. Electronic Journal of Statistics, 3:1360–1392, 2009.
[25] Roman Vershynin. Introduction to the non-asymptotic analysis of random matrices. arXiv preprint
arXiv:1011.3027, 2010.
[26] Nicolas Verzelen et al. Minimax risks for sparse regressions: Ultra-high dimensional phenomenons.
Electronic Journal of Statistics, 6:38–90, 2012.

12

[27] Martin J Wainwright. Information-theoretic limits on sparsity recovery in the high-dimensional and
noisy setting. IEEE Transactions on Information Theory, 55(12):5728–5741, 2009.
[28] Yu-Xiang Wang, James Sharpnack, Alex Smola, and Ryan Tibshirani. Trend filtering on graphs. In
Artificial Intelligence and Statistics, pages 1042–1050, 2015.
[29] Fei Ye and Cun-Hui Zhang. Rate minimaxity of the lasso and dantzig selector for the lq loss in lr
balls. Journal of Machine Learning Research, 11(Dec):3519–3540, 2010.
[30] Cun-Hui Zhang and Jian Huang. The sparsity and bias of the lasso selection in high-dimensional
linear regression. The Annals of Statistics, pages 1567–1594, 2008.
[31] Tong Zhang et al. Some sharp performance bounds for least squares regression with l1 regularization.
The Annals of Statistics, 37(5A):2109–2144, 2009.

A

Proof Ingredients

Proposition 8. (Lemma 1 in [12]) Suppose Z has a chi-squared distribution with d degrees of freedom.
Then,
√
√
P[d − 2 dt ≤ Z ≤ d + 2 dt + 2t] ≥ 1 − 2e−t
∀t ≥ 0.
(5)

Proposition 9. (Proposition 2.5 in [18]) Suppose Ωu ∩ Ωv = ∅, and that X ∈ Rn×p satisfies RIP of
order s0 and level δ > 0 with s0 = |Ωu | + |Ωv |. Then,
√
(6)
kXΩTu XΩv k2→2 ≤ δ

Proposition 10. (Equation (5.5) in [25]) Let X be a Gaussian random variable with mean 0, variance
σ. Then,
2
2
P[|X| > t] ≤ 2e−t /2σ ,
t ≥ 1.

B
B.1

Proofs
Proof of Theorem 3

Consider the window estimators
Sj =

1 X
|yi |2
L i∈Ω
j

1 X
|βi + ηi |2
=
L i∈Ω
j

1 X
1 X
1 X
=
|βi |2 +
|ηi |2 + 2
βi ηi .
L i∈Ω
L i∈Ω
L i∈Ω
j

j

j

P

Set Ej := L1 i∈Ωj |ηi |2 . Ej is a sum of L independent squares of N (0, σ 2 ) random variables. Then Ej
concentrates strongly around its expected value,
E(Ej ) = σ 2 .
Note that Ej has a chi-squared distribution with L degrees of freedom, so by (5) with the choice t =
log(p)2 and after a union bound over all p/L windows, we get that with probability at least 1 − p2 ,

1−

5
log(p)



σ 2 ≤ Ej ≤



1+

5
log(p)



σ2,

holds uniformly for all j ∈ {1, 2, . . . , p/L}, assuming that L ≥ log 3 (p).
p
by assumption, the pigeon hole principle implies that at least
Since L ≤ 2s
Ωβ . On any such “good” window k we have kβk:k+L−1 k22 = 0 and hence
|Sk − σ 2 | ≤

13

5σ 2
.
log(p)

p
2L

windows do not overlap

(7)

2

5σ
Thus, if S is the average over a subset of the good windows, then also |S − σ 2 | ≤ log(p)
.
Now,
to
bound
the
estimator
above
on
any
window,
we
need
some
control
on
the cross term
P
i∈Ωj βi ηi . Note that this quantity is just a sum of i.i.d. Gaussians with mean zero and with vari-

ance kβΩj ∩Ωβ k22 σ 2 ; thus, by concentration, we have that with probability at least 1 − 2/p, the following
holds uniformly over all windows:
p
X
2σkβΩj ∩Ωβ k2 log(p)
βi ηi ≤
.
(8)
L
i∈Ω
j

Hence, for any any window,
2 X
1
kβΩj ∩Ωβ k22 + Ej −
βi ηi
L
L i∈Ω ∩Ω
j
β
√


kβΩj ∩Ωβ k2 2σ log p
1
5
√
√
≥ kβΩj ∩Ωβ k22 + 1 −
σ2 −
L
log(p)
L
L


5
σ 2 log(p)
2
≥ 1−
The minimal value of a quadratic x2 + b − ax is b − a2 /4
σ −
log(p)
L


6
σ2,
≥ 1−
log(p)

Sj ≥

where the final inequality holds because log2 (p)) ≤ L.
c2 = 2L Pp/(2L) S . By construction, σ
c2 ≤ S, where S is
Now, consider the surrogate estimator σ
(j)
S
S
j=1
p
the average over any p/(2L) “good” windows. From the above analysis, we have that with probability
exceeding 1 − p4 ,

c2 = (1 +
Thus, for our final estimator, σ
S

c2 − σ 2 | ≤
|σ
S

1
c2 ,
)σ
log(p)

c2 − σ 2 | ≤
|σ

B.2



6
σ2.
log(p)

we have

6
7
+
log(p)
(log(p))2



σ2

Proof of Theorem 4

Recall that ỹ := Z T y ∈ Rp . Consider the window estimate
1 X
Sj =
|ỹi |2
L i∈Ω
j

1 X T 2
2 X T
1 X
|(Z T Xβ)i |2 +
|Zi η| +
(Z Xβ)i (Z T η)i
=
L i∈Ω
L i∈Ω
L i∈Ω
j

j

j

1 T
2 X T
1
T
XΩβ βk22 + kZΩ
ηk22 +
(Z Xβ)i (Z T η)i
= kZΩ
j
j
L
L
L i∈Ω

(9)

j

The first term is small if Ωj and Ωβ have disjoint support, since X has the RIP, the center term gets
close to its expectation σ 2 due to standard concentration inequalities, and the third term is also small
due to standard concentration inequalities. More concretely, if we assume that Sj is a “good” window,
meaning that Ωj and Ωβ have disjoint support, by equation (6)
1
δkβk22
kXΩTj XΩβ βk22 ≤
.
L
L
√
√
All of the diagonal entries of Σj are in the range [ 1 − δ, 1 + δ], hence by (10)
1 T
1+δ
kZΩj XΩβ βk22 ≤
kXΩTj XΩβ βk22
L
L
δ(1 + δ)kβk22
≤
L
2δkβk22
≤
L

14

(10)

(11)

For the center term, note that kZΩj ηk22 = kPL ηk22 where PL is projection onto the first L coordinates.
Next, we know that kPL ηk22 has a chi-squared distribution with L degrees of freedom, so by (5) with
t = log(p2 ),
h
p
i
2
P |kPL ηk22 − Lσ 2 | ≤ 2σ 2
L log(p) + log(p) ≥ 1 − 2 .
p

Hence by a union bound, with probability at least 1 − p2 , the following holds uniformly over all windows:

T
|kZΩ
ηk22 /L − σ 2 | = |kPL ηk22 /L − σ 2 |
j
r
log(p)
log(p)
2
≤ 2σ
+ 2σ 2
L
L
5σ 2
≤
(12)
log(p)
P
For the final term in 9, note that L2 i∈Ωj (Z T Xβ)i (Z T η)i is a Gaussian random variable with variance
2σkZΩj Xβ βk2 /L. Thus, by Proposition 10 and (11), the following holds uniformly over all windows with
probability at least 1 − p1 :

p
4σkZΩj Xβ βk2 log(p)
2 X T
(Z Xβ)i (Z T η)i ≤
L i∈Ω
L
j
p
√
8 δσkβk2 log(p)
≤
,
L

(13)

(14)

Thus, averaging over any set of p/2L “good” windows, using (11) (12) and (14) we have
√
√
5σ 2
8 δσkβk2 log p
2L X
2δkβk22
2
+
+
Sj − σ ≤
p j
L
log p
L
c2 =
with probability at least 1 − p4 . Thus, by construction, the estimator σ
S

√
√
5σ 2
8 δσkβk2 log p
2δkβk22
2
c
2
+
+
σS ≤ σ +
L
log p
L

2L
p

P

(15)

j

S(j) also satisfies

c2 cannot be too small. The inequalities (13) and (12)
It remains to show that the window estimator σ
S
hold uniformly over all windows, not just good windows; hence, for any window Sj ,
1
1 T
2 T
T
kZΩ
XΩβ βk22 + kZΩ
ηk22 − kZΩ
XΩβ βk2 kXΩTj ηk2
j
j
j
L
L
L
√
8σkXΩTj XΩβ βk2 log p
1
5σ 2
≥ kXΩTj XΩβ βk22 + σ 2 −
−
L
log(p)
L

Sj ≥

≥ σ2 −

4σ 2 log(p)
5σ 2
−
log(p)
L

The minimal value of a quadratic x2 + b − ax is b − a2 /4

Combining the bounds,
−

√
√
4σ 2 log(p)
2L X
2δkβk22
5σ 2
8 δσkβk2 log p
5σ 2
−
≤
S(j) − σ 2 ≤
+
+
log(p)
L
p j
L
log p
L

c2 = (1 +
For our final estimator σ
c2 − σ 2 | ≤
|σ



1+

1
log(p)

1
c2 ,
)σ
S
log(p)

we have




p
√
kβk2
6σ 2
1
2δ
+
+ max 4σ 2 log(p), 8 δσkβk2 log(p)
L
log(p)
L

15

