Convex Regression in Multidimensions: Suboptimality of Least Squares Estimators

The main goal of this paper is to show that nonparametric Least Squares Estimators (LSEs) associated with the constraint of convexity can be minimax suboptimal when the underlying dimension is 5 or larger. Specifically, we consider regression over a variety of different classes of convex functions (and covariate designs) and find lower bounds on the rates of convergence of the corresponding least squares estimator over each class, and show that, when $d\geq 5$, these rates of convergence do not match the minimax rate for the given class.

We work in the standard nonparametric regression setting for estimating an unknown convex function $f_{0}:\Omega\rightarrow{\mathbb{R}}$ defined on a known full-dimensional compact convex domain $\Omega\subseteq{\mathbb{R}}^{d}$ ($d\geq 1$) from observations $(X_{1},Y_{1}),\dots,(X_{n},Y_{n})$ generated via the model:

where $\xi_{1},\ldots,\xi_{n}$ are i.i.d. errors having the $N(0,\sigma^{2})$ distribution, and design points $X_{1},\dots,X_{n}\in\Omega$ that may be fixed or random. This problem is known as convex regression and has a long history in statistics and related fields. Standard references are [29, 28, 22, 21, 44, 31, 35, 3], and applications of convex regression can be found in [49, 50, 2, 38, 1, 30, 46].

The Least Squares Estimator (LSE) over a class ${\cal F}$ of functions on $\Omega$ is defined as any minimizer of the least squares criterion over ${\cal F}$:

In order to introduce the specific LSEs studied in this paper, consider the following four function classes:

1. 1.

   ${\mathcal{C}}_{L}^{B}(\Omega)$: class of all convex functions on $\Omega$ which are uniformly Lipschitz with Lipschitz constant $L$ and uniformly bounded by $B$.

2. 2.

   ${\mathcal{C}}_{L}(\Omega)$: class of all convex functions on $\Omega$ that are uniformly Lipschitz with Lipschitz constant $L$. There is no uniform boundedness assumption on functions in this class.

3. 3.

   ${\mathcal{C}}^{B}(\Omega)$: class of all convex functions on $\Omega$ that are uniformly bounded by $B$. There is no Lipschitz assumption on functions in this class.

4. 4.

   ${\mathcal{C}}(\Omega)$: class of all convex functions on $\Omega$.

Throughout the paper, we assume that $L$ and $B$ are positive constants not depending on the sample size $n$. We study LSEs over each of the above function classes and we use the following terminology for these estimators:

1. 1.

   $\hat{f}_{n}({\mathcal{C}}_{L}^{B}(\Omega))$: Bounded Lipschitz Convex LSE,

2. 2.

   $\hat{f}_{n}({\mathcal{C}}_{L}(\Omega))$: Lipschitz Convex LSE,

3. 3.

   $\hat{f}_{n}({\mathcal{C}}^{B}(\Omega))$: Bounded Convex LSE,

4. 4.

   $\hat{f}_{n}({\mathcal{C}}(\Omega))$: Unrestricted Convex LSE.

Strictly speaking, $\hat{f}_{n}({\mathcal{C}}_{L}^{B}(\Omega))$ should be called the “Uniformly Bounded Uniformly Lipschitz Convex LSE” but we omit the word “uniformly” throughout for brevity. Table 1 gives references to some previous papers that studied each of these estimators. For the unrestricted convex LSE, Table 1 lists only the papers focusing on the multivariate case (univariate investigations are in [28, 17, 22, 14, 20, 23, 12, 6, 11]).

This paper studies the performance of these LSEs and compares their rates of convergence with the corresponding minimax rates. For our results, we assume throughout the paper that the convex body $\Omega$ (which is the domain of the unknown function $f_{0}$) is translated and scaled so that

where ${\mathfrak{B}_{d}}$ is the unit ball in ${\mathbb{R}}^{d}$ and $r_{d}$ is a positive constant depending on $d$ alone. In particular, this assumption implies that $\Omega$ has diameter at most 2 and also volume at most 1. Also the diameter and volume of $\Omega$ are bounded from below by a constant depending on $d$ alone. The classical John’s theorem (see e.g., [5, Lecture 3]) shows that for every convex body $\Omega$, there exists an affine transformation such that (3) holds for the transformed body with $r_{d}=1/d$. In the rest of the paper, whenever we say that a result holds for “any convex body $\Omega$” (such as in Tables 2, 3, 4, 5 and 6), we mean “any convex body $\Omega$ satisfying (3)”.

We focus first on the random design setting where the design points $X_{1},\dots,X_{n}$ are assumed to be independent having the uniform distribution ${\mathbb{P}}$ on $\Omega$ (see Subsection 5.3 for a discussion of more general design assumptions in random design), and work with the loss function

The minimax rate for a function class ${\cal F}$ in the above setting is:

where ${\mathbb{E}}_{f_{0}}$ denotes expectation with respect to the joint distribution of all the observations $(X_{1},Y_{1}),\dots,(X_{n},Y_{n})$ when $f_{0}$ is the true regression function (see (1)), and the infimum is over all estimators $\breve{f}_{n}$.

Minimax rates for the above function classes 1–4 are presented in Table 2. These results are known except perhaps for ${\mathcal{C}}_{L}(\Omega)$ which we state and prove in this paper as Proposition 2.1. These minimax rates can also be derived from the corresponding metric entropy rates which are given in Table 3. Recall that the $\epsilon$-metric entropy of a function class ${\cal F}$ with respect to a metric $\ell$ is defined as the logarithm of the smallest number $N(\epsilon,{\cal F},\ell)$ of closed balls of radius $\epsilon$ whose union contains ${\cal F}$. Classical results from Yang and Barron [51] imply that the minimax rate (for the loss $\ell_{{\mathbb{P}}}^{2}$) for a nonparametric function class ${\cal F}$ satisfying some natural assumptions (see [51, Subsection 3.2] for a listing of these assumptions which are satisfied for ${\mathcal{C}}_{L}^{B}(\Omega)$ and ${\mathcal{C}}^{B}(\Omega)$) equals $\epsilon^{2}_{n}$ where $\epsilon_{n}$ solves the equation:

For ${\mathcal{C}}_{L}^{B}(\Omega)$, the minimax rate is $n^{-4/(d+4)}$ which is the square of the solution for $\epsilon$ in (5): $n\epsilon^{2}=\epsilon^{-d/2}$. For ${\mathcal{C}}_{L}(\Omega)$, the metric entropy is infinite as this class contains every constant function and hence is not totally bounded. But the minimax rate is still $n^{-4/(d+4)}$ and can be easily derived as a consequence of the minimax rate result for ${\mathcal{C}}_{L}^{L}(\Omega)$ (the argument is given in the proof of Proposition 2.1).

For ${\mathcal{C}}^{B}(\Omega)$ with $d\geq 2$, as shown by Han and Wellner [26] and Gao and Wellner [19] respectively, the minimax rate and the metric entropy rate change depending on whether $\Omega$ is polytopal or has smooth boundary. Here (and in the rest of the paper), when we say that a specific set $S$ is a “polytope”, we assume that the number of facets or extreme points of $S$ is bounded by a constant not depending on the sample size $n$. If the number of facets is allowed to grow with $n$, then polytopes can approximate general convex bodies arbitrarily well in which case it is not meaningful to give separate rates for polytopes and other convex bodies.

When $\Omega$ is a polytope, the minimax rate for ${\mathcal{C}}^{B}(\Omega)$ equals $n^{-4/(d+4)}$ corresponding to the $\epsilon$-entropy of $\epsilon^{-d/2}$. In contrast, when $\Omega$ is the unit ball, the minimax rate is $n^{-2/(d+1)}$ corresponding to the $\epsilon$-entropy of $\epsilon^{-(d-1)}$. The larger metric entropy for the ball (compared to polytopes) is driven by the curvature of the boundary and the lower bound is proved by considering indicator-like functions of spherical caps (see [19, Subsection 2.10] for the metric entropy proofs). More high level intuition on the differences between the polytopal case and the ball case is provided in Subsection 5.1.

We now ask whether the LSEs (as defined via (2)) achieve the corresponding minimax rates. The supremum risk of an estimator $\breve{f}$ on a function class ${\cal F}$ is defined as

Existing results for the convex regression LSEs are described in Table 4 where the rate $r_{n,d}$ is given by

The upper bound of $r_{n,d}$ in the first three rows of Table 4 is derived from standard upper bounds on the performance of LSEs [7, 47] which state that the LSE risk is bounded from above by $\delta_{n}^{2}$ where $\delta_{n}$ solves the equation:

Here $N_{[\,]}(\epsilon,{\cal F},\ell_{{\mathbb{P}}})$ denotes the $\epsilon$-bracketing number of ${\cal F}$ under the $\ell_{{\mathbb{P}}}$ metric which is defined as the smallest number of $\epsilon$-brackets:

needed to cover ${\cal F}$. The logarithm of $N_{[\,]}(\epsilon,{\cal F},\ell_{{\mathbb{P}}})$ denotes the $\epsilon$-bracketing entropy of ${\cal F}$. Gao and Wellner [19] (see also Doss [15]) proved bracketing entropy numbers for the convex function classes; these are given in Table 3 and they coincide with the metric entropy rates. Plugging in $\epsilon^{-d/2}$ for the bracketing entropy in (7) and solving the resulting equation in $\delta$ leads to the rate $r_{n,d}$ (although ${\mathcal{C}}_{L}(\Omega)$ has infinite entropy, results for this class can be derived by working instead with ${\mathcal{C}}_{L}^{L}(\Omega)$ which has entropy $\epsilon^{-d/2}$). The split in $r_{n,d}$ for $d\leq 4$ and $d\geq 5$ occurs because $\int_{\delta^{2}}^{\delta}\epsilon^{-d/4}d\epsilon$ depends differently on $\delta$ for $d\leq 4$ compared to $d\geq 5$. The same dimension-dependent split in the upper bounds for the rate can be seen in [7] and [47, Chapter 9] for certain other LSEs on function classes with smoothness restrictions.

It is interesting to note that the last row in Table 4 (corresponding to the case when $\Omega$ is the unit ball) does not have any dimension-dependent split in the rate (it equals $n^{-2/(d+1)}$ for every $d\geq 2$). This result is due to Kur et al. [32]. More details on this case are given in Subsection 5.1.

We are now ready to describe the main results of this paper. Tables 2 and 4 together indicate that the minimax rate for the first three rows equals $n^{-4/(d+4)}$ which also coincides with $r_{n,d}$ for $d\leq 4$. This immediately implies minimax optimality (up to log factors) of the LSE over each function class in the first three rows of Table 2 when $d\leq 4$. However for $d\geq 5$, there is a significant gap between the minimax rate $n^{-4/(d+4)}$ and $r_{n,d}=n^{-2/d}$. The question of whether this gap is real (implying suboptimality of the LSEs) or merely an artifact of some loose argument in the derived upper bounds was previously open. We settle this question in this paper by proving that the LSE for each function class in the first three rows of Table 2 is minimax suboptimal for $d\geq 5$. More precisely, in Theorem 3.1, we prove that $n^{-2/d}(\log n)^{-4(d+1)/d}$ is a lower bound (up to a multiplicative factor independent of $n$) on each supremum risk in the first three rows of Table 4 for $d\geq 5$.

Table 5 has a summary of the minimax rates and rates exhibited by the LSE (along with minimax optimality and suboptimality of the LSE) in our random design setting.

Our proof of minimax suboptimality for the LSE is based on the following idea. Suppose $f_{0}(x):=\|x\|^{2}$ and let $\tilde{f}_{k}$ be a piecewise affine convex approximation to $f_{0}$ with $k$ affine pieces as described in the statment of Lemma 3.2. It is then well-known (see e.g., [4, Lemma 4.1] or [8]) that the $\ell_{{\mathbb{P}}}$ distance between $f_{0}$ and $\tilde{f}_{k}$ is of order $k^{-2/d}$. If we now set $k$ to be of order $\sqrt{n}$, then $\ell_{{\mathbb{P}}}(\tilde{f}_{k},f_{0})\sim n^{-1/d}$, or equivalently $\ell_{{\mathbb{P}}}^{2}(\tilde{f}_{k},f_{0})\sim n^{-2/d}$. Note that $n^{-2/d}$ is much larger than the minimax rate $n^{-4/(d+4)}$ for $d\geq 5$. We study the behavior of each convex LSE $\hat{f}_{n}$ (in the first three rows of Table 2) when the true function is $\tilde{f}_{k}$ for $k\sim\sqrt{n}$, and show (in Theorem 3.1) that ${\mathbb{E}}\ell_{{\mathbb{P}}}^{2}(\hat{f}_{n},\tilde{f}_{k})$ is bounded below by a constant multiple of $n^{-2/d}(\log n)^{-4(d+1)/d}$. In other words, we prove that the LSEs are at a squared distance from the true function $\tilde{f}_{k}$ that is much larger than the minimax rate when $d\geq 5$. Our proof techniques are outlined in Section 4. We are unable to definitively say why the LSEs are behaving this way when the true function is piecewise affine with $\sqrt{n}$ pieces. We believe this is due to overfitting, probably the LSEs have many more affine pieces than $k\sim\sqrt{n}$ in this case and are actually closer to the quadratic function $f_{0}$.

It has previously been observed that LSEs and related Empirical Risk Minimization procedures can be suboptimal for function classes with large entropy integrals. For example, [7, Section 4] designed certain pathological large function classes where the LSE is provably minimax suboptimal. The fact that minimax suboptimality also holds for natural function classes such as ${\mathcal{C}}_{L}^{B}(\Omega)$ and ${\mathcal{C}}_{L}(\Omega)$ (for $d\geq 5$) is noteworthy.

In contrast to our suboptimality results, the last row of Table 4 and Table 5 correspond to the class ${\mathcal{C}}^{B}({\mathfrak{B}_{d}})$ (here ${\mathfrak{B}_{d}}$ is the unit ball) for which the LSE is minimax optimal for all $d\geq 2$ (this result is due to Kur et al. [32]). As already mentioned, ${\mathcal{C}}^{B}({\mathfrak{B}_{d}})$ has much larger metric entropy (compared to ${\mathcal{C}}^{B}(\Omega)$ for polytopal $\Omega$) that is driven by the curvature of ${\mathfrak{B}_{d}}$; one can understand this from the proofs of the metric entropy lower bounds on ${\mathcal{C}}^{B}({\mathfrak{B}_{d}})$ and ${\mathcal{C}}^{B}(\Omega)$ for polytopal $\Omega$ in [19] (specifically see [19, proof of Theorem 1(i) and proof of Theorem 4]). High-level intuition on the differences between the ball and polytopal cases is given in Subsection 5.1. Isotonic regression is another shape constrained regression problem where the LSE is minimax optimal for all $d$ (see Han et al. [25]). The class of coordinatewise monotone functions on $[0,1]^{d}$ is similar to ${\mathcal{C}}^{B}({\mathfrak{B}_{d}})$ in that its metric entropy is driven by well-separated subsets of the domain $[0,1]^{d}$ (see [18, Proof of Proposition 2.1]). Other examples of such classes where the LSE is optimal for all dimensions can be found in Han [24].

In this section, we provide more details for the minimax rate results mentioned in the Introduction.