應(yīng)用數(shù)值數(shù)學(xué) 62 (2012) Pages 1790–1803.pdf

Applied Numerical Mathematics 62 2012 1790 1803 Contents lists available at SciVerse ScienceDirect Applied Numerical Mathematics Approximation error in regularized SVD based Fourier continuations Mark Lyon Department of Mathematics and Statistics University of New Hampshire 33 Academic Way Durham NH 03861 United States a r t i c l ei n f oa b s t r a c t Article history Received 20 December 2011 Received in revised form 5 April 2012 Accepted 24 June 2012 Available online 11 July 2012 Keywords Fourier series Fourier continuation Fourier extension We present an analysis of the convergence of recently developed Fourier continuation techniques that incorporates the required truncation of the Singular Value Decomposi tion SVD Through the analysis the convergence of SVD based continuations are related to the convergence of any Fourier approximation of similar form demonstrating the ef fi ciency and accuracy of the numerical method The analysis determines that the Fourier continuation approximation error can be bounded by a key value that depends only on the parameters of the Fourier continuation and on the points over which it is applied For any given distribution of points a fi nite number of calculations can be performed to ob tain this important value Our numerical computations on evenly spaced points show that as the number of points increases this quantity converges to a fi xed value allowing for broad conclusions on the convergence of Fourier continuations calculated with truncated SVDs We conclude that Fourier continuations can obtain super algebraic or even exponen tial convergence on evenly spaced points for non periodic functions until the convergence is limited by a parameter normally chosen near the machine precision accuracy threshold 2012 IMACS Published by Elsevier B V All rights reserved 1 Introduction Given the values of a function f with k 0 continuous derivatives evaluated at special sets of discrete points i e the Chebyshev points spectrally accurate methods for approximating functions from the discrete data are known If how ever evenly spaced or randomly spaced points are required the approximation algorithms are much more complex and generally yield signifi cantly slower convergence rates The exception being periodic functions where Fourier methods have outstanding properties and exponential rates of convergence for suffi ciently smooth functions When standard Fourier se ries methods are applied to non periodic functions a jump discontinuity is imposed in the approximation at the domain boundaries which gives rise to an adverse oscillation This is known as the Gibbs phenomenon which cripples the accuracy of the Fourier series near this artifi cial singularity A wide range of methods have been introduced to mitigate the ill effects of these oscillations A signifi cant early break through came in the Gegenbauer polynomials 19 20 27 and the subsequent utilization of Pad approximations 14 18 is also noteworthy Recently methods based on Fourier continuation FC 3 6 8 24 25 also termed Fourier extension 1 4 5 22 have demonstrated highly accurate approximations for non periodic functions The convergence of these FC meth ods will be the focus of this paper with the error of the resulting approximations quantifi ed through a combination of analytical and numerical study Fourier continuation uses a least squares approach to create an approximating Fourier series with a period that is larger than the domain of the original function The result is a Fourier series that effectively continues or extends the non periodic function An example is given in Fig 1 The function f x log 3x 1 on the interval 0 1 is shown as the solid line with a periodic Fourier continuation denoted by the dashed dotted line on the interval 0 2 From Fig 1 it is clear that E mail address mark lyon unh edu 0168 9274 36 00 2012 IMACS Published by Elsevier B V All rights reserved http dx doi org 10 1016 j apnum 2012 06 032 M Lyon Applied Numerical Mathematics 62 2012 1790 18031791 Fig 1 The function f solid defi ned by f x log 3x 1 is shown on the interval 0 1 and a periodic FC approximation to f is shown dashed on the interval 0 2 Note that f and its approximation are indistinguishable on 0 1 the these Fourier continuations are not analytic continuations and no connection to analyticity is meant to be implied by the terminology The method is a special case of Whitney s Extension Problem 21 29 that continues to be studied in much more general settings e g 15 16 and included references Fourier continuations have been shown through numerous numerical examples see 4 7 24 to produce super algebraically convergent error that decays faster than the number of Fourier modes to any real negative power Fourier series approximations throughout the domain of original function provided the original function is suffi ciently smooth The convergence is obtained despite the ill conditioning of the least squares system by using a Singular Value Decom position SVD and truncating singular values below a certain limit near machine precision This process is detailed in Section 2 The ill conditioning can be viewed as a consequence of the fact that the system has two complete bases which leads to non uniqueness in the solution as shown by Huybrechs in 22 through the use of frames When possible the system is split using even and odd symmetries to reduce the overall computational cost Fourier continuations have been successful in a variety of contexts including multi dimensional applications with un evenly spaced points such as those presented by Bruno et al in 7 In 8 25 Bruno and Lyon present an alternate but related continuation method FC Gram that relies upon a limited number of pre computed continuations of certain or thogonal polynomials This approach was implemented to reduce the computational cost of the SVD based FC approach of 4 6 7 and achieved a computational cost essentially equivalent to the cost of a single Fast Fourier Transform FFT O N log N The alternating direction Partial Differentiation Equation PDE solvers presented in 8 25 critically depend on the FC method s ability to produce accurate approximations from evenly spaced points The FC Gram approach was also extended by Albin and Bruno in 3 and applied to the compressible Navier Stokes equations in general spatial domains A method for high effi ciency parallel computing of FC based PDE solutions was also developed and successfully tested in 3 With the FC Gram technique the order of convergence corresponds to a degree of certain polynomials that are used in the continuation and therefore the resulting algorithms are not super algebraically convergent but high order accurate either fi fth or sixth order in practice Subsequently an algorithm has been introduced by Lyon in 24 to accelerate the computation of the SVD based FC methods toO N log N 2 or essentially the cost ofO log N FFTs The increased speed is obtained by a special splitting of the operator and uses FFT algorithms as well as randomized algorithms that have recently been introduced to compute decompositions for such low rank systems e g 12 13 17 23 30 These randomized algorithms have signifi cant similari ties to the compressed sensing methods that have provided great advances in applied mathematics over the past decade see 9 10 While slightly slower than the FC gram approximation the super algebraic convergence of SVD based method is obtained with the accelerated algorithm Further certain preliminary calculations need not be performed multiple times allowing effectively FFT speeds for any application requiring the repeated approximation from the same discrete grid such as would be the case for the PDE methods above In 22 Huybrechs showed under certain conditions on the function to be approximated and through the use of a specifi c set of discrete data points that are similar in nature to the Chebyshev points the exponential convergence of Fourier approximations that extend non periodic but analytic functions The result does not apply however to the many applications for which evenly spaced or alternate point distribution are required including the PDEs solvers in 3 8 25 In 7 Bruno et al bounds the L error of the FC approximations using certain numerically obtained constants These constants are effectively fi lling the role of a generalized Lebesgue constant for the approximation The constants however increase with the number of points at an unknown rate as the number of points and subsequently the number of Fourier modes are increased The present contribution can be viewed as an extended variation on the analysis in 7 Here the error of the FC approximations will be bound in L2instead of L resulting in numerically computed constants that remain bounded as the size of the system increases and allowing for broader conclusions concerning the convergence of the FC techniques The analysis presented here is also unique in that it will explicitly take the truncation of the singular values into account in the analysis adding insight into the stability and accuracy of the method 1792M Lyon Applied Numerical Mathematics 62 2012 1790 1803 Table 1 Quick reference table for important symbols used throughout the text ffunction to be approximated by Fourier continuation domain of the function f C f Fourier continuation of f N of data points in half of M of sine modes and cosine modes used inC f bperiod of the functionC f larger than the length of cdomain of a single period ofC f containing cut off for the regularization of the SVDs K of truncated singular values fM any Fourier series with M sine and M cosine modes 1 calculated value defi ned in Eq 20 2 calculated value defi ned in Eq 29 1upper bound on the value 1 Cconstant not bigger than limM 1 2M In the next section certain constants are derived that bound the L2error of FC approximations and then in Section 3 these parameters are evaluated and quantifi ed as the number of points increases From that analysis we conclude that super algebraic convergence and even in some cases exponential convergence can be obtained by FC methods that is until a near machine precision accuracy threshold which is determined by the truncation of the SVD is reached An important recent contribution by Platte et al proves that fast and stable exponentially accurate approximations of analytic but non periodic functions from evenly spaced points is impossible 26 With FC methods by picking appropriate continuation parameters and truncating the SVD of least squares systems the stability of the method is controlled The cost of the stability there must be a cost for the stability per 26 is not necessarily the loss of the rate of convergence but rather that full convergence to absolute zero cannot be obtained as the method will only obtain accuracies to a certain threshold Thus while exponential convergence can be observed with FC methods it is certainly not in violation of the theorem in 26 since the convergence is limited to a certain predetermined precision and the full class of functions for which exponential convergence can be obtained by FC methods has not yet been determined forthcoming work by Adcock et al 2 will address this issue If it were possible removing the accuracy threshold alone would result in exponential ill conditioning In practice however as this threshold can be chosen near machine precision FC methods are very attractive and should continue to gain traction for all practical computations Section 4 then gives numerical examples of functions that are non periodic over the domain they are sampled in and the resulting exponential convergence of the FC method to within an accuracy threshold of approximately 14 digits of accuracy is shown Section 5 briefl y addresses insights from the analysis on the effect of measurement and round off error in the computation of the FC approximations Finally we note that Fourier continuations depend on several parameters including the continuation length This paper does not attempt to determine optimal parameters as all tests to date indicate that optimal parameters depend on the function to be approximated as concluded by Bruno et al in 7 while the analysis here does not The contribution of this paper is in demonstrating the numerical stability of SVD regularized Fourier continuations That is the error produced by the Fourier continuation algorithm can be bounded in terms of the error produced by any similarly defi ned Fourier series approximation generated by any means computational or analytic For quick reference of important symbols see Table 1 2 Analysis of the split SVD based FC algorithm Let f be a function that possesses for some positive integer k a continuous kth derivative on the unit length interval 0 5 0 5 which for brevity we will denote by and let xj N j 1 be a collection of distinct points contained in that interval The Fourier continuation problem is then to defi ne using only the discrete data f xj for j 1 N a periodic approximation C f x m t M me 2 i b mx 1 to the function f with period b 1 where t M is the set of integers M 2 1 M 2 2 M 2 1 M 2 for even values of M or M 1 2 M 1 2 1 M 1 2 1 M 1 2 for odd values of M The coeffi cients m m t M can be determined e g 5 7 by solving the least square problem min m N j 1 m t M me 2 i b mxj f xj 2 2 M Lyon Applied Numerical Mathematics 62 2012 1790 18031793 using an SVD decomposition As mentioned in Section 1 the problem is often split into separate problems for sine and cosine coeffi cients by exploiting even and odd symmetries e g 4 5 7 24 The analysis presented here is performed using the split system while the analysis of the full system would be expected to produce similar results through a simplifi ed version of the ensuing analysis For a given function f we defi ne f e f x f x 2 and fo f x f x 2 then for a set of distinct points x j N j 1 on the interval 0 0 5 a Fourier continuation can be defi ned by C f x c1 2 M m 2 cmcos 2 b m 1 x M m 1 smsin 2 b mx 3 with the coeffi cients resulting from the minimizations min cm N j 1 c1 2 M m 2 cmcos 2 b m 1 x j f e x j 2 4 and min sm N j 1 M m 1 smsin 2 b mx j f o x j 2 5 This alternative formulation of the problem can be computed four times faster than solving the single equation 2 Oth erwise the precise formulation and manner of solution of the Fourier continuation has been observed to have relatively insignifi cant impact on the properties of the solution e g 7 The total number of modes in the split system is 2M and the function would need to be evaluated at either 2N or 2N 1 values depending on whether or not x j 0 for some j 2 1 Algorithm for Fourier continuation With fe f e x 1 f e x N T fo fo x 1 fo x N T c c1 cM T s s1 sM T Lej1 1 2 for j 1 N 6 Lejm cos 2 b m 1 x j for j 1 N m 2 M 7 and Lo jm sin 2 b mx j for j 1 N m 1 M 8 Eqs 4 and 5 are equivalent to the least squares solutions to the linear systems Le c fe andLo s fo 9 The resulting matrices Leand Loare ill conditioned for even modest values of N and M and require truncation of all singular values less than a given cut off to obtain an accurate solution in the presence of round off error Let UeSeVeTbe the full SVD decomposition see 28 of the matrix Le Le UeSeV eT 10 where the matrix Ueis N N and unitary Seis an N M diagonal matrix of the singular values and Veis an M M uni tary matrix and VeTdenotes the complex conjugate transpose of Ve The truncated SVD least squares solution is obtained by c V e Se U eT fe 11 where the M N diagonal matrix Se is formed by setting all values to zero except Se jj 1 Sejj for all j such that Se jj The analogous equation for the coeffi cients s s Vo So U oT fo 12 is obtained with the SVD decomposition Lo UoSoVoT As defi ned above there are four user choosable parameters for the continuation M N b and 1794M Lyon Applied Numerical Mathematics 62 2012 1790 1803 2 2 Preliminary analysis Our goal is to bound the error in the approximation of a given function f byC f in terms of the error of any arbitrary Fourier approximation of f We let fM be any function of the form fM c0 2 M m 2 cmcos 2 b m 1 x M m 1 smsin 2 b mx 13 and note that f C f L2 f fM L2 fM C fM L2 C f fM L2 14 In order to obtain a bound for f C f L2 we proceed by bounding the three simpler quantities f fM L2 fM C fM L2 and C f fM L2 in terms of the quantities f fM L and fM L c for any given Fourier series fM Since the bound will hold for all functions fM and Fourier series can easily be obtained which rapidly converge to a given suffi ciently differential function see Remark 2 1 in what follows we show that the numerically calculated Fourier continuation will also rapidly converge to within a certain accuracy threshold Also since we have taken to a unit length interval the trivial bound f fM L2 f fM L 15 holds and we need only consider the terms fM C fM L2 and C f fM L2 further Often when such bounds are constructed for other approximation methods the term corresponding to fM C fM L2 would be negligible since one might expect that a 2M mode Fourier series approximation of a 2M mode Fourier series with an identical basis would be perfect except for round off errors However this is not the case for SVD regularized Fourier continuations since certain combinations of modes are neglected by the regularization and this term will not in general be zero and the required analysis follows 2 3 Bound for fM C fM L2 The value fM C fM L2 will only be zero if the system is suffi ciently small so that no singular values in the SVDs are below the truncation cut off Yet the singular values of Leand Lodecay so rapidly as M increases that the truncation of certain modes is unavoidable for even modest sized computations The columns of Veand Voprovide alternate bases for analysis of the error Let 1 M T V eT c and 1 M T V oT s 16 then the coeffi cients jand jgive an alternate representation of the Fourier series for fM fM x M j 1 jvej x jvoj x 17 where the functions ve j x are each an up