Problems Implementing the Remez Algorithm

Question

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
    // This code is licensed under the BSD 3-Clause License
    // http://opensource.org/licenses/BSD-3-Clause

    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);
};

Answer 1

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.