slide1 l.
Download
Skip this Video
Loading SlideShow in 5 Seconds..
wavelet++ General-Purpose Biorthogonal Wavelet Library in C++ William Garber, Dmitri Volja, Wei Ku PowerPoint Presentation
Download Presentation
wavelet++ General-Purpose Biorthogonal Wavelet Library in C++ William Garber, Dmitri Volja, Wei Ku

Loading in 2 Seconds...

play fullscreen
1 / 71

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


  • 293 Views
  • Uploaded on

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

loader
I am the owner, or an agent authorized to act on behalf of the owner, of the copyrighted work described.
capcha
Download Presentation

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


An Image/Link below is provided (as is) to download presentation

Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author.While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server.


- - - - - - - - - - - - - - - - - - - - - - - - - - E N D - - - - - - - - - - - - - - - - - - - - - - - - - -
Presentation Transcript
slide1

wavelet++

General-Purpose Biorthogonal Wavelet Library in C++

William Garber, Dmitri Volja, Wei Ku

slide2

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
slide3

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

slide4

Compact Support

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

simultaneously

slide5

Biorthogonality: Primal Basis Determines Dual Basis Uniquely

completeness:

biorthogonality:

Duals of wavelets are also wavelets.

slide6

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.

slide7

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

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
slide9

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

}

slide10

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

slide11

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

slide12

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

slide13

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

slide14

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

slide16

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

slide18

Two Scale Relationship

Haar Scaling Function

Haar Scaling Function

Haar Wavelet

  • compact support
  • detail added by higher scale: more accuracy
slide20

Multiple Dimensions: Tensor Wavelets

two dimensional

wavelet

two dimensional

scaling function

Cubic Spline

slide21

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
slide22

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)
slide23

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

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

slide25

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

slide26

Timing: Laplacian Operator

  • Sparse wavelets faster than Dense
  • Handles larger problem with Same Memory
slide27

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
slide28

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:

slide29

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
slide31

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
slide34

WaveletDefFilter: Basic Convolution

blitz stencils for speed

UPDATE

PREDICT

slide35

Lifting

Haar wavelet

slide36

CDF(4,4) ScalingFunctions

CDF(4,4) Wavelets

multiple scales shown

slide37

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)

slide38

Tensor Wavelets

two dimensional

wavelet

two dimensional

scaling function

Cubic Spline

slide39

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)

slide40

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

slide43

Two Dimensional Cubic Spline Image Averaging

Original 512x512

Averaged 256x256

Averaged 128x128

Averaged 64x64

Averaged 32x32

Averaged 16x16

slide44

Lifting

CDF(2,2) wavelet

Automatic generation of duals

slide49

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

slide50

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

slide52

WaveletDef: Define Transformation

transFwd and transInv

transFwd

LSTEP 0

... etc ...

LSTEP N-1

LSTEP SCALES

LSTEP SCALED

transInv

slide53

Wavelet Constructors: Lifting Steps

typedef WaveletDefLiftStep<double> LSTEP;

typedef WaveletDef<double> WDEF;

const double ScaleS = 1.0;

const double ScaleD = 1.0;

// WaveletDef

WDEFCDF44("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));

slide54

Wavelet Constructors: Convolutions

// construct wavelet with vector of lifting steps

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 convolutionlifting step

std::vector<LSTEP> v(1);

v[0] = LSTEP(h, g, h, g);

WDEF Daubechies4a("Daubechies_4a_CUSTOM", 1.0, 1.0, v);

slide55

Vector Space Library

dense and sparse versions

  • Major classes:
  • VectorSpace
  • Bivector
  • TransInvMat2
  • OuterProduct
  • GlobalFunctions
  • InnerProduct
  • DualConj
  • ConvertToSparse
  • ConvertToDense
  • AddMult
  • Product
  • Vector Space Data:
  • mask for InnerProduct
  • Overlap Matrix for DualConj
slide56

Vector Space Library

Constructors

dense

sparse

slide57

Minimization scheme for electronic structure calculations within

Bi-orthogonal wavelet basis

Dmitri Volja, William Garber, Wei Ku

Condensed Matter Physics & Materials Science Department, BNL

*Physics Dept. SUNY @ Stony Brook

January 17, 2007

outline
Outline
  • Current framework of DFT tools, limitations
  • Minimization approach
  • Applications
  • Basis Importance
  • Examples
slide59

Kohn-Sham scheme :

Auxiliary orbitals of non-interacting electrons

Current framework

Energy of interacting system is unique functional of density :

DFT :

Eigenvalue problem at each iteration : O(N3): Expensive!

additional problems
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
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

Linear scaling for compact basis O(M)

More stable than iterative Kohn-Sham scheme in minima search

* William Garber and Wei Ku

algorithm of cg minimization
Algorithm of CG minimization

O(NM2)

  • Resemblance to starting states
  • Order independence

O(N2M)

Symmetric Orthonormalization:

application to non local potentials minimization w r t truncated dm
Application to non-local potentials:Minimization w.r.t. truncated DM

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

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

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

Need for smart basis choice !

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

slide64

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 !

slide65

Scaling functions and wavelets of primal and dual spaces

Average of signal: V space

Details of signal: W space

basis properties
Basis properties
  • Compactness in Real space
  • Multiresolution
  • Compactness in Fourier space
  • Vanishing moment property
  • Interpolating property
generic algorithm
Generic Algorithm

CG

Supporting Objects

Functional

Constraint

Boundary

Convergence

  • Information on the constraints used
  • Lagrange matrix
  • apply()
  • Information on the boundaries within unit cell
  • Information on the crystal periodicity
  • Information on the energy functional used
  • Easy implementation of new fucntionals
  • Information on the convergence criteria
  • get_gradient()
  • get_dE_2nd_order_corr()
  • apply()
  • apply()
  • modify_gradient()
slide68

CG: Actual scheme

while( ! convergence.apply( hs) ) {// main loop

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

slide69

y

x

One electron under various boundary conditions

a

b

c

d

slide70

y

x

Two electrons in a nano-wire cross-section

slide71

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