Least Squares Isotonic REgression on Partially Ordered Sets

John Spouge, Honghui Wan, and W. John Wilbur

National Center for Biotechnology Information

National Library of Medicine

National Institutes of Health, Bethesda, MD, U.S.A.

Correspondence to: W. J. Wilbur, National Library of Medicine, Bldg.
38A, Rm. 8S806, 8600 Rockville Pike, Bethesda, MD 20894, U.S.A. Phone
301-496-2475. Fax 301-480-9241. Email wilbur@ncbi.nlm.nih.gov

Running title: Order Constrained Least Squares

# Abstract

Given a finite partially ordered set with a positive weighting function
defined on its points, it is well known that any real valued function
defined on the set has a unique best order preserving approximation in
the weighted least squares sense. Many algorithms have been given for
the solution of this isotonic regression problem. Most such algorithms
are either clearly not polynomial algorithms or they are of unknown time
complexity. It has been asserted in the literature that the general
isotonic regression problem can be solved in time *n*^4^ where *n* is
the number of points in the set but the proof given in the literature is
incorrect.

## Introduction.

The least squares isotonic regression problem has a long history and
methods for its solution have numerous applications to order restricted
statistical inference. The area is well summarized in \[Barlow et al.
1972; Qian 1994; Robertson et al. 1985; Robertson et al. 1988\]. A large
number of different algorithms have been put forward as solutions to the
general isotonic regression problem. Most of these algorithms are either
of unknown time complexity or are not polynomial. In the special case
that the *n* points of the set form a two dimensional grid and the
ordering is the natural two dimensional ordering, Block, Qian, and
Sampson \[Block et al. 1994\] give an algorithm with time complexity
*n*^2.5^ and space complexity *n*^2^. Best and Chakravarti \[Best and
Chakravarti 1990\] refer to Maxwell and Muckstadt \[Maxwell and
Muckstadt 1985\] and Roundy \[Roundy 1986\] as providing an algorithm of
time complexity *n*^4^ for the general isotonic regression problem.
Roundy, however, relies on the algorithm as described by Maxwell and
Muckstadt \[Maxwell and Muckstadt 1985\]. On examining the purported
proof in \[Maxwell and Muckstadt 1985\] we find the sketch of a proof
that relies on the computation of a maximal cut for a directed graph
with edge capacities defined and with a source and sink defined. The
problem of finding a maximal cut for an acyclic directed graph has been
shown to be NP-complete by Farley and Proskurowski \[Farley and
Proskurowski 1982\]. Their proof is for an acyclic graph without source
and sink points, but it can be trivially converted to a proof for the
case with source and sink points. One simply attaches to any acyclic
directed graph without source and sink points distinct points *s* and
*t* that are not already in the graph and that have edges directed from
*s* to all points in the original graph and from all points in the
original graph to *t*. All these new edges are given weights of zero.
The result is an acyclic graph with source and sink points and with the
property that a maximal cut for this derived graph provides also a
maximal cut for the original graph by simply dropping *s* and *t*. Thus
the maximal cut problem for a directed graph with weighted edges and
source and sink must also be an NP-complete problem.

It is our purpose here to give a proof that the general isotonic
regression problem for a partially ordered set with *n* points is solved
by solving a series of at most *n* minimal cut problems on acyclic
directed graphs with source, sink, weighted edges, and a simple
bipartite structure. This leads to the conclusion \[Galil and Naamad
1980\] that the general isotonic regression problem can be solved by an
algorithm of time complexity *n*^4^ and space complexity *n*^2^.

## II. Preliminaries.

Our concepts and notation will be based on the presentation in
\[Robertson et al. 1985\] where possible. Here we wish to describe the
setting and list the basic results from \[Robertson et al. 1985\] that
will be important for our discussion. We begin by defining the objects
of study.

**Definition 1.** Let *X* denote a set that has defined on it a relation
![](media/image1.wmf). Consider:

i)  ![](media/image2.wmf)(reflexivity).

ii) ![](media/image3.wmf)(symmetry).

iii) ![](media/image4.wmf)(antisymmetry).

iv) ![](media/image5.wmf) (transitivity).

> If ![](media/image6.wmf) is reflexive and transitive it is called a
> [quasi-order]{.underline}. A quasi-order that is antisymmetric is
> called a [partial order]{.underline}. A quasi-order that is symmetric
> is known as an [equivalence relation]{.underline}.

Our interest is in a finite nonempty quasi-ordered set *X* and two real
valued functions *g* and *w* defined on *X*. While the function *g* is
unrestricted the function *w* must be strictly positive and is known as
the weighting function. We seek to approximate *g* on *X* by functions
on *X* that preserve the ordering.

**Definition 2.** A real valued function *f* defined on *X* is
[isotonic]{.underline} iff when ![](media/image7.wmf) and
![](media/image8.wmf) then ![](media/image9.wmf) (where
"![](media/image10.wmf)" has its ordinary meaning for real numbers).

**Definition 3.** Given *X*, *g*, and *w* as described, an [isotonic
regression for *g*]{.underline} is a real valued function *g*\* that is
isotonic and among all isotonic functions minimizes the sum

![](media/image11.wmf)

The set of all isotonic functions is a closed convex set and it is
evident that *g*\* is the perpendicular projection of
![](media/image12.wmf)onto this subset under the inner product
![](media/image13.wmf) and norm ![](media/image14.wmf). From this one
obtains

**Theorem 1.** There always exists a unique isotonic regression *g*\*.

**Definition 4.** Given *X*, *g*, and *w* as described, and any nonempty
set ![](media/image15.wmf) set

![](media/image16.wmf)

We will follow the common notation that for any function *f* the set of
points \[*f* = *a*\] consists of points in *X* that satisfy
*f*(*x*)=*a*. Similarly \[*f* \> *a*\] will denote the subset of points
from *X* which satisfy *f*(*x*)\>*a* and likewise for the other possible
order relations.

**Theorem 2.** For any real number *r* satisfying ![](media/image17.wmf)
we have ![](media/image18.wmf)

Proof: See \[Robertson et al. 1988\] p. 18, Thm. 1.3.5.

**Corollary 1.** For any real number *c*

![](media/image19.wmf)

**Definition 5.** A subset *U* of *X* is called an [upper
set]{.underline} if ![](media/image20.wmf). A [lower set]{.underline} is
defined in a dual manner. The collection of all upper sets for *X* will
be denoted by ![](media/image21.wmf) and likewise the collection of all
lower sets will be denoted by ![](media/image22.wmf).

We will follow \[Robertson et al. 1985\] by generally denoting an
arbitrary member of ![](media/image23.wmf) by *L* and an arbitrary
member of ![](media/image24.wmf) by *U*. As pointed out in \[Robertson
et al. 1985\], the upper and lower sets possess a number of simple but
important properties. Both collections are closed under arbitrary
intersections and unions over members and both contain the empty set.
Further if *f* is any isotonic function on *X* then for any real number
a, ![](media/image25.wmf) is a lower set and ![](media/image26.wmf) is
an upper set.

Before moving on we would like to make one simplification in the setting
in which we work. If ![](media/image27.wmf) is a quasi-order one can
define an equivalence relation ![](media/image28.wmf) in terms of
![](media/image29.wmf), namely,

![](media/image30.wmf)

Under the equivalence relation ![](media/image31.wmf) the set *X* breaks
up into a collection ![](media/image32.wmf) of nonempty equivalence
classes. It is evident that an isotonic function on *X* is constant on
any equivalence class defined by ![](media/image33.wmf). Further the
quasi-order ![](media/image34.wmf) on *X* carries over in a canonical
way to a partial order on ![](media/image35.wmf). Then for any
equivalence class ![](media/image36.wmf) we may define

![](media/image37.wmf)

In this way the original regression problem for *g* and *w* on *X* is
transformed to the corresponding regression problem for
![](media/image38.wmf) and ![](media/image39.wmf) on
![](media/image40.wmf)and a solution for the latter trivially becomes a
solution for the former. With this in mind we will restrict attention
for the remainder of the paper to the case when ![](media/image41.wmf)
is a partial order.

## III. Decomposition 

We continue based on the content of the previous section. Our goal is to
show how the problem of computing the isotonic regression *g*\* of a
given function *g* on *X* can be reduced to the problem of finding a
special kind of partition of *X* into nonempty upper and lower sets *U*
and *L.* For such special sets we will find that *g*\* on *L* is the
isotonic regression ![](media/image42.wmf) of *g* restricted to *L* and
*g*\* on *U* is the isotonic regression ![](media/image43.wmf). The
decomposition can then be repeated in a recursive manner with a
reduction in size and complexity at each division until the problem of
computing the regression becomes trivial on the pieces. Thus the problem
is solved by finding the appropriate partitions.

**Definition 6.** A) A pair ![](media/image44.wmf) such that *U* is an
upper set and ![](media/image45.wmf) is called a [complementary
pair]{.underline}. The pair is [trivial]{.underline} if either upper or
lower set is empty. B) A nontrivial complementary pair is a [projection
pair]{.underline} for *g* if and only if there is a real number *c* such
that

> ![](media/image46.wmf)

**Theorem 3.** If ![](media/image47.wmf) is a projection pair for *g*
then

![](media/image48.wmf)

Proof: Let ![](media/image49.wmf). Because ![](media/image50.wmf),
![](media/image51.wmf) is isotonic. Using norm notation we have

![](media/image52.wmf)

It follows from the uniqueness of ![](media/image53.wmf) that
![](media/image54.wmf). QED.

The next task is to develop an approach to finding projection pairs. Our
approach is inspired by work of Maxwell and Muckstadt \[Maxwell and
Muckstadt 1985\] on a scheduling problem in economics. We will
essentially adapt their Theorem 1 to our setting. Let *b* be any real
number and define the real valued function *a* on *X* by

![](media/image55.wmf).

Then for any set ![](media/image56.wmf) we define

![](media/image57.wmf).

**Definition 6.** A) A [maximal upper set]{.underline} *U* is an upper
set that satisfies ![](media/image58.wmf). A [minimal lower
set]{.underline} is given a dual definition. B) A complementary pair
![](media/image59.wmf) in which *U* is a maximal upper set (or
equivalently *L* is a minimal lower set) is called a
[partition]{.underline} for *g*.

**Theorem 4.** Every nontrivial partition for *g* is a projection pair
for *g*.

Proof: Let ![](media/image60.wmf) be a nontrivial partition for *g*. Set

![](media/image61.wmf)

Now note that *U* may be written as the disjoint union
![](media/image62.wmf) and since by the maximality of *U* we have
![](media/image63.wmf), we must have

![](media/image64.wmf).

But then it follows that

![](media/image65.wmf).

Equation has a dual equation applicable to the lower set *L* and the
number ![](media/image66.wmf). When these results are combined with
Corollary 1 we obtain

![](media/image67.wmf).

It follows that ![](media/image68.wmf) is a projection pair for *g*.
QED.

**Corollary 2.** If ![](media/image69.wmf) is a nontrivial partition for
*g* and ![](media/image70.wmf), ![](media/image71.wmf) on *U*.

Proof: The proof of the theorem applies to show that

![](media/image72.wmf).

But ![](media/image73.wmf) implies ![](media/image74.wmf).
Also![](media/image75.wmf) and ![](media/image76.wmf) must have the same
average over *U*. It is then necessary that ![](media/image77.wmf). QED.

**Corollary 3.** If with ![](media/image78.wmf), *g* has a trivial
partition then *g*\* is constant on *X*.

Proof: First note that with this value of *b*, ![](media/image79.wmf).
Because there is a trivial partition it must be true that *X* itself is
a maximal upper set. We may then replace *U* by *X* in the proof of
Theorem 4 and obtain from the right side of equation the relation
![](media/image80.wmf). The dual argument likewise gives
![](media/image81.wmf). QED.

The strategy for finding the isotonic regression should now be evident.
If we can find maximal upper sets in an economical manner then we can
recursively partition the space until the pieces have maximal upper sets
with zero ![](media/image82.wmf) values and trivial solutions to the
isotonic regression problem.

## IV. The General Maximal Upper Set Problem

Here we continue with a consideration of the finite poset *X* and a real
valued function *a* defined on *X* but make no other requirement on *a*.
We make use of and Definition 6 to define a maximal upper set. The
problem then is to find a maximal upper set in this general setting. To
begin we remove upper and lower sets that present no difficulty to the
construction of a maximal upper set.

**Definition 7.** Given the poset *X* and the real valued function *a*,
the [*P*-layer]{.underline} is defined by

![](media/image83.wmf)

and the [*N*-layer]{.underline} by

![](media/image84.wmf).

The [residual set]{.underline} ![](media/image85.wmf).

From the transitivity of the order relation it is immediately evident
that the *P*-layer is an upper set and the *N*-layer is a lower set. If
the residual set is empty then the *P*-layer is trivially a maximal
upper set. Thus our interest will be in those cases when
![](media/image86.wmf) is nonempty and we shall assume this is the case
in the remainder of this section. We wish to convert the important
elements in ![](media/image87.wmf) into a bipartite poset \[Trotter
1992\].

Given a nonempty ![](media/image88.wmf), set

![](media/image89.wmf)

and define a partial order on *Z* by ![](media/image90.wmf) in *X*.

Then it is clear that ![](media/image91.wmf) is a partial order on *Z*
and *Z* is a bipartite set. The function *a* is already defined on the
elements of *Z*. For any subset *S* of ![](media/image92.wmf) set

![](media/image93.wmf).

**Definition 8.** Call an upper set *U* in *Z* [canonical]{.underline}
if and only if

![](media/image94.wmf)

**Lemma 1.** Every maximal upper set in *Z* is canonical.

Proof: Let *U* be a maximal upper set in *Z*. Then
![](media/image95.wmf) is also an upper set and has a greater sum
exactly when ![](media/image96.wmf) is nonempty. Since it cannot have a
greater sum this difference must be empty and satisfied. QED.

If *U* is an upper set in *Z* define

![](media/image97.wmf)

**Lemma 2.** If ![](media/image98.wmf)is a maximal upper set for *Z*,
then ![](media/image99.wmf) is a maximal upper set for *X*.

Proof: First we argue that ![](media/image100.wmf) is an upper set for
*X*. Let ![](media/image101.wmf) and let ![](media/image102.wmf) where
![](media/image103.wmf). Clearly *y* cannot be in
![](media/image104.wmf) because *z* is not. If![](media/image105.wmf)
there is nothing to prove. We are left with the case
![](media/image106.wmf). Now ![](media/image107.wmf) implies the
existence of an ![](media/image108.wmf) with ![](media/image109.wmf).
Then ![](media/image110.wmf). If ![](media/image111.wmf) this insures
that ![](media/image112.wmf) and ![](media/image113.wmf). If
![](media/image114.wmf) then we note that any ![](media/image115.wmf)
with ![](media/image116.wmf) already satisfies ![](media/image117.wmf)
and hence ![](media/image118.wmf). From the maximality of
![](media/image119.wmf) then we conclude that
necessarily![](media/image120.wmf). Finally if ![](media/image121.wmf),
then from it follows that ![](media/image122.wmf). Thus
![](media/image123.wmf) is an upper set.

To see that ![](media/image124.wmf) is a maximal upper set we note that
if *M* is a maximal upper set in *X*, then ![](media/image125.wmf) is
also so that we may assume that ![](media/image126.wmf). Further if
![](media/image127.wmf) and ![](media/image128.wmf)
with![](media/image129.wmf), then ![](media/image130.wmf) so that
![](media/image131.wmf). It follows that ![](media/image132.wmf) is an
upper set in *Z*. Thus we can write

![](media/image133.wmf)

This establishes the maximality of ![](media/image134.wmf). QED.

Our approach to finding a canonical maximal upper set for *Z* is to
construct a network flow problem based on the bipartite structure of
*Z*. We add source and sink points ![](media/image135.wmf) to *Z* to
form the set ![](media/image136.wmf) and construct the directed graph
*G* with node set ![](media/image137.wmf) and with edges and capacities:

i)  For each ![](media/image138.wmf) an edge from
    ![](media/image139.wmf) to *x* with ![](media/image140.wmf).

ii) For each ![](media/image141.wmf) an edge from *x* to *y* with
    ![](media/image142.wmf).

iii) For each ![](media/image143.wmf) an edge from *y* to
     ![](media/image144.wmf) with ![](media/image145.wmf).

There are a number of algorithms that solve a network flow problem of
this type to give a maximal flow. We make use of the max-flow min-cut
theorem and assume the solution has the form of a minimal cut for *G*.
If ![](media/image146.wmf) the [capacity]{.underline}
![](media/image147.wmf) equals the sum of all the capacities of edges
that begin in *A* and end in *B*. A [cut]{.underline}
![](media/image148.wmf) for *G* is defined by any set
![](media/image149.wmf) with ![](media/image150.wmf) and
![](media/image151.wmf) (Note that ![](media/image152.wmf)). The
[capacity of the cut]{.underline} is ![](media/image153.wmf). A [minimal
cut]{.underline} is a cut with minimal capacity.

![](media/image154.wmf)**Lemma 3.** If ![](media/image155.wmf) with
![](media/image156.wmf) and ![](media/image157.wmf) and
![](media/image158.wmf), then

![](media/image159.wmf)

Proof: By direct calculation

![](media/image160.wmf)

Here the first equation follows from

![](media/image161.wmf)

because no edge can begin in ![](media/image162.wmf) and end outside of
![](media/image163.wmf). The second equality in relies on the fact that
![](media/image164.wmf) adds up the capacities ![](media/image165.wmf)
while ![](media/image166.wmf) adds up the values
![](media/image167.wmf). QED.

**Lemma 4.** If *S* defines a minimal cut for *G*, then
![](media/image168.wmf) is a maximal upper set for *Z*.

Proof: First we claim that

![](media/image169.wmf).

If there is a ![](media/image170.wmf) then there must be some
![](media/image171.wmf) with ![](media/image172.wmf). But then
![](media/image173.wmf) is an edge that begins in *S* and ends in
![](media/image174.wmf) and ![](media/image175.wmf), which is
impossible. We conclude that ![](media/image176.wmf). From this it
follows that

![](media/image177.wmf).

Let us set

![](media/image178.wmf).

Direct calculation shows that

![](media/image179.wmf)

Clearly if ![](media/image180.wmf) is not empty then
![](media/image181.wmf) is not a minimal cut. Thus is established.

Now it is evident from that ![](media/image182.wmf) is an upper set in
*Z* and because of , Lemma 3 applies and is a consequence. Suppose
![](media/image183.wmf) is not maximal and let *M* be a maximal upper
set in *Z*. By Lemma 1 and we have

![](media/image184.wmf).

Then we may set ![](media/image185.wmf) and note that Lemma 3 applies to
yield

![](media/image186.wmf)

By comparing and we see that we must have ![](media/image187.wmf). This
contradicts the minimality of ![](media/image188.wmf) and shows that
![](media/image189.wmf) must be a maximal upper set in *Z*. QED.

It is evident that a minimal cut for *G* provides a maximal upper set
for the original partially ordered set *X* by combining the results of
the Lemmas 2 and 4. This result combined with the results of the
previous section and with general results on network flow problems
\[Galil and Naamad 1980\] show that the general isotonic regression
problem has a time complexity not greater than *n*^4^ and the space
complexity not greater than *n*^2^. The simple bipartite structure of
the network Z\* raises the question whether it might be solvable even
more efficiently. To this question we do not know the answer.

# References

Barlow, R. E., Bartholomew, D. J., Bremner, J. M., and Brunk, H. D.
(1972). *Statistical inference under order restrictions. The theory and
application of isotonic regression.* John Wiley & Sons, New York.

Best, M. J. and Chakravarti, N. (1990). Active set algorithms for
isotonic regression; a unifying framework. *Mathematical Programming*
*47,* 425-439.

Block, H., Qian, S., and Sampson, A. (1994). Structure algorithms for
partially ordered isotonic regression. *Journal of Computational and
Graphical Statistics* *3,* 3, 285-300.

Farley, A. M. and Proskurowski, A. (1982). Directed maximal-cut
problems. *Information Processing Letters* *15,* 5, 238-241.

Galil, Z. and Naamad, A. (1980). An O(EVlog^2^V) algorithm for the
maximal flow problem. *Journal of Computer and System Sciences* *21,*
203-217.

Maxwell, W. L. and Muckstadt, J. A. (1985). Establishing consistent and
realistic reorder intervals in production-distribution systems.
*Operations Research* *33,* 6, 1316-1341.

Qian, S. (1994). Generalization of least-square isotonic regression.
*Journal of Statistical Planning and Inference* *38,* 389-397.

Robertson, T., Wright, F. T., and Dykstra, R. L. (1985). Advances in
Order Restricted Statistical Inference. Lecture Notes in Statistics, D.
Brillinger, S. Feinberg, J. Gani, J. Hartigan, and K. Krickeberg,
(Eds.), Springer-Verlag, Berlin.

Robertson, T., Wright, F. T., and Dykstra, R. L. (1988). *Order
Restricted Statistical Inference.* John Wiley & Sons, New York.

Roundy, R. (1986). A 98%-effective lot-sizing rule for a multi-product,
multi-stage production/inventory system. *Mathematics of Operations
Research* *11,* 4, 699-727.

Trotter, W. T. (1992). *Combinatorics and partially ordered sets:
dimension theory.* Johns Hopkins University Press, Baltimore.