Isotonic Regression: Theory & Extensions

• Isotonic Regression is a shape-constrained nonparametric method that forces the fitted function to be nondecreasing with respect to a given order.
• It employs exact algorithms like the Pool-Adjacent-Violators Algorithm (PAVA) and active-set methods to achieve optimal Euclidean projections onto isotonic cones.
• Extensions incorporate alternative loss functions, high-dimensional sparsity, and Bayesian priors, thereby broadening its application in complex statistical settings.

Searching arXiv for recent and foundational papers on isotonic regression and its extensions. Isotonic regression is a shape-constrained nonparametric regression problem in which the fitted function is required to be increasing with respect to a total or partial order on the covariates. In its classical least-squares form, it is the Euclidean projection of the data vector onto an isotone cone, or more generally onto an isotonic regression cone determined by order relations on the design points (Németh et al., 2015). The method is notable for imposing monotonicity without global curvature constraints, for admitting exact algorithms in both univariate and partially ordered settings, and for supporting a broad family of extensions involving alternative loss functions, distributional targets, sparsity, latent permutations, Bayesian priors, and multivariate outputs (Kuosmanen et al., 14 May 2026).

1. Classical formulation and order structure

For observations $(x_i,y_i)$, $i=1,\dots,n$, with $x_i\in\mathbb R^s$ and $y_i\in\mathbb R$, the standard model assumes

$$y_i=f(x_i)+\varepsilon_i,\qquad \mathbb E[\varepsilon_i]=0,\qquad \mathrm{Var}(\varepsilon_i)=\sigma^2,$$

with $f$ component-wise monotonic increasing: $x_i\le x_j\implies f(x_i)\le f(x_j).$ Writing $\alpha_i$ for the fitted value at $x_i$, isotonic least squares solves

$$\hat\alpha^{\rm INLS} =\arg\min_{\{\alpha_i\}}\sum_i (y_i-\alpha_i)^2 \quad\text{s.t.}\quad \alpha_i\le \alpha_j\ \text{whenever }x_i\le x_j$$

(Kuosmanen et al., 14 May 2026).

In the discrete univariate setting, the parameter space is the monotone cone

$i=1,\dots,n$0

and the isotonic estimator is the metric projection of the ordered response vector onto $i=1,\dots,n$1 (Guyader et al., 2013). More generally, if $i=1,\dots,n$2, $i=1,\dots,n$3 is a loop-free directed graph, and $i=1,\dots,n$4 are positive weights, the isotonic regression cone is

$i=1,\dots,n$5

(Németh et al., 2015).

A defining feature of the least-squares solution is that it is blockwise constant. The fitted sequence is constant on pooled observations and jumps at block boundaries, so the classical isotonic regression estimator is an increasing step function (Kuosmanen et al., 14 May 2026). This stepwise character is central to both the algorithmics and the later methodological extensions.

2. Projection geometry and exact algorithms

In the univariate case, the standard computational method is the Pool-Adjacent-Violators Algorithm (PAVA). Starting from singleton blocks with values $i=1,\dots,n$6, PAVA scans adjacent blocks and merges any pair violating monotonicity, replacing them by their weighted average,

$i=1,\dots,n$7

and repeats until all adjacent blocks are ordered (Guyader et al., 2013). PAVA terminates in at most $i=1,\dots,n$8 pooling steps, has worst-case complexity $i=1,\dots,n$9, and computes the unique Euclidean projection onto the monotone cone (Guyader et al., 2013).

For partial orders, isotonic regression is a quadratic program over a polyhedral cone, and one solves it by active-set or network-flow methods; commercial solvers such as CPLEX and Gurobi are explicitly mentioned for the general quadratic-program formulation (Kuosmanen et al., 14 May 2026). The cone viewpoint gives a unified treatment of weighted monotone regression, graph-isotonic regression, and projection algorithms (Németh et al., 2015).

The geometric structure is especially explicit in the projection-theoretic literature. A closed convex set $x_i\in\mathbb R^s$0 is $x_i\in\mathbb R^s$1-isotonic if coordinate-wise ordering is preserved by metric projection: $x_i\in\mathbb R^s$2 A central consequence is that every isotonic regression cone $x_i\in\mathbb R^s$3 is an $x_i\in\mathbb R^s$4-isotonic projection set, so its metric projection preserves the coordinate-wise order (Németh et al., 2015).

Several equivalent characterizations sharpen this picture. In univariate least squares, the estimator also has the well-known min–max formula

$x_i\in\mathbb R^s$5

and can be described as the slope sequence of the greatest convex minorant of the cumulative-sum diagram (Jordan et al., 2019). Under exchangeable increments, the slopes of that greatest convex minorant are distributed as the order statistics of the running averages, yielding an exact non-asymptotic risk formula for constant signals (Soloff et al., 2018).

The projection operator itself has additional structure beyond Euclidean nonexpansiveness. A necessary and sufficient condition for isotonic projection to be contractive under a norm is the “nonincreasing under neighbor-averaging” property, and every permutation-invariant norm satisfies it; consequently, the isotonic projection is nonexpansive in every $x_i\in\mathbb R^s$6 norm, $x_i\in\mathbb R^s$7 (Yang et al., 2017).

3. Statistical properties and risk behavior

Under monotone truth, classical isotonic regression is a consistent nonparametric estimator. In the univariate model with non-atomic design and bounded noise, the estimator is consistent in $x_i\in\mathbb R^s$8 and pointwise, and at fixed points it converges at rate $x_i\in\mathbb R^s$9 with a non-Gaussian limit law (Guyader et al., 2013). The same paper also gives a generalized $y_i\in\mathbb R$0-consistency result when the effective regression target is not monotone, by projecting the true regression function onto the cone of monotone functions in $y_i\in\mathbb R$1 (Guyader et al., 2013).

Risk theory changes substantially in higher dimensions. For block-increasing functions on $y_i\in\mathbb R$2, the least-squares estimator achieves the minimax rate of order $y_i\in\mathbb R$3 in empirical $y_i\in\mathbb R$4 loss, up to poly-logarithmic factors; for fixed lattice design, the worst-case rate is $y_i\in\mathbb R$5 up to a logarithmic factor when $y_i\in\mathbb R$6 (Han et al., 2017). The same work establishes a sharp oracle inequality showing that when the true function is piecewise constant on $y_i\in\mathbb R$7 hyperrectangles, the estimator adapts at rate $y_i\in\mathbb R$8, again up to poly-logarithmic factors (Han et al., 2017).

A distinct line of work studies alternative estimators in multi-dimensional spaces and graphs. The block estimator, defined through rectangular blocks rather than all upper and lower sets, attains the minimax rate $y_i\in\mathbb R$9 in $$y_i=f(x_i)+\varepsilon_i,\qquad \mathbb E[\varepsilon_i]=0,\qquad \mathrm{Var}(\varepsilon_i)=\sigma^2,$$0 risk for $$y_i=f(x_i)+\varepsilon_i,\qquad \mathbb E[\varepsilon_i]=0,\qquad \mathrm{Var}(\varepsilon_i)=\sigma^2,$$1 when the range of $$y_i=f(x_i)+\varepsilon_i,\qquad \mathbb E[\varepsilon_i]=0,\qquad \mathrm{Var}(\varepsilon_i)=\sigma^2,$$2 is bounded, and achieves the near parametric adaptation rate

$$y_i=f(x_i)+\varepsilon_i,\qquad \mathbb E[\varepsilon_i]=0,\qquad \mathrm{Var}(\varepsilon_i)=\sigma^2,$$3

when $$y_i=f(x_i)+\varepsilon_i,\qquad \mathbb E[\varepsilon_i]=0,\qquad \mathrm{Var}(\varepsilon_i)=\sigma^2,$$4 is $$y_i=f(x_i)+\varepsilon_i,\qquad \mathbb E[\varepsilon_i]=0,\qquad \mathrm{Var}(\varepsilon_i)=\sigma^2,$$5-piecewise constant (Deng et al., 2018). The same estimator also has an oracle property in variable selection: when $$y_i=f(x_i)+\varepsilon_i,\qquad \mathbb E[\varepsilon_i]=0,\qquad \mathrm{Var}(\varepsilon_i)=\sigma^2,$$6 depends only on a subset $$y_i=f(x_i)+\varepsilon_i,\qquad \mathbb E[\varepsilon_i]=0,\qquad \mathrm{Var}(\varepsilon_i)=\sigma^2,$$7 of variables, its $$y_i=f(x_i)+\varepsilon_i,\qquad \mathbb E[\varepsilon_i]=0,\qquad \mathrm{Var}(\varepsilon_i)=\sigma^2,$$8 risk automatically achieves, up to a poly-logarithmic factor, the minimax rate based on oracular knowledge of $$y_i=f(x_i)+\varepsilon_i,\qquad \mathbb E[\varepsilon_i]=0,\qquad \mathrm{Var}(\varepsilon_i)=\sigma^2,$$9 (Deng et al., 2018).

For the constant-signal case, isotonic least squares has an exact distribution-free risk identity under exchangeable noise. If $f$0, then

$f$1

and for i.i.d. zero-mean unit-variance noise,

$f$2

(Soloff et al., 2018). This exact formula is non-asymptotic and does not require Gaussianity.

4. Beyond least squares: functionals, losses, and secondary objectives

A distinctive feature of isotonic regression is that optimal solutions can be characterized through statistical functionals rather than a single loss. If a target functional $f$3 is defined via an identification function, then there exists a single isotonic estimator that is simultaneously optimal for every loss function consistent with that functional. The paper on optimal solutions develops this theory for the expectation, quantile, and expectile functionals, extends it from total to partial orders, and shows that any solution resulting from the pool-adjacent-violators algorithm is optimal in the total-order case (Jordan et al., 2019). The same work emphasizes that simultaneous optimality is unattainable in unimodal regression, despite its close connection (Jordan et al., 2019).

Generalized isotonic regression replaces squared error by a separable convex differentiable loss,

$f$4

subject to isotonicity constraints. The Generalized Isotonic Recursive Partitioning algorithm provides a recursive partitioning method for such losses, subsuming the $f$5 case, the Barlow–Brunk framework, negative Poisson log-likelihood, Huber loss, and other convex differentiable objectives (Luss et al., 2011). Every intermediate model along the partitioning path remains isotonic, and the algorithm converges to the global optimum (Luss et al., 2011).

Model complexity control leads to reduced isotonic regression, in which the fitted monotone step function is constrained to have exactly $f$6 blocks. For weighted one-dimensional data, exact dynamic-programming algorithms compute the optimal $f$7-step isotonic regression in $f$8 time, where $f$9 is the number of steps of the unconstrained isotonic regression; the same algorithms also determine optimal $x_i\le x_j\implies f(x_i)\le f(x_j).$0-means clustering of weighted one-dimensional data (Hardwick et al., 2014). This addresses the longstanding criticism that unconstrained isotonic regression may overfit by using as many as $x_i\le x_j\implies f(x_i)\le f(x_j).$1 steps (Hardwick et al., 2014).

Another nonclassical objective is $x_i\le x_j\implies f(x_i)\le f(x_j).$2 isotonic regression, or monotonic relabeling, where the goal is to minimize the Hamming distance to an isotonic labeling. Because there may be exponentially many optimal relabelings, secondary criteria become important. Algorithms are given for arbitrary ordinal labels, for real-valued labels with secondary $x_i\le x_j\implies f(x_i)\le f(x_j).$3 objectives, for penalized criteria combining $x_i\le x_j\implies f(x_i)\le f(x_j).$4 and weighted $x_i\le x_j\implies f(x_i)\le f(x_j).$5 error, and for multidimensional coordinate-wise orderings; the paper reports a reduction from previous $x_i\le x_j\implies f(x_i)\le f(x_j).$6 algorithms to $x_i\le x_j\implies f(x_i)\le f(x_j).$7 time (Stout, 2021).

5. High-dimensional, latent-structure, and regularized variants

Multivariate isotonic regression is statistically and computationally more difficult than the univariate case. Isotonic Recursive Partitioning constructs a regularized path of isotonic models by recursively solving “best cut” subproblems on progressively smaller groups. It converges in a finite number of splits to the exact global isotonic regression solution, while intermediate models often have better predictive performance because of complexity control (Luss et al., 2011). The method was motivated in part by higher-dimensional overfitting and was applied to gene–gene interaction search and epistasis (Luss et al., 2011).

Sparsity offers another route through the high-dimensional regime. In sparse isotonic regression, the unknown monotone function depends only on an active subset $x_i\le x_j\implies f(x_i)\le f(x_j).$8 with $x_i\le x_j\implies f(x_i)\le f(x_j).$9. The mixed-integer formulation enforces isotonicity only along those active coordinates, while a two-stage linear-programming procedure first recovers the active set and then solves ordinary isotonic regression in the recovered dimension (Gamarnik et al., 2019). The paper gives VC-entropy bounds, a statistical consistency regime in which $\alpha_i$0 and $\alpha_i$1, and support-recovery guarantees for the linear-programming stage (Gamarnik et al., 2019).

Two latent-structure problems replace the standard paired observation model. In uncoupled isotonic regression, one observes only the unordered sets $\alpha_i$2 and $\alpha_i$3. The problem can be reformulated through push-forward measures $\alpha_i$4 and solved by minimum-Wasserstein deconvolution,

$\alpha_i$5

leading to the minimax rate $\alpha_i$6, which is exponentially slower than the classical $\alpha_i$7 rate (Rigollet et al., 2018).

In isotonic regression with unknown permutations, the grid is permuted along each coordinate by latent permutations $\alpha_i$8. For $\alpha_i$9, the bounded permuted class has worst-case minimax risk of order $x_i$0 up to poly-logarithmic factors, and the Mirsky partition estimator is minimax optimal while also achieving the smallest adaptivity index possible for polynomial-time procedures (Pananjady et al., 2020). The same work shows that a statistical-computational gap emerges for adaptation under an average-case hardness conjecture (Pananjady et al., 2020).

6. Smoothing, distributional, Bayesian, and multivariate-output extensions

The step-function nature of standard isotonic regression is both a strength and a limitation. Because the fitted function is piece-constant, its gradient is zero within each block and undefined at jumps, so one cannot compute meaningful marginal quantities such as $x_i$1, shadow prices, or elasticities (Kuosmanen et al., 14 May 2026). To address this, a piece-wise linear smoothing framework fits a continuous, monotonic, piece-wise linear function to the distinct isotonic levels by constructing anchor points and enforcing “conditional convexity” through a bilevel optimization problem, reformulated as a single-level MILP (Kuosmanen et al., 14 May 2026). Monte Carlo simulations report that the method can reduce mean squared error in both convex and non-convex settings, with univariate low-noise RMSE improving from approximately $x_i$2 for INLS to $x_i$3 for CC-INLS as $x_i$4 grows, and MSE reductions exceeding $x_i$5–$x_i$6 in multivariate high-noise settings (Kuosmanen et al., 14 May 2026). In an application to agglomeration economies in Finnish municipalities, out-of-sample RMSE was $x_i$7 for CC-INLS(vertices) and $x_i$8 for CC-INLS(centroids) (Kuosmanen et al., 14 May 2026).

Isotonic ideas also extend from conditional means to full conditional distributions. Isotonic Distributional Regression estimates conditional CDFs under stochastic-order isotonicity,

$x_i$9

and minimizes average CRPS over isotonic collections of CDFs (Henzi et al., 2019). The estimator is simultaneously optimal under all proper scoring rules that are mixtures of elementary quantile or probability scores, and isotonic quantile regression and isotonic binary regression emerge as special cases (Henzi et al., 2019). The method is unique, calibrated, admits interpolation at new covariates, and can be combined with subagging for smoother fits and computational gains (Henzi et al., 2019).

For multivariate outputs, classical isotonic regression is no longer directly applicable because monotonicity is not readily extendable. Brenier isotonic regression replaces scalar monotonicity by cyclic monotonicity, using the fact that the gradient of a convex potential is cyclically monotone and that optimal transport couplings are supported on $$\hat\alpha^{\rm INLS} =\arg\min_{\{\alpha_i\}}\sum_i (y_i-\alpha_i)^2 \quad\text{s.t.}\quad \alpha_i\le \alpha_j\ \text{whenever }x_i\le x_j$$0-cyclically monotone sets (Bao et al., 11 Mar 2026). The resulting bi-level program fits a discrete Brenier potential and produces a barycentric map that is cyclically monotone; the paper demonstrates applications to multiclass probability calibration and generalized linear models (Bao et al., 11 Mar 2026).

Bayesian isotonic regression introduces prior structure on the positive first differences of the monotone signal. A locally adaptive Bayesian method assigns half shrinkage priors, especially the half-horseshoe, to those differences and uses a fast Gibbs sampler based on truncated-normal and scale-mixture updates (Okano et al., 2022). The theoretical analysis proves that posterior mean estimators are robust to large differences and that asymptotic risk for unchanged points can be improved; empirically, the half-horseshoe variant is reported to be especially effective in jump scenarios and to recover the known structural break in the Nile River flow series (Okano et al., 2022).

A broad pattern emerges across these extensions. Classical isotonic regression remains the canonical projection-based estimator for monotone structure, but recent work shows that the same order-restricted core can support marginal-effect recovery, distributional prediction, cyclically monotone multi-output regression, regularized high-dimensional fitting, and Bayesian local adaptivity without abandoning the underlying shape constraint (Kuosmanen et al., 14 May 2026).