Here is a table of the fastest known isotonic regression algorithms for several metrics and orderings. I've cited the first paper to give a correct algorithm with the given time bound, to the best of my knowledge. In some cases I've listed 2 if they appeared nearly contemporaneously. I've only include algorithms for exact calculations (to within numerical error), not approximations. Feel free to contact me about corrections or updates.

This is not an historical survey of the results leading up to the ones listed here, but rather just a list of the best known now. Many of the papers listed have extensive historical references within them. I've also omitted related topics such as unimodal regression, integer-valued regression, regressions with constraints on the number of level sets or the differences between adjacent ones, etc. I might include some of these in the future if there is sufficient interest.

A directed acyclic graph (DAG) G with n vertices V = {v₁,...,v_n} and m edges E defines a partial order over the vertices, where v_i precedes v_j if and only if there is a path from v_i to v_j. It is assumed that G is connected, and hence m ≥ n-1. If it isn't connected then the algorithms would be applied to each component independently of the others. Given G with a weighted value at each vertex (also called data), (v_i,y_i,w_i): v_i ∈ V, y_i ∈ ℜ, w_i ≥ 0}, an L_p isotonic regression of the data is a set {ý_i: ý_i ∈ ℜ, 1 ≤ i ≤ n } that minimizes

(∑i=1n wi |yi - ýi|p)1/p	if 1 ≤ p < ∞
maxi=1∞ wi |yi - ýi|	if p = ∞

among all sets satisfying the isotonic constraint that if v_i precedes v_j then ý_i ≤ ý_j The regression error is the value of this expression.

The orderings listed are linear, tree, points in multidimensional space, and general (i.e., an algorithm that applies to all orderings). For points in multidimensional space (the ``dim'' orderings), a point p=(p₁,...,p_d) in d-dimensional space is dominated by point q=(q₁,...,q_d) if p_i ≤ q_i for all 1 ≤ i ≤ d. If p is dominated by q then p precedes q, thus the regression value at p must be no larger than the value at q. The multidimensional points are further subdivided into grids and points in arbitrary positions, and into dimension 2 and dimension ≥ 3. The constants in the O-notational analysis depend on d but are not explicitly determined.

The metrics considered are L₁, L₂, and L_∞. These metrics also go under a variety of other names: L₁ is also known as median regression or Manhatten distance; L₂ is squared error regression; and L_∞ is also known as minimax optimization, uniform metric, Chebyshev distance, or least absolute deviation.

The table lists the best time known for the given ordering and metric, and citations to the references where the algorithm is described, or to a relevant note.

1. For linear orders the ``pair adjacent violators'', PAV, approach has been repeatedly rediscovered. Apparently the first paper that used this is by Ayer et al. [1]. For the L₂ metric it is trivial to implement in linear time, while for the others appropriate data structures are needed, increasing the time by a logarithmic factor.

2. The best general algorithm for L₂ is Spouge, Wan, and Wilbur's [8] correction to Maxwell and Muckstadt's [6] algorithm. Maxwell and Muckstadt claimed that their algorithm took Θ(n⁴) time, but Spouge, Wan, and Wilbur showed that it used an NP-complete subproblem. They showed how to modify the algorithm to obtain quartic time. The time is Θ(n F(n,m)), where F(n,m) is the time to compute a maximal flow on a graph with n vertices and m edges. The Sleator and Tarjan max flow algorithm has F(n,m) = Θ(nm log n), which gives the best time known when m = o(n² log n).

3. The best general algorithm for L_∞ is the improvement in [10] to the algorithm of Kaufman and Tamir [4]. This replaced a linear programming approach with one that could be computed via topological sorting, changing the time from Θ(m log n + n log² n) to Θ(m log n), which is faster for sparse graphs where m = o(n log n).

4. For 2-dimensional points with arbitrary placement, [11] shows how to simulate the 2-dimensional grid algorithms for the L₁ and L₂ metrics to obtain the indicated results.

5. For points in arbitrary position and d ≥ 3 (and d=2 for the L_∞ metric), the results are based on constructing a larger set of points for which domination ordering can be specified with a succinct set of edges, and then applying the best result for general sparse DAGS to this new set of points. This was first used in the original version of [10], and subsequently in [11]. However, the revised version of [10] replaces this approach with another which avoids parametric search by using more direct construction. The sparse DAG construction and application to all 3 metrics was moved to a short paper [12].

   The DAG constructed has Θ(n log^d-1 n) edges and vertices, hence the number of edges is linear in the number of vertices. This is true for the other orderings as well, which is why I have mentioned results for sparse graphs in the entries of algorithms for general orders. Without using this embedding, the natural representation of domination ordering for points in arbitrary position has m=Θ(n²), and so the bounds for general dense orderings would apply.

6. For L_∞ the approach used in [4,10] is based on using parametric search and solving the inverse problem. A direct approach for linear orders appears in [5,10], and for tree orders appears in [10], as do the results for dimensional orders with arbitrary placement. There is no known direct construction for grid orderings, d ≥ 2, which achieves the time bound obtained by the algorithm for general orders.

7. For L_∞ and unweighted data, the time reduces to Θ(m). Unweighted data does not improve the best time known for the other two metrics.

8. It is straightforward to show that if the graph is a star, then the time is linear for all of the metrics.

9. For the arbitrary DAG entries I included the crude bound, solely in terms of n, to make it easier to see its relationship to the entries above it, as well as the bound which relies on n and m, since that is the relevant version for the sparse graphs.

References:

1. Ahuja, RK and Olin, JB (2001), ``A fast scaling algorithm for minimizing separable convex functions subject to chain constraints'', Operations Research 49, pp. 784-789
2. Angelov, A, Harb, B, Kannan, S, and Wang, L-S (2006), ``Weighted isotonic regression under the L₁ norm'', Symp. Discrete Algorithms (SODA), pp. 783-791.
3. Ayer, M, Brunk, HD, Ewing, GM, Reid, WT, and Silverman, E (1955), ``An empirical distribution function for sampling with incomplete information'', Annals of Math. Stat. 5, pp. 641-647.
4. Kaufman, Y and Tamir, A (1993), ``Locating service centers with precedence constraints'', Discrete Applied Math. 47, pp. 251-261.
5. Lin, T-X, Kuo, C-C, Hsieh, Y-H, and Wang, B-F (2009), ``Efficient algorithms for the inverse sorting problem with bound constraints under the L_∞-norm and the Hamming distance'', J. Comp. and Sys. Sci., to appear.
6. Maxwell, WL and Muckstadt, JA (1985), ``Establishing consistent and realistic reorder intervals in production-distribution systems'', Operations Research 33, pp. 1316-1341.
7. Pardalos, PM and Xue, G (1999), ``Algorithms for a class of isotonic regression problems'', Algorithmica 23, pp. 211-222.
8. Spouge J, Wan H, and Wilbur WJ (2003), ``Least squares isotonic regression in two dimensions'', J. Optimization Theory and Apps. 117, pp. 585-605.
9. Stout, QF (2008), ``Unimodal regression via prefix isotonic regression'', Comp. Stat. and Data Anal. 53, pp. 289-297. A preliminary version appeared in ``Optimal algorithms for unimodal regression'', Computing and Stat. 32, 2000, and the supplemental portions were posted on the web in 2003.
10. Stout, QF (2009), ``Algorithms for L_∞ isotonic regression'', submitted. This is a revision of the original version (2008).
11. Stout, QF (2008), ``Isotonic median regression via scaling'', submitted.
12. Stout, QF (2009), ``An approach to computing multidimensional isotonic regressions'', submitted.

Copyright © 2009 Quentin F. Stout.