Nonparametric Shape-restricted Regression

arXiv:1709.05707v1 [] 17 Sep 2017

Adityanand Guntuboyina∗ and Bodhisattva Sen†
Department of Statistics
University of California at Berkeley
423 Evans Hall
Berkeley, CA 94720
e-mail: aditya@stat.berkeley.edu
Department of Statistics
Columbia University
1255 Amsterdam Avenue
New York, NY 10027
e-mail: bodhi@stat.columbia.edu
Abstract: We consider the problem of nonparametric regression under shape constraints. The main examples include isotonic regression (with respect to any partial order), unimodal/convex regression, additive shape-restricted regression, and constrained
single index model. We review some of the theoretical properties of the least squares
estimator (LSE) in these problems, emphasizing on the adaptive nature of the LSE. In
particular, we study the risk behavior of the LSE, and its pointwise limiting distribution theory, with special emphasis to isotonic regression. We survey various methods
for constructing pointwise confidence intervals around these shape-restricted functions.
We also briefly discuss the computation of the LSE and indicate some open research
problems and future directions.
Keywords and phrases: adaptive risk bounds, bootstrap, Chernoff’s distribution,
convex regression, isotonic regression, likelihood ratio test, monotone function, order
preserving function estimation, projection on a closed convex set, tangent cone.

1. Introduction
In nonparametric shape-restricted regression the observations {(xi , yi ) : i = 1, . . . , n} satisfy
yi = f (xi ) + εi ,

for i = 1, . . . , n,

(1)

where x1 , . . . , xn are design points in some space (e.g., Rd , d ≥ 1), ε1 , . . . , εn are unobserved
mean-zero errors, and the real-valued regression function f is unknown but obeys certain
known qualitative restrictions like monotonicity, convexity, etc. Let F denote the class
of all such regression functions. Letting θ∗ := (f (x1 ), . . . , f (xn )), Y := (y1 , . . . , yn ) and
ε := (ε1 , . . . , εn ), model (1) may be rewritten as
Y = θ∗ + ε,

(2)

and the problem is to estimate θ∗ and/or f from Y , subject to the constraints imposed by
the properties of F. The constraints on the function class F translate to constraints on θ∗
of the form θ∗ ∈ C, where
C := {(f (x1 ), . . . , f (xn )) ∈ Rn : f ∈ F}
∗
†

(3)

Supported by NSF CAREER Grant DMS-16-54589.
Supported by NSF Grants DMS-17-12822 and AST-16-14743.
1
imsart-generic ver. 2014/10/16 file: ShapeRegArxiv.tex date: September 19, 2017

A. Guntuboyina and B. Sen/Shape-restricted Regression

2

is a subset of Rn (in fact, in most cases, C will be a closed convex set). In the following we
give some examples of shape-restricted regression.
Example 1.1 (Isotonic regression). Probably the most studied shape-restricted regression
problem is that of estimating a monotone (nondecreasing) regression function f when x1 <
. . . < xn are univariate design points. In this case, F is the class of all nondecreasing
functions on the interval [x1 , xn ], and the constraint set C reduces to
I := {(θ1 , . . . , θn ) ∈ Rn : θ1 ≤ . . . ≤ θn },

(4)

which is a closed convex cone in Rn (I is defined through n − 1 linear constraints). The
above problem is typically known as isotonic regression and has a long history in statistics;
see e.g., [20, 4, 124].
Example 1.2 (Order preserving regression on a partially ordered set). Isotonic regression
can be easily extended to the setup where the covariates take values in a space X with a
partial order -1 ; see e.g., [106, Chapter 1]. A function f : X → R is said to be isotonic (or
order preserving) with respect to the partial order - if for every pair u, v ∈ X ,
u - v ⇒ f (u) ≤ f (v).
For example, suppose that the predictors take values in R2 and the partial order - is defined
as (u1 , u2 ) - (v1 , v2 ) if and only if u1 ≤ v1 and u2 ≤ v2 . This partial order leads to a
natural extension of isotonic regression to two-dimensions; see e.g., [67, 105, 26]. One can
also consider other partial orders; see e.g., [117, 118] and the references therein for isotonic
regression with different partial orders. We will introduce and study yet another partial order
in Section 6.
Given data from model (1), the goal is to estimate the unknown regression function
f : X → R under the assumption that f is order preserving (with respect to the partial order
-). The restrictions imposed by the partial order - constrain θ∗ to lie in a closed convex
cone C which may be expressed as
{(θ1 , . . . , θn ) ∈ Rn : θi ≤ θj for every i, j such that xi - xj }.
Example 1.3 (Convex regression). Suppose that the underlying regression function f :
Rd → R (d ≥ 1) is known to be convex, i.e., for every u, v ∈ Rd ,
f (αu + (1 − α)v) ≤ αf (u) + (1 − α)f (v),

for every α ∈ (0, 1).

(5)

Convexity appears naturally in many applications; see e.g., [69, 79, 35] and the references
therein. The convexity of f constrains θ∗ to lie in a (polyhedral) convex set C ⊂ Rn which,
when d = 1 and the xi ’s are ordered, reduces to


θn − θn−1
n θ2 − θ 1
K := (θ1 , . . . , θn ) ∈ R :
≤ ... ≤
,
(6)
x2 − x1
xn − xn−1
whereas for d ≥ 2 the characterization of C is more complex; see e.g., [109].
1

A partial order is a binary relation - that is reflexive (x - x for all x ∈ X ), transitive (u, v, w ∈ X , u - v
and v - w imply u - w), and antisymmetric (u, v ∈ X , u - v and v - u imply u = v).
imsart-generic ver. 2014/10/16 file: ShapeRegArxiv.tex date: September 19, 2017

A. Guntuboyina and B. Sen/Shape-restricted Regression

3

Observe that when d = 1, convexity is characterized by nondecreasing derivatives (subgradients). This observation can be used to generalize convexity to k-monotonicity (k ≥ 1):
a real-valued function f is said to be k-monotone if its (k − 1)’th derivative is monotone;
see e.g., [85, 25]. For equi-spaced design points in R, this restriction constrains θ∗ to lie in
the set
{θ ∈ Rn : ∇k θ ≥ 0}
where ∇ : Rn → Rn is given by ∇(θ) := (θ2 − θ1 , θ3 − θ2 , . . . , θn − θn−1 , 0) and ∇k represents
the k-times composition of ∇. Note that the case k = 1 and k = 2 correspond to isotonic
and convex regression, respectively.
Example 1.4 (Unimodal regression). In many applications f , the underlying regression
function, is known to be unimodal; see e.g., [45, 24] and the references therein. Let Im , 1 ≤
m ≤ n, denote the convex set of all unimodal vectors (first decreasing and then increasing)
with mode at position m, i.e.,
Im := {(θ1 , . . . , θn ) ∈ Rn : θ1 ≥ . . . ≥ θm ≤ θm+1 ≤ . . . ≤ θn }.
Then, the unimodality of f constrains θ∗ to belong to U := ∪nm=1 Im . Observe that now U is
not a convex set, but a union of n convex cones.
Example 1.5 (Shape-restricted additive model). In an additive regression model one assumes that f : Rd → R (d ≥ 1) depends on each of the predictor variables in an additive
fashion, i.e., for (u1 , . . . , ud ) ∈ Rd ,
f (u1 , . . . , ud ) =

d
X

fi (ui ),

i=1

where fi ’s are one-dimensional functions and fi captures the influence of the i’th variable.
Observe that the additive model generalizes (multiple) linear regression. If we assume that
each of the fi ’s are shape-constrained, then one obtains a shape-restricted additive model; see
e.g., [5, 88, 95, 32] for a study of some possible applications, identifiability and estimation
in such a model.
Example 1.6 (Shape-restricted single index model). In a single index regression model one
assumes that the regression function f : Rd → R takes the form
f (x) = m(x> β ∗ ),

for all x ∈ Rd ,

where m : R → R and β ∗ ∈ Rd are unknown. Single index models are popular in many
application areas, including econometrics and biostatistics (see e.g., [103, 83]), as they
circumvent the curse of dimensionality encountered in estimating the fully nonparametric regression function by assuming that the link function depends on x only through a
one-dimensional projection, i.e., x> β ∗ . Moreover, the coefficient vector β ∗ provides interpretability. Observe that single index models extend generalized linear models (where the link
function m is assumed known). Moreover, as most known link functions are nondecreasing,
the monotone single index model (where m is assumed unknown but nondecreasing) arises
naturally in applications; see e.g., [98, 53, 8].
imsart-generic ver. 2014/10/16 file: ShapeRegArxiv.tex date: September 19, 2017

A. Guntuboyina and B. Sen/Shape-restricted Regression

4

Observe that all the aforementioned problems fall under the general area of nonparametric regression. However, it turns out that in each of the above problems one can use
classical techniques like least squares and/or maximum likelihood (without additional explicit regularization/penalization) to readily obtain tuning parameter-free estimators that
have attractive theoretical and computational properties. This makes shape-restricted regression different from usual nonparametric regression, where likelihood based methods are
generally infeasible. In this paper we try to showcase some of these attractive features of
shape-restricted regression and give an overview of the major theoretical advances in this
area.
Let us now introduce the estimator of θ∗ (and f ) that we will study in this paper.
The least squares estimator (LSE) θ̂ of θ∗ in shape-restricted regression is defined as the
projection of Y onto the set C (see (3)), i.e.,
θ̂ := arg min kY − θk2 ,
θ∈C

(7)

where k · k denotes the usual Euclidean norm in Rn . If C is a closed convex set then θ̂ ∈ C
is unique and is characterized by the following condition:
hY − θ̂, θ − θ̂i ≤ 0,

for all θ ∈ C,

(8)

where h·, ·i denotes the usual inner product in Rn ; see [17, Proposition 2.2.1]. It is easy to
see now that the LSE θ̂ is tuning parameter-free, unlike most nonparametric estimators.
However it is not generally easy to find a closed-form expression for θ̂. As for estimating f ,
any fˆn ∈ F that agrees with θ̂ at the data points xi ’s will be considered as a LSE of f .
In this paper we mainly review the main theoretical properties of the LSE θ̂ with special
emphasis on its adaptive nature. The risk behavior of θ̂ (in estimating θ∗ ) is studied in
Sections 2 and 3 — Section 2 mainly deals with the isotonic LSE in detail whereas Section 3
summarizes the main results for other shape-restricted problems. In Section 4 we study the
pointwise asymptotic limiting behavior of the LSE fˆn , in the case of isotonic and convex
regression, focusing on methods for constructing (pointwise) confidence intervals around f .
In the process of this review we highlight the main ideas and techniques used in the proofs
of the theoretical results; in fact, we give (nearly) complete proofs in many cases.
The computation of the LSE θ̂, in the various problems outlined above, is discussed
in Section 5. In Section 6 we mention a few open research problems and possible future
directions. Some of the detailed proofs are relegated to Section 7. Although the paper
mostly summarizes known results, we also present some new results — Theorems 2.1, 2.2
and 6.1, and Lemma 3.1 are new.
There are indeed many other important applications and examples of shape-restricted
regression beyond those highlighted so far. We briefly mention some of these below. Shape
constrained functions also arise naturally in interval censoring problems (e.g., in the current
status model; see [61, 71]), in survival analysis (e.g., in estimation of monotone/unimodal
hazard rates [70]), and in regression models where the response, conditional on the covariate,
comes from a regular parametric family (e.g., monotone response models [11]). It also arises
in the study of many inverse problems; e.g., deconvolution problems (see e.g., [61, 72]) and
the classical Wicksell’s corpuscle problem (see e.g., [54, 111]).
imsart-generic ver. 2014/10/16 file: ShapeRegArxiv.tex date: September 19, 2017

A. Guntuboyina and B. Sen/Shape-restricted Regression

5

In this paper we will mostly focus on estimation of the underlying shape-restricted function using the method of least squares. Although this produces tuning parameter-free estimators, the obtained LSEs are not “smooth”. There is also a line of research that combines
shape constraints with smoothness assumptions — see e.g., [97, 84, 60] (and the references therein) where kernel-based methods have been combined with shape-restrictions,
and see [87, 94, 98, 77] where splines are used in conjunction with the shape constraints.
2. Risk bounds for the isotonic LSE
In this section we attempt to answer the following question: “How good is θ̂ as an estimator
of θ∗ ?”. To quantify the accuracy of θ̂ we first need to fix a loss function. Arguably the
most natural loss function here is the squared error loss: kθ̂ − θ∗ k2 /n. As the loss function
is random, we follow the usual approach and study its expectation:
n

R(θ̂, θ∗ ) :=

X

 1
2
1
Eθ∗ kθ̂ − θ∗ k2 = Eθ∗
θ̂i − θi∗
n
n

(9)

i=1

which we shall refer to as the risk of the LSE θ̂. We focus on the risk in this paper. It
may be noted that upper bounds derived for the risk usually hold on the loss kθ̂ − θ∗ k2 /n
as well, with high probability. When ε ∼ Nn (0, σ 2 In ), this is essentially because kθ̂ − θ∗ k
concentrates around its mean; see [122] and [16] for more details on high probability results.
One can also try to study the risk under more general `p -loss functions. For p ≥ 1, let
n

R

(p)

X

 1
1
(θ̂, θ ) := Eθ∗ kθ̂ − θ∗ kpp = Eθ∗
θ̂i − θi∗
n
n
∗

p

(10)

i=1

1/p
P
n
p
, for u = (u1 , . . . , un ) ∈ Rn . We shall mostly focus on the
where kukp :=
j=1 |uj |
risk for p = 2 in this paper but we shall also discuss some results for p 6= 2.
In this section, we focus on the problem of isotonic regression (Example 1.1) and describe
bounds on the risk of the isotonic LSE. As mentioned in the Introduction, isotonic regression
is the most studied problem in shape-restricted regression where the risk behavior of the
LSE is well-understood. We shall present the main results here. The results described in this
section will serve as benchmarks to which risk bounds for other shape-restricted regression
problems (see Section 3) can be compared.
Throughout this section, θ̂ will denote the isotonic LSE (which is the minimizer of
kY −θk2 subject to the constraint that θ lies in the closed convex cone I described in (4)) and
θ∗ will usually denote an arbitrary vector in I (in some situations we deal with misspecified
risks where θ∗ is an arbitrary vector in Rn not necessarily in I).
The risk, R(θ̂, θ∗ ), essentially has two different kinds of behavior. As long as θ∗ ∈ I and
V (θ∗ ) := θn∗ − θ1∗ (referred to as the variation of θ∗ ) is bounded from above independently
of n, the risk R(θ̂, θ∗ ) is bounded from above by a constant multiple of n−2/3 . We shall refer
to this n−2/3 bound as the worst case risk bound mainly because it is, in some sense, the
maximum possible rate at which R(θ̂, θ∗ ) converges to zero. On the other hand, if θ∗ ∈ I
is piecewise constant with not too many constant pieces, then the risk R(θ̂, θ∗ ) is bounded
from above by the parametric rate 1/n up to a logarithmic multiplicative factor. This rate
imsart-generic ver. 2014/10/16 file: ShapeRegArxiv.tex date: September 19, 2017

A. Guntuboyina and B. Sen/Shape-restricted Regression

6

is obviously much faster compared to the worst case rate of n−2/3 which means that the
isotonic LSE is estimating piecewise constant nondecreasing sequences at a much faster rate.
In other words, the isotonic LSE is adapting to piecewise constant nondecreasing sequences
with not too many constant pieces. We shall therefore refer to this log n/n risk bound as
the adaptive risk bound.
The worst case risk bounds for the isotonic LSE will be explored in Section 2.1 while the
adaptive risk bounds are treated in Section 2.2. Proofs will be provided in Section 7. Before
proceeding to risk bounds, let us first describe some basic properties of the isotonic LSE.
An important fact about the isotonic LSE is that θ̂ = (θ̂1 , . . . , θ̂n ) can be explicitly
represented as (see [106, Chapter 1]):
Pv
j=u yj
θ̂j = min max
, for j = 1, . . . , n.
(11)
v≥j u≤j v − u + 1
This is often referred to as the min-max formula for isotonic regression. The isotonic LSE
is, in some sense, unique among shape-restricted LSEs because it has the above explicit
characterization. It is this characterization that allows for a precise study of the properties
of θ̂.
The above characterization of the isotonic LSE shows that θ̂ is piecewise constant, and
in each “block” (i.e., region of constancy) it is the average of the response values (within
the block); see [106, Chapter 1]. However, the blocks, their lengths and their positions, are
chosen adaptively by the algorithm, the least squares procedure. If θi∗ = f (xi ) for some
design points 0 ≤ x1 < · · · < xn ≤ 1, then we can define the isotonic LSE of f as the
piecewise constant function fˆn : [0, 1] → R which has jumps only at the design points and
such that fˆn (xi ) = θ̂i for each i = 1, . . . , n. Figure 1 shows three different scatter plots,
for three different regression functions f , with the fitted isotonic LSEs fˆn . Observe that for
the leftmost plot the block-sizes (of the isotonic LSE) vary considerably with the change
in slope of the underlying function f — the isotonic LSE, fˆn , is nearly constant in the
interval [0.3, 0.7] where f is relatively flat whereas fˆn has many small blocks towards the
boundary of the covariate domain where f has large slope. This highlights the adaptive
nature of the isotonic LSE fˆn and also provides some intuition as to why the isotonic LSE
adapts to piecewise constant nondecreasing functions with not too many constant pieces.
Moreover, in some sense, fˆn can be thought of as a kernel estimator (with the box kernel)
or a ‘regressogram’ ([120]), but with a varying bandwidth/window.
2.1. Worst case risk bound
The worst case risk bound for the isotonic LSE is given by the following inequality. Under
the assumption that the errors ε1 , . . . , εn are i.i.d. with mean zero and variance σ 2 , the risk
of the isotonic LSE satisfies the bound (see [129]):
 2
2/3
σ V (θ∗ )
σ 2 log(en)
R(θ̂, θ∗ ) ≤ C
+C
(12)
n
n
where V (θ∗ ) = θn∗ − θ1∗ denotes the variation of θ∗ ∈ I and C > 0 is a universal constant.
Let us try to understand each of the terms on the right side of (12). As long as the
variation V (θ∗ ) is not small the risk R(θ̂, θ∗ ) is given by (σ 2 V (θ∗ )/n)2/3 , up to a constant
imsart-generic ver. 2014/10/16 file: ShapeRegArxiv.tex date: September 19, 2017

A. Guntuboyina and B. Sen/Shape-restricted Regression

0.0

0.2

0.4

0.6

0.8

1.0

0.8
0.0

0.2

0.4

0.6

0.8
0.6
0.4
0.2
0.0

0.0

0.2

0.4

0.6

0.8

1.0

Piecewise constant

1.0

Constant

1.0

Cubic

7

0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.2

0.4

0.6

0.8

1.0

Fig 1: Plots of Y (circles), θ̂ (red), and θ∗ (blue) for three different f ’s: (i) cubic polynomial
(left plot), (ii) constant (middle plot), and (iii) piecewise constant. Here n = 60, and
ε ∼ Nn (0, σ 2 In ) with σ = 0.1. Here In denotes the identity matrix of order n.

multiplicative factor. This shows that the rate of estimating any monotone function (under
the `2 -loss) is n−2/3 . Moreover, (12) gives the explicit dependence of the risk on the variation
of θ∗ (and on σ 2 ).
The second term on the right side of (12) is also interesting — when V (θ∗ ) = 0, i.e., θ∗
is a constant sequence, (12) shows that the risk of the isotonic LSE scales like log n/n. This
is a consequence of the fact that θ̂ chooses its blocks (of constancy) adaptively depending
on the data. When θ∗ is the constant sequence, θ̂ has fewer blocks (in fact, it has of the
order of log n blocks; see e.g., [92, Theorem 1]) and some of the blocks will be very large (see
e.g., the middle plot of Figure 1), so that averaging the responses
within the large blocks
P
would yield a value very close to the grand mean Ȳ = ( ni=1 yi )/n (which has risk σ 2 /n
in this problem). Thus (12) already illustrates the adaptive nature of the LSE — the risk
of the LSE θ̂ changes depending on the “structure” of the true θ∗ . In the next subsection
(see (13)) we further highlight this adaptive nature of the LSE.
Remark 2.1. To the best of our knowledge, inequality (12) first appeared in [92, Theorem 1]
who proved it under the assumption that the errors ε1 , . . . , εn are i.i.d. N (0, σ 2 ). Zhang [129]
proved (12) for much more general errors including the case when ε1 , . . . , εn are i.i.d. with
mean zero and variance σ 2 . The proof we give (in Section 7.2) follows the arguments of
[129]. Another proof of an inequality similar to (12) for the case of normal errors has been
given recently by [23] who proved it as an illustration of a general technique for bounding
the risk of LSEs.
Remark 2.2. The LSE over bounded monotone functions also satisfies the bound (12) and
has been observed by many authors including [99, 121, 36]. Proving this result is easier,
however, because of the presence of the uniform bound on the function class (such a bound
is not present for the isotonic LSE). It must also be kept in mind that the bounded isotonic
LSE comes with a tuning parameter that needs to chosen by the user.
Remark 2.3. Inequality (12) also implies that the isotonic LSE achieves the risk (σ 2 V /n)2/3
for θ∗ ∈ IV := {θ ∈ I : θn − θ1 ≤ V } (as long as V is not too small) without any knowledge
of V . It turns out that the minimax risk over IV is of the order (σ 2 V /n)2/3 as long as V is
imsart-generic ver. 2014/10/16 file: ShapeRegArxiv.tex date: September 19, 2017

A. Guntuboyina and B. Sen/Shape-restricted Regression

8

√
in the range σ/ n . V . σn (see e.g., [25, Theorem 5.3]). Therefore, in this wide range of
V , the isotonic LSE is minimax (up to constant multiplicative factors) over the class IV .
This is especially interesting because the isotonic LSE does not require any knowledge of V .
This illustrates another kind of adaptation of the isotonic LSE; further details on this can
be found in [28].
2.2. Adaptive risk bounds
As the isotonic LSE fit is piecewise constant, it may be reasonable to expect that when θ∗
is itself a piecewise constant (with not too many pieces), the risk of θ̂ would be small. The
rightmost plot of Figure 1 corroborates this intuition. This leads us to our second type of
risk bound for the LSE. For θ ∈ I, let k(θ) ≥ 1 denote the number of constant blocks of θ,
i.e., k(θ) is the integer such that k(θ) − 1 is the number of inequalities θi ≤ θi+1 that are
strict, for i = 1, . . . , n − 1 (the number of jumps of θ).
Theorem 2.1. Under the assumption that ε1 , . . . , εn are i.i.d. with mean zero and variance
σ 2 we have


1 ∗
4σ 2 k(θ)
en
∗
2
R(θ̂, θ ) ≤ inf
kθ − θk +
log
(13)
θ∈I n
n
k(θ)
for every θ∗ ∈ Rn .
Note that θ∗ in Theorem 2.1 can be any arbitrary vector in Rn (it is not required that
∈ I). An important special case of inequality (13) arises when θ is taken to be θ∗ in
order to obtain:
en
4σ 2 k(θ∗ )
R(θ̂, θ∗ ) ≤
log
.
(14)
n
k(θ∗ )
θ∗

It makes sense to compare (14) with the worst case risk bound (12). Suppose, for example, θj∗ = 1{j > n/2} (here 1 denotes the indicator function) so that k(θ∗ ) = 2 and
V (θ∗ ) = 1. Then the risk bound in (12) is essentially (σ 2 /n)2/3 while the right side of (14)
is (8σ 2 /n) log(en/2) which is much smaller than (σ 2 /n)2/3 . More generally, if θ∗ is piecewise
constant with k blocks then k(θ∗ ) = k so that inequality (14) implies that the risk is given
by the parametric rate kσ 2 /n with a logarithmic multiplicative factor of 4 log(en/k) — this
is a much stronger bound compared to (12) when k is small.
Inequality (14) is an example of an oracle inequality. This is because of the following.
Let θ̂OR denote the oracle piecewise constant estimator of θ∗ which estimates θ∗ by the
mean of Y in each constant block of θ∗ (note that θ̂OR uses knowledge of the locations of
the constant blocks of θ∗ and hence is an oracle estimator). It is easy to see then that the
risk of θ̂OR is given by
σ 2 k(θ∗ )
R(θ̂OR , θ∗ ) =
.
n
As a result, inequality (14) can be rewritten as


en
∗
R(θ̂, θ ) ≤ 4 log
R(θ̂OR , θ∗ ).
(15)
k(θ∗ )
Because this involves a comparison of the risk of the LSE θ̂ with that of the oracle estimator θ̂OR , inequality (14) is referred to as an oracle inequality. Inequality (15) shows that
imsart-generic ver. 2014/10/16 file: ShapeRegArxiv.tex date: September 19, 2017

A. Guntuboyina and B. Sen/Shape-restricted Regression

9

the isotonic LSE, which uses no knowledge of k(θ∗ ) and the positions of the blocks, has
essentially the same risk performance as the oracle piecewise-constant estimator (up to the
multiplicative logarithmic factor 4 log(en/k(θ∗ ))). This is indeed remarkable!
Inequality (13) can actually be viewed as a more general oracle inequality where the
behavior of the isotonic LSE is compared with oracle piecewise constant estimators even
when θ∗ ∈
/ I. We would like to mention here that, in this context, (13) is referred to as a
sharp oracle inequality because the leading constant in front of the kθ∗ − θk2 /n term on the
right-hand side of (13) is equal to one. We refer to [16] for a detailed explanation of oracle
and sharp oracle inequalities.
Based on the discussion above, it should be clear to the reader that the adaptive risk
bound (13) complements the worst case bound (12) as it gives much finer information about
how well any particular θ∗ (depending on its ‘complexity’) can be estimated by the LSE θ̂.
Remark 2.4 (Model misspecification). As already mentioned, the sharp oracle inequality (13) needs no assumption on θ∗ (which can be any arbitrary vector in Rn ), i.e., the
inequality holds true even when θ∗ ∈
/ I. See [25, Section 6] for another way of handling
model misspecification, where θ̂ is compared with the “closest” element to θ∗ in I (and not
θ∗ ).
Remark 2.5. To the best of our knowledge, an inequality of the form (13) first explicitly
appeared in [25, Theorem 3.1] where it was proved that


4σ 2 k(θ)
en
1 ∗
∗
2
R(θ̂, θ ) ≤ 4 inf
kθ − θk +
log
(16)
θ∈I n
n
k(θ)
under the additional assumption that θ∗ ∈ I. The proof of this inequality given in [25] is
based on ideas developed in [129]. Note the additional constant factor of 4 in the above
inequality compared to (13).
Under the stronger assumption ε ∼ Nn (0, σ 2 In ), Bellec [16, Theorem 3.2] improved (16)
and proved that


1 ∗
σ 2 k(θ)
en
∗
2
R(θ̂, θ ) ≤ inf
kθ − θk +
log
,
(17)
θ∈I n
n
k(θ)
for every θ∗ ∈ Rn . A remarkable feature of this bound is that the multiplicative constants
involved are all tight, which implies, in particular, that


1 ∗
σ 2 k(θ)
en
∗
2
R(θ̂, θ ) ≤ inf
kθ − θk + C
log
θ∈I n
n
k(θ)
cannot hold for every θ∗ if C < 1. This follows from the P
fact that when θ∗ = (0, 0, . . . , 0) ∈ I
n
2
∗
2
2
and ε ∼ Nn (0, σ In ), the risk R(θ̂, θ ) exactly equals σ
j=1 1/j  σ log n; see [16] for an
explanation. It must be noted that this implies, in particular, that the logarithmic term in
these adaptive risk bounds cannot be removed.
Note that inequality (13) has an additional factor of 4 compared to (17) on the second
term in the right-hand side. This is because the errors ε1 , . . . , εn can be non-Gaussian in
Theorem 2.1.
Remark 2.6. One may attempt to prove (12) from the adaptive risk bound (13) by approximating arbitrary θ∗ ∈ I via θ ∈ I with a bound on k(θ). However, it is likely that such
imsart-generic ver. 2014/10/16 file: ShapeRegArxiv.tex date: September 19, 2017

A. Guntuboyina and B. Sen/Shape-restricted Regression

10

an approach will lead to additional logarithmic terms on the right hand side of (12) (see
e.g., [25, Theorem 4.1]).
2.2.1. Adaptive Risk Bounds for R(p) (θ̂, θ∗ ).
The risk bound (13) (or more specifically (15)) implies that the isotonic LSE pays a logarithmic price in risk compared to the oracle piecewise constant estimator. This fact is
strongly tied to the fact that the risk is measured via squared error loss (as in (9)). The
story will be different if one measures risk under `p -metrics for p 6= 2. To illustrate this, we
shall describe adaptive bounds for the risk R(p) (θ̂, θ∗ ) defined in (10).
The following result bounds the risk R(p) (θ̂, θ∗ ) assuming that θ∗ ∈ I. The risk bounds
involve a positive constant Cp that depends on p alone. Explicit expressions for Cp can be
gleaned from the proof of Theorem 2.2 (in Section 7.4).
Theorem 2.2. Assume that the errors ε1 , . . . , εn are i.i.d. N (0, σ 2 ). Fix θ∗ ∈ I and let
p ≥ 1, p 6= 2. Let k denote the number of constant blocks of θ∗ and let the lengths of the
blocks be denoted by n1 , . . . , nk . We then have
 min(p,2)/2
k
σ p X (2−p)+ /2
(p)
∗
p k
R (θ̂, θ ) ≤ Cp
(18)
ni
≤ Cp σ
n
n
i=1

where Cp is a positive constant that depends on p alone.
Remark 2.7. As stated, Theorem 2.2 appears to be new even though its conclusion is
implicit in the detailed risk calculations of [129] for isotonic regression. Zhang’s [129] arguments also seem to indicate that the bound (18) should be tight up to constant multiplicative
terms involving p. We have assumed that ε1 , . . . , εn are normal in Theorem 2.2 but it is
possible to allow non-Gaussian errors by imposing suitable moment conditions.
Let us now compare the isotonic LSE to the oracle piecewise constant estimator θ̂OR
(introduced in the previous subsection) in terms of the `p -risk. It is easy to verify that the
risk of θ̂OR under the `p -loss is given by
k

R(p) (θ̂OR , θ∗ ) = (E|η|p )σ p

1 X (2−p)/2
ni
n

(19)

i=1

for every p > 0 where η := ε1 /σ is standard normal.
Comparing (18) and (19), we see that the isotonic LSE performs at the same rate (up
to constant multiplicative factors) as the oracle piecewise constant estimator for 1 ≤ p < 2
(there is not even a logarithmic price for these values of p). When p = 2, as seen from (15),
the isotonic LSE pays a logarithmic price of 4 log(en/k(θ∗ )). For p > 2 however, there is
a significant price that is paid. For example, if all the constant blocks have roughly equal
size, then the oracle estimator’s risk, when p > 2, equals (k/n)p/2 while the bound in (18)
is of order k/n.
3. Risk bounds for the LSE in other shape-restricted regression problems
In this section we consider the problems of convex regression (Example 1.3), isotonic regression on a partially ordered set (Example 1.2), unimodal regression (Example 1.4) and shape
imsart-generic ver. 2014/10/16 file: ShapeRegArxiv.tex date: September 19, 2017

A. Guntuboyina and B. Sen/Shape-restricted Regression

11

restricted additive models (Example 1.5). In each of these problems, we describe results
related to the performance of the LSEs. The reader will notice that the risk results are not
as detailed as compared to the isotonic regression results of the previous section.
3.1. Convex Regression
Let us consider Example 1.3 where the goal is to estimate a convex function f : [0, 1] → R
from regression data as in (1). The convex LSE θ̂ is defined as the projection of Y onto the
closed convex cone K (see (6)). This estimator was first proposed in [69] for the estimation
of production functions and Engel curves. It can be shown that θ̂ is piecewise affine with
knots only at the design points; see [58, Lemma 2.6]. The accuracy of the LSE, in terms
of the risk R(θ̂, θ∗ ) (defined in (9)), was first studied in [64] followed by [25, 16, 24]. These
results are summarized below. Earlier results on the risk under a supremum loss can be
found in [66, 40].
Suppose that ε ∼ Nn (0, σ 2 In ). In [24], the following worst case risk bound for θ̂ was
given (when xi = i/n are the ordered design points):
∗

R(θ̂, θ ) ≤ C

σ2

!4/5
p
T (θ∗ )
σ2
+ C 4/5
n
n

(20)

where C > 0 is a universal constant and T (θ∗ ) is a constant depending on θ∗ (like V (θ∗ )
in (12) for isotonic regression). Roughly speaking, T (θ∗ ) measures the “distance” of θ∗ from
the set of all affine (functions) sequences. Formally, Let L denote the subspace of Rn spanned
by the constant vector (1, . . . , 1) and the vector (1, 2, . . . , n); i.e., L is the linear subspace
of affine sequences. Let PL denote the orthogonal projection matrix onto the subspace L
and let β ∗ := (In − PL )θ∗ . Then T (θ∗ ) := max1≤i≤n βi∗ − min1≤i≤n βi∗ . Observe that when
θ∗ itself is an affine sequence (which is also a convex sequence), then T (θ∗ ) = 0.
The risk bound (20) shows that the risk of the convex LSE is bounded above by n−4/5 .
Inequality (20) improved a result in [64], which had a similar bound but with an additional
multiplicative logarithmic factor (in n). Comparing with (12), it is natural to conjecture that
the second term in (20) can be improved to Cσ 2 (log(en))/n but this has not been proved
so far. Another feature of (20) is that the errors are assumed to be Gaussian; it might
be possible to extend them to sub-Gaussian errors but this is still a strong assumption
compared to the corresponding result for isotonic regression (see (12)) which holds without
distributional assumptions.
The proof of (20) (and other worst case risk bounds like (20) for shape-restricted regression problems under Gaussian/sub-Gaussian errors) involves tools from the theory of
Gaussian processes like chaining and Dudley’s entropy bound and crucially relies on an accurate ‘size’ measure of the underlying class (e.g., ‘local’ balls of K) as captured by its metric
entropy; see Section 3.5 for a broad outline of the proof strategy. Although the main idea
of the proof is simple, deriving appropriate bounds on the metric entropy of the underlying
class can be challenging.
As with the isotonic LSE, the convex LSE θ̂ exhibits adaptive behavior. As the convex
LSE θ̂ is piecewise affine it may be expected that the risk of θ̂ would be nearly parametric if
the true θ∗ is (well-approximated by) piecewise affine. Indeed this is the case. For θ ∈ K let
imsart-generic ver. 2014/10/16 file: ShapeRegArxiv.tex date: September 19, 2017

A. Guntuboyina and B. Sen/Shape-restricted Regression

12

q(θ) ≥ 1 denote the number of affine pieces of θ; i.e., q(θ) is an integer such that q(θ) − 1 is
the number of inequalities in (6) that are strict. This adaptive behavior can be illustrated
through the following risk bound:


1 ∗
8σ 2 q(θ)
en
∗
2
R(θ̂, θ ) ≤ inf
.
(21)
kθ − θk +
log
θ∈K n
n
q(θ)
This inequality has been proved by [16, Section 4] improving earlier results of [64, 25] which
had superfluous multiplicative constants. Note that this bound holds for ε ∼ N (0, σ 2 In ). It is
not known if the bounds holds for non-Gaussian errors (compare this with the corresponding
inequality (13) for isotonic regression which holds without distributional assumptions on
the errors). Let us also note that risk bounds for the LSE under the R(p) (θ̂, θ∗ ) risk (defined
in (10)) are not available for convex regression.
3.2. Isotonic regression on a partially ordered set
We now turn our attention to Example 1.2 where the covariates are partially ordered and the
goal is to estimate the order preserving (isotonic) regression function. The book Robertson
et al. [106, Chapter 1] gives a nice overview of the characterization and computation of LSEs
in such problems along with their applications in statistics. However, not much is known
in terms of rates of convergence for these LSEs beyond the example of coordinate-wise
nondecreasing ordering introduced in Example 1.2.
In this subsection we briefly review the main results in [26] which considers estimation of
a bivariate (d = 2) coordinate-wise nondecreasing regression function. An interesting recent
paper [65] has extended these results to all dimensions d ≥ 2 (see Remark 3.1). Estimation
of bivariate coordinate-wise nondecreasing functions has applications and connections to the
problem of estimating matrices of pairwise comparison probabilities arising from pairwise
comparison data ([29, 113]) and to seriation ([43]).
As the distribution of the design points xi complicate the analysis of shape-restricted
LSEs, especially when d > 1, for simplicity, we consider the regular uniform grid design.
∗ ) ∈ Rn1 ×n2 from
This reduces the problem to estimating an isotonic ‘matrix’ θ∗ := (θij
observations
∗
yij = θij
+ εij ,
for i = 1, . . . , n1 , j = 1, . . . , n2 ,
where θ∗ is constrained to lie in
M := {θ ∈ Rn1 ×n2 : θij ≤ θkl whenever i ≤ k and j ≤ l},
and the random errors εij ’s are i.i.d. N (0, σ 2 ), with σ 2 > 0 unknown. We refer to any
matrix in M as an isotonic matrix. Letting Y := (yij ) denote the matrix (of order n1 × n2 ;
n := n1 n2 ) of the observed responses, the LSE θ̂ is defined as the minimizer of the squared
Frobenius norm, kY − θk2 , over θ ∈ M, i.e.,
n1 X
n2
X
θ̂ := arg min
(yij − θij )2 .
θ∈M

(22)

i=1 j=1

As M is a closed convex cone in Rn1 ×n2 , the LSE θ̂ exists uniquely.
imsart-generic ver. 2014/10/16 file: ShapeRegArxiv.tex date: September 19, 2017

A. Guntuboyina and B. Sen/Shape-restricted Regression

13

The goal now is to formulate both the worst case and adaptive risk bounds for the matrix
isotonic LSE θ̂ in estimating θ∗ . In [26, Theorem 2.1] it was shown that
!
r
2
2 D(θ ∗ )
σ
σ
R(θ̂, θ∗ ) ≤ C
(log n)4 + (log n)8
(23)
n
n
∗ )2 . The above bound shows that
for a universal constant C > 0, where D(θ∗ ) := (θn∗ 1 n2 − θ11
when the ‘variation’ D(θ∗ ) of θ∗ is non-zero, the risk of θ̂ decays at the rate n−1/2 , while
when D(θ∗ ) = 0 (i.e., θ∗ is a constant), the risk is (almost) parametric. The above bound
probably has superfluous logarithmic factors (see [65]).
To describe the adaptive risk bound for the matrix isotonic LSE we need to introduce
some notation. A subset A of {1, . . . , n1 } × {1, . . . , n2 } is called a rectangle if A = {(i, j) :
k1 ≤ i ≤ l1 , k2 ≤ j ≤ l2 } for some 1 ≤ k1 ≤ l1 ≤ n1 and 1 ≤ k2 ≤ l2 ≤ n2 . A rectangular
partition of {1, . . . , n1 } × {1, . . . , n2 } is a collection of rectangles π = (A1 , . . . , Ak ) that are
disjoint and whose union is {1, . . . , n1 } × {1, . . . , n2 }. The cardinality of such a partition,
|π|, is the number of rectangles in the partition. The collection of all rectangular partitions
of {1, . . . , n1 } × {1, . . . , n2 } will be denoted by P. For θ ∈ M and π = (A1 , . . . , Ak ) ∈ P,
we say that θ is constant on π if {θij : (i, j) ∈ Al } is a singleton for each l = 1, . . . , k. We
define k(θ), for θ ∈ M, as the “number of rectangular blocks” of θ, i.e., the smallest integer
k for which there exists a partition π ∈ P with |π| = k such that θ is constant on π. In [26,
Theorem 2.4] the following adaptive risk bound was stated:
 ∗

kθ − θk2 Cσ 2 k(θ)
∗
8
R(θ̂, θ ) ≤ inf
+
(log n) .
(24)
θ∈M
n
n

where C > 0 is a universal constant.
In [26] the authors also established a property of the LSE that they termed ‘variable’
∗ ) ∈ I has the property
adaptation. Let In1 := {θ ∈ Rn1 : θ1 ≤ · · · ≤ θn1 }. Suppose θ∗ = (θij
∗ only depends on i, i.e., there exists θ ∗∗ ∈ I
∗
∗∗
that θij
n1 such that θij = θi for every i and
∗
j. If we knew this fact about θ , then the most natural way of estimating it would be
to perform
row-averages ȳ := (ȳ1 , . . . , ȳn1 ), where
P 2 vector isotonic estimation based on the
yij /n2 , resulting in an estimator θ̆ of θ∗∗ . Note that the construction of θ̆ requires
ȳi := nj=1
the knowledge that all rows of θ∗ are constant. As a consequence of the adaptive risk bound
(24), it was shown in [26, Theorem 2.4] that the matrix isotonic LSE θ̂ achieves the same
risk bounds as θ̆, up to additional logarithmic factors. This is remarkable because θ̂ uses no
special knowledge on θ∗ ; it automatically adapts to intrinsic dimension of θ∗ .
Remark 3.1 (Extension to d ≥ 2). The recent paper Han et al. [65] has studied ddimensional isotonic regression for general d ≥ 1 and proved versions of inequalities (23)
and (24). Their results imply, in particular, that the worst case risk behavior as well as
adaptive risk behavior of the isotonic LSE changes significantly for d ≥ 3 compared to the
behavior for d ≤ 2. They also obtained results for the isotonic LSE under random design
settings.
Let us reiterate that the bounds (23) and (24) are established under the assumption
that the errors εi,j are i.i.d. N (0, σ 2 ) (sub-Gaussianity). Nothing is known if they still hold
under general non-Gaussian errors. Also risk bounds in other loss functions (such as those
in appropriate `p -metrics) are not available.
imsart-generic ver. 2014/10/16 file: ShapeRegArxiv.tex date: September 19, 2017

A. Guntuboyina and B. Sen/Shape-restricted Regression

14

3.3. Unimodal Regression
In this subsection we summarize the two kinds of risk bounds known for the LSE in unimodal
(decreasing and then increasing) regression, introduced in Example 1.4. The unimodal LSE θ̂
is defined as any projection of Y onto U, a finite union of the closed convex cones described in
Example 1.4. It is known that θ̂ is piecewise constant with possible jumps only at the design
points. Once the mode of the fitted LSE is known (and fixed), θ̂ is just the nonincreasing
(isotonic) LSE fitted to the points to the left of the mode and nondecreasing (isotonic) LSE
fitted to the points on the right of the mode.
As in isotonic regression, the unimodal LSE θ̂ exhibits adaptive behavior. In fact, the
risk bounds for the unimodal LSE θ̂ are quite similar to those obtained for the isotonic LSE.
The two kinds of risk bounds are given below (under the assumption that ε ∼ Nn (0, σ 2 In )):


∗

R(θ̂, θ ) ≤ C
and



∗

R(θ̂, θ ) ≤ C inf

θ∈U

σ 2 S(θ∗ )
n

2/3
+C

σ2
,
n2/3

where θ∗ ∈ U,

σ 2 (k(θ) + 1)
en
1 ∗
kθ − θk2 + C
log
n
n
k(θ) + 1

(25)


(26)

where k(θ) is the number of constant blocks of θ, S(θ∗ ) := maxi,j |θi∗ − θj∗ | is the range of
θ∗ and C > 0 is a universal constant.
The worst case risk bound (25) is given in [27, Theorem 2.1] while the adaptive risk
bound (26) is a consequence of [16, Theorem A.4] (after integrating the tail probability).
The proof of (25) (given in [27, Theorem 2.1]) is based on the general theory of least squares
outlined in Section 3.5; also see [23, Theorem 2.2]. It shows that a unimodal regression
function can also be estimated at the same rate as a monotone function. The adaptive risk
bound (26), although being similar in spirit to that of the isotonic LSE, is weaker than (17)
(obtained for the isotonic LSE). Note that inequality (26) is not sharp (i.e., the leading
constant on the right side of (26) is not 1); in fact it is not known whether a sharp oracle
inequality can be constructed for R(θ̂, θ∗ ) (see [16]). The proof of the adaptive risk bound is
also slightly more involved than that of Theorem 2.1; the fact that the underlying parameter
space U is non-convex complicates the analysis.
3.4. Shape-restricted additive models
Given observations (x1 , y1 ), . . . , (xn , yn ) where {xi = (xij , 1 ≤ j ≤ d)}ni=1 are d-dimensional
design points and y1 , . . . , yn are real-valued, the additive model (see e.g., [68, 86]) assumes
that
d
X
∗
yi = µ +
fj∗ (xij ) + εi
for i = 1, . . . , n
j=1

µ∗

where
∈ R is an unknown intercept term, f1∗ , . . . , fd∗ are unknown univariate functions
satisfying
n
1X ∗
fj (xij ) = 0, for every j = 1, . . . , d,
(27)
n
i=1

imsart-generic ver. 2014/10/16 file: ShapeRegArxiv.tex date: September 19, 2017

A. Guntuboyina and B. Sen/Shape-restricted Regression

15

and ε1 , . . . , εn are unobserved mean-zero errors. Assumption (27) is necessary to ensure
the identifiability of f1∗ , . . . , fd∗ . We focus our attention to shape-restricted additive models
where it is assumed that each fj∗ obeys a known qualitative restriction such as monotonicity
or convexity which is captured by the assumption that fj∗ ∈ Fj for a known class of functions
Fj . One of the main goals in additive modeling is to recover each individual function fj∗ ∈ Fj
for j = 1, . . . , d.
The LSEs µ̂, fˆj of µ∗ , fj∗ , for j = 1, . . . , d are defined as minimizers of the sum of squares
criterion, i.e.,
n 
d
2
X
X
(µ̂, fˆ1 , . . . , fˆd ) := arg min
yi − µ −
fj (xij )
(28)
i=1

j=1

Pn

under the constraints µ ∈ R, fj ∈ Fj , i=1 fj (xij ) = 0 for j = 1, . . . , d. It is natural to
compare the performance of these LSEs to the corresponding oracle estimators defined in
the following way. For each k = 1, . . . , d, the oracle estimator fˆkOR is defined as
fˆkOR := arg min
fk

n 
X
i=1

yi − µ∗ −

X

2
fj∗ (xij ) − fk (xik ) ,

(29)

j6=k

Pn

where fk ∈ Fk and satisfies i=1 fk (xik ) = 0. In other words, fˆkOR assumes knowledge of
fj∗ , for j 6= k, and µ∗ , and performs least squares minimization only over fk ∈ Fk .
A very important aspect about shape-restricted additive models is that it is possible
for the LSE fˆk to be close to the oracle estimator fˆkOR , for each k = 1, . . . , d. Indeed,
this property was proved by Mammen and Yu [89] under certain assumptions for additive
isotonic regression where each function fj is assumed to be monotone. Specifically, [89]
worked with a random design setting where the design points are assumed to be i.i.d. from a
Lipschitz density that is bounded away from zero and infinity on [0, 1]d (this is a very general
setting which allows for non-product measures). They also assumed that each function
fj is differentiable and strictly increasing. Although the design restrictions in this result
are surprisingly minimal, we believe that the assumptions on the fj ’s can be relaxed. In
particular, this result should hold when fj ’s are piecewise constant and even under more
general shape restrictions such as convexity.
Our intuition is based on the following (trivial) observation that there exist design configurations where the LSE fˆk is remarkably close to fˆkOR for each k = 1, . . . , d under almost
no additional assumptions. The simplest such instance is when the design points take values
in a uniform grid. In this case, it is easy to see that fˆk is exactly equal to fˆkOR as shown
below. For simplicity, assume that d = 2 and that the design points are coming from a grid
in [0, 1]2 . It is convenient then to index the observations by (i1 , i2 ) where i1 = 1, . . . , n1 and
i2 = 1, . . . , n2 (n = n1 × n2 ). The design points will be given by xi1 ,i2 = (i1 /n1 , i2 /n2 ) and
the observation model can be written as
yi1 ,i2 = µ∗ + f1∗ (i1 /n1 ) + f2∗ (i2 /n2 ) + εi1 ,i2 ,

(30)

for i1 = 1, . . . , n1 , i2 = 1, . . . , n2 . The following result is proved in Section 7.5.
Lemma 3.1. Consider model (30) where fk∗ ∈ Fk , for k = 1, 2. Suppose that fˆ1 and fˆ2
denote the LSEs of f1∗ and f2∗ respectively, as defined in (28). Also, let the oracle estimators
be as defined fˆ1OR and fˆ2OR ; see (29). Then fˆ1 = fˆ1OR and fˆ2 = fˆ2OR .
imsart-generic ver. 2014/10/16 file: ShapeRegArxiv.tex date: September 19, 2017

A. Guntuboyina and B. Sen/Shape-restricted Regression

16

Note that we have made no assumptions at all on F1 and F2 . Thus under the regular grid
design in [0, 1]d , the LSEs of f1 and f2 are exactly equal to the oracle estimates. For general
design configurations, it might be much harder to relate the LSEs to the corresponding
oracle estimators. Nevertheless, the aforementioned phenomenon for gridded designs allows
us to conjecture that the closeness of fˆk to fˆkOR must hold in much greater generality than
has been observed previously in the literature.
It may be noted that the risk behavior of fˆkOR is easy to characterize. For example, when
fk∗ is assumed to be monotone, fˆkOR will satisfy risk bounds similar to those described in
Section 2. Likewise, when fk∗ is assumed to be convex, then fˆkOR will satisfy risk bounds
described in Subsection 3.1. Thus, when fˆk is close to fˆkOR (which we expect to happen
under a broad set of design configurations), it is natural to expect that fˆk will satisfy such
risk bounds as well.
3.5. General theory of LSEs
In this section, we collect some general results on the behavior of the LSEs that are useful for
proving the risk bounds described in the previous two sections. These results apply to LSEs
that are defined by (7) for a closed convex constraint set C. Convexity of C is crucial here (in
particular, these results do not directly apply to unimodal regression where the constraint
set is non-convex; see Section 3.3). We assume that the observation vector Y = θ∗ + ε for
a mean-zero random vector ε. Except in Lemma 3.4, we assume that ε ∼ Nn (0, σ 2 In ).
The first result reduces the problem of bounding R(θ̂, θ∗ ) to controlling the expected
supremum of an appropriate Gaussian process. This result was proved by Chatterjee [23]
(see [30, 122] for extensions to penalized LSEs).
Lemma 3.2 (Chatterjee). Consider the LSE (7) for a fixed closed convex set C. Assume
that Y = θ∗ + ε where ε ∼ Nn (0, σ 2 In ) and θ∗ ∈ C. Let us define the function gθ∗ : R+ → R
as
"
#
t2
gθ∗ (t) := E
sup
hε, θ − θ∗ i − .
(31)
2
θ∈C:kθ−θ∗ k≤t
Let tθ∗ be the point in [0, ∞) where t 7→ gθ∗ (t) attains its maximum (existence and uniqueness of tθ∗ are proved in [23, Theorem 1.1]). Then there exists a universal positive constant
C such that

C
R(θ̂, θ∗ ) ≤ max t2θ∗ , σ 2 .
(32)
n
Remark 3.2. Chatterjee [23] actually proved that kθ̂ − θ∗ k is tightly concentrated around
tθ∗ . The bound (32) is an easy consequence of this concentration result.
Lemma 3.2 reduces the problem of bounding R(θ̂, θ∗ ) to that of bounding tθ∗ . For this
latter problem, [23, Proposition 1.3] observed that
tθ∗ ≤ t∗∗

whenever t∗∗ > 0 and gθ∗ (t∗∗ ) ≤ 0.

In order to bound tθ∗ , one therefore seeks t∗∗ > 0 such that gθ∗ (t∗∗ ) ≤ 0. This now requires
a bound on the expected supremum of the Gaussian process in the definition of gθ∗ (t)
in (31). A simple upper bound for this expected Gaussian supremum is given by Dudley’s
imsart-generic ver. 2014/10/16 file: ShapeRegArxiv.tex date: September 19, 2017

A. Guntuboyina and B. Sen/Shape-restricted Regression

17

entropy bound (see e.g., [119, Chapter 2]) which is given below. This bound involves covering
numbers. For a subset K ⊆ Rn and  > 0, let N (, K) denote the -covering number of K
under the Euclidean metric k · k (i.e., N (, K) is the minimum number of closed balls of
radius  required to cover K). The logarithm of N (, K) is known as the -metric entropy
of K. Also, for each θ∗ ∈ C and t > 0, let
B(θ∗ , t) := {θ ∈ C : kθ − θ∗ k ≤ t}
denote the ball of radius t around θ∗ . Observe that the supremum in the definition in (31)
is over all θ ∈ B(θ∗ , t). Dudley’s entropy bound leads to the following upper bound for the
expected Gaussian supremum appearing in the definition of gθ∗ (t).
Lemma 3.3 (Chaining). For every θ∗ ∈ C and t > 0,
"
#
 Z 2t

p
√
∗
12
E
sup hε, θ − θ i ≤ σ inf
log N (, B(θ∗ , t)) d + 4δ n .
θ∈B(θ∗ ,t)

0<δ≤2t

δ

Remark 3.3. Dudley’s entropy bound is not always sharp. More sophisticated generic
chaining arguments exist which gives tight bounds (up to universal multiplicative constants)
for suprema of Gaussian processes; see [119].
Lemma 3.2 and Lemma 3.3 present one way of bounding R(θ̂, θ∗ ). This involves controlling the metric entropy of subsets of the constraint set C of the form B(θ∗ , t). This method
is useful but works only for the case of Gaussian/sub-Gaussian errors.
Let us now present another result which is useful for proving adaptive risk bounds under
misspecification. We shall now work with general error distributions for ε that are not
necessarily Gaussian (we only assume that E(ε) = 0). This result essentially states for
bounding R(θ̂, θ∗ ), it is possible to work with tangent cones associated with C instead of C.
It is easier to deal with cones as opposed to general closed convex sets which leads to the
usefulness of this result.
For a closed convex set C and θ ∈ C, the tangent cone of C at θ is defined as
TC (θ) := Closure{t(η − θ) : t ≥ 0, η ∈ C}.
Informally, TC (θ) represents all directions in which one can move from θ and still remain in
C. The following lemma relates the risk R(θ̂, θ∗ ) to tangent cones.
Lemma 3.4. Let C be a closed convex set in Rn . Let θ∗ ∈ Rn and suppose that Y = θ∗ + σZ
for some mean-zero random vector Z with EkZk2 < ∞. Then,
n



o
E kθ̂ − θ∗ k2 ≤ inf kθ∗ − θk2 + σ 2 E kΠTC (θ) (Z)k2 ,
(33)
θ∈C

where ΠTC (θ) (Z) denotes the projection of Z onto the closed convex cone TC (θ).
Some remarks on this lemma are given below.
Remark 3.4 (Statistical dimension). When Z ∼ Nn (0, In ) and K is a closed convex cone
in Rn , the quantity
"
#

2




δ(K) := E kΠK (Z)k2 = E hZ, ΠK (Z)i = E
sup hZ, θi
,
θ∈K:kθk≤1

imsart-generic ver. 2014/10/16 file: ShapeRegArxiv.tex date: September 19, 2017

A. Guntuboyina and B. Sen/Shape-restricted Regression

18

has been termed the statistical dimension of K by Amelunxen et al. [2]. Therefore, when Z ∼
Nn (0, In ),inequality (33) bounds the risk R(θ̂, θ∗ ) of the LSE via the statistical dimension
of the tangent cones TC (θ).
Remark 3.5 (No distributional assumptions). There are no distributional assumptions on
Z for (33) to hold. In particular, the components of Z can be arbitrarily dependent and
non-Gaussian (as long as EkZk2 < ∞). This result first appeared in [16, Proposition 2.1]
who stated it assuming that Z ∼ Nn (0, In ) but an examination of the proof reveals that the
Gaussian assumption is not necessary.
Remark 3.6. When θ∗ ∈ C, then one can take θ = θ∗ in the right side of (33) to deduce
that
Ekθ̂ − θ∗ k2 ≤ σ 2 EkΠTC (θ∗ ) (Z)k2 .
(34)
This inequality (34) was first proved by [101]. Bellec [16] extended it to the case when θ∗ ∈
/C
by proving Lemma 3.4.
Remark 3.7 (Tightness). A remarkable fact proved by Oymak and Hassibi [101] is that
lim
σ↓0

1
Ekθ̂ − θ∗ k2 = EkΠTC (θ∗ ) (Z)k2
σ2

when θ∗ ∈ C.

(35)

Analogues of this inequality when θ∗ ∈
/ C have been recently proved in [42]. The equality
in (35) implies that if rn (θ), for θ ∈ C, is any rate term controlling the adaptive behavior
of the LSE in the following sense:

Ekθ̂ − θ∗ k2 ≤ inf kθ∗ − θk2 + σ 2 rn (θ)
for every θ∗ ∈ Rn
θ∈C

then it necessarily must happen that
rn (θ) ≥ EkΠTC (θ) (Z)k2

for every θ ∈ C.

Thus one loses nothing when working with tangent cones (i.e., focussing on bounding EkΠTC (θ) (Z)k2 )
to prove adaptive risk bounds.
We shall show how to apply Lemma 3.4 to prove the adaptive risk bound (13) in Section
7.1. Lemma 3.4 is also crucially used in [16] to prove the adaptive risk bound (21) for convex
regression. Lemma 3.4 also has applications beyond shape-restricted regression. It has been
recently used to prove risk bounds for total variation denoising and trend filtering (see [63]).
4. Inference with shape-restricted functions
Till now we have focused our attention on (global) risk properties of shape-restricted LSEs.
In this section we investigate the pointwise limiting behavior of the estimators. By the
pointwise behavior we mean the distribution of the LSE fˆn at a fixed point (say t), properly
normalized. Developing asymptotic distribution theory for the LSEs turns out to be rather
non-trivial, mainly because there is no closed form simple expression for the LSEs; all the
properties of the estimator have to be teased out from the general characterization (8).
The LSEs exhibit non-standard asymptotics: The limiting distributions that arise are
non-normal (and the rates of convergence are slower than n−1/2 ) and involve many nuisance
imsart-generic ver. 2014/10/16 file: ShapeRegArxiv.tex date: September 19, 2017

A. Guntuboyina and B. Sen/Shape-restricted Regression

19

parameters (that are difficult to estimate). As before, analyzing the isotonic LSE is probably
the simplest, and we will work with this example in Section 4.1. In Section 4.2 we develop
bootstrap and likelihood based methods for constructing (asymptotically) valid pointwise
confidence intervals, for the isotonic regression function f , that bypass estimation of nuisance
parameters. Section 4.3 deals with the case when f is convex — we sketch a proof of the
pointwise limiting distribution of the convex LSE. Not much is known in this area beyond
d = 1 for any of the shape-restricted LSEs discussed in the Introduction.
4.1. Pointwise limit theory of the LSE in isotonic regression
Let us recall the setup in (1) where f is now an unknown nondecreasing function. Further,
for simplicity, let xi = i/n, for i = 1, . . . , n, be the ordered design points and we assume that
ε1 , . . . , εn are i.i.d. mean zero errors with finite variance σ 2 > 0. The above assumptions can
be relaxed substantially, e.g., we can allow for dependent, heteroscedastic errors and the
xi ’s can be any sequence whose empirical distribution converges to a probability measure
on [0, 1]; see e.g., [3], [123, Section 3.2.15].
We start with another useful characterization of the isotonic LSE ([15, Theorem 1.1]).
Define the cumulative sum diagram (CSD) as the continuous piecewise affine function Fn :
[0, 1] → R (with possible knots only at i/n, for i = 1, . . . , n) for which
Fn (0) := 0,

and Fn

i

i

1X
:=
yj ,
n
n

for i = 1, . . . , n.

(36)

j=1

For any function g : I → R, where I ⊂ R is an interval, we denote by g̃ the greatest convex
minorant (GCM) of g (on I), i.e., g̃ is the largest convex function sitting below g. Thus,
F̃n denotes the GCM of Fn (on the interval [0, 1]). Let fˆn : (0, 1] → R be defined as the
left-hand derivative of the GCM of the CSD; i.e.,
fˆn := [F̃n ]0 ≡ F̃n0 ,
the left-hand slope of F̃n . Then, it can be shown that (see e.g., [106, Chapter 1]) the isotonic
LSE θ̂ is given by θ̂i = fˆn (i/n), for i = 1, . . . , n. Figure 2 illustrates these concepts from a
simple simulation.
Fix 0 < t < 1 and suppose that f has a positive continuous derivative f 0 on some neighborhood of t. The goal is to find the asymptotic distribution of fˆn (t), properly normalized.
We will show that
d
∆n := n1/3 {fˆn (t) − f (t)} → κC
(37)
where C := arg minh∈R {W(h) + h2 } has Chernoff’s distribution (here W(·) is a two-sided
Brownian motion starting from 0) and κ := [4σ 2 f 0 (t)]1/3 ; see e.g., [21, 126, 50, 51]. The
first result of this type was derived in [104] for the Grenander estimator — the maximum
likelihood estimator of a nonincreasing density in [0, ∞) (see [49]). Note that the Chernoff’s
random variable C is pivotal and its quantiles are known; see e.g., [33, 62].
4.1.1. Outline of a proof of (37).
In the following we sketch an outline of a proof of (37). Our proof technique directly
appeals to the characterization of the isotonic LSE as described at the beginning of this
imsart-generic ver. 2014/10/16 file: ShapeRegArxiv.tex date: September 19, 2017

A. Guntuboyina and B. Sen/Shape-restricted Regression
CSD and GCM

0.0

0.0

0.2

0.4

0.4

0.8

0.6

Scatter plot

20

0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.2

0.4

0.6

0.8

1.0

Fig 2: The left panel shows the scatter plot with the fitted function fˆn (in red) and the true f
(in blue) while the right panel shows the CSD (dashed) along with its GCM (in red). Here
n = 10, f (x) = x and ε ∼ Nn (0, σ 2 In ) with σ = 0.5.

section; see [123, Section 3.2.15] for an alternative proof technique that uses the switching
relationship, due to Groeneboom [51]. Let us further assume that the i.i.d. errors εi ’s have
a finite moment generating function near 0. This assumption lets us avoid the use of heavy
empirical process machinery and, we hope, will make the main technical arguments simple
and accessible to a broader audience.
We consider the stochastic process
Zn (h) := n2/3 [Fn (t + n−1/3 h) − Fn (t) − n−1/3 hf (t)],
for h ∈ [−tn1/3 , (1 − t)n1/3 ]. We regard stochastic processes as random elements in D(R),
the space of right continuous functions on R with left limits, equipped with the projection
σ-field and the topology of uniform convergence on compacta; see [102, Chapters IV and V]
for background.
Noting that if u is a bounded function and v is affine then u]
+ v = ũ + v, we have
Z̃n (h) := n2/3 [F̃n (t + n−1/3 h) − Fn (t) − n−1/3 hf (t)],
for h ∈ [−tn1/3 , (1 − t)n1/3 ]. By taking the left-derivative of the above process at h = 0 we
get (w.p. 1),
∆n = [Z̃n ]0 (0).
(38)
The above relation is crucial, as it relates ∆n , the quantity of interest, to a functional of the
process Zn . We study the process Zn (and show its convergence) and apply a (version of)
‘continuous’ mapping theorem (see e.g., [74, pp. 217-218]) to derive the limiting distribution
of ∆n .
Let F̌n : [0, 1] → R be the continuous piecewise affine function (with possible knots only
at i/n, for i = 1, . . . , n) with
F̌n

i
n

i

:=

1 X j 
f
,
n
n

for i = 0, . . . , n,

j=1

imsart-generic ver. 2014/10/16 file: ShapeRegArxiv.tex date: September 19, 2017

A. Guntuboyina and B. Sen/Shape-restricted Regression

21

and let F : [0, 1] → R be defined as
x

Z
F (x) :=

f (s)ds.
0

To study the stochastic process Zn we decompose Zn into the sum of the following three
terms:
Zn,1 (h) := n2/3 [Fn (t + n−1/3 h) − F̌n (t + n−1/3 h) − Fn (t) + F̌n (t)],
Zn,2 (h) := n2/3 [F̌n (t + n−1/3 h) − F (t + n−1/3 h) − F̌n (t) + F (t)],
Zn,3 (h) := n2/3 [F (t + n−1/3 h) − F (t) − n−1/3 hf (t)],
Observe that Fn − F̌n is just the partial sum process, properly normalized. By the Hungarian
embedding theorem (see e.g., [75]) we know that the partial sum process is approximated
by a Brownian motion process such that
Fn (x) − F̌n (x) = n−1/2 σBn (x) + Rn (x),

(39)

where Bn is a Brownian motion on [0,1] and
sup |Rn (x)| = O

 log n 
n

x

w.p. 1.

Thus,
Zn,1 (·) = σWn (·) + op (1),
where the process Wn is defined as Wn (h) := n1/6 {Bn (t + n−1/3 h) − Bn (t)}, h ∈ R, and
Wn ∼ W with W being distributed as a two-sided Brownian motion (starting at 0). This
shows that the process Zn,1 converges in distribution to W.
To study Zn,2 , observe that as f (·) is continuously differentiable in a neighborhood N
around t, we have (by a simple interpolation bound)
sup |F̌n (x) − F (x)| = O(n−1 ).
x∈N

Thus, Zn,2 converges to the zero function. By a simple application of Taylor’s theorem, we
can show that Zn,3 converges, uniformly on compacta, to the function D(h) := h2 f 0 (t)/2.
Combining the above results, we obtain that Zn converges in distribution to the process
Z(h) := σW(h) + h2 f 0 (t)/2, i.e.,
d

Zn → Z
in the topology of uniform convergence on compacta. Thus, it is reasonable to expect that
d

∆n = [Z̃n ]0 (0) → [Z̃]0 (0).
However, a rigorous proof of the convergence in distribution of ∆n involves a little more
than an application of a continuous mapping theorem. The convergence of Zn to Z is only
under the metric of uniform convergence on compacta. However, the GCM near the origin
might be determined by values of the process far away from the origin; the convergence
Zn to Z itself does not imply the convergence of [Z̃n ]0 (0) to [Z̃]0 (0). We need to show that
∆n is determined by values of Zn (h) for h in an Op (1) neighborhood of 0; see e.g., [74,
pp. 217-218] for such a result with a detailed proof.
d

It can be shown that [Z̃]0 (0) = κC (see e.g., [55, Chapter 3.2] and [55, Exercise 3.12])
which completes the proof sketch of (37).
imsart-generic ver. 2014/10/16 file: ShapeRegArxiv.tex date: September 19, 2017

A. Guntuboyina and B. Sen/Shape-restricted Regression

22

4.1.2. Other asymptotic regimes.
Observe that the assumption f 0 (t) 6= 0 is crucial in deriving the limiting distribution in (37).
One may ask, what if f 0 (t) = 0? Or even simply, what if f is a constant function on [0,1]?
In the latter case, we can easily show that, as a stochastic process in [0, 1],
√

d

n{fˆn (t) − f (t)} → σ[B̃]0 (t),

where B is the standard Brownian motion on [0, 1]. The above holds because of the following
√
√
observations. First note that n{fˆn (t)−f (t)} is the left-hand slope of the GCM of n(Fn −
√
F ) at t (as F is now linear). As n(Fn − F ) converges in distribution to the process σB on
D[0, 1], we have
√
√
d
0
0
n{fˆn (t) − f (t)} = n[F^
n − F ] (t) → σ[B̃] (t).
The above heuristic can be justified rigorously; see e.g., [55, Section 3.2]. In the related (nonincreasing) density estimation problem, [51, 22] showed that if f (t) lies on a flat stretch of
the underlying function f then the LSE (which is also the nonparametric maximum likelihood estimator, usually known as the Grenander estimator) converges to a non-degenerate
limit at rate n−1/2 , and they characterized the limiting distribution.
If one assumes that f (j) (t) = 0, for j = 1, . . . , p−1, and f (p) (t) 6= 0 (for p ≥ 1), where f (j)
denotes the j’th derivative of f , then one can derive the limiting distribution of fˆn (t), which
now converges at the rate n−p/(2p+1) ; see e.g., [126, 81]. Note that all the above scenarios
illustrate that the rate of convergence of the isotonic LSE fˆn (t) crucially depends on the
the behavior of f around t; this demonstrates the adaptive behavior of the isotonic LSE
from a pointwise asymptotics standpoint.
4.2. Constructing asymptotically valid pointwise confidence intervals
Although (37) gives the asymptotic distribution of the isotonic LSE at the point t, it is not
immediately clear how it can be used to construct a confidence interval for f (t) — the limiting distribution involves the nuisance parameter f 0 (t) that needs to be estimated. A naive
approach would suggest plugging in an estimator of f 0 (t) in the limiting distribution in (37)
to construct an approximate confidence interval. However, as fˆn is a piecewise constant
function, fˆn0 is either 0 or undefined and cannot be used to estimate f 0 (t) consistently. This
motivates the use of bootstrap and likelihood ratio based methods to construct confidence
intervals for f (t). In the following we just assume that ε1 , . . . , εn i.i.d. mean zero errors with
finite variance.
4.2.1. Bootstrap based inference.
Let us revisit (37) and consider the problem of bootstrapping fˆn to estimate the distribution
of ∆n ∼ Hn (say). Suppose that Ĥn is an approximation of Hn (which will be obtained
from bootstrap in this subsection) that can be computed. Then, an approximate 1 − α
(0 < α < 1) confidence interval for f (t) would be
[fˆn (t) − q̂1−α/2 n−1/3 , fˆn (t) − q̂α/2 n−1/3 ],
imsart-generic ver. 2014/10/16 file: ShapeRegArxiv.tex date: September 19, 2017

A. Guntuboyina and B. Sen/Shape-restricted Regression

23

where q̂α denotes the α’th quantile of Ĥn .
In a regression setup there are two main bootstrapping techniques: ‘bootstrapping pairs’
and ‘bootstrapping residuals’. Bootstrapping pairs refers to drawing with replacement samples from the data {(xi , yi ) : i = 1, . . . , n}; it is more natural when we have i.i.d. bivariate
data from a joint distribution. The residual bootstrap procedure fixes the design points xi ’s
and draws
yi∗ := fˇn (xi ) + ε∗i ,
i = 1, . . . , n
(the ∗ indicates a data point in the bootstrap sample) where fˇn is a natural estimator of
f in the model, and ε∗i ’s are i.i.d. (conditional on the data) having the distribution of the
(centered) residuals {yi − fˇn (xi ) : i = 1, . . . , n}. Let fˆn∗ denote the isotonic LSE computed
from the bootstrap sample. The bootstrap counterpart of ∆n (cf. (37)) is
∆∗n := n1/3 {fˆn∗ (t) − fˇn (t)}.
We now approximate Hn by Ĥn , the conditional distribution of ∆∗n , given the data. Note
that a natural candidate for fˇn in isotonic regression is the LSE fˆn .
Will this bootstrap approximation (by Ĥn ) work? This brings us to the notion of consistency of the bootstrap. Let d denote the Levy metric or any other metric metrizing weak
convergence of distributions. We say that Ĥn is weakly consistent if d(Hn , Ĥn ) → 0 in probability. If the convergence holds with probability 1, then we say that the bootstrap is strongly
consistent. If Hn has a weak limit H, then consistency requires Ĥn to converge weakly to
H, in probability; and if H is continuous, consistency requires supx∈R |Ĥn (x) − H(x)| → 0
in probability.
It is well-known that both the above bootstrap schemes — bootstrapping pairs and
bootstrapping residuals with fˇn = fˆn — yield inconsistent estimators of Hn ; see [1, 110, 76,
112, 52]. Intuitively, the inconsistency of the residual bootstrap procedure can be attributed
to the lack of smoothness of fˆn . Indeed a version of the residual bootstrap where one
considers fˇn as a smoothed version of fˆn (that can approximate the nuisance parameter
f 0 (t) consistently) can be shown to be consistent; see e.g., [112]. Specifically, suppose that
fˇn is a sequence of estimators such that
lim sup fˇn (x) − f (x) = 0,

n→∞ x∈I

(40)

almost surely, where I ⊂ [0, 1] is an open neighborhood of t, and
lim sup n1/3 fˇn (t + n−1/3 h) − fˇn (t) − f 0 (t)n−1/3 h = 0

n→∞ h∈K

(41)

almost surely for any compact set K ⊂ R. It can be shown, using arguments similar to
those in the proof of [112, Theorem 2.1], that if (40) and (41) hold then, conditional on the
data, the bootstrap estimator ∆∗n converges in distribution to κC, as defined in (37), almost
surely. Thus, this bootstrap scheme is strongly consistent.
A natural question that arises now is: Can we construct a smooth fˇn such that (40)
and (41) hold w.p. 1? We briefly describe such a smoothed bootstrap scheme. Let k(·)
be a differentiable symmetric Rdensity (kernel) with compact support (e.g., k(x) ∝ (1 −
x
x2 )2 1[−1,1] (x)) and let K(x) := −∞ k(s) ds be the corresponding distribution function. Let h
imsart-generic ver. 2014/10/16 file: ShapeRegArxiv.tex date: September 19, 2017

A. Guntuboyina and B. Sen/Shape-restricted Regression

24

be a smoothing parameter. Note that h may depend on the sample size n but, for notational
convenience, we write h instead of hn . Let kh (x) := k(x/h)/h and Kh (x) := K(x/h). Then
the smoothed isotonic LSE of f is defined as (cf. [60])
Z
ˇ
ˇ
fn (x) ≡ fn,h (x) := Kh (x − s) dfˆn (s), x ∈ [0, 1].
It can be easily seen that fˇn is a nondecreasing function (if t2 > t1 , then Kh (t2 −s) ≥ Kh (t1 −
s) for all s). Observe that fˇn is a smoothed version of the step function fˆn . In [112] it is shown
that the obtained bootstrap procedure is strongly consistent, i.e., ∆∗n = n1/3 {fˆn∗ (t) − fˇn (t)}
converges weakly to κC, conditional on the data, almost surely.
It is natural to conjecture that a (suitably) smoothed bootstrap procedure would also
yield (asymptotically) valid pointwise confidence intervals for other shape-restricted regression functions (e.g., convex regression). Moreover, it can be expected that the naive ‘with
replacement’ bootstrap and the residual bootstrap using the LSE would lead to inconsistent
procedures. However, as far as we are aware, there is no work that rigorously proves these
claims.
4.2.2. Likelihood ratio based inference.
Banerjee and Wellner [13] proposed a novel method for constructing pointwise confidence
intervals for a monotone function (e.g., f ) that avoids the need to estimate nuisance parameters; also see [11, 56]. Specifically, the strategy is to consider the testing problem
H0 : f (t) = φ0 versus H1 : f (t) 6= φ0 , where φ0 ∈ R is a known constant, using the likelihood ratio statistic (LRS), constructed under the assumption of i.i.d. Gaussian errors. If
one could find the limiting distribution of the LRS under the null hypothesis and show that
the limit is pivotal (as is the case in parametric models where the limiting distribution turns
out to be χ2 ) then that would provide a convenient way to construct a confidence interval
for f (t) via the method of inversion: an asymptotic level 1 − α confidence set would be given
by the set of all φ0 ’s for which the null hypothesis H0 : f (t) = φ0 is accepted.
To study the form of the LRS, we first need to understand the constrained isotonic LSE.
Consider the setup introduced in the beginning of Subsection 4.1 and suppose that l := bntc,
so that l/n ≤ t < (l + 1)/n. Under H0 : f (t) = φ0 , the constrained isotonic LSE fˆn0 is given
by
n
X
 0
fˆn (i/n) : i = 1, . . . , n :=
arg min
(Yi − θi )2 .
θ∈Rn :θ1 ≤···≤θl ≤φ0 ≤θl+1 ≤···≤θn i=1

Note that both functions fˆn and fˆn0 are identified only at the design points. By convention,
we extend them as left-continuous piecewise constant functions defined on the entire interval
(0, 1]. The hypothesis test is based on the following LRS:
Ln :=

n
X

n
2 X
2
0
ˆ
Yi − fn (i/n) −
Yi − fˆn (i/n) .

i=1

i=1

As shown in [11, 12] (in the setting of random design, which can be easily generalized to
cover the uniform grid design; see [6]), if f (t) = φ0 and f 0 (t) 6= 0, then
d

Ln → σ 2 L,
imsart-generic ver. 2014/10/16 file: ShapeRegArxiv.tex date: September 19, 2017

A. Guntuboyina and B. Sen/Shape-restricted Regression

25

where L is a nonnegative random variable expressible as a functional of two-sided Brownian
motion plus quadratic drift {W(h) + h2 : h ∈ R}, and σ 2 is the common variance of the
errors. An important feature of this limiting distribution is that it is pivotal — free of
the parameters of the problem. This readily yields confidence sets for f (t) (obtained by
the method of inversion) that do not need estimation of the nuisance parameter f 0 (t) — a
challenging quantity to estimate in practice. However, an estimate of σ 2 is required, which
can be easily obtained: The natural estimator kY − θ̂k2 /n of σ 2 is asymptotically normal
with mean σ 2 and variance 2σ 4 /n (see [92, Proposition 3]). This methodology has been
applied successfully in several monotone function estimation problems; see e.g., [13, 11, 56].
The method, and extensions thereof, also applies to both short- and long-range dependence
regimes for the errors; see [6]. Also see [14, 6] for illustrations and examples of the superior
performance of the LR based method over plug-in methods for constructing confidence
intervals for f (t), especially when the estimation of the derivative f 0 (t) is difficult.
Not much is known about the (asymptotic) distribution of the LRS beyond monotone
function estimation problems. However, in the recent papers [37, 38] the authors study
the LRS for testing the location of the mode of a log-concave density f using the unconstrained/constrained maximum likelihood estimator of a log-concave density (also see [10])
and show that, under the null hypothesis which fixes the value of the mode of f (and
assumes strict curvature of − log f at the mode), the LRS is asymptotically pivotal.
4.3. Pointwise limit theory of the LSE in convex regression
We assume that we have data from (1) where f : [0, 1] → R is now assumed to be convex (see
Example 1.3). We study the pointwise asymptotic theory for the convex LSE. For simplicity,
as before, we consider equi-spaced design points. Let us first describe the characterization
of the convex LSE that will drive the asymptotic analysis. Given the convex LSE θ̂, let
Θ̂ =P(Θ̂1 , . . . , Θ̂n ) denote the vector of its cumulative sums (divided by n), i.e., Θ̂i :=
n−1 ij=1 θ̂j , for i = 1, . . . , n, and recall Fn , as defined in (36). Then, θ̂ has the following
characterization: θ̂ is the unique vector such that Θ̂n = Fn (1) and
( P
j−1
X
≥ j−1
Fn (i/n) for j = 2, . . . , n,
Θ̂i
(42)
Pi=1
j−1
=
F
(i/n)
if
θ̂
has
a
kink
at
j/n
or
j
=
n;
n
i=1
i=1
this follows from the characterization of projection on the closed convex set K (as defined
in (6)); see [58, Lemma 2.6] for a complete proof. We define the convex LSE fˆn : [1/n, 1] → R
of f as the piecewise linear interpolation of the points {(i/n, θ̂i ) : i = 1, . . . , n}.
Fix t ∈ (0, 1) and consider the estimation of f (t) using the convex LSE fˆn (t) under
the assumption that f 00 is continuous and nonzero in a neighborhood of t. If the errors are
i.i.d. sub-Gaussian with mean zero and f 00 (t) 6= 0, the rate of convergence of fˆn (t) is known
to be n−2/5 (see [85]). The pointwise asymptotic distribution of the convex LSE (properly
normalized) is derived in [58]. In Groeneboom et al. [58] the authors show that
d

∆n := n2/5 {fˆn (t) − f (t)} → H00 (0),

(43)

where H is the “invelope” of integrated Brownian motion with quartic drift (+h4 ), and H00 (0)
is the second derivative of H at 0 (which exists w.p. 1). The invelope is a cubic spline lying
imsart-generic ver. 2014/10/16 file: ShapeRegArxiv.tex date: September 19, 2017

A. Guntuboyina and B. Sen/Shape-restricted Regression

26

above and touching integrated Brownian motion +h4 ; compare this with the “envelope” of
Brownian motion with a parabolic drift (+h2 ) that appears when analyzing the isotonic
LSE (see Section 4.1.1). Although a rigorous proof of the above weak convergence is long
and delicate (see [58, Theorem 6.3]), the main intuition for such a limit can be gotten
from looking at the characterization given in (42). We describe some of the main ideas
below. The first step is to show that the characterization in (42) can be ‘localized’ in an
appropriate sense. Then we show that the right side of the inequality in (42), appropriately
localized and normalized, converges to a limiting process involving integrated Brownian
motion +h4 . Then, a continuous mapping-like result, where we look at the limiting version
of the localized (42), yields the convergence of ∆n .
Let Θ̂n : [0, 1] → R be the continuous piecewise affine function (with possible knots only
at i/n) such that Θ̂n (i/n) = Θ̂i , for i = 1, . . . , n, and Θ̂n (0) = 0. Let Θ̃n : [0, 1] → R be
defined as
Z x
Θ̃n (x) :=
fˆn (s) ds.
0

Note that Θ̂n and Θ̃n are asymptotically the same, but it is easier to study Θ̃n . Let Gn
denote the empirical distribution function of the design points {i/n : i = 1, . . . , n}. Further,
loc
let us define the stochastic processes Hloc
n and H̃n (on R) as
Hloc
n (h)

:= n

4/5

t+n−1/5 h

Z
t

H̃loc
n (h)

:= n

4/5

Z

v


Z
Θ̂n (v) − Θ̂n (t) −



0
f (t) − (u − t)f (t) dGn (u) dv + An t + Bn ,

t
t+n−1/5 h


Z
Θ̃n (v) − Θ̃n (t) −

t

v




f (t) − (u − t)f (t) du dv + An t + Bn ,
0

t

where
An := n

3/5

{Θ̂n (t) − Fn (t)},

and Bn := n

4/5

Z

t

{Θ̂n (s) − Fn (s)} ds.
0

The process Hloc
n can be thought of as the ‘localization’ of the left side of (42). The process
H̃loc
is
important
for the following reason. As Θ̃0n (x) = fˆn (x), for x ∈ (0, 1), by differentiating
n
loc
the process H̃n twice we get
00
2/5 ˆ
(H̃loc
{fn (t + n−1/5 h) − f (t) − f 0 (t)n−1/5 h}.
n ) (h) = n

Thus, the quantity of interest ∆n is related to the process H̃loc
n as
00
∆n = (H̃loc
n ) (0),
loc
and this motivates the study of the process H̃loc
n (cf. (38)). Further, one can show that H̃n
loc
and Hn are asymptotically the same; see [58, p. 1696].
The following process can be thought of as the localization of the right side of (42):

4/5
Zloc
n (h) := n

Z
t

t+n−1/5 h


Z
Fn (v) − Fn (t) −

v



f (t) − (u − t)f 0 (t) dGn (u) dv.

t

imsart-generic ver. 2014/10/16 file: ShapeRegArxiv.tex date: September 19, 2017

A. Guntuboyina and B. Sen/Shape-restricted Regression

27

Moreover, as we will show, the process Zloc
n converges to a limiting distribution and is
loc ) as, for all h ∈ R,
related to Hloc
(and
thus
to
H̃
n
n
Hloc
n (h)

−

Zloc
n (h)

= n

4/5

Z

t+n−1/5 h

{Θ̂n (v) − Fn (v)}dv + Bn
t

= n

4/5

Z

t+n−1/5 h

{Θ̂n (v) − Fn (v)}dv ≥ 0,
0

with equality if t + n−1/5 h is a kink point. As the process Zloc
n involves empirical averages
we can use empirical process techniques (or the Hungarian embedding theorem; see (39))
and Taylor’s expansion to show that
Z h
1
d p
Zloc
(h)
→
f
(t)
W(s)ds + f 00 (t)h4 := Z(h), h ∈ [−K, K],
n
24
0
as a stochastic process in D[−K, K], for any K > 0, under the metric of uniform convergence
on [−K, K], where W is a two-sided Brownian motion; see [58, pp. 1694–1696].
By [58, Theorem 6.1] there exists an almost surely uniquely defined random continuous
process H, called an invelope of the process Z, such that (i) H(h) ≥ Z(h), for each h ∈ R;
(ii) H has a convexRsecond derivative, and, with probability 1, H is three times differentiable
at h = 0; and (iii) {H(h) − Z(h)}dH(3) (h) = 0 (which signifies that H = Z on the set where
H(3) has a jump). Further, it is shown in [58, pp. 1689–1692] that along with the process
loc
Zloc
n , the “invelope” Hn converges in such a way that the second and third derivatives of
Hloc
n at zero converges in distribution to the corresponding quantities of H; also see [57].
loc
loc
Thus, H̃loc
n converges weakly to H (as H̃n and Hn are asymptotically equivalent) and
00
00
∆n = (H̃loc
n ) (0) (see (43)) converges weakly to H (0).
Remark 4.1 (Multi-scale inference in shape-restricted problems). The pointwise asymptotic theory for isotonic and convex regression is developed under suitable smoothness assumptions on f , e.g., (37) needs f 0 (t) 6= 0 whereas the weak convergence of (43) assumes
f 00 (t) 6= 0. In [39], utilizing suitable multiscale tests, the author constructs confidence bands
for f that are locally adaptive in a certain sense (to the underlying smoothness in f ) and
have guaranteed coverage, assuming that f is isotonic or convex. These confidence bands
are computationally feasible and are also shown to be asymptotically sharp optimal in an
appropriate sense.
5. Computation of the LSE
In this section we discuss the computation of the LSE θ̂ in nonparametric shape-restricted
regression problems. Note that in most cases (see e.g., Examples 1.1–1.5) the LSE θ̂ is the
projection of Y onto C, a (finite union of) closed convex set(s) in Rn . If C is a polyhedral
convex set, then the computation of θ̂ involves solving a quadratic program with a bunch of
linear constraints. Many off-the-shelf solvers — e.g., CPLEX, MOSEK, Gurobi — can solve
these quadratic programs easily even for moderately large sample sizes (e.g., n ≈ 105 ). In
the following we consider the main examples in the Introduction and discuss some problem
specific algorithms that are computationally more efficient.
imsart-generic ver. 2014/10/16 file: ShapeRegArxiv.tex date: September 19, 2017

A. Guntuboyina and B. Sen/Shape-restricted Regression

28

Isotonic and unimodal regression. For the monotone regression problem [15] presented a graphical interpretation of the isotonic LSE (defined in (11)) in terms of the GCM
of the CSD; see Section 4.1. The method of successive approximation to the GCM can
be described algebraically as the pool-adjacent-violators algorithm (PAVA); see e.g., [106,
p. 9-10]. Roughly speaking, PAVA works as follows. We start with y1 on the left. We move
to the right until we encounter the first violation yi > yi+1 . Then we replace this pair by
their average, and back-average to the left as needed, to get monotonicity. We continue
this process to the right, until finally we reach yn . If skillfully implemented, PAVA has a
computational complexity of O(n); see [117] for a comparison of various algorithms to solve
isotonic regression in `p -metrics, for p ≥ 1. The isoreg command in the stats package in the
R programming language implements the isotonic LSE. Further, see [116] for an efficient
(requiring only O(n) time) computation of the unimodal LSE (Example 1.4).
Order preserving regression on a partially ordered set. Given a partial order on the design points xi ’s we can compute the LSE θ̂ of the isotonic (order preserving) θ∗
by solving (7). Here C, the space where Y is projected onto to obtain θ̂, is a closed convex
cone and can be represented as
C := {(θ1 , . . . , θn ) ∈ Rd : θi ≤ θj if xi - xj , for some i 6= j}.
Thus, the computation of θ̂ involves solving a quadratic program with O(n2 ) linear constraints (although for some special situations, like isotonic regression in d = 1, C can be
represented by O(n) linear constraints). The computation of the order preserving LSE on a
partially ordered set (Example 1.2) has received quite a bit of attention recently; see e.g., [80,
118] and the references therein. In particular, https://github.com/sachdevasushant/Isotonic
gives an implementation of the isotonic LSE using interior point methods. The special case
of the matrix isotonic LSE defined in (22) can be computed efficiently by an iterative algorithm (see e.g., [47] and [106, Chapter 1]).
Once we obtain θ̂ = (θ̂1 , . . . , θ̂n ) we can then easily construct an estimate fˆn of the order
preserving f at any x (not necessarily a design point) by taking a maximum over a selected
number of coordinates of θ̂: We can define fˆn as
fˆn (x) := sup θ̂j ,

for x ∈ Rd ,

j:xj - x

where we take the convention that sup(∅) = −∞. Note that fˆn is indeed order preserving
— for u, v with u - v we have fˆn (u) = supj:xj - u θ̂j ≤ supj:xj - v θ̂j = fˆn (v), as in the right
side the supremum is taken over a bigger set.
Convex regression. Algorithms for the computation of the convex LSE (Example 1.3)
when d = 1 can be found in [41], [44], [59] and the references therein. When d > 1, the
problem is substantially harder: Due to the lack of a natural ordering of points in Rd (for
d > 1), the constraint set C is not easy to express (cf. (6)). In fact, in this case C can be
expressed as the projection of the higher-dimensional polyhedron
n
o
(ξ, θ) ∈ Rdn+n : ξ = [ξ1> , . . . , ξn> ]> , θj + ξj> (xi − xj ) ≤ θi , ∀ i, j = 1, . . . , n
onto the space of θ; see [31]. The above characterization can be seen as a consequence
of the subgradient inequality for convex functions; see [107, Theorem 25.1, p. 242]. Thus
imsart-generic ver. 2014/10/16 file: ShapeRegArxiv.tex date: September 19, 2017

A. Guntuboyina and B. Sen/Shape-restricted Regression

29

the computation of the convex LSE θ̂ involves solving a quadratic program with n(d + 1)
variables and n(n−1) linear constraints; see [109] for the characterization, computation and
consistency of the convex LSE where off-the-shelf interior point solvers (e.g., cvx, MOSEK,
etc.) were used to compute θ̂. However, these off-the-shelf solvers do not scale well and
become prohibitively expensive for n ≥ 300 mainly due to the presence of O(n2 ) linear
constraints. In [91], exploiting problem specific structure, the authors propose a scalable
algorithmic framework based on the augmented Lagrangian method to compute the convex
LSE θ̂.
Shape constrained additive models. The computation of the additive shape-restricted
(Example 1.5) LSE is discussed in [96, 32]. As this reduces to solving a quadratic program
with O(n) linear constraints, off-the-shelf solvers can be effectively used for computing the
LSE. The shape-restricted LSE can also be computed efficiently by the back-fitting algorithm ([19]) — a simple iterative procedure used to fit a generalized additive model —
which involves fitting one function at a time (out of the d many univariate nonparametric
functions).
Shape-restricted single index model. Let us now look at Example 1.6. Here interest focuses on estimating the nonparametric (shape-restricted) function m and the finitedimensional parameter β ∗ . Although single index models are well-studied in the statistical
literature (see e.g., [103], [82], [34] and the references therein), estimation and inference in
shape-restricted single index models are not very well-developed, despite their numerous
applications. The LSE in the monotone single index model is defined as
(m̂, β̂) := arg min
ψ,β

n
X
2
(yi − ψ(x>
i β)) ,

(44)

i=1

where the minimization is over all nondecreasing functions ψ : R → R and over β ∈ Rd
(with kβk = 1, for identifiability). As the above LSE solves a non-convex problem, its
computation is non-trivial. A version of the following alternating minimization scheme is
typically applied to compute the LSE. For a fixed β, the sum-of-squared errors can be easily
minimized over all nondecreasing functions (as this reduces to the problem to univariate
isotonic regression). However, the minimization of the profiled least squares criterion (over
β), for a fixed ψ, is non-smooth and non-convex; see [73, 53] and the references therein
for some strategies to find β̂. A similar strategy is employed in computing the maximum
likelihood estimate in the related problem of current status regression; see e.g., [53]. Also
see [77] and the R package simest for the related computation in a convex single index
model.
6. Some open problems
In this section we state and motivate a few open research problems and some possible future
directions.
Beyond i.i.d. Gaussian errors. We have mentioned in Section 2 that risk bounds for
isotonic regression do not assume that the errors ε1 , . . . , εn are Gaussian. Indeed, the worst
case risk bound (12) as well as the adaptive risk bound in Theorem 2.1 work under no
imsart-generic ver. 2014/10/16 file: ShapeRegArxiv.tex date: September 19, 2017

A. Guntuboyina and B. Sen/Shape-restricted Regression

30

distributional assumptions on the errors (it is only assumed that the errors are i.i.d. with
mean zero and finite variance σ 2 ; even the i.i.d. assumption can be relaxed considerably; see
[129]). However the risk bounds for other shape-restricted regression problems (including
convex regression, isotonic regression on partially ordered sets such as multivariate isotonic
regression, unimodal regression, etc.) assume Gaussian (or sub-Gaussian) errors. Based on
the results for the univariate isotonic LSE, we believe that the assumption of Gaussianity
should not really be necessary for these other problems as well. However the existing proof
techniques for these risk bounds strongly rely on the assumption of sub-Gaussianity. It will
be very interesting to prove risk bounds in these problems without Gaussianity. We believe
that new techniques will need to be developed for this.
Beyond `2 -loss. Most of the risk results available in the shape constrained literature
apply only to the `2 -loss. A notable exception is the case of isotonic regression where risk
bounds are available under the `p -loss for every p ≥ 1 (as already described in Subsection
2.2.1). It will be interesting to develop risk results for `p -losses in problems such as convex
regression and multivariate shape-restricted regression. The results for isotonic regression
(see Theorem 2.2 and the following discussion) indicate that adaptation risk bounds for
LSEs have a different relationship with oracle risk bounds for p 6= 2. For example, for p < 2,
the isotonic LSE is suboptimal only by a constant factor in comparison to the oracle while
for p > 2, the isotonic LSE is significantly suboptimal. We believe that it is quite non-trivial
to study risk of the LSEs under `p -loss functions for p 6= 2. The existing abstract theory for
studying LSEs seems to give risk results only under the `2 -loss function.
Minimax Results. The risk of a LSE R(θ̂, θ∗ ) in a shape-restricted regression problem
usually varies quite significantly as θ∗ varies over the parameter space. For example, in
isotonic regression, the risk behaves as n−2/3 when V (θ∗ ) is bounded and as (k/n) log(n/k)
when θ∗ is piecewise constant having k constant pieces. Minimax lower bounds over these
parameter classes allow the assessment of optimality of the LSE compared to other estimators. We mentioned in Remark 2.3 that the isotonic LSE is minimax optimal over
{θ ∈ I : θn − θ1 ≤ V } for a wide range of values of V . In an interesting recent paper [46],
the authors characterized the minimax risk over the class of all monotone vectors with at
most k constant pieces. Their results imply that the risk (k/n)(log(n/k)) achieved by the
isotonic LSE over this class is only suboptimal by a factor of (log n)/(log log n) in comparison to the minimax risk. Minimax lower bounds exist for other shape-restricted regression
problems (see e.g., [64, 25, 16, 65, 28]) which suggest that the LSE is nearly minimax optimal
but some of these results are not as tight as the corresponding results for univariate isotonic
regression. It will be interesting to develop tight minimax results for other shape-restricted
regression problems which will allow a precise evaluation of the minimaxity properties of
the LSEs.
Estimation of other shape constrained regression functions. In the recent years
there has been quite a bit of interest in studying different shape-restricted regression functions, beyond d = 1. We have already seen a few such examples in this paper (e.g., Examples 1.2, 1.3 and 1.5). What are other useful shape-restrictions in multi-dimension? In the
following we mention a few such shape constraints (that have many real applications): (i)
unordered weak majorization, (ii) quasiconvexity, and (iii) supermodularity.

imsart-generic ver. 2014/10/16 file: ShapeRegArxiv.tex date: September 19, 2017

A. Guntuboyina and B. Sen/Shape-restricted Regression

31

Unordered weak majorization. In Example 1.2 we discussed the problem of estimating an
order preserving regression function (with respect to a partial order). Robertson et al. [106,
Chapter 1] gives a nice overview of the properties of the LSE in this problem. As an example
we introduced a generalization of monotonicity beyond d = 1, namely, coordinate-wise
monotonicity. In the following we introduce and characterize another related (and slightly
stronger) notion of monotonicity in multi-dimensions that is closely tied to the concept of
majorization and Schur-convexity (see e.g., [90]). We define the unordered weak majorization
partial order - as
(u1 , . . . , ud ) - (v1 , . . . , vd ) if and only if

i
X
k=1

uk ≤

i
X

vk for i = 1, 2, . . . , d.

k=1

The following result characterizes all functions that preserve the ordering -.
Theorem 6.1. Let f : Rd → R be a continuously differentiable function. Then f preserves
the partial order - if and only if for any z ∈ Rd ,
f(1) (z) ≥ f(2) (z) ≥ . . . ≥ f(d) (z) ≥ 0,

(45)

where f(i) denotes the partial derivate of f with respect to the i’th coordinate.
The above result shows that in a regression setup if f can be assumed to obey (45), i.e.,
the influence of the predictor variables is ordered, LS estimation under the unordered weak
majorization partial order can be used to estimate f . Constraints like (45) appear quite
often in econometrics; see e.g., [128] and the references therein.
Quasiconvexity. A function f : Rd → R is quasiconvex if and only if its sub-level sets
Sα (f ) = {u ∈ Rd : f (u) ≤ α} are convex for every α ∈ R; see [18, Section 3.4]. Alternatively,
a function f is quasiconvex if and only if f (αu + (1 − α)v) ≤ max{f (u), f (v)}, for all
u, v ∈ Rn , and α ∈ [0, 1] (cf. (5)). Quasiconvex functions extend the notion of unimodality to
multi-dimensions and have applications in mathematical optimization and economics. The
computation of f using the method of least squares is likely to be a non-convex problem.
Supermodularity. A function f : Rd → R is supermodular if f (u ∨ v) + f (u ∧ v) ≥ f (u) +
f (v) for all u, v ∈ Rd , where u ∨ v and u ∧ v denote the component-wise maximum and minimum of u and v, respectively, i.e., (u1 , . . . , ud )∨(v1 , . . . , vd ) := (max{u1 , v1 }, . . . , max{ud , vd })
and (u1 , . . . , ud ) ∧ (v1 , . . . , vd ) := (min{u1 , v1 }, . . . , min{ud , vd }). The concept of supermodularity is used in the social sciences (economics and game theory). If f is twice continuously
∂ 2 f (u)
differentiable, then supermodularity is equivalent to the condition ∂u
≥ 0 for all i 6= j;
i ∂uj
see e.g., [114].
In all the above problems, computation of the LSE, its theoretical properties (consistency,
rates of convergence, etc.) are unknown.
Connection to nonnegative least squares. In many shape-restricted regression problems the LSE θ̂ is defined as the projection of Y onto a closed convex polyhedral cone C;
see Examples 1.1–1.3 and 1.5. As every closed convex polyhedral cone C in Rn can be represented in terms of its generators (i.e., there exists a finite subset of C, whose elements
are referred to as generators of C, such that every vector θ ∈ C is a nonnegative linear
combination of the generators; see, for example, [108, Corollary 7.1a]), θ̂ can be thought
imsart-generic ver. 2014/10/16 file: ShapeRegArxiv.tex date: September 19, 2017

A. Guntuboyina and B. Sen/Shape-restricted Regression

32

of as solving a nonnegative least squares problem. Moreover, in the examples mentioned
above, C is generated by at least O(n) vectors (which form the design matrix) that are
highly correlated; see e.g., [93] for the exact form of the generators of some of the examples discussed. In fact, for isotonic and convex regression in d = 1 there are exactly n + 1
and n + 2 generators, respectively. However, for their higher dimensional analogues (i.e.,
d > 1) it is not clear what the generators are. We think this is an open problem. More
generally, one can ask how does one construct the design matrix (or the generators) corresponding to any closed convex polyhedral cone expressed in terms of linear inequalities,
e.g., C := {θ ∈ Rn : Aθ ≤ 0} where A is an m × n matrix (m being the number of linear
constraints) and the ‘≤’ is interpreted coordinate-wise.
As the number of generators (or the columns of the design matrix) is increasing with n, we
are essentially solving a ‘high-dimensional’ nonnegative least squares problem; see e.g., [115].
A general open question is: Can a theory be developed on the estimation accuracy of the
LSE θ̂ based solely on the properties of the design matrix? It may be noted here that the
generators can be highly correlated.
Boundary behavior of shape-restricted LSEs. It is well known that the isotonic
LSE is inconsistent at the boundary of the covariate domain, i.e., fˆn (0+) does not consistently estimate f (0+) (see e.g., [125, 78, 9] for detailed discussions on the properties of the
LSE for a nonincreasing density near 0). Intuitively, this inconsistency is because there are
very few ‘constraints’ near the boundary (of the covariates). This phenomenon is expected to
persist for other shape constrained LSEs, especially in multi-dimensional problems. However
not much is known about the boundary behavior of these LSEs. Even in one-dimensional
convex regression, as far as we are aware, whether fˆn (0+) is a Op (1) random variable is
not known; see [48] for some results on fˆn (0+) and its derivative (also see [7]). This has
motivated the study of bounded/penalized shape-restricted LSEs; see e.g., [31, 127, 77].
Shape-restricted single index models. Although several smoothing based methods
have been proposed and investigated in single index models (see e.g., [98, 53] and the
√
references therein) to obtain n-consistent and efficient estimators of β ∗ (see [100] for a
brief overview of the notion of semiparametric efficiency), not much is known for just shaperestricted single index models. Durot et al. [8] studied the LSE in a monotone single index
model (see (44)) and showed the n1/3 -consistency of the LSEs of m and β ∗ ; also see [98]
and [53]. However many open questions remain. The limiting distribution of the LSE β̂ of
√
β ∗ is unknown; in fact, it is not known whether β̂ is n-consistent.
In a convex single index model (i.e., m is convex), [77] shows that the Lipschitz constrained convex LSE (where we minimize the least squares criterion over the class of all
L-Lipschitz convex functions, for L fixed) yields a semiparametrically efficient estimator of
the index parameter β ∗ . However, the behavior of the convex LSE (without the Lipschitz
assumption) is unknown.
As mentioned above, the computation of the shape-restricted LSEs is non-trivial. Usually
an alternating minimization scheme is used to compute the LSEs. However, no convergence
guarantees (to a local optimum) exist for such an alternating minimization procedure.

imsart-generic ver. 2014/10/16 file: ShapeRegArxiv.tex date: September 19, 2017

A. Guntuboyina and B. Sen/Shape-restricted Regression

33

Acknowledgements
We would like to thank Moulinath Banerjee, Sabyasachi Chatterjee and Probal Chaudhuri
for helpful discussions.
7. Proofs
7.1. Proof of Theorem 2.1.
To prove the adaptive risk bound (13) for isotonic regression we shall apply inequality (33)
with Z = ε/σ which reduces the task to proving that
EkΠTI (θ) (Z)k2 ≤ 4k(θ) log

en
k(θ)

for every θ ∈ I.

(46)

To prove this, we obviously need to characterize the tangent cone TI (θ). Fix θ ∈ I and let
k = k(θ). This means that θ has k constant pieces (and P
k − 1 jumps). Let the lengths of
the k constant pieces be n1 , . . . , nk (so that nj ≥ 1 and kj=1 nj = n). It is now an easy
exercise to verify that the tangent cone TI (θ) is given by
TI (θ) = In1 × · · · × Ink

(47)

where × denotes Cartesian product and, for an integer ` ≥ 1,
n
o
I` := (x1 , . . . , x` ) ∈ Rl : x1 ≤ · · · ≤ x` .
The equality (47) immediately implies that
2

EkΠTI (θ) (Z)k =

k
X

EkΠInj (Z (j) )k2

(48)

j=1

where Z (j) denotes the part of the vector Z corresponding to the j th constant piece of θ.
The key now is to prove that


1
1
2
for every n ≥ 1.
(49)
EkΠIn (Z)k ≤ 4 1 + + · · · +
2
n
P
This, along with the standard fact that nj=1 (1/j) ≤ log(en), implies that
EkΠInj (Z (j) )k2 ≤ 4 log(enj ).
Thus, we obtain, via (48),
EkΠTI (θ) (Z)k2 ≤

k
X
j=1

4 log(enj ) ≤ 4k log

en
k

imsart-generic ver. 2014/10/16 file: ShapeRegArxiv.tex date: September 19, 2017

A. Guntuboyina and B. Sen/Shape-restricted Regression

34

where the last inequality follows from the concavity of the logarithm function. This proves
(46) and consequently (13). To prove (49) let
U = (U1 , . . . , Un ) := ΠIn (Z).
The representation (11) gives
for every 1 ≤ j ≤ n

Uj = min max Z̄uv
v≥j u≤j

where Z̄uv represents the mean of Zu , . . . , Zv . It is clear from here that
Uj ≤ max Z̄un
u≤j

so that
E(Uj )2+ ≤ E max(Z̄un )2+ .
u≤j

Here x+ := max(x, 0) and a2+ := (a+ )2 . The key now is to realize that the sequence of random variables Z̄1n , . . . , Z̄jn is a martingale (easy to verify) so that by Doob’s submartingale
maximal inequality, we deduce that
E(Uj )2+ ≤ 4E(Z̄jn )2+ .
By an analogous argument, we can also deduce that
E(Uj )2− ≤ 4E(Z̄1j )2−
where x− := max(0, −x). Putting together, we deduce
EkΠIn (Z)k2 =

n
X

EUj2 ≤ 4

j=1

n
X


E(Z̄jn )2+ + E(Z̄1j )2− .

j=1

Now because Z1 , . . . , Zn are i.i.d. (Z1 , . . . , Zn ) has the same distribution as (Zn , . . . , Z1 ) so
that
EkΠIn (Z)k2 ≤ 4

n
X

n
X


E(Z̄1,n−j+1 )2+ + E(Z̄n−j+1,n )2− = 4
E(Z̄1j )2+ + E(Z̄jn )2− .

j=1

j=1

Adding the above two inequalities, we obtain
2

2EkΠIn (Z)k ≤ 4

n
X

2

E(Z̄jn ) + E(Z̄1j )

2



=4

j=1

n 
X
j=1

1
1
+
n−j+1 j


=8

n
X
1
j=1

j

.

This proves (49) and completes the proof of (13).

imsart-generic ver. 2014/10/16 file: ShapeRegArxiv.tex date: September 19, 2017

A. Guntuboyina and B. Sen/Shape-restricted Regression

35

7.2. Proof of inequality (12).
The initial steps of this proof follow the same steps as in the beginning of the proof of
Theorem 2.2 (the proof of Theorem 2.2 is provided in Subsection 7.4). Specifically, we use
inequality (55) for p = 2 which says that

 2


2

∗
∗
∗
¯
Eθ∗ θ̂j − θj
≤ Eθ∗ θ j,j+m − θj + max ε̄u,j+m
.
u≤j

+

+

√
for every 1 ≤ j ≤ n and 0 ≤ m ≤ n−j. The elementary inequality E(c+X)2 ≤ (c+ EX 2 )2
for c ≥ 0 applied to c = θ¯∗ j,j+m − θj∗ ≥ 0 and X := (maxu≤j ε̄u,j+m )+ now gives


Eθ∗ θ̂j − θj∗

2
+

≤

!2
r
2
θ¯∗ j,j+m − θj∗ + E max (ε̄u,j+m )+ .
u≤j

Doob’s maximal submartingale inequality applied to the expectation above results in
E max (ε̄u,j+m )2+ ≤ 4E (ε̄u,j+m )2+ ≤ 4 (ε̄u,j+m )2 =
u≤j

4σ 2
.
m+1

We have thus proved that


Eθ∗ θ̂j −

θj∗

2
+


2
2σ
∗
∗
¯
≤ min
θ j,j+m − θj + √
.
0≤m≤n−j
m+1

(50)

2

2σ
∗
∗
¯
.
≤ min
θj − θ j−m,j + √
0≤m≤j−1
m+1

(51)

Similarly, we can prove that


Eθ∗ θ̂j −

θj∗

2
−

For each 1 ≤ j ≤ n, let us define


2σ
∗
∗
¯
m1 (j) := max 0 ≤ m ≤ n − j : θ j,j+m − θj ≤ √
m+1
and



2σ
∗
∗
¯
m2 (j) := max 0 ≤ m ≤ j − 1 : θj − θ j−m,j ≤ √
.
m+1

Note that these choices of m1 (j) and m2 (j) are different from those made in the proof of
Theorem 2.2. From (50) and (51), we obtain


2
1
1
Eθ∗ θ̂j − θj∗ ≤ 16σ 2
+
m1 (j) + 1 m2 (j) + 1
so that the risk of θ̂ is bounded by


n
n
2
X
X
16σ
1
1

.
R(θ̂, θ∗ ) ≤
+
n
m1 (j) + 1
m2 (j) + 1
j=1

(52)

j=1

imsart-generic ver. 2014/10/16 file: ShapeRegArxiv.tex date: September 19, 2017

A. Guntuboyina and B. Sen/Shape-restricted Regression

36

We shall bound j (m1 (j) + 1)−1 below. The bound for the term involving m2 (j) will be
similar.
For each m ≥ 0, let
P

ρ(m) :=

n
X

I {m1 (j) = m}

and

l(m) :=

j=1

n
X

I {m1 (j) < m} .

j=1

Clearly ρ(m) = l(m + 1) − l(m) and hence
n
X
j=1

X ρ(m)
X l(m + 1) − l(m)
X
l(m + 1)
1
=
=
=
m1 (j) + 1
m+1
m+1
(m + 1)(m + 2)
m≥0

m≥0

(53)

m≥0

because l(0) = 0. For 0 ≤ m ≤ n − 1, we can write
l(m + 1) =

n
X

I{m1 (j) < m + 1} ≤

j=1

n−m−1
X

I{m1 (j) < m + 1} + m + 1.

j=1

When 1 ≤ j ≤ n − m − 1 and m1 (j) < m + 1, the definition of m1 (j) implies that
2σ
θ¯∗ j,j+m+1 − θj∗ > √
m+2
so that
I {m1 (j) < m + 1} <

√

m+2

θ¯∗ j,j+m+1 − θj∗
2σ

for 1 ≤ j ≤ n − m − 1.

and
√
l(m + 1) ≤ m + 1 +

n−m−1

m + 2 X ¯∗
θ j,j+m+1 − θj∗
2σ

for 0 ≤ m ≤ n − 1.

j=1

We now use Lemma 7.1 to bound the right hand side above in terms of V := V (θ∗ ). Indeed,
Lemma 7.1 (applied with l = m + 2) gives
√
(m + 1) m + 2
V.
l(m + 1) ≤ m + 1 +
2σ
Also, the trivial upper bound l(m + 1) ≤ n always holds. We thus obtain
√


(m + 1) m + 2
l(m + 1) ≤ min n, m + 1 +
V
for all m ≥ 0.
2σ
This inequality, together with (53), gives
n
X
j=1

√


X
1
1
(m + 1) m + 2
≤
min n, m + 1 +
V .
m1 (j) + 1
(m + 1)(m + 2)
2σ
m≥0

imsart-generic ver. 2014/10/16 file: ShapeRegArxiv.tex date: September 19, 2017

A. Guntuboyina and B. Sen/Shape-restricted Regression

Now let
s := min

 
σn 2/3
V

37


,n

and write (C below stands for a positive constant whose value may change from appearance
to appearance; it does not depend on n, V or σ)
n
X
j=1

√


X
1
(m + 1) m + 2
1
≤
min n, m + 1 +
V
m1 (j) + 1
(m + 1)(m + 2)
2σ
m≥0
√
X
X (m + 1) + (m + 1) m + 2(V /(2σ))
1
+n
≤
(m + 1)(m + 2)
(m + 1)(m + 2)
m>s
0≤m≤s

≤

X
m≤n

1
+
m+2

X
m≤(σn/V )2/3

X
1
V
√
+n
(m + 1)(m + 2)
2σ m + 2
m>s

V  σn 1/3 Cn
+
σ V
s
 2/3 !
2/3
1/3
1
V n
V
,
+ Cn max
= C log(en) + C
2/3
σn
n
σ
≤ C log(en) + C

≤ C log(en) + C

V 2/3 n1/3
.
σ 2/3

One can analogously argue that
n
X
j=1

1
V 2/3 n1/3
≤ C log(en) + C
.
m2 (j) + 1
σ 2/3

The above pair of inequalities, when combined with (52), complete the proof of inequality
(12).
The following auxiliary lemma was used in the proof of inequality (12).
Lemma 7.1. For every nondecreasing sequence a1 , . . . , an and l ≥ 1, we have
n−l+1
X

(āi,i+l−1 − ai ) ≤

i=1

l−1
(an − a1 )
2

(54)

where āi,j for i ≤ j is the mean of ai , . . . , aj .
Proof. Write bi = ai − ai−1 for i ≥ 2 and b1 = a1 . It is elementary to check that
l−1

āi,i+l−1

1X
=
(l − j)bi+j
l

for i = 1, . . . , n − l + 1.

j=1

As a result, the right hand side of (54) can be easily seen to be a linear combination of
b2 , . . . , bn which are all nonnegative because a1 , . . . , an is nondecreasing. The maximum
value of the coefficient of any bj is (l − 1)/2 which completes the proof.

imsart-generic ver. 2014/10/16 file: ShapeRegArxiv.tex date: September 19, 2017

A. Guntuboyina and B. Sen/Shape-restricted Regression

38

7.3. Proof of (17).
It is now assumed that ε1 , . . . , εn are i.i.d. N (0, σ 2 ). Let Zi := εi /σ, i = 1, . . . , n so that
Z1 , . . . , Zn are i.i.d. N (0, 1). We will use the same ideas as in the proof of Theorem 2.1.
However, when Z ∼ Nn (0, In ), (49) holds with equality and without the multiplicative
factor 4, for every n ≥ 1. This can be proved via symmetry arguments based on the theory
of finite reflection groups. A sketch of this argument can be found in [2, Subsection D.4].
This observation therefore completes the proof of (17).
7.4. Proof of Theorem 2.2
We only need to prove the first inequality in (18). The second inequality is a consequence
of Jensen’s inequality (note that x 7→ x(2−p)+ /2 is concave on (0, ∞)).
We use the representation (11) of the isotonic LSE. Fix 1 ≤ j ≤ n and let 0 ≤ m ≤ n − j.
By (11), we have

θ̂j = min max Ȳuv ≤ max Ȳu,j+m = max θ¯∗ u,j+m + ε̄u,j+m .
u≤j

v≥j u≤j

u≤j

Because θ∗ ∈ I, we have θ¯∗ u,j+m ≤ θ¯∗ j,j+m for all u ≤ j. We therefore have

θ̂j − θj∗ ≤ θ¯∗ j,j+m − θj∗ + max ε̄u,j+m .
u≤j

Taking positive parts on both sides and then raising to power p and taking expectations on
both sides, we derive

p

p

∗
∗
∗
¯
Eθ∗ θ̂j − θj
≤ Eθ∗ θ j,j+m − θj + max ε̄u,j+m
.
(55)
u≤j

+

+

Let us now introduce some notation. Let k = k(θ∗ ) and let n1 , . . . , nk denote the lengths of
the k constant blocks of θ∗ . Let s0 := 0 and let si := n1 + · · · + ni for 1 ≤ i ≤ k. For each
j = 1, . . . , n, let us define two integers m1 (j) and m2 (j) in the following way: m1 (j) = si − j
and m2 (j) = j − 1 − si−1 when si−1 + 1 ≤ j ≤ si . The key is to realize that θ¯∗ j,j+m1 (j) = θj∗
for every j. As a result, inequality (55) with m = m1 (j) gives

p

p

p
p
∗
Eθ∗ θ̂j − θj
≤ E max ε̄u,j+m1 (j)
= σ E max Z̄u,j+m1 (j)
(56)
+

u≤j

+

u≤j

+

where, as before, Zj := εj /σ for j = 1, . . . , n. The fact that Z̄1,j+m1 (j) , . . . , Z̄j,j+m1 (j) is
a martingale allows us to deduce, via Doob’s Lp maximal inequality for nonnegative submartingales, that

p

p
p
p
E max Z̄u,j+m1 (j)
≤
E Z̄j,j+m1 (j) +
(57)
u≤j
p−1
+


p
p/2
p
1
p
=
E (η)+
p−1
m1 (j) + 1
imsart-generic ver. 2014/10/16 file: ShapeRegArxiv.tex date: September 19, 2017

A. Guntuboyina and B. Sen/Shape-restricted Regression

39

where η ∼ N (0, 1) is a standard normal random variable. This inequality requires p > 1 so
we assume now that p > 1. The argument for controlling the left hand side above for p = 1
will be given subsequently.
We have therefore proved that, for p > 1,


Eθ∗ θ̂j −

θj∗

p
+

≤σ

p



p
p−1

p

E (η)p+



1
m1 (j) + 1

p/2

for every 1 ≤ j ≤ n. A similar argument gives


Eθ∗ θ̂j −

θj∗

p
−

≤σ

p



p
p−1

p

E (η)p−



1
m2 (j) + 1

p/2
.

Putting the above two inequalities together, we obtain
(

p
p/2 
p/2 )
p
p
1
1
σp
∗
p
Eθ∗ θ̂j − θj ≤
E|η|
+
2 p−1
m1 (j) + 1
m2 (j) + 1
for every j = 1, . . . , n. Note now that
n 
X
j=1

1
m1 (j) + 1

p/2
=

k
X

P

j (m2 (j)

n

R



i=1 j=si−1 +1

and the same bound holds for
(p)

si
X

X
1
(θ̂, θ ) = Eθ∗
θ̂j − θj∗
n
∗

1
si − j + 1

p/2
=

ni  p/2
k X
X
1
i=1 j=1

j

+ 1)−p/2 as well. This gives
p

i=1

σp
≤
n



p
p−1

p

p

E|η|

ni  p/2
k X
X
1
i=1 j=1

j

.

(58)

We now consider the two cases 1 < p < 2 and 2 < p < ∞ separately. The simple bound
n
X

j −p/2 ≤

j=1

2
n1−(p/2)
2−p

for 1 ≤ p < 2

gives
R

(p)

2 σp
(θ̂, θ ) ≤
2−p n
∗



p
p−1

p

p

E|η|

k
X

1−(p/2)

ni

i=1

which proves (18) for 1 < p < 2.
P
P P i  1 p/2
−p/2 in (58). This completes the
For p > 2, we simply use ki=1 nj=1
≤k ∞
j=1 j
j
proof of (18) for the case when p > 1, p 6= 2.
For p = 1, we shall bound the expectation in the right hand side of (56) in the following
way. Let τ 2 := 1/(m1 (j) + 1) be the variance of Z̄j,j+m1 (j) .


Z
E max Z̄u,j+m1 (j)
≤τ+
u≤j

+

τ

∞


P max Z̄u,j+m1 (j)
u≤j


+


≥ t dt

imsart-generic ver. 2014/10/16 file: ShapeRegArxiv.tex date: September 19, 2017

A. Guntuboyina and B. Sen/Shape-restricted Regression

40

Because Z̄u,j+m1 (j) , u = 1, . . . , j is a martingale, the sequence (Z̄u,j+m1 (j) )+ is a nonnegative
submartingale and hence Doob’s maximal inequality gives


oi

 n

1 h
P max Z̄u,j+m1 (j) + ≥ t ≤ E Z̄u,j+m1 (j) + Z̄u,j+m1 (j) + ≥ t
u≤j
t
so that


Z ∞n
o dt


E max Z̄u,j+m1 (j)
Z̄u,j+m1 (j) + ≥ t
≤ τ + E Z̄u,j+m1 (j) +
u≤j
t
τ
+

 Z ∞
Z ∞
dt
dt
2
{η+ ≥ t}
≤ τ + τ Eη+
= τ + τ E η+
t
t2
1
1
3τ
3
≤
=
2
2(m1 (j) + 1)
where, as before, η is a standard normal random variable. Using the above inequality in
place of (57) proves (18) for p = 1 thereby completing the proof of Theorem 2.2.
7.5. Proof of Lemma 3.1
Assumption (27) becomes
n1
n2
1 X
1 X
∗
f1 (i1 /n1 ) =
f2∗ (i2 /n2 ) = 0.
n1
n2
i1 =1

i2 =1

It is easy to see then that for every pair of functions f1 and f2 satisfying
n1
n2
1 X
1 X
f1 (i1 /n1 ) =
f2 (i2 /n2 ) = 0,
n1
n2
i1 =1

(59)

i2 =1

and µ ∈ R, the following quantity
n1 X
n2
X

yi1 ,i2 − µ − f1 (i1 /n1 ) − f2 (i2 /n2 )

2

i1 =1 i2 =1

is equal to
n1 X
n2
X

(yi1 ,i2 − ȳi1 ,· − ȳ·,i2 + ȳ·,· )2 + n2

i1 =1 i2 =1

n1
X

(ȳi1 ,· − ȳ·,· − f1 (i1 /n1 ))2

i1 =1

+ n1

n2
X

(ȳ·,i2 − ȳ·,· − f2 (i2 /n2 ))2 + n1 n2 (ȳ·,· − µ)2

i2 =1

where
ȳi1 ,·

n2
1 X
yi1 ,i2
:=
n2
i2 =1

and

ȳ·,i2

n1
1 X
:=
yi1 ,i2
n1
i1 =1

and

n2
n1 X
1 1 X
ȳ·,· :=
yi1 ,i2 .
n1 n2
i1 =1 i2 =1

imsart-generic ver. 2014/10/16 file: ShapeRegArxiv.tex date: September 19, 2017

A. Guntuboyina and B. Sen/Shape-restricted Regression

41

This nice decomposition of the least squares criterion implies that
(n
)
n1
1
X
X
2
fˆ1 = arg min
(ȳi ,· − ȳ·,· − f1 (i1 /n1 )) : f1 ∈ F1 ,
f1 (i1 /n1 ) = 0
1

i1 =1

and

(
fˆ2 = arg min

n2
X

i1 =1

2

(ȳ·,i2 − ȳ·,· − f2 (i2 /n2 )) : f2 ∈ F2 ,

n2
X

)
f2 (i2 /n2 ) = 0

i2 =1

i2 =1

and also that µ̂ = ȳ·,· .
The above decomposition of the least squares criterion holds for every µ, f1 and f2
satisfying (59). In particular, it holds when f1 (resp. f2 ) is replaced by f1∗ (resp. f2∗ ). It
follows therefore that fˆ1 = fˆ1OR and fˆ2 = fˆ2OR .
7.6. Proof of Theorem 6.1
The proof of the above result uses the following lemma, which we prove in Section 7.7.
Lemma 7.2. Let f be a continuous real-valued function defined on Rd (d ≥ 1). Define the
transformation Tk, : Rn → Rn for k ∈ {1, 2, . . . , d} and  ∈ R as
Tk, (z1 , . . . , zd ) := (z1 , . . . , zk−1 , zk + , zk+1 − , zk+2 , . . . , zd ).
Then, f is order preserving with respect to the partial order - if and only if for every z ∈ Rd ,
and k = 1, . . . , d, f (Tk, (z)) is nondecreasing in .
Fix z ∈ Rd . Since f is differentiable and order preserving, by Lemma 7.2 this is equivalent
to
d
f (Tk, (z)) ≥ 0,
d
i.e., f(k) (Tk, (z)) − f(k+1) (Tk, (z)) ≥ 0,
for all k ∈ {1, 2, . . . , d − 1}. Letting  go to zero yields the desired result.
7.7. Proof of Lemma 7.2
Suppose that f is order preserving, i.e., if u - v, where u, v ∈ Rd , then f (u) ≤ f (v). Let
z ∈ Rd and  < 0 ∈ R. We want to show that f (Tk, (z)) ≤ f (Tk,0 (z)), for k = 1, . . . , d. Fix
k ∈ {1, 2, . . . , d} and define
u := Tk, (z)

and

v := Tk,0 (z).
P
Pj
P
Pk
P
Therefore, ji=1 ui =
vi , for j = 1, . . . , k − 1, ki=1 ui <
vi , and ji=1 ui =
i=1
i=1
Pj
i=1 vi , for j = k + 1, . . . , d. Thus, u - v and we have f (u) ≤ f (v), i.e., f (Tk, (z)) ≤
f (Tk,0 (z)).
Now suppose that f (Tk, (z)) is nondecreasing in , for every z ∈ Rd . Take u ≡ (u1 , . . . , ud ) P
(v1 , . . . , vd ) ≡ v in Rd . Let k := ki=1 (vi − ui ), for k = 1, . . . , d. By definition, k ≥ 0, for
all k. Define
sk := Tk,k (sk−1 ),
imsart-generic ver. 2014/10/16 file: ShapeRegArxiv.tex date: September 19, 2017

A. Guntuboyina and B. Sen/Shape-restricted Regression

42

for k = 1, . . . , d, where s0 ≡ u. Note that sk ∈ Rd for every k = 1, . . . , d and sd =
Td,d (sd−1 ) = v. Then from the nondecreasing property of f (Td, (·)) (and using the fact
that k ≥ 0 for all k),
f (u) ≡ f (T1,0 (s0 )) ≤ f (T1,1 (s0 )) = f (T2,0 (s1 )) ≤ f (T2,2 (s1 )) = f (T3,0 (s3 ))
≤ · · · ≤ f (Td,0 (sd−1 )) ≤ f (Td,d (sd−1 )) ≡ f (v),
yielding the desired result.
References
[1] J. Abrevaya and J. Huang. On the bootstrap of the maximum score estimator. Econometrica, 73(4):1175–1204, 2005.
[2] D. Amelunxen, M. Lotz, M. B. McCoy, and J. A. Tropp. Living on the edge: phase
transitions in convex programs with random data. Inf. Inference, 3(3):224–294, 2014.
[3] D. Anevski and O. Hössjer. A general asymptotic scheme for inference under order
restrictions. Ann. Statist., 34(4):1874–1930, 2006.
[4] M. Ayer, H. D. Brunk, G. M. Ewing, W. T. Reid, and S. E. An empirical distribution
function for sampling with incomplete information. Ann. Math. Statist., 26:641–647,
1955.
[5] P. Bacchetti. Additive isotonic models. Journal of the American Statistical Association, 84(405):289–294, 1989.
[6] P. Bagchi, M. Banerjee, and S. A. Stoev. Inference for monotone functions under
short- and long-range dependence: confidence intervals and new universal limits. J.
Amer. Statist. Assoc., 111(516):1634–1647, 2016.
[7] F. Balabdaoui. Consistent estimation of a convex density at the origin. Math. Methods
Statist., 16(2):77–95, 2007.
[8] F. Balabdaoui, C. Durot, and H. Jankowski. Least squares estimation in the monotone
single index model. arXiv preprint arXiv:1610.06026, 2016.
[9] F. Balabdaoui, H. Jankowski, M. Pavlides, A. Seregin, and J. Wellner. On the Grenander estimator at zero. Statist. Sinica, 21(2):873–899, 2011.
[10] F. Balabdaoui, K. Rufibach, and J. A. Wellner. Limit distribution theory for maximum
likelihood estimation of a log-concave density. Ann. Statist., 37(3):1299–1331, 2009.
[11] M. Banerjee. Likelihood based inference for monotone response models. Ann. Statist.,
35(3):931–956, 2007.
[12] M. Banerjee. Inference in exponential family regression models under certain shape
constraints using inversion based techniques. In Advances in multivariate statistical
methods, volume 4 of Stat. Sci. Interdiscip. Res., pages 249–271. World Sci. Publ.,
Hackensack, NJ, 2009.
[13] M. Banerjee and J. A. Wellner. Likelihood ratio tests for monotone functions. Ann.
Statist., 29(6):1699–1731, 2001.
[14] M. Banerjee and J. A. Wellner. Confidence intervals for current status data. Scand.
J. Statist., 32(3):405–424, 2005.
[15] R. E. Barlow, D. J. Bartholomew, J. M. Bremner, and H. D. Brunk. Statistical
inference under order restrictions. The theory and application of isotonic regression.
imsart-generic ver. 2014/10/16 file: ShapeRegArxiv.tex date: September 19, 2017

A. Guntuboyina and B. Sen/Shape-restricted Regression

[16]
[17]
[18]
[19]

[20]
[21]

[22]
[23]
[24]
[25]
[26]
[27]
[28]
[29]
[30]

[31]

[32]
[33]
[34]
[35]

43

John Wiley & Sons, London-New York-Sydney, 1972. Wiley Series in Probability and
Mathematical Statistics.
P. C. Bellec. Sharp oracle inequalities for least squares estimators in shape restricted
regression. arXiv preprint arXiv:1510.08029, 2015.
D. P. Bertsekas. Convex analysis and optimization. Athena Scientific, Belmont, MA,
2003. With Angelia Nedić and Asuman E. Ozdaglar.
S. Boyd and L. Vandenberghe. Convex optimization. Cambridge University Press,
Cambridge, 2004.
L. Breiman and J. H. Friedman. Estimating optimal transformations for multiple
regression and correlation. J. Amer. Statist. Assoc., 80(391):580–619, 1985. With
discussion and with a reply by the authors.
H. D. Brunk. Maximum likelihood estimates of monotone parameters. Ann. Math.
Statist., 26:607–616, 1955.
H. D. Brunk. Estimation of isotonic regression. In Nonparametric Techniques in
Statistical Inference (Proc. Sympos., Indiana Univ., Bloomington, Ind., 1969), pages
177–197. Cambridge Univ. Press, London, 1970.
C. Carolan and R. Dykstra. Asymptotic behavior of the Grenander estimator at
density flat regions. Canad. J. Statist., 27(3):557–566, 1999.
S. Chatterjee. A new perspective on least squares under convex constraint. Ann.
Statist., 42(6):2340–2381, 2014.
S. Chatterjee. An improved global risk bound in concave regression. Electron. J.
Stat., 10(1):1608–1629, 2016.
S. Chatterjee, A. Guntuboyina, and B. Sen. On risk bounds in isotonic and other
shape restricted regression problems. Ann. Statist., 43(4):1774–1800, 2015.
S. Chatterjee, A. Guntuboyina, and B. Sen. On matrix estimation under monotonicity
constraints. Bernoulli (to appear; arXiv preprint arXiv:1506.03430), 2017.
S. Chatterjee and J. Lafferty. Adaptive risk bounds in unimodal regression. arXiv
preprint arXiv:1512.02956, 2015.
S. Chatterjee and J. Lafferty.
Denoising flows on trees.
arXiv preprint
arXiv:1602.08048, 2017.
S. Chatterjee and S. Mukherjee. On estimation in tournaments and graphs under
monotonicity constraints. arXiv preprint arXiv:1603.04556, 2017.
X. Chen, A. Guntuboyina, and Y. Zhang. A note on the approximate admissibility of regularized estimators in the gaussian sequence model. arXiv preprint
arXiv:1703.00542, 2017.
X. Chen, Q. Lin, and B. Sen. On degrees of freedom of projection estimators with applications to multivariate shape restricted regression. arXiv preprint
arXiv:1509.01877, 2015.
Y. Chen and R. J. Samworth. Generalized additive and index models with shape
constraints. J. R. Stat. Soc. Ser. B. Stat. Methodol., 78(4):729–754, 2016.
H. Chernoff. Estimation of the mode. Ann. Inst. Statist. Math., 16:31–41, 1964.
X. Cui, W. K. Härdle, and L. Zhu. The EFM approach for single-index models. Ann.
Statist., 39(3):1658–1688, 2011.
I. Demetriou and P. Tzitziris. Infant mortality and economic growth: Modeling by
increasing returns and least squares. In Proceedings of the World Congress on Engiimsart-generic ver. 2014/10/16 file: ShapeRegArxiv.tex date: September 19, 2017

A. Guntuboyina and B. Sen/Shape-restricted Regression

44

neering, volume 2, 2017.
[36] D. L. Donoho. Gelfand n-widths and the method of least squares. Technical report
282, 1990.
[37] C. R. Doss and J. A. Wellner. Inference for the mode of a log-concave density. arXiv
preprint arXiv:1611.10348, 2016.
[38] C. R. Doss and J. A. Wellner. Mode-constrained estimation of a log-concave density.
arXiv preprint arXiv:1611.10335, 2016.
[39] L. Dümbgen. Optimal confidence bands for shape-restricted curves. Bernoulli,
9(3):423–449, 2003.
[40] L. Dümbgen, S. Freitag, and G. Jongbloed. Consistency of concave regression with
an application to current-status data. Math. Methods Statist., 13(1):69–81, 2004.
[41] R. L. Dykstra. An algorithm for restricted least squares regression. J. Amer. Statist.
Assoc., 78(384):837–842, 1983.
[42] B. Fang and A. Guntuboyina. On the risk of convex-constrained least squares estimators under misspecification. arXiv preprint arXiv:1706.04276, 2017.
[43] N. Flammarion, C. Mao, and P. Rigollet. Optimal rates of statistical seriation. arXiv
preprint arXiv:1607.02435, 2016.
[44] D. A. S. Fraser and H. Massam. A mixed primal-dual bases algorithm for regression
under inequality constraints. Application to concave regression. Scand. J. Statist.,
16(1):65–74, 1989.
[45] M. Frisén. Unimodal regression. The Statistician, pages 479–485, 1986.
[46] C. Gao, F. Han, and C.-H. Zhang. Minimax risk bounds for piecewise constant models.
arXiv preprint arXiv:1705.06386, 2017.
[47] F. Gebhardt. An algorithm for monotone regression with one or more independent
variables. Biometrika, 57(2):263–271, 1970.
[48] P. Ghosal and B. Sen. On univariate convex regression. Sankhya A, pages 1–39, 2016.
[49] U. Grenander. On the theory of mortality measurement. II. Skand. Aktuarietidskr.,
39:125–153 (1957), 1956.
[50] P. Groeneboom. The concave majorant of Brownian motion. Ann. Probab.,
11(4):1016–1027, 1983.
[51] P. Groeneboom. Estimating a monotone density. In Proceedings of the Berkeley
conference in honor of Jerzy Neyman and Jack Kiefer, Vol. II (Berkeley, Calif., 1983),
Wadsworth Statist./Probab. Ser., pages 539–555, Belmont, CA, 1985. Wadsworth.
[52] P. Groeneboom and K. Hendrickx. Confidence intervals for the current status model.
arXiv preprint arXiv:1611.08299, 2016.
[53] P. Groeneboom and K. Hendrickx. Current status linear regression. ArXiv e-prints
arXiv:1601.00202, 2016.
[54] P. Groeneboom and G. Jongbloed. Isotonic estimation and rates of convergence in
Wicksell’s problem. Ann. Statist., 23(5):1518–1542, 1995.
[55] P. Groeneboom and G. Jongbloed. Nonparametric estimation under shape constraints,
volume 38 of Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press, New York, 2014. Estimators, algorithms and asymptotics.
[56] P. Groeneboom and G. Jongbloed. Nonparametric confidence intervals for monotone
functions. Ann. Statist., 43(5):2019–2054, 2015.
[57] P. Groeneboom, G. Jongbloed, and J. A. Wellner. A canonical process for estimation
imsart-generic ver. 2014/10/16 file: ShapeRegArxiv.tex date: September 19, 2017

A. Guntuboyina and B. Sen/Shape-restricted Regression

[58]
[59]

[60]

[61]
[62]
[63]
[64]
[65]
[66]
[67]
[68]
[69]
[70]
[71]

[72]
[73]

[74]
[75]

[76]

45

of convex functions: the “invelope” of integrated Brownian motion +t4 . Ann. Statist.,
29(6):1620–1652, 2001.
P. Groeneboom, G. Jongbloed, and J. A. Wellner. Estimation of a convex function:
characterizations and asymptotic theory. Ann. Statist., 29(6):1653–1698, 2001.
P. Groeneboom, G. Jongbloed, and J. A. Wellner. The support reduction algorithm for
computing non-parametric function estimates in mixture models. Scand. J. Statist.,
35(3):385–399, 2008.
P. Groeneboom, G. Jongbloed, and B. I. Witte. Maximum smoothed likelihood estimation and smoothed maximum likelihood estimation in the current status model.
Ann. Statist., 38(1):352–387, 2010.
P. Groeneboom and J. A. Wellner. Information bounds and nonparametric maximum
likelihood estimation, volume 19 of DMV Seminar. Birkhäuser Verlag, Basel, 1992.
P. Groeneboom and J. A. Wellner. Computing Chernoff’s distribution. J. Comput.
Graph. Statist., 10(2):388–400, 2001.
A. Guntuboyina, D. Lieu, S. Chatterjee, and B. Sen. Spatial adaptation in trend
filtering. arXiv preprint arXiv:1702.05113, 2017.
A. Guntuboyina and B. Sen. Global risk bounds and adaptation in univariate convex
regression. Probab. Theory Related Fields, 163(1-2):379–411, 2015.
Q. Han, T. Wang, S. Chatterjee, and R. J. Samworth. Isotonic regression in general
dimensions. arXiv preprint arXiv:1708.09468, 2017.
D. L. Hanson and G. Pledger. Consistency in concave regression. Ann. Statist.,
4(6):1038–1050, 1976.
D. L. Hanson, G. Pledger, and F. T. Wright. On consistency in monotonic regression.
Ann. Statist., 1:401–421, 1973.
T. J. Hastie and R. J. Tibshirani. Generalized additive models, volume 43 of Monographs on Statistics and Applied Probability. Chapman and Hall, Ltd., London, 1990.
C. Hildreth. Point estimates of ordinates of concave functions. J. Amer. Statist.
Assoc., 49:598–619, 1954.
J. Huang and J. A. Wellner. Estimation of a monotone density or monotone hazard
under random censoring. Scand. J. Statist., 22(1):3–33, 1995.
J. Huang and J. A. Wellner. Interval censored survival data: A review of recent
progress. Proceedings of the First Seattle Symposium in Biostatistics: Survival Analysis, Ed. D. Y. Lin and T. R. Fleming, pages 123–169, 1997.
G. Jongbloed and F. H. van der Meulen. Estimating a concave distribution function
from data corrupted with additive noise. Ann. Statist., 37(2):782–815, 2009.
S. M. Kakade, V. Kanade, O. Shamir, and A. Kalai. Efficient learning of generalized linear and single index models with isotonic regression. In Advances in Neural
Information Processing Systems, pages 927–935, 2011.
J. Kim and D. Pollard. Cube root asymptotics. Ann. Statist., 18(1):191–219, 1990.
J. Komlós, P. Major, and G. Tusnády. An approximation of partial sums of independent RV’s, and the sample DF. II. Z. Wahrscheinlichkeitstheorie und Verw. Gebiete,
34(1):33–58, 1976.
M. R. Kosorok. Bootstrapping in Grenander estimator. In Beyond parametrics in
interdisciplinary research: Festschrift in honor of Professor Pranab K. Sen, volume 1
of Inst. Math. Stat. (IMS) Collect., pages 282–292. Inst. Math. Statist., Beachwood,
imsart-generic ver. 2014/10/16 file: ShapeRegArxiv.tex date: September 19, 2017

A. Guntuboyina and B. Sen/Shape-restricted Regression

46

OH, 2008.
[77] A. K. Kuchibhotla, R. K. Patra, and B. Sen. Efficient estimation in convex single
index models. arXiv preprint arXiv:1708.00145, 2017.
[78] V. N. Kulikov and H. P. Lopuhaä. The behavior of the NPMLE of a decreasing density
near the boundaries of the support. Ann. Statist., 34(2):742–768, 2006.
[79] T. Kuosmanen. Representation theorem for convex nonparametric least squares.
Econometrics J., 11(2):308–325, 2008.
[80] R. Kyng, A. Rao, and S. Sachdeva. Fast, provable algorithms for isotonic regression
in all l p-norms. In Advances in Neural Information Processing Systems, pages 2719–
2727, 2015.
[81] S. Leurgans. Asymptotic distributions of slope-of-greatest-convex-minorant estimators. Ann. Statist., 10(1):287–296, 1982.
[82] K.-C. Li and N. Duan. Regression analysis under link violation. Ann. Statist.,
17(3):1009–1052, 1989.
[83] Q. Li and J. S. Racine. Nonparametric econometrics. Princeton University Press,
Princeton, NJ, 2007. Theory and practice.
[84] E. Mammen. Estimating a smooth monotone regression function. Ann. Statist.,
19(2):724–740, 1991.
[85] E. Mammen. Nonparametric regression under qualitative smoothness assumptions.
Ann. Statist., 19(2):741–759, 1991.
[86] E. Mammen, O. Linton, and J. Nielsen. The existence and asymptotic properties of
a backfitting projection algorithm under weak conditions. Ann. Statist., 27(5):1443–
1490, 1999.
[87] E. Mammen and C. Thomas-Agnan. Smoothing splines and shape restrictions. Scand.
J. Statist., 26(2):239–252, 1999.
[88] E. Mammen and K. Yu. Additive isotone regression. IMS Lecture Notes-Monograph
Series Asymptotics: Particles, Processes, and Inverse Problems, 55:179–195, 2007.
[89] E. Mammen and K. Yu. Additive isotone regression. In Asymptotics: particles, processes and inverse problems, pages 179–195. Institute of Mathematical Statistics, 2007.
[90] A. W. Marshall, I. Olkin, and B. C. Arnold. Inequalities: theory of majorization and
its applications. Springer Series in Statistics. Springer, New York, second edition,
2011.
[91] R. Mazumder, A. Choudhury, G. Iyengar, and B. Sen. A computational framework for
multivariate convex regression and its variants. J. Amer. Statist. Assoc. (accepted);
arXiv preprint arXiv:1509.08165, 2015.
[92] M. Meyer and M. Woodroofe. On the degrees of freedom in shape-restricted regression.
Ann. Statist., 28(4):1083–1104, 2000.
[93] M. C. Meyer. An extension of the mixed primal-dual bases algorithm to the case of
more constraints than dimensions. J. Statist. Plann. Inference, 81(1):13–31, 1999.
[94] M. C. Meyer. Inference using shape-restricted regression splines. Ann. Appl. Stat.,
2(3):1013–1033, 2008.
[95] M. C. Meyer. Semi-parametric additive constrained regression. J. Nonparametr. Stat.,
25(3):715–730, 2013.
[96] M. C. Meyer. A simple new algorithm for quadratic programming with applications
in statistics. Comm. Statist. Simulation Comput., 42(5):1126–1139, 2013.
imsart-generic ver. 2014/10/16 file: ShapeRegArxiv.tex date: September 19, 2017

A. Guntuboyina and B. Sen/Shape-restricted Regression

47

[97] H. Mukerjee. Monotone nonparameteric regression. Ann. Statist., 16(2):741–750,
1988.
[98] S. A. Murphy, A. W. van der Vaart, and J. A. Wellner. Current status regression.
Math. Methods Statist., 8(3):407–425, 1999.
[99] A. S. Nemirovski, B. T. Polyak, and A. B. Tsybakov. Convergence rate of nonparametric estimates of maximum-likelihood type. Problemy peredachi informatsii,
21(4):17–33, 1985.
[100] W. K. Newey. Semiparametric efficiency bounds. Journal of applied econometrics,
5(2):99–135, 1990.
[101] S. Oymak and B. Hassibi. Sharp MSE bounds for proximal denoising. Found. Comput.
Math., 16(4):965–1029, 2016.
[102] D. Pollard. Convergence of stochastic processes. Springer Series in Statistics. SpringerVerlag, New York, 1984.
[103] J. L. Powell, J. H. Stock, and T. M. Stoker. Semiparametric estimation of index
coefficients. Econometrica, 57(6):1403–1430, 1989.
[104] B. L. S. Prakasa Rao. Estkmation of a unimodal density. Sankhyā Ser. A, 31:23–36,
1969.
[105] T. Robertson and F. T. Wright. Consistency in generalized isotonic regression. Ann.
Statist., 3:350–362, 1975.
[106] T. Robertson, F. T. Wright, and R. L. Dykstra. Order restricted statistical inference.
John Wiley & Sons, 1988.
[107] R. T. Rockafellar. Convex analysis. Princeton Mathematical Series, No. 28. Princeton
University Press, Princeton, N.J., 1970.
[108] A. Schrijver. Theory of linear and integer programming. Wiley-Interscience Series
in Discrete Mathematics. John Wiley & Sons, Ltd., Chichester, 1986. A WileyInterscience Publication.
[109] E. Seijo and B. Sen. Nonparametric least squares estimation of a multivariate convex
regression function. Annals of Statistics, 39:1633–1657, 2011.
[110] B. Sen, M. Banerjee, and M. Woodroofe. Inconsistency of bootstrap: the Grenander
estimator. Ann. Statist., 38(4):1953–1977, 2010.
[111] B. Sen and M. Woodroofe. Bootstrap confidence intervals for isotonic estimators in a
stereological problem. Bernoulli, 18(4):1249–1266, 2012.
[112] B. Sen and G. Xu. Model based bootstrap methods for interval censored data. Comput.
Statist. Data Anal., 81:121–129, 2015.
[113] N. B. Shah, S. Balakrishnan, A. Guntuboyina, and M. J. Wainwright. Stochastically
transitive models for pairwise comparisons: statistical and computational issues. IEEE
Trans. Inform. Theory, 63(2):934–959, 2017.
[114] D. Simchi-Levi, X. Chen, and J. Bramel. Convexity and supermodularity. The Logic
of Logistics: Theory, Algorithms, and Applications for Logistics and Supply Chain
Management, pages 13–32, 2005.
[115] M. Slawski and M. Hein. Non-negative least squares for high-dimensional linear models: consistency and sparse recovery without regularization. Electron. J. Stat., 7:3004–
3056, 2013.
[116] Q. F. Stout. Unimodal regression via prefix isotonic regression. Comput. Statist. Data
Anal., 53(2):289–297, 2008.
imsart-generic ver. 2014/10/16 file: ShapeRegArxiv.tex date: September 19, 2017

A. Guntuboyina and B. Sen/Shape-restricted Regression

48

[117] Q. F. Stout. Fastest isotonic regression algorithms, 2014.
[118] Q. F. Stout. Isotonic regression for multiple independent variables. Algorithmica,
71(2):450–470, 2015.
[119] M. Talagrand. Upper and lower bounds for stochastic processes: modern methods and
classical problems, volume 60. Springer Science & Business Media, 2014.
[120] J. W. Tukey. Curves as parameters, and touch estimation. In Proc. 4th Berkeley Sympos. Math. Statist. and Prob., Vol. I, pages 681–694. Univ. California Press, Berkeley,
Calif., 1961.
[121] S. Van de Geer. Estimating a regression function. Annals of Statistics, 18:907–924,
1990.
[122] S. van de Geer and M. J. Wainwright. On concentration for (regularized) empirical
risk minimization. Sankhya A, pages 1–42, 2017.
[123] A. W. van der Vaart and J. A. Wellner. Weak convergence and empirical processes.
Springer Series in Statistics. Springer-Verlag, New York, 1996. With applications to
statistics.
[124] C. van Eeden. Maximum likelihood estimation of ordered probabilities. Nederl. Akad.
Wetensch. Proc. Ser. A. 59, Indag. Math., 18:444–455, 1956.
[125] M. Woodroofe and J. Sun. A penalized maximum likelihood estimate of f (0+) when
f is nonincreasing. Statist. Sinica, 3(2):501–515, 1993.
[126] F. T. Wright. The asymptotic behavior of monotone regression estimates. Ann.
Statist., 9(2):443–448, 1981.
[127] J. Wu, M. C. Meyer, and J. D. Opsomer. Penalized isotonic regression. J. Statist.
Plann. Inference, 161:12–24, 2015.
[128] A. Yatchew. Semiparametric regression for the applied econometrician. Cambridge
University Press, 2003.
[129] C.-H. Zhang. Risk bounds in isotonic regression. Ann. Statist., 30(2):528–555, 2002.

imsart-generic ver. 2014/10/16 file: ShapeRegArxiv.tex date: September 19, 2017

