Remez’s algorithm and best uniform approximation

Chebyshev showed that the polynomial of degree \(n\) closest to \(f\) in the maximum norm is characterized by an error curve that reaches its peak magnitude \(n + 2\) times with alternating sign. Evgeny Remez’s 1934 exchange algorithm finds it by repeatedly levelling the error on a set of reference points. This script computes the best degree-6 approximation to \(e^x \sin(3x)\) on \([-1, 1]\) and compares its error with Chebyshev interpolation of the same degree.

import matplotlib.pyplot as plt
import numpy as np

from mathematicskit.numerical_analysis import ChebyshevInterpolant, remez_minimax

The equioscillating error curve

f = lambda x: np.exp(x) * np.sin(3.0 * x)
n = 6
best = remez_minimax(f, n)
cheb = ChebyshevInterpolant(f, n + 1)  # degree n through n + 1 Chebyshev points

x = np.linspace(-1.0, 1.0, 1000)
err_best = f(x) - best.evaluate(x)
err_cheb = f(x) - cheb(x)

fig, ax = plt.subplots(figsize=(7, 4.5))
ax.plot(x, err_cheb, color="steelblue", label="Chebyshev interpolant")
ax.plot(x, err_best, color="firebrick", label="Remez minimax")
ax.scatter(best.reference, f(best.reference) - best.evaluate(best.reference), color="firebrick", zorder=3, label="reference points")
ax.axhline(best.max_error, color="gray", ls="--", lw=0.8)
ax.axhline(-best.max_error, color="gray", ls="--", lw=0.8)
ax.set_title(f"Degree {n}: the minimax error equioscillates at {n + 2} points")
ax.legend(fontsize=8)
fig.tight_layout()

print(f"Remez: {best.iterations} exchanges, max error {best.max_error:.3e}")
print(f"Chebyshev interpolation: max error {np.max(np.abs(err_cheb)):.3e}")

plt.show()

Remez: 4 exchanges, max error 3.860e-03
Chebyshev interpolation: max error 8.326e-03