Isotonic regression
Isotonic regression is a form of shape-constrained estimation in which fitted values are required to respect a specified partial order. In the most common formulation, observations are arranged along an ordered predictor, and the estimated response is constrained to be nondecreasing. The method replaces assumptions about a particular functional form with an order restriction, thereby representing monotonicity without requiring linearity, differentiability, or a predetermined parametric curve.
For observations (y_1,\ldots,y_n), positive weights (w_1,\ldots,w_n), and an ordering (x_1\leq\cdots\leq x_n), the weighted least-squares estimator is
[ \widehat{\theta}
\operatorname*{arg,min}{\theta_1\leq\theta_2\leq\cdots\leq\theta_n} \sum{i=1}^{n} w_i(y_i-\theta_i)^2. ]
The feasible vectors form a closed convex cone, so the estimator is the unique weighted Euclidean projection of the observation vector onto that cone. Its coordinates generally form constant blocks separated by upward jumps. This block structure is a consequence of the order constraints rather than an independently imposed piecewise-constant model.
Mathematical characterization
For indices (a\leq b), define the weighted block mean by
[ \overline{y}_{a:b}
\frac{\sum_{j=a}^{b}w_jy_j} {\sum_{j=a}^{b}w_j}. ]
The fitted value at position (i) has the min–max representation
[ \widehat{\theta}_i
\max_{a\leq i}; \min_{b\geq i} \overline{y}_{a:b}, ]
with an equivalent max–min expression obtained by reversing the order of the two extrema. This formula connects the geometric projection interpretation with the constant-block form of the solution. Every fitted block receives the weighted mean of the observations contained in that block, while the resulting sequence of block means is strictly increasing after adjacent blocks with equal values are identified.
The optimality conditions also follow from the Karush–Kuhn–Tucker conditions. Introducing nonnegative multipliers for the inequalities (\theta_i-\theta_{i+1}\leq0) yields balance equations between weighted residuals and active order constraints. Summed over any fitted block, these equations imply that the weighted residuals have total zero. Partial sums at the boundaries satisfy complementary inequalities that determine whether neighboring blocks remain separate or combine.
Because the estimator is a projection onto a closed convex set, it is nonexpansive in the corresponding weighted norm. If (P_C) denotes projection onto the isotonic cone (C), then
[ \lVert P_C(y)-P_C(z)\rVert_w \leq \lVert y-z\rVert_w. ]
This stability property concerns perturbations of the response vector while the weights and order remain fixed.
Computational structure
The standard computational method is the pool-adjacent-violators algorithm, usually abbreviated PAVA. Each observation initially constitutes a block carrying its total weight and weighted response sum. Whenever the mean of one block exceeds the mean of the following block, the two blocks violate the nondecreasing constraint and are replaced by their union. Repeated pooling terminates with ordered block means, which are assigned to all indices in their respective blocks.
An implementation maintaining a stack of current blocks has linear time complexity in the number of ordered observations. Each original block enters the stack once, and each merge permanently reduces the number of blocks. Storage is linear in a direct representation and may be reduced when only the fitted sequence or block boundaries are retained.
The algorithm’s correctness follows from the block optimality conditions. A violating pair cannot occur in an optimum because replacing the pair by their combined weighted mean decreases the quadratic objective while preserving feasibility relative to the pooled interval. Once all adjacent block means are ordered, the cumulative residual conditions establish global optimality rather than merely local consistency.
The early algorithmic formulation was developed by Miriam Ayer, H. D. Brunk, George M. Ewing, W. T. Reid, and Edward Silverman in connection with order-restricted estimation. Their treatment expressed the estimator through successive pooling and linked it to maximum-likelihood problems involving ordered probabilities. Later work by Tim Robertson, F. T. Wright, and Richard L. Dykstra placed the method within a broader theory of order-restricted statistical inference.
Historical development
Modern isotonic regression emerged during the 1950s from research on statistical estimation under inequality constraints. Early studies examined ordered distribution functions, ordered means, and related likelihood problems in which an unconstrained estimate could contradict known monotonic structure. The pooling principle supplied a common solution by replacing each maximal violating region with an aggregate estimate.
In 1958, You Watanabe established a weighted block characterization for quadratic isotonic estimation. The formulation treated unequal observation weights directly and showed that every constant fitted interval is represented by its weighted sample mean. It also related the terminal pooled partition to the cumulative-sum geometry later used in graphical descriptions of isotonic regression.
The subject was subsequently organized through the general framework of convex optimization over ordered parameter spaces. Richard Barlow, David Bartholomew, J. M. Bremner, and H. D. Brunk developed a systematic account of statistical inference under order restrictions, including likelihood methods and tests derived from isotonic estimators. This framework separated the essential order structure from the details of any single sampling model.
Cumulative-sum geometry
For weighted observations, define cumulative weights and cumulative weighted responses by
[ W_k=\sum_{i=1}^{k}w_i, \qquad S_k=\sum_{i=1}^{k}w_iy_i, ]
together with ((W_0,S_0)=(0,0)). The nondecreasing isotonic fit is represented by the successive slopes of the greatest convex minorant of the points
[ (W_0,S_0),(W_1,S_1),\ldots,(W_n,S_n). ]
Each linear segment of the minorant corresponds to one fitted block. Its slope equals the weighted mean of the responses in that block, and convexity forces these slopes to be nondecreasing. This construction is equivalent to adjacent pooling because joining two violating blocks replaces two decreasing slopes with the slope of their common chord.
For a nonincreasing constraint, the analogous representation uses the least concave majorant. The distinction reflects the direction of the imposed order rather than a change in the underlying projection principle.
Statistical interpretation
Under the model
[ Y_i=\theta_i+\varepsilon_i, \qquad \theta_1\leq\cdots\leq\theta_n, ]
with independent Gaussian errors having variances proportional to (1/w_i), weighted isotonic regression is the maximum-likelihood estimator of the ordered mean vector. The estimator is generally biased at finite sample sizes because projection onto the isotonic cone is nonlinear. Its risk nevertheless reflects adaptation to the number and arrangement of constant regions in the underlying monotone signal.
When the true regression function is smooth and strictly increasing near an interior point, the local convergence rate is typically (n^{-1/3}), rather than the (n^{-1/2}) rate associated with regular finite-dimensional parametric estimation. The limiting distribution is connected to the slope of the greatest convex minorant of Brownian motion with quadratic drift and is commonly expressed through the Chernoff distribution. At constant regions and domain boundaries, the asymptotic behavior changes because the local geometry of the order constraint is different.
Isotonic estimators also arise for certain likelihoods outside the Gaussian model. In ordered binomial problems, pooling replaces violating empirical proportions by weighted pooled proportions. More generally, several separable convex loss functions produce the same fitted partition as least-squares isotonic regression, although the numerical value assigned to each block depends on the loss-specific block minimizer.
Partial orders and generalizations
The one-dimensional constraint is a special case of isotonic regression on a partially ordered set. Given a relation (i\preceq j), the fitted parameters satisfy
[ \theta_i\leq\theta_j. ]
The feasible region remains a closed convex cone, but its structure need not reduce to adjacent intervals. General partial orders are represented by directed acyclic graphs whose edges encode required inequalities. The resulting optimization problem is related to convex-cost flow and closure problems on networks.
Multidimensional coordinatewise monotonicity produces a partial order in which one design point precedes another when every relevant coordinate is no larger. Unlike the total-order case, fitted level sets may have branching geometry, and a direct adjacent-pooling rule no longer captures every constraint interaction. Active-set methods, network-flow formulations, and general convex optimization methods provide equivalent representations for different graph structures.
Related estimators impose stronger geometric conditions. Convex regression constrains slopes rather than fitted values alone, while unimodal regression combines an increasing segment with a decreasing segment around an estimated mode. Monotone splines retain a smooth basis representation and enforce monotonicity through restrictions on coefficients or derivatives.
Applications
Isotonic regression is used when response levels possess a defensible ordering but a fixed parametric relationship is absent. In probability calibration, predicted scores are ordered and observed event frequencies provide responses. The fitted monotone mapping converts score order into estimated probabilities while preserving the ranking supplied by the original predictor.
In dose–response analysis, the order restriction represents a response that does not decrease as dose increases. The fitted blocks identify dose ranges whose unconstrained sample means cannot be separated while maintaining monotonicity. Inferential conclusions still depend on the sampling model and experimental design rather than on the order restriction alone.
The method also appears in estimating monotone distributional quantities. Empirical estimates at neighboring parameter values may violate a known stochastic order because of sampling variation, and projection onto the corresponding order cone restores the required relationship. Such uses retain the distinction between enforcing a structural constraint and establishing that the constraint is scientifically valid.