l1fit

Fast constrained L1 curve fitting for Python and NumPy. l1fit solves

\min_{x \in \mathbb{R}^n} \; \lVert A x - b \rVert_1
\quad \text{s.t.} \quad C x = d, \;\; E x \le f

with A of size k x n, optional equality constraints C x = d (l rows) and optional inequality constraints E x <= f (m rows), and optional sign restrictions on x and on the residuals. The L1 norm makes the fit robust against outliers: where a least-squares line is pulled towards them, the L1 line ignores them.

The solver is the modified simplex method of Barrodale and Roberts (ACM TOMS algorithm 552, CL1), written in modern Fortran. It works on a single dense tableau, so it needs no LP modelling layer and is much faster than general-purpose solvers on this problem class (see Benchmark). The Python package is a thin, typed nanobind binding that releases the GIL.

Installation

l1fit is available on PyPI:

uv add l1fit        # or: pip install l1fit

There is one wheel per platform for Linux (x86_64, aarch64), macOS (arm64, x86_64) and Windows (x86_64); each covers Python 3.12 to 3.15 (stable ABI). They need nothing but NumPy: the Fortran code uses no runtime library, so there is no libgfortran to bundle. The wheels are built by .github/workflows/wheels.yml.

To build from source (other platforms, or for development) you need a Fortran compiler (gfortran; on Windows the MinGW-w64 one, next to MSVC for the C++), a C++17 compiler and CMake. uv builds the extension with scikit-build-core:

uv sync     # development environment (builds the extension in editable mode)
uv build    # wheel and sdist in dist/

Usage

import numpy as np
from l1fit import Status, solve_l1

# Fit a line through 7 points, one of them an outlier of size 10: the L1 fit ignores it.
t = np.arange(7.0)
b = 1.0 + 2.0 * t
b[3] += 10.0
A = np.column_stack([np.ones(7), t])

result = solve_l1(A, b)
assert result.status is Status.OPTIMAL
print(result.x)  # [1. 2.]
print(result.objective)  # 10.0 (up to rounding) = the sum of the absolute residuals

# Optional equality constraints C x = d, inequality constraints E x <= f and sign restrictions:
result = solve_l1(
    A,
    b,
    C=np.array([[1.0, 1.0]]),
    d=np.array([4.0]),  # x0 + x1 = 4
    E=np.array([[0.0, 1.0]]),
    f=np.array([1.5]),  # x1 <= 1.5
    x_sign=np.array([0, 1]),  # x1 >= 0
)

solve_l1(A, b, C, d, E, f, *, x_sign, residual_sign, tol, max_iter) takes float64 numpy arrays (other dtypes raise TypeError) and never modifies them. (The inequality matrix is called E, as in the paper.) It returns a frozen L1Result with x, residual (b - A x), equality_residual, inequality_slack, objective, status and iterations. Infeasibility, rounding errors and reaching max_iter are reported through result.status (Status.INFEASIBLE, ROUNDING_ERRORS, MAX_ITERATIONS), not as exceptions; invalid arguments raise TypeError or ValueError. NaN and infinity are rejected. The binding releases the GIL, and the solver is reentrant, so threads can solve independent problems in parallel.

The solver works in double precision. The default tolerance is 1e-10 and the default iteration limit 10 * (k + l + m).

Benchmark against CVXPY

CVXPY is the usual way to solve this problem in Python without a dedicated routine. The benchmark (benchmarks/bench_cvxpy.py, instances from tests/instances) uses the natural formulation:

x = cp.Variable(n)
problem = cp.Problem(cp.Minimize(cp.norm(b - A @ x, 1)), [C @ x == d, E @ x <= f])
problem.solve(solver=...)

CVXPY rewrites it as a linear program with one epigraph variable per row (|b_i - A_i x| <= t_i: n + k variables and 2k + m inequality rows) and passes it to a solver. Three were used:

• HiGHS dual simplex (solver="HIGHS", option solver="simplex") - the same algorithm family as l1fit, on the larger LP;
• HiGHS interior point (solver="HIGHS", option solver="ipm");
• Clarabel (solver="CLARABEL"), an interior-point solver for conic problems and CVXPY's default for LPs.

Wall-clock time of the whole call, minimum of 20 runs for l1fit and of 3 for CVXPY (a new cp.Problem for every run, so the times include CVXPY's canonicalization). The factor is the time relative to l1fit:

Apple M1 Pro, Python 3.12.13, NumPy 2.5.3, CVXPY 1.9.3 (HiGHS 1.15.1 through highspy, Clarabel 0.11.1), default solver tolerances, l1fit with its defaults (tol=1e-10). Raw numbers, versions and per-solver details are in benchmarks/results/cvxpy.json; reproduce with

uv run --group bench python benchmarks/bench_cvxpy.py --json benchmarks/results/cvxpy.json

l1fit is 9 to 127 times faster. Even when only the time that the solver itself reports is compared (leaving out CVXPY's modelling overhead), it stays 7 to 124 times faster: the dedicated simplex works on a k x n tableau, while the general LP is much larger. The optimal objectives agree to 1e-12 (HiGHS) and about 1e-9 (Clarabel, whose interior-point tolerance is looser). The comparison is between a specialised and a general-purpose solver on dense random problems (five with a planted feasible point and Student-t noise, cvxpy_script_1000x250 with pure Gaussian data); the factors will differ on other data, and simplex methods slow down on very large or degenerate problems. CVXPY, of course, solves far more than this one problem class.

Background: a modernized CALGO 552

The algorithm is not new: this is a refactored and modernized version of CALGO552.f, the Fortran 77 routine CL1 published in 1980 by I. Barrodale and F. D. K. Roberts (Algorithm 552: Solution of the constrained l1 linear approximation problem, ACM Transactions on Mathematical Software 6(2), 231-235). The original is a single subroutine with 18 arguments, about 60 numeric labels, GOTOs and single precision. The copy this project started from also had a few garbled lines that did not compile; they were repaired first. What was done:

• Source form and types: fixed form to free-form Fortran 2018 modules (implicit none (type, external), private by default); single precision to double precision (one kind parameter wp, the accumulator dp).
• Control flow: every GOTO, arithmetic IF and label removed. The single routine was split step by step into small procedures (pivot_tableau, select_entering_column, select_leaving_row, decide_at_optimum, set_up_phase1_costs, ...), each verified before the next one.
• Interface: assumed-shape contiguous arrays, so four dimension arguments (KLMD, KLM2D, NKLMD, N2D) and the three workspace arrays (CU, IU, S) are gone: the workspace is allocated inside, and 11 arguments are left. Arguments are validated, and failures are reported as status codes (0 optimal, 1 infeasible, 2 rounding errors, 3 iteration limit, 4 invalid input, 5 allocation failed) instead of undefined behaviour. A bind(C) entry point (l1_cl1) serves the Python extension.
• Behaviour kept: the algorithm itself was not changed. Compiled in single precision, the refactored code reproduced the original bit for bit (status, iterations, solution, residuals) on about 2600 test problems, and it runs at the same speed. Double precision is the only numerical difference: it makes the solver reliable on problems where the single-precision original stops with rounding errors (for example a 1000 x 250 problem with 200 inequality rows).
• Packaging: the Python extension, wheels for all major platforms, and a test suite that checks every answer against an independent HiGHS reference.

The history of the repository contains the original and the intermediate steps. The algorithm and the original Fortran are the work of Barrodale and Roberts and are credited as such.

Development

uv sync                                  # rebuilds when pyproject.toml, CMakeLists.txt or src/ change
uv run pytest tests                      # tests (about 6 s), see tests/README.md
uv run ruff check && uv run ruff format --check
uv run pyrefly check
uv run python -m nanobind.stubgen -m l1fit._core -o python/l1fit/_core.pyi   # after changing src/_core.cpp

Ruff and pyrefly check the Python package only (python/); Fortitude lints the Fortran (src/); see pyproject.toml and fortitude.toml.

License

BSD 3-Clause, see LICENSE.