# Problems Implementing the Remez Algorithm

So first off: *** This code is not being used in production software.
It is a personal project of mine, trying to understand approximation theory and advanced curve fitting.

In other words, I'm trying to understand how it works, not trying to get a currently existing solution.

So I have been trying to implement the Remez algorithm for polynomial approximation. I sort of/maybe/kind of have it working (not really).

My current solution generate ok polynomials, but a) the coefficients are not converging & b) while monitoring the coefficients at each stage, I've noticed that the x-values seem to slip past each other.

I'll give some examples to show what I mean.
My base function I'm trying to model is the square root function with a 4 degree polynomial on the domain [0.25, 1]

Round 1
X[0] = 0.25
X[1] = 0.4
X[2] = 0.55
X[3] = 0.7
X[4] = 0.85
X[5] = 1

Round 2
X[0] = 0.25
X[1] = 0.595076928220583
X[2] = 0.493988453622788
X[3] = 0.714333640596557
X[4] = 1.14135676154991
X[5] = 1

Round 3
X[0] = 0.25
X[1] = 0.638393337463021
X[2] = 0.63752199068821
X[3] = 0.538600997945798
X[4] = 1.07101841739164
X[5] = 1

Round 4
X[0] = 0.25
X[1] = 0.559423143598625
X[2] = 0.560673538304964
X[3] = 0.580378820375143
X[4] = 1.04454592077508
X[5] = 1


So here's a look at my actual code.

template <typename func_t>
type_t MiniMax(size_t Degree, const type_t& LowerLimit, const type_t& UpperLimit, unsigned char Iterations, func_t F0, func_t F1, func_t F2)
{
// (C) Jacob Wells 2015

if((Degree < 1) || (Iterations < 1))
{
return (type_t)NAN;
}

const type_t ONE(1), NEG1(-1), ZERO(0);
matrix_t<type_t> M(Degree + 2, Degree + 3);
vector<type_t> XVal;
type_t Delta, Sign, Pow, Err;
type_t D1, D2;
size_t I, J;

Delta = (UpperLimit - LowerLimit) / (Degree + 1);
XVal.resize(Degree + 2);
Coef.resize(Degree + 1);
Sign = ONE;

for(I = 0; I < (Degree + 2); I++)               // Generate our initial x-values
{
XVal[I] = (I * Delta) + LowerLimit;
}

do
{
Sign = NEG1;
for(I = 0; I < (Degree + 2); I++)
{
Pow = ONE;
M[I][Degree + 1] = Sign;            // Enters the alternating error sign
M[I][Degree + 2] = F0(XVal[I]);     // Enters the f(x) value
Sign *= NEG1;
for(J = 0; J <= Degree; J++)        // Evaluates the polynomial for each power
{
M[I][J] = Pow;
Pow *= XVal[I];
}
}

rref(Degree + 2, M);                    // Use Row Reduction Echelon Form to find the polynomial coefficients

Err = M[Degree + 1][Degree + 2];

for(I = 0; I <= Degree; I++)            // Copy the coefficients into the polynomial class' array
{
Coef[I] = M[I][Degree + 2];
}

if(Iterations > 1)
{
for(I = 1; I <= Degree; I++)        // Use Newton's method to find our new x-values
{
D1 = nth_deriv(XVal[I], 1) - F1(XVal[I]);
D2 = nth_deriv(XVal[I], 2) - F2(XVal[I]);

if((D1 != ZERO) && (D2 != ZERO))
{
XVal[I] -= (D1 / D2);
}
}
}

}while(--Iterations != 0);

return Err;
}


A quick little guide to some of my code:

F0, F1, & F2 are, respectively, the square root functions, it's first derivative, and it's second derivative.

nth_deriv calculate the Nth Derivative of the current polynomial.

rref reduces the Matrix to Row Reduction Echelon Form

matrix_t is a bare bones matrix class I came up with.

Now I have done a lot of testing on these functions, and I haven't found a single error with them, so I feel very confident that the problem is in the MiniMax Function.

EDIT:

My barebones matrix class

template <typename type_t>
class matrix_t
{
public:
~matrix_t() {}

matrix_t(size_t ColL, size_t RowL)
{
Arr.resize(ColL * RowL);
RowLen = RowL;
}

type_t* operator [] (size_t I)
{
return &Arr[I * RowLen];
}

const type_t* operator [] (size_t I) const
{
return &Arr[I * RowLen];
}

type_t At(size_t I)
{
return Arr[I];
}

private:
size_t RowLen;
vector<type_t> Arr;

matrix_t();
matrix_t(const matrix_t& m);
void operator = (const matrix_t& m);
};

• I've had occasion to recommend implementing the Remez algorithm in a spreadsheet to get a feel for how the exchange strategies work with polynomial interpolation. – hardmath Sep 3 '15 at 14:32
• OK, but the problem is I'm not sure if I'm implementing the algorithm correctly in the first place. – Mandalf The Beige Sep 3 '15 at 14:49
• Edited my original post. When I said barebones matrix, I meant nothing but what I needed. But to get back on track, from what I've read, my code is correct, but my program doesn't work, ergo I'm wrong. – Mandalf The Beige Sep 3 '15 at 14:57
• The "slipping past" is the evidence that there is a problem. What I mean by is that if you look at my starting x-values, you'll see that their all sorted, i.e. x[0] < x[1], x[1] < x[2], etc. In every example of the Remez algorithm I've seen even after several iterations in which the x-values are exchanged, they are still sorted. Mine don't stay sorted, and they don't seem to converge. Can I just ask, do you see anything obviously wrong with my code? I've been working on this for the last 6 months on my own, and haven't been able to figure anything out. – Mandalf The Beige Sep 3 '15 at 16:18
• I guess you are trying to do a replacement of more than one point at a time, but you should review the logic for doing so in Step 3 of the Wikipedia article. Using Newton's method as your code comments suggest can cause iterations that jump outside the prescribed intervals, and this seems to be your "sliding past" problem. – hardmath Sep 3 '15 at 19:33

Let's use the example (approximate the square root function on $[0.25,1.0]$ with a quartic polynomial) to step through your calculations. I suspect that the code is going to work with only modest changes.

Having chosen an initial set of six interpolation points:

$$x_k = 0.25, 0.4, 0.55, 0.7, 0.85, 1.0\;\; (k=0,\ldots,5)$$

we proceed to interpolate the function values at these points with a linear combination of a quartic polynomial and values $\pm 1$ that alternate at these interpolation points:

$$p(x_k) + (-1) e = c_4 x_k^4 + c_3 x_k^3 + c_2 x_k^2 + c_1 x + c_0 + (-1)^k e = \sqrt{x_k}$$

This interpolation problem amounts to solving a linear system:

$$\begin{pmatrix} 1 & x_0 & x_0^2 & x_0^3 & x_0^4 & +1 \\ 1 & x_1 & x_1^2 & x_1^3 & x_1^4 & -1 \\ 1 & x_2 & x_2^2 & x_2^3 & x_2^4 & +1 \\ 1 & x_3 & x_3^2 & x_3^3 & x_3^4 & -1 \\ 1 & x_4 & x_4^2 & x_4^3 & x_4^4 & +1 \\ 1 & x_5 & x_5^2 & x_5^3 & x_5^4 & -1 \end{pmatrix} \begin{pmatrix} c_0 \\ c_1 \\ c_2 \\ c_3 \\ c_4 \\ e \end{pmatrix} = \begin{pmatrix} \sqrt{x_0} \\ \sqrt{x_1} \\ \sqrt{x_2} \\ \sqrt{x_3} \\ \sqrt{x_4} \\ \sqrt{x_5} \end{pmatrix}$$

Of course we will resolve this problem each time we change the interpolation points.

Let's look at the solution for your initial set of points:

C0 =  0.195191580257
C1 =  1.506696296725
C2 = -1.365826192268
C3 =  0.938340584388
C4 = -0.274493119408
E = -9.0850E-5


You would solve for these using your "bare bones matrix class". I computed them using a spreadsheet, putting the $x_k$'s in one column A5:A10, the function values $\sqrt{x_k}$ in a parallel column C5:C10, then the matrix of coefficients above in a block E5:J10, and finally a matrix formula for a column beyond that something like:

   { = MMULT(MINVERSE(E5:J10),C5:C10) }


I am then able to change the interpolation points and (since everything else is given by a spreadsheet formula) quickly observe how the approximation error $e$ changes.

Intuitively you might think that as we converge our interpolation points to the relative extrema of approximation error on the interval, that the error $e$ will decrease in absolute value. But the opposite will occur.

The oscillating error $e$ that we are able to interpolate at a non-optimal set of points is less than the best error of approximation on the interval as a whole. This is in effect why we go looking in the interval for relative extrema of the error function, the exact function minus the approximating quartic polynomial:

$$\text{error at } x = \sqrt{x} - p(x)$$

We want to identify points inside the interval where the approximation error is bigger (in absolute value) than what we achieved at the interpolation points, namely $\pm e$.

It is typical (but not inevitable) that the endpoints of the interval will remain among the interpolation points throughout the computation. Such is the case here, and you kept these two endpoints in your computation. One of the points you brought into the solution in the second iteration was $x=1.14135676154991$, but that lies outside the interval and should be ignored.

Let's take the three extrema you found inside the interval and decide what these should be swapped out with among the original points. The basic rule is to look at the original points that bracket the new extremum, and swap with the point whose error of approximation is the same sign as that of the new point.

Taking the first extremum you found, $x_*=0.595076928220583$, we notice that it lies between our original points $0.55$ and $0.7$. Therefore we will swap it for whichever of those has the same sign error value. Because $e\approx -9.085\times 10^{-5}$, the error at $x_0 = 0.25$ is negative and alternates in sign with successive interpolation points. Thus the error was again negative at $x_2 = 0.55$ and positive at $x_3 = 0.7$.

We compute the error of approximation at the new point $x_* = 0.595076928220583$, something that is easily added to our spreadsheet:

\begin{align*} \text{error at } x_* &= \sqrt{x_*} - p(x_*) \\ &\approx -3.0107\times 10^{-5} \end{align*}

Based solely on the criterion of matching signs of errors, we would replace old point $0.55$ with $0.595076928220583$, since both have a negative error value. But we see that something is wrong! This new point $x_*$ is clearly not an extremum of the error function; the error there is smaller in absolute value than at all the existing interpolation points.

Okay, let's back up and figure out what the actual extrema of the error values are. For this purpose a graphing calculator application is helpful. I had never used it before, but GeoGebra gets a fair amount of Community ad spots at Math.SE, and it was pretty easy to plot this:

Visually (because high accuracy is not vital at this stage) the six points labelled A,B,C,D,E,F above are the local extrema, and they become our new interpolation points:

$$x_k = 0.25, 0.31, 0.51, 0.71, 0.92, 1.0\;\; (k=0,\ldots,5)$$

Once these replaced the original points in the spreadsheet, the new value of "equi-oscillating" error value is found, $e \approx -1.6737\times 10^{-4}$. Note that this is, as theory predicts, larger in absolute value than the "level error" of interpolation for the first set of points:

C0 =  0.192730711777
C1 =  1.527365848218
C2 = -1.421661988977
C3 =  0.999188048845
C4 = -0.297789990497
E = -1.6737E-04


Lather, rinse, repeat.

• I've been staring at this for the last week trying to figure out what to say, to no avail. So I'll just be straight with you: I'm sorry, but this doesn't help me. you summarize the Remez algorithm, and I know that will help someone else looking into it, but the one part I needed help with, the Multi-Point Exchange, the solution you present is to use a graph. I appreciate you trying to help me, but I have some other contacts I'm going to try. Again, thank you for trying. – Mandalf The Beige Sep 16 '15 at 22:02
• I took your example and showed you got the revised interpolation points (extrema of the error function) wrong. Since you did not include the coefficients of your quartic polynomial, I have no way of knowing whether your interpolation of the original (equally spaced) points was wrong, or your search for extrema failed. So I included all the details of the first iteration so that you can make this comparison. I would be curious to know. – hardmath Sep 17 '15 at 0:32