Nearly-Isotonic Regression
Ryan J. Tibshirani∗,

Holger Hoefling†,

Robert Tibshirani‡

Abstract
We consider the problem of approximating a sequence of data points with a “nearly-isotonic”,
or nearly-monotone function. This is formulated as a convex optimization problem that yields
a family of solutions, with one extreme member being the standard isotonic regression fit. We
devise a simple algorithm to solve for the path of solutions, which can be viewed as a modified
version of the well-known pool adjacent violators algorithm, and computes the entire path inO(n log n) operations, (n being the number of data points). In practice, the intermediate fits
can be used to examine the assumption of monotonicity. Nearly-isotonic regression admits anice property in terms of its degrees of freedom: at any point along the path, the number ofjoined pieces in the solution is an unbiased estimate of its degrees of freedom. We also extend
the ideas to provide “nearly-convex” approximations.
Key words: isotonic regression, pool adjacent violators, path algorithm, degrees of freedom

1Introduction

Isotonic regression solves the following problem: given a sequence of n data points y1 , . . . yn , how
can we best summarize this by a monotone sequence β1 , . . . βn ? Formally, the problem is to findβ̂ iso = argmin
β∈Rn

n
X

(yi − βi )2 subject to β1 ≤ . . . ≤ βn .

(1)

i=1

Here and throughout we assume that “monotone” means monotone nondecreasing: an analogous
problem statement exists for monotone nonincreasing.
A unique solution to problem (1) exists and can be obtained using the pool adjacent violators
algorithm (PAVA) (Barlow et al. 1972). Roughly speaking, PAVA works as follows. We start withy1 on the left. We move to the right until we encounter the first violation yi > yi+1 . Then we
replace this pair by their average, and back-average to the left as needed, to get monotonicity. We
continue this process to the right, until finally we reach yn . If skillfully implemented, PAVA has acomputational complexity of O(n) (Grotzinger & Witzgall 1984).
There is quite a large literature on isotonic regression. The aforementioned book (Barlow et al.
1972) is a classic reference, along with Robertson et al. (1988), although there have been many
earlier references to the same idea (Brunk 1955, Ayer et al. 1955, Miles 1959, Bartholomew 1959a,
Bartholomew 1959b). A recent paper by de Leeuw et al. (2009) gives a nice overview of the problem’s
history and computational aspects. Luss et al. (2010) propose a recursive partitioning algorithm
that produces a series of less smooth isotonic fits, culminating in the usual isotonic regression.
In this paper we consider a different problem, that of approximating the data with a “nearlyisotonic” function. We formulate a convex optimization problem that yields a family of solutions,with the saturated fit at one end (β̂ = y) and the isotonic regression fit (β̂ = β̂ iso ) at the other. The∗ Dept. of Statistics, Stanford University, ryantibs@stanford.edu† Molecular Diagnostics, Novartis, hhoeflin@gmail.com‡ Depts. of Health, Research & Policy, and Statistics, Stanford University, tibs@stanford.edu

1entire path of solutions can be computed in O(n log n) operations, using a simple algorithm thatcan be thought of as a modified version of PAVA. Furthermore, we show that at any point on thepath, an unbiased estimate of the solution’s degrees of freedom is achieved by simply counting thenumber of joined pieces.

2

Nearly-isotonic fitting

Given a sequence y1 , . . . yn , we seek a nearly-monotone approximation, and so we consider the
problemnn−1X1X
β̂λ = argmin(yi − βi )2 + λ(βi − βi+1 )+ ,
(2)
2 i=1
β∈Rn
i=1
with x+ indicating the positive part, x+ = x · 1(x > 0). This is a convex optimization problem foreach fixed λ ≥ 0. The penalty term penalizes adjacent pairs that violate the monotonicity property,
that is, having βi > βi+1 . When λ = 0, the solution interpolates the data, and letting λ → ∞,we obtain the solution to the isotonic regression problem (1). For intermediate values of λ we get
nonmonotone solutions which trade off monotonicity with goodness-of-fit.
We have implicitly assumed in (2) that the data points are measured along an equally spaced
grid. If this is not the case, then it would make sense to change the penalty to
λ

n−1
X

(βi − βi+1 )+,xi+1 − xi
i=1

where xi denotes the value where the measurement yi is taken. Also, we note that the criterionin (2) is closely related to the criterion of the fused lasso signal approximator (FLSA) (Friedman
et al. 2007), which is
n−1
n
n
X
X
1X
|βi − βi+1 |.
|βi | + λ2
(yi − βi )2 + λ1
2 i=1
i=1
i=1
The first term above encourages sparsity, while the second term penalizes the jumps in adjacentcoefficients. In nearly-isotonic regression, we omit the sparsity term and penalize only nonmonotonicities.
Figure 1 shows a toy example with n = 7 data points. The solution β̂λ is shown for four valuesof λ, with the bottom right panel showing the full isotonic regression. In fact, these four values of λrepresent the knots in the entire path of solutions, with each coordinate β̂λ,i a linear function of λ inbetween these knots. The arrows indicate places where the adjacent pairs β̂λ,i , β̂λ,i+1 are joined andset equal to a common value. It is not hard to believe that, in general, if two adjacent coordinates
are joined at some value of λ, then they remain joined for all larger values of λ. This is indeed true,and it leads to a simple algorithm to compute the path of solutions, discussed in the next section.

3

The path of solutions

In this section we derive an algorithm for computing the entire path of solutions for problem (2)as a function of λ. This algorithm is a kind of “modified” version of PAVA that starts not at theleft end of the data, but rather joins adjacent points when needed, to produce the sequence of bestnear-isotonic approximations. In contrast, the sequence of intermediate fits in PAVA are not of
special interest, only the final isotonic fit.
We begin by stating a lemma suggested at the end of the last section. Essentially, the lemmasays that the pieces in the solution can only be joined together, not split apart, as λ increases. This
fact is analogous to that for the FLSA (see Section 4 of Friedman et al. 2007), and the proof givenin the Appendix, Section A.1.
2

λ=0

λ = 0.25

λ = 0.7

λ = 0.77

Figure 1: Nearly-isotonic fits for toy example. An interpolating function is shown in the top left
panel. There are three joining events (indicated by the arrows) shown in the remaining panels, withthe usual isotonic regression appearing in the bottom right panel.

Lemma 1. Suppose that for some λ, we have two adjacent coordinates of the solution satisfying
β̂λ,i = β̂λ,i+1 . Then β̂λ0 ,i = β̂λ0 ,i+1 for all λ0 > λ.
This fact greatly simplifies the construction of the nearly-isotonic solution path, and our treatment is similar to that for the FLSA in Hoefling (2010). Suppose that at a parameter value λ, there
are Kλ joined pieces, which we represent by groups of coordinates A1 , . . . AKλ . Then we can rewrite
the criterion in (2) as
Kλ X
KX
λ −1
1X
(3)
(βAi − βAi+1 )+ .
(yj − βAi )2 + λ
2 i=1
i=1
j∈Ai

To find the optimal levels β̂λ,Ai of each group Ai , we examine the subgradient equations of (3):X−
(4)
yj + |Ai |β̂λ,Ai + λ(si − si−1 ) = 0 for i = 1, . . . Kλ .
j∈AiHere si = 1(β̂λ,Ai − β̂λ,Ai+1 > 0), and for notational convenience, we let s0 = sKλ = 0.
Now suppose that the groups A1 , . . . AKλ do not change for an interval of increasing λ. Then we
can differentiate (4) with respect to λ, treating the groups Ai , and hence the values si , as constants.
This yields
si−1 − sidβ̂λ,Ai
=
,
(5)
dλ
|Ai |
3

which is itself a constant, meaning that each β̂λ,Ai is a (locally) linear function of λ.
As λ increases, equation (5) will continue to give the slope of the solution until the groups
A1 , . . . AKλ change. According to Lemma 1, this can only happen when two groups merge, that is,their solution paths intersect. Using the slopes mi = dβ̂λ,Ai /dλ derived in (5), and doing a little
algebra, we find that groups Ai and Ai+1 will merge at
ti,i+1 =β̂λ,Ai+1 − β̂λ,Ai+ λ,mi − mi+1(6)

for each i = 1, . . . Kλ − 1. Therefore we can move all the way until the “next” value of λ
λ∗ =

min

ti,i+1 ,(7)i∗ = argmin ti,i+1 .

(8)

i: ti,i+1 >λand merge groups Ai∗ and Ai∗ +1 , wherei: ti,i+1 >λNote that there may be more than one pair of groups that achieve this minimum, in which case wemerge all of the appropriate groups. The minimizations (7) and (8) are only taken over the values of
ti,i+1 that are larger than λ. If none of the ti,i+1 are larger than λ, then none of the existing groups
will ever merge [in fact, the slopes in (5) are zero], so we know the solution for the rest of the pathand the algorithm can terminate.
For the sake of clarity, we lay out the steps of the algorithm. The algorithm finds the set ofcritical points or knots {λ1 , . . . λT } at which the groups merge, and also computes the solution β̂λat each critical point.
Algorithm 1 (Modified Pool Adjacent Violators).
• Start with λ = 0, Kλ = n, and Ai = {i} for each i = 1, . . . n. The solution is β̂λ,i = yi foreach i.
• Repeat:
1. Compute the slopes mi of each group, according to (5).
2. Compute the collision times ti,i+1 of each pair of adjacent groups, according to (6).
3. If each ti,i+1 ≤ λ, terminate.
4. Compute the critical point λ∗ as in (7), and update the solution based on the slopes:β̂λ∗ ,Ai = β̂λ,Ai + mi · (λ∗ − λ)for each i = 1, . . . Kλ . Merge the appropriate groups as in (8) (so Kλ∗ = Kλ − 1), andfinally, set λ = λ∗ .
It is not difficult to verify that the path visited by the algorithm satisfies the subgradient equations(4) at each λ, and hence is indeed the solution path to nearly-isotonic regression (2). It is also
not difficult to see that the total number of critical points (or the total number of iterations) is
T = n − K∞ , where K∞ is the number of joined pieces in the full isotonic regression fit. Since thesolution path is linear between the critical points, a practical implementation could, for example,
store the solution at each of these points and then use linear interpolation to return the solution at
any other λ. See Hoefling (2010) for a more efficient tree-based implementation; when applied tothe current situation this gives a total running time of O(T log n) = O(n log n).

4

Rewriting (4) gives us the following formula for the fitted value β̂λ,Ai :P
j∈Ai yj − λ(si − si−1 )β̂λ,Ai =
,
|Ai |

(9)

where recall si = 1(β̂Ai − β̂Ai+1 > 0). Hence in any monotone increasing “stretch” of data points, theestimate is simply an average of yi values, just like in PAVA. In general, there are two key differences
between our approach and PAVA:
1. our algorithm chooses to join groups based on the tradeoff between goodness-of-fit with monotonicity, governed by the parameter λ, rather than sequentially enforcing monotonicity likePAVA;
2. our algorithm averages a shrunken version of the data points yi , as in (9), so that nonmonotone
pieces are pulled together (for example, see the top right panel in Figure 1).

4

Degrees of freedom

Suppose that the data vector y is drawn from the normal modely ∼ N (µ, σ 2 I),(10)

where I is the n × n identity matrix. The degrees of freedom of a fitted vector ŷ of y can be definedasn
1 X
Cov(ŷi , yi ).
(11)df(ŷ) = 2
σ i=1This is discussed in Efron (1986) and Hastie et al. (2008), for example. For isotonic regression, it
turns out that the number of joined pieces in β̂ iso is an unbiased estimate of this quantity. This
follows from Proposition 1 of Meyer & Woodroofe (2000), and was pointed out in Luss et al. (2010).
Fortuitously, this same property holds for the entire nearly-isotonic path. In particular, if Kλ isthe number of joined pieces in β̂λ , the solution of (2) at λ, thenndf(β̂λ ) =1 X
Cov(β̂λ,i , yi ) = E(Kλ ).
σ 2 i=1(12)

A proof of this appears in the Appendix, Section A.2.
One might ask: why is this property useful? Degrees is freedom is a measure of model complexity:
it describes the effective number of parameters that are used in the fit. In this sense, it is reassuring
to know that the number of parameters used by nearly-isotonic regression is not drastically morethan the number of joined pieces in the fit, but is actually equal to the number of joined pieces,on average. Moreover, an estimate of degrees of freedom allows us to use different model selection
criteria like Cp , AIC, or BIC. The Cp statistic, for example, is given by
Cp (λ) =nX(yi − β̂λ,i )2 − nσ 2 + 2σ 2 df(β̂λ ),

(13)

i=1

and is an unbiased estimate of the true risk E
Ĉp (λ) =

n
X

 Pn



2
i=1 (µi − β̂λ,i ) . We can define Ĉp (λ) as(yi − β̂λ,i )2 − nσ 2 + 2σ 2 Kλ ,

i=1

5

(14)

replacing df(β̂λ ) in (13) with its unbiased estimate Kλ . The modified statistic Ĉp (λ) is still unbiased
as an estimate of the true risk, which suggests choosing λ to minimize Ĉp (λ).
Note that Kλ is a step function with respect to λ, with steps at the critical points {λ1 , . . . λT }.
Pn
Meanwhile, the residual sum of squares i=1 (yi − β̂λ,i )2 is monotone nondecreasing for λ in between
critical points (this is proved in the Appendix, Section A.3). This means that Ĉp (λ) must achieve its
minimum at some λk ∈ {λ1 , . . . λT }, and hence Algorithm 1 can be used to simultaneously computethe path of solutions and select a value of the tuning parameter. We provide an example of this in
the next section.

5

Global warming data

Here we look at data on annual temperature anomalies from 1856 to 1999, relative to the 1961-1990
mean, studied in Wu et al. (2001). The data, along with an interpolating function, are shown in the
top left of Figure 2. Temperature seems to increase monotonically with possible decreases around1900 and 1950. The entire nearly-isotonic path has 139 knots, four of which are shown in Figure 2,
with the usual isotonic regression shown in the bottom right panel.
In order to try to determine whether the nonmonotonicity was real, we carried out 10-fold crossvalidation for λ, on a grid of 20 equally spaced values. This produced the CV error curve on theleft-hand side of Figure 3. The curve and its standard error bars indicate that the isotonic regression
fit (the largest value of λ on the right) has significantly worse error than the near-isotonic fits. The
one standard error rule chooses the value λ = 0.66. We also computed the modified Cp statistic, Ĉp ,given in (14). This was evaluated at each of the critical points {λ1 , . . . λT } (in the current example
we also had to use an estimate of σ 2 , which was derived from the full isotonic regression fit). The Ĉpcurve is shown on the right-hand side of Figure 3, and is minimized at λ = 0.44. This corresponds tothe top right panel in Figure 2, and is a less regularized model than that chosen by cross-validation.

6

Nearly-convex fitting

We can apply an idea similar to nearly-isotonic fitting to obtain a “nearly-convex” approximation
to a sequence of data. Namely, given y1 , . . . yn , this is
n−2
n
X
1X
2
β̂λ = argmin(βi − 2βi+1 + βi+2 )+ .
(yi − βi ) + λ
2 i=1
β∈Rn
i=1

(15)

The quantity −βi + 2βi+1 − βi+2 can be viewed as an approximation to the second derivative at
βi+1 , and so in a sense, the penalty in problem (15) encourages the sequence β to have a negative
second derivative. As another way of looking at the penalty, note that the ith term is positive when
βi+1 < (βi + βi+2 )/2, that is, βi+1 lies below the line through its neighbors. For a sufficiently large
value of λ, problem (15) produces the best convex fit to the data. For smaller values of λ it yields
a path of nearly-convex approximations that lie close to the data.
The path algorithm for this problem is more complicated than that for nearly-isotonic regression.
As λ increases, some of the differences βi − 2βi+1 + βi+2 that are negative (and hence violate
convexity), are set to zero. This is the analogue of the joining events that occur in nearly-isotonic
regression. However such differences that are zero for one value of λ can become nonzero again for
λ0 > λ, unlike the more orderly behavior guaranteed by Lemma 1 for nearly-isotonic regression.
Despite this, a path algorithm can be derived, by letting αi = βi − βi+1 for i = 1, . . . n − 1 andαn = βn , rewriting problem (15) asα̂λ = argmin
α∈Rn

n−2
n
n
2
X
X
1 X
(αi − αi+1 )+ ,αj + λ
yi −
2 i=1
i=1
j=i

6

0.2
−0.6

−0.4

−0.2

0.0

0.2
0.0
−0.2
−0.6

−0.4

Temperature anomalies

0.4

lam=0.44, rss=0.81, pieces=32

0.4

lam=0, rss=0, pieces=160

1850

1900

1950

2000

1850

1900

1950

2000

Year

−0.4

−0.2

0.0

0.2

0.4

lam=1.79, rss=1.39, pieces=22

−0.6

−0.6

−0.4

−0.2

0.0

0.2

0.4

lam=1.02, rss=1.15, pieces=26

1850

1900

1950

2000

1850

1900

1950

2000

Figure 2: Global warming data: temperature anomalies by year. An interpolating function is shownin the top left, along with the fit for three larger values of λ, and the usual isotonic regression shownin
P the bottom2 right panel. The title above each panel gives the value of λ, the residual sum of squares(yi − β̂λ,i ) , and the number of joined pieces Kλ .

7

1.5
1.0

0.012

0.5

Ĉp

0.010
0.006

0.0

0.008

CV error

0.0

0.5

1.0

1.5

0.0

λ

0.5

1.0

1.5

λ

Figure 3: Left panel: CV error curve for the global warming data, with ±1 standard error bars. The
vertical dashed line indicates the value of λ chosen by the one standard error rule. Right panel: Ĉpcurve for the same dataset, with the vertical dashed line marking the minimizing value of λ.

and using the tools developed in Hoefling (2010). Figure 4 shows a toy example.

7

Discussion

We have proposed a method for fitting a path of nearly-isotonic approximations to a sequence ofdata that culminates in the usual isotonic regression fit. We proved that the degrees of freedom ateach point in the path is simply the average number of joined pieces in the fitted function. We havealso extended the method to find nearly-convex approximations.
There are other ways that this work could be generalized. For two-dimensional data on a equally
spaced m × n grid, we can modify the nearly-isotonic regression criterion to
m−1
n
m n−1
m n
i
hX
XX
X
1 XX
(βi,j − βi+1,j )+ ,(βi,j − βi,j+1 )+ +(yi,j − βi,j )2 + λ
2 i=1 j=1
i=1 j=1i=1 j=1

which encourages monotonicity in both the horizontal (i index) and vertical (j index) directions. In
another generalization, we could control the overall smoothness of the fitted function, in addition toits montonicity. For this, we could use the criterion
n−1
n−1
n
X
X
1X
|βi − βi+1 |,(βi − βi+1 )+ + λ2
(yi − βi )2 + λ1
2 i=1
i=1
i=1

where λ1 and λ2 are both tuning parameters. For example, letting λ1 → ∞ gives a fully monotone fit,
whose smoothness can be controlled by varying λ2 . In this case (λ1 → ∞), the choice λ2 = 0 gives the

8

λ=0

λ = 0.5

λ = 0.8

λ = 3.3

λ = 6.7

λ = 10.6

Figure 4: Nearly-convex fits for a toy example. An interpolating function is shown in the top left.
The other panels show members of the nearly-convex path where a nonconvex stretch (indicated by
the arrow) is flattened into a linear segment. The best convex fit is shown in the bottom right panel.

9standard isotonic regression, while values larger than zero yield smoother monotone approximations.
This provides a different way of achieving the same goal as that discussed in Luss et al. (2010).
An R package “neariso” that implements our path algorithm for nearly-isotonic regression will
be made available on the CRAN website (R Development Core Team 2008).

Acknowledgements
The authors would like to thank Jonathan Taylor for his insights on the dual problem. The authors
also thank Saharon Rosset and his coauthors for sending us a preprint of their related work. The
third author was supported by National Science Foundation Grant DMS-9971405 and National
Institutes of Health Contract N01-HV-28183.

A

Appendix

A.1

Proof of Lemma 1

Similar to (4), but disregarding group structure, we can rewrite the subgradient equations of problem(2) as−yi + β̂λ,i + λ(si − si−1 ) = 0 for i = 1, . . . n,
(16)
where si = 1(β̂λ,i − β̂λ,i+1 > 0) if β̂λ,i − β̂λ,i+1 is > 0 or < 0, and otherwise si ∈ [0, 1]. Following
an argument of Friedman et al. (2007), suppose that we have a stretch of joined coordinates β̂λ,j =. . . = β̂λ,j+k . We show that, as λ increases, this group of remains intact until it merges with aneighboring group.
An important point is that, at the value λ, both sj−1 and sj+k are neccessarily in {0, 1}. As λincreases, sj−1 and sj+k will be constant as long as the group of coordinates {j, . . . j + k} doesn’tmerge with an adjacent one. Hence we consider a subset of the subgradient equations (16) corresponding to i = j, . . . j + k, and show that as λ increases, it has a solution with β̂λ,j = . . . = β̂λ,j+k
and sj , . . . sj+k−1 ∈ [0, 1]. This is sufficient to prove the lemma.
Taking pairwise differences between equations, and using the fact that β̂λ,j = . . . = β̂λ,j+k , weget1
As = Dy + c,
λ
where




−1
1 0 ...
0 0
2 −1
0 ... 0
0 0

 −1
0 0 
2 −1 . . . 0
0 0 , , D =  0 −1 1 . . .
A=



 ...
...
0
0 0 . . . −1 100
0 . . . 0 −1 2and s = (sj , . . . sj+k−1 ), y = (yj , . . . yj+k ), and c = (sj−1 , 0, . . . sj+k ). But A is invertible, so we canagain rewrite this as1
s = A−1 Dy + A−1 c.
(17)
λWe only need to check that s will remain in [0, 1] coordinatewise. As λ increases, the first term above
only gets smaller in magnitude. Therefore, if A−1 c is in [0, 1] coordinatewise, then the right-handside of (17) will stay in [0, 1] for increasing λ, completing the proof. As A is tridiagonal, its inverse
has an explicit form (Schlegel 1970). All we need to know is that (A−1 )i1 = (n − i + 1)/(n + 1) and(A−1 )in = i/(n + 1), and then one can check directly that A−1 c is coordinatewise in [0, 1] for all
three possibilities for c.

10

A.2

Degrees of freedom proof

In many cases, the most straightforward calculation for degrees of freedom
 comes from Stein’s
formula (Stein 1981), which states that for a function g(y) = g1 (y), . . . gn (y) of y
n
n
hX
∂gi i
1 X
Cov(g
,
y
)
=
E
.
i
i
σ 2 i=1∂yi
i=1

But Stein’s result only holds when g is continuous and almost differentiable.
Recall that if Ai , i = 1 . . . Kλ are the joined pieces in the solution β̂λ , then we can express thefitted value of each joined piece asP
j∈Ai yj − λ(si − si−1 )β̂λ,Ai =,|Ai |where si = 1(β̂λ,Ai − β̂λ,Ai+1 > 0). Suppose we take the derivative of the above expression withrespect to yj for j ∈ Ai , and simply treat the Ai ’s and si ’s as constants. Doing so would give∂ β̂λ,j /∂yj = 1/|Ai |, and assuming that we can apply Stein’s formula to β̂λ ,df(β̂λ ) = E

Kλ X
hX

i=1 j∈Ai

1 i= E(Kλ ).
|Ai |

(18)

The following lemma states that β̂λ is continuous and almost differentiable, and that for almost
every y the groups Ai and values si , i = 1, . . . Kλ , are locally constant. These conditions allow us
to precisely conclude (18).
Lemma 2. For fixed λ:(i) the solution β̂λ is continuous and almost differentiable as a function of y;(ii) for almost every y, there is a ball B containing y such that Kλ and Ai , si , i = 1, . . . Kλ are
the same for any y0 ∈ B.
Proof. The easiest proof comes from looking at the dual of problem (2). Letting D be the (n − 1) × n
matrix
−1
1 0 ...
0 0
 0 −1 1 . . .
0 0 
,
D=


...
0
0 0 . . . −1 1
we can rewrite minimization (2) as
1
β̂λ = argmin ky − βk22 + λ1T (Dβ)+ .
β∈Rn 2(19)

Here we use x+ for the positive part applied coordinate-wise to a vector x. Convex analysis tells usthat the dual of (19) is1ûλ = argmin ky − DT uk22 subject to 0 ≤ u ≤ λ.
2
n−1
u∈R

(20)

The above inequality represents a coordinate-wise inequality, 0 ≤ ui ≤ λ for each i = 1, . . . n − 1.
The solutions of (19) and (20) are related byβ̂λ = y − DT ûλ .
11

We will write β̂λ (y) and ûλ (y) to emphasize that these are functions of the data y. Studying(20), the fit DT ûλ (y) is the projection of y onto the convex set Pλ = {DT u : 0 ≤ u ≤ λ}. But a
projection onto a convex set is a contraction, and hence β̂λ (y) is Lipschitz with constant 2, sincekβ̂λ (y) − β̂λ (y0 )k2 ≤ ky − y0 k2 + kDT ûλ (y) − DT ûλ (y0 )k2 ≤ 2ky − y0 k2 .
This immediately means that β̂λ (y) is continuous and almost differentiable [for example, see Theorem
2 in Section 3.2 of Evans & Gariepy (1992)], which establishes part (i) of the lemma.
As for part (ii), it turns out that the number of pieces Kλ , as well as the sets Ai and values si ,i = 1, . . . Kλ , are determined entirely by the face of the polytope Pλ onto which y projects. For
almost every y ∈ Rn , there is a ball around y whose interior projects to a single face. The points
that don’t share this property necessarily lie on a ray that emanates orthogonally from one of thecorners of Pλ , but the union of these rays has measure zero. For more details, the reader can refer
to Section 10 of Tibshirani & Taylor (2010) or Section 2 of Meyer & Woodroofe (2000).

A.3

P
Proof that ni=1 (yi − β̂λ,i ) is monotone nondecreasing for λ in between
critical points

Pn
Consider the function f (λ) = i=1 (yi − β̂λ,i )2 , where β̂λ is the solution to nearly-isotonic regression(2). We show that f (λ) is nondecreasing for λ in between any two critical points in the solution
path. First define, for j ∈ Ai , rj = si = 1(β̂λ,Ai − β̂λ,Ai+1 > 0). Then we can rewrite (9) asP
Pj∈Ai (rj − rj−1 )j∈Ai yj − λβ̂λ,Ai =,|Ai |because the second sum above is telescoping. Therefore we can write the whole solution in vector
notation asβ̂λ = P (y − λDr),(21)
where P is the matrix that averages within each group A1 , . . . AKλ , and D is the matrix that gives
adjacent differences. As λ varies, equation (21) continues to give the solution as long as the groups
Ai and signs si do not change, which happens at the critical points. We havey − β̂λ = (I − P )y + λP Dr,
and using the fact that P is a projection matrix,ky − β̂λ k22 = y T (I − P )y + λ2 rT DT P Dr.
The coefficient multiplying λ2 is nonnegative, which completes the proof.

References
Ayer, M., Brunk, H. D., Ewing, G. M., Reid, W. T. & Silverman, E. (1955), ‘An empirical distribution function for sampling with incomplete information’, Annals of Mathematical Statistics
26(4), 641–647.
Barlow, R. E., Bartholomew, D., Bremner, J. M. & Brunk, H. D. (1972), Statistical inference underorder restrictions; the theory and application of isotonic regression, Wiley, New York.
Bartholomew, D. J. (1959a), ‘A test for homogeneity for ordered alternatives’, Biometrika 46(1), 36–
48.
12

Bartholomew, D. J. (1959b), ‘A test for homogeneity for ordered alternatives II’, Biometrika
46(3), 328–335.
Brunk, H. D. (1955), ‘Maximum likelihood estimates of monotone parameters’, Annals of Mathematical Statistics 26(4), 607–616.
de Leeuw, J., Hornik, K. & Mair, P. (2009), ‘Isotone optimization in R: Pool-adjacent-violators(PAVA) and active set methods’, Journal of Statistical Software 32(5), 1–24.
Efron, B. (1986), ‘How biased is the apparent error rate of a prediction rule?’, Journal of theAmerican Statistical Association 81(394), 461–470.
Evans, L. & Gariepy, R. (1992), Measure theory and fine properties of functions, CRC Press, Boca
Raton.
Friedman, J., Hastie, T., Hoefling, H. & Tibshirani, R. (2007), ‘Pathwise coordinate optimization’,
Annals of Applied Statistics 1(2), 302–332.
Grotzinger, S. J. & Witzgall, C. (1984), ‘Projections onto simplices’, Applied Mathematics andOptimization 12(1), 247–270.
Hastie, T., Tibshirani, R. & Friedman, J. (2008), The Elements of Statistical Learning; Data Mining,Inference and Prediction (2nd edition), Springer Verlag, New York.
Hoefling, H. (2010), A path algorithm for the fused lasso signal approximator. Unpublished.
URL: http: // www. holgerhoefling. com/ Articles/ FusedLasso. pdf
Luss, R., Rosset, S. & Shahar, M. (2010), Isotonic recursive partioning. Submitted.
Meyer, M. & Woodroofe, M. (2000), ‘On the degrees of freedom in shape-restricted regression’,
Annals of Statistics 28(4), 1083–1104.
Miles, R. E. (1959), ‘The complete amalgamation into blocks, by weighted means, of a finite set ofreal numbers’, Biometrika 46(3), 317–327.
R Development Core Team (2008), R: A Language and Environment for Statistical Computing, RFoundation for Statistical Computing, Vienna, Austria. ISBN 3-900051-07-0.
URL: http: // www. R-project. org
Robertson, T., Wright, F. T. & Dykstra, R. L. (1988), Order restricted statistical inference, Wiley,
New York.
Schlegel, P. (1970), ‘The explicit inverse of a tridiagonal matrix’, Mathematics of Computation
24(111), 665–665.
Stein, C. (1981), ‘Estimation of the mean of a multivariate normal distribution’, Annals of Statistics
9(6), 1135–1151.
Tibshirani, R. J. & Taylor, J. (2011), ‘The solution path of the generalized lasso’, Annals of Statistics
39(3), 1335–1371.
Wu, W., Woodroofe, M. & Mentz, G. (2001), ‘Isotonic regression; another look at the changepoint
problem’, Biometrika 88(3), 794–804.

13