Draft Lattice Builder A Software Tool for Constructing Lattice - - PowerPoint PPT Presentation

▶
draft
SMART_READER_LITE
LIVE PREVIEW

Draft Lattice Builder A Software Tool for Constructing Lattice - - PowerPoint PPT Presentation

Jan 20, 2023 231 likes •929 views

1 Draft Lattice Builder A Software Tool for Constructing Lattice Rules Pierre LEcuyer and David Munger Informatique et Recherche Op erationnelle, Universit e de Montr eal 1 Draft Lattice Builder A Software Tool for Constructing

slide-1
SLIDE 1

Draft

1

Lattice Builder

A Software Tool for Constructing Lattice Rules

Pierre L’Ecuyer and David Munger

Informatique et Recherche Op´ erationnelle, Universit´ e de Montr´ eal

slide-2
SLIDE 2

Draft

1

Lattice Builder

A Software Tool for Constructing Lattice Rules

Pierre L’Ecuyer and David Munger

Informatique et Recherche Op´ erationnelle, Universit´ e de Montr´ eal Outline:

  • 1. Setting: integration by (randomly-shifted) lattice rules.
  • 2. Want a tool to search for good lattices adapted to problem at hand,

for arbitrary dimension, size, figure of merit, ...

  • 3. Selection criteria, construction methods, weights on projections, ...
  • 4. Introducing Lattice Builder.
slide-3
SLIDE 3

Draft

2

Monte Carlo integration

Want to estimate µ =

  • [0,1)s f (u) du = E[f (U)]

where f : [0, 1)s → R and U is a uniform r.v. over [0, 1)s. Standard Monte Carlo:

◮ Generate n independent realizations of U, say U1, . . . , Un; ◮ estimate µ by ˆ

µn = 1

n

n

i=1 f (Ui).
slide-4
SLIDE 4

Draft

2

Monte Carlo integration

Want to estimate µ =

  • [0,1)s f (u) du = E[f (U)]

where f : [0, 1)s → R and U is a uniform r.v. over [0, 1)s. Standard Monte Carlo:

◮ Generate n independent realizations of U, say U1, . . . , Un; ◮ estimate µ by ˆ

µn = 1

n

n

i=1 f (Ui).

Variance: Var[ˆ µn] = σ2/n where σ2 =

  • [0,1)s f 2(u) du − µ.

Central limit theorem: √n(ˆ µn − µ)/σ ⇒ N(0, 1) when n → ∞.

slide-5
SLIDE 5

Draft

3

Quasi-Monte Carlo (QMC)

Replace the random points Ui by a set of deterministic points Pn = {u0, . . . , un−1} that cover [0, 1)s more evenly. That is, low discrepancy between empirical distribution of Pn and uniform distribution.

slide-6
SLIDE 6

Draft

3

Quasi-Monte Carlo (QMC)

Replace the random points Ui by a set of deterministic points Pn = {u0, . . . , un−1} that cover [0, 1)s more evenly. That is, low discrepancy between empirical distribution of Pn and uniform distribution. Main construction methods: lattice rules and digital nets (Korobov, Hammersley, Halton, Sobol’, Faure, Niederreiter, etc.)

slide-7
SLIDE 7

Draft

4

Randomized quasi-Monte Carlo (RQMC) estimator

ˆ µn,rqmc = 1 n

n−1
  • i=0

f (Ui), with Pn = {U0, . . . , Un−1} ⊂ (0, 1)s an RQMC point set: (i) each point Ui has the uniform distribution over (0, 1)s; (ii) Pn as a whole is a low-discrepancy point set. E[ˆ µn,rqmc] = µ (unbiased).

Var[ˆ µn,rqmc] = Var[f (Ui)] n + 2 n2
  • i<j
Cov[f (Ui), f (Uj)].

We want to make the last sum as negative as possible.

slide-8
SLIDE 8

Draft

5

Integration lattice

Ls =   v =

s
  • j=1

zjvj such that each zj ∈ Z    , where v1, . . . , vs ∈ Rs are linearly independent over R and where Ls contains Zs. Lattice rule: Take Pn = {u0, . . . , un−1} = Ls ∩ [0, 1)s.

slide-9
SLIDE 9

Draft

5

Integration lattice

Ls =   v =

s
  • j=1

zjvj such that each zj ∈ Z    , where v1, . . . , vs ∈ Rs are linearly independent over R and where Ls contains Zs. Lattice rule: Take Pn = {u0, . . . , un−1} = Ls ∩ [0, 1)s. Lattice rule of rank 1: ui = iv1 mod 1 for i = 0, . . . , n − 1. nv1 = a = (a1, . . . , as) ∈ Zs

n.

Korobov rule: a = (1, a, a2 mod n, . . . , as−1 mod n).

slide-10
SLIDE 10

Draft

5

Integration lattice

Ls =   v =

s
  • j=1

zjvj such that each zj ∈ Z    , where v1, . . . , vs ∈ Rs are linearly independent over R and where Ls contains Zs. Lattice rule: Take Pn = {u0, . . . , un−1} = Ls ∩ [0, 1)s. Lattice rule of rank 1: ui = iv1 mod 1 for i = 0, . . . , n − 1. nv1 = a = (a1, . . . , as) ∈ Zs

n.

Korobov rule: a = (1, a, a2 mod n, . . . , as−1 mod n). For any u ⊂ {1, . . . , s}, the projection Ls(u) of Ls is also a lattice, with point set Pn(u).

slide-11
SLIDE 11

Draft

5

Integration lattice

Ls =   v =

s
  • j=1

zjvj such that each zj ∈ Z    , where v1, . . . , vs ∈ Rs are linearly independent over R and where Ls contains Zs. Lattice rule: Take Pn = {u0, . . . , un−1} = Ls ∩ [0, 1)s. Lattice rule of rank 1: ui = iv1 mod 1 for i = 0, . . . , n − 1. nv1 = a = (a1, . . . , as) ∈ Zs

n.

Korobov rule: a = (1, a, a2 mod n, . . . , as−1 mod n). For any u ⊂ {1, . . . , s}, the projection Ls(u) of Ls is also a lattice, with point set Pn(u). Random shift modulo 1: generate a single point U uniformly over (0, 1)s and add it to each point of Pn, modulo 1, coordinate-wise: Ui = (ui + U) mod 1. Each Ui is uniformly distributed over [0, 1)s.

slide-12
SLIDE 12

Draft

6

Randomly-shifted lattice: Two-dim. example

1 1 ui,1 ui,2

slide-13
SLIDE 13

Draft

6

Randomly-shifted lattice: Two-dim. example

1 1 ui,1 ui,2

slide-14
SLIDE 14

Draft

6

Randomly-shifted lattice: Two-dim. example

1 1 ui,1 ui,2 U

slide-15
SLIDE 15

Draft

6

Randomly-shifted lattice: Two-dim. example

1 1 ui,1 ui,2

slide-16
SLIDE 16

Draft

6

Randomly-shifted lattice: Two-dim. example

1 1 ui,1 ui,2

slide-17
SLIDE 17

Draft

7

Example of a poor lattice: s = 2, n = 101, a = 51

1 1 ui,1 ui,2 Good uniformity in one dimension, but not in two!

slide-18
SLIDE 18

Draft

8

Variance expression

Suppose f has Fourier expansion f (u) =

  • h∈Zs

ˆ f (h)e2π√−1htu. For a randomly shifted lattice, the exact variance is Var[ˆ µn,rqmc] =

  • 0=h∈L∗
s

|ˆ f (h)|2, where L∗

s = {h ∈ Rs : htv ∈ Z for all v ∈ Ls} ⊆ Zs is the dual lattice.

From the viewpoint of variance reduction, an optimal lattice for f minimizes D2

f (Pn) = Var[ˆ

µn,rqmc]. But not a practical criterion...

slide-19
SLIDE 19

Draft

9

Periodic smooth functions

If f has square-integrable mixed partial derivatives up to order α/2, and the periodic continuations of its derivatives up to order α/2 − 1 are continuous across the unit cube boundaries, then |ˆ f (h)|2 = O((max(1, |h1|) · · · max(1, |hs|))−α). Moreover, there is a vector v1 = v1(n) such that Pα

def

=

  • 0=h∈L∗
s

(max(1, |h1|) · · · max(1, |hs|))−α = O(n−α+δ). This Pα has been proposed long ago as a figure of merit, often with α = 2. It is the variance for a worst-case f having |ˆ f (h)|2 = (max(1, |h1|) · · · max(1, |hs|))−α.

slide-20
SLIDE 20

Draft

10

If α/2 is a positive integer, this worst-case f is f (u) =

  • u⊆{1,...,s}
  • j∈u

(2π)α/2 (α/2)! Bα/2(uj). where Bα/2 is the Bernoulli polynomial of degree α/2. This worst-case function is not necessarily representative of what happens in applications. Moreover, the hidden factor in O increases quickly with s, so this result is not very useful for large s. To get a bound that is uniform in s, the Fourier coefficients must decrease faster with the dimension and “size” of vectors h; that is, f must be “smoother” in high-dimensional projections. This is typically what happens in applications where RQMC is really effective! [cf., tractability theory; W´

  • zniakowski, Novak, Sloan, Dick, etc.]
slide-21
SLIDE 21

Draft

11

A very general weighted Pα

Pα can be generalized by giving different weights w(h) to the vectors h: ˜ Pα

def

=

  • 0=h∈L∗
s

w(h)(max(1, |h1|) · · · max(1, |hs|))−α. But how do we choose these weights? There are too many! The optimal weights to minimize the variance are: w(h) = (max(1, |h1|) · · · max(1, |hs|))α|ˆ f (h)|2.

slide-22
SLIDE 22

Draft

12

ANOVA decomposition

A cruder expansion of f (u) = f (u1, . . . , us): f (u) =

  • u⊆{1,...,s}

fu(u) = µ +

s
  • i=1

f{i}(ui) +

s
  • i,j=1

f{i,j}(ui, uj) + · · · where fu(u) =

  • [0,1)|¯
u| f (u) du¯ u −
  • v⊂u

fv(uv). The Monte Carlo variance decomposes as σ2 =

  • u⊆{1,...,s}

σ2

u,

where σ2

u = Var[fu(U)].

The σ2

u’s can be estimated by MC or RQMC.

Heuristic intuition: Make sure the projections Pn(u) are very uniform for the important subsets u (i.e., with larger σ2

u).
slide-23
SLIDE 23

Draft

13

Weighted Pγ,α with projection-dependent weights γu

Denote u(h) = u(h1, . . . , hs) the set of indices j for which hj = 0. Pγ,α =
  • 0=h∈L∗
s γu(h)(max(1, |h1|) · · · max(1, |hs|))−α. For α/2 integer > 0, if ui = (ui,1, . . . , ui,s) = iv1 mod 1, Pγ,α =
  • ∅=u⊆{1,...,s}
1 n n−1
  • i=0
γu −(−4π2)α/2 (α)! |u| j∈u Bα(ui,j) (finite sum) and the corresponding variation is V 2 γ (f ) =
  • ∅=u⊆{1,...,s}
1 γu(4π2)α|u|/2
  • [0,1]|u|
  • ∂α|u|/2
∂uα/2 fu(u)
  • 2
du, for f : [0, 1)s → R smooth enough. Then, Var[ˆ µn,rqmc] =
  • u⊆{1,...,s}
Var[ˆ µn,rqmc(fu)] ≤ V 2 γ (f )Pγ,α .
slide-24
SLIDE 24

Draft

14

This Pγ,α is the RQMC variance for the worst-case function f ∗(u) =

  • u⊆{1,...,s}

√γu

  • j∈u

(2π)α/2 (α/2)! Bα/2(uj), whose square Fourier coefficients are |ˆ f ∗(h)|2 = γu(h)(max(1, |h1|) · · · max(1, |hs|))−α. For this function, we have

σ2 u = γu
  • Var[Bα/2(U)] (4π2)α/2
((α/2)!)2 |u| = γu
  • |Bα(0)|(4π2)α/2
(α)! |u| .

For α = 2, this gives γu = (3/π2)|u|σ2

u ≈ (0.30396)|u|σ2 u.

For α = 4, this gives γu = [45/π4]|u|σ2

u ≈ (0.46197)|u|σ2 u.

For α → ∞, we have γu → (0.5)|u|σ2

u.

Note: The correct weights are not proportional to the variances σ2

u.
slide-25
SLIDE 25

Draft

15

Heuristics for choosing the weights

For f ∗, we should take γu = ρ|u|σ2

u for some constant ρ.

But there are still 2s − 1 subsets u to consider! One could define a simple parametric model for the square variations and then estimate the parameters by matching the ANOVA variances σ2

u (e.g.,

Wang and Sloan 2006, L. and Munger 2012). For example, product weights: γu =

j∈u γj for some constants γj ≥ 0.

Order-dependent weights: γu depends only on |u|. Example: γu = 1 for |u| ≤ d and γu = 0 otherwise. Wang (2007) suggests this with d = 2. Mixture: POD weights (Kuo et al. 2011). Note that all one-dimensional projections (before random shift) are the

  • same. So the weights γu for |u| = 1 are irrelevant.
slide-26
SLIDE 26

Draft

16

Weighted Rγ,α

Take D2

u(Pn) = 1

n

n−1
  • i=0
  • j∈u

 

⌊n/2⌋
  • h=−⌊(n−1)/2⌋

max(1, |h|)−αe2πιhui,j − 1   . Upper bounds on Pα can be computed in terms of Rα. Can be computed for any α > 0 (finite sum). We compute it using FFT.

slide-27
SLIDE 27

Draft

17

Figure of merit based on the spectral test

Compute the shortest vector ℓu(Pn) in dual lattice for each projection u and normalize by an upper bound ℓ∗

|u|(n) (with Euclidean length):

Du(Pn) = ℓ∗

|u|(n)

ℓu(Pn) ≥ 1.

slide-28
SLIDE 28

Draft

17

Figure of merit based on the spectral test

Compute the shortest vector ℓu(Pn) in dual lattice for each projection u and normalize by an upper bound ℓ∗

|u|(n) (with Euclidean length):

Du(Pn) = ℓ∗

|u|(n)

ℓu(Pn) ≥ 1.

  • L. and Lemieux (2000), etc., maximize

Mt1,...,td = min   min

2≤r≤t1

ℓ{1,...,r}(Pn) ℓ∗

r (n)

, min

2≤r≤d

min

u={j1,...,jr}⊂{1,...,s} 1=j1<···<jr≤tr

ℓu(Pn) ℓ∗

r (n)

  . Computing time of ℓu(Pn) is almost independent of n, but exponential in |u|. Poor lattices can be eliminated quickly. Can use a different norm, compute shortest vector in primal lattice, etc.

slide-29
SLIDE 29

Draft

18

Search methods

Korobov lattices. Search over all admissible a, for a = (1, a, a2, . . . , ...). Random Korobov. Try r random values of a.

slide-30
SLIDE 30

Draft

18

Search methods

Korobov lattices. Search over all admissible a, for a = (1, a, a2, . . . , ...). Random Korobov. Try r random values of a. Rank 1, exhaustive search. Pure random search. Try admissible vectors a at random.

slide-31
SLIDE 31

Draft

18

Search methods

Korobov lattices. Search over all admissible a, for a = (1, a, a2, . . . , ...). Random Korobov. Try r random values of a. Rank 1, exhaustive search. Pure random search. Try admissible vectors a at random. Component by component (CBC) construction. (Sloan, Kuo, etc.). Let a1 = 1; For j = 2, 3, . . . , s, find z ∈ {1, . . . , n − 1}, gcd(z, n) = 1, such that (a1, a2, . . . , aj = z) minimizes Dq(Pn({1, . . . , j})). Fast CBC construction for Pγ,α: use FFT. (Nuyens, Cools).

slide-32
SLIDE 32

Draft

18

Search methods

Korobov lattices. Search over all admissible a, for a = (1, a, a2, . . . , ...). Random Korobov. Try r random values of a. Rank 1, exhaustive search. Pure random search. Try admissible vectors a at random. Component by component (CBC) construction. (Sloan, Kuo, etc.). Let a1 = 1; For j = 2, 3, . . . , s, find z ∈ {1, . . . , n − 1}, gcd(z, n) = 1, such that (a1, a2, . . . , aj = z) minimizes Dq(Pn({1, . . . , j})). Fast CBC construction for Pγ,α: use FFT. (Nuyens, Cools). Randomized CBC construction. Let a1 = 1; For j = 2, . . . , s, try r random z ∈ {1, . . . , n − 1}, gcd(z, n) = 1, and retain (a1, a2, . . . , aj = z) that minimizes Dq(Pn({1, . . . , j})). Can add filters to eliminate poor lattices more quickly.

slide-33
SLIDE 33

Draft

19

Embedded lattices

Pn1 ⊂ Pn2 ⊂ . . . Pnm with n1 < n2 < · · · < nm, for some m > 0. Usually: nk = bc+k for integers c ≥ 0 and b ≥ 2, typically with b = 2, ak = ak+1 mod nk for all k < m, and the same random shift. We need a measure that accounts for the quality of all m lattices. We standardize the merit at all levels k so they have a comparable scale: Eq(Pn) = Dq(Pn)/D∗

q(n),

where D∗

q(n) is a normalization factor, e.g., a bound on Dq(Pn) or a

bound on its average over all (a1, . . . , as) under consideration. For Pγ,α, bounds by Sinescu and L. (2012) and Dick et al. (2008). For CBC, we do this for each coordinate j = 1, . . . , s (replace s by j). Also used as filters. Then we can take as a global measure (with sum or max): ¯ Eq,m(Pn1, . . . , Pnm) q =

m
  • k=1

wk [Eq(Pnk)]q .

slide-34
SLIDE 34

Draft

20

Existing tools

Construction: Nuyens (2012) provides Matlab code for fast-CBC construction of lattice rules based on Pγ,α, with product and

  • rder-dependent weights.

Precomputed tables for fixed criteria: Maisonneuve (1972), Sloan and Joe (1994), L. and Lemieux (2000), Kuo (2012), etc. Software for using (randomized) lattice rules in simulations is also available in many places (e.g., in SSJ).

slide-35
SLIDE 35

Draft

21

Lattice Builder

Implemented as C++ library, modular object-oriented design, accessible from a program via API. Various choices of figures of merit, arbitrary weights, construction methods, etc. Easily extensible. For better run-time efficiency, uses static polymorphism, via templates, rather than dynamic polymorphism. Several other techniques to reduce computations and improve speed. Offers a pre-compiled program with Unix-like command line interface. Also graphical interface. Available for download on GitHub, with source code, documentation, and precompiled executable codes for Linux or Windows, in 32-bit and 64-bit versions.

Show graphical interface

slide-36
SLIDE 36

Draft

22

Examples of applications

slide-37
SLIDE 37

Draft

23

Example: Playing with the Weights

To see the effect of weights selection on RQMC variance, when choosing a lattice rule, we shall integrate the worst-case function f ∗

α (u) =
  • u⊆{1,...,s}

√vu

  • j∈u

(2π)α/2 (α/2)! Bα/2(uj), whose RQMC variance is Pγ,α. The ideal weights are γu = vu. In our experiments, we measure the inflation factor of RQMC variance when we use a lattice rule constructed via fast CBC with different weights γu = vu. We start with s = 10 and vu = Γ|u| for |u| ≤ k, for a given integer k and a given constant Γ > 0. We select the weights as γu = ˜ Γ|u| for |u| ≤ ˜ k, where ˜ Γ and ˜ k may differ from Γ and k.

slide-38
SLIDE 38

Draft

24

Ratio of RQMC variances for modified vs ideal weights n A1 A2 B1 B2 C1 C2 28 1.11 1.21 1.13 4.08 3.82 6.80 29 1.21 1.10 1.42 10.5 2.93 7.25 210 1.36 1.38 2.04 4.64 2.86 5.94 211 1.24 1.43 2.40 6.18 2.15 5.14 212 1.42 1.66 3.79 13.2 2.47 5.94 213 1.30 2.38 5.51 9.09 2.66 5.97 214 1.51 2.54 30.5 8.66 9.11 29.1 215 1.46 1.93 25.6 13.3 3.52 9.71 216 1.80 2.55 3.13 12.9 2.73 10.2 A1: k = ˜ k = s = 10, Γ2 = 0.1 and ˜ Γ2 = 0.001. Not too bad. A2: k = ˜ k = s = 10, Γ2 = 0.001 and ˜ Γ2 = 0.1. More impact.

slide-39
SLIDE 39

Draft

25

B1: Γ2 = ˜ Γ2 = 0.1, k = 4 and ˜ k = 2. Criterion is blind on projections of

  • rder 3 and 4. This has unpredictable (sometimes dramatic) impact on

RQMC variance. B2: Γ2 = ˜ Γ2 = 0.5, k = 2 and ˜ k = 4. Gives weight to irrelevant

  • projections. Stronger degradation on average.
slide-40
SLIDE 40

Draft

26

C1: Γ2 = ˜ Γ2 = 0.1 and k = ˜ k = 4, but increase the variation of f by replacing v2

u with v2 u + ˜

v2

u , with

˜ v2

u = 1.0 for u = {1, 3}, {3, 5}, {5, 7}, {7, 9},

˜ v2

u = 0.5 for u = {2, 3, 4}, {4, 5, 6}, {6, 7, 8}, {8, 9, 10},

˜ v2

u = 0.25 for u = {1, 2, 3, 4}, {4, 5, 6, 7}, {7, 8, 9, 10},

and ˜ vu = 0 for all other u. Important projections are not given enough weight relative to others. C2: Like C1, but v2

u is replaced with only ˜

v2

u , as defined above. Many

irrelevant projections u, with ˜ vu = 0, now have weights γu > 0. In all cases, the variance ratio increases (non-monotonically) with n.

slide-41
SLIDE 41

Draft

27

Competing Projections for spectral and P2 criteria

Plots of order 2 vs order 3 contributions, for 1000 random lattices with n = 220, with s = 3 (top) and 10 (bottom). 100 101 100 101
  • rder-2
  • rder-3
spectral, s = 3 10−9 10−6 10−3 10−8 10−6 10−4
  • rder-2
  • rder-3
P2, s = 3 101.5 102 102.4 102.6
  • rder-3
spectral, s = 10 10−9 10−6 10−3 10−6 10−4 10−2
  • rder-3
P2, s = 10
slide-42
SLIDE 42

Draft

28

Construction via CBC vs random CBC

Pγ,2 for product weights, γ2

j = 0.1, s = 10 (left) and s = 5 (right).

rz = min{r : Erandom CBC[Pγ,2] ≤ (1 + z/100) Pγ,2(CBC)}. rz as a function of n for z = 5 ( ), z = 10 ( ) and z = 20 ( ).

26 210 214 22 26 210 n rz 26 210 214 22 26 210 n rz
slide-43
SLIDE 43

Draft

29

Quantiles of figure of merit

We computed Pγ,2 with product weights with γ2

j = 0.3 for all j, for all

admissible vectors a ∈ {1} × Us−1

n

, for n = 2e, ..., 219. For s = 2, a linear regression of log Pγ,2 vs log n for 212 ≤ n ≤ 219 gives decreasing rates of n−1.92 for the best, and n−1.87 and n−1.77 for the 10 % and 90 % quantiles. The mean decreases as n−1 an the worst-case as n0 (it is near 0.1948 for all n, obtained with a = (1, 1)).

27 210 213 216 219 10−10 10−5 100 n merit value s = 2 worst mean 90 % median 10 % best
slide-44
SLIDE 44

Draft

30

For s = 3, a linear regression log Pγ,2 vs log n for 210 ≤ n ≤ 215 gives decreasing rates of n−1.76 for the best, and n−1.64 and n−1.39 for the 10 % and 90 % quantiles. The mean decreases as n−1 and the worst-case as n0 (it is near 0.6393).

27 210 213 216 219 10−6 10−3 100 n merit value s = 3 worst mean 90 % median 10 % best
slide-45
SLIDE 45

Draft

31

Convergence of CBC and random CBC with r = 10 and r ≈ log n

10−10 10−6 10−2 merit value s = 2 10 % R-CBC (r = 10) R-CBC (r ≈ log n) CBC best 27 210 213 216 219 10−7 10−4 10−1 n merit value s = 3 10 % R-CBC (r = 10) R-CBC (r ≈ log n) CBC best
slide-46
SLIDE 46

Draft

32

Example: a stochastic activity network

Each arc j has random length Vj = F −1

j

(Uj). Let T = f (U1, . . . , U13) = length of longest path from node 1 to node 9. Want to estimate q(x) = P[T > x] for a given constant x.

1 source 2 V1 3 V2 V3 4 V4 5 V8 6 V5 V6 V10 7 V7 8 V9 V12 9 sink V11 V13
slide-47
SLIDE 47

Draft

33

To estimate q(x) by MC, we generate n independent realizations of T, say T1, . . . , Tn, and (1/n) n

i=1 I[Ti > x].

For RQMC, we replace the n realizations of (U1, . . . , U13) by the n points

  • f a randomly-shifted lattice.
slide-48
SLIDE 48

Draft

33

To estimate q(x) by MC, we generate n independent realizations of T, say T1, . . . , Tn, and (1/n) n

i=1 I[Ti > x].

For RQMC, we replace the n realizations of (U1, . . . , U13) by the n points

  • f a randomly-shifted lattice.

CMC estimator. Generate the Vj’s only for the 8 arcs that do not belong to the cut L = {5, 6, 7, 9, 10}, and replace I[T > x] by its conditional expectation given those Vj’s, P[T > x | {Vj, j ∈ L}]. This makes the integrand continuous in the Uj’s.

slide-49
SLIDE 49

Draft

33

To estimate q(x) by MC, we generate n independent realizations of T, say T1, . . . , Tn, and (1/n) n

i=1 I[Ti > x].

For RQMC, we replace the n realizations of (U1, . . . , U13) by the n points

  • f a randomly-shifted lattice.

CMC estimator. Generate the Vj’s only for the 8 arcs that do not belong to the cut L = {5, 6, 7, 9, 10}, and replace I[T > x] by its conditional expectation given those Vj’s, P[T > x | {Vj, j ∈ L}]. This makes the integrand continuous in the Uj’s. Illustration: Vj ∼ Normal(µj, σ2

j ) for j = 1, 2, 4, 11, 12, and

Vj ∼ Exponential(1/µj) otherwise. The µj: 13.0, 5.5, 7.0, 5.2, 16.5, 14.7, 10.3, 6.0, 4.0, 20.0, 3.2, 3.2, 16.5.

slide-50
SLIDE 50

Draft

34 1 2 V1 3 V2 V3 4 V4 5 V8 6 V5 V6 V10 7 V7 8 V9 V12 9 V11 V13
slide-51
SLIDE 51

Draft

35

ANOVA Variances for the Stochastic Activity Network

20 40 60 80 100 x = 64 x = 100 CMC, x = 30 CMC, x = 64 CMC, x = 100 % of total variance Order 1 Order 2 Order 3 Order 4 Order 5 Order 6 Order 7
slide-52
SLIDE 52

Draft

36

Lattices of Rank 1 with CBC

26 28 210 212 214 10−7 10−6 10−5 10−4 10−3 n variance Stochastic Activity Network (x = 64) MC Sobol M13,13,13,13,13,13 P2 order 2D P2 product v 2 u = (3/π2)|u|σ2 u P2 product v 2 u = (45/π4)|u|σ2 u P2 product Wang & Sloan (2006) n−2
slide-53
SLIDE 53

Draft

37

Lattices of Rank 1 with CBC

26 28 210 212 214 10−7 10−6 10−5 10−4 n variance Stochastic Activity Network (x = 100) MC Sobol M13,13,13,13,13,13 P2 order 2D P2 product v 2 u = (3/π2)|u|σ2 u P2 product v 2 u = (45/π4)|u|σ2 u P2 product Wang & Sloan (2006) n−2
slide-54
SLIDE 54

Draft

38

Lattices of Rank 1 with CBC

26 28 210 212 214 10−10 10−9 10−8 10−7 10−6 10−5 n variance Stochastic Activity Network (CMC x = 30) MC Sobol M13,13,13,13,13,13 weighted M13,13,13,13,13,13 P2 order 2D P2 order v 2 u = (3/π2)|u|σ2 u P2 product v 2 u = (3/π2)|u|σ2 u P2 product v 2 u = (3/π2)|u|σ2 u (no baker) P2 product v 2 u = (45/π4)|u|σ2 u P2 product Wang & Sloan (2006) n−2
slide-55
SLIDE 55

Draft

39

Lattices of Rank 1 with CBC

26 28 210 212 214 10−9 10−8 10−7 10−6 10−5 10−4 n variance Stochastic Activity Network (CMC x = 64) MC Sobol M13,13,13,13,13,13 weighted M13,13,13,13,13,13 P2 order 2D P2 order v 2 u = (3/π2)|u|σ2 u P2 product v 2 u = (3/π2)|u|σ2 u P2 product v 2 u = (3/π2)|u|σ2 u (no baker) P2 product v 2 u = (45/π4)|u|σ2 u P2 product Wang & Sloan (2006) n−2
slide-56
SLIDE 56

Draft

40

Lattices of Rank 1 with CBC

26 28 210 212 214 10−8 10−7 10−6 10−5 n variance Stochastic Activity Network (CMC x = 100) MC Sobol M13,13,13,13,13,13 weighted M13,13,13,13,13,13 P2 order 2D P2 order v 2 u = (3/π2)|u|σ2 u P2 product v 2 u = (3/π2)|u|σ2 u P2 product v 2 u = (3/π2)|u|σ2 u (no baker) P2 product v 2 u = (45/π4)|u|σ2 u P2 product Wang & Sloan (2006) n−2
slide-57
SLIDE 57

Draft

41

Random vs. Full CBC

26 28 210 212 214 10−10 10−9 10−8 10−7 10−6 n variance Stochastic Activity Network (CMC x = 30) Full CBC (P2 product) Random CBC (P2 product)
slide-58
SLIDE 58

Draft

42

Random vs. Full CBC

26 28 210 212 214 10−8 10−7 10−6 10−5 10−4 n variance Stochastic Activity Network (CMC x = 64) Full CBC (P2 product) Random CBC (P2 product)
slide-59
SLIDE 59

Draft

43

Random vs. Full CBC

26 28 210 212 214 10−8 10−7 10−6 10−5 n variance Stochastic Activity Network (CMC x = 100) Full CBC (P2 product) Random CBC (P2 product)
slide-60
SLIDE 60

Draft

44

Prime vs. Power-of-2 Number of Points

26 28 210 212 214 10−10 10−9 10−8 10−7 10−6 n variance Stochastic Activity Network (CMC x = 30) prime (P2 product) power of 2 (P2 product)
slide-61
SLIDE 61

Draft

45

Prime vs. Power-of-2 Number of Points

26 28 210 212 214 10−8 10−7 10−6 10−5 10−4 n variance Stochastic Activity Network (CMC x = 64) prime (P2 product) power of 2 (P2 product)
slide-62
SLIDE 62

Draft

46

Prime vs. Power-of-2 Number of Points

26 28 210 212 214 10−8 10−7 10−6 10−5 n variance Stochastic Activity Network (CMC x = 100) prime (P2 product) power of 2 (P2 product)
slide-63
SLIDE 63

Draft

47

Korobov vs. CBC

26 28 210 212 214 10−10 10−9 10−8 10−7 10−6 10−5 n variance Stochastic Activity Network (CMC x = 30) CBC (P2 product) Korobov (P2 product)
slide-64
SLIDE 64

Draft

48

Korobov vs. CBC

26 28 210 212 214 10−8 10−7 10−6 10−5 10−4 n variance Stochastic Activity Network (CMC x = 64) CBC (P2 product) Korobov (P2 product)
slide-65
SLIDE 65

Draft

49

Korobov vs. CBC

26 28 210 212 214 10−8 10−7 10−6 10−5 n variance Stochastic Activity Network (CMC x = 100) CBC (P2 product) Korobov (P2 product)
slide-66
SLIDE 66

Draft

50

Histograms

0.5 1 0.2 0.4 0.6 0.8 1 probability single MC draw (x = 100) 6 7 ·10−2 5 · 10−2 0.1 0.15 probability MC estimator (x = 100) 6.5 7 ·10−2 5 · 10−2 0.1 probability RQMC estimator (x = 100)
slide-67
SLIDE 67

Draft

51

Histograms

0.5 1 0.1 0.2 0.3 probability single MC draw (CMC x = 100) 6 6.5 7 ·10−2 5 · 10−2 0.1 0.15 probability MC estimator (CMC x = 100) 6.4 6.5 6.6 6.7 ·10−2 5 · 10−2 0.1 0.15 probability RQMC estimator (CMC x = 100)
slide-68
SLIDE 68

Draft

52

Conclusion and future work

– We (and you!) have a flexible software tool to search for good lattice rules based on various criteria, arbitrary weights, with several choices of method, arbitrary dimension and number of points. – Can be handy for applications, as well as for empirical investigations on lattice rules behavior. Future: – Consider other figures of merit, e.g., based on other function spaces. – Design and implement a similar software tool for digital nets.

slide-69
SLIDE 69

Draft

53