Total Variation Classes Beyond 1d: Minimax Rates, and the
Limitations of Linear Smoothers

arXiv:1605.08400v1 [math.ST] 26 May 2016

Veeranjaneyulu Sadhanala∗

Yu-Xiang Wang∗

Ryan J. Tibshirani

Carnegie Mellon University
( These authors contributed equally)
∗

Abstract
We consider the problem of estimating a function defined over n locations on a d-dimensional grid
(having all side lengths equal to n1/d ). When the function is constrained to have discrete total variation
bounded by Cn , we derive the minimax optimal (squared) `2 estimation error rate, parametrized by n and
Cn . Total variation denoising, also known as the fused lasso, is seen to be rate optimal. Several simpler
estimators exist, such as Laplacian smoothing and Laplacian eigenmaps. A natural question is: can these
simpler estimators perform just as well? We prove that these estimators, and more broadly all estimators
given by linear transformations of the input data, are suboptimal over the class of functions with bounded
variation. This extends fundamental findings of Donoho and Johnstone [1998] on 1-dimensional total
variation spaces to higher dimensions. The implication is that the computationally simpler methods cannot
be used for such sophisticated denoising tasks, without sacrificing statistical accuracy. We also derive
minimax rates for discrete Sobolev spaces over d-dimensional grids, which are, in some sense, smaller
than the total variation function spaces. Indeed, these are small enough spaces that linear estimators can
be optimal—and a few well-known ones are, such as Laplacian smoothing and Laplacian eigenmaps, as
we show. Lastly, we investigate the problem of adaptivity of the total variation denoiser to these smaller
Sobolev function spaces.

1

Introduction

Let G = (V, E) be a d-dimensional grid graph, i.e., a lattice graph, with equal side lengths. Label the nodes
as V = {1, . . . , n}, and edges as E = {e1 , . . . , em }. Consider data y = (y1 , . . . , yn ) ∈ Rn observed over
the nodes, from a model
yi ∼ N (θ0,i , σ 2 ), i.i.d., for i = 1, . . . , n,
(1)
where θ0 = (θ0,1 , . . . , θ0,n ) ∈ Rn is an unknown mean parameter to be estimated, and σ 2 > 0 is the marginal
noise variance. It is assumed that θ0 displays some kind of regularity over the grid G, e.g., θ0 ∈ Td (Cn ) for
some Cn > 0, where

Td (Cn ) = θ : kDθk1 ≤ Cn ,
(2)
and D ∈ Rm×n is the edge incidence matrix of G. This has `th row D` = (0, . . . , −1, . . . , 1, . . . , 0), with
a −1 in the ith location, and 1 in the jth location, provided that the `th edge is e` = (i, j) with i < j.
Equivalently, L = DT D is the graph Laplacian matrix of G, and thus
X
X
kDθk1 =
|θi − θj |,
and
kDθk22 = θT Lθ =
(θi − θj )2 .
(i,j)∈E

(i,j)∈E

We will refer to the class in (2) as a discrete total variation (TV) class, and to the quantity kDθ0 k1 as the
discrete total variation of θ0 , though for simplicity we will often drop the word “discrete”.
1

The problem of estimating θ0 given a total variation bound as in (2) is of great importance in both
nonparametric statistics and signal processing, and has many applications, e.g., changepoint detection for 1d
grids, and image denoising for 2d and 3d grids. There has been much methodological and computational
work devoted to this problem, resulting in practically efficient estimators in dimensions 1, 2, 3, and beyond.
However, theoretical performance, and in particularly optimality, is only really well-understood in the 1dimensional setting. This paper seeks to change that, and offers theory in d-dimensions that parallel more
classical results known in the 1-dimensional case.
Estimators under consideration. Central role to our work is the total variation (TV) denoising or fused
lasso estimator (e.g., Rudin et al. [1992], Vogel and Oman [1996], Tibshirani et al. [2005], Chambolle and
Darbon [2009], Hoefling [2010], Tibshirani and Taylor [2011], Sharpnack et al. [2012], Barbero and Sra
[2014]), defined by the convex optimization problem
θ̂TV = argmin ky − θk22 + λkDθk1 ,

(3)

θ∈Rn

where λ ≥ 0 is a tuning parameter. Another pair of methods that we study carefully are Laplacian smoothing
and Laplacian eigenmaps, which are most commonly seen in the context of clustering, dimensionality
reduction, and semi-supervised learning, but are also useful tools for estimation in a regression setting like
ours (e.g., Belkin and Niyogi [2002, 2003], Smola and Kondor [2003], Zhu et al. [2003], Belkin and Niyogi
[2004], Zhou et al. [2005], Belkin et al. [2005], Belkin and Niyogi [2005], Ando and Zhang [2006], Sharpnack
and Singh [2010]). The Laplacian smoothing estimator is given by
θ̂LS = argmin ky − θk22 + λkDθk22 ,

i.e., θ̂LS = (I + λL)−1 y,

(4)

θ∈Rn

for a tuning parameter λ ≥ 0, where in the second expression we have written θ̂LS in closed-form (this is
possible since it is the minimizer of a convex quadratic). For Laplacian eigenmaps, we must introduce the
eigendecomposition of the graph Laplacian, L = V ΣV T , where Σ = diag(ρ1 , . . . , ρn ) with 0 = ρ1 < ρ2 ≤
. . . ≤ ρn , and where V = [V1 , V2 , . . . , Vn ] ∈ Rn×n has orthonormal columns. The Laplacian eigenmaps
estimator is
T
θ̂LE = V[k] V[k]
y, where V[k] = [V1 , V2 , . . . , Vk ] ∈ Rn×k ,
(5)
where now k ∈ {1, . . . , n} acts as a tuning parameter.
Laplacian smoothing and Laplacian eigenmaps are appealing because they are (relatively) simple: they
are just linear transformations of the data y. Indeed, as we are considering G to be a grid, both estimators in
(4), (5) can be computed very quickly, in nearly O(n) time, since the columns of V here are discrete cosine
transform (DCT) basis vectors when d = 1, or Kronecker products thereof, when d ≥ 2 (e.g., Conte and
de Boor [1980], Godunov and Ryabenkii [1987], Kunsch [1994], Ng et al. [1999], Wang et al. [2008]). The
TV denoising estimator in (3), on the other hand, cannot be expressed in closed-form, and is much more
difficult to compute, especially when d ≥ 2, though several advances have been made over the years (see the
references above, and in particular Barbero and Sra [2014] for an efficient operator-splitting algorithm and
nice literature survey). Importantly, these computational difficulties are often worth it: TV denoising often
practically outperforms `2 -regularized estimators like Laplacian smoothing (and also Laplacian eigenmaps)
in image denoising tasks, as it is able to better preserve sharp edges and object boundaries (this is now widely
accepted, early references are, e.g., Acar and Vogel [1994], Dobson and Santosa [1996], Chambolle and
Lions [1997]). See Figure 1 for an example, using the often-studied “cameraman” image.
In the 1d setting, classical theory from nonparametric statistics draws a clear distinction between the
performance of TV denoising and estimators like Laplacian smoothing and Laplacian eigenmaps. Perhaps
surprisingly, this theory has not yet been fully developed in dimensions d ≥ 2. Arguably, the comparison
2

Noisy image

Laplacian smoothing

TV denoising

Figure 1: Comparison of Laplacian smoothing and TV denoising for the common“cameraman” image. TV denoising
provides a more visually appealing result, and also achieves about a 35% reduction in MSE compared to Laplacian
smoothing (MSE being measured to the original image). Both methods were tuned optimally.

between TV denoising and Laplacian smoothing and Laplacian eigenmaps is even more interesting in higher
dimensions, because the computational gap between the methods is even larger (the former method being
much more expensive, say in 2d and 3d, than the latter two). Shortly, we review the 1d theory, and what is
known in d-dimensions, for d ≥ 2. First, we introduce notation.
Notation. For deterministic (nonrandom) sequences an , bn we write an = O(bn ) to denote that an /bn is
−1
upper bounded for all n large enough, and an  bn to denote that both an = O(bn ) and a−1
n = O(bn ).
Also, for random sequences An , Bn , we write An = OP (Bn ) to denote that An /Bn is bounded in probability.
We abbreviate a ∧ b = min{a, b} and a ∨ b = max{a, b}. For an estimator θ̂ of the parameter θ0 in (1), we
define its mean squared error (MSE) to be
1
kθ̂ − θ0 k22 .
n
The risk of θ̂ is the expectation of its MSE, and for a set K ⊆ Rn , we define the minimax risk and minimax
linear risk to be




R(K) = inf sup E MSE(θ̂, θ0 ) and RL (K) = inf sup E MSE(θ̂, θ0 ) ,
MSE(θ̂, θ0 ) =

θ̂

θ0 ∈K

θ̂ linear θ0 ∈K

respectively, where the infimum on in the first expression is over all estimators θ̂, and in the second expression
over all linear estimators θ̂, meaning that θ̂ = Sy for a matrix S ∈ Rn×n . We will also refer to linear
estimators as linear smoothers. Note that both Laplacian smoothing in (4) and Laplacian eigenmaps in (5) are
linear smoothers, but TV denoising in (3) is not. Lastly, in somewhat of an abuse of nomenclature, we will
often call the parameter θ0 in (1) a function, and a set of possible values for θ0 as in (2) a function space; this
comes from thinking of the components of θ0 as the evaluations of an underlying function over n locations
on the grid. This embedding has no formal importance, but it is convenient notationally, and matches the
notation in nonparametric statistics.
Review: TV denoising in 1d. The classical nonparametric statistics literature [Donoho et al., 1990, Donoho
and Johnstone, 1998, Mammen and van de Geer, 1997] provides a more or less complete story for estimation
under total variation constraints in 1d. See also Tibshirani [2014] for a translation of these results to a setting
more consistent (notationally) to that in the current paper. Assume that d = 1 and Cn = C > 0, a constant
(not growing with n). The results in Donoho and Johnstone [1998] imply that
R(T1 (C))  n−2/3 .
3

(6)

Further, Mammen and van de Geer [1997] showed that the TV denoiser θ̂TV in (3), with λ  n1/3 , satisfies
MSE(θ̂TV , θ0 ) = OP (n−2/3 ),

(7)

for all θ0 ∈ T1 (C), and is thus minimax rate optimal over T1 (C). (In assessing rates here and throughout, we
do not distinguish between convergence in expectation versus convergence in probability.) Wavelet denoising,
under various choices of wavelet bases, also achieves the minimax rate. However, many simpler estimators
do not. To be more precise, it is shown in Donoho and Johnstone [1998] that
RL (T1 (C))  n−1/2 .

(8)

Therefore, a substantial number of commonly used nonparametric estimators—such as running mean
estimators, smoothing splines, kernel smoothing, Laplacian smoothing, and Laplacian eigenmaps, which are
all linear smoothers—have a major deficiency when it comes to estimating functions of bounded variation.
Roughly speaking, they will require many more samples to estimate θ0 within the same degree of accuracy as
an optimal method like TV or wavelet denoising (on the order of −1/2 times more samples to achieve an
MSE of ). Further theory and empirical examples (e.g., Donoho and Johnstone [1994a, 1998], Tibshirani
[2014]) offer the following perspective: linear smoothers cannot cope with functions in T (C) that have
spatially inhomogeneous smoothness, i.e., that vary smoothly at some locations and vary wildly at others.
Linear smoothers can only produce estimates that are smooth throughout, or wiggly throughout, but not a mix
of the two. They can hence perform well over smaller, more homogeneous function classes like Sobolev or
Holder classes, but not larger ones like total variation classes (or more generally, Besov and Triebel classes),
and for these, one must use more sophisticated, nonlinear techniques. A motivating question: does such a gap
persist in higher dimensions, between optimal nonlinear and linear estimators, and if so, how big is it?
Review: TV denoising in multiple dimensions. Recently, Wang et al. [2016] established rates for TV
denoising over various graph models, including grids, and Hutter and Rigollet [2016] made improvements,
particularly in the case of d-dimensional grids with d ≥ 2. We can combine Propositions 4 and 6 of Hutter
and Rigollet [2016] with Theorem 3 of Wang et al. [2016] to give the following result: if d ≥ 2, and Cn
is an arbitrary sequence (potentially unbounded with n), then the TV denoiser θ̂TV in (3) satisfies, over all
θ0 ∈ Td (Cn ),


 √

Cn log n
Cn log n
TV
TV
MSE(θ̂ , θ0 ) = OP
for d = 2, and MSE(θ̂ , θ0 ) = OP
for d ≥ 3, (9)
n
n
√
with λ  log n for d = 2, and λ  log n for d ≥ 3. Note that, at first glance, this is a very different result
from the 1d case. We expand on this next.

2

Summary of results

A gap in multiple dimensions. For estimation of θ0 in (1) when d ≥ 2, consider, e.g., the simplest possible
linear smoother: the mean estimator, θ̂mean = ȳ1 (where 1 = (1, . . . , 1) ∈ Rn , the vector of all 1s). Lemma
4, given below, implies that over θ0 ∈ Td (Cn ), the MSE of the mean estimator is bounded in probability by
Cn2 log n/n for d = 2, and Cn2 /n for d ≥ 3. Compare this to (9). When Cn = C > 0 is a constant, i.e., when
the TV of θ0 is assumed to be bounded (which is assumed for the 1d results in (6), (7), (8)), this means that
the TV denoiser and the mean estimator converge to θ0 at the same rate, basically (ignoring log terms), the
“parametric rate” of 1/n, for estimating a finite-dimensional parameter! That TV denoising and such a trivial
linear smoother perform comparably over 2d and 3d grids could not be farther from the story in 1d, where
TV denoising is separated by an unbridgeable gap from all linear smoothers, as shown in (6), (7), (8).
Our results in Section 3 clarify this conundrum, and can be summarized by three points.
4

• We argue in Section 3.1 that there is a proper “canonical” scaling for the TV class defined in (2). E.g.,
√
when d = 1, this yields Cn  1, a constant, but when d = 2, this yields Cn  n, and Cn also
diverges with n for all d ≥ 3. Sticking with d = 2 as an interesting example, we see that under such a
scaling, the MSE rates achieved by TV denoising and the mean estimator respectively, are drastically
different; ignoring log terms, these are
Cn2
 1,
(10)
n
√
respectively. Hence, TV denoising has an MSE rate of 1/ n, in a setting where the mean estimator
has a constant rate, i.e., a setting where it is not even known to be consistent.
Cn
1
√
n
n

and

• We show in Section 3.3 that our choice to study the mean estimator here is not somehow “unlucky” (it
is not a particularly bad linear smoother, nor is the upper bound on its MSE loose): the minimax linear
risk over Td (Cn ) is on the order Cn2 /n, for all d ≥ 2. Thus, even the best linear smoothers have the
same poor performance as the mean over Td (Cn ).
• We show in Section 3.2 that the TV estimator is (essentially) minimax optimal over Td (Cn ), as the
minimax risk over this class scales as Cn /n (ignoring log terms).

100

100

10-1

10-1

MSE

MSE

To summarize, these results reveal a significant gap between linear smoothers and optimal estimators like
TV denoising, for estimation over Td (Cn ) in d dimensions, with d ≥ 2, as long as Cn scales appropriately.
Roughly speaking, the TV classes encompass a challenging setting for estimation because they are very
broad, containing a wide array of functions—both globally smooth functions, said to have homogeneous
smoothness, and functions with vastly different levels of smoothness at different grid locations, said to have
heterogeneous smoothness. Linear smoothers cannot handle heterogeneous smoothness, and only nonlinear
methods can enjoy good estimation properties over the entirety of Td (Cn ). To reiterate, a telling example is
√
√
d = 2 with the canonical scaling Cn  n, where we see that TV denoising achieves the optimal 1/ n rate
√
(up to log factors), meanwhile, the best linear smoothers have max risk that is constant over T2 ( n). See
Figure 2 for an illustration.
√
Trivial scaling, Cn  1
Canonical scaling, Cn  n

10-2

10-3

10-4 2
10

10-3

TV denoising (-tted slope -0.88)
Laplacian smoothing (-tted slope -0.99)
Mean estimator (-tted slope -1.01)
Trivial rate: n!1

103

104

10-2

10-4 2
10

105

n

TV denoising (-tted slope -0.84)
Laplacian smoothing (-tted slope -0.01)
Mean estimator (-tted slope 0.00)
Minimax rate: n!1=2

103

104

105

n

√

Figure 2: MSE curves for estimation over a 2d grid, under two very different scalings of Cn : constant and n. The
parameter θ0 was a “one-hot” signal, with all but one component equal to 0. For each n, the results were averaged
over 5 repetitions, and Laplacian smoothing and TV denoising were tuned for optimal average MSE.

5

Minimax rates over smaller function spaces, and adaptivity. Sections 4 and 5 are focused on different
function spaces, discrete Sobolev spaces, which are `2 analogs of discrete TV spaces as we have defined them
in (2). Under the canonical scaling of Section 3.1, Sobolev spaces are contained in TV spaces, and the former
can be roughly thought of as containing functions of more homogeneous smoothness. The story now is more
optimistic for linear smoothers, and the following is a summary.
• In Section 4, we derive minimax rates for Sobolev spaces, and prove that linear smoothers—in
particular, Laplacian smoothing and Laplacian eigenmaps—are optimal over these spaces.
• In Section 5, we discuss an interesting phenomenon, a phase transition of sorts, at d = 3 dimensions.
When d = 1 or 2, the minimax rates for a TV space and its inscribed Sobolev space match; when d ≥ 3,
they do not, and the inscribed Sobolev space has a faster minimax rate. Aside from being an interesting
statement about the TV and Sobolev function spaces in high dimensions, this raises an important
question of adaptivity over the smaller Sobolev function spaces. As the minimax rates match for d = 1
and 2, any method optimal over TV spaces in these dimensions, such as TV denoising, is automatically
optimal over the inscribed Sobolev spaces. But the question remains open for d ≥ 3—does, e.g., TV
denoising adapt to the faster minimax rate over Sobolev spaces? We present empirical evidence to
suggest that this may be true, and leave a formal study to future work.
Other considerations and extensions. There are many problems related to the one that we study in this
paper. Clearly, minimax rates for the TV and Sobolev classes over general graphs, not just d-dimensional
grids, are of interest. Our minimax lower bounds for TV classes actually apply to generic graphs with
bounded max degree, though it is unclear whether to what extent they are sharp beyond grids; a detailed
study will be left to future work. Another related topic is that of higher-order smoothness classes, i.e., classes
contaning functions whose derivatives are (say) of bounded variation. The natural extension of TV denoising
here is called trend filtering, defined via the regularization of discrete higher-order derivatives. In the 1d
setting, minimax rates, the optimality of trend filtering, and the suboptimality of linear smoothers is already
well-understood [Tibshirani, 2014]. Trend filtering has been defined and studied to some extent on general
graphs [Wang et al., 2016], but no notions of optimality have been investigated beyond 1d. This will also be
left to future work. Lastly, it is worth mentioning that there are other estimators (i.e., other than the ones
we study in detail) that attain or nearly attain minimax rates over various classes we consider in this paper.
E.g., wavelet denoising is known to be optimal over TV classes in 1d [Donoho and Johnstone, 1998]; and
comparing recent upper bounds from Needell and Ward [2013], Hutter and Rigollet [2016] with the lower
bounds in this work, we see that wavelet denoising is also nearly minimax in 2d (ignoring log terms).

3
3.1

Analysis over TV classes
Canonical scaling for TV and Sobolev classes

We start by establishing what we call a “canonical” scaling for the radius Cn of the TV ball Td (Cn ) in (2), as
well as the radius Cn0 of the Sobolev ball Sd (Cn0 ), defined as

Sd (Cn0 ) = θ : kDθk2 ≤ Cn0 .
(11)
Proper scalings for Cn , Cn0 will be critical for properly interpreting our new results in d dimensions, in a way
that is comparable to known results for d = 1 (which are usually stated in terms of the 1d scalings Cn  1,
√
Cn0  1/ n). To study (2), (11), it helps to introduce a third function space,
n
o
Hd (1) = θ : θi = f (i1 /` . . . , id /`), i = 1, . . . , n, for some f ∈ Hdcont (1) .
(12)
6

Above, we have mapped each location i on the grid to a multi-index (i1 , . . . , id ) ∈ {1, . . . , `}d , where
` = n1/d , and Hdcont (1) denotes the (usual) continuous Holder space on [0, 1]d , i.e., functions that are 1Lipschitz with respect to the `∞ norm. We seek an embedding that is analogous to the embedding of
continuous Holder, Sobolev, and total variation spaces in 1d functional analysis, namely,
Hd (1) ⊆ Sd (Cn0 ) ⊆ Td (Cn ).

(13)

Our first lemma provides a choice of Cn , Cn0 that makes the above true. Its proof, as with all proofs in this
paper, can be found in the appendix.
Lemma 1. For d ≥ 1, the embedding in (13) holds with choices Cn  n1−1/d and Cn0  n1/2−1/d . Such
choices are called the canonical scaling for the function classes in (2), (11).
As a sanity check, both the (usual) continuous Holder and Sobolev function spaces in d dimensions are
known to have minimax risks that scale as n−2/(2+d) , in a standard nonparametric regression setup (e.g.,
Gyorfi et al. [2002]). Under the canonical scaling Cn0  n1/2−1/d , our results in Section 4 show that the
discrete Sobolev class Sd (n1/2−1/d ) also admits a minimax rate of n−2/(2+d) .

3.2

Minimax rates over TV classes

The following is a lower bound for the minimax risk of the TV class Td (Cn ) in (2).
Theorem 2. Assume n ≥ 2, and denote dmax = 2d. Then, for constants c > 0, ρ1 ∈ (2.34, 2.35),

p
√

σCn 1 + log(σdmax n/Cn )
√


if Cn ∈ [σdmax log n, σdmax n/ ρ1 ]

dmax n
√
R(Td (Cn )) ≥ c ·
.
2 /n
2 /(d2
n)
∨
σ
if
C
<
σd
C
log n

n
max
max
n


√
σ 2 /ρ
if Cn > σdmax n/ ρ1
1

(14)

The proof uses a simplifying reduction of the TV class, via Td (Cn ) ⊇ B1 (Cn /dmax ), the latter set
denoting the `1 ball of radius Cn /dmax in Rn . It then invokes a sharp characterization of the minimax risk in
normal means problems over `p balls due to Birge and Massart [2001]. Several remarks are in order.
Remark 1. The first line on the right-hand side in (14) often provides the most useful lower bound. To see this,
recall that under the canonical
scaling for TV classes, we have Cn = n1−1/d . For all d ≥ 2, this certainly
√
√
implies Cn ∈ [σdmax log n, σdmax n/ ρ1 ], for large n.
Remark 2. Even though its construction is very simple, the lower bound
risk in (14) is
√ on the minimax
√
sharp or nearly sharp
in
many
interesting
cases.
Assume
that
C
∈
[σd
log
n,
σd
n/
ρ
n
max
max
1 ]. The lower
p
bound rate is Cn log(n/Cn )/n. When d = 2, we see that this is very close to the upper bound rate of
Cn log n/n achieved by the TV denoiser, as stated in (9). These two differ by at most a log n factor (achieved
when
upper bound rate of
√ Cn  n). When d ≥ 3, we see that the lower bound rate is even closer to the √
Cn log n/n achieved by the TV denoiser, as in (9). These two now differ by at most a log n factor (again
achieved when Cn  n). We hence conclude that the TV denoiser is essentially minimax optimal in all
dimensions d ≥ 2.
√
Remark 3. When d = 1, and (say) Cn  1, the lower bound rate of log n/n given by Theorem 2 is not
sharp; we know from Donoho and Johnstone [1998] (recall (6)) that the minimax rate over T1 (1) is n−2/3 .
The result in the theorem (and also Theorem 3) in fact holds more generally, beyond grids: for an arbitrary
graph G, its edge incidence matrix D, and Td (Cn ) as defined in (2), the result holds for dmax equal to the
max degree of G. It is unclear to what extent this is sharp, for different graph models.

7

3.3

Minimax linear rates over TV classes

We now turn to a lower bound on the minimax linear risk of the TV class Td (Cn ) in (2).
Theorem 3. Recall the notation dmax = 2d. Then
σ 2 Cn2
σ2
1
RL (Td (Cn )) ≥ 2
∨
≥
2
2
Cn + σ dmax n
n
2

Cn2
d2max n

!
∧ σ2

∨

σ2
.
n

(15)

The proof relies on an elegant meta-theorem on minimax rates from Donoho et al. [1990], which uses
the concept of a “quadratically convex” set, whose minimax linear risk is the same as that of its hardest
rectangular subproblem. An alternative proof can be given entirely from first principles. See the appendix.
√
Remark 4. When Cn2 grows with n, but not too fast (scales as n, at most), the lower bound rate in (15) will
be Cn2 /n. Compared to the Cn /n minimax rate from Theorem 2 (ignoring log terms), we see a clear gap
between optimal nonlinear and linear estimators. In fact, under the canonical scaling Cn  n1−1/d , for any
d ≥ 2, this gap is seemingly huge: the lower bound for the minimax linear rate will be a constant, whereas
the minimax rate from Theorem 2 (ignoring log terms) will be n−1/d .
We now show that the lower bound in Theorem 3 is essentially tight, and remarkably, it is certified by
analyzing two trivial linear estimators: the mean estimator and the identity estimator.
Lemma 4. Let Mn denote the largest column norm of D† . For the mean estimator θ̂mean = ȳ1,

 σ 2 + Cn2 Mn2
,
E MSE(θ̂mean , θ0 ) ≤
n
θ0 ∈Td (Cn )
sup

√
From Proposition 4 in Hutter and Rigollet [2016], we have Mn = O( log n) when d = 2 and Mn = O(1)
when d ≥ 3.
The risk of the identity estimator θ̂id = y is clearly σ 2 . Combining this logic with Lemma 4 gives the
upper bound RL (Td (Cn )) ≤ (σ 2 + Cn2 Mn2 )/n ∧ σ 2 . Comparing this with the lower bound described in
Remark 4, we see that the two rates basically match, modulo the Mn2 factor in the upper bound, which only
provides an extra log n factor when d = 2. The takeaway message: in the sense of max risk, the best linear
smoother does not perform much better than the trivial estimators.
Additional empirical experiments, similar to those shown in Figure 2, are given in the appendix.

4

Analysis over Sobolev classes

Our first result here is a lower bound on the minimax risk of the Sobolev class Sd (Cn0 ) in (11).
Theorem 5. For a universal constant c > 0,
R(Sd (Cn0 )) ≥

 σ2
2
2d
c
(nσ 2 ) d+2 (Cn0 ) d+2 ∧ nσ 2 ∧ n2/d (Cn0 )2 + .
n
n

Elegant tools for minimax analysis from Donoho et al. [1990], which leverage the fact that the ellipsoid
Sd (Cn0 ) is orthosymmetric and quadratically convex (after a rotation), are used to prove the result.
The next theorem gives upper bounds, certifying that the above lower bound is tight, and showing that
Laplacian eigenmaps and Laplacian smoothing, both linear smoothers, are optimal over Sd (Cn0 ) for all d and
for d = 1, 2, or 3 respectively.

8

n−2/3

Dimension 2
√
n−1/2 log n

Dimension d ≥ 3
√
n−1/d log n

n−2/3

n−1/2

n− 2+d

Function class

Dimension 1

TV ball Td (n1−1/d )
Sobolev ball Sd (n1/2−1/d )

2

Table 1: Summary of rates for canonically-scaled TV and Sobolev spaces.
Linear signal in 2d

Linear signal in 3d

100

100
TV denoising (-tted slope -0.54)
Laplacian smoothing (-tted slope -0.62)
TV-ball minimax rate: n!1=2
Sobolev-ball minimax rate: n!1=2

MSE

10-1

MSE

10-1

10-2

10-2
TV denoising (-tted slope -0.44)
Laplacian smoothing (-tted slope -0.50)
TV-ball minimax rate: n!1=3
Sobolev-ball minimax rate: n!2=5

10-3 2
10

103

104

10-3 2
10

105

n

103

104

105

n

Figure 3: MSE curves for estimating a “linear” signal, a very smooth signal, over 2d and 3d grids. For each n, the
results were averaged over 5 repetitions, and Laplacian smoothing and TV denoising were tuned for best average MSE
performance. The signal was set to satisfy kDθ0 k2  n1/2−1/d , matching the canonical scaling.

Theorem 6. For Laplacian eigenmaps, θ̂LE in (5), with k  ((n(Cn0 )d )2/(d+2) ∨ 1) ∧ n, we have
 cσ 2


2
2d
c
sup E MSE(θ̂LE , θ0 ) ≤
(nσ 2 ) d+2 (Cn0 ) d+2 ∧ nσ 2 ∧ n2/d (Cn0 )2 +
,
n
n
0)
θ0 ∈Sd (Cn
for a universal constant c > 0, and n large enough. When d = 1, 2, or 3, the same bound holds for Laplacian
smoothing θ̂LS in (5), with λ  (n/(Cn0 )2 )2/(d+2) (and a possibly different constant c).
Remark 5. As shown in the proof, Laplacian smoothing is nearly minimax rate optimal over Sd (Cn0 ) when
d = 4, just incurring an extra log factor. It is unclear to us whether this method is still (nearly) optimal when
d ≥ 5; based on insights from our proof technique, we conjecture that it is not.

5

A phase transition, and adaptivity

The TV and Sobolev classes in (2) and (11), respectively, display a curious relationship. We reflect on
Theorems 2 and 5, under the canonical scaling with Cn  n1−1/d and Cn0  n1/2−1/d , that, recall, guarantees
Sd (Cn0 ) ⊆ Td (Cn ). (This is done for concreteness; statements could also be made outside of this case, subject
√
to an appropriate relationship with Cn /Cn0  n.) When d = 1, both the TV and Sobolev classes have a
minimax rate of n−2/3 (this TV result is actually due to Donoho and Johnstone [1998], as stated in (6), not
−1/2 , the
Theorem 2). When d = 2, both the TV and Sobolev
√ classes again have the same minimax rate of n
caveat being that the rate for TV class has an extra log n factor. But for all d ≥ 3, the rates for the canonical
TV and Sobolev classes differ, and the smaller Sobolev spaces have faster rates than their inscribing TV
spaces. This may be viewed as a phase transition at d = 3; see Table 1.
We may paraphrase to say that 2d is just like 1d, in that expanding the Sobolev ball into a larger TV ball
does not hurt the minimax rate, and methods like TV denoising are automatically adaptive, i.e., optimal over
9

both the bigger and smaller classes. However, as soon as we enter the 3d world, it is no longer clear whether
TV denoising can adapt to the smaller, inscribed Sobolev ball, whose minimax rate is faster, n−2/5 versus
n−1/3 (ignoring log factors). Theoretically, this is an interesting open problem that we do not approach in
this paper and leave to future work.
We do, however, investigate the matter empirically: see Figure 3, where we run Laplacian smoothing
and TV denoising on a highly smooth “linear” signal θ0 . This is constructed so that each component θi is
proportional to i1 + i2 + . . . + id (using the multi-index notation (i1 , . . . , id ) of (12) for grid location i), and
the Sobolev norm is kDθ0 k2  n1/2−1/d . Arguably, these are among the “hardest” types of functions for TV
denoising to handle. The left panel, in 2d, is a case in which we know that TV denoising attains the minimax
rate; the right panel, in 3d, is a case in which we do not, though empirically, TV denoising surely seems to be
doing better than the slower minimax rate of n−1/3 (ignoring log terms) that is associated with the larger TV
ball.
Even if TV denoising is shown to be minimax optimal over the inscribed Sobolev balls when d ≥ 3, note
that this does not necessarily mean that we should scrap Laplacian smoothing in favor of TV denoising, in
all problems. Laplacian smoothing is the unique Bayes estimator in a normal means model under a certain
Markov random field prior (e.g., Sharpnack and Singh [2010]); statistical decision theory therefore tells that
it is admissible, i.e., no other estimator—TV denoising included—can uniformly dominate it.

6

Discussion

We conclude with a quote from Albert Einstein: “Everything should be made as simple as possible, but no
simpler”. In characterizing the minimax rates for TV classes, defined over d-dimensional grids, we have
shown that simple methods like Laplacian smoothing and Laplacian eigenmaps—or even in fact, all linear
estimators—must be passed up in favor of more sophisticated, nonlinear estimators, like TV denoising, if one
wants to attain the optimal max risk. Such a result was previously known when d = 1; our work has extended
it to all dimensions d ≥ 2. We also characterized the minimax rates over discrete Sobolev classes, revealing
an interesting phase transition where the optimal rates over TV and Sobolev spaces, suitably scaled, match
when d = 1 and 2 but diverge for d ≥ 3. It is an open question as to whether an estimator like TV denoising
can be optimal over both spaces, for all d.

A

Proofs

We present proofs of all results, according to the order in which they appear in the paper.

A.1

Proof of Lemma 1 (canonical scaling)

Suppose that θ ∈ Hd (1) that is a discretization of a 1-Lipshitz function f , i.e., θi = f (i1 /` . . . , id /`),
i = 1, . . . , n. We first we compute and bound its squared Sobolev norm
X
X
2
kDθk22 =
(θi − θj )2 =
f (i1 /`, . . . , id /`) − f (j1 /`, . . . , jd /`)
(i,j)∈E

(i,j)∈E

≤

X

(i1 /`, . . . , id /`) − (j1 /`, . . . , jd /`)

2
∞

(i,j)∈E

= m/`2 ,
where, recall, we denote by m = |E| the number of edges in the grid. In the second line we used the
1-Lipschitz property of f , and in the third we used that multi-indices corresponding to adjacent locations
10

√
on the grid are exactly 1 apart, in `∞ distance. Thus we see that setting Cn0 = m/` gives the desired
containment Sd (Cn0 ) ⊇ Hd (1). It is always true that m  n for a d-dimensional grid (though the constant
√
may depend on d), so that m/`  n1/2−1/d . This completes the proof for the Sobolev class scaling.
√
As for TV class scaling, the result follows from the simple fact that kxk1 ≤ mkxk2 for any x ∈ Rm ,
√
so that we may take Cn = mCn0 = n1−1/d .

A.2

Proof of Theorem 2 (minimax rates over TV classes)

Here and henceforth, we use the notation Bp (r) = {x : kxkp ≤ r} for the `p ball of radius r, where p, r > 0
(and the ambient dimension will be determined based on the context).
We begin with a very simple lemma, that embeds an `1 ball inside the TV ball Td (Cn ).
Lemma 7. Let G be a graph with maximum degree dmax , and let D ∈ Rm×n be its incidence matrix. Then
for any r > 0, it holds that B1 (r/dmax ) ⊆ Td (r).
Proof. Write Di for the ith column of D. The proof follows from the observation that, for any θ,
n
n


X
X
kDθk1 =
D i θi ≤
kDi k1 |θi | ≤
max kDi k1 kθk1 = dmax kθk1 .
i=1

1

i=1,...,n

i=1

To prove Theorem 2, we will rely on a result from Birge and Massart [2001], which gives a lower bound
for the risk in a normal means problem, over `p balls. Another related, earlier result is that of Donoho and
Johnstone [1994b]; however, the Birge and Massart result places no restrictions on the radius of the ball in
question, whereas the Donoho and Johnstone result does. Translated into our notation, the Birge and Massart
result is as follows.
Lemma 8 (Proposition 5 of Birge and Massart [2001]). Assume i.i.d. observations yi ∼ N (θ0,i , σ 2 ), i =
1, . . . , n, and n ≥ 2. Then the minimax risk over the `p ball Bp (rn ), where 0 < p < 2, satisfies


 p 1−p/2
√

σ n
√
2−p
p

σ rn 1 + log
if σ log n ≤ rn ≤ σn1/p / ρp
p

rn
√
.
n · R(Bp (rn )) ≥ c ·
2
r
if rn < σ log n

n


√
 2
σ n/ρp
if rn > σn1/p / ρ
Here c > 0 is a universal constant, and ρp > 1.76 is the unique solution of ρp log ρp = 2/p.
Finally, applying Lemma 8 to B1 (Cn /dmax ) almost gives the lower bound as stated in Theorem 2.
However, note that the minimax risk in question is trivially lower bounded by σ 2 /n, because
1
inf
sup
Ekθ̂ − θ0 k22 ≥ inf
θ̂ θ0 ∈Td (Cn ) n
θ̂

n

sup
θ0 :θ0,1 =...=θ0,n

1X
E(θ̂i − θ0,1 )2
n
i=1
2

= inf sup E(θ̂1 − θ0,1 )
θ̂1

=

θ0,1

σ2

.
n
In the second to last line, the problem is to estimate a 1-dimensional mean parameter θ0,1 , given the
observations yi ∼ N (θ0,1 , σ 2 ), i.i.d., for i = 1, . . . , n; this has a well-known minimax risk of σ 2 /n. What
this means for our TV problem: to derive a lower bound for the minimax rate over Td (Cn ), we may take the
maximum of the result of applying Lemma 8 to B1 (Cn /dmax ) and σ 2 /n. √
One can see that the term σ 2 /n
only plays a role for small Cn , i.e., it effects the case when Cn < σdmax log n, where the lower bound
becomes Cn2 /(d2max n) ∨ σ 2 /n.
11

A.3

Proof of Theorem 3 (minimax linear rates over TV classes)

First we recall a few definitions, from Donoho et al. [1990]. Given a set A ⊆ Rk , its quadratically convex
hull qconv(A) is defined as

qconv(A) = (x1 , . . . , xk ) : (x21 , . . . , x2k ) ∈ conv(A2+ ) , where

A2+ = (a21 , . . . , a2k ) : a ∈ A, ai ≥ 0, i = 1, . . . , k .
(Here conv(B) denotes the convex hull of a set B.) Furthermore, the set A is called quadratically convex
provided that qconv(A) = A. Also, A is called orthosymmetric provided that (a1 , . . . , ak ) ∈ A implies
(σ1 a1 , . . . , σk ak ) ∈ A, for any choice of signs σ1 , . . . , σk ∈ {−1, 1}.
Now we proceed with the proof. Following from equation (7.2) of Donoho et al. [1990],

qconv B1 (Cn /dmax ) = B2 (Cn /dmax ).
Theorem 11 of Donoho et al. [1990] states that, for orthosymmetric, compact sets, such as B1 (Cn /dmax ),
the minimax linear risk equals that of its quadratically convex hull. Moreover, Theorem 7 of Donoho et al.
[1990] tells us that for sets that are orthosymmetric, compact, convex, and quadratically convex, such as
B2 (Cn /dmax ), the minimax linear risk is the same as the minimax linear risk over the worst rectangular
√
√
subproblem. We consider B∞ (Cn /(dmax n)), and abbreviate rn = Cn /(dmax n). It is fruitful to study
rectangles because the problem separates across dimensions, as in
inf

sup

θ̂ linear θ0 ∈B∞ (rn )


 X
n
n 
1X
1
2
(θ̂i − θ0,i ) =
inf
E
n
n
θ̂i linear
i=1

i=1

= inf

sup

sup

2



E(θ̂i − θ0,i )

|θ0,i |≤rn

E(θ̂1 − θ0,1 )2 .

θ̂1 linear |θ0,1 |≤rn

Thus it suffices to compute the minimax linear risk over the 1d class {θ0,1 : |θ0,1 | ≤ rn }. It is easily shown
(e.g., see Section 2 of Donoho et al. [1990]) that this is rn2 σ 2 /(rn2 + σ22 ), and so this is precisely the minimax
linear risk for B2 (Cn /dmax ), and for B1 (Cn /dmax ).
To get the first lower bound as stated in the theorem, we simply take a maximum of rn2 σ 2 /(rn2 + σ22 ) and
σ 2 /n, as the latter is the minimax risk for estimating a 1-dimensional mean parameter given n observations
in a normal model with variance σ 2 , recall the end of the proof of Theorem 2. To get the second, we use the
fact that 2ab/(a + b) ≥ min{a, b}. This completes the proof.

A.4

Alternative proof of Theorem 3

Here, we reprove Theorem 3 using elementary arguments. We write y = θ0 + , for  ∼ N (0, σ 2 I). Given
an arbitary linear estimator, θ̂ = Sy for a matrix S ∈ Rn×n , observe that

 1
1
E MSE(θ̂, θ0 ) = Ekθ̂ − θ0 k22 = EkS(θ0 + ) − θ0 k22
n
n
1
1
= EkSk22 + k(S − I)θ0 k22
n
n
σ2
1
=
kSk2F + k(S − I)θ0 k22 ,
n
n

12

(16)

which we may view as the variance and (squared) bias terms, respectively. Now denote by ei the ith standard
basis vector, and consider
!


2
2
σ2
1
σ
C
max k(I − S)ei k22
kSk2F +
sup
k(S − I)θ0 k22 ≥
kSk2F + 2 n
n
n
dmax n i=1,...,n
θ0 :kDθ0 k1 ≤Cn n
n
σ2
Cn2 X
2
≥
kSkF + 2
k(I − S)ei k22
n
dmax n2

=
≥
=

i=1
2
C
kSk2F + 2 n 2 k(I − S)k2F
n
dmax n
n
n
σ2 X 2
C2 X
Sii + 2 n 2
(1 − Sii )2
n
dmax n
i=1
i=1

n 
2
X
Cn
1
2
2 2
(1 − Sii ) .
σ Sii + 2
n
dmax n
i=1

σ2

Here Sii , i = 1, . . . , n denote the diagonal entries of S. To bound each term in the sum, we apply the simple
inequality ax2 + b(1 − x)2 ≥ ab/(a + b) for all x (since a short calculation shows that the quadratic in x
here is minimized at x = b/(a + b)). We may continue on lower bounding the last displayed expression,
giving


σ2
1
σ 2 Cn2
2
2
kSkF +
sup
k(S − I)θ0 k2 ≥ 2
.
n
Cn + σ 2 d2max n
θ0 :kDθ0 k1 ≤Cn n
Lastly, we may take the maximum of this with σ 2 /n in order to derive a final lower bound, as argued in the
proof of Theorem 3.

A.5

Proof of Lemma 4 (mean estimator over TV classes)

For this estimator, the smoother matrix is S = 11T /n and so kSk2F = 1. From (16), we have

 σ2
1
+ kθ0 − θ̄0 1k22 ,
E MSE(θ̂mean , θ0 ) =
n
n
where θ̄0 = (1/n)

Pn

i=1 θ0,i .

Now
sup

θ0 :kDθ0 k1 ≤Cn

1
1
kθ0 − θ̄0 1k22 =
sup
kxk22
n
x∈row(D):kDxk1 ≤Cn n
=

sup
z∈col(D):kzk1 ≤Cn

≤

sup
z:kzk1 ≤Cn

1
kD† zk22
n

1
kD† zk22
n

Cn2
max kDi† k22
n i=1,...,n
C2M 2
≤ n n,
n
=

which establishes the desired bound.

13

A.6

Proof of Theorem 5 (minimax rates over Sobolev classes)

Recall that we denote by L = V ΣV T the eigendecomposition of the graph Laplacian L = DT D, where
Σ = diag(ρ1 , . . . , ρn ) with 0 = ρ1 < ρ2 ≤ . . . ≤ ρn , and where V ∈ Rn×n has orthonormal columns.
Also denote by D = U Σ1/2 V T the singular value decomposition of the edge incidence matrix D, where
U ∈ Rm×n has orthonormal columns.1 First notice that
kDθ0 k2 = kU Σ1/2 V T θ0 k2 = kΣ1/2 V T θ0 k2 .
This suggests that a rotation by V T will further simplify the minimax risk over Sd (Cn0 ), i.e.,
inf

sup

θ̂

0
θ0 :kΣ1/2 V T θ0 k2 ≤Cn

1
Ekθ̂ − θ0 k22 = inf
n
θ̂

sup
0
θ0 :kΣ1/2 V T θ0 k2 ≤Cn

= inf
γ̂

1
EkV T θ̂ − V T θ0 k22
n

1
Ekγ̂ − γ0 k22 ,
n

sup
0
γ0 :kΣ1/2 γ0 k2 ≤Cn

(17)

where we have rotated and now consider the new parameter γ0 = V T θ0 , constrained to lie in
Ed (Cn0 )


=

γ:

n
X

ρi γi2

≤

(Cn0 )2


.

i=2

To be clear, in the rotated setting (17) we observe a vector y 0 = V T y ∼ N (γ0 , σ 2 I), and the goal is to
estimate the mean parameter γ0 . Since there are no constraints along the first dimension, we can separate
out the MSE in (17) into that incurred on the first component, and all other components. Decomposing
γ0 = (α0 , β0 ) ∈ R1×(n−1) , with similar notation for an estimator γ̂,
inf
γ̂

sup
0)
γ0 ∈Ed (Cn

1
1
Ekγ̂ − γ0 k22 = inf sup E(α̂ − α0 )2 + inf
α̂
n
n
α0
β̂
=

σ2
+ inf
n
β̂

sup
0 ))
β0 ∈P−1 (Ed (Cn

1
Ekβ̂ − β0 k22
n

1
Ekβ̂ − β0 k22 ,
0 )) n
β0 ∈P−1 (Ed (Cn
sup

(18)

where P−1 projects onto all coordinate axes but the 1st, i.e., P−1 (x) = (0, x2 , . . . , xn ), and in the second line
we have used the fact that the minimax risk for estimating a 1-dimensional parameter α0 given an observation
z ∼ N (α0 , σ 2 ) is simply σ 2 .
Let us lower bound the second term in (18), i.e., R(P−1 (Ed (Cn0 ))). The ellipsoid P−1 (Ed (Cn0 )) is
orthosymmetric, compact, convex, and quadratically convex, hence Theorem 7 in Donoho et al. [1990] tells
us that its minimax linear risk is the minimax linear risk of its hardest rectangular subproblem. Further,
Lemma 6 in Donoho et al. [1990] then tells us the minimax linear risk of its hardest rectangular subproblem
is, up to a constant factor, the same as the minimax (nonlinear) risk of the full problem. More precisely,
Lemma 6 and Theorem 7 from Donoho et al. [1990] imply
4
R(P−1 (Ed (Cn0 ))) ≥ RL (P−1 (Ed (Cn0 ))) =
sup
RL (H),
5
0 ))
H⊆P−1 (Ed (Cn

(19)

where the supremum above is taken over all rectangular subproblems, i.e., all rectangles H contained in
P−1 (Ed (Cn0 )).
1

When d = 1, we have m = n − 1 edges, and so it is not be possible for U to have orthonormal columns; however, we can just
take its first column to be all 0s, and take the rest as the eigenbasis for Rn−1 , and all the arguments given here will go through.

14

To study rectangular subproblems, it helps to reintroduce the multi-index notation for a location i on the
d-dimensional grid, writing this as (i1 , . . . , id ) ∈ {1, . . . , `}d , where ` = n1/d . For a parameter 2 ≤ τ ≤ `,
we consider rectangular subsets of the form2

H(τ ) = β ∈ Rn−1 : |βi | ≤ ti (τ ), i = 2, . . . , n , where
(
P
Cn0 /( j1 ,...,jd ≤τ ρj1 ,...,jd )1/2 if i1 , . . . , id ≤ τ
ti (τ ) =
, for i = 2, . . . , n.
0
otherwise
P
It is not hard to check that H(τ ) ⊆ {β ∈ Rn−1 : ni=2 ρi βi2 ≤ (Cn0 )2 } = P−1 (Ed (Cn0 )). Then, from (19),
n

R(P−1 (Ed (Cn0 )))

1 X ti (τ )2 σ 2
≥ sup RL (H(τ )) = sup
n
ti (τ )2 + σ 2
τ
τ
i=1

= sup
τ

(τ d − 1)σ 2 (Cn0 )2
1
P
.
n (Cn0 )2 + j1 ,...,jd ≤τ ρj1 ,...,jd

The first equality is due to the fact that the minimax risk for rectangles decouples across dimensions, and
the 1d minimax linear risk is straightforward to compute for an interval, as argued in the proof Theorem 3;
the second equality simply comes from a short calculation following the definition of ti (τ ), i = 2, . . . , n.
Applying Lemma 9, on the eigenvalues of the graph Laplacian matrix L for a d-dimensional grid, we have
that for a constant c > 0,
(τ d − 1)σ 2 (Cn0 )2
(τ d − 1)σ 2 (Cn0 )2
1
σ 2 (Cn0 )2
P
≥
≥
.
(Cn0 )2 + j1 ,...,jd ≤τ ρj1 ,...,jd
2 (Cn0 )2 τ −d + cσ 2 τ 2 /`2
(Cn0 )2 + cσ 2 τ d+2 /`2
We can choose τ to maximize the expression on the right above, given by
∗



τ =

`2 (Cn0 )2
cσ 2



1
d+2

.

When 2 ≤ τ ∗ ≤ `, this provides us with the lower bound on the minimax risk
R(P−1 (Ed (Cn0 ))) ≥ RL (H(τ ∗ )) ≥

2d
2
c1
1
τ d σ 2 (Cn0 )2
(nσ 2 ) d+2 (Cn0 ) d+2 ,
d
4
2d =
2n 2(cσ 2 ) d+2 (C 0 ) d+2 `− d+2
n

(20)

n

for a constant c1 > 0. When τ ∗ < 2, we can use τ = 2 as lower bound on the minimax risk,
R(P−1 (Ed (Cn0 ))) ≥ RL (H(2)) ≥

1
σ 2 `2 (Cn0 )2
c2
≥ `2 (Cn0 )2 ,
2
0
2
−d
2
2
2n ` (Cn ) 2 + cσ 2
n

(21)

for a constant c2 > 0, where in the last inequality, we used the fact that `2 (Cn0 )2 ≤ cσ 2 2d+2 (just a constant)
since we are in the case τ ∗ < 2. Finally, when τ ∗ > `, we can use τ = ` as a lower bound on the minimax
risk,
1
σ 2 (Cn0 )2
R(P−1 (Ed (Cn0 ))) ≥ RL (H(`)) ≥
≥ c3 σ 2 ,
(22)
2n `−d (Cn0 )2 + cσ 2
for a constant c3 > 0, where in the last inequality, we used that cσ 2 ≤ `−d (Cn0 )2 as we are in the case τ ∗ > `.
Taking a minimum of the lower bounds in (20), (21), (22), as a way to navigate the cases, gives us a final
lower bound on R(P−1 (Ed (Cn0 ))), and completes the proof.
Here, albeit unconvential, it helps to index β ∈ H(τ ) ⊆ Rn−1 according to components i = 2, . . . , n, rather than i =
1, . . . , n − 1. This is so that we may keep the index variable i to be in correspondence with positions on the grid.
2

15

A.7

Proof of Theorem 6 (Laplacian eigenmaps and Laplacian smoothing over Sobolev classes)

We will prove the results for Laplacian eigenmaps and Laplacian separately.
T , for a tuning parameter
Laplacian eigenmaps. The smoother matrix for this estimator is Sk = V[k] V[k]
k = 1, . . . , n. From (16),


 σ2
1
E MSE(θ̂LE , θ0 ) =
k + k(I − Sk )θ0 k22 .
n
n
Now we write k = τ d , and analyze the max risk of the second term,
sup
0
θ0 :kDθ0 k2 ≤Cn

1
1
k(I − Sk )θ0 k22 = sup
k(I − Sk )D† zk22
n
0 n
z:kzk2 ≤Cn

(Cn0 )2 2
σmax (I − Sk )D†
n
(Cn0 )2
1
≤
2
n 4 sin (πτ /(2`))
(C 0 )2 4`2
≤ n
.
n π2τ 2
=

Here we denote by σmax (A) the maximum singular value of a matrix A. The last inequality above used the
simple lower bound sin(x) ≥ x/2 for x ∈ [0, π/2]. The earlier inequality used that


T
(I − Sk )D† = (I − V[k] V[k]
)V T (Σ† )1/2 U T = 0, . . . , 0, Vk+1 , . . . , Vn (Σ† )1/2 U T ,
where we have kept the same notation for the singular value decomposition of D as in the proof of Theorem
2 ((I − S )D † ) is the reciprocal of the (k + 1)st smallest eigenvalue ρ
5. Therefore σmax
k
k+1 of the graph
Laplacian L. For any subset A of the eigenvalues Λ of the Laplacian, with |A| = k, ρk+1 ≥ min Λ \ A. So
for a d-dimensional grid,
ρk+1 ≥ min Λ \ {ρi1 ,··· ,id : 1 ≤ ij ≤ τ for all 1 ≤ j ≤ d}
= 4 sin2 (πτ /(2`))
where recall ` = n1/d , as explained by (23), in the proof of Lemma 9.
Hence, we have established

 σ 2 σ 2 d (Cn0 )2 4`2
+ τ +
E MSE(θ̂LE , θ0 ) ≤
.
n
n
n π2τ 2
0
θ0 :kDθ0 k2 ≤Cn
sup

2

Choosing τ to balance the two terms on the right-hand side above results in τ ∗ = (2`Cn0 /(πσ)) d+2 . Plugging
in this choice of τ , while utilizing the bounds 1 ≤ τ ≤ `, very similar to the arguments given at the end of
the proof of Theorem 5, gives the result for Laplacian eigenmaps.
Laplacian smoothing. The smoother matrix for this estimator is Sλ = (I + λL)−1 , for a tuning parameter
λ ≥ 0. From (16),
n

 σ2 X
1
1
LS
E MSE(θ̂ , θ0 ) =
+ k(I − Sλ )θ0 k22 .
2
n
(1 + λρi )
n
i=1

16

When d = 1, 2, or 3, the first term upper is bounded by c1 σ 2 /n + c2 σ 2 /λd/2 , for some constants c1 , c2 > 0,
by Lemma 10. As for the second term,
sup
0
θ0 :kDθ0 k2 ≤Cn

1
k(I − Sλ )θ0 k22 = sup k(I − Sλ )D† zk22
n
0
z:kzk2 ≤Cn

(Cn0 )2 2
σmax (I − Sλ )D†
n

2
1
1
(Cn0 )2
1−
=
max
n i=2,...,n
1 + λρi ρi
0
2
(C )
λρi
= n λ max
i=2,...,n (1 + λρi )2
n
(C 0 )2 λ
≤ n
.
4n

=

In the third equality we have used the fact the eigenvectors of I − Sλ are the left singular vectors of D† ,
and in the last inequailty we have used the simple upper bound f (x) = x/(1 + x)2 ≤ 1/4 for x ≥ 0 (this
function being maximized at x = 1).
Therefore, from what we have shown,

 c1 σ 2 c2 σ 2 (Cn0 )2 λ
+ d/2 +
.
E MSE(θ̂LS , θ0 ) ≤
n
4n
0
λ
θ0 :kDθ0 k2 ≤Cn
sup

Choosing λ to balance the two terms on the right-hand side above gives λ∗ = c(n/(Cn0 )2 )2/(d+2) , for a
constant c > 0. Plugging in this choice, and using upper bounds from the trivial cases λ = 0 and λ = ∞
when Cn0 is very small or very large, respectively, gives the result for Laplacian smoothing.
P
Remark 6. When d = 4, Lemma 10 gives a slightly worse upper bound on ni=1 1/(1 + λρi )2 , with an
“extra” term (nc2 /λd/2 )) log(1 + c3 λ), for constants c2 , c3 > 0. It is not hard to show, by tracing through the
same arguments as given above that we can use this to establish an upper bound on the max risk of
 cσ 2


2
2d
c
(nσ 2 ) d+2 (Cn0 ) d+2 log(n/(Cn0 )2 ) ∧ nσ 2 ∧ n2/d (Cn0 )2 +
,
sup E MSE(θ̂LE , θ0 ) ≤
n
n
0)
θ0 ∈Sd (Cn
only slightly worse than the minimax optimal rate, by a log factor.
When d ≥ 5, our analysis provides a much worse bound for the max risk of Laplacian smoothing, as the
integral denoted I(d) in the proof of Lemma 10 grows very large when d ≥ 5. We conjecture that this not due
to slack in our proof technique, but rather, to the Laplacian smoothing estimator itself, since all inequalities
the proof can be fairly tight.

B

Utility lemmas used in the proofs of Theorems 5 and 6

This section contains some calculations on the partial sums of eigenvalues of the Laplacian matrix L, for
d-dimensions grids. These are useful for the proofs of both Theorem 5 and Theorem 6.
Lemma 9. Let L ∈ Rn×n denote the graph Laplacian matrix of a d-dimensional grid graph, and ρi1 ,...,id ,
(i1 , . . . , id ) ∈ {1, . . . , `}d be its eigenvalues, where ` = n1/d . Then there exists a constant c > 0 (dependent
on d) such that, for any 1 ≤ τ ≤ `,
X

ρi1 ,··· ,id ≤ c

(i1 ,...,id )∈{1,...,τ }d

17

τ d+2
.
`2

Proof. The eigenvalues of L can be written explicitly as
ρi = 4 sin2

 π(i − 1) 
 π(i − 1) 
1
d
+ . . . + 4 sin2
,
2`
2`

(i1 , . . . , id ) ∈ {1, . . . , `}d .

(23)

This follows from known facts about the eigenvalues for the Laplacian matrix of a 1d grid, and the fact
that the Laplacin matrix for higher-dimensional grids can be expressed in terms of a Kronecker sum of the
Laplacian matrix of an appropriate 1d grid (e.g., Conte and de Boor [1980], Kunsch [1994], Ng et al. [1999],
Wang et al. [2008, 2016], Hutter and Rigollet [2016]). We now use the fact that sin(x) ≤ x for all x ≥ 0,
which gives us the upper bound
X

ρi1 ,··· ,id ≤

(i1 ,...,id )∈{1,...,τ }d

≤
≤
=

π2
`2



X

(i1 − 1)2 + . . . + (id − 1)2



(i1 ,...,id )∈{1,...,τ }d
τ
π 2 d d−1 X
τ
(i − 1)2
`2
i=1
π 2 d d−1 3
τ
τ
`2

π 2 d d+2
τ
,
`2

as desired.
Lemma 10. Let L ∈ Rn×n denote the graph Laplacian matrix of a d-dimensional grid graph, and ρi ,
i = 1, . . . , n be its eigenvalues. Let λ ≥ 0 be arbitrary. For d = 1, 2, or 3, there are constants c1 , c2 > 0
such that
n
X
1
n
≤ c1 + c2 d/2 .
2
(1 + λρi )
λ
i=1

For d = 4, there are constants c1 , c2 , c3 > 0 such that
n
X
i=1


1
n 
≤
c
+
c
1
+
log(1
+
c
λ)
.
1
2 d/2
3
(1 + λρi )2
λ

Proof. We will use the explicit form of the eigenvalues as given in the proof of Lemma 9. In the expressions
below, we use c > 0 to denote a constant whose value may change from line to line. Using the inequality
sin x ≥ x/2 for x ∈ [0, π/2],
n
X
i=1

1
≤
(1 + λρi )2

X
(i1 ,...,id )∈{1,...,`}d

1+

Z
≤1+
[0,`]d 1
Z ` √d

=1+c
0

=1+c

n

Z

λd/2

|0

+

2
λ π4

π2
λ 4`
2

1
Pd

j=1 (ij

1
Pd


2 2 2
j=1 xj /`

1

2 r
2
1 + λ π4 r2 /`2
√
π
λd
d−1
2

18

d−1

u
du .
(1 + u2 )2
{z
}
I(d)

− 1)2
dx

dr

2

In the second inequality, we used the fact that the right-endpoint Riemann sum is always an underestimate for
the integral of a function that is monotone nonincreasing in each coordinate. In the third, we made a change
to spherical coordinates, and suppressed all of the angular variables, as they contribute at most a constant
factor. It remains to compute I(d), which can be done by symbolic integration:
√
 √  1 π
π d
1
−1 π
I(1) =
+
tan
λd ≤ + ,

2
2
2
4 4
4 1 + π λd
4

1
1
1
I(2) = −
,
≤

2
2 2 1 + π λd
2
4
 √  π
1
−1 π
I(3) = tan
λd ≤ , and
2
2
4


1
π2 
1
1
1
π2  1
I(4) = log 1 + λd +
−
≤
log
1
+
λd + .

2
2
4
2
2
4
2
2 1 + π λd
4

This completes the proof.

C

Additional experiments comparing TV denoising and Laplacian smoothing for piecewise constant functions
Piecewise constant signal in 2d

Piecewise constant signal in 3d
100

10-1

10-1

MSE

MSE

100

10-2

10-2

10-3

TV denoising (-tted slope -0.68)
Laplacian smoothing (-tted slope -0.36)
Minimax rate: n!1=2
Minimax linear rate: constant

102

103

104

10-3
102

105

TV denoising (-tted slope -0.57)
Laplacian smoothing (-tted slope -0.27)
Minimax rate: n!1=3
Minimax linear rate: constant

103

104

105

n

n

Figure 4: MSE curves for estimating a “piecewise constant” signal, having a single elevated region, over 2d and
3d grids. For each n, the results were averaged over 5 repetitions, and the Laplacian smoothing and TV denoising
estimators were tuned for best average MSE performance. We set θ0 to satisfy kDθ0 k1  n1−1/d , matching the
canonical scaling. Note that all estimators achieve better performance than that dictated by their minimax rates.

References
Robert Acar and Curtis R. Vogel. Analysis of total variation penalty methods. Inverse Problems, 10:
1217–1229, 1994.
Rie Ando and Tong Zhang. Learning on graph with Laplacian regularization. Advances in Neural Information
Processing Systems, 9, 2006.

19

Alvero Barbero and Suvrit Sra. Modular proximal optimization for multidimensional total-variation regularization. arXiv: 1411.0589, 2014.
Mikhail Belkin and Partha Niyogi. Using manifold structure for partially labelled classification. Advances in
Neural Information Processing Systems, 15, 2002.
Mikhail Belkin and Partha Niyogi. Laplacian eigenmaps for dimensionality reduction and data representation.
Neural Computation, 15(6):1373–1396, 2003.
Mikhail Belkin and Partha Niyogi. Semi-supervised learning on Riemannian manifolds. Machine Learning,
56(1–3):209–239, 2004.
Mikhail Belkin and Partha Niyogi. Towards a theoretical foundation for Laplacian-based manifold methods.
Conference on Learning Theory (COLT-05), 18, 2005.
Mikhail Belkin, Partha Niyogi, and Vikas Sindhwani. On manifold regularization. Proceedings of the
International Conference on Artificial Intelligence and Statistics, 8, 2005.
Lucien Birge and Pascal Massart. Gaussian model selection. Journal of the European Mathematical Society,
3(3):203–268, 2001.
Antonin Chambolle and Jerome Darbon. On total variation minimization and surface evolution using
parametric maximum flows. International Journal of Computer Vision, 84:288–307, 2009.
Antonin Chambolle and Pierre-Louis Lions. Image recovery via total variation minimization and related
problems. Numerische Mathematik, 76(2):167–188, 1997.
Samuel Conte and Carl de Boor. Elementary Numerical Analysis: An Algorithmic Approach. McGraw-Hill,
New York, 1980. International Series in Pure and Applied Mathematics.
David Dobson and Fadil Santosa. Recovery of blocky images from noisy and blurred data. SIAM Journal on
Applied Mathematics, 56(4):1181–1198, 1996.
David Donoho and Iain Johnstone. Ideal spatial adaptation by wavelet shrinkage. Biometrika, 81(3):425–455,
1994a.
David Donoho and Iain Johnstone. Minimax risk over `p -balls for `q -error. Probability Theory and Related
Fields, 99(2):277–303, 1994b.
David Donoho and Iain Johnstone. Minimax estimation via wavelet shrinkage. Annals of Statistics, 26(8):
879–921, 1998.
David Donoho, Richard Liu, and Brenda MacGibbon. Minimax risk over hyperrectangles, and implications.
Annals of Statistics, 18(3):1416–1437, 1990.
S. Godunov and V. Ryabenkii. Difference Schemes: An Introduction to the Underlying Theory. Elsevier,
Amsterdam, 1987. Number 19 in Studies in Mathematics and Its Applications.
Laszlo Gyorfi, Michael Kohler, Adam Krzyzak, and Harro Walk. A Distribution-Free Theory of Nonparametric Regression. Springer, New York, 2002.
Holger Hoefling. A path algorithm for the fused lasso signal approximator. Journal of Computational and
Graphical Statistics, 19(4):984–1006, 2010.

20

Jan-Christian Hutter and Philippe Rigollet. Optimal rates for total variation denoising. In Conference on
Learning Theory (COLT-16), 2016. to appear.
Hans Kunsch. Robust priors for smoothing and image restoration. Annals of the Institute of Statistical
Mathematics, 46(1):1–19, 1994.
Enno Mammen and Sara van de Geer. Locally apadtive regression splines. Annals of Statistics, 25(1):
387–413, 1997.
Deanna Needell and Rachel Ward. Stable image reconstruction using total variation minimization. SIAM
Journal on Imaging Sciences, 6(2):1035–1058, 2013.
Michael Ng, Raymond Chan, and Wun-Cheung Tang. A fast algorithm for deblurring models with Neumann
boundary conditions. SIAM Journal on Scientific Computing, 21(3):851–866, 1999.
Leonid Rudin, Stanley Osher, and Emad Faterni. Nonlinear total variation based noise removal algorithms.
Physica D: Nonlinear Phenomena, 60:259–268, 1992.
James Sharpnack and Aarti Singh. Identifying graph-structured activation patterns in networks. Advances in
Neural Information Processing Systems, 13, 2010.
James Sharpnack, Alessandro Rinaldo, and Aarti Singh. Sparsistency of the edge lasso over graphs. Proceedings of the International Conference on Artificial Intelligence and Statistics, 15:1028–1036, 2012.
Alexander Smola and Risi Kondor. Kernels and regularization on graphs. Proceedings of the Annual
Conference on Learning Theory, 16, 2003.
Robert Tibshirani, Michael Saunders, Saharon Rosset, Ji Zhu, and Keith Knight. Sparsity and smoothness via
the fused lasso. Journal of the Royal Statistical Society: Series B, 67(1):91–108, 2005.
Ryan J. Tibshirani. Adaptive piecewise polynomial estimation via trend filtering. Annals of Statistics, 42(1):
285–323, 2014.
Ryan J. Tibshirani and Jonathan Taylor. The solution path of the generalized lasso. Annals of Statistics, 39
(3):1335–1371, 2011.
Curtis R. Vogel and M. E. Oman. Iterative methods for total variation denoising. SIAM Journal on Scientific
Computing, 17(1):227–238, 1996.
Yilun Wang, Junfeng Yang, Wotao Yin, and Yin Zhang. A new alternating minimization algorithm for total
variation image reconstruction. SIAM Journal on Imaging Sciences, 1(3):248–272, 2008.
Yu-Xiang Wang, James Sharpnack, Alex Smola, and Ryan J. Tibshirani. Trend filtering on graphs. Journal
of Machine Learning Research, 2016. To appear.
Dengyong Zhou, Jiayuan Huang, and Bernhard Scholkopf. Learning from labeled and unlabeled data on a
directed graph. Proceedings of the International Conference on Machine Learning, 22, 2005.
Xiaojin Zhu, Zoubin Ghahramani, and John Lafferty. Semi-supervised learning using Gaussian fields and
harmonic functions. International Conference on Machine Learning (ICML-03), 20, 2003.

21

