An isotonic regression (https://en.wikipedia.org/wiki/Isotonic_regression) computes a fit (vector y_hat) which minimizes a specified loss to an observations vector y subject to isotonicity (non-decreasing) constraints.

Isotonic regression usually refers to the sum of squares loss (L2 norm), but can be generalized to other losses as well (sometimes referred to as generalized isotonic regression). For example, sum of absolute values (L1 norm). In additions, when a weight vector is specified the weighted loss is minimized. The total order indicates the constraints form a line graph, while partial order takes a generad Direct Acyclic Graph (DAG) as input. Sklearn has existing implementation for total order L2-isotonic regression.

Are there any open source python implementations for total order L1 isotonic regression?

Formnally, y_hat = argmin_z sum w_i |y_i - z_i|, subject to z_i <= z_j for i < j.


1 Answer 1

Copyright 2020 Google LLC.
SPDX-License-Identifier: Apache-2.0

import heapq
import numpy as np

def isotonic_regression_l1_total_order(y, w):
  """Finds a non-decreasing fit for the specified `y` under L1 norm.

  The O(n log n) algorithm is described in:
  "Isotonic Regression by Dynamic Programming", Gunter Rote, SOSA@SODA 2019.

    y: The values to be fitted, 1d-numpy array.
    w: The loss weights vector, 1d-numpy array.

    An isotonic fit for the specified `y` which minimizies the weighted
    L1 norm of the fit's residual.
  h = []  # max heap of values
  p = np.zeros_like(y)  # breaking position
  for i in range(y.size):
    a_i = y[i]
    w_i = w[i]
    heapq.heappush(h, (-a_i, 2 * w_i))
    s = -w_i
    b_position, b_value = h[0]
    while s + b_value <= 0:
      s += b_value
      b_position, b_value = h[0]
    b_value += s
    h[0] = (b_position, b_value)
    p[i] = -b_position
  z = np.flip(np.minimum.accumulate(np.flip(p)))  # right_to_left_cumulative_min
  return z

Example use:

isotonic_regression_l1_total_order(np.array([1, 3, 5]), np.ones((3,)))
isotonic_regression_l1_total_order(np.array([5, 3, 1]), np.ones((3,)))

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.