The Remez Algorithm

Question

The Remez algorithm is a well-known iterative routine to approximate a function by a polynomial in the minimax norm. But, as Nick Trefethen [1] says about it:

Most of these [implementations] go back many years and in fact, most of them do not solve the general best approximation problem as posed above but variants involving discrete variables or digital filtering. One can find a few other computer programs in circulation, but overall, it seems that there is no widely-used program at present for computing best approximations.

One can compute the minimax solution also by applying least-squares or convex optimization, for example using Matlab and the free CVX toolbox applied to the Runge function on [-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

The approximation with Chebyshev polynomials has a minimax norm of 0.1090 while this approach here reaches a minimum of 0.0654, the same value that gets computed with the Remez algorithm in the Matlab chebfun toolbox.

Is there any advantage in applying the more complicated Remez algorithm if you can compute the minimax solution faster and more accurate with an optimization solver? Are there any reports/articles comparing these two approaches on some difficult problems or test cases?

--
[1] R. Pachon and L. N. Trefethen. BIT Numerical Mathematics (2008) Vol. 46.

Answer 1

The "right" answer strongly depends on what you need your approximant for. Do you really need the best approximation for some error bound? Or just a good approximation? Or just a good approximation in the minmax sense?

Nick Trefethen recently gave a nice example where Remez approximation is a bad idea since it minimizes the maximum error irrespective of the average error over the entire interval, which may not be what you want. Of course, the maximum error may be large, but this is bounded for smooth functions.

Update

Following the discussion in the comments below, I downloaded the CVX Toolbox and did a direct comparison with the Chebfun Remez algorithm (disclaimer: I am part of the Chebfun development team):

% 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' );

After a lot of output, I get, on my laptop with Matlab 2012a, CVX version 1.22 and Chebfun's latest SVN Snapshot:

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

Note that the Chebfun f used to measure the error is accurate to 15 digits. So Chebfun's Remez takes twice as long, but gets a smaller global error. It should be pointed out that CVX uses compiled code for the optimization whereas Chebfun is 100% native Matlab. The minimum error of 0.00654 is the minimum error 'on the grid', off that grid, the error can be up to 0.00659. Increasing the grid size to m = 1001 I get

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

i.e. almost the same speed, but the discrete optimization is still worse as of the fourth decimal digit. Finally, ncreasing the grid size further to m = 10001 I get

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

i.e. the discrete optimization is now more than ten times slower and is still worse as of the sixth digit.

The bottom line is that Remez will get you the globally optimal result. While the discrete analog may be fast on small grids, it will not give a correct result.