Problems with Exchange Procedure in 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 about 6 months ago I started implementing the Remez algorithm for polynomial approximation. It sorts of works, but I can't get the X-Values to converge during the procedure. I've been toying around with it for the last few months, but have finally given up trying to solve it on my own.

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.e. the x-values start off sorted, but after several iterations of the exchange process, they become unsorted.

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

And 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.