Articulo de referencia

Algoritmo de Remez

El algoritmo de Remez o algoritmo de intercambio de Remez , publicado por Evgeny Yakovlevich Remez en 1934, es un algoritmo iterativo que se utiliza para encontrar aproximacione...

El algoritmo de Remez o algoritmo de intercambio de Remez , publicado por Evgeny Yakovlevich Remez en 1934, es un algoritmo iterativo que se utiliza para encontrar aproximaciones simples de funciones, específicamente, aproximaciones de funciones en un espacio de Chebyshev que son las mejores en el sentido de la norma uniforme L . [ 1 ] A veces se le denomina algoritmo de Remes o algoritmo de Reme . [ 2 ]

Un ejemplo típico de espacio de Chebyshev es el subespacio de polinomios de Chebyshev de orden n en el espacio de funciones continuas reales en un intervalo C [ a , b ]. El polinomio de mejor aproximación dentro de un subespacio dado se define como aquel que minimiza la máxima diferencia absoluta entre el polinomio y la función. En este caso, la forma de la solución se precisa mediante el teorema de equioscilación .

Procedimiento

El algoritmo de Remez comienza con la funciónF{\displaystyle f}ser aproximado y un conjuntoincógnita{\displaystyle X}denorte+2{\displaystyle n+2}puntos de muestraincógnita1,incógnita2,...,incógnitanorte+2{\displaystyle x_{1},x_{2},...,x_{n+2}}En el intervalo de aproximación, generalmente los extremos del polinomio de Chebyshev se mapean linealmente al intervalo. Los pasos son:

  • Resuelve el sistema de ecuaciones lineales.
b0+b1incógnitai+...+bnorteincógnitainorte+(1)imi=F(incógnitai){\displaystyle b_{0}+b_{1}x_{i}+...+b_{n}x_{i}^{n}+(-1)^{i}E=f(x_{i})}(dóndei=1,2,...norte+2{\displaystyle i=1,2,...n+2}),
para lo desconocidob0,b1...bnorte{\displaystyle b_{0},b_{1}...b_{n}}y E.
  • Utilice elbi{\displaystyle b_{i}}como coeficientes para formar un polinomioPAGnorte{\displaystyle P_{n}}.
  • Encuentra el conjuntoMETRO{\displaystyle M}de puntos de error máximo local|PAGnorte(incógnita)F(incógnita)|{\displaystyle |P_{n}(x)-f(x)|}.
  • Si los errores en cadametroMETRO{\displaystyle m\in M}son de igual magnitud y alternan en signo, entoncesPAGnorte{\displaystyle P_{n}}es el polinomio de aproximación minimax. Si no, reemplácelo.incógnita{\displaystyle X}conMETRO{\displaystyle M}y repita los pasos anteriores.

El resultado se denomina polinomio de mejor aproximación o algoritmo de aproximación minimax .

W. Fraser ofrece una revisión de los aspectos técnicos de la implementación del algoritmo de Remez. [ 3 ]

Elección de inicialización

Los nodos de Chebyshev son una opción común para la aproximación inicial debido a su papel en la teoría de la interpolación polinómica . Para la inicialización del problema de optimización para la función f mediante el interpolante de Lagrange L n ( f ), se puede demostrar que esta aproximación inicial está acotada por

FLnorte(F)(1+Lnorte)infpagPAGnorteFpag{\displaystyle \lVert f-L_{n}(f)\rVert _{\infty }\leq (1+\lVert L_{n}\rVert _{\infty })\inf _{p\in P_{n}}\lVert fp\rVert }

siendo la norma o constante de Lebesgue del operador de interpolación de Lagrange L n de los nodos ( t 1 , ..., t n  +  1 )

Lnorte=Λ¯norte(T)=máximo1incógnita1λnorte(T;incógnita),{\displaystyle \lVert L_{n}\rVert _{\infty }={\overline {\Lambda }}_{n}(T)=\max _{-1\leq x\leq 1}\lambda _{n}(T;x),}

siendo T los ceros de los polinomios de Chebyshev y las funciones de Lebesgue las

λnorte(T;incógnita)=j=1norte+1|lj(incógnita)|,lj(incógnita)=iji=1norte+1(incógnitati)(tjti).{\displaystyle \lambda _{n}(T;x)=\sum _{j=1}^{n+1}\left|l_{j}(x)\right|,\quad l_{j}(x)=\prod _{\stackrel {i=1}{i\neq j}}^{n+1}{\frac {(x-t_{i})}{(t_{j}-t_{i})}}.}

Theodore A. Kilgore, [ 4 ] Carl de Boor y Allan Pinkus [ 5 ] demostraron que existe un único t i para cada L n , aunque no se conoce explícitamente para polinomios (ordinarios). De manera similar,Λ_norte(T)=min1incógnita1λnorte(T;incógnita){\displaystyle {\underline {\Lambda }}_{n}(T)=\min _{-1\leq x\leq 1}\lambda _{n}(T;x)}y la optimalidad de una elección de nodos se puede expresar comoΛ¯norteΛ_norte0.{\displaystyle {\overline {\Lambda }}_{n}-{\underline {\Lambda }}_{n}\geq 0.}

Para los nodos de Chebyshev, que proporcionan una elección subóptima, pero analíticamente explícita, el comportamiento asintótico se conoce como [ 6 ].

Λ¯norte(T)=2πregistro(norte+1)+2π(γ+registro8π)+αnorte+1{\displaystyle {\overline {\Lambda }}_{n}(T)={\frac {2}{\pi }}\log(n+1)+{\frac {2}{\pi }}\left(\gamma +\log {\frac {8}{\pi }}\right)+\alpha _{n+1}}

(donde γ es la constante de Euler-Mascheroni ) con

0<αnorte<π72norte2{\displaystyle 0<\alpha _{n}<{\frac {\pi }{72n^{2}}}}paranorte1,{\displaystyle n\geq 1,}

y límite superior [ 7 ]

Λ¯norte(T)2πregistro(norte+1)+1{\displaystyle {\overline {\Lambda }}_{n}(T)\leq {\frac {2}{\pi }}\log(n+1)+1}

Lev Brutman [ 8 ] obtuvo la cota paranorte3{\displaystyle n\geq 3}, yT^{\displaystyle {\hat {T}}}siendo los ceros de los polinomios de Chebyshev expandidos:

Λ¯norte(T^)Λ_norte(T^)<Λ¯316cunaπ8+π641pecado2(3π/16)2π(γregistroπ)0,201.{\displaystyle {\overline {\Lambda }}_{n}({\hat {T}})-{\underline {\Lambda }}_{n}({\hat {T}})<{\overline {\Lambda }}_{3}-{\frac {1}{6}}\cot {\frac {\pi }{8}}+{\frac {\pi }{64}}{\frac {1}{\sin ^{2}(3\pi /16)}}-{\frac {2}{\pi }}(\gamma -\log \pi )\approx 0.201.}

Rüdiger Günttner [ 9 ] obtuvo a partir de una estimación más precisa paranorte40{\displaystyle n\geq 40}

Λ¯norte(T^)Λ_norte(T^)<0,0196.{\displaystyle {\overline {\Lambda }}_{n}({\hat {T}})-{\underline {\Lambda }}_{n}({\hat {T}})<0.0196.}

Discusión detallada

Esta sección proporciona más información sobre los pasos descritos anteriormente. En esta sección, el índice i va de 0 a n + 1.

Paso 1: Dadoincógnita0,incógnita1,...incógnitanorte+1{\displaystyle x_{0},x_{1},...x_{n+1}}, resolver el sistema lineal de n + 2 ecuaciones

b0+b1incógnitai+...+bnorteincógnitainorte+(1)imi=F(incógnitai){\displaystyle b_{0}+b_{1}x_{i}+...+b_{n}x_{i}^{n}+(-1)^{i}E=f(x_{i})}(dóndei=0,1,...norte+1{\displaystyle i=0,1,...n+1}),
para lo desconocidob0,b1,...bnorte{\displaystyle b_{0},b_{1},...b_{n}}y E.

Debe quedar claro que(1)imi{\displaystyle (-1)^{i}E}en esta ecuación solo tiene sentido si los nodosincógnita0,...,incógnitanorte+1{\displaystyle x_{0},...,x_{n+1}}están ordenados , ya sea estrictamente crecientes o estrictamente decrecientes. Entonces este sistema lineal tiene una solución única. (Como es bien sabido, no todos los sistemas lineales tienen una solución). Además, la solución se puede obtener con soloO(norte2){\displaystyle O(n^{2})}operaciones aritméticas, mientras que un solucionador estándar de la biblioteca tardaríaO(norte3){\displaystyle O(n^{3})}operaciones. He aquí la prueba sencilla:

Calcular el interpolante estándar de grado npag1(incógnita){\displaystyle p_{1}(x)}aF(incógnita){\displaystyle f(x)}en los primeros n + 1 nodos y también el interpolante estándar de grado npag2(incógnita){\displaystyle p_{2}(x)}a las ordenadas(1)i{\displaystyle (-1)^{i}}

pag1(incógnitai)=F(incógnitai),pag2(incógnitai)=(1)i,i=0,...,norte.{\displaystyle p_{1}(x_{i})=f(x_{i}),p_{2}(x_{i})=(-1)^{i},i=0,...,n.}

Para ello, utilice en cada caso la fórmula de interpolación de Newton con las diferencias divididas de orden0,...,norte{\displaystyle 0,...,n}yO(norte2){\displaystyle O(n^{2})}operaciones aritméticas.

El polinomiopag2(incógnita){\displaystyle p_{2}(x)}tiene su i -ésimo cero entreincógnitai1{\displaystyle x_{i-1}}yincógnitai, i=1,...,norte{\displaystyle x_{i},\ i=1,...,n}y por lo tanto no hay más ceros entreincógnitanorte{\displaystyle x_{n}}yincógnitanorte+1{\displaystyle x_{n+1}}:pag2(incógnitanorte){\displaystyle p_{2}(x_{n})}ypag2(incógnitanorte+1){\displaystyle p_{2}(x_{n+1})}tienen el mismo signo(1)norte{\displaystyle (-1)^{n}}.

La combinación lineal pag(incógnita):=pag1(incógnita)pag2(incógnita)mi{\displaystyle p(x):=p_{1}(x)-p_{2}(x)\!\cdot \!E}es también un polinomio de grado n y

pag(incógnitai)=pag1(incógnitai)pag2(incógnitai)mi = F(incógnitai)(1)imi,    i=0,,norte.{\displaystyle p(x_{i})=p_{1}(x_{i})-p_{2}(x_{i})\!\cdot \!E\ =\ f(x_{i})-(-1)^{i}E,\ \ \ \ i=0,\ldots ,n.}

Esto es lo mismo que la ecuación anterior parai=0,...,norte{\displaystyle i=0,...,n}y para cualquier elección de E. La misma ecuación para i = n + 1 es

pag(incógnitanorte+1) = pag1(incógnitanorte+1)pag2(incógnitanorte+1)mi = F(incógnitanorte+1)(1)norte+1mi{\displaystyle p(x_{n+1})\ =\ p_{1}(x_{n+1})-p_{2}(x_{n+1})\!\cdot \!E\ =\ f(x_{n+1})-(-1)^{n+1}E}y requiere un razonamiento especial: resuelto para la variable E , es la definición de E :
mi := pag1(incógnitanorte+1)F(incógnitanorte+1)pag2(incógnitanorte+1)+(1)norte.{\displaystyle E\ :=\ {\frac {p_{1}(x_{n+1})-f(x_{n+1})}{p_{2}(x_{n+1})+(-1)^{n}}}.}

Como se mencionó anteriormente, los dos términos en el denominador tienen el mismo signo: E y por lo tantopag(incógnita)b0+b1incógnita++bnorteincógnitanorte{\displaystyle p(x)\equiv b_{0}+b_{1}x+\ldots +b_{n}x^{n}}Siempre están bien definidos.

El error en los n + 2 nodos ordenados dados es positivo y negativo alternativamente porque

pag(incógnitai)F(incógnitai) = (1)imi,  i=0,...,norte+1.{\displaystyle p(x_{i})-f(x_{i})\ =\ -(-1)^{i}E,\ \ i=0,...,n\!+\!1.}

El teorema de equioscilación establece que bajo esta condición no existe ningún polinomio de grado n con error menor que E. De hecho, si existiera tal polinomio, llamémoslopag~(incógnita){\displaystyle {\tilde {p}}(x)}, entonces la diferencia pag(incógnita)pag~(incógnita)=(pag(incógnita)F(incógnita))(pag~(incógnita)F(incógnita)){\displaystyle p(x)-{\tilde {p}}(x)=(p(x)-f(x))-({\tilde {p}}(x)-f(x))}seguiría siendo positivo/negativo en los nodos n + 2incógnitai{\displaystyle x_{i}}y por lo tanto tienen al menos n + 1 ceros, lo cual es imposible para un polinomio de grado n . Por lo tanto, este E es una cota inferior para el error mínimo que se puede lograr con polinomios de grado n .

El paso 2 cambia la notación de b0+b1incógnita+...+bnorteincógnitanorte{\displaystyle b_{0}+b_{1}x+...+b_{n}x^{n}}apag(incógnita){\displaystyle p(x)}.

El paso 3 mejora los nodos de entrada.incógnita0,...,incógnitanorte+1{\displaystyle x_{0},...,x_{n+1}}y sus errores±mi{\displaystyle \pm E}como sigue.

En cada región P, el nodo actualincógnitai{\displaystyle x_{i}}se reemplaza con el maximizador localincógnita¯i{\displaystyle {\bar {x}}_{i}}y en cada región Nincógnitai{\displaystyle x_{i}}se reemplaza con el minimizador local. (Se esperaincógnita¯0{\displaystyle {\bar {x}}_{0}}en A , elincógnita¯i{\displaystyle {\bar {x}}_{i}}cercaincógnitai{\displaystyle x_{i}}, yincógnita¯norte+1{\displaystyle {\bar {x}}_{n+1}}en B .) No se requiere alta precisión aquí, la búsqueda lineal estándar con un par de ajustes cuadráticos debería ser suficiente. (Ver [ 10 ] )

Dejarzi:=pag(incógnita¯i)F(incógnita¯i){\displaystyle z_{i}:=p({\bar {x}}_{i})-f({\bar {x}}_{i})}Cada amplitud|zi|{\displaystyle |z_{i}|}es mayor o igual que E. El teorema de La Vallée Poussin y su demostración también se aplican az0,...,znorte+1{\displaystyle z_{0},...,z_{n+1}}conmin{|zi|}mi{\displaystyle \min\{|z_{i}|\}\geq E}como el nuevo límite inferior para el mejor error posible con polinomios de grado n .

Además,máximo{|zi|}{\displaystyle \max\{|z_{i}|\}}Resulta útil como límite superior obvio para ese mejor error posible.

Paso 4: Conmin{|zi|}{\displaystyle \min \,\{|z_{i}|\}}ymáximo{|zi|}{\displaystyle \max \,\{|z_{i}|\}}como límite inferior y superior para el mejor error de aproximación posible , se tiene un criterio de parada fiable: repetir los pasos hastamáximo{|zi|}min{|zi|}{\displaystyle \max\{|z_{i}|\}-\min\{|z_{i}|\}}es suficientemente pequeño o ya no disminuye. Estos límites indican el progreso.

Variantes

En la literatura se encuentran algunas modificaciones del algoritmo. [ 11 ] Estas incluyen:

  • Sustituir más de un punto de muestreo por las ubicaciones de las diferencias absolutas máximas cercanas.
  • Reemplazar todos los puntos de muestra en una sola iteración con las ubicaciones de todas las diferencias máximas, alternando el signo. [ 12 ]
  • Utilizar el error relativo para medir la diferencia entre la aproximación y la función, especialmente si la aproximación se va a utilizar para calcular la función en un ordenador que utiliza aritmética de punto flotante ;
  • Incluyendo restricciones de punto de error cero. [ 12 ]
  • La variante de Fraser-Hart, utilizada para determinar la mejor aproximación racional de Chebyshev. [ 13 ]

Véase también

Referencias

  1. ^ Remez, E. Ya. (1934). "Sobre la determinación de polinômes d'approximation de gré donnée". Com. Soc. Matemáticas. Jarkov . 10 : 41. (1934). "Sobre un procedimiento convergente de aproximaciones sucesivas para determinar los polinomios de aproximación" . compt. Desgarrar. Acad. Ciencia. (en francés). 198 : 2063-5 . (1934). "Sobre el cálculo efectivo de los polinomes de aproximación de Tschebyschef" . compt. Desgarrar. Acad. Ciencia. (en francés). 199 : 337–340 .
  2. Chiang, Yi-Ling F. (noviembre de 1988). "Un algoritmo Remes modificado" . SIAM Journal on Scientific and Statistical Computing . 9 (6): 1058– 1072. doi : 10.1137/0909072 . ISSN 0196-5204 . 
  3. Fraser, W. (1965). "Un estudio de los métodos para calcular aproximaciones polinómicas minimax y casi minimax para funciones de una sola variable independiente" . J. ACM . 12 (3): 295–314 . doi : 10.1145/321281.321282 . S2CID 2736060 . 
  4. Kilgore, TA (1978). "Una caracterización de la proyección interpoladora de Lagrange con norma de Tchebycheff mínima". J. Approx. Theory . 24 (4): 273– 288. doi : 10.1016/0021-9045(78)90013-8 .
  5. de Boor, C.; Pinkus, A. (1978). "Demostración de las conjeturas de Bernstein y Erdös sobre los nodos óptimos para la interpolación polinomial" . Journal of Approximation Theory . 24 (4): 289– 303. doi : 10.1016/0021-9045(78)90014-X .
  6. Luttmann, FW; Rivlin, TJ (1965). "Algunos experimentos numéricos en la teoría de la interpolación polinomial". IBM J. Res. Dev . 9 (3): 187– 191. doi : 10.1147/rd.93.0187 .
  7. Rivlin, TJ (1974). "Las constantes de Lebesgue para la interpolación polinomial" . En Garnir, HG; Unni, KR; Williamson, JH (eds.). Análisis funcional y sus aplicaciones . Lecture Notes in Mathematics. Vol. 399. Springer. pp. 422–437 . doi : 10.1007/BFb0063594 . ISBN   978-3-540-37827-3.
  8. Brutman, L. (1978). "Sobre la función de Lebesgue para la interpolación polinomial". SIAM J. Numer. Anal . 15 (4): 694– 704. Bibcode : 1978SJNA...15..694B . doi : 10.1137/0715046 .
  9. Günttner, R. (1980). "Evaluación de las constantes de Lebesgue". SIAM J. Numer. Anal . 17 (4): 512– 520. Bibcode : 1980SJNA...17..512G . doi : 10.1137/0717043 .
  10. Luenberger, DG; Ye, Y. (2008). «Métodos básicos de descenso» . Programación lineal y no lineal . Serie internacional en investigación operativa y ciencias de la gestión. Vol. 116 (3.ª ed.). Springer. pp. 215–262 . doi : 10.1007/978-0-387-74503-9_8 . ISBN    978-0-387-74503-9.
  11. ^ Egidi, Nadaniela; Fatone, Lorella; Misici, Luciano (2020), "Un nuevo algoritmo tipo Remez para la mejor aproximación polinómica" , en Sergeyev, Yaroslav D.; Kvasov, Dmitri E. (eds.), Cálculos numéricos: teoría y algoritmos , vol. 11973, Cham: Springer, págs. 56 a 69, doi : 10.1007/978-3-030-39081-5_7 , ISBN   978-3-030-39080-8, S2CID 211159177 
  12. 1 2 Temes, GC; Barcilon, V.; Marshall, FC (1973). "La optimización de sistemas con ancho de banda limitado". Actas del IEEE . 61 (2): 196– 234. Bibcode : 1973IEEEP..61..196T . doi : 10.1109/PROC.1973.9004 . ISSN 0018-9219 . 
  13. Dunham, Charles B. (1975). "Convergencia del algoritmo de Fraser-Hart para la aproximación racional de Chebyshev" . Matemáticas de la Computación . 29 (132): 1078– 1082. doi : 10.1090/S0025-5718-1975-0388732-9 . ISSN 0025-5718 . 
  • Aproximaciones minimax y el algoritmo de Remez , capítulo introductorio en la documentación de Boost Math Tools, con enlace a una implementación en C++.
  • Introducción al DSP Archivado el 23/04/2014 en Wayback Machine
  • Aarts, Ronald M .; Vínculo, Carlos; Mendelsohn, Phil y Weisstein, Eric W. "Algoritmo Remez" . MundoMatemático .