Algorytm Remeza

14

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

Hans W.
źródło

Odpowiedzi:

4

„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.

Pedro
źródło
N. Trefethen podkreślał to samo i podał podobny przykład w cytowanym przeze mnie artykule. Pytanie nie dotyczyło najlepszego przybliżenia, ale: Jaka jest zaleta algorytmu Remeza (obecnie), jeśli można uzyskać ten sam wynik za pomocą rozsądnego wypukłego rozwiązania ?
Hans W.
1
@HansWerner, przepraszam, źle odczytałem twoje pytanie. Twój wypukły solver nie daje tego samego wyniku, przynajmniej nie dla wszystkich cyfr. Jeśli dobrze rozumiem twój wypukły kod, minimalizujesz maksymalny błąd w stosunku do dyskretnego zestawu punktów. Jest to dobre przybliżenie - ale nie to samo - minimalizowanie globalnego błędu maksymalnego.
Pedro
Wypukły solver w tym przypadku dał lepszy wynik. Pomyśl o tym, algorytm Remeza jest procedurą iteracyjną, dość podobną do procedury optymalizacji i również nie zwróci dokładnego wyniku. W konkretnym przypadku powyżej rozwiązanie z optymalizacji było lepsze, tj. Miało ogólnie mniejszą maksymalną normę, niż wynik najlepszej implementacji Remeza, jaką znam. Pytanie jest nadal otwarte .
Hans W.
@HansWerner, jak zmierzyłeś maksymalny błąd rozwiązania uzyskany za pomocą wypukłego solvera? Algorytm Remeza chebfunpowinien iterować, dopóki nie zostanie osiągnięte minimum precyzji maszyny (w pewnym sensie).
Pedro
Niekoniecznie; istnieją opcje takie jak „tol” (czyli tolerancja względna) lub „maxiter” dla chebfun/remez, ale istnieją podobne opcje dla wszystkich rodzajów solverów optymalizacyjnych. Można powiedzieć, że Remez jest procedurą optymalizacji specjalizującą się w określonym zadaniu. A pytanie brzmi: czy solwery ogólnego przeznaczenia nadrobiły zaległości i nie ma już potrzeby tak wyspecjalizowanego solvera? Oczywiście mogę się mylić.
Hans W.