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!

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.