Created
April 25, 2023 21:27
-
-
Save polynomialherder/85219bbb5163a55eca15ddd735828b0e to your computer and use it in GitHub Desktop.
An implementation of the Remez exchange algorithm
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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]) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment