El algoritmo de Remez

14

El algoritmo Remez es una rutina iterativa bien conocida para aproximar una función por un polinomio en la norma minimax. Pero, como Nick Trefethen [1] dice al respecto:

La mayoría de estas [implementaciones] se remontan a muchos años y, de hecho, la mayoría de ellas no resuelven el problema general de mejor aproximación como se plantea anteriormente, sino variantes que involucran variables discretas o filtrado digital. Uno puede encontrar algunos otros programas de computadora en circulación, pero en general, parece que actualmente no hay un programa ampliamente utilizado para calcular las mejores aproximaciones.

También se puede calcular la solución minimax mediante la aplicación de mínimos cuadrados u optimización convexa, por ejemplo, usando Matlab y la caja de herramientas CVX gratuita aplicada a la función Runge en [-1, 1]:

m = 101; n = 11;            % 101 points, polynomial of degree 10
xi = linspace(-1, 1, m);    % equidistant points in [-1, 1]
ri = 1 ./ (1+(5*xi).^2);    % Runge function

tic                         % p is the polynomial of degree (n-1)
cvx_begin                   % minimize the distance in all points
    variable p(n);
    minimize( max(abs(polyval(p, xi) - ri)) );
cvx_end
toc                         % 0.17 sec for Matlab, CVX and SeDuMi

La aproximación con los polinomios de Chebyshev tiene una norma 0.1090mínima de 0.0654, mientras que este enfoque aquí alcanza un mínimo del mismo valor que se calcula con el algoritmo Remez en la chebfuncaja de herramientas de Matlab .

¿Hay alguna ventaja en aplicar el algoritmo Remez más complicado si puede calcular la solución minimax más rápido y más preciso con un solucionador de optimización? ¿Hay informes / artículos que comparen estos dos enfoques sobre algunos problemas difíciles o casos de prueba?

-
[1] R. Pachon y LN Trefethen. BIT Numerical Mathematics (2008) vol. 46)

Hans W.
fuente

Respuestas:

4

La respuesta "correcta" depende en gran medida de para qué necesita su aproximación. ¿Realmente necesita la mejor aproximación para algún error vinculado? O simplemente una buena aproximación? ¿O simplemente una buena aproximación en el sentido minmax?

Nick Trefethen recientemente dio un buen ejemplo donde la aproximación de Remez es una mala idea ya que minimiza el error máximo independientemente del error promedio en todo el intervalo, que puede no ser lo que desea. Por supuesto, el error máximo puede ser grande, pero esto está limitado para funciones suaves.

Actualizar

Después de la discusión en los comentarios a continuación, descargué el CVX Toolbox e hice una comparación directa con el algoritmo Chebfun Remez (descargo de responsabilidad: soy parte del equipo de desarrollo de Chebfun):

% Do the convex optimization bit.
m = 101; n = 11;            % 101 points, polynomial of degree 10
xi = linspace(-1, 1, m);    % equidistant points in [-1, 1]
ri = 1 ./ (1+(5*xi).^2);    % Runge function

tic                         % p is the polynomial of degree (n-1)
cvx_begin                   % minimize the distance in all points
    variable p(n);
    minimize( max(abs(polyval(p, xi) - ri)) );
cvx_end
toc_or = toc                % 0.17 sec for Matlab, CVX and SeDuMi

% Extract a Chebfun from the result
x = chebfun( [-1,1] );
A = [ chebfun(1) , x ];
for k=3:n, A(:,k) = A(:,k-1).*x; end
or = A * flipud(p)

% Make a chebfun of Runge's function
f = chebfun( @(x) 1 ./ ( 1 + 25*x.^2 ) )

% Get the best approximation using Remez
tic, cr = remez( f , 10 ); toc_cr = toc

% Get the maximum error in each case
fprintf( 'maximum error of convex optimization: %e (%f s)\n' , norm( f - or , inf ) , toc_or );
fprintf( 'maximum error of chebfun remez: %e (%f s)\n' , norm( f - cr , inf ) , toc_cr );

% Plot the two error curves
plot( [ f - cr , f - or ] );
legend( 'chebfun remez' , 'convex optimization' );

Después de una gran cantidad de resultados, obtengo, en mi computadora portátil con Matlab 2012a, CVX versión 1.22 y la última instantánea SVN de Chebfun:

maximum error of convex optimization: 6.665479e-02 (0.138933 s)
maximum error of chebfun remez: 6.592293e-02 (0.309443 s)

Tenga en cuenta que el Chebfun futilizado para medir el error tiene una precisión de 15 dígitos. Entonces, Remez de Chebfun toma el doble de tiempo, pero obtiene un error global más pequeño. Cabe señalar que CVX utiliza código compilado para la optimización, mientras que Chebfun es Matlab 100% nativo. El error mínimo de 0.00654es el error mínimo 'en la cuadrícula', fuera de esa cuadrícula, el error puede ser hasta 0.00659. Aumentar el tamaño de la cuadrícula para m = 1001obtener

maximum error of convex optimization: 6.594361e-02 (0.272887 s)
maximum error of chebfun remez: 6.592293e-02 (0.319717 s)

es decir, casi la misma velocidad, pero la optimización discreta es aún peor a partir del cuarto dígito decimal. Finalmente, aumentando el tamaño de la cuadrícula más para m = 10001obtener

maximum error of convex optimization: 6.592300e-02 (5.177657 s)
maximum error of chebfun remez: 6.592293e-02 (0.312316 s)

es decir, la optimización discreta ahora es más de diez veces más lenta y aún es peor a partir del sexto dígito.

La conclusión es que Remez le dará el resultado óptimo a nivel mundial. Si bien el análogo discreto puede ser rápido en redes pequeñas, no dará un resultado correcto.

Pedro
fuente
Y N. Trefethen hizo hincapié en lo mismo y dio un ejemplo similar en el artículo que estaba citando. La pregunta no era acerca de la mejor aproximación, sino: ¿Cuál es la ventaja del algoritmo Remez (hoy en día) si puede obtener el mismo resultado con un solucionador convexo razonable ?
Hans W.
1
@HansWerner, lo siento, leí mal tu pregunta. Su solucionador convexo no le está dando el mismo resultado, al menos no a todos los dígitos. Si entiendo su código convexo correctamente, está minimizando el error máximo en un conjunto discreto de puntos. Esta es una buena aproximación de, pero no es lo mismo que, minimizar el error máximo global.
Pedro
El solucionador convexo en este caso dio un mejor resultado. Piénselo, el algoritmo Remez es un procedimiento iterativo, bastante similar a una rutina de optimización, y tampoco devolverá un resultado exacto. En el caso concreto anterior, la solución de la optimización fue mejor, es decir, tenía una norma máxima general más pequeña que el resultado de la mejor implementación de Remez que conozco. La pregunta aún está abierta .
Hans W.
@HansWerner, ¿cómo midió el error máximo de la solución obtenida con el solucionador convexo? El algoritmo de Remez chebfundebería iterar hasta que se logre el mínimo de precisión de la máquina (en cierto sentido).
Pedro
No necesariamente; hay opciones como 'tol' (que es tolerancia relativa) o 'maxiter' para chebfun/remez, pero hay opciones similares para todo tipo de solucionadores de optimización. En cierto modo, se podría decir que Remez es una rutina de optimización especializada para una determinada tarea. Y la pregunta es: ¿Han resuelto los solucionadores de propósito general y ya no hay necesidad de un solucionador tan especializado? Por supuesto, puedo estar equivocado.
Hans W.