Durand-Kerner Roots - Python

fast polynomial roots.ipynb

notebook.ipynb

Performance improvement

Performing

polynom = polynom / polynom[0]

i.e., by ensuing the coefficient of the highest power of x to be 1 significantly improves the numerical stability.

For example, the above algorithm fails for:

# 1 - fail example:
polynom = np.array([4, 0, -1])
benchmark(polynom)
# 2 - fail example
polynom = np.array([2, 13, -7])
benchmark(polynom)

but the problem gets resolved when executing as such:

# 1 - fail example - resolved:
polynom = np.array([4, 0, -1])
polynom = polynom / polynom[0]
benchmark(polynom)
# 2 - fail example - resolved
polynom = np.array([2, 13, -7])
polynom = polynom / polynom[0]
benchmark(polynom)

Thank you so much for doing this! I was a bit lost after reading the Wikipedia article, but this write-up made it so much clearer. I would also like to propose an alternative implementation that, even if not 100% equivalent, makes the computation a lot more SIMD-friendly (the actual SIMD implementation is missing, though):

import numpy as np


def durand_kerner(
        p: np.ndarray,
        max_iter: int = 1000,
        tol: float = 1e-8
) -> np.ndarray:
    # Polynomial scaling (roots remain unchanged)
    p = p / p[0]

    # Deterministic init using the complex unit circle
    deg = len(p) - 1
    roots = np.exp(np.arange(deg) * 2j * np.pi / deg)

    f0, cnt = np.polyval(p, roots), 0
    avg_delta = tol + 1  # do-while enforcing
    while avg_delta > tol and cnt < max_iter:
        cnt += 1

        # Simultaneous distance calculation using advanced indexing
        dens = np.prod(roots[:, None] - roots[None, :] + np.eye(deg), axis=1)
        roots = roots - f0 / dens

        f1 = np.polyval(p, roots)
        # 2-norm of the change in polyval
        avg_delta = np.sqrt(np.sum((f1 - f0)**2) / deg)
        f0 = f1

    return roots