G.l.p., optimal coefficients, rank-1 lattice rules, ... Dirk - - PowerPoint PPT Presentation

Lattice rules Lattice sequences Cos space + tent transform Conclusions

G.l.p., optimal coefficients, rank-1 lattice rules, ...

Dirk Nuyens

Department of Computer Science KU Leuven, Belgium RICAM Special Semester on Algebra and Number Theory Workshop 1: Uniform distribution and quasi-Monte Carlo methods Austrian Academy of Sciences, Linz October 14–18, 2013

Lattice rules — Dirk Nuyens (KU Leuven) 1/41

Lattice rules Lattice sequences Cos space + tent transform Conclusions Definition

Rank-1 lattice rules

For a given generating vector z ∈ Zd

N take

QN(f ; z) := 1 N

N−1

• k=0

f zk mod N N

• ≈
• [0,1)d f (x) dx.

Korobov (1959), Sloan, Joe, Niederreiter, Kuo, Dick, ...

Components zj are taken relatively prime to N (→ permutations of k/N). Here: N = 17, z = (1, 5). Looks good.

♣ how to generate these points: N=17; z=[1; 5]; P = mod(z*(0:N-1), N)/N;

But, how to pick “good” z?

♣ goodandbad.m; try z = [1; 1];

Lattice rules — Dirk Nuyens (KU Leuven) 2/41

Lattice rules Lattice sequences Cos space + tent transform Conclusions Definition

Rank-1 lattice rules

For a given generating vector z ∈ Zd

N take

QN(f ; z) := 1 N

N−1

• k=0

f zk mod N N

• ≈
• [0,1)d f (x) dx.

Korobov (1959), Sloan, Joe, Niederreiter, Kuo, Dick, ...

Components zj are taken relatively prime to N (→ permutations of k/N). Here: N = 17, z = (1, 5). Looks good.

♣ how to generate these points: N=17; z=[1; 2]; P = mod(z*(0:N-1), N)/N;

But, how to pick “good” z?

♣ goodandbad.m; try z = [1; 1];

Lattice rules — Dirk Nuyens (KU Leuven) 2/41

Lattice rules Lattice sequences Cos space + tent transform Conclusions Definition

Rank-1 lattice rules

For a given generating vector z ∈ Zd

N take

QN(f ; z) := 1 N

N−1

• k=0

f zk mod N N

• ≈
• [0,1)d f (x) dx.

Korobov (1959), Sloan, Joe, Niederreiter, Kuo, Dick, ...

Components zj are taken relatively prime to N (→ permutations of k/N). Here: N = 17, z = (1, 5). Looks good.

♣ how to generate these points: N=17; z=[1; 5]; P = mod(z*(0:N-1), N)/N;

But, how to pick “good” z?

♣ goodandbad.m; try z = [1; 1];

Lattice rules — Dirk Nuyens (KU Leuven) 2/41

Lattice rules Lattice sequences Cos space + tent transform Conclusions Worst-case error

A Koskma–Hlawka error bound

Suppose f can be represented as f (x) =

• h∈Zd

ˆ f (h) exp(2πi h · x) where

• h∈Zd

|ˆ f (h)| < ∞, then QN(f ; z) − I(f ) = 1 N

N−1

• k=0
• h∈Zd

ˆ f (h) exp(2πi h · z k/N) − ˆ f (0)

Lattice rules — Dirk Nuyens (KU Leuven) 3/41

Lattice rules Lattice sequences Cos space + tent transform Conclusions Worst-case error

A Koskma–Hlawka error bound

Suppose f can be represented as f (x) =

• h∈Zd

ˆ f (h) exp(2πi h · x) where

• h∈Zd

|ˆ f (h)| < ∞, then QN(f ; z) − I(f ) =

• 0=h∈Zd

ˆ f (h)

• 1

N

N−1

• k=0

exp(2πi (h · z) k/N)

• Lattice rules — Dirk Nuyens (KU Leuven)

3/41

Lattice rules Lattice sequences Cos space + tent transform Conclusions Worst-case error

A Koskma–Hlawka error bound

Suppose f can be represented as f (x) =

• h∈Zd

ˆ f (h) exp(2πi h · x) where

• h∈Zd

|ˆ f (h)| < ∞, then QN(f ; z) − I(f ) =

• 0=h∈Zd

ˆ f (h) 1h·z≡0

Lattice rules — Dirk Nuyens (KU Leuven) 3/41

Lattice rules Lattice sequences Cos space + tent transform Conclusions Worst-case error

A Koskma–Hlawka error bound

Suppose f can be represented as f (x) =

• h∈Zd

ˆ f (h) exp(2πi h · x) where

• h∈Zd

|ˆ f (h)| < ∞, then QN(f ; z) − I(f ) =

• 0=h∈Zd

ˆ f (h) rα,γ(h) r−1

α,γ(h)

• =1

1h·z≡0

Lattice rules — Dirk Nuyens (KU Leuven) 3/41

Lattice rules Lattice sequences Cos space + tent transform Conclusions Worst-case error

A Koskma–Hlawka error bound

Suppose f can be represented as f (x) =

• h∈Zd

ˆ f (h) exp(2πi h · x) where

• h∈Zd

|ˆ f (h)| < ∞, then |QN(f ; z) − I(f )| ≤

• 0=h∈Zd

|ˆ f (h)| rα,γ(h) r−1

α,γ(h) 1h·z≡0

Lattice rules — Dirk Nuyens (KU Leuven) 3/41

Lattice rules Lattice sequences Cos space + tent transform Conclusions Worst-case error

A Koskma–Hlawka error bound

Suppose f can be represented as f (x) =

• h∈Zd

ˆ f (h) exp(2πi h · x) where

• h∈Zd

|ˆ f (h)| < ∞, then |QN(f ; z) − I(f )| ≤  

h∈Zd

|ˆ f (h)|prp

α,γ(h)

 

1/p

• =:f p,α,γ

   

• 0=h∈Zd

h·z≡0

r−q

α,γ(h)

   

1/q

• =:e(z,N;·p,α,γ)

for some positive function rα,γ(h).

Lattice rules — Dirk Nuyens (KU Leuven) 3/41

Lattice rules Lattice sequences Cos space + tent transform Conclusions Worst-case error

The worst-case error

The worst-case error in using Q for f ∈ F to approximate I(f ) e(Q; F) := sup

f ∈F f F≤1

• Q(f ) − I(f )
• ⇒ |Q(f ) − I(f )| ≤ f F e(Q; F).

For the Korobov space, with weights γu, f p

p,α,γ =

• h∈Zd

|ˆ fh|p  γ−1

u(h)

• j∈u(h)

|hj|α  

p

=

• h∈Zd

|ˆ fh|p rp

α,γ(h) < ∞.

Riesz representer of the error for lattice rule when p > 1: ξ(x) =

• 0=h∈Zd

r−q

α,γ(h) exp(2πi h · x).

Then e(z, N; · p,α,γ) = QN(ξ; z)1/q.

Lattice rules — Dirk Nuyens (KU Leuven) 4/41

Lattice rules Lattice sequences Cos space + tent transform Conclusions Worst-case error

Calculating the worst-case error

◮ Korobov space with product weights γu = j∈u γj:

χ(x) = −1 +

d

• j=1
• 1 + γq

j ω(xj)

• ,

where ω(x) =

• 0=h∈Z

exp(2πi hx) |h|αq . Need α > 1/q.

◮ The R-criterion → star-discrepancy:

χ(x) = · · · where ω(x) =

• 0=h∈
• − N

2 , N 2

• exp(2πi hx)

|h|αq . Need α > 0. This corresponds to truncated Fourier series.

Lattice rules — Dirk Nuyens (KU Leuven) 5/41

Lattice rules Lattice sequences Cos space + tent transform Conclusions Component-by-component

The component-by-component algorithm

Brute force search for good z ∈ Zd

N would cost exponential in d.

Therefore (Sloan, Joe, Kuo, ... and also Korobov): for s = 1 to d do for all zs ∈ Z×

N do

calculate es(z1, . . . , zs−1, zs) end for zs = argmin

z∈Z×

N

es(z1, . . . , zs−1, z) end for Cost of O(d2N 2).

◮ Rules with optimal rate of convergence in weighted

Korobov space, i.e., O(N −α+δ), δ > 0.

◮ Rules with optimal rate of convergence in weighted

Sobolev space, i.e., O(N −1+δ), δ > 0.

Kuo (2003), Dick (2004)

Lattice rules — Dirk Nuyens (KU Leuven) 6/41

Lattice rules Lattice sequences Cos space + tent transform Conclusions Fast component-by-component construction

Fast component-by-component construction

Component-by-component by using the recursion es((z∗

1, . . . , z∗ s−1, zs), N; || · ||p,α,γ)q

= es−1((z∗

1, . . . , z∗ s−1), N; || · ||p,α,γ)q + θs(zs).

Write ω(k · zj) := ω(xk,j) = ω({kzj/N}) = ω((kzj mod N)/N) =

• 0=h∈Z

r−q

α (h) exp(2πi hkzj/N),

Lattice rules — Dirk Nuyens (KU Leuven) 7/41

Lattice rules Lattice sequences Cos space + tent transform Conclusions Fast component-by-component construction

For prime N

... where θs(zs) =

• u⊆{1:s−1}

γu∪{s} 1 N

N−1

• k=0

Yu(k) ω(k · zs) = 1 N

N−1

• k=0

Ys(k) ω(k · zs) = 1 N Ys(0) ω(0) + 1 N

N−1

• k=1

Ys(k) ω(k · zs)

• convolution

with Yu(k) =

• hu∈Z|u|

∗

• j∈u

r−q

α (hj) exp(2πi hjz∗ j k/N) =

• j∈u

ω(k · z∗

j ),

Ys(k) =

• u⊆{1:s−1}

γu∪{s} Yu(k).

Lattice rules — Dirk Nuyens (KU Leuven) 8/41

Lattice rules Lattice sequences Cos space + tent transform Conclusions Fast component-by-component construction

Cyclic structure → circulant matrix

When N is prime ZN = Z×

N ∪ {0} ... and Ωn is basically just

homomorphic to the multiplication table of the group Z×

N.

E.g., for N = 7: · 1 2 3 4 5 6 1 1 2 3 4 5 6 2 2 4 6 1 3 5 3 3 6 2 5 1 4 4 4 1 5 2 6 3 5 5 3 1 6 4 2 6 6 5 4 3 2 1

Lattice rules — Dirk Nuyens (KU Leuven) 9/41

Lattice rules Lattice sequences Cos space + tent transform Conclusions Fast component-by-component construction

Cyclic structure → circulant matrix

When N is prime ZN = Z×

N ∪ {0} ... and Ωn is basically just

homomorphic to the multiplication table of the group Z×

N.

E.g., for N = 7: · 1 2 3 4 5 6 1 1 2 3 4 5 6 2 2 4 6 1 3 5 3 3 6 2 5 1 4 4 4 1 5 2 6 3 5 5 3 1 6 4 2 6 6 5 4 3 2 1 Multiplication modulo N can be much easier using a representation in powers of a generator g... gα gβ ≡ gα+β (mod ϕ(N)) (mod N).

Lattice rules — Dirk Nuyens (KU Leuven) 9/41

Lattice rules Lattice sequences Cos space + tent transform Conclusions Fast component-by-component construction

Cyclic structure → circulant matrix

When N is prime ZN = Z×

N ∪ {0} ... and Ωn is basically just

homomorphic to the multiplication table of the group Z×

N.

E.g., for N = 7: · 1 2 3 4 5 6 1 1 2 3 4 5 6 2 2 4 6 1 3 5 3 3 6 2 5 1 4 4 4 1 5 2 6 3 5 5 3 1 6 4 2 6 6 5 4 3 2 1 → · 1 5 4 6 2 3 1 1 5 4 6 2 3 3 3 1 5 4 6 2 2 2 3 1 5 4 6 6 6 2 3 1 5 4 4 4 6 2 3 1 5 5 5 4 6 2 3 1 Multiplication modulo N can be much easier using a representation in powers of a generator g... gα gβ ≡ gα+β (mod ϕ(N)) (mod N).

Lattice rules — Dirk Nuyens (KU Leuven) 9/41

Lattice rules Lattice sequences Cos space + tent transform Conclusions Fast component-by-component construction

Cyclic structure → circulant matrix

When N is prime ZN = Z×

N ∪ {0} ... and Ωn is basically just

homomorphic to the multiplication table of the group Z×

N.

E.g., for N = 7: · 1 2 3 4 5 6 1 1 2 3 4 5 6 2 2 4 6 1 3 5 3 3 6 2 5 1 4 4 4 1 5 2 6 3 5 5 3 1 6 4 2 6 6 5 4 3 2 1 → · 1 5 4 6 2 3 1 1 5 4 6 2 3 3 3 1 5 4 6 2 2 2 3 1 5 4 6 6 6 2 3 1 5 4 4 4 6 2 3 1 5 5 5 4 6 2 3 1 Multiplication modulo N can be much easier using a representation in powers of a generator g... gα gβ ≡ gα+β (mod ϕ(N)) (mod N).

Lattice rules — Dirk Nuyens (KU Leuven) 9/41

Lattice rules Lattice sequences Cos space + tent transform Conclusions Fast component-by-component construction

This is just a permutation on the rows and columns:

1 41

◮ The CBC algorithm multiplies with this matrix in each step. ◮ Matrix-vector multiplication with a circulant matrix Cm of

size m × m takes O(m log(m)) (instead of O(m2)): Cm x = IFFT(FFT(c) .* FFT(x)). ⇒ Fast construction of good lattice rule in time O(sN log(N)).

(N. & Cools, 2006a)

Lattice rules — Dirk Nuyens (KU Leuven) 10/41

Lattice rules Lattice sequences Cos space + tent transform Conclusions Fast component-by-component construction

The effect

1 min 10 mins 1 hour 1 day 1 month 1 year 8 years 120 years 102 103 104 105 106 107 108 n 10−3 100 103 106 109 total time for s = 20 (in secs) fastrank1 spmvrank1 rank1 slowrank1

Figure: Timings generated on a P4 2.4GHz ht, 2GB RAM for 20 dimensions

Note: this is known as Rader factorization for prime number FFTs, but is here done for general functions.

Lattice rules — Dirk Nuyens (KU Leuven) 11/41

Lattice rules Lattice sequences Cos space + tent transform Conclusions Non-prime N

Non-prime number of points: Chinese remainder

The technique for prime N makes use of a generator so...

◮ The multiplicative group Z× N = {z ∈ ZN : gcd(z, N) = 1}, is

cyclic whenever N = 2, 4, pk or 2pk, with p an odd prime. A generator for the cyclic group Z×

N is

called a primitive root modulo N.

◮ Every group Z× N, N = N1N2 · · · Nr, all relatively prime, can

be written as Z×

N ≃

• Z×

N1 ⊕ Z× N2 ⊕ · · · ⊕ Z× Nr

if 8 ∤ N, (−52k ⊕ 52k) ⊕ Z×

N2 ⊕ · · · ⊕ Z× Nr

if 8 | N, N1 = 2k, k ≥ 3, with every part a cyclic multiplicative group.

(N. & Cools, 2006b)

Lattice rules — Dirk Nuyens (KU Leuven) 12/41

Lattice rules Lattice sequences Cos space + tent transform Conclusions Non-prime N

1 2 3 6 9 18 5 10 15 30 45 90 1 2 3 6 9 18 5 10 15 30 45 90

Lattice rules — Dirk Nuyens (KU Leuven) 13/41

Lattice rules Lattice sequences Cos space + tent transform Conclusions Prime power sequences

Prime powers and sequences

Now take N = pM, a prime power, and set Ppm = z k mod pm pm : k ∈ Zpm

• ,

0 ≤ m ≤ M, then the point set is embedded P1 ⊂ Pp ⊂ Pp2 ⊂ · · · ⊂ PpM . We can construct such that (Cools, Kuo, N., 2006):

• 1. each embedded point set, Ppm, is a good lattice rule;
• 2. every subsequence anchored at and spanning a power of

the base, xk pm, xk pm+1, . . . , x(k+1) pm−1, is a good lattice.

Lattice rules — Dirk Nuyens (KU Leuven) 14/41

Lattice rules Lattice sequences Cos space + tent transform Conclusions Prime power sequences Lattice rules — Dirk Nuyens (KU Leuven) 15/41

1 3 9 27 81 1 3 9 27 81

Examples for a normal prime power case with N = 34 = 81 and the “odd” even prime power case with N = 27 = 128.

Lattice rules Lattice sequences Cos space + tent transform Conclusions Prime power sequences

Lattice sequences...

Base 3 sequence with 9, 64 (not a power of 3) and 81 points:

References on construction: N. (2005), N. & Cools (2006a), N. & Cools (2006b), Cools, Kuo & N. (2006), Hickernell, Kritzer, Kuo, & N. (2011), ... also constructions over finite fields for polynomial lattice rules and constructions for (weighted) degree of exactness...

Lattice rules — Dirk Nuyens (KU Leuven) 16/41

Lattice rules Lattice sequences Cos space + tent transform Conclusions Stop anywhere

Lattice sequences at intermediate N...

That is: A generating vector z ∈ Zs which is “good” for all values of N of the form bm. This gives a sequence of embedded lattice rules Pb0 ⊂ Pb1 ⊂ Pb2 ⊂ · · · ⊂ Pb10 ⊂ · · · ⊂ Pb20 ⊂ · · · We could stop anywhere using points in a permuted order: x(k) = ϕb,m(k) z mod bm bm = ϕb(k) z mod 1, m sufficiently large. Where ϕb,m gives a permutation on the integers 0, . . . , bm − 1 which keeps the embedding. E.g., ϕb could be van der Corput’s radical inverse function in base b.

Hickernell & Niederreiter (2003), Hickernell, Hong, L’Écuyer & Lemieux (2000), Cools, Kuo & N. (2006), Dick, Pillichshammer & Waterhouse (2008)

Question: what happens at intermediate N?

Lattice rules — Dirk Nuyens (KU Leuven) 17/41

Lattice rules Lattice sequences Cos space + tent transform Conclusions Stop anywhere

A classical result from Kuipers and Niederreiter (1974)

Theorem: Suppose you have M point sets Pi with Ni points, which have discrepancy D∗

Ni(Pi) each, then the combined

discrepancy satisfies D∗

N(P) ≤ M

• i=1

Ni N D∗

Ni(Pi),

with P =

i Pi and N = i Ni.

This is a remarkable result. There are no relational conditions

• n the different point sets!

Of course, you want to keep M small.

Lattice rules — Dirk Nuyens (KU Leuven) 18/41

Lattice rules Lattice sequences Cos space + tent transform Conclusions Stop anywhere

A worst-case error variation to the theme

The previous result is easily modified for the worst-case error: Theorem [HKKN2011]: Given M equal-weight rules Qi using point sets Pi, which have worst-case error eNi(Pi, K ), in a reproducing kernel Hilbert space H(K ), then the worst-case error for the equal-weight rule using all the points satisfies eN(P, K ) ≤

M

• i=1

Ni N eNi(Pi, K ), with P =

i Pi and N = i Ni.

Lattice rules — Dirk Nuyens (KU Leuven) 19/41

Lattice rules Lattice sequences Cos space + tent transform Conclusions Stop anywhere

What can be obtained by this?

Given a “good” lattice sequence in base b up to N points, where N =

L−1

• ℓ=0

nℓ bℓ, L = ⌊logb N⌋ + 1, and chop it up in blocks of powers of b: nb of blocks

• f size

with worst-case error . . . . . . . . . nℓ bℓ eℓ = O(b−ℓ(log bℓ)s−1) . . . . . . . . . Then, if each block has eℓ = O(b−ℓ(log bℓ)s−1), the total worst-case error behaves like eN(P, K ) = O (b − 1)(log N)s N

• .

Lattice rules — Dirk Nuyens (KU Leuven) 20/41

Lattice rules Lattice sequences Cos space + tent transform Conclusions Stop anywhere

O(N −1) anywhere for lattice sequences

A lattice sequence constructed in a shift-invariant reproducing kernel Hilbert space has exactly the right properties. E.g., for the first 22 points of a sequence in base 3:

S22

• S21
• S18
• S9
• x0, x1, x2, x3, . . . , x8
• P2,1

, x9, x10, x11, . . . , x17

• P2,2

, x18, x19, x20

• P1,1

, x21

• P0,1

the worst-case error of P2,2 is exactly that of the first 9 points, the worst-case error of P1,1 is exactly that of the first three points, etc...

Lattice rules — Dirk Nuyens (KU Leuven) 21/41

Lattice rules Lattice sequences Cos space + tent transform Conclusions Stop anywhere

But we want more...

The previous result showed that one can stop anywhere with a good lattice sequence. The error bound of O(N −1+ǫ) automatically holds for all N, not only powers of the base. However for some function spaces the worst-case error bounds are eN(SN, K ) = O(N −α+ǫ), N = bm, m = 0, 1, . . . , with α > 1, and ǫ > 0 hides the log N term. But using eN(P, K ) ≤

M

• i=1

Ni N eNi(Pi, K ), can at best only give O(N −1+ǫ). We need something else...

Lattice rules — Dirk Nuyens (KU Leuven) 22/41

Lattice rules Lattice sequences Cos space + tent transform Conclusions Stop anywhere

Equal-weight rules are stuck with O(N −1)

Theorem [HKKN2011]: Given a quasi-Monte Carlo rule Q(f ) = 1 N

N−1

• k=0

f (x(k)), then we have for the error for N + 1 points and N points: EN+1(f ) = I(f ) − 1 N + 1

N

• k=0

f (x(k)), EN(f ) = I(f ) − 1 N

N−1

• k=0

f (x(k)),

Lattice rules — Dirk Nuyens (KU Leuven) 23/41

Lattice rules Lattice sequences Cos space + tent transform Conclusions Stop anywhere

Equal-weight rules are stuck with O(N −1)

Theorem [HKKN2011]: Given a quasi-Monte Carlo rule Q(f ) = 1 N

N−1

• k=0

f (x(k)), then we have for the error for N + 1 points and N points: (N + 1) EN+1(f ) = (N + 1) I(f ) −

N

• k=0

f (x(k)), N EN(f ) = N I(f ) −

N−1

• k=0

f (x(k)),

Lattice rules — Dirk Nuyens (KU Leuven) 23/41

Lattice rules Lattice sequences Cos space + tent transform Conclusions Stop anywhere

Equal-weight rules are stuck with O(N −1)

Theorem [HKKN2011]: Given a quasi-Monte Carlo rule Q(f ) = 1 N

N−1

• k=0

f (x(k)), then we have for the error for N + 1 points and N points: (N + 1) EN+1(f ) = (N + 1) I(f ) −

N

• k=0

f (x(k)), N EN(f ) = N I(f ) −

N−1

• k=0

f (x(k)), such that using their difference gives |N EN(f )| + |(N + 1) EN+1(f )| ≥ |I(f ) − f (x(N))| = O(1). Suppose the error is O(N −β), then it follows that β ≤ 1.

Lattice rules — Dirk Nuyens (KU Leuven) 23/41

Lattice rules Lattice sequences Cos space + tent transform Conclusions Stop anywhere

Adding in weights

Theorem [HKKN2011]: If the weights are chosen as in Q(f ) :=

M

• i=1

wi Qi(f ), wi := N a

i

N a

1 + · · · + N a M

, a > 0, and the worst-case error for Qi in H(K ) satisfies eNi(Qi, K ) ≤ C N −α

i

, then the worst-case error for the compound rule Q satisifies eN(Q, K ) ≤ C M max(min(α,a),1) N min(α,a) . If a ≥ α and α > 1 then we have O(M α/N α).

Lattice rules — Dirk Nuyens (KU Leuven) 24/41

Lattice rules Lattice sequences Cos space + tent transform Conclusions Korobov’s sequence

Korobov’s sequence

If we are happy with O(N −1) then there is a different trick. Korobov: Construct a good lattice rule of d + 1-dimensions w.r.t. R (i.e., the star discrepancy) with generating vector z = (1, z1, . . . , zd), then leave out the first dimension. ...

Lattice rules — Dirk Nuyens (KU Leuven) 25/41

Lattice rules Lattice sequences Cos space + tent transform Conclusions Cosine space

Periodicity and lattice rules?

Up till now everything was periodic:

0.2 0.4 0.6 0.8 1.0 1.0 0.5 0.5 1.0

1 ∪ √ 2 cos(2πkx)

• k≥1

∪ √ 2 sin(2πkx)

• k≥1

Let’s change that:

0.2 0.4 0.6 0.8 1.0 1.0 0.5 0.5 1.0

1 ∪ √ 2 cos(πkx)

• k≥1

Lattice rules — Dirk Nuyens (KU Leuven) 26/41

Lattice rules Lattice sequences Cos space + tent transform Conclusions Cosine space

Cosine space

Can represent linear functions: x = 1 2 − 4 π2

∞

• k=1

k odd

cos(πkx) k2 . In fact, for the unanchored Sobolev space with r = 1, f , gK sob

1,γ =

1 f (x) dx 1 g(x) dx + 1 γ 1 f ′(x)g′(x) dx the reproducing kernel coincides (Dick, N., Pillichshammer 2013): K sob

1,γ (x, y) = 1 + γB1(x)B1(y) + γ B2(|x − y|)

2 = 1 + γ π2

∞

• k=1

1 k2 2 cos(kπx) cos(kπy) = K cos

1,γπ−2(x, y).

Lattice rules — Dirk Nuyens (KU Leuven) 27/41

Lattice rules Lattice sequences Cos space + tent transform Conclusions Tent transformation

Tent transformation

Tent transformed lattice points Pϕ(z, N) (a multiset):

◮ Tent transformation, ϕ : [0, 1] → [0, 1] is given by

ϕ(x) := 1 − |2x − 1|.

◮ Applied to each coordinate separately. ◮ Brings symmetry to the sampling points. ◮ For any k ∈ Zd +:

√ 2

|k|0 d j=1 cos(πkjϕ(xj)) =

√ 2

|k|0 d j=1 cos(2πkjxj).

Also “Bakers” transform in Hickernell (2000).

Lattice rules — Dirk Nuyens (KU Leuven) 28/41

Lattice rules Lattice sequences Cos space + tent transform Conclusions Tent transformation

Tent transformation

0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 Rank-1 lattice points with M = 34, z = [1, 21] 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 Tent transformed lattice points

Lattice rules — Dirk Nuyens (KU Leuven) 29/41

Lattice rules Lattice sequences Cos space + tent transform Conclusions Tent transformation

Why it works

0.2 0.4 0.6 0.8 1 −1 1 √ 2 cos(3πx) sampled with a 1d lattice, N = 15 0.2 0.4 0.6 0.8 1 −1 1 Tent transformed points 0.5 1 1.5 2 −1 1 Symmetry 0.2 0.4 0.6 0.8 1 −1 1 √ 2 cos(3πx) sampled with a 1d lattice, N = 16 0.2 0.4 0.6 0.8 1 −1 1 Tent transformed points 0.5 1 1.5 2 −1 1 Symmetry

Lattice rules — Dirk Nuyens (KU Leuven) 30/41

Lattice rules Lattice sequences Cos space + tent transform Conclusions Tent transformation

Results for cosine space

From (Dick, N., Pillichshammer 2013): e2(H(K cos

α,γ,s); Pϕ(z, N)) =

• h∈Zd

h·z≡0 (mod N)

rα,γ,s(h) = e2(H(K kor

α,γ,s); P(z, N)).

→ No random shifting needed for O(N −1) with non-periodic functions. Also: Construction: just construct for Korobov space!

Lattice rules — Dirk Nuyens (KU Leuven) 31/41

Lattice rules Lattice sequences Cos space + tent transform Conclusions Collocation

From integration to collocation

Consider the Poisson partial differential equation ∇2u(x) =

d

• j=1

∂2u(x) ∂x2

j

= f (x), ∀x ∈ Ω = (0, 1)d, with homogeneous Neumann boundary conditions u

• ∂Ω = 0.

Li and Hickernell (2003) studied this for periodic functions. Munthe-Kaas and Sorevik (2012) also for periodic functions.

Here: Assume u and f are expressed by cosine series, then u(x) = ˆ u0 +

• 0=k∈Zd

+

−ˆ f (k) π2k2

2

( √ 2)|k|0

d

• j=1

cos(πkjxj). Coefficients ˆ f (k) are calculated from samples on tent transformed rank-1 lattice points. → Limited to Ad(m).

Lattice rules — Dirk Nuyens (KU Leuven) 32/41

Lattice rules Lattice sequences Cos space + tent transform Conclusions Collocation

Weighted degrees

In Cools, Kuo, N. (2010) for integration in periodic case:

◮ Weighted trigonometric degree:

Ad,β(m) =

• h ∈ Zd :

d

• j=1

|hj| βj ≤ m

• .

◮ Weighted Zaremba degree (hyperbolic cross):

Ad,β(m) =

• h ∈ Zd :

d

• j=1

max

• 1, |hj|

βj

• ≤ m
• .

◮ Weighted product degree:

Ad,β(m) =

• h ∈ Zd : max

1≤j≤d

|hj| βj ≤ m

• .

Lattice rules — Dirk Nuyens (KU Leuven) 33/41

Lattice rules Lattice sequences Cos space + tent transform Conclusions Collocation

And their controlled sizes

◮ Weighted Zaremba degree (hyperbolic cross):

2m+1 ≤ |Ad,β(m)| ≤ mτ

d

• j=1
• 1 + 2ζ(τ)βτ

j

• for all

τ > 1.

◮ Weighted product degree:

|Ad,β(m)| =

d

• j=1
• 1 + 2⌊βjm⌋
• ≤ exp
• 2m

d

• j=1

βj

• .

◮ Weighted trigonometric degree (see product degree).

→ These sets are bounded in size independently of d if

∞

• j=1

βj < ∞.

Lattice rules — Dirk Nuyens (KU Leuven) 34/41

Lattice rules Lattice sequences Cos space + tent transform Conclusions Collocation

And their 2D Fourier spectra

Take m = 5, β1 = 1, β2 = 0.6 (and d = 2):

Trigonometric degree Zaremba degree Product degree

Lattice rules — Dirk Nuyens (KU Leuven) 35/41

Lattice rules Lattice sequences Cos space + tent transform Conclusions Collocation

What do we need for reconstruction?

We want the cosine coefficients ∀k ∈ Ad,γ(m) ˆ f (k) =

• [0,1]d f (x) (

√ 2)|k|0

d

• j=1

cos(πkjxj) dx. Using a tent-transformed lattice rule Pϕ(z, N): ˆ fa(k) = 1 N

• x∈Pϕ(z,N)

f (x) ( √ 2)|k|0

d

• j=1

cos(πkjxj) =

• ℓ∈Ad,γ(m)

ˆ f (ℓ) ( √ 2)|ℓ|0+|k|0 N

• x∈Pϕ(z,N)

d

• j=1

cos(πℓjxj) cos(πkjxj)

• =δk,ℓ

.

Lattice rules — Dirk Nuyens (KU Leuven) 36/41

Lattice rules Lattice sequences Cos space + tent transform Conclusions Collocation

What do we need for reconstruction

Discrete orthogonality:

Lemma

For exact evaluation of ˆ f (k), ∀k ∈ Ad,γ(m) we need 1 N

• x∈P(z,N)

( √ 2)|ℓ|0+|k|0 22d

• σ,σ′∈{±1}d

exp(2πi (σ(ℓ) − σ′(k)) · x) = δk,ℓ. Note: this uses the untransformed lattice points!

Lattice rules — Dirk Nuyens (KU Leuven) 37/41

Lattice rules Lattice sequences Cos space + tent transform Conclusions Collocation

Related work

Kämmerer, Kunis & Potts (2012) reconstructs trigonometric

polynomials with support on hyperbolic cross Ad,γ(m) ⊂ Zd.

◮ For the exact reconstruction it is necessary

1 N

• x∈P(z,N)

exp(2πi (ℓ − k) · x) = δk,ℓ.

◮ A difference set

Dd,γ

T

is then defined

• Dd,γ

T

:=

• ℓ − k ∈ Zd : ℓ, k ∈

Ad,γ(m)

• .

⇒ We can define a similar difference set Dd,γ

T

which collides with Dd,γ

T

.

Lattice rules — Dirk Nuyens (KU Leuven) 38/41

Lattice rules Lattice sequences Cos space + tent transform Conclusions Collocation

The difference sets

2 4 6 8 10 12 14 16 18 20 1 2 3 4 5 6 7 8 9 10

H

2,{ 1

2 , 1 4 }

40

−40 −30 −20 −10 10 20 30 40 −20 −15 −10 −5 5 10 15 20

D

2,{ 1

2 , 1 4 }

40

−20 −15 −10 −5 5 10 15 20 −10 −8 −6 −4 −2 2 4 6 8 10

• H

2,{ 1

2 , 1 4 }

40

−40 −30 −20 −10 10 20 30 40 −20 −15 −10 −5 5 10 15 20

• D

2,{ 1

2 , 1 4 }

40

Figure : Hyperbolic crosses and the difference sets for the cosine and Fourier series.

Lattice rules — Dirk Nuyens (KU Leuven) 39/41

Lattice rules Lattice sequences Cos space + tent transform Conclusions Results

Results so far

◮ We can construct these point sets:

using modified algorithm of Cools, Kuo, N. (2010), similar as the algorithm from Kämmerer, Kunis & Potts (2012).

◮ We have modified algorithms for FFT based reconstruction

and evaluation, similar as results from Kämmerer

◮ The weights are important in the complexity of the

algorithms! Manuscript in preparation with Cools and Suryanarayana.

Lattice rules — Dirk Nuyens (KU Leuven) 40/41

Lattice rules Lattice sequences Cos space + tent transform Conclusions

Conclusions

◮ Fast CBC constructions: O(dN log(N) + dT). ◮ Tent transformed lattice rule give you optimal rate O(N −1)

in the unanchored Sobolev space.

◮ Collocation with tent-transformed lattice rules (in

progress). Software available:

◮ Application of QMC point sets:

http://www.cs.kuleuven.be/~dirkn/qmc-generators/

◮ Construction of lattice rules:

http://www.cs.kuleuven.be/~dirkn/fast-cbc/

Lattice rules — Dirk Nuyens (KU Leuven) 41/41