---
title: 'Isotonic Regression: Theory & Extensions'
url: https://www.emergentmind.com/topics/isotonic-regression
type: topic
---

# Isotonic Regression: Theory & Extensions

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 [1503.08458]. 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 [2605.14943].

## 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
\]
[2605.14943].

In the discrete univariate setting, the parameter space is the monotone cone
\[
\mathcal C_n^+=\{u\in\mathbb R^n:u_1\le u_2\le \cdots \le u_n\},
\]
and the isotonic estimator is the metric projection of the ordered response vector onto \(\mathcal C_n^+\) [1303.4288]. More generally, if \(V=\{1,\dots,m\}\), \(E\subseteq V\times V\) is a loop-free directed graph, and \(w^1,\dots,w^m\) are positive weights, the isotonic regression cone is
\[
K_E^w=\{x\in\mathbb R^m:\sqrt{w^i}\,x^i\le \sqrt{w^j}\,x^j\ \text{for every }(i,j)\in E\}
\]
[1503.08458].

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 [2605.14943]. 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 \(y_{(i)}\), PAVA scans adjacent blocks and merges any pair violating monotonicity, replacing them by their weighted average,
\[
v_j'=\frac{|B_j|\,v_j+|B_{j+1}|\,v_{j+1}}{|B_j|+|B_{j+1}|},
\]
and repeats until all adjacent blocks are ordered [1303.4288]. PAVA terminates in at most \(n\) pooling steps, has worst-case complexity \(O(n)\), and computes the unique Euclidean projection onto the monotone cone [1303.4288].

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 [2605.14943]. The cone viewpoint gives a unified treatment of weighted monotone regression, graph-isotonic regression, and projection algorithms [1503.08458].

The geometric structure is especially explicit in the projection-theoretic literature. A closed convex set \(C\subseteq \mathbb R^m\) is \(\mathbb R_+^m\)-isotonic if coordinate-wise ordering is preserved by metric projection:
\[
x\le_{\mathbb R_+^m} y \implies P_Cx\le_{\mathbb R_+^m} P_Cy.
\]
A central consequence is that every isotonic regression cone \(K_E^w\) is an \(\mathbb R_+^m\)-isotonic projection set, so its metric projection preserves the coordinate-wise order [1503.08458].

Several equivalent characterizations sharpen this picture. In univariate least squares, the estimator also has the well-known min–max formula
\[
\hat g(z_\ell)
=
\min_{j\ge \ell}\max_{i\le j}\frac{1}{j-i+1}\sum_{k=i}^j y_k
=
\max_{i\le \ell}\min_{j\ge i}\frac{1}{j-i+1}\sum_{k=i}^j y_k,
\]
and can be described as the slope sequence of the greatest convex minorant of the cumulative-sum diagram [1904.04761]. 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 [1812.04249].

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 \(\ell_p\) norm, \(1\le p\le \infty\) [1706.01852].

## 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 \(L_2\) and pointwise, and at fixed points it converges at rate \(n^{-1/3}\) with a non-Gaussian limit law [1303.4288]. The same paper also gives a generalized \(L_2\)-consistency result when the effective regression target is not monotone, by projecting the true regression function onto the cone of monotone functions in \(L_2(\mu)\) [1303.4288].

Risk theory changes substantially in higher dimensions. For block-increasing functions on \([0,1]^d\), the least-squares estimator achieves the minimax rate of order \(n^{-\min\{2/(d+2),1/d\}}\) in empirical \(L_2\) loss, up to poly-logarithmic factors; for fixed lattice design, the worst-case rate is \(n^{-1/d}\) up to a logarithmic factor when \(d\ge 3\) [1708.09468]. The same work establishes a sharp oracle inequality showing that when the true function is piecewise constant on \(k\) hyperrectangles, the estimator adapts at rate \((k/n)^{\min(1,2/d)}\), again up to poly-logarithmic factors [1708.09468].

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 \(n^{-1/d}\) in \(\ell_2\) risk for \(d\ge 3\) when the range of \(f\) is bounded, and achieves the near parametric adaptation rate
\[
(K/n)\{1\vee\log(n/K)\}^{d}
\]
when \(f\) is \(K\)-piecewise constant [1812.08944]. The same estimator also has an oracle property in variable selection: when \(f\) depends only on a subset \(S\) of variables, its \(\ell_2\) risk automatically achieves, up to a poly-logarithmic factor, the minimax rate based on oracular knowledge of \(S\) [1812.08944].

For the constant-signal case, isotonic least squares has an exact distribution-free risk identity under exchangeable noise. If \(\theta^*=c\mathbf 1_n\), then
\[
R(\widehat\theta,\theta^*)
=\frac{\sigma^2}{n}\,\delta_n(\mu),
\]
and for i.i.d. zero-mean unit-variance noise,
\[
\delta_n(\mu)=H_n=1+\frac12+\cdots+\frac1n
\]
[1812.04249]. 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 \(T\) 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 [1904.04761]. The same work emphasizes that simultaneous optimality is unattainable in unimodal regression, despite its close connection [1904.04761].

Generalized isotonic regression replaces squared error by a separable convex differentiable loss,
\[
f(\hat y)=\sum_{i=1}^n f_i(\hat y_i),
\]
subject to isotonicity constraints. The Generalized Isotonic Recursive Partitioning algorithm provides a recursive partitioning method for such losses, subsuming the \(\ell_2\) case, the Barlow–Brunk framework, negative Poisson log-likelihood, Huber loss, and other convex differentiable objectives [1104.1779]. Every intermediate model along the partitioning path remains isotonic, and the algorithm converges to the global optimum [1104.1779].

Model complexity control leads to reduced isotonic regression, in which the fitted monotone step function is constrained to have exactly \(b\) blocks. For weighted one-dimensional data, exact dynamic-programming algorithms compute the optimal \(b\)-step isotonic regression in \(\Theta(n+bm)\) time, where \(m\) is the number of steps of the unconstrained isotonic regression; the same algorithms also determine optimal \(k\)-means clustering of weighted one-dimensional data [1412.2844]. This addresses the longstanding criticism that unconstrained isotonic regression may overfit by using as many as \(n\) steps [1412.2844].

Another nonclassical objective is \(L_0\) 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 \(L_p\) objectives, for penalized criteria combining \(L_p\) and weighted \(L_0\) error, and for multidimensional coordinate-wise orderings; the paper reports a reduction from previous \(\Theta(n^3)\) algorithms to \(o(n^{3/2})\) time [2106.00279].

## 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 [1102.5496]. The method was motivated in part by higher-dimensional overfitting and was applied to gene–gene interaction search and epistasis [1102.5496].

Sparsity offers another route through the high-dimensional regime. In sparse isotonic regression, the unknown monotone function depends only on an active subset \(A\subseteq\{1,\dots,d\}\) with \(|A|=s\). 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 [1907.01715]. The paper gives VC-entropy bounds, a statistical consistency regime in which \(n=\omega(s\ln d)\) and \(n=\exp\{\omega(s^2)\}\), and support-recovery guarantees for the linear-programming stage [1907.01715].

Two latent-structure problems replace the standard paired observation model. In uncoupled isotonic regression, one observes only the unordered sets \(\{x_1,\dots,x_n\}\) and \(\{y_1,\dots,y_n\}\). The problem can be reformulated through push-forward measures \(\pi_g=\frac1n\sum_i\delta_{g(x_i)}\) and solved by minimum-Wasserstein deconvolution,
\[
\hat f\in\arg\min_{g\in F_V}W_2^2(\pi_g*D,\hat\pi),
\]
leading to the minimax rate \((\log\log n)/(\log n)\), which is exponentially slower than the classical \(n^{-1/3}\) rate [1806.10648].

In isotonic regression with unknown permutations, the grid is permuted along each coordinate by latent permutations \(\pi_j\). For \(d\ge 3\), the bounded permuted class has worst-case minimax risk of order \(n^{-1/d}\) 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 [2009.02609]. The same work shows that a statistical-computational gap emerges for adaptation under an average-case hardness conjecture [2009.02609].

## 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 \(\partial f/\partial x_r\), shadow prices, or elasticities [2605.14943]. 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 [2605.14943]. 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 \(0.0736\to 0.0476\to 0.0374\) for INLS to \(0.0481\to 0.0424\to 0.0369\) for CC-INLS as \(n\) grows, and MSE reductions exceeding \(40\)–\(50\%\) in multivariate high-noise settings [2605.14943]. In an application to agglomeration economies in Finnish municipalities, out-of-sample RMSE was \(2.5191\) for CC-INLS(vertices) and \(2.5651\) for CC-INLS(centroids) [2605.14943].

Isotonic ideas also extend from conditional means to full conditional distributions. Isotonic Distributional Regression estimates conditional CDFs under stochastic-order isotonicity,
\[
x_i\preceq x_j \Longrightarrow F_{x_i}(y)\ge F_{x_j}(y)\quad\forall y,
\]
and minimizes average CRPS over isotonic collections of CDFs [1909.03725]. 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 [1909.03725]. The method is unique, calibrated, admits interpolation at new covariates, and can be combined with subagging for smoother fits and computational gains [1909.03725].

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 \(c\)-cyclically monotone sets [2603.10452]. 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 [2603.10452].

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 [2208.05121]. 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 [2208.05121].

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 [2605.14943].

Source: https://www.emergentmind.com/topics/isotonic-regression