A stability barrier for reconstructions from Fourier samples

0.87 MB·English

by  Ben Adcock

A stability barrier for reconstructions from Fourier samples
Ben Adcock
Department of Mathematics
Purdue University
2 150 N. University Street
1
West Lafayette, IN 47907
0
2 USA

t
c Anders C. Hansen
O
DAMTP, Centre for Mathematical Sciences
9 University of Cambridge
2 Wilberforce Rd, Cambridge CB3 0WA

United Kingdom
]
A
Alexei Shadrin
N
DAMTP, Centre for Mathematical Sciences
.
h University of Cambridge
t
a Wilberforce Rd, Cambridge CB3 0WA
m
United Kingdom
[

1 October 31, 2012
v
1
3 Abstract
8
We prove that any stable method for resolving the Gibbs phenomenon – that is, recovering
7
high-orderaccuracyfromthefirstmFouriercoefficientsofananalyticandnonperiodicfunction
.
0 –canconvergeatbestroot-exponentiallyfastinm. Anymethodwithfasterconvergencemust
1
alsobeunstable,andinparticular,exponentialconvergenceimpliesexponentialill-conditioning.
2
ThisresultisanalogoustoarecenttheoremofPlatte,Trefethen&Kuijlaarsconcerningrecovery
1
from pointwise function values on an equispaced m-grid. The main step in our proof is an
:
v estimate for the maximal behaviour of a polynomial of degree n with bounded m-term Fourier
i series. This bound is related to a conjecture of Hrycak & Gro¨chenig. In the second part of the
X
paper we discuss the implications of our main theorem to polynomial-based interpolation and
r
least-squaresapproachesforovercomingtheGibbsphenomenon. Finally,weproposetheuseof
a
so-calledFourierextensionsasanattractivealternativeforthisproblem. Wepresentnumerical
results demonstrating rapid convergence of the resulting approximation in a stable manner.
1 Introduction
The Fourier series
F f(x)= √1 (cid:88) fˆeijπx, fˆ = √1 (cid:90) 1 f(x)e−ijπxdx,
m j j
2 2
−1
|j|≤m
of an analytic and periodic function f : [−1,1] → R converges exponentially fast in the truncation
parameter m. However, such rapid convergence is destroyed once periodicity is no longer present.
1
Fornonperiodicf,theseriesF f(x)convergesonlylinearlyincompactsubsetsof(−1,1),andthere
m
isnouniformconvergenceon[−1,1]. Neartheendpointsx=±1onewitnessesthewell-knownGibbs
phenomenon [34, 49].
Although the Gibbs phenomenon has a long history [38], and is completely understood mathe-
matically [34, 49], it remains a significant hurdle in practical applications of Fourier series. Indeed,
itisatestamenttoitsimportancethatthequestionoftheresolutionoftheGibbsphenomenon–i.e.
restoring high-order convergence given only the partial Fourier series F f of a function – remains
m
an area of active inquiry.
Although there are many existing methods for this problem (see Section 2 for a review), no one
stands out as being inherently superior. In particular, many exponentially convergent strategies
appear to suffer from some sort of ill-conditioning, and whilst there are other methods that resolve
the Gibbs phenomenon in a stable manner, these typically result in slower orders of convergence.
The purpose of this paper is to explain these observations.
Specifically, we show that any exponentially-convergent method must also be exponentially ill-
conditioned, and in general, if a method has a convergence rate of ρ−mτ for some ρ > 1 and
τ ∈ (1,1], then it must also possess ill-conditioning of order ρm2τ−1. In particular, this result
2
implies the following fundamental stability barrier: the best possible convergence rate for a stable
method for resolving the Gibbs phenomenon is root-exponential in m.
A theorem of this type is not new. Our main result is a direct analogue of a theorem of Platte,
Trefethen & Kuijlaars for the problem of overcoming the Runge phenomenon in equispaced poly-
nomial interpolation. In [45] it was proved that any method for recovering high accuracy from the
pointwisevaluesofafunctiononanequispacedgridmustalsoexhibittheaforementionedinstability
behaviour. The problem we consider in this paper, namely recovering high accuracy from the first
2m+1 Fourier coefficients of f, can be considered a continuous analogue of this problem. Indeed,
the problem of recovery from equispaced pointwise values is equivalent to that of recovery from
discrete Fourier coefficients.
In proving our main result we follow a similar argument to that of [45]. The key step therein
is the use of an estimate of Coppersmith & Rivlin [24] concerning the maximal behaviour of an
algebraic polynomial of degree n which is bounded on an equispaced grid of m points. Our main
theoretical contribution is an analogous result for Fourier series: namely, we estimate the maximal
behaviour of a polynomial p of degree n with bounded m-term Fourier series F p. In doing so, we
m
provide a partial answer to a conjecture of Hrycak & Gr¨ochenig [39] (see Section 5 for a discussion).
Our main theorem on the behaviour of polynomials with bounded Fourier sums has an impor-
tant consequence for a particular method that is sometimes used in practice. Our result implies
that the so-called inverse polynomial reconstruction method (IPRM) [41, 42, 44] is exponentially
unstable. Moreover, although it is possible to stabilize this method via a least squares procedure
(henceforth referred to as the polynomial least squares method), our main theorem demonstrates
that this necessarily decreases the convergence rate to root-exponential.
Although root-exponential convergence is the best possible permitted, our theorem says nothing
about spectral convergence. Nor does is apply to methods for which convergence occurs only down
to some finite tolerance (cid:15) . In the final part of this paper we propose the use of so-called Fourier
tol
extensions [7,8,16,21,40]forthisproblem. Aswedemonstratebynumericalexample,thisapproach
gives rapid convergence – sometimes exponential, but always spectral – but only down to a finite
tolerance on the order of machine precision. Moreover, we show that it typically outperforms the
aforementioned polynomial least squares method. We conclude that Fourier extensions present an
attractive method for this problem.
The outline of the remainder of this paper is as follows. In Section 2 we give an overview
of common methods for the resolution of the Gibbs phenomenon, and discuss the key issues of
convergence and stability. The aforementioned estimate is obtained in Section 3, and in Section 4
wepresentthemainresultofthepaper. InSection5weconsiderpolynomialinterpolationandleast
squares, and Fourier extensions are introduced in Section 6.
2
2 The Gibbs phenomenon and its resolution
The Gibbs phenomenon has a long history dating back to Wilbraham in 1848 [50] (to acknowledge
thecontributionofWilbraham,thenameGibbs–Wilbrahamphenomenonisalsooccasionallyused).
Forgotten for half a century, this phenomenon was rediscovered by Michelson [43], with the ensuing
debate regarding convergence, or lack thereof, between Michelson and Love (carried out in Nature)
being eventually settled by Gibbs [29, 30] in 1899, after the arbitration of Poincar´e. Gibbs contri-
bution to this problem was first recognised by Bˆocher in 1906 [13], who introduced the term the
Gibbs phenomenon. A detailed and fascinating review of the Gibbs phenomenon and its history is
provided in [38], with shorter summaries also presented in [22, 34].
Many methods have been proposed to ameliorate or resolve the Gibbs phenomenon. Of these,
perhaps the earliest to appear were filters and mollifiers. Here the Gibbs phenomenon is viewed as
noise polluting the high-order Fourier coefficients, which can therefore be mitigated by premultipli-
cation with a rapidly decaying function [34, 49]. Unfortunately, standard filters do not lead to high
uniform accuracy: they only ensure faster convergence in regions of [−1,1] away from the endpoints
x = ±1. More recently, attention has focused on adaptive filters and mollifiers (see [49] and refer-
ences therein), which can be constructed to obtain exponential convergence in (−1,1) in a stable
manner, with typically polynomial accuracy at the endpoints. Note that this does not contradict
the main result of this paper, since the rate of uniform convergence is not exponential. As a general
principle, exponential convergence away from x=±1 can be obtained without ill-conditioning.
An alternative approach (which can also be viewed as a type of mollifier [49]) is the technique of
spectralreprojection,introducedanddevelopedbyGottliebetal[32,33,34,35,36]. Heretheslowly
convergentFourierseriesisreprojectedontoasuitablebasis;theso-calledGibbscomplementarybasis.
OftenGegenbauerpolynomialsareusedforthisbasis. Providedthevariousparametersarecarefully
tuned, exponential convergence, uniformly on [−1,1], can be restored. Unfortunately, the original
Gegenbauer reconstruction procedure of [35] has been shown to be rather sensitive to the choice
of parameters [17, 28], with the wrong parameters giving potentially divergent approximations.
To mitigate this effect, a substantially more robust procedure, based on Freud polynomials, was
introducedin[28]. Nonetheless,ourmainresultstatesthatspectralreprojectionmusteitherexhibit
exponentialinstabilityornotbetrulyexponentiallyconvergent. Wenoteatthispointthatastability
analysisofspectralreprojectionhasnotyetbeencarriedout. Nevertheless,thisapproachhasfound
application in a number of areas, including image processing [9, 10] and the spectral approximation
of PDEs with discontinuous solutions [37].
Spectralreprojectionissometimesreferredtoadirect technique,sinceitdoesnotrequiresolution
of a linear system. The most obvious inverse method is to seek to ‘interpolate’ the first 2m+1
Fourier coefficients of f with an algebraic polynomial of degree 2m + 1. This technique, which
requires solution of a linear system of equations, is sometimes referred to as the inverse polynomial
reconstructionmethod (IPRM)inliterature[41,42,44]. However,interpolationofFouriercoefficients
canbeseenasacontinuousanalogueofpolynomialinterpolationonequispacednodes. Itshouldcome
as little surprise, therefore, that there are substantial issues with both convergence and stability. In
particular,aRunge-typephenomenoniswitnessed. See[4,39],aswellasSection5,foradiscussion.
Since‘interpolating’Fouriercoefficientsmaynotwork,onecanalsousealowerdegreepolynomial
in combination with a least squares fit (so-called polynomial least squares). This was first discussed
in detail in [39], and later in [4]. Unfortunately, as we shall later prove, to ensure stability one
√
requires that the degree n scale like m. This corresponds to only root-exponential convergence in
m, consistently with the stability barrier we establish in Section 4. Nonetheless, despite this slower
convergence, we remark that this approach often outperforms spectral reprojection in practice. For
a comparison, see [4].
As an alternative to lowering the polynomial degree, one may also try using a higher degree
polynomial. Underdetermined least squares leads to a poor approximation in this case. However,
better accuracy can be restored by using l1-minimization instead. To the best of knowledge no
3
analysiscurrentlyexistsforthisapproach. Forarelateddiscussioninthecaseofequispacedfunction
values, see [19, 45]. One may also consider Sobolev norm minimization, such has been considered in
the equispaced case in [23].
In [4] a general framework was introduced for stable reconstructions in Hilbert spaces. Given
Fourier coefficients, one can reconstruct in any other basis of functions, with one example being the
polynomial least squares method discussed above. An alternative to polynomials involves the use
of splines. As discussed in [6], fixed-order splines result in algebraic convergence in a stable manner
(see also [51]). One may also consider variable-order splines, but stability becomes an issue.
A different approach to overcoming the Gibbs phenomenon is to smooth the function f via
subtraction so as to make it periodic up to a given order, and then compute its Fourier series. This
idea dates back to Krylov and Lanczos, amongst others, and was later studied by Lax and Gottlieb
& Orszag – see [1, 2, 46] and references therein. Such smoothing can be carried out implicitly,
via extrapolation on the high-order Fourier coefficients; an approach sometimes known as Eckhoff’s
methodinliterature[26]. Asshownin[46], thismethodconvergesalgebraicallyfastinm. However,
there are also issues with instability [1]. A hybrid approach, combining Gegenbauer reconstruction
and Eckhoff’s method, was also developed in [27].
Alternative methods for overcoming the Gibbs phenomenon arise from sequence extrapolation
techniques. See [12, 20]. In a similar spirit, Driscoll & Fornberg introduced Pad´e-based method in
[25]. This approach gives exponential convergence [11], however in view of our theorem, must also
be exponentially unstable. Such instability was noted in [25].
There are numerous other methods for resolving the Gibbs phenomenon, and we have not pre-
sented a complete list. The reader is referred to [18, 34, 49] and references therein for further
information. We also mention that many methods designed for the related problem of recovering
high accuracy from function values on equispaced grids (i.e. overcoming the Runge phenomenon)
can potentially be adapted to the Fourier coefficient problem. See [19, 45] for a comprehensive list
of such methods.
One such technique that has recently been successfully applied to overcome the Runge phe-
nomenonisthatofFourierextensions. Inthefinalpartofthispaper(Section6)weproposetheuse
of this approach to overcome the Gibbs phenomenon. As we show, the corresponding method can
be extremely effective for this problem.
3 Maximal behaviour of an algebraic polynomial with bounded
Fourier series
Thefirststeptowardsthestabilitytheoremin[45]isanestimateduetoCoppersmith&Rivlin[24].
This concerns the behaviour of the quantity
(cid:110) (cid:111)
A =sup (cid:107)p(cid:107) :p∈P ,(cid:107)p(cid:107) =1 ,
n,m ∞ n m,∞
where (cid:107)p(cid:107) = sup |p(x)|, (cid:107)p(cid:107) = max |p(x )| and {x ,...,x } is a grid of m
∞ x∈[−1,1] m,∞ j=1,...,m j 1 m
equispaced points in [−1,1]. In [24] it was shown that
(c1)nm2 ≤An,m ≤(c2)nm2. (3.1)
for constants c ≥ c > 1. This estimate determines how many equispaced gridpoints are required
2 1
to control the behaviour of a polynomial of degree n. Observe that if m=o(n2) then (3.1) implies
that there exists a polynomial which is bounded on the grid {x ,...,x }, but which grows large in
1 m
between grid points. On the other hand, A is bounded as n,m → ∞ if and only if m = O(cid:0)n2(cid:1)
n,m
(that the scaling m = O(cid:0)n2(cid:1) is sufficient for boundedness is an older result which dates back to
Sch¨onhage [48]).
4
The problem of the behaviour of A is a classical one in approximation theory (see [45] for
n,m
a summary of its history). However, A involves equispaced grid values. To prove a stability
n,m
theorem for overcoming the Gibbs phenomenon, we need to study a related quantity involving
Fourier coefficients. To this end, we define
B =sup{(cid:107)p(cid:107) :p∈P ,(cid:107)p(cid:107) =1}, (3.2)
n,m 2 n m
(cid:113)
where (cid:107)p(cid:107) = (cid:82)1 |p(x)|2dx is the L2-norm on [−1,1] and
2 −1
(cid:107)p(cid:107)2 = (cid:88) |pˆ |2 =(cid:107)F p(cid:107)2, (3.3)
m j m
|j|≤m
is the l2-norm of the first 2m+1 Fourier coefficients of p.
WeremarkthatthequantityB differsfromA intworespects. First,itinvolvescontinuous
n,m n,m
Fourier coefficients, as opposed to discrete Fourier coefficients (i.e. equispaced grid values). Second,
rather than using the uniform norm, we consider the L2 (respectively l2)-norm. This is natural,
since Fourier series satisfy Parseval’s relation in the L2/l2-norms. Indeed, the second equality in
(3.3) follows directly from this relation.
With this aside, note that
lim B =1,
n,m
m→∞
for any fixed n ∈ N. This follows from strong convergence of the operators F → I on L2(−1,1),
m
and the fact that the space P is finite dimensional. Moreover, much as in the case of A , it can
n n,m
be shown that the scaling m = O(cid:0)n2(cid:1) is sufficient for boundedness of B . This result was first
n,m
proved by Hrycak & Gr¨ochenig [39, Thm. 4.1] (see also [4, Lem. 3.1]).
Our main result in this section shows that B admits a similar lower bound to that found in
n,m
(3.1). We have
Theorem 3.1. Let B be as in (3.2). Then there exists a constant c>1 such that
n,m
Bn,m ≥cnm2, ∀n,m∈N.
Specifically, we have the lower bound
(Bn,m)2 ≥1+ 8nm + 16nmdnm2, (3.4)
where d> 9.
4
Proof. Let p∈P be arbitrary. Integrating by parts, we find that
n
(cid:88)n (−1)j
p = b , j ∈Z\{0},
(cid:98)j k jk
k=1
where
1 (cid:104) (cid:105)
b =−√ p(k−1)(1)−p(k−1)(−1) .
k
2(iπ)k
Therefore, we can write
n
p =(−1)jp(1), j ∈Z\{0}, p(t):=(cid:88)b tk. (3.5)
(cid:98)j (cid:101) j (cid:101) k
k=1
5
Note that p ∈ P is uniquely defined by the values pˆ ,b ,...,b . Therefore, (3.5) defines a one-to-
n 0 1 n
one correspondence between polynomials p=p(x) with pˆ =0, say, and polynomials p˜(t) satisfying
0
p˜(0)=0. Hence, using Parseval’s relation, we have
B ≥sup{(cid:107)p(cid:107) :p∈P ,pˆ =0,(cid:107)p(cid:107) =1}= sup B(n,m,p),
n,m 2 n 0 m
p∈P∗
n
where P∗ ={p∈P :p(0)=0} and
n n
(cid:118)
(cid:117)(cid:80) |p(1)|2
B(n,m,p)=(cid:117)(cid:116) 1≤|j|<∞ j .
(cid:80) |p(1)|2
1≤|j|≤m j
To prove the theorem it is sufficient to find a particular P ∈ P∗ such that B(n,m,P) admits the
n
bound (3.4).
Let n=4q+1 for some q ∈N. Consider the polynomial
P(x):=xT∗(x2)A (x2),
q q
where
q
A (x2):= (cid:89)(x2− 1 ),
q j2
j=1
(cid:104) (cid:105)
and T∗(x) denotes the Chebyshev polynomial of degree q on the interval 1 , 1 . Then, by
q m2 (q+1)2
definition,
1
P( )=0, 1<|j|≤q,
j
and therefore we have
(cid:80) |P(1)|2 (cid:80) |P(1)|2 (cid:80) |P(1)|2
(B(n,m,P))2 = 1≤|j|<∞ j = q<|j|<∞ j =1+ m<j<∞ j (3.6)
(cid:80) |P(1)|2 (cid:80) |P(1)|2 (cid:80) |P(1)|2
1≤|j|≤m j q<|j|≤m j q<j≤m j
Note that A (·) has all its zeros in the interval [ 1 ,1], hence |A | is monotonically decreasing on the
q q2 q
interval [0, 1 ]. Therefore,
q2
1 1 1
0<|A ( )|≤|A ( )|<|A ( )|, q <j ≤m<j .
q j2 q m2 q j2 1 2
1 2
So, putting expression for P into (3.6), we obtain
(cid:80) 1 |T∗( 1 )|2|A ( 1 )|2
(B(n,m,P))2 =1+ m<j<∞ j2 q j2 q j2
(cid:80) 1 |T∗( 1 )|2|A ( 1 )|2
q<j≤m j2 q j2 q j2
1
(cid:80) 1 |T∗( 1 )|2|A ( 1 )|2
>1+ m<j<∞ j2 q j2 q m2
(cid:80) 1 |T∗( 1 )|2|A ( 1 )|2
q<j≤m j2 q j2 q m2
(cid:80) 1 |T∗( 1 )|2
=1+ m<j<∞ j2 q j2
(cid:80) 1 |T∗( 1 )|2
q<j≤m j2 q j2
N
:=1+
D
1 1
Since |T |≤1 on the interval [ , ], for the denominator D we have the estimate
q m2 (q+1)2
D = (cid:88) 1 |T∗( 1 )|2 ≤ (cid:88) 1 ,
j2 q j2 j2
q<j≤m q<j<∞
6
so that D < 1. As for the numerator N, we split it into two parts:
q
N := (cid:88) 1 |T∗( 1 )|2, N := (cid:88) 1 |T∗( 1 )|2.
1 j2 q j2 2 j2 q j2
m<j<2m 2m≤j<∞
By definition, T∗ is the Chebyshev polynomial of degree q on the interval [ 1 , 1 ], so all of
q m2 (q+1)2
its zeros are located inside that interval. Hence |T∗| is decreasing on [0, 1 ] towards the value
q m2
|T∗( 1 )|=1. Therefore
q m2
N >(cid:12)(cid:12)T∗(cid:16) 1 (cid:17)(cid:12)(cid:12)2 (cid:88) 1 > 1 ,
1 (cid:12) q m2 (cid:12) j2 2m
m<j<2m
N >(cid:12)(cid:12)T∗(cid:16) 1 (cid:17)(cid:12)(cid:12)2 (cid:88) 1 > 1 (cid:12)(cid:12)T∗(cid:16) 1 (cid:17)(cid:12)(cid:12)2 .
2 (cid:12) q (2m)2 (cid:12) j2 2m(cid:12) q (2m)2 (cid:12)
2m≤j<∞
i.e.
1 1 (cid:12) (cid:16) (cid:17)(cid:12)2
N > , N > (cid:12)T∗ 1 (cid:12) .
1 2m 2 2m (cid:12) q (2m)2 (cid:12)
Let us evaluate the value |T∗(x∗)|, where x∗ := 1 . Take the affine mapping
q 0 0 (2m)2
M :{I∗ =[ 1 , 1 ]}→{I =[−1,1]},
m2 (q+1)2
so that
T∗(x)=T (M(x)) ⇒ T∗(x∗)=T (M(x∗))=T (x ),
q q q 0 q 0 q 0
where T (x) is the standard Chebyshev polynomial on I = [−1,1]. The length of the interval I∗ is
q
less than 1 , so mapping M onto the interval I with the length 2 uses the length magnifying factor
q2
λ>2q2. The point x∗ = 1 lies at the distance δ∗ = 3 from the left endpoint 1 of I∗, so it
0 (2m)2 0 4m2 m2
will be mapped to the point x <−1 given by
0
−x =1+δ , δ >2q2δ∗ = 3q2 .
0 0 0 0 2m2
For x=1+δ, we have
T (x)= 1(cid:16)(x+(cid:112)x2−1)q+(x−(cid:112)x2−1)q(cid:17)> 1(x+(cid:112)x2−1)q > 1(1+√2δ)q,
q 2 2 2
√ √
and since 2δ = 3q > q , we have
0 m m
T∗( 1 )=T∗(x∗)=T (x )> 1(cid:16)1+ q (cid:17)q = 1(cid:104)(cid:0)1+ q (cid:1)mq (cid:105)q2/m := 1γq2/m. (3.7)
q (2m)2 q 0 q 0 2 m 2 m 2
Hence
1 m
N > γ2q2/m, γ :=(1+ 1)r, r := . (3.8)
2 4m r q
Combing these results together, and using the fact that q ≈ n, we obtain
4
N N N
(B(n,m,P))2 >1+ =1+ 1 + 2
D D D
q q
>1+ + γ2q2/m
2m 4m
n n
>1+ + γn2/8m.
8m 16m
Since m>n/2 and q <n/4, we have r = m >2 so that a lower bound for γ is γ =(1+ 1)r >9/4.
q r
This completes the proof.
7
108 105
1015
104
1012 106
1000
109
104
106 100
100
1000 10
20 40 60 80 100 20 40 60 80 100 20 40 60 80 100
0.6 1.2
0.3 0.5 1.0
0.4 0.8
0.2
0.3 0.6
0.2 0.4
0.1
0.1 0.2
20 40 60 80 100 20 40 60 80 100 20 40 60 80 100
(α,β)=(1,1) (α,β)=(1,5) (α,β)=(1,3)
2 4 4 8 2
Figure 1: Toprow: thequantityB (squares)andthelowerboundB∗ (circles)againstn. Bottom
n,αnβ n,αnβ
row: thescaledquantitiesnβ−2log(B )andnβ−2log(B∗ ). ComputationswerecarriedoutinMath-
n,αnβ n,αnβ
ematica using additional precision. Note that the quantity B can be computed, since concides with the
n,m
minimum singular value of a particular matrix.
In Figure 1 we confirm this theorem by plotting the quantity B and the lower bound
n,m
(cid:118)
(cid:117)(cid:117) n n (cid:18)9(cid:19)nm2
B∗ =(cid:116)1+ + .
n,m 8m 16m 4
Note that B∗ not only provides a lower bound, it also appears to correctly predict the behaviour
n,m
of B . We conjecture that an upper bound of the form
n,m
Bn,m ≤cnm2, ∀n,m∈N, (3.9)
also holds, much as in the case of the quantity A .
n,m
4 An impossibility theorem for the resolution of the Gibbs
phenomenon
We now turn our attention to the main result of this paper. First, we require some notation.
Following[45],let{φm}m∈N beafamilyofmappingL2(−1,1)→L2(−1,1)suchthatφm(f)depends
onlyonthevalues{fˆ} . Notethatthemappingsφ canbebothlinearornonlinear. Wedefine
j |j|≤m m
the condition number for φ by
m
(cid:107)φ (f +g)−φ (f)(cid:107)
κ =suplim sup m m 2.
m f (cid:15)→0 g: (cid:107)g(cid:107)m
0<(cid:107)g(cid:107) ≤(cid:15)
m
For a compact set E ⊆C we shall also let B(E) be the Banach space of functions continuous on E
and analytic in its interior, with norm (cid:107)f(cid:107) =sup |f(z)|.
E z∈E
Our main result is as follows:
8
Theorem 4.1. Let a compact set E ⊆C containing [−1,1] in its interior be fixed, and suppose that
{φm}m∈N are such that (i) for any f, the quantity φm(f) depends only on the values {fˆj}|j|≤m, and
(ii) for some M <∞, σ >1 and τ ∈(1,1], we have
2
(cid:107)f −φ (f)(cid:107) ≤Mσ−mτ(cid:107)f(cid:107) , ∀f ∈B(E),m∈N. (4.1)
m 2 E
Then the condition numbers κ satisfy
m
κ ≥cm2τ−1,
m
for some c>1 and all sufficiently large m.
Proof. Weproceedasin[45]. WemaywithoutlossofgeneralityassumethatE isaBernsteinellipse
E(ρ) for some ρ>1. We may then replace (4.1) by
1
(cid:107)f −φ (f)(cid:107) ≤ ρ−αmτ(cid:107)f(cid:107) , m≥m , ∀f ∈B(E(ρ)). (4.2)
m 2 2 E 0
wherem issufficientlylargeandα>0isfixed. Letp∈P . AninequalityofBernstein[45,Lemma
0 n
1] implies that
(cid:107)p(cid:107) ≤ρn(cid:107)p(cid:107) .
E [−1,1]
Now, since p∈P , we have that (cid:107)p(cid:107) ≤cn(cid:107)p(cid:107) for some constant c>0. Thus,
n [−1,1] 2
(cid:107)p(cid:107) ≤cnρn(cid:107)p(cid:107) . (4.3)
E 2
Setting p=f in (4.2) now gives
(cid:16) (cid:17)
(cid:107)φ (p)(cid:107) ≥(cid:107)p(cid:107) − 1ρ−αmτ(cid:107)p(cid:107) ≥ 1− 1cnρn−αmτ (cid:107)p(cid:107) .
m 2 2 2 E 2 2
Suppose that n≤ 1αmτ. Then
2
cnρn−αmτ ≤ 21cαmτρ−21αmτ <1, ∀m≥m1,
where m is sufficiently large. Set m∗ :=max{m ,m }. If m≥m∗ this now gives
1 0 1
(cid:107)φ (p)(cid:107) ≥ 1(cid:107)p(cid:107) , ∀m≥m∗.
m 2 2 2
Since φ (0)=0 (this follows from (4.1)) we have
m
(cid:107)φ ((cid:15)p)−φ (0)(cid:107) (cid:107)φ ((cid:15)p)(cid:107) 1 (cid:107)(cid:15)p(cid:107) 1 (cid:107)p(cid:107)
m n 2 = m 2 ≥ 2 = 2 .
(cid:107)(cid:15)p(cid:107) (cid:107)(cid:15)p(cid:107) 2(cid:107)(cid:15)p(cid:107) 2(cid:107)p(cid:107)
m m m m
Therefore, since p∈P is arbitrary, we obtain
n
(cid:26) (cid:27)
(cid:107)p(cid:107)
κ ≥ 1 sup 2 = 1B(n,m), ∀n≤ 1αmτ,
m 2p∈Pn (cid:107)p(cid:107)m 2 2
p(cid:54)=0
where B(n,m) is given by (3.2). Using Theorem 3.1 with n= 1αmτ, this now yields
2
κm ≥ 12c14α2m2τ−1,
as required.
This theorem implies that any method for overcoming the Gibbs phenomenon that results in a
convergence rate of σ−mτ for all analytic functions f ∈ B(E) must also possess ill-conditioning of
order cm2τ−1. In particular, the best convergence rate that can be achieved with a stable method
is root-exponential in the number of Fourier coefficients m. As commented in [45], we stress that,
despite the use of polynomials in the proof, this theorem is not about polynomials or polynomial-
based approximation procedures. It holds for all (linear or nonlinear) mappings φ satisfying a
m
bound of the form (4.1).
9
5 Fourier coefficient interpolation and least squares
As discussed in Section 2, an obvious way to seek to overcome the Gibbs phenomenon is by inter-
polating the first 2m+1 Fourier coefficients of f with a polynomial of degree 2m. In other words,
we construct an approximation f satisfying
m
f(cid:99)mj =fˆj, |j|≤m, fm ∈P2m. (5.1)
It can be shown that such a polynomial exists uniquely for any m [39]. However, as we shall see
in a moment, this approach is both exponentially unstable and divergent. A simple modification
involves computing an overdetermined least squares with (n≤m):
(cid:40) (cid:41)
f =argmin (cid:88) |fˆ −pˆ |2 . (5.2)
n,m j j
p∈P2n |j|≤m
Note that f , as defined by (5.1), coincides with f when n = m. As mentioned, (5.1) is often
m n,m
referred to as the inverse polynomial reconstruction method (IPRM) [41, 42]. The modification
(5.2), referred to as polynomial least squares, was introduced and analysed in [39], and developed
further in [4] (see also [3]).
In [6] it was shown that f satisfies the sharp bound
n,m
(cid:107)f −f (cid:107)≤B inf (cid:107)f −p(cid:107),
n,m 2n,m
p∈P2n
whereB istheconstantdefinedin(3.2)(previous,butnon-sharp,estimatesweregivenin[4,39]).
n,m
Moreover, it was also shown that the condition number κ = κ of the mapping f (cid:55)→ f is
n,m n,m
precisely
κ =B , (5.3)
n,m 2n,m
where B is defined by (3.2). Hence, Theorem 3.1 allows us to explain both the convergence and
n,m
stability of this approach. In particular, (5.3) and Theorem 3.1 give that
κn,m ≥cnm2. (5.4)
√
and therefore polynomial least squares is unstable whenever n grows faster than m. In the partic-
ular case n=m, this shows that the IPRM method (5.1) is exponentially ill-conditioned. Note that
such exponential growth was previously observed, although not analysed, in [6, 42, 39].
We remark that Hrycak & Gr¨ochenig have conjectured that an upper bound of the form (5.4)
also holds. In other words, the condition numbers can grow no worse than cnm2 for some c>1. This
is of course equivalent to the conjecture (3.9) concerning the constant B .
n,m
Regarding the convergence of fn,m, it is useful to recall that the error infp∈P2n(cid:107)f −p(cid:107) decays
like ρ−2n, where ρ > 1 is the parameter of the largest Bernstein ellipse within which f is analytic
[47]. Thus, we have the bound
(cid:107)f −f (cid:107)≤c B ρ−2n, (5.5)
n,m f n,m
wherec >0issomeconstantdependingonf only. Thisindicatesthatf mayfailtoconvergeto
f n,m
f if the rate of growth of B exceeds that of ρ2n. In particular, if n=O(m), in which case B
n,m n,m
is exponentially large, there will always be functions (analytic within only a small ellipse E(ρ)) for
which the right-hand side of (5.5) diverges. Thus, the polynomial least squares method (5.2), and
in particular the IPRM (5.1), may well suffer from a Runge-type phenomenon – i.e. divergence of
f for some nontrivial family of analytic functions – whenever n = O(m). Although we have no
n,m
proof of this fact (the upper bound (5.5) need not be sharp for an individual function f), there is
substantial numerical evidence to support this conjecture [39, 42].
√
Ontheotherhand,whenn=O( m)itisknownthatB =O(1)[4,39],andthereforewehave
n,m
stability, as well as guaranteed convergence of the approximation. Unfortunately, the convergence
rate is only root-exponential, and Theorem 4.1 demonstrates that it can be no faster.
10