Fast Adaptive Algorithms in the Non-Standard Form for Multidimensional Problems
aa r X i v : . [ m a t h . NA ] A ug Fast Adaptive Algorithms in the Non-StandardForm for Multidimensional ProblemsGregory Beylkin, Vani Cheruvu and Fernando PérezDepartment of Applied Mathemati s, University of Colorado, Boulder, CO80309-0526, United StatesAbstra tWe present a fast, adaptive multiresolution algorithm for applying integral operatorswith a wide lass of radially symmetri kernels in dimensions one, two and three.This algorithm is made e(cid:30) ient by the use of separated representations of the kernel.We dis uss operators of the lass ( − ∆ + µ I ) − α , where µ ≥ and < α < / ,and illustrate the algorithm for the Poisson and S hrödinger equations in dimensionthree. The same algorithm may be used for all operators with radially symmetri kernels approximated as a weighted sum of Gaussians, making it appli able a rossmultiple (cid:28)elds by reusing a single implementation.This fast algorithm provides ontrollable a ura y at a reasonable ost, ompa-rable to that of the Fast Multipole Method (FMM). It di(cid:27)ers from the FMM bythe type of approximation used to represent kernels and has an advantage of beingeasily extendable to higher dimensions.Key words: Separated representation; multiwavelets; adaptive algorithms; integraloperators.1 Introdu tionFor a number of years, the Fast Multipole Method (FMM) [1,2,3℄ has beenthe method of hoi e for applying integral operators to fun tions in dimen-sions d ≤ . On the other hand, multiresolution algorithms in wavelet andmultiwavelet bases introdu ed in [4℄ for the same purpose were not e(cid:30) ient This resear h was partially supported by DOE grant DE-FG02-03ER25583,DOE/ORNL grant 4000038129 and DARPA/ARO grant W911NF-04-1-0281.Preprint submitted to Elsevier O tober 25, 2018nough to be pra ti al in more than one dimension. Re ently, with the in-trodu tion of separated representations [5,6,7℄, pra ti al multiresolution algo-rithms in higher dimensions [8,9,10,11℄ be ame available as well. In this paperwe present a new fast, adaptive algorithm for applying a lass of integral op-erators with radial kernels in dimensions d = 1 , , , and we brie(cid:29)y dis uss itsextension to higher dimensions.In physi s, hemistry and other applied (cid:28)elds, many important problems maybe formulated using integral equations, typi ally involving Green's fun tionsas their kernels. Often su h formulations are preferable to those via partialdi(cid:27)erential equations (PDEs). For example, evaluating the integral expressingthe solution of the Poisson equation in free spa e (the onvolution of theGreen's fun tion with the mass or harge density) avoids issues asso iated withthe high ondition number of a PDE formulation. Integral operators appearin (cid:28)elds as diverse as ele trostati s, quantum hemistry, (cid:29)uid dynami s andgeodesy; in all su h appli ations fast and a urate methods for evaluatingoperators on fun tions are needed.The FMM and our approa h both employ approximate representations ofoperators to yield fast algorithms. The main di(cid:27)eren e lies in the type ofapproximations that are used. For example, for the Poisson kernel /r indimension d = 3 , the FMM [3℄ uses a plane wave approximation starting fromthe integral r = 12 π Z ∞ e − λ ( z − z ) Z π e iλ (( x − x ) cos α +( y − y ) sin α ) dαdλ, (1)where r = q ( x − x ) + ( y − y ) + ( z − z ) . The elegant approximation de-rived from this integral in [3℄ is valid within a solid angle, and thus requiressplitting the appli ation of an operator into dire tional regions; the numberof su h regions grows exponentially with dimension. For the same kernel, ourapproa h starts with the integral r = 2 √ π Z ∞−∞ e − r e s + s ds, (2)and its dis retization with (cid:28)nite a ura y ǫ yields a spheri ally symmetri approximation as a weighted sum of gaussians. Other radial kernels an besimilarly treated by a suitable hoi e of integrals. The result is a separatedrepresentation of kernels and, therefore, an immediate redu tion in the ostof their appli ation. This di(cid:27)eren e in the hoi e of approximation di tatesthe di(cid:27)eren es in the orresponding algorithms. In dimension d ≤ both ap-proa hes are pra ti al and yield omparable performan e. The key advantageof our approa h is its straightforward extensibility to higher dimensions [6,7℄.Given an arbitrary a ura y ǫ , we e(cid:27)e tively represent kernels by a set of ex-ponents and weights des ribing the terms of the gaussian approximation of2ntegrals like in (2). The number of terms in su h sum is roughly proportionalto log( ǫ − ) , or a low power of log( ǫ − ) , depending on the operator. Sin e op-erators are fully des ribed up to an a ura y ǫ by the exponents and weightsof the sum of gaussians, a single algorithm applies all su h operators. Thesein lude operators su h as ( − ∆ + µ I ) − α , where µ ≥ and < α < / ,and ertain singular operators su h as the proje tor on divergen e-free fun -tions. Sin e many physi ally signi(cid:28) ant operators depend only on the distan ebetween intera ting obje ts, our approa h is dire tly appli able to problemsinvolving a wide lass of operators with radial kernels.We ombine separated and multiresolution representations of kernels and usemultiwavelet bases [12℄ that provide inter alia a method for dis retizing inte-gral equations, as is the ase in quantum hemistry [8,9,11,10℄. This hoi e ofmultiresolution bases a ommodates integral and di(cid:27)erential operators as wellas a wide variety of boundary onditions, without degrading the order of themethod [13,14℄. Multiwavelet bases retain the key desired properties of waveletbases, su h as vanishing moments, orthogonality, and ompa t support. Due tothe vanishing moments, wide lasses of integro-di(cid:27)erential operators have ane(cid:27)e tively sparse matrix representation, i.e., they di(cid:27)er from a sparse matrixby an operator with small norm. Some of the basis fun tions of multiwaveletbases are dis ontinuous, similar to those of the Haar basis and in ontrast towavelets with regularity (see e.g. [15,16℄). The usual hoi es of s aling fun -tions for multiwavelet bases are either the s aled Legendre or interpolatingpolynomials. Sin e these are also used in the dis ontinuous Galerkin and dis- ontinuous spe tral elements methods, our approa h may also be seen as anadaptive extension of these methods.The algorithm for applying an operator to a fun tion starts with omputingits adaptive representation in a multiwavelet basis, resulting in a d -tree withblo ks of oe(cid:30) ients at the leaves. Then the algorithm adaptively applies the(modi(cid:28)ed) separated non-standard form [4℄ of the operator to the fun tion byusing only the ne essary blo ks as di tated by the fun tion's tree represen-tation. We note that in higher dimensions, d ≫ , fun tions need to be ina separated representation as well, sin e the usual onstru tions via bases orgrids are prohibitive (see [6,7℄).We start in Se tion 2 by re alling the basi notions of multiresolution analy-sis, non-standard operator form and adaptive representation of fun tions un-derlying our development. We then onsider the separated representation forradially symmetri kernels in Se tion 3, and use it to e(cid:30) iently extend themodi(cid:28)ed ns-form to multiple spatial dimensions in Se tion 4. We pay parti -ular attention to omputing the band stru ture of the operator based on onedimensional information. We use this onstru tion in Se tion 5 to introdu ethe adaptive algorithm for appli ation of multidimensional operators in themodi(cid:28)ed ns-form, and illustrate its performan e in Se tion 6. We onsider3wo examples: the Poisson equation in free spa e and the ground state of theHydrogen atom. We on lude with a brief dis ussion in Se tion 7.2 Preliminary onsiderationsThis se tion and Appendix are provided for the onvenien e of the reader inorder to keep this paper reasonably self- ontained. We provide ba kgroundmaterial and introdu e ne essary notation.The essen e of our approa h is to de ompose the operator using proje torson a Multiresolution Analysis (MRA), and to e(cid:30) iently apply its proje tionsusing a separated representation. We use the de omposition of the operatorinto the ns-form [4℄, but we organize it di(cid:27)erently (thus, modi(cid:28)ed ns-form) toa hieve greater e(cid:30) ien y. This modi(cid:28) ation be omes important as we extendthis algorithm to higher dimensions.In this se tion we introdu e notation for MRA, des ribe the adaptive repre-sentation of fun tions and asso iated data stru tures, introdu e the modi(cid:28)edns-form and an algorithm for its adaptive appli ation in dimension d = 1 asba kground material for the multidimensional ase.2.1 Multiresolution analysisLet us onsider the multiresolution analysis as a de omposition of L ([0 , d ) into a hain of subspa es V ⊂ V ⊂ V ⊂ · · · ⊂ V n ⊂ . . . , so that L ([0 , d ) = ∪ ∞ j =0 V j . We note that our indexing of subspa es (in reas-ing towards (cid:28)ner s ales) follows that in [13℄, and is the reverse of that in [4,15℄.On ea h subspa e V j , we use the tensor produ t basis of s aling fun tions ob-tained using the fun tions φ jkl ( x ) ( k = 0 , . . . , p − whi h we brie(cid:29)y des ribein Appendix.The wavelet subspa es W j are de(cid:28)ned as the orthogonal omplements of V j in V j − , thus V n = V ⊕ nj =0 W j . Introdu ing the orthogonal proje tor on V j , P j : L ([0 , d ) → V j and on-sidering an operator T : L ([0 , d ) → L ([0 , d ) , we de(cid:28)ne its proje tion T j : V j → V j as T j = P j TP j . We also onsider the orthogonal proje tor Q j : L ([0 , d ) → W j , de(cid:28)ned as Q j = P j +1 − P j .4.2 Adaptive representation of fun tionsLet us des ribe an adaptive re(cid:28)nement strategy for onstru tion multiresolu-tion representations of fun tions f : B → B , where B = [0 , d . We pro eedby re ursive binary subdivision of the box B , so the basi stru ture represent-ing our fun tions is a d − tree with arrays of oe(cid:30) ients stored at the leaves(terminal nodes) and no data stored on internal nodes. On ea h box obtainedvia this subdivision, our basis is a tensor produ t of orthogonal polynomials ofdegree k = 0 , . . . , p − in ea h variable, as des ribed in Appendix 8. Therefore,the leaves arry d -dimensional arrays of p d oe(cid:30) ients whi h may be used toapproximate fun tion values anywhere in the box orresponding to the spatialregion overed by it, via (30) or its equivalent for higher values of d . For on- iseness, we will often refer to these d -dimensional arrays of oe(cid:30) ients storedat tree nodes as fun tion blo ks.This adaptive fun tion de omposition algorithm is similar to that used in [17℄.Su h onstru tion formally works in any dimension d . However, sin e its om-plexity s ales exponentially with d , its pra ti al use is restri ted to fairly lowdimension, e.g. d . . In higher dimensions, alternate representation strate-gies for fun tions su h as [6,7℄ should be onsidered. In high dimensions, thesestrategies deal with the exponential growth of omplexity by using ontrolledapproximations that have linear ost in d .For simpli ity, we will des ribe the pro edure for the one-dimensional asesin e the extension to dimensions d = 2 , is straightforward. Sin e we annot a(cid:27)ord to onstru t our representation by starting from a (cid:28)ne s ale (es-pe ially in d = 3 ), we pro eed by su essive re(cid:28)nements of an initial oarsesampling. This approa h may result in a situation where the initial samplingis insu(cid:30) ient to resolve a rapid hange in a small volume; however, in pra ti alappli ations su h situations are rare and may be avoided by an appropriate hoi e of the initial sampling s ale.Let B jl = [2 − j l, − j ( l + 1)] , l = 0 , . . . , j − , represent a binary subinterval ons ale j . We denote by f jl = n f jl ( x k ) o p − k =0 the ve tor of values of the fun tion f on the Gaussian nodes in B jl . From these values we ompute the oe(cid:30)- ients n s jkl o p − k =0 (see (32) in Appendix) and interpolate f ( x ) for any x ∈ B jl by using (30). We then subdivide B jl into two hild intervals, B j +12 l and B j +12 l +1 ,and evaluate the fun tion f on the Gaussian nodes in B j +12 l and B j +12 l +1 . Wethen interpolate f by using the oe(cid:30) ients n s jkl o p − k =0 from their parent inter-val and denote by e f j +12 l and e f j +12 l +1 the ve tors of interpolated values on thetwo subintervals. Now, if for a given toleran e ǫ either (cid:13)(cid:13)(cid:13) f j +12 l − e f j +12 l (cid:13)(cid:13)(cid:13) > ǫ or (cid:13)(cid:13)(cid:13) f j +12 l +1 − e f j +12 l +1 (cid:13)(cid:13)(cid:13) > ǫ , we repeat the pro ess re ursively for both subintervals,5 j +12 l and B j +12 l +1 ; otherwise, we keep the oe(cid:30) ients n s jkl o p − k =0 to represent thefun tion on the entire interval B jl . This interval then be omes a leaf in ourtree.At this stage we use the ℓ ∞ norm, thus onstru ting an approximation ˜ f tothe original fun tion f su h that (cid:13)(cid:13)(cid:13) f − e f (cid:13)(cid:13)(cid:13) ∞ < ǫ , whi h immediately impliesthat (cid:13)(cid:13)(cid:13) f − e f (cid:13)(cid:13)(cid:13) < ǫ . This estimate learly extends to any dimension. On ethe approximation with ℓ ∞ norm is onstru ted, the orresponding tree maybe pruned if an appli ation only requires the approximation to be valid inthe ℓ norm. We start this pro ess on the (cid:28)nest s ale and simply remove allblo ks whose umulative ontribution is below ǫ . Other norms, su h as H , an be a ommodated by appropriately weighing the error toleran e with as ale-dependent fa tor in the initial ( oarse to (cid:28)ne) de omposition pro ess.The omplete de omposition algorithm pro eeds by following the above re ipe,starting with an initial oarse s ale (typi ally j = 0 ) and ontinuing re ursivelyuntil the stopping riterion is met for all subintervals. In pra ti e, we hoose astopping s ale j max, beyond whi h the algorithm will not attempt to subdivideany further. Rea hing j max means that the fun tion has signi(cid:28) ant variationswhi h are not a urately resolved over an interval of width − j max using a basisof order p . A pseudo- ode listing of this pro ess is presented as Algorithm 1.Algorithm 1 Adaptive Fun tion De omposition.Start at a oarse s ale, typi ally j = 0 .Re ursively, for all boxes b j on s ale j , pro eed as follows:Constru t the list C of d hild boxes b j +1 on s ale j + 1 .Compute the values of the fun tion f ( b j ) at the p d Gauss-Legendre quadra-ture nodes in b j .for all boxes b j +1 ∈ C doFrom f ( b j ) , interpolate to the Gauss-Legendre quadrature nodes in b j +1 ,produ ing values ˜ f ( b j +1 ) .Compute the values of f ( b j +1 ) at the Gauss-Legendre nodes of b j +1 , bydire t evaluation.if (cid:13)(cid:13)(cid:13) f ( b j +1 ) − ˜ f ( b j +1 ) (cid:13)(cid:13)(cid:13) ∞ > ǫ thenRe ursively repeat the entire pro ess for all boxes b j +1 ∈ C .end ifend for b j in the fun tiontree. 6 . . . . . . − . − . . . . (a) Fun tion. . . . . . . − − − − − − − − − − | Error | Cutoff (b) Pointwise error. N nod = 8 ,ǫ = 1 . e − ,N blocks = 17 ( ) Adaptive tree N nod = 8 ,ǫ = 1 . e − ,N blocks = 33 (d) Redundant treeFigure 1. The fun tion f ( x ) = sin(16 πx ) shown in (a), is de omposed with p = 8 and ǫ = 10 − ; (b) shows the pointwise approximation error. ( ) is the resulting adap-tive tree, where smaller subdivisions are required in regions with higher frequen y ontent. (d) is the redundant tree asso iated with this adaptive de omposition, whereall internal nodes have been (cid:28)lled with data.2.2.1 Tree stru tures for representing fun tionsThe de omposition Algorithm 1 naturally leads to a tree data stru ture torepresent fun tions, with the leaves of the tree orresponding to the spatialintervals over whi h the multiwavelet basis provides a su(cid:30) iently a uratelo al approximation. By using (30) or its higher-dimensional extensions, theonly data needed to approximate f ( x ) anywhere in B is the array of basis oe(cid:30) ients on these leaves. Thus we use a tree stru ture where the leavesstore these oe(cid:30) ients and the internal nodes do not ontain any data (andare e(cid:27)e tively removed sin e we use hash tables for storage). We will referto this stru ture as an adaptive tree. Ea h level in the tree orresponds to as ale in the MRA, with the root node orresponding to the oarsest proje tion f ∈ V .Now, in order to apply the modi(cid:28)ed non-standard form of an operator to a7un tion, we will show in the next se tion that we also need the basis oe(cid:30)- ients orresponding to internal nodes of the tree. Hen e, from the adaptivetree data stru ture we will ompute a similar tree but where we do keep the oe(cid:30) ients of s aling fun tions on all nodes (leaves and internal).The oe(cid:30) ients on the internal nodes are redundant sin e they are omputedfrom the fun tion blo ks stored in the leaves. We will thus refer to the tree ontaining oe(cid:30) ients on all nodes as a redundant tree. It is onstru ted re- ursively starting from the leaves, by proje ting the s aling oe(cid:30) ients fromall sibling nodes onto their parent node, using the de omposition (41).Figure 1 shows both the adaptive and the redundant trees for a sample fun -tion. This (cid:28)gure displays the oarsest s ales at the bottom and progressively(cid:28)ner ones further up, with (cid:28)lled boxes representing nodes where s aling o-e(cid:30) ients are stored and empty boxes indi ating nodes with no data in them(these do not need to be a tually stored in the implementation).2.3 Modi(cid:28)ed ns-formThe non-standard form [4℄ (see also [13℄ for the version spe ialized for multi-wavelets) of the operator T is the olle tion of omponents of the teles opi expansion T n = ( T n − T n − ) + ( T n − − T n − ) + · · · + T = T + n − X j =0 ( A j + B j + C j ) , (3)where A j = Q j TQ j , B j = Q j TP j , and C j = P j TQ j . The main propertyof this expansion is that the rate of de ay of the matrix elements of the op-erators A j , B j and C j away from the diagonal is ontrolled by the numberof vanishing moments of the basis and, for a (cid:28)nite but arbitrary a ura y ǫ ,the matrix elements outside a ertain band an be set to zero resulting in anerror of the norm less than ǫ . Su h behavior of the matrix elements be omes lear if we observe that the derivatives of kernels of Calderón-Zygmund andpseudo-di(cid:27)erential operators de ay faster than the kernel itself. If we use theTaylor expansion of the kernel to estimate the matrix elements away from thediagonal, then the size of these elements is ontrolled by a high derivative ofthe kernel sin e the vanishing moments eliminate the lower order terms [4℄.We note that for periodi kernels the band is measured as a periodi distan efrom the diagonal, resulting in (cid:28)lled-in ` orners' of a matrix representation.Let us introdu e notation to show how the teles opi expansion (3) is usedwhen applying an operator to a fun tion. If we apply the proje tion of theoperator T j − not on its (cid:16)natural s ale(cid:17) j − , but on the (cid:28)ner s ale j , wedenote its upsampled version as ↑ ( T j − ) . In the matrix representation of8 j − , this operation results in the doubling of the matrix size in ea h dire tion.This upsampling ↑ ( · ) and downsampling ↓ ( · ) notation will also be used forproje tions of fun tions.With this notation, omputing g = T f via (3) splits a ross s ales, ˆ g = T f ˆ g = [ T − ↑ ( T )] f ˆ g = [ T − ↑ ( T )] f (4) . . . . . . . . . ˆ g j = [ T j − ↑ ( T j − )] f j . . . . . . . . . where f j = P j f .As in the appli ation of the usual ns-form in [4℄, to obtain g n after buildingthe set { ˆ g , ˆ g , . . . , ˆ g n } , we have to ompute g n = ˆ g n + ↑ (ˆ g n − + ↑ (ˆ g n − + ↑ (ˆ g n − + . . . + ( ↑ ˆ g ) . . . ))) . (5)The order of the parentheses in this expression is essential, as it indi ates theorder of the a tual operations whi h are performed starting on the oarsestsubspa e V . For example, if the number of s ales n = 4 , then (5) yields g = ˆ g + ↑ (ˆ g + ↑ (ˆ g + ↑ (ˆ g + ( ↑ ˆ g )))) , des ribing the sequen e of ne essaryoperations.Unfortunately, the sparsity of the non-standard form indu ed by the vanishingmoments of bases is not su(cid:30) ient for fast pra ti al algorithms in dimensionsother than d = 1 . For algorithms in higher dimensions, we need an additionalstru ture for the remaining non-zero oe(cid:30) ients of the representation. We willuse separated representations (see Se tion 3) introdu ed in [6,18℄ and (cid:28)rstapplied in a multiresolution setting in [8,9,10,11℄. Within the retained bands,the omponents of the non-standard form are stored and applied in a separatedrepresentation and, as a result, the numeri al appli ation of operators be omese(cid:30) ient in higher dimensions.2.4 Modi(cid:28)ed ns-form in 1DLet us des ribe a one-dimensional onstru tion for operators on L ([0 , tointrodu e all the features ne essary for a multidimensional algorithm. Sin ewe use banded versions of operators, we need to introdu e the ne essary book-keeping. 9he (cid:16)template(cid:17) for the band stru ture on s ale j omes from the band on theprevious s ale j − . For ea h blo k on s ale j − , the upsampling operation ↑ ( T j − ) reates four blo ks (all ombinations of even/odd row and olumnindi es). We insist on maintaining the stri t orresponden e between these fourblo ks and those of T j . For this reason the des ription of the retained blo ksof T j involves the parity of their row and olumn indi es. Let us denote theblo ks in the matrix representing T j by t j ; ll ′ , where l, l ′ = 0 . . . j − . Individualelements within these blo ks are indexed as t j ; ll ′ ii ′ , where i, i ′ = 0 , . . . , p − , and p is the order of the multiwavelet basis. For a given width of the band b j , wekeep the operator blo ks t j ; ll ′ with indi es satisfying l − b j + 1 ≤ l ′ ≤ l + b j , for even l , l − b j ≤ l ′ ≤ l + b j − , for odd l. (6)We denote the banded operators where we keep only blo ks satisfying (6) as T b j j and ↑ ( T j − ) b j . If we downsample the operator ↑ ( T j − ) b j ba k to itsoriginal s ale j − , then (6) leads to the band des ribed by the ondition l − ⌊ b j / ⌋ ≤ l ′ ≤ l + ⌊ b j / ⌋ , (7)where ⌊ b j / ⌋ denotes the integer part of b j / . We denote the banded operatoron s ale j − as T ⌊ b j / ⌋ j − , where we retain blo ks satisfying (7).If we now rewrite (4) keeping only blo ks within the bands on ea h s ale, weobtain ˆ g = T f ˆ g = [ T b − ↑ ( T ) b ] f = T b f − ↑ ( T ) b f ˆ g = [ T b − ↑ ( T ) b ] f = T b f − ↑ ( T ) b f (8) . . . . . . . . . ˆ g j = [ T b j j − ↑ ( T j − ) b j ] f j = T b j j f j − ↑ ( T j − ) b j f j . . . . . . . . . For any arbitrary but (cid:28)nite a ura y, instead of applying the full [ T j − ↑ ( T j − )] , we will only apply its banded approximation.A simple but important observation is that ↓ ([ ↑ ( T j − )] f j ) = T j − f j − , (9)whi h follows from the fa t that Q j P j = P j Q j = 0 , sin e these are orthogonalproje tions. Thus, we observe that ↓ ( ↑ ( T j − ) b j ) = T ⌊ b j / ⌋ j − ; so instead of ap-plying ↑ ( T j − ) b j f j on s ale j , we an obtain the same result using (9), so that10igure 2. Modi(cid:28)ed non-standard form of the onvolution operator in (11) in a mul-tiwavelet basis, with white representing 0 and bla k representing large values. Thetop matrix is the proje tion of this operator on V , resulting in a dense matrix. Thelower half depi ts the multiresolution representation in (10) with only the blo ksthat are a tually retained for a given a ura y. We will all the leftmost matrix inthis series a whole band matrix and all others outer band matri es. The two emptyouter band matri es on s ales j = 0 , are explained in the main text. ↑ (cid:16) T ⌊ b j / ⌋ j − f j − (cid:17) = ↑ ( T j − ) b j f j . Therefore, we will ompute only T ⌊ b j / ⌋ j − f j − ons ale j − and ombine it with omputing T b j − j − f j − . In orporating this into(5), we arrive at g n = T b n n f n + ↑ h(cid:16) T b n − n − − T ⌊ b n / ⌋ n − (cid:17) f n − + (10) ↑ h(cid:16) T b n − n − − T ⌊ b n − / ⌋ n − (cid:17) f n − + . . . + h ↑ h(cid:16) T − T ⌊ b / ⌋ (cid:17) f ii . . . ii . Using this expression yields an e(cid:30) ient algorithm for applying an operator,as on ea h s ale j , T b j j − T [ b j +1 / j is a sparse obje t, due to the an ellationwhi h o urs for most of the blo ks. In parti ular, T b j j − T [ b j +1 / j is missingthe blo ks near the diagonal, and we will refer to it as an outer band matrix.We will all T b j j a whole band matrix as it ontains both the inner and outerbands. 11he stru ture of these matri es is illustrated in Figure 2. The two emptys ales j = 0 , arise due to the omplete an ellation of blo ks on these s ales.We note that the modi(cid:28)ed non-standard form is onstru ted adaptively inthe number of s ales ne essary for a given fun tion. For just two s ales, this onstru tion will have the s ale j = 1 non-empty. Given an adaptive de om-position of a fun tion on n s ales, we pre ompute the modi(cid:28)ed non-standardform (depi ted in Figure 2 for (cid:28)ve s ales) on all s ales n, n − , . . . , . Formatri es requiring n · n blo ks on the (cid:28)nest s ale n , we need to keep andapply only O (2 n ) blo ks, as with the original non-standard form in [4℄.2.4.1 Adaptive appli ationLet us show how to use the multis ale representation in (10) to apply theoperator T to a fun tion f with ontrolled a ura y ǫ . We des ribe an adaptiveappli ation of the operator to a fun tion, where we assume (as is often the ase) that the tree stru ture of the input is su(cid:30) ient to adequately des ribethe output with a ura y ǫ . This assumption will be removed later.Our algorithm uses the stru ture shown in Figure 2 in an adaptive fashion.We opy the stru ture of the redundant tree for the input fun tion, and usethat as a template to be (cid:28)lled for the output g . For ea h node of the outputtree we determine whether it is a leaf or an internal node: for leaves, we mustapply a whole band matrix on the s ale of that node, su h as the leftmostmatrix for j = 5 in the example shown in Figure 2. For internal nodes, weapply an outer band matrix (for that s ale). We note that our onstru tionof the operator produ es both whole and outer band matri es for all s ales,and we simply hoose the appropriate kind for ea h node of output as needed.Upon ompletion of this pro ess, we apply the proje tion (5) to onstru t the(cid:28)nal adaptive tree representing the output.Algorithm 2 returns an adaptive tree representing the fun tion g . This tree ontains su(cid:30) ient information to evaluate g at arbitrary points by interpola-tion and may be used as an input in further omputations.We note that Step 2 in Algorithm 2 naturally resolves the problem that isusually addressed by mortar methods, see e.g. [19,20,21,22℄. Sin e adaptiverepresentations have neighboring blo ks of di(cid:27)erent sizes, they en ounter dif-(cid:28) ulties when applying non-diagonal operators, as they require blo ks whi hdo not exist on that s ale. Our approa h simply onstru ts these as neededand a hes them for reuse, without requiring any additional onsideration onthe part of the user. 12lgorithm 2 Adaptive non-standard form operator appli ation in d = 1 , g = T f Initialization: Constru t the redundant tree for f and opy it as skeletontree for g (see Se tion 2.2.1).for all s ales j = 0 , . . . , n − dofor all fun tion blo ks g jl in the tree for g at s ale j do T jll ′ (see Se tion 2.3):if g jl belongs to a leaf thenRead operator blo ks T jll ′ from row l of whole band matrix T b j j .elseRead operator blo ks T jll ′ from row l of outer band matrix T b j j − T [ b j +1 / j .end if f jl ′ of the input fun tion f :if fun tion blo k f jl ′ exists in the redundant tree for f thenRetrieve it.elseCreate it by interpolating from a oarser s ale and a he for reuse.end if ˆ g jl = P l ′ T jll ′ f jl ′ , where the operation T jll ′ f jl ′ indi ates a regular matrix-ve tormultipli ation.end forend for ˆ g jl on all s ales into a proper adap-tive tree by using Eq. (5).Dis ard from the resulting tree unne essary fun tion blo ks at the requesteda ura y.Return: the fun tion g represented by its adaptive tree.2.4.2 Numeri al exampleLet us brie(cid:29)y illustrate the appli ation of the modi(cid:28)ed ns-form with an exam-ple of a singular onvolution on the unit ir le, the operator with the kernel K ( x ) = cot( πx ) , ( Cf )( y ) = p.v. Z cot( π ( y − x )) f ( x ) dx, (11)a periodi analogue of the Hilbert transform. In order to (cid:28)nd its representationin multiwavelet bases, we ompute 13 j ; lii ′ = 2 − j Z − K (2 − j ( x + l )) Φ ii ′ ( x ) dx = 2 − j Z − cot( π − j ( x + l )) Φ ii ′ ( x ) dx , (12)where Φ ii ′ ( x ) , i, i ′ = 0 , . . . , k − are ross- orrelation fun tions des ribedin Appendix 8.4 and l = 0 , ± , ± , . . . j − . We ompute r j ; lii ′ using the onvergent integrals r j ; lii ′ = 2 − j i ′ + i X k = i ′ − i c kii ′ Z Φ + k, ( x ) (cid:16) cot( π − j ( x + l )) + ( − i + i ′ cot( π − j ( − x + l )) (cid:17) dx, where Φ + k, is a polynomial des ribed in Appendix 8.4. In our numeri al ex-periment, we apply (11) to the periodi fun tion on [0 , , f ( x ) = X k ∈ Z e − a ( x + k − / , whi h yields ( Cf )( y ) = − X k ∈ Z e − a ( y + k − / Er(cid:28) [ √ a ( y + k − / i r πa X n ∈ Z sign ( n ) e − n π /a e πiny , (13)where e − y Er(cid:28) ( y ) = √ π R y e s − y ds. Expression (13) is obtained by (cid:28)rst observ-ing that the Hilbert transform of e − ax is − e − ay Er(cid:28) ( √ ay ) , and then evaluatingthe latti e sum, noting that (see [23, formula 4.3.91℄) cot( πx ) = 1 π ( 1 x + ∞ X k =1 xx − k ) . Table 1 summarizes the numeri al onstru tion of this solution for a = 300 ,at various requested pre isions. Optimal performan e is obtained by adjustingthe order of the basis p as a fun tion of the requested pre ision, to ensure thatthe operator remains a banded matrix with small band, and that the adaptiverepresentation of the input fun tion requires a moderate number of s ales.The resulting numeri al error (as ompared to the exa t analyti al solution),measured in the ℓ norm, is shown in the last olumn. Figure 3 shows theinput and results for this example, as well as the point-wise error for the asewhere ǫ req = 10 − and p = 14 (the last row in the table).3 Separated representations of integral kernelsThe approa h we've dis ussed so far does not e(cid:30) iently generalize to theappli ation of non-separable multidimensional integral kernels. Sin e several14 S ales N blo ks ǫ E − . · − − . · −
11 [2,4,5℄ 14 − . · −
14 [3,4,5℄ 16 − . · − Table 1Results from evaluating (13) with our algorithm. The order of the basis p is adjustedas a fun tion of the requested pre ision ǫ . The se ond olumn indi ates s ales presentin the adaptive tree for the input. The third olumn shows the total number ofblo ks of oe(cid:30) ients in this tree. The last olumn ( E ) shows the a tual error of the omputed solution in the ℓ norm. . . . . . . − . . . . f ( x ) = e ( − a ( x − ) ) [cot( f )]( x ) . . . . . . − − − − − − − − Figure 3. Results of applying the otangent kernel to a periodized Gaussian usingbasis of order p = 14 (the last row in Table 1). The pointwise error is shown on theright for a requested a ura y of ǫ = 10 − .physi ally important kernels belong to this ategory (e.g. the Poisson kernelsin d = 2 and d =3), additional tools are needed. We now des ribe the key ideathat allows us to perform this generalization to d > .We use the separated representation of operators introdu ed in [6,24℄ to re-du e the omputational ost of the straightforward generalization of the mul-tiresolution approa h in [4℄. Su h representations are parti ularly simple for onvolution operators and are based on approximating kernels by a sum ofGaussians [24,8,9,11,10℄. This approximation has a multiresolution hara terby itself and requires a remarkably small number of terms. In fa t, our algo-rithm uses the oe(cid:30) ients and the exponents of the Gaussian terms as the onlyinput from whi h it sele ts the ne essary terms, s ale by s ale, a ording thedesired a ura y threshold ǫ . Therefore, our algorithm works for all operatorswith kernels that admit approximation by a sum of Gaussians. Examples ofsu h operators in lude the Poisson and the bound state Helmholtz Green'sfun tions, the proje tor on divergen e free fun tions, as well as regular and15ra tional derivative operators. Let us onsider a parti ular family of operators ( − ∆ + µ I ) − α , where µ ≥ and < α < / . The kernel of this operator K µ,α ( r ) = 2 − + α · C α · ( µr ) − α K − α ( µr ) , where K − α is the modi(cid:28)ed Bessel fun tion, r = || x − y || and C α = 2 · − α π − / Γ( α ) , has an integral representation K µ,α ( r ) = C α Z ∞−∞ e − r e s e − µ e − s +(3 − α ) s ds. (14)Using the trapezoidal rule, we onstru t an approximation valid over a rangeof values δ ≤ r ≤ R with a ura y ǫ , of the form (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) K µ,α ( r ) − M X m =1 w m e − τ m r (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ≤ ǫK ,α ( r ) = ǫ Γ(3 / − α ) · C α r − α , (15)where τ m = e s m , w m = h C α e − µ e − sm / − α ) s m , h = ( B − A ) /M and s m = A + mh . The limits of integration, A , B and the step size h are sele ted asindi ated in [24℄, where it is shown that for a (cid:28)xed a ura y ǫ the number ofterms M in (15) is proportional to log( Rδ − ) . Although it is possible to sele t δ and R following the estimates in [25℄ and optimize the number of terms fora desired a ura y ǫ , in this paper we start with an approximation that has anobviously ex essive range of validity and thus, an ex essive number of terms.An example of su h approximation is shown in Figure 4. For a requested toler-an e of ǫ = 10 − , roughly 300 terms are enough to provide a valid approxima-tion over a range of de ades. We then let the algorithm hoose the ne essaryterms, s ale by s ale, to satisfy the user-supplied a ura y requirement ǫ . Thisapproa h may end up with a few extra terms on some s ales in omparisonwith that using a nearly optimal number of terms [8,10,11,9℄. Whereas the ost of applying a few extra terms is negligible, we gain signi(cid:28) antly in havinga mu h more (cid:29)exible and general algorithm.We note that approximation in (15) learly redu es the problem of applyingthe operator to that of applying a sequen e of Gauss transforms [26,27℄, oneby one. From this point of view, the algorithm that we present may be on-sidered as a pro edure for applying a linear ombination of Gauss transformssimultaneously.In order to represent the kernel K of the operator in multiwavelet bases, weneed to ompute the integrals, 16 − − − − − − − − − − − − − − − − − − − | RelativeError | Figure 4. Relative error of the Gaussian approximation for the Poisson kernel in 3dimensions. This unoptimized expansion uses 299 terms to over a dynami rangeof roughly 15 de ades with ǫ = 10 − relative a ura y. r j ; ℓi i ′ ,i i ′ ,i i ′ = M X m =1 − j Z B w m e − τ m k − j ( x + ℓ ) k Φ i i ′ ( x ) Φ i i ′ ( x ) Φ i i ′ ( x ) d x , (16)where Φ ii ′ ( x ) are the ross- orrelations of the s aling fun tions (see Appendix).We obtain r j ; ℓi i ′ ,i i ′ ,i i ′ = M X m =1 w m F j ; m,l i i ′ F j ; m,l i i ′ F j ; m,l i i ′ , (17)where F j ; m,lii ′ = 12 j Z − e − τ m ( x + l ) / j Φ ii ′ ( x ) dx. (18)Sin e the Gaussian kernel is not homogeneous, we have to ompute integrals(18) for ea h s ale. Although in prin iple l ∈ Z , in the next se tion we explainhow to restri t it to a limited range on ea h s ale, for a given a ura y ǫ .4 Modi(cid:28)ed ns-form of a multidimensional operatorIn this se tion, we des ribe how the separated representation approximationsof Se tion 3 an be used to onstru t a multidimensional extension of the ns-form representation from Se tion 2.3, using only one-dimensional quantitiesand norm estimates. This makes our approa h viable for d > . We use themodi(cid:28)ed ns-form as in the one-dimensional ase des ribed in Se tion 2.3. We17nd the ns-form essential for adaptive algorithms in more than one dimension,sin e:(1) S ales do not intera t as the operator is applied. All intera tions betweens ales are a ounted for by the (inexpensive) proje tion at the (cid:28)nal stepof the algorithm.(2) For the same reason, the subdivision of spa e at di(cid:27)erent s ales natu-rally maps into the supporting data stru tures. We note that one of themain di(cid:30) ulties in developing adaptive algorithms is in organizing om-putations with blo ks of an adaptive de omposition of a fun tion fromdi(cid:27)erent s ales but with a ommon boundary. The methods for su h om-putations are known as mortar elements methods. In our approa h thisissue does not present any obsta le, as all relevant intera tions are natu-rally a ounted for by the data stru tures.The key feature that makes our approa h e(cid:30) ient in dimensions d ≥ isthe separated stru ture of the modi(cid:28)ed ns-form. Namely, the blo ks of T j − ↑ ( T j − ) are of the form (for d = 3 ) T ℓj − ↑ ( T ℓj − ) = M X m =1 w m F j ; ml F j ; ml F j ; ml − ↑ M X m =1 w m F j − ml F j − ml F j − ml ! = M X m =1 w m F j ; ml F j ; ml F j ; ml (19) − M X m =1 w m ↑ ( F j − ml ) · ↑ ( F j − ml ) · ↑ ( F j − ml ) . As in the ase d = 1 , the norm of the operator blo ks of T ℓj − ↑ ( T ℓj − ) de aysrapidly with k ℓ k , ℓ = ( l , l , l ) , and the rate of de ay depends on the numberof vanishing moments of the basis [4℄. Moreover, we limit the range of shiftindi es k ℓ k using only one-dimensional estimates of the di(cid:27)eren es F j ; m, lii ′ − ↑ F j ; m,lii ′ and F j ; m, l +1 ii ′ − ↑ F j ; m,lii ′ (20)of operator blo ks omputed via (18). The norms of individual blo ks F j ; m,lii ′ are illustrated in Figure 5 (a), for s ale j = 1 . By sele ting the number of vanishing moments for a given a ura y, it is su(cid:30)- ient to use k ℓ k ∞ ≤ in pra ti al appli ations that we have en ountered. Also,not all terms in the Gaussian expansion of an operator need to be in ludedsin e, depending on the s ale j , their ontribution may be negligible for agiven a ura y, as shown in Figure 5 (b). Below we detail how we sele t terms18f the Gaussian expansion on a given s ale as well as the signi(cid:28) ant blo ks of T ℓj − ↑ ( T ℓj − ) . This pro edure establishes the band stru ture of the operator.We then proje t the banded operator ↑ ( T b j j − ) ba k to the s ale j − andthen ombine blo ks on the natural s ale for ea h proje tion in order to applythe operator e(cid:30) iently as was explained in Se tion 2.3 for the one-dimensional ase.We note that in de iding whi h terms to keep in (19), we do not ompute thedi(cid:27)eren e between the full three dimensional blo ks as it would arry a high omputational ost; instead we use estimates based on the one dimensionalblo ks of the separated representation. We note that sin e the resulting bandstru ture depends only on the operator and the desired a ura y of its ap-proximation, one of the options is to store su h band information as it is likelyto be reused.In order to e(cid:30) iently identify the signi(cid:28) ant blo ks in T ℓj − ↑ ( T ℓj − ) as afun tion of ℓ , we develop norm estimates based only on the one-dimensionalblo ks. The di(cid:27)eren e between two terms of the separated representation, say F F F − G G G , may be written as F F F − G G G = ( F − G ) F F + G ( F − G ) F + G G ( F − G ) . We average six di(cid:27)erent ombinations of the three terms to in lude all dire -tions in a symmetri manner, whi h results in the norm estimate k F F F − G G G k ≤ sym [ k F − G kk F kk F k + (21) k G kk F − G kk F k + k G kk G kk F − G k ] , where the symmetrization is over the three dire tions and generates terms.For the rotationally symmetri operators with Gaussian expansion as in (15) omputing the right hand side in this estimate involves just three types of onedimensional blo ks and their norms, N j ; m ; l dif = (cid:13)(cid:13)(cid:13) F j ; m ; l − ↑ ( F j − m ; l ) (cid:13)(cid:13)(cid:13) ,N j ; m ; lF = (cid:13)(cid:13)(cid:13) F j ; m ; l (cid:13)(cid:13)(cid:13) , (22) N j ; m ; l ↑ F = (cid:13)(cid:13)(cid:13) ↑ ( F j − m ; l ) (cid:13)(cid:13)(cid:13) , where index j indi ates the s ale, m the term in the Gaussian expansion (15),and l the position of the blo k in a given dire tion.These estimates allow us to dis ard blo ks whose norm falls below a giventhreshold of a ura y, namely, for ea h multi-index ℓ , we estimate19
60 180 200 220 240 260 280 300
Separation index − − − − − − − − − − − − − − − N o r m -2-10123Cutoff (a) Norms of ea h one-dimensionalblo k omputed via (18) for s ale j = 1 ,as a fun tion of the index m in the sep-arated representation.
160 180 200 220 240 260 280 300
Separation index − − − − − − − − − − − − − − − N o r m -2-10123Cutoff (b) Norm estimates (23) for s ale j = 1 as a fun tion of the index m in the sep-arated representation. Based on theseestimates, only terms above the uto(cid:27)are a tually applied.Figure 5. Comparison of norms of matrix blo ks generated by the Gaussian ap-proximation for the Poisson kernel in dimension d = 3 . In ea h pi ture, the urves orrespond to the di(cid:27)erent o(cid:27)sets l for whi h blo ks are generated. Figure (b) illus-trates the estimate in (23) (see main text for details). (cid:13)(cid:13)(cid:13) T ℓj − ↑ ( T ℓj − ) (cid:13)(cid:13)(cid:13) ≤ M X m =1 w m sym h N j ; m ; l dif N j ; m ; l F N j ; m ; l F + (23) N j ; m ; l ↑ F N j ; m ; l dif N j ; m ; l F + N j ; m ; l ↑ F N j ; m ; l ↑ F N j ; m ; l dif i . For ea h s ale j and ea h blo k T ℓj − ↑ ( T ℓj − ) labeled by the multi-index ℓ = ( l , l , l ) , we ompute all terms of the sum in (23) and identify the range [ m , m ] whi h we need to keep for that blo k, by dis arding from the sumterms whose umulative ontribution is below ǫ . If the entire sum falls below ǫ ,this range may be empty and the entire T ℓj − ↑ ( T ℓj − ) is dis arded. The range [ m , m ] di(cid:27)ers signi(cid:28) antly depending whether or not the blo k is a(cid:27)e tedby the singularity of the kernel as is illustrated in Figure 5. In Figure 5 (a)the rate of de ay for the blo ks with shift | l | = 2 , is signi(cid:28) antly faster thanfor the blo ks with | l | ≤ a(cid:27)e ted by the singularity. We note that all blo ksof the (cid:28)rst 150 terms in the separated representation (17) have norm 1 (andrank 1 as matri es) and are not shown in Figure 5 (a).Sin e the di(cid:27)eren e in (23) involves blo ks upsampled from a oarser s ale, allshifts | l | ≤ are a(cid:27)e ted by the singularity. Figure 5 (b) shows the r.h.s. ofthe estimate in (23) for di(cid:27)erent shifts along one of the dire tions, where theblo ks along the other two dire tions are estimated by the maximum normover all possible shifts.After dis arding blo ks with norms less that ǫ using the estimate in (23), wedownsample the remaining blo ks of ↑ ( T ℓj − ) ba k to the original s ale. This20eaves only blo ks of T ℓj on the s ale j and we remove additional blo ks of T ℓj for the shifts | l | = 2 , where the de ay is su(cid:30) iently fast to make their ontribution less than ǫ .This leads us to arrange the blo ks on ea h s ale into several subsets by thee(cid:27)e t the singularity of the kernel has on them and (cid:28)nd the appropriate range [ m , m ] for ea h subset. There are three su h sets in dimension d = 2 andfour sets in dimension d = 3 . For ea h index l, we will say that the indexbelongs to the ore if l = − , , and to the boundary otherwise. The oreindi es orrespond to one-dimensional blo ks whose de(cid:28)ning integrals in ludethe singularity of the kernel. We then divide all possible values of the multi-index ℓ = ( l , l , l ) , a ording to the number of ore indi es it has. In d = 3 this gives us four sets: • Core: all indi es ( l , l , l ) belong to the ore. • Boundary-1: two of the indi es belong to the ore and one to the boundary • Boundary-2: one of the indi es belongs to the ore, the other two to theboundary • Boundary-3: all indi es ( l , l , l ) belong to the boundary.We then (cid:28)nd the range [ m , m ] for ea h subset and apply blo ks of ea hsubset separately, thus avoiding unne essary omputations with blo ks whose ontribution is negligible. This range analysis only needs to be done on e peroperator and the desired a ura y, and the results may be saved for repeateduse.5 Multidimensional adaptive appli ation of ns-formIn this se tion we present an algorithm for applying the modi(cid:28)ed non-standardform whi h is an extension of (2) (based on (8) and (10)) to higher dimensions.We are now seeking to ompute g ( x ) = [ T f ]( x ) = Z K ( y − x ) f ( y ) d y , where x , y ∈ R d for d = 2 , . The separated approximation (19) redu es the omplexity of applying the operator by allowing partial fa torization of thenested loops in ea h s ale indi ated by the order of summation and illustratedfor d = 3 , 21 j ; l l l = X m w m X l ′ F j ; m ; l − l ′ X l ′ F j ; m,l − l ′ X l ′ F j ; m,l − l ′ −↑ X l ′ F j − m ; l − l ′ X l ′ F j − m,l − l ′ X l ′ F j − m,l − l ′ f j ; l ′ l ′ l ′ . As des ribed in the previous se tion, this evaluation is done by regions of in-di es. These regions of indi es are organized so that (for a given a ura y) thenumber of retained terms of the separated representation is roughly the samefor all blo ks within ea h region. Thus, we perform the summation over theterms of the separated representation last, applying only the terms that a tu-ally ontribute to the result above the requested a ura y threshold, a ordingto estimate (23). Therefore, we avoid introdu ing he ks per individual blo kand the resulting loss of performan e.Just as in the one-dimensional ase, we use (19) in a `natural s ale' manner.That is, blo ks belonging to s ale j are only applied on that s ale. As in one-dimensional ase, the intera tion between s ales is a hieved by the proje tion(10) that redistributes blo ks a umulated in this manner properly betweenthe s ales to obtain the adaptive tree for the resulting fun tion. The overallapproa h is the same as des ribed in (2). We note that, as expe ted, the sep-arated representation requires more detailed bookkeeping when onstru tingthe data stru tures for the operator.Remark 1 Our multiresolution de omposition orresponds to the geometri- ally varying re(cid:28)nement in (cid:28)nite element methods. In this ase the adjoiningboxes do not ne essarily share ommon verti es, forming what orresponds tothe so- alled non- onforming grid. In (cid:28)nite element methods su h situationrequires additional onstru tion provided by the mortar element methods.Mortar element methods were introdu ed by Patera and his asso iates, seee.g. [19,20,21,22℄. These methods permit oupling dis retizations of di(cid:27)erenttypes in non-overlapping domains. Su h methods are fairly ompli ated andinvolve, for example, the introdu tion of interfa e onditions through an L minimization. In our approa h we do not fa e these issues at all and do nothave to introdu e any additional interfa e onditions. The proper onstru -tion for adjoining boxes is taken are by the redundant tree data stru tureand Step 2 of Algorithm 3 for applying the kernel, whi h generates the ne es-sary missing boxes on appropriate s ales.Remark 2 Although Algorithm 3 applies onvolution operators, only minor hanges are needed to use it for non- onvolutions. Of ourse in su h ase, theseparated representation for the modi(cid:28)ed ns-form should be onstru ted by adi(cid:27)erent approa h. 22lgorithm 3 Adaptive non-standard form operator appli ation in multipledimensions (illustrated for d = 2 ), g = T f Initialization: Constru t the redundant tree for f and opy it as skeletontree for g (see Se tion 2.2.1).for all s ales j = 0 , . . . , n − dofor all fun tion blo ks g jl l in the tree for g at s ale j do F j ; m ; l − l ′ , F j ; m ; l − l ′ (see Se tion 4):if g jl l belongs to a leaf thenRead operator blo ks F j ; m ; l − l ′ , F j ; m ; l − l ′ for Core, Boundary-1 andBoundary-2, and their weights w m and orresponding ranges fromthe separated representation.elseRead operator blo ks F j ; m ; l − l ′ , F j ; m ; l − l ′ for Boundary-1 andBoundary-2, and their weights w m and orresponding ranges fromthe separated representation.end if f jl ′ l ′ of the input fun tion f :if fun tion blo k f jl ′ l ′ exists in the redundant tree for f thenRetrieve it.elseCreate it by interpolating from a oarser s ale and a he for reuse.end if S of indi es determined in Step 1 (Core, Boundary-1,Boundary-2) and the orresponding ranges of terms in the separatedrepresentation, ompute the sum ˆ g j ; Sl l = X m w m X l ′ F j ; m ; l − l ′ X l ′ F j ; m ; l − l ′ f jl ′ l ′ , Add all omputed sums to obtain ˆ g jl l .end forend for ˆ g jl l on all s ales into a properadaptive tree by using Eq. (5).Dis ard from the resulting tree unne essary fun tion blo ks at the requesteda ura y.Return: the fun tion g represented by its adaptive tree.23.1 Operation ount estimatesAdaptive de omposition of fun tionsThe ost of adaptively de omposing a fun tion in d dimensions is essentiallythat of an adaptive wavelet transform. Spe i(cid:28) ally, it takes O ( N blo ks · p d +1 ) operations to ompute su h representation, where N blo ks is the (cid:28)nal numberof signi(cid:28) ant blo ks in the representation and p is the order of multiwaveletbasis hosen. In omparison with the usual wavelet transform, it appears tobe signi(cid:28) antly more expensive. However, these O ( p d +1 ) operations pro ess O ( p d ) points, thus in ounting signi(cid:28) ant oe(cid:30) ients as it is done in the usualadaptive wavelet transform, we end up with O ( p ) operations per point.Operator appli ationThe ost of applying an operator in the modi(cid:28)ed ns-form is O ( N blo ks M p d +1 ) ,where N blo ks is the number of blo ks in the adaptive representation of theinput fun tion, M is the separation rank of the kernel in (15) and p is theorder of the multiwavelet basis. For a given desired a ura y ǫ , we typi allysele t p ∝ log ǫ − ; M has been shown to be proportional to (log ǫ − ) ν , where ν depends on the operator [24℄. In our numeri al experiments, M is essentiallyproportional to log ǫ − , sin e we never use the full separated representation,as dis ussed in Se tion 4 and illustrated in Figure 5.This operation ount an be potentially redu ed to O ( N blo ks M p d ) by usingthe stru ture of the matri es in 18, and we plan to address this in the future.Final proje tionAfter the operator has been applied to a fun tion in a s ale-independent fash-ion, a (cid:28)nal proje tion step is required as dis ussed in Se tion 2.3. This steprequires O ( N blo ks · p d ) operations, the same as in the original fun tion de om-position. In pra ti e, this time is negligible ompared to the a tual operatorappli ation. 24 Numeri al examples6.1 The Poisson equationWe illustrate the performan e of the algorithm by solving the Poisson equationin d = 3 ∇ φ ( r ) = − ρ ( r ) (24)with free spa e boundary onditions, φ ( r ) → and ∂φ/∂r → as r → ∞ . Wewrite the solution as φ ( r ) = 14 π Z R | r − r ′ | ρ ( r ′ ) d r ′ and adaptively evaluate this integral. We note that our method an equallybe used for d = 2 , sin e the orresponding Green's fun tion an also be rep-resented as a sum of Gaussians, and the operator appli ation algorithm hasbeen implemented in for both d = 2 and d = 3 .For our test we sele t φ ( r ) = X i =1 e − α | r − r ′ | , so that we solve the Poisson equation with ρ ( r ) = −∇ φ ( r ) = − X i =1 (4 α | r − r i | − α ) e − α | r − r i | . Our parameters are hosen as follows: α = 300 , r = (0 . , . , . , r =(0 . , . , . and r = (0 . , . , . . These ensure that ρ ( r ) is well belowour requested thresholds on the boundary of the omputational domain. Allnumeri al experiments were performed on a Pentium-4 running at 2.8 GHz,with 2 GB of RAM. The results are summarized in Table 2.In order to gauge the speed of algorithm in reasonably omputer-independentterms, we use a similar approa h to that of [17℄ and also provide timings ofthe Fast Fourier Transform (FFT). Spe i(cid:28) ally, we display timings for twoFFTs as an estimate of the time needed to solve the Poisson equation witha smooth right hand side and periodi boundary onditions in a ube. As in[17℄, we ompute the rate that estimates the number of pro essed points perse ond. We observe that for our adaptive algorithm su h rate varies between . · and . · (see Table 2), whereas for the FFTs it is around (seeTable 3). We note that our algorithm is not fully optimized, namely, we do notuse the stru ture of the matri es in (18) and the symmetries a(cid:27)orded by theradial kernels. We expe t a substantial impa t on the speed by introdu ingthese improvements and will report them separately.25equested ǫ = 10 − p E E ∞ Time (s) Rate (pts/s)6 . · − . · − . . · . · − . · − .
51 7 . · . · − . · − .
68 1 . · Requested ǫ = 10 − p E E ∞ Time (s) Rate (pts/s)10 . · − . · − . · . · − . · − . · . · − . · − . · Requested ǫ = 10 − p E E ∞ Time (s) Rate (pts/s)16 . · − . · − . · . · − . · − . · . · − . · − . · Table 2A ura y and timings for the adaptive solution of the Poisson equation in (24).size Time (s) Rate (pts/s) . · . · . · . · Table 3Timings of two 3D FFTs to estimate the speed of a non-adaptive, periodi Poissonsolver on a ube for smooth fun tions.We note that the multigrid method (see e.g. [28,29℄) is frequently used as atool for solving the Poisson equation (and similar problems) in di(cid:27)erentialform. The FFT-based gauge suggested in [17℄ is useful for omparisons withthese algorithms as well. 26igure 6. A two-dimensional sli e of the three-dimensional subdivision of spa e bythe s aling fun tions and an illustration of the sour e term for the Poisson equation(24).6.2 The ground state of the hydrogen atomA simple example of omputing the ground state of the hydrogen atom illus-trates the numeri al performan e of the algorithms developed in this paper,and their utility for onstru ting more omplex odes. The eigenfun tions ψ for the hydrogen atom satisfy the time-independent S hrödinger equation(written in atomi units and spheri al oordinates), −
12 ∆ ψ − r ψ = Eψ, (25)where E is the energy eigenvalue. For the ground state, E = − / and the(unnormalized) wave fun tion is ψ = e − r . Following [30℄, we write φ = − G µ V φ, (26)where G µ = ( − ∆ + µ I ) − is the Green's fun tion for some µ and V = − /r is the nu lear potential. For µ = √− E the solution φ of (26) has k φ k = 1 and oin ides with that of (25). We solve (26) by a simple iteration startingfrom some value µ and hanging µ to obtain the solution with k φ k = 1 . Thealgorithm pro eeds as follows:(1) Initialize with some value µ and fun tion φ . The number of iterations ofthe algorithm is only weakly sensitive to these hoi es.(2) Compute the produ t of the potential V and the fun tion φ .(3) Apply the Green's fun tion G µ to the produ t V φ via the algorithm of27igure 7. Convergen e of the iteration to obtain the ground state of the hydrogenatom omputed via formulation in (26) for the non-relativisti S hrödinger equation.The requested a ura y in applying the Green's fun tion is set to − .this paper to ompute φ new = − G µ V φ. (4) Compute the energy for φ new , E new = h∇ φ new , ∇ φ new i + h V φ new , φ new ih φ new , φ new i . (5) Set µ = √− E new , φ = φ new / k φ new k and return to Step 2.The iteration is terminated as the hange in µ and k φ k − falls below thedesired a ura y. The progress of the iteration is illustrated in Figure 7. The omputations in Steps 2 and 4 use the three dimensional extension of theapproa h des ribed in [13℄, to ompute point-wise multipli ations of adaptivelyrepresented fun tions and weak di(cid:27)erential operators of the same.This example illustrates an appli ation of our algorithm to problems in quan-tum hemistry. Multiresolution quantum hemistry developed in [8,9,10,11℄also uses separated representations. The main te hni al di(cid:27)eren e with [8,9,10,11℄is that we use the modi(cid:28)ed non-standard form and apply operator blo ks ontheir (cid:16)natural(cid:17) s ale, thus produ ing a fully adaptive algorithm. We are ur-rently using this algorithm as part of a new method for solving the multipar-ti le S hrödinger equation and will report the results elsewhere.7 Dis ussion and on lusionsWe have shown that a ombination of separated and multiresolution repre-sentations of operators yields a new multidimensional algorithm for applying28 lass of integral operators with radial kernels. We note that the same algo-rithm is used for all su h operators as they are approximated by a weightedsum of Gaussians. This fa t makes our approa h appli able a ross multiple(cid:28)elds, where a single implementation of the ore algorithm an be reused fordi(cid:27)erent spe i(cid:28) problems. The algorithm is fully adaptive and avoids issuesusually addressed by mortar methods.The method of approximation underlying our approa h is distin t from that ofthe FMM, has similar e(cid:30) ien y and has the advantage of being more readilyextendable to higher dimensions. We also note that semi-analyti approxima-tions via weighted sums of Gaussians provide additional advantages in someappli ations. Although we des ribed the appli ation of kernels in free spa e,there is a simple extension to problems with radial kernels subje t to periodi ,Diri hlet or Neumann boundary onditions on a ube that we will des ribeseparately.The algorithm may be extended to lasses of non- onvolution operators, e.g.,the Calderon-Zygmund operators. For su h extensions the separated represen-tation may not be available in analyti form, as it is for the operators of thispaper, and may require a numeri al onstru tion. The separated representa-tion of the kernel permits further generalization of our approa h to dimensions d ≫ for applying operators to fun tions in separated representation.A notable remaining hallenge is an e(cid:30) ient, high order extension of this te h-nique to the appli ation of operators on domains with ompli ated geometriesand surfa es.8 Appendix8.1 S aling fun tionsWe use either the Legendre polynomials P , . . . , P p − or the interpolating poly-nomials on the Gauss-Legendre nodes in [ − , to onstru t an orthonormalbasis for ea h subspa e V j [12,13℄.Let us brie(cid:29)y des ribe some properties of the Legendre s aling fun tions φ k , k = 0 , . . . , p − , de(cid:28)ned as φ k ( x ) = √ k + 1 P k (2 x − , x ∈ [0 , , x / ∈ [0 , , (27)and forming a basis for V . The subspa e V j is spanned by j p fun tions29btained from φ , . . . , φ p − by dilation and translation, φ jkl ( x ) = 2 j/ φ k (2 j x − l ) , k = 0 , . . . , p − , l = 0 , . . . , j − . (28)These fun tions have support on [2 − j l, − j ( l + 1)] and satisfy the orthonor-mality ondition Z ∞−∞ φ jkl ( x ) φ jk ′ l ′ ( x ) dx = δ kk ′ δ ll ′ . (29)A fun tion f , de(cid:28)ned on [0 , , is represented in the subspa e V j by its nor-malized Legendre expansion f ( x ) = j − X l =0 p − X k =0 s jkl φ jkl ( x ) , (30)where the oe(cid:30) ients s jkl are omputed via s jkl = Z − j ( l +1)2 − j l f ( x ) φ jkl ( x ) dx. (31)As long as f ( x ) is smooth enough and is well approximated on [2 − j l, − j ( l +1)] by a polynomial of order up to p − , we may use Gauss-Legendre quadraturesto al ulate the s jkl via s jkl = 2 − j/ p − X i =0 f (2 − j ( x i + l )) φ k ( x i ) w i , (32)where x , . . . , x p − are the roots of P p (2 x − and w , . . . , w p − are the orre-sponding quadrature weights.In more than one dimension, the above formulas are extended by using a tensorprodu t basis in ea h subspa e. For example, in two dimensions equation (30)be omes f ( x, x ′ ) = j − X l =0 p − X k =0 2 j − X l ′ =0 p − X k ′ =0 s jkk ′ ll ′ φ jkl ( x ) φ jk ′ l ′ ( x ′ ) . (33)8.2 MultiwaveletsWe use pie ewise polynomial fun tions ψ , . . . , ψ p − as an orthonormal basisfor W [12,13℄, Z ψ i ( x ) ψ j ( x ) dx = δ ij . (34)Sin e W ⊥ V , the (cid:28)rst p moments of all ψ , . . . , ψ p − vanish: Z ψ i ( x ) x i dx = 0 , i, j = 0 , , . . . , p − . (35)30he spa e W j is spanned by j p fun tions obtained from ψ , . . . , ψ p − bydilation and translation, ψ jkl ( x ) = 2 j/ ψ k (2 j x − l ) , k = 0 , . . . , p − , l = 0 , . . . , j − , (36)and supported in the interval I jl = [2 − j l, − j ( l + 1)] . A fun tion f ( x ) de(cid:28)nedon [0 , is represented in the multiwavelet basis on n s ales by f ( x ) = p − X k =0 s k, φ k ( x ) + n − X j =0 2 j − X l =0 p − X k =0 d jkl ψ jkl ( x ) (37)with the oe(cid:30) ients d jkl omputed via d jkl = Z − j ( l +1)2 − j l f ( x ) ψ jkl ( x ) dx. (38)8.3 Two-s ale relationsThe relation between subspa es, V ⊕ W = V , is expressed via the two-s aledi(cid:27)eren e equations, φ k ( x ) = √ p − X k ′ =0 (cid:16) h (0) kk ′ φ k ′ (2 x ) + h (1) kk ′ φ k ′ (2 x − (cid:17) , k = 0 , . . . , p − , (39) ψ k ( x ) = √ p − X k ′ =0 (cid:16) g (0) kk ′ φ k ′ (2 x ) + g (1) kk ′ φ k ′ (2 x − (cid:17) , k = 0 , . . . , p − , (40)where the oe(cid:30) ients h (0) ij , h (1) ij and g (0) ij , g (1) ij depend on the type of polyno-mial basis used (Legendre or interpolating) and its order p . The matri es of oe(cid:30) ients H (0) = { h (0) kk ′ } , H (1) = { h (1) kk ′ } , G (0) = { g (0) kk ′ } , G (1) = { g (1) kk ′ } are the multiwavelet analogues of the quadrature mirror (cid:28)lters in the usualwavelet onstru tion, e.g., [15℄. These matri es satisfy a number of importantorthogonality relations and we refer to [13℄ for omplete details, in luding the onstru tion of the H, G matri es themselves. Let us only state how thesematri es are used to onne t the s aling s jkl and wavelet d jkl oe(cid:30) ients onneighboring s ales j and j + 1 . The de omposition pro edure ( j + 1 → j ) isbased on 31 jkl = p − X k ′ =0 (cid:16) h (0) kk ′ s j +1 k, l + h (1) kk ′ s j +1 k, l +1 (cid:17) , (41) d jkl = p − X k ′ =0 (cid:16) g (0) kk ′ s j +1 k, l + g (1) kk ′ s j +1 k, l +1 (cid:17) ; (42)the re onstru tion ( j → j + 1 ) is based on s j +1 k, l = p − X k ′ =0 (cid:16) h (0) kk ′ s jkl + g (0) kk ′ d jkl (cid:17) , (43) s j +1 k, l +1 = p − X k ′ =0 (cid:16) h (1) kk ′ s jkl + g (1) kk ′ d jkl (cid:17) . (44)8.4 Cross- orrelation of the s aling fun tionsFor onvolution operators we only need to ompute integrals with the ross- orrelation fun tions of the s aling fun tions, Φ ii ′ ( x ) = Z ∞−∞ φ i ( x + y ) φ i ′ ( y ) dy. (45)Sin e the support of the s aling fun tions is restri ted to [0 , , the fun tions Φ ii ′ are zero outside the interval [ − , and are polynomials on [ − , and [0 , of degree i + i ′ + 1 , Φ ii ′ ( x ) = Φ + ii ′ ( x ) , ≤ x ≤ , Φ − ii ′ ( x ) , − ≤ x < , , < | x | , (46)where i, i ′ = 0 , . . . , p − and Φ + ii ′ ( x ) = Z − x φ i ( x + y ) φ i ′ ( y ) dy , Φ − ii ′ ( x ) = Z − x φ i ( x + y ) φ i ′ ( y ) dy . (47)We summarize relevant properties of the ross- orrelation fun tions Φ ii ′ inProposition 3(1) Transposition of indi es: Φ ii ′ ( x ) = ( − i + i ′ Φ i ′ i ( x ) ,(2) Relations between Φ + and Φ − : Φ − i,i ′ ( − x ) = ( − i + i ′ Φ + i,i ′ ( x ) for ≤ x ≤ ,323) Values at zero: Φ ii ′ (0) = 0 if i = i ′ , and Φ ii (0) = 1 for i = 0 , . . . , p − ,(4) Upper bound: max x ∈ [ − , | Φ ii ′ ( x ) | ≤ for i, i ′ = 0 , . . . , p − ,(5) Conne tion with the Gegenbauer polynomials: Φ +00 ( x ) = C ( − / (2 x −
1) + and Φ + l ( x ) = √ l + 1 C ( − / l +1 (2 x − ,for l = 1 , , . . . , where C ( − / l +1 is the Gegenbauer polynomial,(6) Linear expansion: if i ′ ≥ i then we have Φ + ii ′ ( x ) = i ′ + i X l = i ′ − i c lii ′ Φ + l ( x ) , (48)where c lii ′ = l ( l + 1) R Φ + ii ′ ( x )Φ + l ( x )(1 − (2 x − ) − dx, i ′ > i, l ( l + 1) R (Φ + ii ( x ) − Φ +00 ( x )) Φ + l ( x )(1 − (2 x − ) − dx, i ′ = i, (49)for l ≥ and c ii ′ = δ ii ′ .(7) Vanishing moments: we have R − Φ ( x ) dx = 1 and R − x k Φ ii ′ ( x ) dx = 0 for i + i ′ ≥ and ≤ k ≤ i + i ′ − .Proof of these properties an be found in [25℄.8.5 Separated representations of radial fun tionsAs an example, onsider approximating the fun tion /r α by a olle tion ofGaussians. The number of terms needed for this purpose is mer ifully small.We have [24℄Proposition 4 For any α > , < δ ≤ , and < ǫ ≤ min n , α o , thereexist positive numbers τ m and w m su h that (cid:12)(cid:12)(cid:12)(cid:12) r − α − M X m =1 w m e − τ m r (cid:12)(cid:12)(cid:12)(cid:12) ≤ r − α ǫ, for all δ ≤ r ≤ (50)with M = log ǫ − [ c + c log ǫ − + c log δ − ] , (51)where c k are onstants that only depend on α . For (cid:28)xed power α and a ura y ǫ , we have M = O (log δ − ) . 33he proof of Proposition 4 in [24℄ is based on using the trapezoidal rule to dis- retize an integral representation of /r α . Similar estimates may be obtainedfor more general radial kernels using their integral representations as in (14).We note that approximations of the fun tion /r via sums of Gaussians havebeen also onsidered in [31,32,33℄.8.6 EstimatesBy sele ting appropriate ǫ and δ in the separated representations of a radialkernel K (as in Proposition 4), we obtain a separated approximation for the oe(cid:30) ients t j ; ℓii ′ ,jj ′ ,kk ′ = 2 − j Z [ − , K (2 − j ( x + ℓ )) Φ ii ′ ( x ) Φ jj ′ ( x ) Φ kk ′ ( x ) d x . (52)Sin e the number of terms, M , depends logarithmi ally on ǫ and δ , we a hieveany (cid:28)nite a ura y with a very reasonable number of terms. For example, forthe Poisson kernel K ( r ) = 1 / πr , we have the following estimate [25℄,Proposition 5 For any ǫ > and < δ ≤ the oe(cid:30) ients t j ; ℓii ′ ,jj ′ ,kk ′ in (52)have an approximation with a low separation rank, r j ; ℓii ′ ,jj ′ ,kk ′ = M X m =1 w m F m,l ii ′ F m,l jj ′ F m,l kk ′ , (53)su h that if max i | l i | ≥ , then | t j ; ℓii ′ ,jj ′ ,kk ′ − r j ; ℓii ′ ,jj ′ ,kk ′ | ≤ c − j ǫ, (54)and if max i | l i | ≤ , then | t j ; ℓii ′ ,jj ′ ,kk ′ − r j ; ℓii ′ ,jj ′ ,kk ′ | ≤ − j ( c δ + c ǫ ) (55)where ǫ , δ , M , τ m , w m , m = 1 , . . . , M are des ribed in Proposition 4 for α = 1 and c and c are (small) onstants.As des ribed in Se tion 4, our adaptive algorithm sele ts only some of theterms, as needed on a given s ale for the desired a ura y ǫ .34.7 Evaluation of integrals with the ross- orrelation fun tionsWe need to ompute integrals in (18), where the ross- orrelation fun tions Φ ii ′ are given in (45). We note that using (48), it is su(cid:30) ient to ompute F j ; m,li = 12 j Z − e − τ m ( x + l ) / j Φ i ( x ) dx for ≤ i ≤ p − rather than p integrals in (18). Using the relation between Φ − and Φ + in Proposition 3, we have Z − Φ − i ( x ) e − τ m ( x + l ) dx = Z Φ − i ( − x ) e − τ m ( − x + l ) dx = ( − i Z Φ + i ( x ) e − τ m ( x − l ) dx, so that F j ; m,li = Z [ e − τ m ( x + l ) + ( − i e − τ m ( x − l ) ]Φ + i ( x ) dx. L for the sparse representation of integraloperators, SIAM J. Math. Anal 24 (1) (1993) 246(cid:21)262.[13℄ B. Alpert, G. Beylkin, D. Gines, L. Vozovoi, Adaptive solution of partialdi(cid:27)erential equations in multiwavelet bases, J. Comput. Phys. 182 (1) (2002)149(cid:21)190.[14℄ G. Beylkin, Approximations and fast algorithms, in: A.Laine, M.Unser,A.Aldroubi (Eds.), Pro . SPIE: Wavelets: Appli ations in Signal and ImagePro essing IX, Vol. 4478, 2001, pp. 1(cid:21)9.[15℄ I. Daube hies, Ten Le tures on Wavelets, CBMS-NSF Series in AppliedMathemati s, SIAM, 1992.[16℄ C. Chui, An Introdu tion to Wavelets, A ademi Press, 1992.[17℄ F. Ethridge, L. Greengard, A new fast-multipole a elerated Poisson solver intwo dimensions, SIAM J. S i. Comput. 23 (3) (2001) 741(cid:21)760 (ele troni ).[18℄ G. Beylkin, M. J. Mohlenkamp, Algorithms fornumeri al analysis in high dimensions, APPM preprint /x by exponential sums in [1 , ∞ ))