# Presentation: 576263

## wavelet++

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

## 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

## What is a Wavelet Representation?

- basis for L2(R) formed by
-  dilating ψ  by 2j
-   translating ψ  by integer k

- φ on coarsest scale

- ψ on progressively finer scales

- Coeff in Wavelet Basis

- sJ-3

- dJ-3

- dJ-2

- dJ-1

## Compact Support

-  compact support in x
-  localized in k.
-  Potential V(x) sparse
-  Kinetic Energy T(k) sparse

- simultaneously

## Biorthogonality:  Primal Basis Determines Dual Basis Uniquely

- completeness:
- biorthogonality:

- Duals of wavelets are also wavelets.

## Interpolating Property Useful for Nonlinear Operations

- φ  has zeros at xk  so
- coefficienk sk equals value of function at gridpoints at this scale.

- zeros at xk

## 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.

## 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

## 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);
- }

## Original 512x512

- Level 4   256x256

- Level 3   128x128

- Level 2   64x64

- Level 1   32x32

- Level 0  16x16

- wavelet++  library is Easy to Use:
- Two Dimensional Cubic Spline Forward Transform

## 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);

## 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);

## 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);

## Vector Space Library is Easy To Use:  

- Laplacian Operator

- Input Data:  Gaussians

- Output:  Laplacian

- Input Geometry

- Size:         fine scale j=2  1024x1024
- Timing:     1.4 sec
- Runs on single core of
- athlon 64x2 3600+   1.9GHz    2GB ram

## Mathematical Introduction to Wavelets

## Wavelet Transform for One Dimensional Data

- s J

- sJ-1

- dJ-1

- even

- odd

- sJ-3

- dJ-3

- sJ-2

- dJ-2

- average

- detail

- Result:  Coefficients of Data in Wavelet Basis

- sJ-3

- dJ-3

- dJ-2

- dJ-1

## Adding Levels of Detail

## Two Scale Relationship

- Haar Scaling Function

- Haar Scaling Function

- Haar Wavelet

-  compact support
-  detail added by higher scale:  more accuracy

## Wavelets             vs.         Fourier Series

## Cubic Spline

- two dimensional
- scaling function

- two dimensional
- wavelet

- Multiple Dimensions:  Tensor Wavelets

## Multiple Dimensions:  

- Nonstandard Form Wavelet Representation

- Nonstandard Form:
- Forward Transform X;
- Forward Transform Y;
- Recur on V V
- Advantage:
-  more sparse
-  can be used for NS operator multiply

- W V

- W W

- V W

- W V

- W W

- V W

- V V

- W V

- WW

-  V W

- Standard Form:
- Forward Transform X;
- Forward Transform Y;
- Recur on whole row/col
- Disadvantage:
-  mix scales

- ψ0

- φ0

- ψ2

- ψ1

- φ0

- ψ0

- ψ1

- ψ2

## Operator Matrices: Nonstandard Form Matrix Multiply

- Advantage:
-  do not mix scales
-  progressive refinement
-  for translationally invariant matrix,
-     blocks are simple filters:  O(N)

- Data  =  Matrix  x  Data
- one dimensional data shown

## 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.

## Refresh Vx Algorithm

- 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

## Nonlinear Function Composition Algorithm

-       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

## Timing:  Laplacian Operator

| 
j | 
Size, m x m | Time, sec |  | Speed = m2 /T  (106/s) |  | Speed,
Sparse/
Dense |
| --- | --- | --- | --- | --- | --- | --- |
|  |  | Dense | Sparse | Dense | Sparse |  |
| 2 | 1024x1024 | 2.9 | 1.4 | 0.36 | 0.75 | 2.1 x faster |
| 3 | 2048x2048 | 11.5 | 4.1 | 0.36 | 1.02 | 2.8 x faster |
| 4 | 4096x4096 | 83 | 21.5 | 0.20 | 0.78 | 3.9 x faster |
| 5 | 8192x8192 | 837 (swaps) | 98 | 0.080 | 0.68 | 8.5 x faster |

-  Sparse wavelets faster than Dense
-  Handles larger problem with Same Memory

## 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

## 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:

## Multiple Dimensions:  

- Nonstandard Form Wavelet Transform

- Nonstandard Form:
- Forward Transform X;
- Forward Transform Y;
- Recur on V V
- Advantage:
-  more sparse
-  can be used for NS operator multiply

- W V

- W W

- V W

- W V

- W W

- V W

- V V

- W V

- WW

-  V W

- ψ0

- φ0

- ψ2

- ψ1

- φ0

- ψ0

- ψ1

- ψ2

- Standard Form:
- Forward Transform X;
- Forward Transform Y;
- Recur on whole row/col
- Disadvantage:
-  mix scales

## Interpolating Properties Are Useful:

- Multipole Expansion

##  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

- 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.

## Goal:  Represent Terms in Kohn-Sham Equation

## Convolution

## WaveletDefFilter:  Basic Convolution

- blitz stencils for speed

- UPDATE

- PREDICT

## Lifting

- Haar wavelet

## CDF(4,4) Scaling Functions

- CDF(4,4) Scaling Functions

- CDF(4,4) Wavelets

- multiple scales shown

## Sampled Function

- Transformed data
- s(0,k), d(0,k),d(1,k)...,d(J,k)

- High Freq Signal d(k)

- Peak Signal d(k)

- Low Freq Signal s(k)

- peak signal

- low freq signal

- high freq signal

## Cubic Spline

- two dimensional
- scaling function

- two dimensional
- wavelet

- Tensor Wavelets

## Operator Matrix in Nonstandard Form

- Blocks on j>0 obtained from V V on j=0
- (only a few numbers)

- Advantage:  Do Not Mix Scales

## 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

## Timing:  Laplacian Operator

|  |  | Time Dense, sec |  | Speed Dense,
 n2 /T  (106/s) | Time Sparse, sec
 |  | Speed Sparse,
n2 /T (106/s) | Speed,
Sp/De |
| --- | --- | --- | --- | --- | --- | --- | --- | --- |
| j | Size, nxn | 1st iter | 2nd iter | both iter | 1st iter | 2nd iter | both iter |  |
| 2 | 1024x1024 | 2.8 | 0.11 | 0.36 | 0.70 | 0.70 | 0.75 | 2.1 x |
| 3 | 2048x2048 | 11.1 | 0.37 | 0.36 | 2.1 | 2.0 | 1.02 | 2.8 x |
| 4 | 4096x4096 | 82 | 1.32 | 0.20 | 14.9 | 6.6 | 0.78 | 3.9 x |
| 5 | 8192x8192 | 511 | 326 (swap) | 0.080 | 68 | 30. | 0.68 | 8.5 x |

## blitz++ optimized matrix algebra

## Two Dimensional Cubic Spline Image Averaging

- Original 512x512

- Averaged 256x256

- Averaged 128x128

- Averaged 64x64

- Averaged 32x32

- Averaged 16x16

## Lifting

- CDF(2,2) wavelet

- Automatic generation of duals

## Generators for Matrices          by Beylkin’s Method 

## Implementation:  WaveletRepDense

## Slide 47

## Error of Wavelet Representation

## 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

## 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

## Implementation:  WaveletRepSparse

## WaveletDef:  Define Transformation

- transFwd and transInv

- LSTEP 0

- LSTEP N-1

- LSTEP SCALED

- LSTEP SCALES

- ... etc ...

- transFwd

- transInv

## Wavelet Constructors:  Lifting Steps

-   typedef WaveletDef<double>             WDEF;
-   const double ScaleS = 1.0;
-   const double ScaleD = 1.0;
-   // WaveletDef
-   WDEF CDF44("CDF_4_4_CUSTOM", ScaleS, ScaleD,
- 	           LSTEP(LS_PREDICT, 2, sq5,  1/16.,     -9/16.,     -9/16.,      1/16.),
- 	           LSTEP(LS_UPDATE,  2, sq5, -1/sq232, 9/sq232, 9/sq232, -1/sq232));

##    // single convolution lifting step

-   blitz::Array<double,1> h(blitz::Range(0,3));
-   blitz::Array<double,1> g(blitz::Range(-2,1));
-    h =  (1-sq3)/(4*sq2),
-           (3-sq3)/(4*sq2),
-           (3+sq3)/(4*sq2),
-           (1+sq3)/(4*sq2);
-     // g(k) = (-1)^k  h(1-k)
-     g = h(3),
-          -h(2),
-           h(1),
-          -h(0);
-    // single convolution lifting step
-     std::vector<LSTEP> v(1);
-     v[0] = LSTEP(h, g, h, g);
-     WDEF Daubechies4a("Daubechies_4a_CUSTOM", 1.0, 1.0, v);

- Wavelet Constructors:  Convolutions

## Vector Space Library

- Major classes:
-  VectorSpace
-  Bivector
-  TransInvMat2
-  OuterProduct

- GlobalFunctions
-  InnerProduct
-  DualConj
-  ConvertToSparse
-  ConvertToDense
-  AddMult
-  Product

- dense and sparse versions

- Vector Space Data:
-  mask for InnerProduct
-  Overlap Matrix for DualConj

## Vector Space Library

- Constructors

- dense

- sparse

## Dmitri Volja, William Garber, Wei Ku
Condensed Matter Physics & Materials Science Department, BNL
*Physics Dept. SUNY @ Stony Brook

*Condensed Matter Physics & Materials Science Department, BNL*
**Physics Dept. SUNY @ Stony Brook*

- January 17, 2007

- Minimization scheme for electronic structure calculations within
- Bi-orthogonal wavelet basis

## Outline

-  Current framework of DFT tools, limitations
-  Minimization approach
-  Applications
-  Basis Importance
-  Examples

## Current framework

**DFT :**

**Energy of interacting system is unique functional of density :**

**   ****Eigenvalue problem at each iteration  :  O(N****3****): ****Expensive!**

## Additional problems

-  Need for certain approximations for exchange-correlation
-  Difficulty in introduction of new functionals
-   Orbital dependent functionals: OEP, HF, OPM
-  Kohn-Sham non-interacting approximation
-   Need for density matrix functional theory
-  Broken symmetries (orbital, charge, inversion)
-   are not treated properly

## Solution

-  Direct functional minimization without KS equations

- Still iterative scheme: Basic steps
-  Acting on the state with Hamiltonian O(NM2)
-   Orthonormalization of states O(N2M)
- N-number of states, M basis size, N<<M

- * William Garber and Wei Ku

- Linear scaling for compact basis O(M)
-  More stable than iterative Kohn-Sham scheme in minima search

## Algorithm of CG minimization

- Symmetric Orthonormalization:

**O(NM****2****)**

**  ****O(N****2****M)**

-  Resemblance to starting states
-  Order independence

## Application to non-local potentials:Minimization w.r.t. truncated DM

- * T.L. Gilbert, Phys. Rev. B **12** 2111 (1975)   ** W. Kohn, Phys. Rev. **115**, 809 (1959)

- 1)  Generalization of Kohn-Sham theory to non-local potentials
- 2)  Beyond Kohn-Sham: exact Ekin
- 3)  Exponential decay of density matrix in insulators, semiconductors.
- 4)  Single parameter :Rc cutoff: O(N)
- 5)  For insulators g.s. solution are Wannier states

- Any ground-state observable can be written in terms of 1-RDM

- Need for smart basis choice !

-   In wavelet basis DM is sparse even for itinerant (Metal) systems

## Basis requirements  

-  Simultaneous description of short/long scales
- - Core/valence, localized/delocalized, smooth/sharp etc.
-  Compactness, Large system size, order O(N) scaling
-  Systematic improvement of numerical accuracy
-  Elimination of basis error

- Possible Solution:  Wavelet basis !

## Scaling functions and wavelets of primal and dual spaces

- Average of signal: **V** space

- Details of signal:  **W** space

## Basis properties

-  Compactness in Real space

-  Compactness in Fourier space

-  Vanishing moment property

-  Multiresolution

-  Interpolating property

## Generic Algorithm

- CG

- Functional

- Constraint

- Boundary

-  Information on the energy functional used
-  Easy implementation of new fucntionals

- Convergence

-  Information on the constraints used
-  Lagrange matrix
-  apply()

-  Information on the boundaries within unit cell
-  Information on the crystal periodicity

-  Information on the convergence criteria

-  get_gradient()
-  get_dE_2nd_order_corr()

-  modify_gradient()

-  apply()

-  apply()

**Supporting Objects**

## CG: Actual scheme 

-     constraint.apply( states);                        //  **modify states according constriants**
-     functional.get_gradient( gs, states);       // **get gradient of search**
-     constraint.modify_gradient( gs, states); // **add Lagrange multipliers**
-     boundary.apply( gs);                              // **apply boundary condition**
-     small_vec::iterator norm2p( norm2_gs.begin());
-  for( St_iterator gp = gs.begin(), hp = hs.begin(); gp != gs.end();gp++,hp++,norm2p++) {
-       const double norm2_g( norm2( *gp));
-     CG_MINI::mult_subtract_assign( *hp, (norm2_g)/(*norm2p) , (*gp));
-       (*norm2p)  = norm2_g;
-     }                                                               // **construct true direction of search**

-  functional.get_dE_2nd_order_corr( states, hs, gs, constraint,  dE_2nd);
-                                                                     // **perform 1D minimization**
-      small_vec::iterator  dEp = dE_2nd.begin();
-  for( St_iterator sp = states.begin(), hp = hs.begin(); sp != states.end(); sp++,  hp++,dEp++)
-           { CG_MINI::add_mult( *sp, *hp, *dEp );}
-                                                                      // **actual step to a new position**
-     cycle++;
-   **}  //end of while loop**

- CG: Actual scheme

## One electron under various boundary conditions

- b

- a

- c

- d

## Two electrons in a nano-wire cross-section

## Summary: part 1 

-  Possibility to describe broken symmetry states
-  High mobility in introduction of orbital functionals
- -HF, EXX , OEP etc.
-  Can go beyond Kohn-Sham:
- -Density Matrix Functional Theory
-  Order O(N) calculations