On derivatives of smooth functions represented in multiwavelet bases

Contents lists available atScienceDirect

Journal of Computational Physics: X

www.elsevier.com/locate/jcpx

On derivatives of smooth functions represented in multiwavelet bases ^✩

Joel Anderson

^a

, Robert J. Harrison

^a,∗

, Hideo Sekino

^a

, Bryan Sundahl

^a

, Gregory Beylkin

^b

, George I. Fann

^c

, Stig R. Jensen

^d

, Irina Sagert

^e

aInstituteforAdvancedComputationalScience,StonyBrookUniversity,StonyBrook,NY11794-5250,USA bDepartmentofAppliedMathematics,UniversityofColoradoatBoulder,UCB526,Boulder,CO80309-0526,USA

cComputerScienceandMathematicsDivision,OakRidgeNationalLaboratory,MS6211,1BethelValleyRoad,OakRidge,TN37831-6367,USA dDepartmentofChemistry,HylleraasCentreforQuantumMolecularSciences,UiTTheArcticUniversityofNorway,Mailbox6050,Langnes N-9037,Tromsø,Norway

eCenterforTheoreticalAstrophysics,LosAlamosNationalLaboratory,LosAlamos,NM,87545,USA

a r t i c l e i n f o a b s t r a c t

Articlehistory:

Received12December2018 Receivedinrevisedform28May2019 Accepted29May2019

Availableonline24July2019

Keywords:

Multiwavelets Multiresolution Derivatives Numerical Discontinuous Bandlimitedexponentials

We construct high-order derivative operators for smooth functions represented via discontinuous multiwavelet bases.The need for suchoperators arises in orderto avoid artifacts when computing functionals involving high-order derivatives of solutions of integral equations. Previously high-order derivatives had to be formed by repeated application of a first-derivative operator that, while uniquely defined, has a spectral norm that grows quadratically with polynomial order and, hence, greatly amplifies numericalnoise(truncationerror)inthemultiwaveletcomputation.Thenewconstructions proceedvialeast-squaresprojectionontosmoothbasesandprovidesubstantiallyimproved numerical propertiesas well aspermitting directconstruction ofhigh-order derivatives.

Weemployeitherb-splinesorbandlimitedexponentialsastheintermediatesmoothbasis, withtheformermaintainingtheconceptofapproximationorderwhilethelatterpreserves the pureimaginary spectrum of the first-derivative operator and providesmore direct control overthe bandlimitand accuracyofcomputation.Wedemonstratethe properties ofthesenewoperatorsviaseveralnumericaltestsaswellasapplicationtoaproblemin nuclearphysics.

©^{2019TheAuthors.PublishedbyElsevierInc.ThisisanopenaccessarticleundertheCC} BYlicense(http://creativecommons.org/licenses/by/4.0/).

1. Introduction

Inthispaperwerevisit theissueofcomputinghigh-order derivativesofsmoothfunctionsrepresentedinmultiwavelet bases. During the development of various flavors of wavelet bases, multiwavelet bases were introduced perhaps as an afterthought. Initially,the wholethrustofwaveletconstructionssuch asDaubechies wavelets[12] or splinewavelets(see e.g. [11])was to ensure smoothnessofa basis.Multiwavelets [1] appeared followinga discrete wavelet-likeconstruction in [2] (it turns out that the same bases were sketched earlier in [15] but never used for numerical purposes). While

✩ Anderson,HarrisonandSundahlweresupportedinpartbytheNSFundergrantACI-1450344.JensenwassupportedbytheNorwegianResearchCouncil throughtheCoEHylleraasCentreforQuantumMolecularSciencesGrantNo.262695.SagertwassupportedinpartbyDOEgrantsDE-FG02-87ER40365and DE-SC0018083.FannwassupportedbytheDepartmentofEnergyunderContractNo.DE-AC05-00OR22725.

*

Correspondingauthor.

E-mailaddress:robert.harrison@stonybrook.edu(R.J. Harrison).

https://doi.org/10.1016/j.jcpx.2019.100033

2590-0552/©^{2019TheAuthors.PublishedbyElsevierInc.ThisisanopenaccessarticleundertheCCBYlicense} (http://creativecommons.org/licenses/by/4.0/).

generalizing theHaar basis[17],multiwaveletskept the“unpleasant”featurethatsome basisfunctionsarediscontinuous.

Counter-intuitively, asnotedin[3],thisfeaturemaybeconsidered“a blessingindisguise”.Inparticular,thefactthat the scaling functionsaresupported onnon-overlapping intervalsallowshigh-orderschemesonboundedintervalsand, also,it isrelativelyeasytoperformnonlinearoperations.Fromanalgorithmicpointofview,itisalsorelativelyeasytomaintaina sparse representationoffunctionsinhigherdimensionswhencomparedwithrequiredbookkeepingforoverlappingbases.

The observationthatdiscontinuousbasis functionsallowhigh-orderschemesisnotuniquetomultiresolutionbasesasthe same hasbeenobserved in,e.g., thediscontinuous Galerkinmethods.Consequently,since representationsinmultiwavelet bases have manyusefulpropertiesandresult inalgorithms well suited formodern computer processors, they becamea successful practical tool for highperformance computing andare used in quantum chemistry [18,20,45,46] and nuclear physics[34,14,31,32],andserveasmathematicalunderpinningsofMADNESS[21].

Yet, thefact that basis functionsare allowed to bediscontinuous leads to anumberof problems.Asdiscussed in[3], the only scale-consistent derivative operator available to usis the first derivative. Operators for the second and higher derivatives do not exist dueto the discontinuous nature of thebasis andare applied asa power of the firstderivative.

Naively,ifthefunctionhashigh-orderderivatives,suchanapproachappearsreasonable.However, duetotheapproximate nature of representation in multiwavelet bases there are

-size discontinuities between representations on neighboring intervals (where

is user-selected) and repeated applications of derivative operators may lead to large errors and the appearance ofspiky artifactsatboundariesofintervals.Whileinmanyinstancesdifferentialequationscanbereplacedby integralequationsasin[18,20,45,46],theevaluationofhighorderderivativesofsolutionsisanecessarystepincomputing functionals that carry physical meaning.An interesting observationwas made in [13]:while the basis functionsare not smooth,givenmultiwaveletcoefficientsonecanlookforafunctionwiththemaximalnumberofderivativesthathasthese coefficients.Inthatapproach,coefficientsremainunchangedandthefunctioninterpretedashavingderivativesinbasesdual tomultiwavelets.

Themainideasofourapproachareasfollows.

1. We assume that the input coefficients are (slightly) corrupted and, for thisreason, first projectthe function onto a smooth basis and then use this projection to differentiate the original function. If we were to take the projection andcompute itsmultiwavelet coefficients, they willslightlydiffer fromtheoriginal ones. This approachpermits the constructionoffirst,second,andhigherderivatives.

2. Twosmooth intermediatebasesare explored.Whilewe cansafelyassumethat splinesare wellknown,we introduce exponential bases for bandlimited functions in greater detail since such bases are less known. We use them as an alternative to projections on splines to emphasize the generic nature of our approach and, also, due to remarkable propertiesoftheresultingoperators.

3. Theconstructionforallbasesretainsthe3-boxstenciloftheoriginalderivativeoperatorandforb-splinesalsoretains thesameorderofapproximation,andhencethenewoperatorscanbeusedasadrop-inreplacementinthesoftware.

Forbandlimited exponentials, theconceptof orderofapproximationis maintainedto theaccuracyofthe underlying quadratureandtheassociatedbandlimit.

Ournumericalresults forup tothird-order derivativesdemonstratethat thenewoperators displaythe expectedorderof approximation andare substantially superiorin controllingboth pointwiseandnormwise errorsin derivatives,especially withnoisyinputfunctions.Alsodemonstratedisimprovedconsistencybetweenthederivativeoperatorsandtheirpseudo- inverse (i.e., convolution with the corresponding Green’s function). We will see that the original operator retains some advantages,notablyincomputingthequantummechanicalkineticenergytohighaccuracy.

Motivating our work are applications inthe structure of nuclear matter [34,14,31,32] andmolecular electronic struc- ture[43,47,45,46,48,9,16,20] composedwithinthe MADNESSframework[21] andthe independentMRChempackage [25].

Both disciplines employ theKohn-Sham formulationof densityfunctionaltheory (DFT)[30] thatexpresses theenergyas a functionalof the density that is inturn parameterized by a setof orthonormal one-particle orbitals (i.e.,one-particle wave functions).Variation ofthe orbitalsto minimize the energy whileenforcing the orthonormalityconstraintyields a non-linear eigenvalueproblem(asetofcoupleddifferentialequations) thatforaccurate andefficient computationwithin themultiwaveletbasisisrecastintoasetofcoupledintegralequations.Inthemolecularproblemtheorbitalscanusually be chosen asreal-valued functionsof the positionof an electron in3D space, andthe occupied(i.e., bound ornegative energy) orbitalswillasymptotically decayexponentiallytosatisfy free-spaceboundaryconditions.Inthe nuclearproblem theorbitalswilltypicallybe2-componentspinorswitheachcomponentbeingacomplex-valuedfunctionofthepositionof a nucleon,andbothoccupiedandunoccupied(unboundorpositiveenergy)statesmustbecomputedforenergyfunction- als that incorporatesuperfluidity. Theunbound statesare notexponentially localizedand other techniquesare necessary to map the computation onto a finite volume [32]. While the DFT energy functionals are quite different for nucleiand electrons, they havemultiple elements in commonincluding the incorporation offirst, second and higherderivatives of the electron density aswell asstrongly non-linear functionsof these. Properties beyondthose of the groundstate (e.g., excitation energies or theresponse to external perturbations such astime-dependent electric ormagneticfields) require even higher derivatives.The amplification of numericalnoise inherent in computinghigher derivativeswith theoriginal scheme, thesensitivityof some non-linearfunctionals to point-wiseratherthan norm-wise errors,andthe sensitivityof some functionalstoregionswherethedensityissmall(andhencepotentiallydominatedbynumericalnoise),combineto

make computationswith functionals that incorporatederivatives moreexpensive comparedto those that do not include derivatives.Thisislargelyduetotheneedtocompute muchmoreaccurately(i.e.,oversample)inordertoobtainaccurate derivativesandfunctiontails.

Anotherissuearisingfromthenormoftheoriginalderivativeoperatorandtheneedtocomputehigherderivativesbyits repeatedapplicationisthatdifferentialoperatorsandtheirpseudo-inverse(i.e.,convolutionwiththecorrespondingGreen’s function)donotnumericallycommutewithintheprecisionofcomputation—i.e.,theorderofoperationsmattersandthere isa lackofconsistency betweenthenumericalandformal calculus. Thishasemerged asasignificant issueforscientists makingfirstuseofMADNESS,whichtypicallyrequiresreformulatingpartialdifferentialequationsasintegralequationswith free-spaceorperiodicboundaryconditions.Nominallyequivalentexpressionsmayhaveverydifferenterrorpropertiesthat can affectboth accuracy andconvergence.Forinstance, inthe electronic structure problemitis appealingto compute a residualfromthestandard eigenvalueproblemin differentialformwhileupdatingthe orbitalsusingthe Green’sfunction tothebound-stateHelmholtz(BSH) equation(whichbuildsindesirablephysicsinadditiontobeingnumericallyefficient;

[19]).However,thisfailstoconvergerobustlytohighaccuracysincetheavailableLaplacianoperatorstronglyamplifieshigh frequency noise that is suppressed by the integral operator. In the results section we will seethat the newdifferential operatorssignificantlyimproveuponpriorbehavior.

Wenotethatthereisalargenumberofpublicationsonconstructingderivativeoperators(intheformofafinitediffer- encescheme) withadditionaldesiredproperties- forexampleTVD[22,36,23] orWENO[24,29].However,inthiscasethe difficulty addressedishowto differentiateaccuratelyandavoidartifactsina schemewhereafunction developsashock, arapidvariation overasmallinterval. Notethatwe areaddressing an“opposite”problemwherefunctionsareknownto be differentiablemanytimesbutare representedviafunctionsthat havediscontinuities. Inmolecularelectronic structure thesolutionsaresmooth exceptforcuspsattheknownlocationsofthenuclearcentersifpointnuclearchargesareused, howeveritisroutinetoeliminatethesecuspseitheranalyticallyorthroughtheuseofeffectivepotentials.Anotherissuein constructingderivative operators addressedin[7,35,26] is howtoavoidspurious eigenvaluesofsuch operators(thatnec- essarilyoccurwhenusingpolynomialbases),whichnegativelyaffecttheoperatornorm.Weobtainasimilarresultinour constructionwhenusingbandlimitedexponentials.

Thisreturnsustothebasicconceptofthispaper—themultiwaveletbasisisveryconvenientforefficientlycomputing withlocally-refinedrepresentationsoffunctionsandoperators,andwecancircumventsomeofthedeficienciesofthemul- tiwaveletbasisby constructingoperator representationsinanunderlying basisselectedto providethedesiredproperties.

Webrieflyintroducethemultiwaveletbasis,review theconstructionandformoftheoriginalderivativeoperator,andthen introduceournewconstructions.Subsequently,wepresentanddiscussnumericalresultsandconclusions.

2. Themultiwaveletbasis

Themultiwaveletscalingfunctions[1] aregivenbyscaledandnormalizedLegendrepolynomialsdefinedon[0,1],

φ

i

(

x

) = √

2i

+

^1Pi

(

2x

−

¹

),

x

∈

^[0

,

1]

0

,

x

∈ / (

0

,

1

) ,

(1)

wherei=⁰,. . . ,k−^{1,kistheorderofthebasis(i.e.,thefirstk}polynomialsareused), Pi(x)aretheLegendrepolynomials, and√

2i+^{1 is a}normalizationfactor.From thesethescalingfunctionsinthel^th boxatlevelnareformedby translation anddilation

φ

ⁿ_il

(

x

) =

^2n/2

φ

i

(

2ⁿx

−

^l

)

(2)

Byconstructionthesefunctionsformanorthonormalbasisatleveln.

Alpert’sconstruction[1] ofthewaveletbasisstartswiththetwo-scalerelationshipbetweenpolynomialsatlevelnand n−¹

φ

ⁿ_i⁻¹

(

x

) = √

2

k−¹

j=⁰

(

h⁽_{i j}⁰⁾

φ

ⁿ_j

(

2x

) +

^h(_{i j}¹⁾

φ

ⁿ_j

(

2x

−

¹

),

i

=

⁰

, . . . ,

k

−

¹

,

(3)

whereH⁽⁰⁾= h⁽_{i j}⁰⁾

andH⁽¹⁾= h⁽_{i j}¹⁾

arethefiltercoefficients,whicharereadilycomputedbyorthogonalprojection.The orthonormalwaveletsatleveln−^{1 aresimilarlyexpandedintermsofthescalingfunctionsatleveln}

ψ

_iⁿ⁻¹

(

x

) = √

2

k−1

j=⁰

(

g⁽_{i j}⁰⁾

φ

ⁿ_j

(

2x

) +

^g_{i j}⁽¹⁾

φ

ⁿ_j

(

2x

−

¹

),

i

=

⁰

, . . . ,

k

−

¹

,

(4)

wherethefiltercoefficientsG⁽⁰⁾= g⁽_{i j}⁰⁾

andG⁽¹⁾= g⁽_{i j}¹⁾

are againcomputedby orthogonalprojection notingthatthe waveletsareorthogonaltothescalingfunctionsatthesameandcoarserlevels,andthatthewaveletsareorthogonaltoeach

other across allscales. Together,theserelations enablethe transformationbetweendifferentscalesofthe multiresolution analysis.Notethatattheirmid-pointthemultiwavelets(4) areeitherdiscontinuousorhaveadiscontinuousderivative.

Theprojectionofafunction f(x)intotheorder-kscalingfunctionbasisatlevelnis

fⁿ

(

x

) =

2

ⁿ−¹ l=⁰

k−¹

j=0

sⁿ_jl

φ

ⁿ_jl

(

x

),

(5)

wheresarethescalingfunctioncoefficients.Assuming f isk-timesdifferentiable,theapproximationerroris[1]

f

−

^fn

≤

^{2−nk 2} 4^kk

!

x∈[^sup0,1]

d^k dx^k f

(

x

)

.

(6) Repeated applicationofthe filters H⁽⁰⁾, H⁽¹⁾, G⁽⁰⁾ andG⁽¹⁾ (i.e., thefastwavelettransform [4])performs the orthogonal transformationintothecorrespondingwaveletrepresentation

fⁿ

(

x

) =

k−1

j=⁰

⎛

⎝

s⁰_j0

φ

j

(

x

) +

n−1

m=⁰ 2

^m−1

l=⁰

d^m_jl

ψ

^m_jl

(

x

)

⎞

⎠ .

Adaptive refinement isaccomplished by discarding small wavelet coefficientsaccording to user-specifiedcriteria. The waveletcoefficientswillbesmallbelowtherefinementlevelatwhichthefunctionlocallybecomeswellapproximatedby the Legendre basis. Thus,the noise orerrorin numericalcomputation has aknown form—it is dominated,forsmooth functions,bywaveletsatthenextfinestscale.

2.1. Derivativeoperator

As discussedinthe introduction,there isjust onescale-consistent derivativeoperator inthisbasis,whichmaybe ob- tainedeitherbyseekingtheuniqueoperatorthatishomogeneousofdegreeone,orbyprojectionandintegrationbyparts asistraditionallydonetoobtaintheweakderivative[3].Freeparametersintheoperatorpermitincorporationofboundary conditionsandconstructionofforward,central,orbackwardforms.Inthispaperwejustconsiderthecentral(orperiodic) form.

Giventheprojection(5) ofafunctionatleveln,thescalingfunctioncoefficients(¯sⁿ_lj)ofthederivativeinboxlatleveln arerelatedtotheinputcoefficientsfromthecorrespondingboxanditsneighbors(l−^{1 andl}+^{1)bythetransitionmatrices} r₋1,r0,andr1 ofa3-blockstencil

¯

sⁿ_jl

=

k−¹

i=0

[r1]_jisⁿ_il−1

+

[r0]_jisⁿ_il

+

[r₋1]_jisⁿ_il+1

,

(7)

withr₋1= −^rt1.Forthisunique,scale-consistentderivativeoperator,wealsohavethatthecentralblockisskew-symmetric, i.e., r₀= −^rt₀ ^{so that the entire operator is} skew-symmetric. This property will not hold for the smoothed derivatives constructedbelow.

3. Constructionofsmoothingdifferentialoperators

The scale-consistent construction ofthe derivative operator [3] proceeds withthe perspective that thefunction being differentiated is the expansion in the discontinuous scaling function (or equivalently multiwavelet) basis. Here, we in- stead changethis perspective and assume that the expansion inscaling functionsis an approximation to some at least p-differentiablefunctionaboutwhichweknowonlyapproximationstothefirstfewgeneralizedmomentswithineachbox

—i.e.,theexpansioncoefficientsareinterpretedastheprojectionofthissmoothfunctionintothelocalLegendrepolynomial basis oforderk,withthesecoefficientsonlybeingknowntotheprecision ofcomputation. Weconstructlocalrepresenta- tions ofthissmoothfunctionby projectingintoasmooth basisthat isatleast p-differentiableandwhich isdefinedover the domainof thederivative operator(a 3-boxstencil).The projection isperformedin aleast-squares sense soasto re- produceasaccuratelyaspossibletheknownmoments(i.e.,expansioncoefficientswithrespecttoscalingfunctions).With this smoothrepresentation inhand,the derivativeis computedandexpansion inthescaling function basisrecovered by orthogonal projection.Inpracticeallthreeoperations(projectionintothesmoothbasis,differentiation, projectionintothe scalingfunctionbasis)arecombinedintoasinglematrixrepresentation.

Byapplying thesmoothinglocallyatthesamescale atwhichthescalingfunction expansionis known,weaccommo- date adaptiverefinement(thatcan introducenoise atanyscale)anddonotenforce aglobalband limit thatwouldlimit applicability.Togetherwithincreasingtheaccuracy/band-limit/orderofthesmoothbasisintandemwiththatofthescaling functionbasis,wepreservetheabilitytocontroltheoverallaccuracyofcomputationthrougheitheradaptiverefinementor useofahigher-orderbasis(largerk). Therefore,nonewparameters areintroducedthat controltheerrorincomputation.

Theconstructionpreservesthe3-boxstencilandscale-invariant operator,andhencecanbeusedasadrop-inreplacement inexistingsoftware(e.g.,MADNESS[21],MRChem[25]).

We now explain our construction in a generalsetting that combines the projection anddifferentiation operators. Let usassumethatwehavetwobases,theoriginalbasisu1(x) ,u2(x) ,. . .un(x)andanauxiliary basisv1(x) ,v2(x) ,. . .vm(x), where m≤^{n and the latter basis spans a subspace of the former (exactly or} approximately within a given accuracy).

Specifically,we are interestedinthe casewhere thebasisfunctions v1,v2,. . .vm allhavederivativesoforder p whereas not all of thefunctions u1(x) ,u2(x) ,. . .un(x)are p-differentiable.As we alreadypointedout, our main exampleis the basiswhereu₁(x) ,u₂(x) ,. . .u_n(x)arepiece-wisepolynomials.

We assume that thefunction f hasthederivative of order p that we wantto evaluate butiswritten inthe original basis,

f

(

x

) =

n k=¹

f_ku_k

(

x

) +

^O

( ε ),

(8)

where

ε

istheaccuracyofcomputation.Notethatifthecoefficients f_kareknownexactlythentheerrorwouldbeentirely inthewaveletspacesatthecurrentandfinerscales,butinpracticethecoefficientswouldonlybeknownapproximately.

Letusstartbydefiningthen×^{mmatrixAwithentries}

α

il

=

u_i

(

x

)

v_l

(

x

)

dx

,

(9)

and,iftheoriginalbasisisnotorthonormal,alsothen×^{nGrammatrixGwithentries} G_ik

=

u_i

(

x

)

u_k

(

x

)

dx

,

i

,

k

=

¹

, . . .

n

.

Ourfirststepistoproject f onthesubspacespannedbyfunctionsv1(x) ,v2(x) ,. . .vm(x)reproducinginaleastsquares sensetheprojectionontotheoriginalbasis.Weseekthecoefficientsg_l sothat

n k=1

f_ku_k

(

x

) ≈

m

l=1

g_lv_l

(

x

).

(10)

Denotingf=(f1,f2, . . .fn),g=(g1,g2, . . .gm)andprojectingfromtheleftwithui(x),weseekgsuchthat A^∗Gf

=

^A∗Ag

.

Thus,weobtaincoefficientsof f intheauxiliarybasisas g

=

A^∗A

₋1

A^∗Gf andwrite

f

(

x

) =

m l=1

g_lv_l

(

x

).

Next,weapplythederivativeoforderptothefunction f andobtain d^p

dx^p f

(

x

) =

m l=¹

g_l d^p dx^pv_l

(

x

) .

Inordertoreturntoarepresentationintheoriginalbasis,weprojectthederivative _dx^dppv_l(x)onfunctionsu_i,

β

_il^(p)

=

u_i

(

x

)

^d

p

dx^pv_l

(

x

)

dx

,

i

=

¹

,

2

, . . . ,

n

,

l

=

¹

,

2

, . . .

m

.

Denotingby B^(p)then×^{m matrixwithentries}β_il^(p),weobtain thecoefficientsofthe pthderivativeoffunction f inthe originalbasisas

f^(p)

=

^G−1B(p)g

or

f^(p)

=

^G−1B(p)

A^∗A

₋1

A^∗Gf anddefinethedifferentiationmatrixas

D^(p)

=

^G−1B(p)

A^∗A

₋₁

A^∗G

.

Iftheoriginalbasisu₁,u₂,. . .u_n isorthonormal,thenGistheidentitymatrixandwehave D^(p)

=

^B(p)

A^∗A

₋₁

A^∗

.

(11)

Inourcasethefirstbasisconsistsofthemultiwaveletscalingfunctions(a.k.a.thenormalizedLegendrepolynomialsor, alternatively, the normalizedLagrange interpolating polynomials) in the three adjacentintervals andthe second basis is eitherthatofsplinesorofband-limitedfunctionsonthecombinedthree-boxintervals.

3.1. DerivativeoperatorsusingB-splines

As an auxiliary basis, a b-splinebasis functions oforder M are formed frompiecewise polynomialsof degree M−¹ andspans M+^{1 knots.}Derivatives up toorder M−^{2 arecontinuous acrossthe domain.Weemploy N knotsuniformly} distributed over [−¹,2] ^{(i.e., the union of the three boxes} [−¹,0]^, [⁰,1] ^and [¹,2]^{), with the location of the ith knot} (i=⁰,. . . ,N−^{1)givenby x}i= _N³ⁱ₋₁−^{1.Using theseknots weconstructthefullsplinebasis oforderM,repeating knots} as necessaryatthe endpoints. Associatedwitheach endpoint thereare M−^{1 splines withrepeatedknots,andin the} interior thereare N−^{M splines withoutrepeatedknots. Thus, thereis a totalof N}+^M−2 functions inthe basis. The b-splines withsupportincludingthecenterregion[⁰,1] ^{correspondtothecardinalsplines, whichhaveFouriertransform} proportionalto(sinc(

ω

))^M−1.Itisthispropertythatprovidestheprimarysourceofsmoothing.

Choosingtheb-splinebasistohavetheminimumordernecessarysuchthatthepth derivativeresultsinsplinesoforder k,weidentifyM=^k+^{p.ThischoiceisconsistentwithprojectionofthederivativebackintotheLegendrepolynomialbasis} oforderk,andis alsotheminimumordernecessarytoobtainnon-zeroapproximationstothefirst orsecond derivatives when computing with the Haar (piecewise constant) basis. Furthermore, Theorem 2 of [42] establishes that, with this choice of M,there exists a b-spline approximation to an appropriately smooth function f(x) such that the errorin the

pth derivativeoftheapproximationis O

c(k,p)2^−nk

,wherenisthelevelofrefinementinthedyadically-refinedadaptive meshsotheboxsizeatlevelnis O

2⁻ⁿ

,andc(k,p)issomeconstantindependentofn.Thiserror,whichisnumerically verifiedintheresultssection,isofthesameorderasthatoftheprojectionofasmoothfunctionintotheorder-k Legendre basis (6) andwiththatoftheoriginalderivativeoperator[3].Remarkably,the pth derivative(p≤^{2k)isexactiftheinput} function f(x)isexactlyapolynomial ofdegreek+^p−^{1 orlessover theentiredomain} [−¹,2] ^despitetheintermediate projection intothediscontinuous scaling functionbasis thatincludesonly polynomialsup todegree k−^{1 withineach of} the sub-intervals[−¹,0]^,[⁰,1]^and[¹,2]^{.Thisisduetothelinear}least-squaresreproductionofthemoments(expansion coefficientsinthescalingfunctionbasis)havingtheexactexpansionofthepolynomialintheb-splinebasis asitsunique, zero-residualsolution.

Inthisworkwechosethenumberofknots N=^k+1 independentof p.Theprimary motivationsforthisarethatitis againtheminimumnumbernecessarytoconstructfirstandsecondderivativesfortheHaarbasis,andforthefirstderivative (p=^{1)thischoiceyields2kdegreesoffreedominthesplinebasistobecomparedwith3kdegreesoffreedomintheinput} basis (i.e.,thefirstk Legendrepolynomialsineach ofthe3boxes).Thus,thereisadditionalsmoothingresultingfromthe fittingto thelower dimensionbasis, whichwas thoughttobe desirablefromtheperspective ofsuppressingnoise.Other choicesare clearlyfeasible, however,thischoice hasprovedsatisfactory andlimitedexperimentation didnotidentifyany clearlysuperiorchoices.

3.2. Derivativeoperatorsusingbandlimitedexponentials

Instead ofusingb-splines asanauxiliary basisintheconstruction ofsmoothingderivativeoperators, wecanuseban- dlimited exponentials (which are not necessarily periodic on an interval).The reasons to use these functionsinstead of piece-wise polynomialsare:first,itfurtherclarifiesourapproachtoconstructing“non-native”derivativeoperators inmul- tiwavelet bases and,second, the resulting operators haveremarkableproperties whichwe describe below. Asa practical matter,theaccuracyoftheresultsobtainedusingsuchderivativesissimilartothoseobtainedusingb-splines.

Polynomial approximationshave a longtradition in numericalanalysis andlead,interalia, tothe notion ofthe order of convergence of numericalschemes. For smooth periodic functions a well known alternative is to use approximations via trigonometricpolynomials(i.e.,viaatruncatedFourierseries).Forbandlimitednon-periodic functionsonanintervala naturalbasis isthatofeigenfunctionsofthe band-and-timelimitingoperatorintroducedina seriesofseminal papersby Slepian,LandauandPollak,[41,27,28,37–40],wherethey observedthattheeigenfunctionsofsuch operatoraretheProlate SpheroidalWaveFunctions (PSWFs)ofclassicalMathematicalPhysics.InsteadofusingPSWFsdirectly(which isapossible choice),weusebandlimitedexponentialsconstructedin[5,6] toprovideusanapproximateauxiliarybasis.

We observethat whilesmooth periodic band-limitedfunctionsmayberepresented viaa truncatedFourierseries,the Fourierseriesisnotefficientforband-limitedbutnon-periodicfunctionsonaninterval.Thismotivatesconsideringalinear spaceoffunctionsadmittingarepresentationviaexponentials

e^ibx

|b|≤c,x∈^[−¹,1], Ec

=

u

∈

^L∞

(

[

−

¹

,

1]

) |

^u

(

x

) =

k∈Z

a_ke^icbkx

: {

^ak

}

k∈Z

∈

^l1

,

b_k

∈

[

−

¹

,

1]

,

wherec(bandlimit)isafixedparameter(see[5,6]).Weseek(anapproximate)basisforthisspace.WhilePSWFscanserve thispurpose,itismoreconvenientinourcasetouseafixedsetofexponentials

e^icθkxM

k=1asdescribedinthesequel.

The use of band-limited functionson an interval for integration and interpolation relies on quadratures constructed in[44,5,33].¹ TheseGaussian-typequadraturesforbandlimitedexponentials differfromtheclassicalGaussian quadratures forpolynomialsinthattheyareapproximate.Indeed,withafinitenumberofnodesitisimpossibletointegrateexactlyan infinitenumberoffunctionsbut,itturnsout,thatallfunctionsEc canbeintegratedwithinanyfiniteuser-selectedaccuracy

.The generalizedGaussianquadratures forexponentialsto integratefunctionsinEc withaccuracy

are summarizedas (see[5,33,6])

Lemma.Forc>0andany >0,thereexistafinitenumberofnodes−¹< θ1< θ2<· · ·< θM<1andcorrespondingweights wk>0,suchthatforanyx∈^[−¹,1],

1

−1

e^ictxdt

−

M k=1

wke^icθkx

< ,

(12) wherethenumberofnodes,M=^c/

π

₊O(logc),is(nearly)optimal.Thenodesandweightsmaintainthenaturalsymmetry,θ_k=

−θM−^k+¹andwk=^wM−^k+¹.

Wenotethattheconstructionofquadraturesin[5,33] ismoregeneralthanformulatedhereandyieldsquadraturesfor band-limitedexponentialsintegratedwithaweightfunction.

Wenow describeinterpolation inEc.Givenafiniteaccuracy

,we seektorepresentfunctionsinEc bya fixedset of exponentials

e^icθkxM

k=1,whereMisassmallaspossible.Itturnsout(see[5])thatbyfindingquadraturenodes{θk}_k^M₌₁^and weights{^wk}^M_k=1 ^forexponentials withbandlimit2c andaccuracy ²,we infact obtainanapproximatebasistorepresent allexponentialsinEc withaccuracy .Asabasis,thesetofexponentials

e^icθkxM

k=¹ issimilartomonomials(aswellasto exponentialsintheFourierseries),whereasthePSWFsarestrikinglysimilartoorthogonalpolynomials.While,asageneral rule,itispreferabletoworkwithorthogonalbases,sinceweuseonlythefinal resultinconstructingderivativeoperators, we work directly with the basis

e^icθkxM

k=¹ anduse extended precision to address ill-conditioning. Once the necessary matricesare computedone can usethemin computationswiththe standarddouble precision. Whenbuilding derivative operators usingexponentials,it ismoreconvenientto firstsplittheinterval [−¹,1] intothreesubintervals,constructthe operators onthese subintervals,andthen rescalethe resulttointervals[−¹,0], [0,1] and [1,2] that havebeenused for b-splines.

Ourfirststepistoexpand thefunctions e^icθkxM

k=¹ on[−¹,1] into thebasisoftheLegendre polynomialsonintervals [−¹,−¹/3],[−¹/3,1/3] and[1/3,1].Wedenotetheorthonormalscalingfunctionsofamultiwaveletbasis(i.e.normalized Legendrepolynomials)onintervals[−¹,−¹/3],[−¹/3,1/3] and[1/3,1] as

φ⁻_i (x)i=^k−¹ i=⁰ ,

φ_i^c(x)i=k−1 i=⁰ and

φ⁺_i (x)i=^k−¹ i=⁰ , andcomputecoefficients

A⁻_im

=

−

1/3

−¹

e^icθmx

φ

_i⁻

(

x

)

dx

A^c_im

=

1/3

−1/3

e^icθmx

φ

_i^c

(

x

)

dx

,

and

A⁺_im

=

1 1/3

e^icθmx

φ

_i⁺

(

x

)

dx

,

where i=⁰,. . .k−^{1 and m}=¹,. . . ,M.Following (9), we then merge thesethree matricesto form3k×^{M matrix A}= {^Aim}i=1,...,3k

m=1,...,M,where3k≥^{M.Wechoosetherelationbetweenthenumberofnodes Mandtheorderofthebasisksothat} thebasis functions

e^icθkxM

k=1 areapproximatedby theLegendrepolynomialsofdegreesuptok−^{1 onthe}sub-intervals

1 Thesequadraturesallowustobreakwiththeconventionalapproachofusingpolynomialsandleadtoimprovementsinperformanceofalgorithms forinterpolation,estimationandsolvingpartialdifferentialequations(seee.g.[7,35,26,8]).Theimprovementshavetodowithbeingabletoapproachthe Nyquistrateforbandlimitedfunctionsonthereallinewhileworkingwithbandlimitedfunctionsonaninterval.

Fig. 1.Density plot of three blocks of the derivative matrix fork=16 with a resulting accuracy of approximately 10⁻⁷.

[−¹,−¹/3],[−¹/3,1/3] and[1/3,1] withaccuracy

,where

isthe accuracyofrepresentingexponentials inEc via the basisfunctions.

Givenabandlimited(orapproximatelybandlimited)functioninEc,

f

(

x

) =

M m=1

g_me^icθmx

,

(13)

we denotethevector ofits coefficientsasg= {^gm}_m^M₌₁^{.Using matrixA,we obtainthecoefficientsofthisfunction inthe} Legendrebasisontheunionofthreesubintervalsas

f

=

^Ag

.

(14)

Nextwecomputethederivativeof f in(13) as

f

(

x

) =

^ic

M m=1

g_m

θ

_me^icθmx

,

notingthatthegeneralizationtohigherorderderivativesisstraightforward.DenotingbySthediagonalmatrix,

S

=

^icdiag

(θ

1

, θ

2

, . . . θ

M

) ,

(15)

wecomputecoefficientsofthederivativeintheLegendrebasisas f

=

^ASg

.

Sincewewanttomapfdirectlyintof,weestimategfrom(14) byusingthenormalequationA^∗Ag=^A∗f, f

=

^AS

A^∗A

₋1

A^∗Ag

=

^AS

A^∗A

₋1

A^∗f

.

The3k×^{3kderivativematrixisthen} D

=

^AS

A^∗A

₋1

A^∗

,

(16)

where (A^∗A)⁻¹A^∗ isthepseudo-inverse ofA.All computationsmustbedone withtheextended precisionto accountfor ill-conditioningofAbut,oncecomputed,theresultingmatrixDcanbeusedinthestandardprecision.

Remarkably, while the resulting matrix D is not anti-symmetric, its eigenvalues are either pure imaginary or zeros.

Indeed,onits rangematrixDissimilartoS asitfollowsfrom(16). ExtractingthreemiddleblocksofthematrixD,each ofthemofsizek×^{k,re-scalingthemby thefactor 2}₃^{,we obtainthederivative operatorwithsmoothingbyprojection on} bandlimitedfunctionswiththebandlimit ²₃c.Thefactor ²₃ accountsforswitchingfromtheinterval[−¹,1] (oflength2)to theinterval[−¹,2] (oflength3)whichweusedinconstructingspline-basedprojections.

The sub-matricescomprisingthe resultingderivative operatorare illustratedinFig. 1.There isnoobserved symmetry betweenorwithintheblocksandtheoff-diagonalblocksarenotrankoneasinthecaseofthemultiwaveletderivative.

Similarly,thesecondderivativeisconstructedas D²

=

^AS

A^∗A

₋1

A^∗AS

A^∗A

₋1

A^∗

=

^AS2

A^∗A

₋1

A^∗

andweuseitsthreemiddleblocks.Again,we observethattheresulting3k×^{3kmatrixD2 isnotsymmetricbuthasonly} negativeorzeroeigenvalues.Higherorderderivativesareobtainedinasimilarmanner(notethatin(11)B^(p)=^ASp).

SinceonitsrangematrixDissimilartoS,thenormofthederivativeoperatorD^piscontrolledbythenormofS^p which doesnotexceedc^p,wherecisthebandlimit(15).Asin[7,35,26],wethusavoidspuriouseigenvaluesthatnecessarilyoccur

Fig. 2.Anexampleofabasisfunction,e^icθ5x,fromtheset e^icθkxM

k=1,whereM=25 displayedontheinterval[−¹,−¹/3].Suchfunctionsareapproximated byalinearcombinationoftheLegendrepolynomialsofdegreesi=⁰,. . .k−^{1,wherein1weusedk}=16 (yieldingaccuracy≈³·^10−7).

whenusingpolynomialbases.Thechoiceofthebandlimitc dependsontheorderofthemultiwaveletbasiswhichinturn dependsontheproblembeingsolvedandthedesiredaccuracy.

Asan example,we usebasis e^icθkxM

k=1,whereM=^{25,andthebandlimitc}=¹⁵.70796326794900 anddisplayone of thefunctions,e^icθ5x,ontheinterval[−¹,−¹/3] inFig.2.Asscalingfunctions,weusetheLegendrepolynomialsofdegrees i=0,. . .k−1,wherek=16,approximatingthebasis

e^icθkxM

k=¹ separately onthesub-intervals [−1,−1/3], [−1/3,1/3]

and[1/3,1] with accuracy ≈³·^{10−7.Asdescribed above,we createa 3k}×^{3kderivative matrixandthenextractthree} k×^{kblockstobeusedasthefirstderivativeoperator.Thesethree16}×16 matricesaredisplayedinFig.1.

4. Demonstrationanddiscussion 4.1. Point-wiseerrors

InFig.3,weexaminethenatureoflocalerrorsinthederivativesbyrepeatedapplicationofthefirstderivativeoperator toa3DGaussianandplottingalongthex-axistheabsoluteerrorfromtheexactresults.Byvaryingthelevelatwhichthe Gaussian is projected intothe scaling function basis we explore from low accuracy andundersampling (level n=^{4) to} highaccuracyandoversampling(n=^{7).Itisapparenthowtheoriginalderivativeoperatoramplifiesthe}discontinuitiesat theedges ofthesamplingboxes (thecomputation isinthecube[−¹⁶,16]^{3 sotheboxesatlevelnhavesize25−n),with} repeatedapplicationtoobtainhigherderivativesmagnifyingthiseffect.Refinementalsomagnifiesthisdefectintheoriginal operatorwithanextrafactorofabout2witheachlevelofrefinement.Toexplainthemainoriginofthiseffect,considerthe function(1+^erf(2ⁿx))

/2 whichisasmoothapproximationtoastepdiscontinuityofheight atscalen—itsderivativeat x=^{0 is2n}/

π

.Themaximumvalueofthederivativenecessarilyincreasesatfinerscalesinordertopreservethemagnitude of its integral. The other two constructions greatly suppressthe impact of these artifactsof the discontinuous basis. At low refinement(n=^{4),the original}first-derivative operator performs slightlybetterthan theb-spline andBLE operators

—presumablytheir filteringisdiscarding essentialinformationinthehigh-frequencies,andsince theBLEoperatorfilters moststrongly(seeFig.9)itperformstheleastwell.Withincreasingaccuracy/refinement(n=⁵,6,7),theoriginaloperator performsincreasinglyworserelativetotheb-splineoperatorthatisabletoattainhighprecision.Atintermediaterefinement (n=^{5), the b-spline andBLE operators behave similarly. However, at deeper refinement (n}=⁶,7), the accuracy of the quadratureusedintheconstructionoftheBLEoperator dominatestheerror(in thiscasearound 10⁻⁹)sinceaccuracy of computationishigherthanthatofthequadrature.Whileitispossibletoconstructahigheraccuracyquadrature,itisnot typicalthatwecomputetogreaterthan10significantdigits.

InFig.4,usingthesametest3Dtestfunction,fortheb-splineandBLEconstructionswecontrastthesizeoflocalerrors obtainedbyrepeatedapplicationofthefirstderivativeoperator(

d dx

p

for p=²,3)withthosefromdirectapplicationof thehigh-order derivative operator (_dx^dpp) constructedasdescribed inSection 3.Forthe secondderivative, repeatedversus direct application have fairly comparable errors. However, for the third derivative (available for the b-spline only) it is apparent that repeatedapplication atdeeper refinement (orhigher accuracy) hasa significant advantage over thedirect application.Both rulesshouldhavethesameorderaccuracy(byconstructioninSection3.1andnumericallydemonstrated inFigs.6and7),andweinsteadattributethistotherepeatedapplicationcorrespondingtoa7-boxstencilandratherthan justthe3-boxstencilofthedirectapplication.

4.2. Norm-wisepropertiesanderrors

Thespectralnormofthethreeoperators(largestsingularvalueoftherectangularmatrixrepresentation)wascomputed asafunctionofk(Fig.5).Asexpected,thenormoftheoriginaloperatorgrowsquadraticallywithk,whereastheothertwo

Fig. 3.Fortheoriginal,b-splineandBLEconstructions,theplotdisplaystheabsolutevalueoftheerroralongthex-axisofthefirst,secondandthird derivativesw.r.t.xcomputedwithrepeatedapplicationofthefirstderivativeoperatorappliedin3-dimensionstof(x,y,z)=_π

2

(3/2)

e^{−(x2+y2+z2)}projected intothek=^{10scalingfunctionbasisinthecube}[−¹⁶,16]^{3.Columnsone,twoandthreecorrespondtothefirst,secondandthird}derivatives,respectively.

Rowsonethroughfourcorrespondtothetestfunctionbeingprojectedatlevels4through7,respectively.

growonlylinearlywithanearunitslope.Furthermore,examinationoftherightsingularvectorassociatedwiththelargest singularvalueoftheoriginaloperatorrevealsastrongresemblancetoamultiwavelet.

Figs. 6and7demonstratethatweobtain thetargetorderofaccuracyandconvergence.InFig.6,we demonstratethat forfixedk=^{7 weobtainanerrorforallthreeb-spline}derivativesthatscalesas O

ω

^kforsmall

ω

whenappliedtoone ofsin(

ω

x)orcos(

ω

x).InFig.7,wedemonstratethatforfixed

ω

=^{4.642weobtainfor}applicationtosin(

ω

x)anerrorthat asymptoticallyscalesasO

h^k

forallthreeb-splinederivatives(includingtheerrorarisingfromprojection).

AllthreefirstderivativeoperatorsareexaminedinFigs.8and9withapplicationtosin(

ω

x)inthek=^{7 basis.InFig.8,} which shouldbecompared withFig. 6,we demonstratethatall threeoperators providetheexpectedscalingoftheerror with

ω

,andthattheBLEandb-splineoperatorsaremoreaccuratethantheoriginaloperatorsexceptathigh

ω

oratlow

ω

(where theaccuracyoftheBLEoperatorisconstrainedbythechosenaccuracyofthequadraturerule).Wecannotstrictly speak ofthefrequencyresponse orFouriertransformofthe operators becausethey are appliedon afinite,non-periodic domain,however,inFig.9weexaminethemagnificationfactor(essentiallythenormoftheresult)foreachoperatorapplied tosin(

ω

x)asafunctionof

ω

.Thelargenormoftheoriginaloperatorisimmediatelyapparent,asisthefilteringproperties of the other two operators, with the BLE construction filtering more strongly than the b-spline construction as already discussedabove.

4.3. Kineticenergyinquantummechanics

Innon-relativisticquantummechanics,thekineticenergyoperatorisessentiallyjustthenegativeoftheLaplacian.Thus, ifthewavefunctiondescribingthestateis f(x),theenergyisgivenby

−

f

(

x

)

^∗

∇

^2f

(

x

)

dx

=

(∇

^f

(

x

))

^†

.∇

^f

(

x

)

dx (17)

withtheequalityholdingforperiodicorfree-spaceboundaryconditions.Similarintegralsariseinothersettings,forinstance as the stiffness matrix in finite element analysis. The weak (second) form is especially interesting because its error is