A minimal naive implementation of the Remez exchange algorithm

remez.py

import numpy as np

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

import warnings

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


def gss(f, a, b, tol=1e-15):
    """ 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
    roots = fsolve(p, [E.min(), E.max()])
    search_endpoints = [E[0]] + list(roots) + [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), p(extrema)


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


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



def remez(p, x, n, w=lambda x: 1, max_iterations=100, tolerance=1e-10):
    """ 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
    """
    E = np.linspace(-1, 1, 100)
    # 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
    iteration = 0
    diff = -1
    A = np.ones((n+2,n+2))
    for j in range(len(A)):
        A[j, n+1] = (-1)**(n-j)
    while continuing(diff, iteration):
        # Construct a vandermonde matrix from the control nodes
        vandermonde = np.polynomial.chebyshev.chebvander(x, n)
        for i in range(n+2):
            wxi = w(x[i])
            for j in range(n+1):
                A[i,j] = vandermonde[i, j]*wxi
        # Evaluate the polynomial on the control nodes
        px = p(x)
        # Solve the system of linear equations
        v = solve(A, 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
        chebcoefs = v[:-1]
        coefficients = np.polynomial.chebyshev.cheb2poly(chebcoefs)[::-1]
        poly = np.poly1d(coefficients)
        # Calculate the extrema
        # We're doing a multi-point exchange, so swap all the control points
        x, y = extrema(lambda x: (p(x) - poly(x))*w(x), E)
        print(len(x))
        miny = np.abs(y).min()
        maxy = np.abs(y).max()
        diff = (maxy - miny)/maxy
        e = v[-1]
        print(diff)
        iteration += 1
    return poly


if __name__ == "__main__":
    z0 = -1/2 + 1j/2
    w = lambda x: np.abs(np.array(x) - z0)**2
    # TODOs: Weighting
    # Improve linsolve stability for higher degrees, investigate https://mpmath.org/ for arbitrary precision float types
    f = lambda x: 1/w(x)
    C = remez(f, np.array([-1, 0, 1]), w=w, n=1)
    #C3 = chebyshev(3)
    #C4 = chebyshev(4)
    #C = chebyshev(50, max_iterations=100)
    # 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])