An implementation of the Remez exchange algorithm

remez.py

import numpy as np

from scipy.linalg import norm, solve
from itertools import pairwise


golden_ratio = (np.sqrt(5) + 1) / 2


def gss(f, a, b, tol=1e-10):
    """ Golden section search implementation
        Python code code here is from https://en.wikipedia.org/wiki/Golden-section_search#Iterative_algorithm
    """
    c = b - (b - a) / golden_ratio
    d = a + (b - a) / golden_ratio
    while abs(b - a) > tol:
        if f(c) < f(d):
            b = d
        else:
            a = c
        c = b - (b - a) / golden_ratio
        d = a + (b - a) / golden_ratio
    return (b + a) / 2


def extrema(p, E):
    """ Compute the extrema of the polynomial p over the discretized domain E
    """
    # We break E into the set of intervals A = [X0, r0], [r1, r2], ...., [r_{n-1}, X1]
    # where X0, X1 are the left and right endpoints of E, and the r_i correspond to
    # roots of p
    # We'll search for an extremum in each interval using golden section search
    search_endpoints = [E[0]] + list(p.r) + [E[-1]]
    intervals = pairwise(search_endpoints)
    extrema = []
    for a, b in intervals:
        # As written, our golden section search algorithm searches for a minimum value
        # so we search for the minimum of the function -|p(x)|
        extrema.append(gss(lambda x: -np.abs(p(x)), a, b))
    return np.array(extrema)


def monomial(n):
    """ Given n, returns the monomial x^n
    """
    return np.poly1d([1] + [0 for _ in range(n)])


def chebyshev(n):
    """ Returns the normalized Chebyshev polynomial of degree n for the interval [-1, 1]
    """
    m = monomial(n)
    E = np.linspace(-1, 1, 1000000)
    critical_points = np.sort([-np.cos(k*np.pi/n) for k in range(n)])
    T = remez(m, critical_points, E)
    TT = norm(T(E), np.inf)
    return T/TT


def remez(p, x, E=np.linspace(-1, 1, 100), max_iterations=20, tolerance=1e-08):
    """ Computes p*, the best approximation to p, using the Remez exchange algorithm
        and returns the error function p - p*

        x is the initial array of control nodes
        E is a discretized domain over which we'll search for extrema
    """
    # We use this function to determine if we've hit a stopping condition
    # The stopping conditions are that we've hit max_iterations, or the error
    # is at a tolerable threshold (configured by the tolerance parameter)
    continuing = lambda e, i: not np.isclose(e, 0, atol=tolerance) and i < max_iterations
    # On our first iteration, we evaluate the polynomial we're approximating
    poly = p
    i = 0
    e = -1
    while continuing(e, i):
        # Construct a vandermonde matrix from the control nodes
        vandermonde = np.vander(x)
        # Roll the vandermonde matrix
        # If v is [-1, 0, 1], then its Vandermonde matrix (decreasing powers)
        # will be
        # [1, -1, 1]
        # [0, 0, 1]
        # [1, 1, 1]
        # Rolling shifts everything columnwise 1 unit to the left, giving
        # [-1, 1, 1]
        # [0, 1, 0]
        # [1, 1, 1]
        rolled = np.roll(vandermonde, -1, axis=1)
        # Update the elements of the last column to contain the error
        # terms (-1)^j
        for j, row in enumerate(rolled):
            row[-1] = (-1)**j
        # Evaluate the polynomial on the control nodes
        px = (poly)(x)
        # Solve the system of linear equations
        v = solve(rolled, px)
        # The last element of v is the error, so we use only the first n+1 elements
        # as the coefficients of our approximating polynomial
        coefficients = v[:-1]
        poly = np.poly1d(coefficients)
        # Calculate the extrema
        # We're doing a multi-point exchange, so swap all the control points
        x = extrema(p - poly, E)
        # The error term lives in the last element of v, so we use that and i to determine
        # whether we're ready to stop
        e = v[-1]
        i += 1
    return p - poly


if __name__ == "__main__":
    # The following return the Chebyshev polynomials for [-1, 1] using the Remez exchange algorithm implementation
    # above as expected
    C2 = chebyshev(2) # returns poly1d([ 2.,  0., -1.])
    C3 = chebyshev(3) # returns poly1d([ 3.99995843e+00,  1.03915152e-05, -2.99996363e+00, -5.19575763e-06])
    C4 = chebyshev(4) # returns poly1d([ 8.0000000e+00, -8.8817842e-16, -8.0000000e+00,  0.0000000e+00, 1.0000000e+00])
    C10 = chebyshev(10) # returns poly1d([ 5.11999708e+02, -5.85122951e-05, -1.27999931e+03,  1.18160954e-04, 1.11999944e+03, -7.66385893e-05, -3.99999814e+02,  1.59482896e-05, 4.99999778e+01,  0.00000000e+00, -9.99999430e-01])
    # This does not return the Chebyshev polynomial due to an ill-conditioned matrix:
    # LinAlgWarning: Ill-conditioned matrix (rcond=1.0329e-18): result may not be accurate.
    C5 = chebyshev(5) # returns poly1d([ 6.16848674,  0.67739591, -6.38875921, -0.91652085,  1.06249059, 0.08134306])