Remez exchange Algorithm for approximation

10 次查看（过去 30 天）
Kevser Cifci
Kevser Cifci 2021-10-13
评论： chun 2024-12-14
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!
  1 个评论
chun
chun 2024-12-14
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: chunwc@umich.edu

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Jan
Jan 2021-10-13
编辑：Jan 2021-10-13
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.

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Surrogate Optimization 的更多信息

标签

产品


版本

R2019a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by