wavelet++ General-Purpose Biorthogonal Wavelet Library in C++ William Garber, Dmitri Volja, Wei Ku

1. wavelet++ General-Purpose Biorthogonal Wavelet Library in C++ William Garber, Dmitri Volja, Wei Ku

2. Introduction to wavelet++ • Why Use Wavelets: • compact support in space x • localized in scale k: • high res detail, low res averages • systematic control of error • sparse representation: • identify, compute with, store only critical data • conservation of moments • interpolating properties • fast O(N) algorithms for • wavelet transform • differential operators • nonlinear operations • products • Applications: • physical problems • biorthogonal bases (bra/ket) • large data sets • high resolution

3. What is a Wavelet Representation? • basis for L2(R) formed by • dilating ψby 2j • translating ψby integer k ψon progressively finer scales φon coarsest scale sJ-3 dJ-3 dJ-2 dJ-1 Coeff in Wavelet Basis

4. Compact Support • compact support in x • localized in k. • Potential V(x) sparse • Kinetic Energy T(k) sparse simultaneously

5. Biorthogonality: Primal Basis Determines Dual Basis Uniquely completeness: biorthogonality: Duals of wavelets are also wavelets.

6. Interpolating Property Useful for Nonlinear Operations zeros at xk φhas zeros at xk so coefficienk sk equals value of function at gridpoints at this scale.

7. Conservation of Moments • CDF(4,4) wavelet ψ • interpolates polynomials up to x3 • has zero moments up to 3rd order. • Coulomb interaction trivial and efficient. • Only 0th moment of φ contributes. • Acts like small number of point charges.

8. Wavelet library Vector Space library • Data Structures: • Filter • basic convolution • LiftingStep • WaveletDef: • define wavelet coeff h,g • provide transform • WaveletRepDense • WaveletRepSparse • store data • Operations: • forward transform • inverse transform • function composition • product • convert to dense • convert to sparse • Data Structures: • VectorSpaceDense • VectorSpaceSparse • Overlap Matrix for finding Duals • Explicit treatment of crystal • translational symmetry • Bivector • Wrapper associating • WaveletRep with VectorSpace • TranslationalyInvariantMatrix • Operations: • Inherit Wavelet operations • DualConj • AddMult: Matrix Multiply • InnerProduct

9. wavelet++ library is Easy to Use: Two Dimensional Cubic Spline Forward Transform WDEF wav = &cubic_spline; BASIS basis(wav); TinyI extent(512,512); WREP wrep(extent, basis); loadPhoto(wrep, fnamePhotoIn); while(nlev-- > 0) { string fname = "photo”; fname += nlev + ".dat"; wrep.transFwd(1); savePhoto(wrep, fnamePhoto); }

10. wavelet++ library is Easy to Use: Two Dimensional Cubic Spline Forward Transform Original 512x512 Level 4 256x256 Level 3 128x128 Level 2 64x64 Level 1 32x32 Level 0 16x16

11. wavelet++ Library is Flexible: Define Your Own Class of Wavelet typedef WaveletDef<double> WDEF; typedef WaveletDefLiftStep<double> LSTEP; // Haar Wavelet with Lifting Steps WDEF haar(“haar”, 1/sq2, sq2, LSTEP(LS_PREDICT, 1, 1, -1.0), LSTEP(LS_UPDATE, 0, 1, 0.5)); // Daubechies Wavelet as Convolution h = (1+sq3)*sq2/8, (3+sq3)*sq2/8, // Filter coefficients (3-sq3)*sq2/8, (1-sq3)*sq2/8; g = h(3), -h(2), h(1), -h(0); std::vector<LSTEP> v; v[0] = LSTEP(h,g,h,g); // convolution step WDEF daubechies(“daubechies”, 1, 1, v);

12. Vector Space Library is Easy to Use: Algebra, Interface dense or sparse: ip = InnerProduct(v1, v2); ip = InnerProductShift(v1, v2, deltaCell); vz = AddMult(vy, LaplacianMatrix, vx); vz = DualConj(vy, vx); vz = Product(v1, v2, v3); vz.FunctionComp(vx, functionToApply); summary: using blitz++ algebra on blitz::Array base class v1 += v2 + Product(v3, v4, v5) + InnerProduct(v3,v4) * v5 + AddMult(v6, mat, v7);

13. Vector Space Library is Easy To Use: Laplacian Operator // constructors BASIS basis(WAV); BOXS geometry(fnameBox); VECSPACE_SPARSE vecspaces(basisp, geometry); VECSPACE_DENSE vecspaced(basisp, geometry.extent()); BIVEC_SPARSE VEC1(vecspaces, VEC_BRA); BIVEC_SPARSE VEC2(vecspaces, VEC_BRA); BIVEC_DENSE vec1(vecspaced, VEC_BRA); BIVEC_DENSE vec2(vecspaced, VEC_BRA); LAPLACIAN mat(vecspaces); // input data storePolyDenseTopLevel(VEC1, vec1, function); // convert to sparse convertToSparse(VEC1, vec1); // VEC2 += mat * VEC1; AddMult(VEC2, mat, VEC1); // convert to dense convertToDense(VEC2, vec1); // plot string fnameOut = "denseout.dat"; plotBox(fnameOut, vec1);

14. Vector Space Library is Easy To Use: Laplacian Operator Input Geometry Size: fine scale j=2 1024x1024 Timing: 1.4 sec Runs on single core of athlon 64x2 3600+ 1.9GHz 2GB ram Input Data: Gaussians Output: Laplacian

15. Mathematical Introduction to Wavelets

16. Wavelet Transform for One Dimensional Data s J average even odd detail sJ-1 dJ-1 sJ-2 dJ-2 sJ-3 dJ-3 Result: Coefficients of Data in Wavelet Basis sJ-3 dJ-3 dJ-2 dJ-1

17. Adding Levels of Detail

18. Two Scale Relationship Haar Scaling Function Haar Scaling Function Haar Wavelet • compact support • detail added by higher scale: more accuracy

19. Wavelets vs. Fourier Series

20. Multiple Dimensions: Tensor Wavelets two dimensional wavelet two dimensional scaling function Cubic Spline

21. Multiple Dimensions: Nonstandard Form Wavelet Representation φ0 ψ0 ψ1 ψ2 φ0 V V V W V W ψ0 W V WW V W ψ1 W V W W ψ2 W V W W • Nonstandard Form: • Forward Transform X; • Forward Transform Y; • Recur on V V • Advantage: • more sparse • can be used for NS operator multiply • Standard Form: • Forward Transform X; • Forward Transform Y; • Recur on whole row/col • Disadvantage: • mix scales

22. Operator Matrices: Nonstandard Form Matrix Multiply Data = Matrix x Data one dimensional data shown • Advantage: • do not mix scales • progressive refinement • for translationally invariant matrix, • blocks are simple filters: O(N)

23. Evaluation of Nonlinear Operations in Wavelet Basis • Given: • f(x) in Wavelet Basis • Goal: • h(x)=g(f(x)) in Wavelet Basis Without Transform to Fine Scale • Algorithm: • Use interpolating properties to compute hk=g(fk) for V data; • forward transform recursively projecting details to lower levels. • Advantage: • FASTER and less memory needed (sparse) • remains in Wavelet Basis.

24. Refresh Vx Algorithm Representation: Start at Vx (j=0) Vi  Vx Wi Wm + neighboring Wm Inverse transform From Vi to Vx recur to j+1 Wm = Wavelet Detail Data level j Vx = Scaling fn Average Data level j+1

25. Nonlinear Function Composition Algorithm RefreshVx; Vx represents averaged-out f(x) apply nonlinear function g to Vx Starting at fine scale, forward transform Vx to Vm/Wm OVERWRITE Vx on level j-1 with Vm on corresponding domain recur to level j-1 The overwrite retains all significant data from finer scales. forward transform Vx(j); OVERWRITE Vx(j-1) with Vm

26. Timing: Laplacian Operator • Sparse wavelets faster than Dense • Handles larger problem with Same Memory

27. Conclusion • Why use wavelets? • compact in x, local in k • conservation of moments • interpolating properties • Why use wavelet++ ? • implement general operators • implement general classes of wavelets • easy to use; flexible; object oriented; templates • Algorithmic Advantages • O(N) algorithms • sparsity saves storage space, computational time • systematic refinement • Benchmarking • speed increase corresponds to degree of sparsity • sparsity greater at finer scales • Feedback welcome • wgarber@bnl.gov, weiku@bnl.gov, dvolja@bnl.gov

28. Coefficients at Finer Scales Add Detail Wavelets in Dual Basis: compact support, zero mean. Wavelet coefficients determine differences or details as wavelet series expansion adds finer accuracy. Compact support means coefficients dk represent a finite region in space at a specific scale:

29. Multiple Dimensions: Nonstandard Form Wavelet Transform φ0 ψ0 ψ1 ψ2 φ0 V V V W V W ψ0 W V WW V W ψ1 W V W W ψ2 W V W W • Standard Form: • Forward Transform X; • Forward Transform Y; • Recur on whole row/col • Disadvantage: • mix scales • Nonstandard Form: • Forward Transform X; • Forward Transform Y; • Recur on V V • Advantage: • more sparse • can be used for NS operator multiply

30. Interpolating Properties Are Useful: Multipole Expansion

31. The efficiency and compactness of wavelets makes them an ideal basis to represent information of general structure. We present the basic concept, implementation and results of wavelet++, our biorthogonal wavelet library in C++. The use of template and metaprogramming technique via blitz++ array library allows simple user-defined wavelet types, and ensures the flexibility and efficiency of this general purpose library. • biorthogonal and orthogonal wavelets • lifting or convolution and conversion • primal or dual and conversion • library of predefined wavelets • built on blitz++ template arrays • operators are O(N) convolutions • evaluation of nonlinear functions • all calculations in wavelet basis not fine scale

32. Goal: Represent Terms in Kohn-Sham Equation

33. Convolution

34. WaveletDefFilter: Basic Convolution blitz stencils for speed UPDATE PREDICT

35. Lifting Haar wavelet

36. CDF(4,4) ScalingFunctions CDF(4,4) Wavelets multiple scales shown

37. Transformed data s(0,k), d(0,k),d(1,k)…,d(J,k) Sampled Function high freq signal low freq signal peak signal High Freq Signal d(k) Peak Signal d(k) Low Freq Signal s(k)

38. Tensor Wavelets two dimensional wavelet two dimensional scaling function Cubic Spline

39. Operator Matrix in Nonstandard Form Advantage: Do Not Mix Scales Blocks on j>0 obtained from V V on j=0 (only a few numbers)

40. Sparse Representation: Small Boxes gui for entry of geometry small boxes centered around atomic cores. boxes form tree; level in tree is level j of wavelet data

41. Timing: Laplacian Operator

42. blitz++ optimized matrix algebra

43. Two Dimensional Cubic Spline Image Averaging Original 512x512 Averaged 256x256 Averaged 128x128 Averaged 64x64 Averaged 32x32 Averaged 16x16

44. Lifting CDF(2,2) wavelet Automatic generation of duals

45. Generators for Matrices by Beylkin’s Method

46. Implementation: WaveletRepDense

48. Error of Wavelet Representation

49. Wavelet Library Implementation Represent Data: WaveletRepDense multidimensional blitz++ array base class. Invoke transformations using WDEF WaveletRepSparse one-dimensional blitz++ array base class. multidimensional boxes referernce subsets of base array. This enables blitz optimized algebra WREP += WREP*WREP + WREP… Define Transformations: WaveletDefFilter Basic filter (convolution) class WaveletDefLiftingStep Predict or Update WaveletDef Convolution (h,g) or Lifting Steps WaveletInstance.cpp Global WDEF library

50. Implementation: WaveletRepSparse gui for entry of geometry small boxes centered around atomic cores. boxes form tree; level in tree is level j of wavelet data