LW = 'linewidth'; lw = 1.6; MS = 'markersize'; FS = 'Fontsize'; fs = 16;
Least-squares approximation
If $f:[-1,1]\rightarrow R$ is an $L^2$-integrable function, then its least-squares or best $L^2$ approximation of degree $n$ is the polynomial $p_n$ of degree at most $n$ such that
$$ | f - p_n |_2 = \mbox{minimum}. $$
A good introduction to $L^2$ approximations can be found in [2]. The polyfit
command in Chebfun returns the best $L^2$ approximation of a given degree to a chebfun:
help chebfun/polyfit
POLYFIT Fit polynomial to a CHEBFUN. F = POLYFIT(Y, N) returns a CHEBFUN F corresponding to the polynomial of degree N that fits the CHEBFUN Y in the least-squares sense. If Y is a global polynomial of degree n then this code has an O(n (log n)^2) complexity. If Y is piecewise polynomial then it has an O(n^2) complexity. F = POLYFIT(X, Y, N, D), where D is a DOMAIN object, returns a CHEBFUN F on the domain D which corresponds to the polynomial of degree N that fits the data (X, Y) in the least-squares sense. X should be a real-valued column vector and Y should be a matrix with size(Y,1) = size(X,1). F = POLYFIT(Y, N) where Y is represented as a periodic TRIGFUN object returns the degree N trigonometric polynomial fit of length 2N+1. Note CHEBFUN/POLYFIT does not not support more than one output argument in the way that MATLAB/POLYFIT does. See also INTERP1.
The coefficients of $p_n$ in the Legendre basis can be computed by truncating the Legendre expansion for $f$ after $n+1$ terms. For example,
n = 5; x = chebfun('x'); f = abs(x); P = legpoly(0:n,[-1,1],'norm'); % Legendre-Vandermonde matrix cleg = P'*f; % compute Legendre coefficients pn = P*cleg; % form chebfun of best L^2 approximation plot(f,LW,lw), hold on, plot(pn,'r',LW,lw) title('Best L^2 approximation to |x| of degree 5','fontsize',16), hold off
This approach works well, but requires $O(n^2)$ operations to compute the Legendre coefficients with a relatively large constant. This was the algorithm using in Chebfun's polyfit
command for many years, but was changed last week.
polyfit
using fast Chebyshev-Legendre transform
Recently, the command cheb2leg
was added in Chebfun, which converts a vector of Chebyshev coefficients (of the first kind) to Legendre coefficients in $O(n(\log n)^2/\log\log n)$ operations [1,3]. The command leg2cheb
is its inverse. To compute the Legendre expansion of a function (accurate to machine precision) we can first compute its Chebyshev expansion and then use cheb2leg
. Chebfun already computes the Chebyshev expansion of a function that is accurate to machine precision. Therefore, here is another way to compute the best approximation of a smooth function via cheb2leg
:
n = 10; f = 1./(1+25*x.^2); % Runge function ccheb = chebcoeffs(f); % get the Chebyshev coefficients of f cleg = cheb2leg(ccheb); % convert Cheb coeffs of f to Leg coeffs cleg = cleg(1:(n+1)); % truncate ccheb = leg2cheb(cleg); % convert them back to form a chebfun pn = chebfun(ccheb,'coeffs'); % form a chebfun plot(f,LW,lw), hold on, plot(pn,'r',LW,lw) title('Best L^2 approx to Runge function of degree 10',FS,14), hold off
This is the algorithm that is used in Chebfun's polyfit
, as of today. So we can obtain the same result from the code:
n = 10; f = 1./(1+25*x.^2); % Runge function pn = polyfit(f,n); plot(f,LW,lw), hold on, plot(pn,'r',LW,lw) title('Best L^2 approx to Runge function of degree 10',FS,14), hold off
High-degree best $L^2$ approximation for smooth functions
The fast algorithms now employed by polyfit
enable us to compute very high degree $L^2$ approximations.
n = 1e4; f = 1./(1+1e6*x.^2); % Runge function s = tic; pn = polyfit(f,n); t = toc(s); fprintf('L^2 error is %1.3e\n',norm(f - pn)) fprintf('L^2 approximation of degree %u in t = %1.3f\n',n,t)
L^2 error is 1.797e-06 L^2 approximation of degree 10000 in t = 1.494
Piecewise smooth functions
Computing the Legendre coefficients for piecewise smooth functions is a little trickier. The Legendre coefficients are computed by quadrature rules and then a chebfun object is constructed via Chebyshev coefficients computed using the cheb2leg
command. The algorithm for piecewise smooth function requires $O(n^2)$ operations, but the implicit constant is much smaller. Here is the best $L^2$-approximation to the piecewise smooth absolute value function.
f = abs(x); nn = 10.^(0:3); j=1; for n = nn pn = polyfit(f, n); err(j) = norm(f - pn); j = j+1; end loglog(nn, err,'k.-',LW,lw,MS,24), hold on loglog(nn,nn.^(-3/2),'k--',LW,lw) legend('|| |x| - p_n ||_2','n^{-3/2}') xlabel('n',FS,fs), ylabel('|| |x| - p_n ||_2',FS,fs) title('Convergence of || |x| - p_n ||_2',FS,fs)
References
-
N. Hale and A. Townsend, A fast, simple, and stable Chebyshev-Legendre transform using an asymptotic formula, SIAM Journal on Scientific Computing, 32 (2014), A148-A167.
-
M. Powell, Approximation Theory and Methods, Cambridge University Press, 1981.
-
A. Townsend and N. Hale, A fast Chebyshev-Legendre transform, Chebfun Example, August 2013.