Isotonic Regression
Fitting the best monotonic curve, via the pool-adjacent-violators algorithm.
On this page
Beginner: sometimes you know something very specific about the shape of the relationship you're trying to fit, without knowing its exact form. A drug's effect should never go DOWN as the dose goes UP; a probability model's calibrated output should never go down as the raw score goes up; a person's height-for-age curve should never decrease as age increases. In every one of these cases, the only thing you're confident about is monotonicity — the curve only ever goes up (or only ever goes down), never both. Isotonic regression takes that one assumption and nothing else, and asks: what is the best-fitting curve, in the usual squared-error sense, among ALL curves that respect it? Unlike every earlier lesson in this module, there's no line, no polynomial, no spline basis being fit — the "model" is simply "any non-decreasing step function," and the data itself, subject to that one constraint, determines the shape completely.
Intermediate: formally, isotonic regression finds the sequence of fitted values f(x₁) ≤ f(x₂) ≤ … ≤ f(xₙ) (for x sorted ascending) that minimizes the sum of squared errors against the observed y values, subject to that chain of ≤ constraints. Without the constraint, the unconstrained minimizer is trivial — just set f(xᵢ) = yᵢ for every point, zero error. The entire problem is what happens when that trivial solution violates monotonicity somewhere: you need to find the "closest" monotone sequence to the raw data, in a precise least- squares sense, and the result is always a step function — flat on stretches where enforcing monotonicity forced several points to be averaged together, jumping up (or down) at the boundaries between those stretches. The classic algorithm that finds this exact optimum, efficiently, is the pool-adjacent-violators algorithm (PAVA), derived in full below.
Advanced: isotonic regression is the natural counterpart to every other method in this module. Ridge, lasso, and elastic net constrain coefficient magnitude; splines and GAMs constrain smoothness (via a roughness penalty or a fixed basis); isotonic regression constrains order instead — and, remarkably, that single hard constraint is enough to make the problem well-posed without any tunable penalty strength at all. There's no λ to choose here: the constraint set (all non-decreasing sequences) is fixed, and the fit is simply the L2 projection of the raw y vector onto that set (made precise in the Going Deeper note below). It is also fully non-parametric in the same spirit as LOESS (2.3.14) — no fixed functional form, no polynomial degree, no number of knots to pick — but where LOESS trades that freedom for local smoothness assumed via a kernel, isotonic regression trades it for the much sparser assumption of monotonicity alone, which is why its fitted curve is a staircase rather than a smooth curve.
The xᵢ are assumed sorted ascending already (if not, sort (x, y) pairs together first — the very first pitfall below). f ranges over every possible assignment of fitted values to the n data points subject only to the non-decreasing chain of constraints — not over any fixed parametric family. For the non-increasing (monotone-decreasing) case, simply flip every ≤ in the constraint to ≥; everything below carries over by symmetry (or equivalently, negate y, fit non-decreasing, negate back).
Step 1 — the unconstrained minimizer is trivial. Drop the monotonicity constraint entirely for a moment. Minimizing Σᵢ(yᵢ − f(xᵢ))² with absolutely no restriction on f is n completely independent one-variable minimizations — the cross terms don't exist, each f(xᵢ) only appears in its own squared term. The minimizer of a single term (yᵢ − f(xᵢ))² is obviously f(xᵢ) = yᵢ, achieving zero error. So the unconstrained optimum is just "return the data unchanged" — and the entire content of isotonic regression is what to do when that unconstrained optimum isn't monotone.
Step 2 — two adjacent violators should be pooled to their weighted average. Suppose f(x_i) > f(x_{i+1}) somewhere in the candidate solution — a violation of the constraint. Claim: at the true optimum, whenever two adjacent fitted values would otherwise want to violate monotonicity, the optimal fix is to force them to be equal, and specifically equal to the weighted average of the y values they represent. Here's why. Consider any two adjacent blocks of points currently fit to values a and b with a > b, covering w_a and w_b points respectively (weight = point count, if every point is unweighted so far). Fixing them at any two DIFFERENT values with a' ≤ b' to satisfy the constraint is provably worse than fixing them at one shared value equal to their weighted average — because for a fixed sum of two numbers, the sum of two weighted squared-error terms w_a(y̅_a − v)² + w_b(y̅_b − v)² (where y̅_a, y̅_b are each block's own average y) is a convex parabola in v, minimized exactly at:
— the weighted average — and any other choice of a single shared value increases the sum. Since the unconstrained per-block optima (a and b separately) already violate the ordering, and any two valid values with a' < b' would only be even further from each block's own unconstrained optimum than the shared value v* is, pooling to the weighted average strictly dominates every alternative that keeps the two blocks apart. This is the one-line engine behind PAVA: scan left to right, and the instant you find f(xᵢ) > f(x_{i+1}), replace both points with their weighted average, merging them into a single block.
Step 3 — merging can create a new violation, so you may need to pool again. After merging blocks i and i+1 into one block at the shared average value, that new value might now be GREATER than the block immediately to its left — a fresh violation that didn't exist before the merge. The algorithm handles this by stepping back and re-checking the boundary just behind the newly merged block, pooling again if needed, and repeating until no adjacent pair violates the ordering anywhere in the sequence. Because every merge strictly reduces the number of blocks (by at least one) and the number of blocks is finite and bounded below by 1, this process is guaranteed to terminate, and it terminates at a sequence with no violations left — i.e., a valid monotone solution. Because every single merge along the way was shown in Step 2 to be the exact, provably optimal fix for that specific violation (never merely a heuristic patch), and because the squared-error objective is convex over a convex constraint set (chains of linear inequalities carve out a convex polytope), this greedy, local procedure reaches the global optimum — not just a good monotone approximation.
Worked numeric trace. Take six points with y = [3, 1, 4, 2, 5, 6] at x = 1..6. Scan left to right: 3 > 1 — violation, pool into [2, 2] (weight 2 each, average 2). Now the sequence reads [2, 2, 4, 2, 5, 6] (the pooled block shown at both original positions) — check the boundary just behind: no block before it, move on. Next, 4 > 2 — violation, pool 4 and 2 into a block of weight 2 at value 3: sequence is now [2, 2, 3, 3, 5, 6]. But now the boundary just behind that new block needs rechecking: is the block at 2 (weight 2) greater than the new block at 3 (weight 2)? No, 2 ≤ 3, no further violation there. Continue scanning right: 3 ≤ 5, fine; 5 ≤ 6, fine. Final result: ŷ = [2, 2, 3, 3, 5, 6] — a valid monotone step function that pools indices 1-2 into value 2 and indices 3-4 into value 3, leaving the last two points untouched because they were already consistent with everything before them.
Where this is used: PAVA runs in O(n) time using a simple stack-based implementation (push each new point, and while the top two blocks on the stack violate the order, pop and merge — precisely the two-pointer/stack pattern implemented in the code tabs below), making isotonic regression cheap enough to run as a routine post-processing step on model outputs, not just a standalone fitting method — exactly the calibration use case in the second diagram below.
Grey dots are the raw, noisy (but genuinely increasing on average) data. The red step function is the current partial fit; red dots mark each block's current pooled value. The amber band highlights the next adjacent pair that violates monotonicity and is about to be pooled into their weighted average. It plays through automatically on load — step through manually with the slider afterward.
A synthetic classifier's raw scores are a monotone but distorted (overconfident) function of the true probability. The red reliability curve shows how often events actually occurred within each raw-score bin — far from the grey diagonal means the raw scores can't be trusted as probabilities. Isotonic regression, fit directly on (raw score, outcome) pairs via PAVA, produces the green calibrated curve, which tracks the diagonal far more closely at every distortion level. Drag the slider to make the raw classifier more or less overconfident.
The from-scratch tabs implement PAVA with an explicit stack of pooled blocks — push a new point, then merge backward while the top two blocks violate the order — exactly the O(n) algorithm derived above, no library involved. The library tab shows the same idea via sklearn.isotonic.IsotonicRegression, including the two options a from-scratch version leaves you to handle yourself: increasing for the non-increasing case, and out_of_bounds for what to do at query points outside the training range.
#include <cstdio>
#include <vector>
#include <cmath>
struct Block {
double value;
double weight;
int size;
};
// Pool Adjacent Violators Algorithm -- O(n) via a stack of pooled blocks.
// Assumes y is already sorted by x.
std::vector<double> isotonicRegression(const std::vector<double>& y) {
std::vector<Block> stack;
for (double yi : y) {
stack.push_back({yi, 1.0, 1});
while (stack.size() > 1 && stack[stack.size() - 2].value > stack.back().value) {
Block b = stack.back(); stack.pop_back();
Block a = stack.back(); stack.pop_back();
double mergedWeight = a.weight + b.weight;
double mergedValue = (a.value * a.weight + b.value * b.weight) / mergedWeight;
stack.push_back({mergedValue, mergedWeight, a.size + b.size});
}
}
std::vector<double> fitted(y.size());
int pos = 0;
for (const Block& blk : stack) {
for (int k = 0; k < blk.size; ++k) fitted[pos++] = blk.value;
}
return fitted;
}
int main() {
std::vector<double> y = {3, 1, 4, 2, 5, 6};
std::vector<double> fitted = isotonicRegression(y);
std::printf("raw: ");
for (double v : y) std::printf("%.2f ", v);
std::printf("\nfitted: ");
for (double v : fitted) std::printf("%.2f ", v);
std::printf("\n");
bool monotone = true;
for (size_t i = 1; i < fitted.size(); ++i)
if (fitted[i] < fitted[i - 1] - 1e-9) monotone = false;
std::printf("monotone non-decreasing: %s\n", monotone ? "true" : "false");
return 0;
}- Calibrating classifier probabilities. A classifier's raw output (a softmax score, an SVM decision-function value, a tree ensemble's leaf average) is often systematically over- or under-confident even when it ranks examples correctly. Isotonic regression fit on (raw score, actual outcome) pairs — exactly the second diagram above — is one of the two standard calibration methods in practice, alongside Platt scaling (logistic calibration).
- Dose-response curves in pharmacology. A drug's measured effect as a function of dose is biologically expected to be monotone (more drug, at least as much effect) but its exact shape — linear, saturating, sigmoidal — is often unknown or expensive to model correctly; isotonic regression estimates the response curve directly from noisy trial data without committing to any of those shapes.
- Click-through-rate (CTR) calibration in ad ranking and search. A ranking model's raw relevance or CTR score is usually monotonically related to true click probability by construction (that's the whole point of the ranking model), but the precise numeric CTR values it outputs are frequently miscalibrated; isotonic regression recalibrates them into usable probabilities for downstream bidding or budget systems.
- Monotonic feature transforms before a linear model. When a feature is known to have a monotone (but non-linear) relationship with the target, fitting an isotonic curve to that one feature and using its output as a new, transformed feature gives a downstream linear model access to the right monotone shape without hand-designing a basis for it.
- Growth curves and other biologically monotone measurements. Cumulative quantities (a child's height over age, a patient's cumulative dose received over time) are monotone by physical necessity even though measurement noise can make raw observations dip locally; isotonic regression is a natural smoother that respects that physical constraint exactly rather than approximately.
- Forgetting to sort by x first. PAVA's entire logic assumes the input sequence is already in
xorder — it pools ADJACENT violators in that order. Feeding ityvalues in the wrong order produces a "monotone" fit that's meaningless with respect to the actualxvalues. Always sort(x, y)pairs together byxbefore running PAVA (library implementations like scikit-learn's handle this internally, but a from-scratch call needs it done explicitly). - Applying it when the true relationship isn't actually monotonic. Isotonic regression will happily return a monotone fit for any data you feed it — it can't warn you that the monotonicity assumption is wrong for your problem. Forcing a genuinely non-monotone relationship (e.g. one with a clear peak or U-shape) through isotonic regression produces a systematically biased fit that flattens out exactly the structure you actually cared about.
- Overfitting via too many distinct step levels. With enough noise and enough distinct
xvalues, unconstrained PAVA can still produce a fit with almost as many blocks as there are data points, each block driven by only a few points — a jagged, overfit staircase with no smoothing at all. There's noλin plain isotonic regression to dial this back with; if that's a concern, either bin/aggregatexfirst, or move to a smooth monotone alternative (see the next point). - Confusing it with monotone smooth fits. Isotonic regression's output is a piecewise-constant step function — flat segments and sharp jumps, never a smooth curve. If a smooth monotone shape is what you actually need, the right tool is a shape-constrained spline or GAM (2.3.7, 2.3.8) fit with a monotonicity constraint, not isotonic regression, which was never trying to be smooth in the first place.
- Extrapolating past the training range. A step function has nothing defined past its last training point in either direction — a from-scratch implementation needs an explicit policy (typically holding at the nearest edge value, as shown in the library code tab's
out_of_boundsoption), and treating an isotonic fit as reliable for genuine extrapolation is a mistake regardless of which policy is chosen.
Going deeper
Isotonic regression has a clean geometric reading: it is the L2 (Euclidean) projection of the raw y vector onto the monotone cone — the set of all vectors in ℝⁿ satisfying v₁ ≤ v₂ ≤ … ≤ vₙ. That set is a closed convex cone (closed under addition and non-negative scaling, and intersection of the halfspaces vᵢ₊₁ − vᵢ ≥ 0), and projecting any point onto a closed convex set is a well-posed, unique problem — precisely why isotonic regression has one single, well-defined optimal answer with no ambiguity, and why it can be framed as a quadratic program (minimize a quadratic objective subject to linear inequality constraints) even though PAVA solves it far more efficiently than a general-purpose QP solver would. This projection view is also what makes the result exactly analogous to ridge regression's picture from earlier in this module: ridge projects (in a softer, shrinkage sense) onto a small ball of coefficient vectors; isotonic regression projects (in a hard, exact sense) onto the monotone cone. Same "closest point in a constrained set" logic, entirely different constraint set.
On the calibration side specifically, isotonic regression is a strictly more flexible alternative to Platt scaling, which calibrates by fitting a single logistic (sigmoid) curve through the raw scores. Platt scaling assumes the miscalibration has a specific sigmoidal shape and estimates only two parameters, which makes it far more data-efficient and less prone to overfitting on small calibration sets; isotonic calibration assumes only monotonicity and can fit any shape of miscalibration at all, at the cost of needing more calibration data to avoid the overfitting pitfall above. Choosing between them in practice is itself a small bias-variance trade-off, one level up from the bias-variance trade-off inside the calibrated model itself.
Finally, on complexity: naive PAVA re-scans from the start after every merge and would cost O(n²) in the worst case, but the stack-based formulation implemented in the code tabs above — where a merge only ever needs to re-check the boundary immediately behind the newly formed block, never re-scanning from the beginning — processes each point in amortized constant work, giving the full algorithm its true O(n) running time (after an initial O(n log n) sort by x, if the input isn't sorted already).
You run PAVA on y = [3, 1, 4, 2, 5, 6] (x = 1..6). Walk through which points get pooled together and give the final fitted values. Then explain, using the L2-projection view, why this result is guaranteed to be the single best possible monotone fit rather than merely a reasonable one.
Scanning left to right: 3 > 1 is a violation, so points 1 and 2 pool into a block of weight 2 at value (3+1)/2 = 2. Next compare that pooled value to point 3 (value 4): 2 <= 4, no violation there, but scanning onward, 4 > 2 (point 4) is a violation, so points 3 and 4 pool into a block of weight 2 at value (4+2)/2 = 3. Rechecking the boundary just behind this new block: is the earlier block's value 2 greater than the new block's value 3? No, 2 <= 3, fine. Continuing right: 3 <= 5 (point 5), fine, and 5 <= 6 (point 6), fine. Final fitted values: [2, 2, 3, 3, 5, 6]. This is guaranteed to be THE single best monotone fit, not merely a good one, because isotonic regression is exactly the Euclidean (L2) projection of the raw y vector onto the monotone cone -- a closed convex set -- and projection onto a closed convex set has a unique minimizer by convexity. Each PAVA pooling step was shown in the derivation to be the exact, locally optimal fix for that violation (via the weighted-average argument), and because the overall objective is convex over a convex constraint set, that sequence of exact local fixes reaches the true global optimum, not just a locally reasonable stopping point.
Isotonic regression fits the best-possible curve under nothing but a monotonicity constraint — no functional form, no smoothness assumption, no tunable λ. The key idea, derived and traced numerically above, is that whenever the unconstrained per-point optimum (just f(xᵢ) = yᵢ) violates the ordering, the exact fix is to pool the violating points into their weighted average — repeating this until no violation remains is the O(n) pool-adjacent-violators algorithm, and it is provably exact because isotonic regression is precisely the L2 projection of the data onto the convex monotone cone. This lesson closes out the Regression module: every method in it, from OLS through ridge and lasso's penalized coefficients, through splines and GAMs' smooth basis expansions, through LOESS's fully local fits, to isotonic regression's pure order constraint, is the same underlying question — what extra structure do you believe about the relationship, and how does encoding that belief change the fit — answered a different way each time.