# Are these any existing implementation of L1 isotonic regression in python?

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.

"""
"""

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.

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

Returns:
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
heapq.heappop(h)
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,)))