Remez exchange Algorithm for approximation

Hello,

I tried to write Remez excahnge algorithm function but I don't know how to improve this code:

function [out] = remez(f,a,b) %% we want at the end of the function x, p et e

%% f(xi) - p(xi) = ((-1)^i)*e for i=0,...,n+1

%% we work on the interval [a,b]
%% a < x0 < ... < x(n) < x(n+1) < b : these are the points we work with in the function  


    x = a:0.05:b
    y = f(x);
    p = polyval(poly_alt_e(f,x)); %% poly_alt_e function return the alternating coefficients of the polynom p 
    %% and it gives also the value of e
    v = y - p;
    
    M = max(abs(v)); %% M is the maximum point of abs(f-p)
    
    %% we want to construct new (n+2) points including M and the (n+1) points we already have 
    %% we have to make sure that the signs of (f-p) alternate and I don't kwow how to make sure of that. 
    
    
    if a <= M < x(2) 
        if sign(v(M)) == sign(v(x(2)))
            x(x == x(2)) = [];
        else x = [M x];
             x(end)=[];
        end 
        
    %%% if M is between two points xj and xj+1 : 
    %%  we reject the point where the signs of (v(xj) or v(xj+1) and (v(M)) are equals
        for i = 1:length(x)
            M = max(abs(v));
            while sign(v(M))~=sign(v(i))
                if i < M < i+1 
                    x(x == i) = [];
                else x(x == i+1) = [];
                end
            end 
        end
        
     if v(end-1)< M <= b 
        if sign(v(M)) == sign(v(x(end-1)))
            x(x == x(end-1)) = [];
        else x = [x M];
             x(end)=[];
        end 
        
        for i = 1:length(x)
            M = max(abs(v));
            while sign(v(M))~=sign(v(i))
                if i < M < i+1 
                    x(x == i) = [];
                else x(x == i+1) = [];
                end
            end 
        end
    %% we have to stop when the value of e is not changing anymore
    %% or when x (here M) is not changing. 
end 

Can anyone tell me how to have the polynom p, the value e and x (here M in my function) ??

Thank you so much!

[chun]

I have done research on this function previously. If you have done the same, you will notice it is not only lacking in comments, but also, often fails to calculate the results. I later revised a similar remez function in python, and I published it on github. https://github.com/chunwangpro/Remez_Python

My code can work on most situations since I did a relaxation in find-root stage. I believe you can easily collect the data you want (the polynom p, the value e and x) in my code. I hope it helps you.

Welcome for questions: [e-mail]

Answer 1

[Jan]

if a <= M < x(2)
if v(end-1)< M <= b
if i < M < i+1

These commands will not do, what you want. The comparison is evaluated from left to right. The first expression a <= M replies TRUE or FALSE, which are treated as 1 or 0. This is the input for the 2nd comparison: 1 < x(2) or 0 < x(2) respectively. I assume you want:

if a <= M && M < x(2)
if v(end-1) < M && M <= b
if i < M && M < i+1

This looks strange:

    for i = 1:length(x)
         M = max(abs(v));
         while sign(v(M))~=sign(v(i))
             if i < M < i+1 
                 x(x == i) = [];
             else x(x == i+1) = [];
             end
         end 
    end    

M = max(abs(v)) is calculated in each iteration with the same input. Move the call before the loop and do it once only to save time. But

M = max(abs(v));
while sign(v(M)) ~= sign(v(i))

might be a bug already: M is the maximum absolute value and it is used as index of v. Perhaps you want the index of the maximum absolute value?

[~, M] = max(abs(v));

The missing comments make it hard to guess, what each line of code should do. Code without clear comments cannot be maintained or debugged efficiently.