aberth-hybrid

Hybrid contour-integral + Ehrlich-Aberth pipeline for finding all zeros of an analytic function (or all eigenvalues of a matrix pencil) in a region of the complex plane. Handles polynomials, transcendentals, bounded regions, unbounded strips, and non-normal matrix nonlinear eigenvalue problems (NEPs).

Reference implementation for the preprint "A Hybrid Contour-Integral / Ehrlich-Aberth Root-Finder for Analytic and Transcendental Functions on Unbounded Regions" (arXiv: tbd).

Highlights

• 2.2× faster and 58× tighter median residual than skzeros on the Bowhay-Nakatsukasa-Zaid combustion DDE benchmark (Python vs Python, same 2024 Apple M-series laptop). Julia port is 2.6× faster than skzeros.
• 17× to 415× faster than cxroots on scalar analytic root-finding in a disk.
• The only tool that solves the Hayes DDE unbounded strip in a single call (all 86 Lambert-$W$ zeros in $|\mathrm{Im}(z)| < 265$, at machine precision).
• 10⁷× better separation resolution than the standard Hermitian Rayleigh polish on close-eigenvalue non-normal matrix pencils, from the biorthogonal block Aberth iteration.
• Passes the standard NLEVP benchmark suite at machine-precision backward error (12 canonical problems from polynomial to transcendental).

Install

pip install aberth-hybrid          # once published to PyPI
# or, faster:
uv pip install aberth-hybrid

From source (for contributors):

git clone https://github.com/peppe-tringali/aberth-hybrid
cd aberth-hybrid
uv sync --extra dev
uv run pytest                      # 27 unit tests

Three canonical examples

Example 1: scalar analytic function in a bounded disk

import numpy as np
from aberth_hybrid import hybrid_solve

# Zeros of z + exp(-z) inside the disk |z| < 8
def f(z):  return z + np.exp(-z)
def fp(z): return 1.0 - np.exp(-z)

result = hybrid_solve(f, fp, center=0+0j, radius=8.0, tol=1e-13)
print(f"{result.n_zeros} zeros recovered")
for z in sorted(result.roots, key=lambda z: z.imag):
    print(f"  {z:+.6f}   |f(z)| = {abs(f(z)):.2e}")

Expected output:

4 zeros recovered
  -2.062278-7.588631j   |f(z)| = 9.93e-16
  -0.318132-1.337236j   |f(z)| = 2.29e-16
  -0.318132+1.337236j   |f(z)| = 2.29e-16
  -2.062278+7.588631j   |f(z)| = 9.93e-16

Example 2: unbounded strip (Hayes DDE)

import numpy as np
from aberth_hybrid import Window, solve_on_windows

def f(z):  return z + np.exp(-z)      # Hayes DDE characteristic function
def fp(z): return 1.0 - np.exp(-z)

# Cover a vertical strip along the imaginary axis with overlapping disks
windows = [Window(center=complex(-2.5, k * 8.8), radius=5.5)
           for k in range(-30, 31)]
result = solve_on_windows(f, fp, windows, tol=1e-13)

print(f"{len(result.roots)} zeros in the strip |Im(z)| < 265")

Expected: 86 zeros, each at machine precision, all Lambert-$W$ branches.

Example 3: matrix nonlinear eigenvalue problem

import numpy as np
from aberth_hybrid.beyn import beyn_extract
from aberth_hybrid.block_aberth import block_aberth_from_beyn

# Simple 3x3 matrix DDE: T(lam) = -lam*I + A + B*exp(-lam)
A = np.diag([-2.0, -3.0, -5.0]) + 0j
B = np.diag([1.5, 2.0, 1.2]) + 0j
def T(z):
    return -complex(z) * np.eye(3, dtype=complex) + A + B * np.exp(-complex(z))
def T_prime(z):
    return -np.eye(3, dtype=complex) - B * np.exp(-complex(z))

# Beyn extraction on a disk
beyn = beyn_extract(T, center=0+0j, radius=3.0, n_quad=1024,
                    n_probes=10, n_moments=10,
                    return_left_eigenvectors=True)

# Refine with prewarm + biorthogonal block Aberth
result = block_aberth_from_beyn(T, T_prime, beyn, tol=1e-14,
                                 prewarm_steps=3)
print(f"{len(result.eigenvalues)} eigenvalues:")
for lam in sorted(result.eigenvalues, key=lambda z: (z.real, z.imag)):
    print(f"  {lam:+.6f}")

For a walkthrough of the theory (7-layer pipeline, prewarm theorem, certificate), see the companion blog series. For a fuller worked notebook, open docs/getting_started.ipynb.

Public API

Top-level (from aberth_hybrid import ...):

Matrix NEPs (from aberth_hybrid.beyn import ... and from aberth_hybrid.block_aberth import ...):

Parameter continuation (from aberth_hybrid.continuation import ...):

Julia port

A performance-focused Julia translation lives in julia_port/. It matches the Python implementation's algorithm and runs the combustion DDE benchmark in 9.8 ms (median of 5 runs), 17% faster than the Python reference. To run:

cd julia_port
julia --project=. test/runtests.jl              # 12 unit tests
julia --project=. benchmarks/combustion_dde.jl  # combustion benchmark

Benchmarks

Reproducible benchmarks under benchmarks/:

uv run python -m benchmarks.nlevp.run_all                            # NLEVP suite
uv run python -m benchmarks.head_to_head.scalar_vs_cxroots           # vs cxroots
uv run python -m benchmarks.head_to_head.matrix_vs_neppack           # vs NEP-PACK
uv run python -m benchmarks.head_to_head.combustion_vs_skzeros       # vs skzeros
uv run python -m benchmarks.head_to_head.vs_big_math_libraries       # vs numpy/scipy
bash benchmarks/run_full_suite.sh                                    # everything

Or in Julia:

julia --project=julia_port julia_port/benchmarks/combustion_dde.jl

Tests

uv run pytest -q          # 27 unit tests, all pass
uv run ruff check .       # lint clean

Citation

If you use aberth-hybrid in your work, please cite:

@misc{TringaliAberthHybrid2026,
  author       = {Tringali, Giuseppe},
  title        = {A Hybrid Contour-Integral / {Ehrlich--Aberth} Root-Finder
                  for Analytic and Transcendental Functions on Unbounded
                  Regions},
  year         = {2026},
  eprint       = {tbd},
  archivePrefix= {arXiv},
  primaryClass = {math.NA}
}

License

MIT.