Review on Isotonic and Convex Regression

1. Introduction

One of the central problems in statistics and simulation is to estimate an unknown function $f_0:\mathrm{\Omega }\subset {\mathbb{R}}^d\to \mathbb{R}$ from noisy observations $({\mathit{X}}_1,Y_1),\dots ,({\mathit{X}}_n,Y_n)$, where we assume

$$Y_i=f_0\left({\mathit{X}}_i\right)+{\epsilon }_i,$$

for $i=1,\dots ,n$, $\mathrm{\Omega }$ is a closed convex subset of ${\mathbb{R}}^d$, and the ${\epsilon }_i$’s are independent and identically distributed (iid) mean zero random variables with finite variance ${\sigma }^2$. This type of problem arises in many settings. For instance, $\mathit{x}$ may represent a parameter in a stochastic system, such as the service rate of a single-server queue, and $f_0\left(\mathit{x}\right)$ a performance metric, such as the steady-state mean waiting time in that single server queue when the service rate equals $\mathit{x}$.

In many applications, there is no-closed form formula for $f_0\left(\mathit{x}\right)$. In such a case, a practical approach is to simulate the system at selected design points $\mathit{x}={\mathit{X}}_1,\dots ,{\mathit{X}}_n$, obtain simulated outputs $Y_1,\dots ,Y_n$ as estimates of $f_0\left({\mathit{X}}_1\right),\dots ,f_0\left({\mathit{X}}_n\right)$, respectively, and use the data $({\mathit{X}}_1,Y_1),\dots ,$$({\mathit{X}}_n,Y_n)$ to approximate $f_0$: This is a typical regression problem. The problem becomes challenging when it is impossible to place the relationship into a simple parametric form. However, when $f_0$ is known to satisfy certain structural properties such as monotonicity or convexity, shape-restricted regression offers powerful estimators.

Given such prior shape information, a natural estimator of $f_0$ fits a function f with such a shape characterization to the data set and solves

$$\begin{array}{cccc}\hfill & \underset{f\in \mathcal{F}}{\mathrm{minimize}}\hfill & \hfill & {\displaystyle \frac{1}{n}}\sum _{i=1}^n{\left(Y_i-f\left({\mathit{X}}_i\right)\right)}^2,\hfill \end{array}$$

(1)

where $\mathcal{F}$ is the class of functions satisfying the assumed shape constraint. For example, if $f_0$ is known to be nondecreasing, we can set $\mathcal{F}={\mathcal{F}}_m$, where

$$\begin{array}{c}{\mathcal{F}}_m=\{f:\mathrm{\Omega }\to \mathbb{R}\mid f\left(\mathbf{v}\right)\le f\left(\mathbf{w}\right)if{v}_i\le w_{i}fori=1,\dots ,d,\hfill \\ \hfill \mathbf{v}=(v_1,\dots ,v_d),\mathbf{w}=(w_1,\dots ,w_d)\}.\end{array}$$

On the other hand, if $f_0$ is known to be convex, then we can set $\mathcal{F}={\mathcal{F}}_c$, where

$$\begin{array}{c}{\mathcal{F}}_c=\{f:\mathrm{\Omega }\to \mathbb{R}\mid f\mathrm{is} \mathrm{convex},\mathrm{i}.\mathrm{e}.,\hfill \\ \hfill f(\alpha \mathbf{v}+(1-\alpha \left)\mathbf{w}\right)\le \alpha f\left(\mathbf{v}\right)+(1-\alpha )f\left(\mathbf{w}\right)\mathrm{for}\alpha \in [0,1]\mathrm{and}\mathbf{v},\mathbf{w}\in \mathrm{\Omega }\}.\end{array}$$

Throughout this paper, we will focus on these two types of problems, the one with $\mathcal{F}={\mathcal{F}}_m$ which leads to isotonic regression and the other with $\mathcal{F}={\mathcal{F}}_c$ which leads to convex regression.

One of the first questions that arises is how to compute the solution to (1) numerically when $\mathcal{F}={\mathcal{F}}_m$ or $\mathcal{F}={\mathcal{F}}_c$. Below, we describe how the solution to (1) can be obtained by solving quadratic programming (QP) problems.

1.1. QP Formulation for Isotonic Regression When $d=1$

When $f_0$ is nondecreasing in $\mathit{x}$, as when price increases with demand, we estimate $f_0$ by fitting a nondecreasing function f that minimizes squared errors ${\sum }_{i=1}^n{(Y_i-f\left({\mathit{X}}_i\right))}^2$. Writing $f_i\triangleq f\left({\mathit{X}}_i\right)$ and assuming the design is ordered so ${\mathit{X}}_1<\dots <{\mathit{X}}_n$, the existence of a nondecreasing fit through $({\mathit{X}}_i,f_i)$ is equivalent to $f_1\le \dots \le f_n$. Thus, the estimator solves the following quadratic programming (QP) problem:

$$\begin{array}{cccc}\hfill & \underset{f_1,\dots ,f_n\in \mathbb{R}}{\mathrm{minimize}}\hfill & \hfill & {\displaystyle \frac{1}{n}}\sum _{i=1}^n{(Y_i-f_i)}^2\hfill \\ \hfill & \mathrm{subject}\mathrm{to}\hfill & \hfill & f_1\le f_2\le \dots \le f_n.\hfill \end{array}$$

(2)

Let ${\widehat{f}}_1^m,\dots ,{\widehat{f}}_n^m$ be a solution to (2). Any nondecreasing function interpolating $(({\mathit{X}}_i,{\widehat{f}}_i^m):1\le i\le n)$ defines the isotonic regression estimator and we will denote it by ${\widehat{f}}_n^m(·)$; see, for example, Brunk [1].

1.2. QP Formulation for Convex Regression When $d=1$

Suppose instead $f_0$ is convex. Writing $f_i\triangleq f\left({\mathit{X}}_i\right)$ and assuming the design is ordered so ${\mathit{X}}_1<\dots <{\mathit{X}}_n$, a convex fit through $({\mathit{X}}_i,f_i)$ exists if and only if

$${\displaystyle \frac{f_2-f_1}{{\mathit{X}}_2-{\mathit{X}}_1}}\le \dots \le {\displaystyle \frac{f_n-f_{n-1}}{{\mathit{X}}_n-{\mathit{X}}_{n-1}}}.$$

Hence, the least-squares convex fit solves

$$\begin{array}{cccc}\hfill & \underset{f_1,\dots ,f_n\in \mathbb{R}}{\mathrm{minimize}}\hfill & \hfill & {\displaystyle \frac{1}{n}}\sum _{i=1}^n{(Y_i-f_i)}^2\hfill \\ \hfill & \mathrm{subject}\mathrm{to}\hfill & \hfill & {\displaystyle \frac{f_2-f_1}{{\mathit{X}}_2-{\mathit{X}}_1}}\le \dots \le {\displaystyle \frac{f_n-f_{n-1}}{{\mathit{X}}_n-{\mathit{X}}_{n-1}}}.\hfill \end{array}$$

(3)

Let ${\widehat{f}}_1^c,\dots ,{\widehat{f}}_n^c$ be a solution to (3). Any convex function interpolating $(({\mathit{X}}_i,{\widehat{f}}_i^c):1\le i\le n)$ is a convex regression estimator and we will denote it by ${\widehat{f}}_n^c(·)$; see, for example, Hildreth [2].

1.3. Key Questions

Isotonic and convex regression are nonparametric methods that avoid any bandwidth or tuning parameters in their basic form yet enjoy strong theoretical guarantees. In particular, as $n\to \infty$, these estimators converge to $f_0$ under standard conditions even when the noise variance $Var\left(Y_i\right)$ does not decrease. This is due to the fact that the estimators are required to satisfy the shape constraints, and hence, just adding more and more data points enforces the estimators to converge to the true underlying function.

Over the past seven decades, extensive work was done to study their statistical properties, computational algorithms, and practical applications. This paper reviews the literature on these topics and highlights unresolved questions. We focus on answering the following questions:

• How can the estimator be computed numerically from data $({\mathit{X}}_1,Y_1),\dots ,({\mathit{X}}_n,Y_n)$?
• What are its statistical properties?
• Are there efficient algorithms that are tailored to the computation of the isotonic or convex regression estimators?
• What applications benefit most from these approaches?
• What are its limitations, and how can we overcome them?
• What are the challenges that still lie in this field?

1.4. Organization of This Paper

The remainder of the paper is organized as follows. Section 2 reviews how the isotonic and convex regression estimators can be computed via QP problems. Section 3 focuses on isotonic regression while Section 4 examines convex regression in greater detail. Section 5 includes some conclusions while Section 6 discusses open challenges and future research directions.

1.5. Notation and Definitions

Throughout this paper, we view vectors as columns and write ${\mathit{x}}^{\mathsf{T}}$ to denote the transpose of a vector $\mathit{x}\in {\mathbb{R}}^d$. For $\mathit{x}\in {\mathbb{R}}^d$, we write its ith component as $x_i$, so $\mathit{x}={(x_1,\dots ,x_d)}^{\mathsf{T}}$. We define $\parallel \mathit{x}\parallel ={(x_1^2+\dots +x_d^2)}^{1/2}$ and ${\parallel \mathit{x}\parallel }_{\infty }={max}_{1\le i\le d}\left|x_i\right|$.

For any convex set $S\subset {\mathbb{R}}^d$ and a convex function $g:S\to \mathbb{R}$, a vector $\beta \in {\mathbb{R}}^d$ is said to be a subgradient of g at $\mathit{x}\in S$ if $g\left(\mathit{y}\right)\ge g\left(\mathit{x}\right)+{\beta }^{\mathsf{T}}(\mathit{y}-\mathit{x})$ for all $\mathit{y}\in S$. We denote a subgradient of f at $\mathit{x}\in S$ by $subgradg\left(\mathit{x}\right)$. The set of all subgradients of g at $\mathit{x}$ is called the subdifferential of g at $\mathit{x}$, denoted by $\partial g\left(\mathit{x}\right)$. When a function $g:S\to \mathbb{R}$ is differentiable at $\mathit{x}\in {[a,b]}^d$, we denote its gradient at $\mathit{x}$ by $\nabla g\left(\mathit{x}\right)$. When $g:\mathbb{R}\to \mathbb{R}$ is differentiable, we denote its gradient at $\mathit{x}$ by $g^{\prime }\left(\mathit{x}\right)$.

For a sequence of random variables $(Z_n:n\ge 1)$ and a sequence of positive real numbers $({\alpha }_n:n\ge 1)$, we say $Z_n={\mathcal{O}}_p\left({\alpha }_n\right)$ as $n\to \infty$ if, for any $ϵ>0$, there exist constants c and $n_0$ such that $\mathbb{P}\left(\right|Z_n/{\alpha }_n|>c)<ϵ$ for all $n\ge n_0$.

For any sequences of real numbers $(a_n:n\ge 1)$ and $(b_n:n\ge 1)$, $a_n=O\left(b\left(n\right)\right)$ if there exists a positive real number M and a real number $n_0$ such that $|a_n|\le Mb\left(n\right)$ for all $n\ge n_0$.

For any sequence of random variables $(A_n:n\ge 1)$ and any random variable B, $A_n\Rightarrow B$ means $A_n$ converges in distribution to B, or $A_n$ converges weakly to B as $n\to \infty$.

2. The Basic QP Formulation

In this section, we describe how the isotonic and convex regression estimators can be found by solving QP problems. Since there are many efficient algorithms for solving QP problems, expressing the isotonic and convex regression estimators as solutions to QP problems offers a basic way to compute them numerically. In $d=1$, (2) and (3) already provide QP formulations for the isotonic and convex regression estimators, respectively. We now present multivariate analogues.

2.1. Isotonic Regression When $d>1$

The main challenge in finding the QP formulation for the multivariate isotonic regression estimator is how to express the monotonicity of a multivariate function in linear inequalities as in the constraints of (2). We will define a non-decreasing function g as follows:

$$g\left(\mathbf{v}\right)\le g\left(\mathbf{w}\right)\mathrm{whenever}\mathbf{v}⪯\mathbf{w}$$

for $\mathbf{v}=(v_1,\dots ,v_d)$ and $\mathbf{w}=(w_1,\dots ,w_d)$ in ${\mathbb{R}}^d$, where $\mathbf{v}⪯\mathbf{w}$ indicates $v_i\le w_i$ for $1\le i\le d$.

Using this definition, the nondecreasing function with the minimum sum of squared errors can be found by solving the following QP problem in decision variables $f_1,\dots ,f_n\in \mathbb{R}$:

$$\begin{array}{cccc}\hfill & \mathrm{minimize}\hfill & \hfill & {\displaystyle \frac{1}{n}}\sum _{i=1}^n{(Y_i-f_i)}^2\hfill \\ \hfill & \mathrm{subject}\mathrm{to}\hfill & \hfill & f_i\le f_{j}if{\mathit{X}}_i⪯{\mathit{X}}_{j}for1\le i,j\le n.\hfill \end{array}$$

(4)

Let ${\widehat{f}}_1^m,\dots ,{\widehat{f}}_n^m$ be a solution to (4). Any nondecreasing function interpolating $(({\mathit{X}}_i,{\widehat{f}}_i^m):1\le i\le n)$ can serve as the isotonic regression estimator and we will denote it by ${\widehat{f}}_n^m(·)$

2.2. Convex Regression When $d>1$

Again, the key question is how to express the existence of a convex function going through $({\mathit{X}}_1,f_1),\dots ,({\mathit{X}}_n,f_n)$ as a set of linear constraints. One should notice that there exists a convex function passing through $({\mathit{X}}_1,f_1),\dots ,({\mathit{X}}_n,f_n)$ if and only if there exist ${\beta }_1,\dots ,{\beta }_n\in {\mathbb{R}}^d$ satisfying

$$f_i\ge f_j+{\beta }_j^{\mathsf{T}}({\mathit{X}}_i-{\mathit{X}}_j)for1\le i,j\le n;$$

see page 338 of Boyd and Vandenberghe [3]. Here, ${\beta }_i$ serves as a subgradient of the convex function f at ${\mathit{X}}_i$. Hence, a convex function with the least squares can be found by solving the following QP problem in decision variables $f_1,\dots ,f_n\in \mathbb{R}$ and ${\beta }_1,\dots ,{\beta }_n\in {\mathbb{R}}^d$:

$$\begin{array}{cccc}\hfill & \mathrm{minimize}\hfill & \hfill & {\displaystyle \frac{1}{n}}\sum _{i=1}^n{(Y_i-f_i)}^2\hfill \\ \hfill & \mathrm{subject}\mathrm{to}\hfill & \hfill & f_i\ge f_j+{\beta }_j^{\mathsf{T}}({\mathit{X}}_i-{\mathit{X}}_j)for1\le i,j\le n.\hfill \end{array}$$

(5)

Let ${\widehat{f}}_1^c,\dots ,{\widehat{f}}_n^c,{\widehat{\beta }}_1,\dots ,{\widehat{\beta }}_n$ be a solution to (5). Any convex function going through $({\mathit{X}}_1,{\widehat{f}}_1^c),\dots ,$$({\mathit{X}}_n,{\widehat{f}}_n^c)$ can serve as the convex regression estimator in the multidimensional case. However, to remove any ambiguity, we will define the convex regression estimator ${\widehat{f}}_n^c(·)$ as follows:

$${\widehat{f}}_n^c\left(\mathit{x}\right)=\underset{1\le i\le n}{max}\{{\widehat{f}}_i^c+{\beta }_i^{\mathsf{T}}(\mathit{x}-{\mathit{X}}_i)\}$$

for $\mathit{x}\in \mathrm{\Omega }$.

In (2), (3), (4), and (5), we established the isotonic and convex regression estimators as solutions to QP problems. Since there are many efficient numerical algorithms to solve QP problems, these QP formulations enable us to compute the isotonic and convex regression estimators numerically.

3. Isotonic Regression

In this section, we review statistical properties, interesting applications, and some challenges of isotonic regression estimator. We also review some algorithms for computing the isotonic regression estimator numerically.

3.1. Statistical Properties

3.1.1. Univariate Case

The isotonic regression estimator in one and multiple dimensions was first introduced by Brunk [1]. In one dimension, Brunk [4] established its strong consistency. More specifically, Brunk [4] proved that if $\mathrm{\Omega }=(a,b)\subset \mathbb{R}$, for any $c>a$ and $d<b$,

$$\underset{c\le \mathit{x}\le d}{max}|{\widehat{f}}_n^m\left(\mathit{x}\right)-f_0\left(\mathit{x}\right)|\to 0$$

as $n\to \infty$ with probability one.

Brunk [5] established the asymptotic distribution theory for the isotonic regression estimator when $d=1$. Brunk [5] proved that when $f_0^{\prime }\left(\mathit{x}\right)\ne 0$,

$$n^{1/3}\left({\widehat{f}}_n^m\left(\mathit{x}\right)-f_0\left(\mathit{x}\right)\right)\Rightarrow {\left[4{\sigma }^2{f}_0^{\prime }\left(\mathit{x}\right)\right]}^{1/3}{argmin}_{t\in \mathbb{R}}\{W\left(t\right)+t^2\}$$

as $n\to \infty$, where $W(·)$ is the two-sided Brownian motion starting from 0.

The $l_2$-risk of the isotonic regression estimator, defined by

$$R(f_0,{\widehat{f}}_n^m)\triangleq \mathbb{E}\left[{\displaystyle \frac{1}{n}}\sum _{i=1}^n{\left({\widehat{f}}_n^m\left({\mathit{X}}_i\right)-f_0\left({\mathit{X}}_i\right)\right)}^2\right],$$

(6)

is proven to have the following bound

$$R(f_0,{\widehat{f}}_n^m)\le C{\sigma }^{4/3}V{\left(f_0\right)}^{2/3}{n}^{-2/3},$$

where $V\left(f_0\right)$ is the total variation of $f_0$ and $C>0$ is a universal constant. Meyer and Woodroofe [6] proved (6) under the assumption that the ${\epsilon }_i$’s are normally distributed, and Zhang [7] proved (6) under a more general assumption that the ${\epsilon }_i$’s are iid with mean zero and variance ${\sigma }^2$.

3.1.2. Multivariate Case

In multiple dimensions, Lim [8] established the strong consistency. Corollary 3.1 of Lim [8] proved that when $\mathrm{\Omega }={(0,1)}^d$, for any $ϵ>0$

$$\underset{\mathit{x}\in {[ϵ,1-ϵ]}^d}{sup}|{\widehat{f}}_n^m\left(\mathit{x}\right)-f_0\left(\mathit{x}\right)|\to 0$$

as $n\to \infty$ with probability one.

Lim [8] also established the rate of convergence and proved that

$$\mathbb{E}\left[{\left({\tilde{f}}_n^m\left(\mathit{X}\right)-f_0\left(\mathit{X}\right)\right)}^2\right]=\left\{\begin{array}{cc}{\mathcal{O}}_p\left(n^{-2/3}\right),\hfill & \mathrm{if}d=1\hfill \\ {\mathcal{O}}_p\left(n^{-1/2}{(logn)}^2\right),\hfill & \mathrm{if}d=2\hfill \\ {\mathcal{O}}_p\left(n^{-{\displaystyle \frac{1}{2(d-1)}}}\right),\hfill & \mathrm{if}d>2,\hfill \end{array}\right.$$

where the expectation is taken with respect to $\mathit{X}$, $\mathit{X}$ has the same distribution as the ${\mathit{X}}_i$’s, and ${\tilde{f}}_n^m(·)$, a variant of the isotonic regression estimator, minimizes the sum of squared errors ${\sum }_{i=1}^n{(Y_i-f\left({\mathit{X}}_i\right))}^2$ over all nondecreasing functions f satisfying $\left|f\right(\mathit{x}\left)\right|\le C$ for some pre-set constant C for all $\mathit{x}\in \mathrm{\Omega }$.

Bagchi and Dhar [9] derived the pointwise limit distribution of the isotonic regression estimator in multiple dimensions. In particular, Bagchi and Dhar [9] proved

$$n^{-1/(d+2)}({\widehat{f}}_n^m\left(\mathit{x}\right)-f_0\left(\mathit{x}\right))\Rightarrow H$$

as $n\to \infty$ for some non-degenerate limit distribution H.

Han et al. [10] showed that the $l_2$-risk of a multivariate isotonic regression estimator is bounded by a constant multiple of $n^{-1/d}{(logn)}^4$ for all $d\ge 1$.

3.2. Computational Algorithms

Formulations (2) and (3) are QPs with linear constraints, and hence, can be solved by off-the-shelf solvers such as CVX, CPLEX, or MOSEK. However, there exist algorithms tailored to the specific setting of the isotonic regression estimator. In $d=1$, Barlow et al. [11] discussed a graphical procedure to compute the isotonic regression estimator by identifying it as the greatest common minorant for a cumulative sum diagram, and proposed a tabular procedure called the pool-adjacent-violators (PAV) algorithm. PAV algorithm computes the isotonic regression estimator in $O\left(n\right)$ steps. Various modifications of the PAV algorithm have been proposed in the literature; see Brunk [1], Mukerjee [12], Mammen [13], Friedman and Tibshirani [14], Ramsay [15], or Hall and Huang [16], among many others.

However, for regression functions of more than one variable, the problem of computing isotonic regression estimators is substantially more difficult, and good algorithms are not available except for special cases. Most of the literature deals with algorithms for multivariate isotonic regression on a grid (see Gebhardt [17], Dykstra and Robertson [18], Lee [19], or Qian and Eddy [20], among others), but little work has been done for the case of random design.

3.3. Applications

In many applications of economics, operations research, biology, and medical sciences, monotonicity is often assumed to be hold or is mathematically proven to hold. For example, in economics, economists often assume the downward slope of the demand as a function of the price. Also, monotonicity applies to indirect utility, expenditure production, profit, and cost functions; see Gallant and Golub [21], Matzkin [22], or Aït-Sahalia and Duarte [23], among many others. In operations research, the steady-state mean waiting time of a customer in a single-server queue is proven to be nonincreasing as a function of the service rate; see Weber [24]. In biology, the weight or height of growing objects over time are known to be nondecreasing over time. In medical sciences, the blood pressure is believed to be a monotone function of the use of tobacco and the body weight (see Moolchan et al. [25]) and the probability of contracting a cancer depends monotonically on certain factors such as smoking frequency, drinking frequency, and weight.

3.4. Challenges

Two common issues in the isotonic regression estimator are (i) non-smoothness (piecewise constant fits) and (ii) boundary overfitting (“spiking”). Figure 1 displays an instance of the true regression function $f_0$, the $({\mathit{X}}_i,Y_i)$’s, and the isotonic regression estimator. The isotonic regression estimator in Figure 1 is piecewise constant, and hence, non-smooth. Also, it overfits near the boundary of its domain.

3.4.1. Non-Smoothness of Isotonic Regression

To overcome the non-smoothness of the isotonic regression estimator, several researchers proposed variants of the isotonic regression estimator. Ramsay [15] considered the univariate case and proposed estimating $f_0$ as a nonnegative linear combination of the I-splines given the knot sequence ${\mathit{X}}_1,\dots ,{\mathit{X}}_n$ in $\mathbb{R}$. The I-splines are monotone by definition, so the corresponding estimator is guaranteed to be monotone. This method, however, can be sensitive to the number and placement of knots and the order of the I-splines. Mammen [13] also considered the univariate case and proposed computing the isotonic regression estimator and then making it smooth by applying kernel regression. He and Shi [26] considered the univariate case and proposed fitting normalized B-splines to the data set given the knot sequence ${\mathit{X}}_1,\dots ,{\mathit{X}}_n$ in $\mathbb{R}$. This method is also sensitive to the number and placement of the knots. See Mammen and Thomas-Agnan [27] and Meyer [28] for other approaches in the univariate case.

In multiple dimensions, Dette and Scheder [29] proposed a procedure that smooths the data points using a kernel-type method and then applies an isotonization step to each coordinate so that the resulting estimate maintains the monotonicity in each coordinate.

3.4.2. Overfitting of Isotonic Regression

The isotonic regression estimator is known to be inconsistent at boundaries; this is called the “spiking” problem. See Figure 1 for an illustration of this phenomenon. Proposition 1 of Lim [30] proved that when $\mathrm{\Omega }=[a,b]$, ${\widehat{f}}_n^m\left(a\right)$ does not converge to $f_0\left(a\right)$ in probability as $n\to \infty$. One popular remedy for this is using a penalized isotonic regression estimator, which minimizes the sum of squared errors plus a penalty term over all monotone functions. A penalized isotonic regression estimator can be expressed as the solution to the following problem:

$$\begin{array}{cccc}\hfill & \underset{f\in {\mathcal{F}}_m}{\mathrm{minimize}}\hfill & \hfill & {\displaystyle \frac{1}{n}}\sum _{i=1}^n{(Y_i-f\left({\mathit{X}}_i\right))}^2+{\lambda }_n{P}_n\left(f\right),\hfill \end{array}$$

(7)

where ${\lambda }_n\ge 0$ is a smoothing parameter and $P_n(·)$ is a penalty term. Wu et al. [31] considered the case where $P_n\left(f\right)$ equals the range of f divided by n. They proved that when $d=1$, the solution to (7) (which is often referred to as the “bounded isotonic regression”) evaluated at $\mathit{x}\in \mathrm{\Omega }$ converges to $f_0\left(\mathit{x}\right)$ in probability as $n\to \infty$. This proves that the bounded isotonic regression estimator is consistent at the boundary unlike the traditional isotonic regression estimator. Wu et al. [31] also proved similar results for the $d=2$ case when the covariates are on a grid. Luss and Rosset [32] proposed an algorithm for solving the bounded isotonic regression for the multivariate case. On the other hand, Lim [30] considered the case where $P_n\left(f\right)={max}_{\mathit{x}\in \mathrm{\Omega }}\left|f^{\prime }\left(\mathit{x}\right)\right|$ for the $d=1$ case and established the strong consistency of the corresponding estimator at $\mathit{x}\in \mathrm{\Omega }$, thereby proving that the estimator is consistent at the boundary of $\mathrm{\Omega }$. Lim [30] also computed the convergence rate of $n^{-1/3}$ for the proposed estimator.

Even though several approaches are proposed for the $d=1$ case, limited work is available for the multivariate case. Isotonic regression estimator’s overfitting problem in the $d>1$ case is a good research topic for future.

4. Convex Regression

In this section, we review statistical properties, interesting applications, and some challenges of convex regression estimator. We also review some algorithms for computing the convex regression estimator numerically.

4.1. Statistical Properties

4.1.1. Univariate Case

The convex regression estimator, along with its QP formulation in (3), was first introduced by Hildreth [2] for the $d=1$ case. Hanson and Pledger [33] established strong consistency on compact subsets of $\mathrm{\Omega }$. They proved that, if $\mathrm{\Omega }=(a,b)$, for any $c>a$ and $d<b$,

$$\underset{c\le \mathit{x}\le d}{max}|{\widehat{f}}_n^c\left(\mathit{x}\right)-f_0\left(\mathit{x}\right)|\to 0$$

as $n\to \infty$ with probability one.

Results on rates of convergence were obtained by Mammen [34] and Groeneboom et al. [35]. In particular, Mammen [34] showed that at any interior point $\mathit{x}\in \mathrm{\Omega }$,

$${\widehat{f}}_n^c\left(\mathit{x}\right)-f_0\left(\mathit{x}\right)={\mathcal{O}}_p\left(n^{-2/5}\right)$$

as $n\to \infty$, provided that the second derivative of $f_0$ at $\mathit{x}$ does not vanish.

Groeneboom et al. [35] identified the limiting distribution of $n^{2/5}({\widehat{f}}_n^c\left(\mathit{x}\right)-f_0\left(\mathit{x}\right))$ at a fixed point $\mathit{x}\in \mathrm{\Omega }$. Groeneboom et al. [35] proved that if the errors are iid sub-Gaussian with mean zero and $f_0^{\prime }\left(\mathit{x}\right)\ne 0$,

$$n^{2/5}\left({\widehat{f}}_n^c\left(\mathit{x}\right)-f_0\left(\mathit{x}\right)\right)\Rightarrow {\mathbb{H}}^{\prime \prime }\left(0\right),$$

where $\mathbb{H}$ denotes the greatest convex minorant (i.e., envelope) of integrated Brownian motion with quartic drift $h^4$, and ${\mathbb{H}}^{\prime \prime }\left(0\right)$ is its second derivative at 0. The envelope can be viewed as a cubic curve lying above and touching the integrated Brownian motion $+h^4$.

More recently, Chatterjee [36] derived an upper bound on the $l_2$-risk of the convex regression estimator. Under the assumption of iid Gaussian errors with mean zero, they showed that

$$R(f_0,{\widehat{f}}_n^c)\triangleq \mathbb{E}\left[{\displaystyle \frac{1}{n}}\sum _{i=1}^n{\left({\widehat{f}}_n^c\left({\mathit{X}}_i\right)-f_0\left({\mathit{X}}_i\right)\right)}^2\right]\le C{n}^{-4/5}$$

for all $n\ge 1$, where C is a constant.

4.1.2. Multivariate Case

The multivariate convex regression estimator and its QP formulation in (5) were first introduced by Allon et al. [37]. Strong consistency over any compact subset of the interior of $\mathrm{\Omega }$ was established by Seijo and Sen [38] and Lim and Glynn [39]. Seijo and Sen [38] proved that for any compact subset A of the interior of $\mathrm{\Omega }$,

$$\underset{\mathit{x}\in A}{max}|{\widehat{f}}_n^c\left(\mathit{x}\right)-f_0\left(\mathit{x}\right)|\to 0$$

as $n\to \infty$ with probability one under some modest assumptions.

Rates of convergence were obtained by Lim [40]. In particular, Lim [40] considered a variant of the convex regression estimator, say ${\tilde{f}}_n^c(·)$, which minimizes the sum of squared errors ${\sum }_{i=1}^n{(Y_i-f\left(X_i\right))}^2$ over the set of convex functions f with subgradients bounded uniformly by a constant. Lim [40] showed that

$${\displaystyle \frac{1}{n}}\sum _{i=1}^n{({\tilde{f}}_n^c\left({\mathit{X}}_i\right)-f_0\left({\mathit{X}}_i\right))}^2=\left\{\begin{array}{cc}{\mathcal{O}}_p\left(n^{-4/(4+d)}\right),\hfill & \mathrm{if}d<4\hfill \\ {\mathcal{O}}_p\left((logn)n^{-1/2}\right),\hfill & \mathrm{if}d=4\hfill \\ {\mathcal{O}}_p\left(n^{-2/d}\right),\hfill & \mathrm{if}d>4\hfill \end{array}\right.$$

as $n\to \infty$.

A bound on the global risk has been established by Han and Wellner [41] as follows:

$$\mathbb{E}\left[{\int }_{\mathrm{\Omega }}{\left({\overline{f}}_n^c\left(\mathit{x}\right)-f_0\left(\mathit{x}\right)\right)}^2\lambda \left(\mathit{x}\right)d\mathit{x}\right]=\left\{\begin{array}{cc}O\left(n^{-4/(4+d)}\right),\hfill & \mathrm{if}d<4\hfill \\ O\left((logn)n^{-1/2}\right),\hfill & \mathrm{if}d=4\hfill \\ O\left(n^{-2/d}\right),\hfill & \mathrm{if}d>4\hfill \end{array}\right.$$

as $n\to \infty$, where $\lambda (·)$ is the density function of ${\mathit{X}}_1$ and ${\overline{f}}_n^c(·)$ minimizes the sum of squared errors ${\sum }_{i=1}^n{(Y_i-f\left(X_i\right))}^2$ over the set of convex functions f bounded uniformly by a constant.

The limiting distribution of $r_n^{-1}({\widehat{f}}_n^c\left(\mathit{x}\right)-f_0\left(\mathit{x}\right))$ for an appropriate rate of convergence $r_n$ has not be established so far and this is still an open problem.

4.2. Computational Algorithms

In the univariate case, Dent [42] introduced a QP formulation of the convex regression estimator as in (3). Dykstra [43] developed an iterative algorithm to compute the convex regression estimator numerically. The procedure proposed by Dykstra [43] is guaranteed to converge to the convex regression estimator as the number of iterations approaches infinity. For other algorithms for solving the QP in (3), see Wu [44] and Fraser and Massam [45].

For the multivariate setting, Allon et al. [37] first presented the QP formulation given in (5). However, this formulation requires $O\left(n^2\right)$ linear constraints, which limits its applicability to data sets of only a few hundred observations under current computational capabilities. To address this, Mazumder et al. [46] proposed an algorithm based on the augmented Lagrangian method. Their approach was able to solve (5) for $n=1000$ and $d=10$ within a few seconds.

To mitigate the computational burden, several authors (Magnani and Boyd [47]; Aguilera et al. [48]; Hannah and Dunson [49]) proposed approximating the convex regression estimator by the maximum of a small set of hyperplanes. Specifically, they solved the optimization problem

$$({\alpha }^{*},{\beta }^{*},K^{*})={argmin}_{\alpha \in \mathbb{R},\beta \in {\mathbb{R}}^d,K\in \mathbb{N}}\sum _{i=1}^n{\left[Y_i-\underset{k=1,\dots ,K}{max}\left({\alpha }_k+{\beta }_k^{\mathsf{T}}{\mathit{x}}_i\right)\right]}^2,$$

where $(\alpha ,\beta )\in \mathbb{R}\times {\mathbb{R}}^d$ defines a hyperplane and $\mathbb{N}$ is the set of positive integers. Their estimator ${\widehat{f}}_n(·)$ is then given by the maximum over these hyperplanes:

$${\widehat{f}}_n\left(\mathit{x}\right)=\underset{k=1,\dots ,K^{*}}{max}\left({\alpha }_k^{*}+{\beta }_k^{*\mathsf{T}}\mathit{x}\right).$$

More recently, to handle larger data sets, Bertsimas and Mundru [50] introduced a penalized convex regression estimator defined as the solution to the following QP problem in decision variables $f_1,\dots ,f_n\in \mathbb{R}$ and ${\beta }_1,\dots ,{\beta }_n\in {\mathbb{R}}^d$:

$$\begin{array}{cccc}\hfill & \mathrm{minimize}\hfill & \hfill & {\displaystyle \frac{1}{n}}\sum _{i=1}^n{(Y_i-f_i)}^2+{\displaystyle \frac{{\lambda }_n}{n}}\sum _{i=1}^n{\parallel {\beta }_i\parallel }^2\hfill \\ \hfill & \mathrm{subject}\mathrm{to}\hfill & \hfill & f_i\ge f_j+{\beta }_j^{\mathsf{T}}({\mathit{X}}_i-{\mathit{X}}_j),1\le i,j\le n,\hfill \end{array}$$

(8)

where ${\lambda }_n>0$ is a smoothing parameter. They proposed a cutting-plane algorithm for solving (8) and reported that it can handle data sets with $n=10000$ and $d=10$ within minutes.

Another strategy to alleviate the computational burden of the convex regression estimator is to reformulate the problem as a linear programming (LP) problem rather than a quadratic program. Luo and Lim [51] proposed estimating the convex function that minimizes the sum of absolute deviations (instead of the sum of squared errors). The resulting estimator can be obtained by solving the following LP problem in decision variables $f_1,\dots ,f_n\in \mathbb{R}$ and ${\beta }_1,\dots ,{\beta }_n\in {\mathbb{R}}^d$:

$$\begin{array}{cccc}\hfill & \mathrm{minimize}\hfill & \hfill & {\displaystyle \frac{1}{n}}\sum _{i=1}^n|Y_i-f_i|\hfill \\ \hfill & \mathrm{subject}\mathrm{to}\hfill & \hfill & f_i\ge f_j+{\beta }_j^{\mathsf{T}}({\mathit{X}}_i-{\mathit{X}}_j)for1\le i,j\le n.\hfill \end{array}$$

(9)

Luo and Lim [51] further showed that the corresponding estimator converges a.s. to the true function as $n\to \infty$ over any compact subset of the interior of $\mathrm{\Omega }$.

4.3. Applications

Convexity and concavity play an important role across various fields, including economics and operations research. In economics, for instance, it is common to assume diminishing marginal productivity with respect to input resources, which implies that the production function is concave in its inputs (Allon et al. [37]). Also, demand (Varian [52]) and utility (Varian [53]) functions are often assumed to be convex/concave. In operations research, the steady-state waiting time in a single-server queue has been shown to be convex in the service rate; see, for example, Weber [24]. An additional example is provided in Demetriou and Tzitziris [54], where infant mortality rates are modeled as exhibiting increasing returns with respect to the gross domestic product (GDP) per capita, leading to the assumption that the infant mortality function is convex in GDP per capita. For further applications of convex regression, see Johnson and Jiang [55].

4.4. Challenges

The convex regression estimator faces two primary limitations: (i) non-smoothness (piecewise linear fits) and (ii) overfitting near the boundary of its domain (“spiking”). Figure 2 illustrates these issues by displaying the true regression function $f_0$, the observed data $({\mathit{X}}_i,Y_i)$, and the convex regression estimator. As shown, the convex regression estimator is piecewise linear and therefore non-smooth. It also exhibits boundary overfitting.

4.4.1. Non-Smoothness of Convex Regression

In the univariate setting, several approaches have been proposed to address the lack of smoothness. Mammen and Thomas-Agnan [27] developed a two-step procedure in which the data are first smoothed using nonparametric smoothing methods, followed by convex regression applied to the smoothed values. Alternatively, Aït-Sahalia and Duarte [23] suggested first estimating the convex regression estimator and then applying local polynomial smoothing to the estimator. Birke and Dette [56] proposed a different strategy: smoothing the data with a nonparametric method such as the kernel, local polynomial, series, or spline methods, computing its derivative, applying isotonic regression to the derivative, and integrating to obtain a convex estimator. Other methods for the univariate case are discussed in Du et al. [57].

In the multivariate case, Mazumder et al. [46] constructed a smooth approximation to the convex regression estimator that remains uniformly close to the original estimator.

4.4.2. Overfitting of Convex Regression

As with isotonic regression, convex regression tends to overfit near domain boundaries, often producing extremely large subgradients and inconsistent estimates of $f_0$ near the boundary. For $d=1$, Ghosal and Sen [58] showed that the convex regression estimator fails to converge in probability to the true value at boundary points and that its subgradients at the boundary are unbounded in probability as $n\to \infty$.

A common strategy to mitigate overfitting is to constrain the subgradients. Lim [40] proposed bounding the ${\parallel ·\parallel }_{\infty }$ norm of the subgradients, imposing $\parallel {\beta }_i{\parallel }_{\infty }\le C$ for all i for some fixed $C>0$ to the constraints of (5). Instead, Mazumder et al. [46] added $\parallel {\beta }_i\parallel \le C$ to the constraints of (5).

Another widely used approach is penalization, which adds a penalty term to the objective function of (1). Chen et al. [59] and Bertsimas and Mundru [50] studied a penalized multivariate convex regression estimator defined as the solution to the following QP problem in decision variables $f_1,\dots ,f_n\in \mathbb{R},{\beta }_1,\dots ,{\beta }_n\in {\mathbb{R}}^d$:

$$\begin{array}{cccc}\hfill & \mathrm{minimize}\hfill & \hfill & {\displaystyle \frac{1}{n}}\sum _{i=1}^n{(Y_i-f_i)}^2+{\displaystyle \frac{1}{n}}\sum _{i=1}^n{\parallel {\beta }_i\parallel }^2\hfill \\ \hfill & \mathrm{subject}\mathrm{to}\hfill & \hfill & f_i\ge f_j+{\beta }_j^{\mathsf{T}}({\mathit{X}}_i-{\mathit{X}}_j)for1\le i,j\le n.\hfill \end{array}$$

(10)

Lim [60] generalized this idea and proposed the following three alternative approaches:

$$\begin{array}{cccccc}\hfill & \left(\mathrm{A}\right)\hfill & \hfill & {\widehat{f}}_n^A={argmin}_{f\in {\mathcal{F}}_c}\hfill & \hfill & {\displaystyle \frac{1}{n}}\sum _{i=1}^n{(Y_i-f\left({\mathit{X}}_i\right))}^2+{\lambda }_{n}J\left(f\right)\hfill \\ \hfill & \left(\mathrm{B}\right)\hfill & \hfill & {\widehat{f}}_n^B={argmin}_{f\in {\mathcal{F}}_c}\hfill & \hfill & {\displaystyle \frac{1}{n}}\sum _{i=1}^n{(Y_i-f\left({\mathit{X}}_i\right))}^2\hfill & \hfill & subjecttoJ\left(f\right)\le u_n\hfill \\ \hfill & \left(\mathrm{C}\right)\hfill & \hfill & {\widehat{f}}_n^C={argmin}_{f\in {\mathcal{F}}_c}\hfill & \hfill & J\left(f\right)\hfill & \hfill & subjectto{\displaystyle \frac{1}{n}}\sum _{i=1}^n{(Y_i-f\left({\mathit{X}}_i\right))}^2\le s_n\hfill \end{array}$$

for some smoothing parameter ${\lambda }_n\ge 0$ and constants $u_n,s_n\ge 0$, where $J\left(f\right)$ is a penalty term measuring the overall magnitude of the subgradient of f.

In particular, Lim [60] considered the case where

$$J\left(f\right)=\underset{1\le i\le d}{max}\underset{\mathit{x}\in \mathrm{\Omega }}{sup}\left|{\displaystyle \frac{\partial f}{\partial x_i}}\left(\mathit{x}\right)\right|$$

and $\partial f/\partial x_i\left(\mathit{x}\right)$ represents the partial derivative of f with respect to the ith component of $\mathit{x}$. Lim [60] established the uniform strong consistency of ${\widehat{f}}_n^A$ and ${\widehat{f}}_n^C$ and their derivatives. Specifically, she showed that

$$\underset{\mathit{x}\in \mathrm{\Omega }}{sup}|{\widehat{f}}_n^A\left(\mathit{x}\right)-f_0\left(\mathit{x}\right)|\to 0$$

(11)

and

$$\underset{\mathit{x}\in \mathrm{\Omega }}{sup}\underset{\beta \in \partial {\widehat{f}}_n^A\left(\mathit{x}\right)}{sup}{\parallel \beta -\nabla f_0\left(\mathit{x}\right)\parallel }_{\infty }\to 0$$

(12)

as $n\to \infty$ with probability one. Analogous results were shown for ${\widehat{f}}_n^C$. These results confirm that penalization successfully eliminates boundary overfitting.

Finally, Lim [60] derived convergence rates for ${\widehat{f}}_n^A$ and ${\widehat{f}}_n^C$. For ${\widehat{f}}_n^A$, she proved

$${\displaystyle \frac{1}{n}}\sum _{i=1}^n{({\widehat{f}}_n^A\left({\mathit{X}}_i\right)-f_0\left({\mathit{X}}_i\right))}^2=\left\{\begin{array}{cc}{\mathcal{O}}_p\left(n^{-{\displaystyle \frac{4}{4+d}}}\right),\hfill & \mathrm{if}d<4\hfill \\ {\mathcal{O}}_p(n^{-1/2}logn),\hfill & \mathrm{if}d=4\hfill \\ {\mathcal{O}}_p\left(n^{-2/d}\right),\hfill & \mathrm{if}d>4\hfill \end{array}\right.$$

as $n\to \infty$, and for ${\widehat{f}}_n^C$, she established

$${\displaystyle \frac{1}{n}}\sum _{i=1}^n{({\widehat{f}}_n^C\left({\mathit{X}}_i\right)-f_0\left({\mathit{X}}_i\right))}^2=\left\{\begin{array}{cc}{\mathcal{O}}_p\left(n^{-1/2}\sqrt{logn}\right),\hfill & \mathrm{if}d<4\hfill \\ {\mathcal{O}}_p(n^{-1/2}logn),\hfill & \mathrm{if}d=4\hfill \\ {\mathcal{O}}_p\left(n^{-2/d}\right),\hfill & \mathrm{if}d>4\hfill \end{array}\right.$$

as $n\to \infty$.

5. Conclusions

Shape-restricted regression offers a way to incorporate qualitative information, such as monotonicity and convexity, into nonparametric function estimation without relying on bandwidths or other tuning parameters. In this review, we have focused on the isotonic and convex regression, and have highlighted how they arise naturally in areas such as statistics, economics, and operations research. By expressing these estimators as solutions to QP problems, we established a unified optimization framework that facilitates their numerical computation.

On the theoretical side, a substantial body of work has established strong consistency, convergence rates, and risk bounds for isotonic and convex regression. At the same time, the literature has identified important limitations of these estimators, including their inherent non-smoothness and boundary overfitting (“spiking”), which can severely distort behavior near the edges of the design space.

These shortcomings have motivated a variety of methodological extensions. For isotonic regression, smoothing via splines, kernel methods, and other post- or pre-processing steps, as well as penalization schemes such as bounded isotonic regression, have been proposed to improve smoothness and boundary behavior. For convex regression, bounding or penalizing subgradients has emerged as an effective strategy to control overfitting while preserving convexity.

6. Future Directions

This section outlines several open problems that, in our view, remain both important and underexplored. Although some of these topics have appeared in prior work, they call for substantially deeper investigation.

6.1. Problem 1. Scalable Algorithms for Large-Scale Data

Modern applications frequently involve data sets with millions of observations. For example, when a hiring manager receives an application through LinkedIn, the platform is expected to infer the applicant’s seniority level using massive amounts of historical data. Problems of this nature essentially reduce to isotonic regression with millions of data points. At present, no computational algorithm exists that can efficiently compute the isotonic regression estimator at such a scale. In the era of “Big Data,” there is an urgent need for scalable algorithms capable of handling large-scale isotonic and convex regression problems.

6.2. Problem 2. Algorithms for Multivariate Isotonic Regression

In the multivariate setting, efficient algorithms for isotonic regression under random designs remain largely unavailable.

6.3. Problem 3. Pointwise Limit Distribution of Multivariate Convex Regression

Although most major theoretical properties of shape-restricted estimators have been established over the past seven decades, one central question remains open: the pointwise limiting distribution of the multivariate convex regression estimator. Specifically, it is unknown whether

$$r_n^{-1}\left({\widehat{f}}_n^c\left(\mathit{x}\right)-f_0\left(\mathit{x}\right)\right)$$

converges weakly to a non-degenerate limit for some appropriate normalization $r_n$, and if so, what the form of that limit should be. Resolving this would close one of the most significant remaining theoretical gaps in convex regression.

6.4. Problem 4. Incorporating Smoothness into Shape-restricted Regression in Multiple Dimensions

Although several studies have explored ways to address the inherent non-smoothness of shape-restricted regression, existing methods remain preliminary. More effective techniques are needed for incorporating smoothness into multivariate isotonic and convex regression estimators.

6.5. Problem 5. Overfitting in Isotonic and Convex Regression

Only a handful of methods exist to mitigate overfitting in isotonic and convex regression. The most common strategies involve (i) adding a penalty term to the least squares objective function in (1), or (ii) imposing a hard bound on the subgradients (in the convex regression case). Both approaches introduce tuning parameters such as the smoothing constants or gradient bounds that must be chosen in practice, often arbitrarily or in a data-dependent ad-hoc manner. This raises two fundamental questions: How can these parameters be optimally determined? Are there alternative strategies for preventing overfitting without relying on externally imposed parameters?

6.6. Problem 6. Theory for Penalized Isotonic and Convex Regression

There remains substantial scope for advancing the theoretical foundations of penalized isotonic and convex regression. For penalized isotonic regression, most existing results address only the univariate case; extending these insights to higher dimensions remains open. For penalized convex regression, many fundamental questions about its theoretical behavior are still unresolved.