Algorytm Remeza

Algorytm Remeza jest dobrze znaną procedurą iteracyjną przybliżającą funkcję wielomianem w normie minimax. Ale, jak mówi o tym Nick Trefethen [1]:

Większość tych [wdrożeń] sięga wielu lat wstecz, a właściwie większość z nich nie rozwiązuje ogólnego problemu najlepszego przybliżenia przedstawionego powyżej, ale warianty obejmujące zmienne dyskretne lub cyfrowe filtrowanie. W obiegu znajduje się kilka innych programów komputerowych, ale ogólnie wydaje się, że obecnie nie ma powszechnie stosowanego programu do obliczania najlepszych przybliżeń.

Można obliczyć rozwiązanie minimax również przez zastosowanie optymalizacji najmniejszych kwadratów lub wypukłej, na przykład za pomocą Matlaba i darmowego zestawu narzędzi CVX zastosowanego do funkcji Runge na [-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

Przybliżenie wielomianów Czebyszewa ma 0.1090minimalną normę, podczas gdy to podejście osiąga minimum 0.0654, tej samej wartości, która jest obliczana za pomocą algorytmu Remeza w chebfunprzyborniku Matlaba .

Czy jest jakaś korzyść ze stosowania bardziej skomplikowanego algorytmu Remeza, jeśli można szybciej i dokładniej obliczyć rozwiązanie minimax za pomocą solvera optymalizacyjnego? Czy są jakieś raporty / artykuły porównujące te dwa podejścia do niektórych trudnych problemów lub przypadków testowych?

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

„Właściwa” odpowiedź silnie zależy od tego, do czego potrzebujesz swojego przybliżenia. Czy naprawdę potrzebujesz najlepszego przybliżenia dla niektórych błędów? Czy po prostu dobre przybliżenie? Czy po prostu dobre przybliżenie w sensie minmax?

Nick Trefethen ostatnio podał ładny przykład, w którym przybliżenie Remeza jest złym pomysłem, ponieważ minimalizuje maksymalny błąd niezależnie od średniego błędu w całym przedziale czasowym, co może nie być tym, czego chcesz. Oczywiście maksymalny błąd może być duży, ale jest to ograniczone dla płynnych funkcji.

Aktualizacja

Po dyskusji w komentarzach poniżej pobrałem CVX Toolbox i wykonałem bezpośrednie porównanie z algorytmem Chebfun Remez (zastrzeżenie: jestem częścią zespołu programistycznego 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' );

Po wielu wynikach otrzymuję na moim laptopie z Matlab 2012a wersję CVX 1.22 i najnowszą migawkę SVN firmy Chebfun:

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

Zauważ, że Chebfun fużyty do pomiaru błędu ma dokładność do 15 cyfr. Remez Chebfuna zajmuje dwa razy więcej czasu, ale dostaje mniejszy błąd globalny. Należy zauważyć, że CVX używa skompilowanego kodu do optymalizacji, podczas gdy Chebfun jest w 100% natywnym Matlabem. Błąd minimalny to błąd 0.00654minimalny „na siatce”, poza tą siatką błąd może być do 0.00659. Zwiększenie rozmiaru siatki do m = 1001uzyskania

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

tj. prawie taka sama prędkość, ale optymalizacja dyskretna jest jeszcze gorsza od czwartej cyfry dziesiętnej. Wreszcie, zwiększając rozmiar siatki dalej do m = 10001mnie

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

tzn. optymalizacja dyskretna jest teraz ponad dziesięć razy wolniejsza i wciąż gorsza od szóstej cyfry.

Najważniejsze jest to, że Remez zapewni ci globalnie optymalny wynik. Chociaż dyskretny analog może być szybki na małych siatkach, nie da poprawnego wyniku.