Tractability Using Periodized Generalized Faure Sequences - - PowerPoint PPT Presentation

Tractability Using Periodized Generalized Faure Sequences

Christiane Lemieux

Department of Statistics and Actuarial Science University of Waterloo, Canada

Information-Based Complexity and Stochastic Computation Workshop ICERM, September 15, 2014

Christiane Lemieux University of Waterloo Tractability using PGFS 1 / 26

Motivation

Faure sequences and their generalizations are often used because they achieve the optimal value of 0 for the quality parameter t. Problem is they are not extensible in the dimension, since the base b must be at least as large as the dimension s to achieve t = 0 Can we make them extensible in the dimension, and if so, for what kind of problems would that work? Based on * C. Lemieux, H. Faure, New perspectives on (0, s)-sequences, in:

• P. L’Ecuyer, A. Owen (Eds.), Monte Carlo and Quasi-Monte Carlo

Methods 2008, Springer-Verlag, 2009, pp. 113–130. * C. Lemieux, Tractability using periodized generalized Faure sequences, to appear in Journal of Complexity, 2014.

Christiane Lemieux University of Waterloo Tractability using PGFS 2 / 26

Outline

1 Problem Setup and Review 2 Periodized Generalized Faure Sequences 3 Tractability Results 4 Numerical Results Christiane Lemieux University of Waterloo Tractability using PGFS 3 / 26

1 – Problem Setup and Review

Goal is to estimate Is(f ) =

• [0,1)s f (x)dx

where f : I s = [0, 1)s → R is a real-valued function. Monte Carlo or (randomized) quasi-Monte Carlo estimate Is(f ) by QN,s(f , PN) = 1 N

N

• i=1

f (xi) where PN = {x1, . . . , xN} ⊆ [0, 1)s. With (r)QMC, PN is a low-discrepancy point set: E(PN, z) = 1 N

N

• i=1

s

• j=1

1xi,j≤zj −

s

• j=1

zj, then low (star-)discrepancy means D∗(PN) = sup

z∈I s |E(PN, z)| ∈ O(N−1 logs N).

Christiane Lemieux University of Waterloo Tractability using PGFS 4 / 26

Example: want to estimate the exp. nb of clients waiting more than 5 minutes in a bank over fixed horizon T; fcts g and h resp. generate inter-arrival and service times Ai, Si, so f (x) = φ(Y) = nb of clients who waited more than 5 minutes Y = (A1 = g(x1), S1 = h(x2), . . . , AN = g(x2N−1), SN = h(x2N )) N = number of clients observed over horizon [0, T] Here s = 2N is not known ahead of time. In such cases, we need PN to be extensible in the dimension, i.e., coordinates for each xi can be added “on the fly”.

Christiane Lemieux University of Waterloo Tractability using PGFS 5 / 26

Review of Faure sequences

Elementary intervals in base b: subsets of I s of the form

s

• j=1

lj brj , lj + 1 brj

• (volume is b−M where M = r1 + . . . + rs)

where 0 ≤ lj < brj and rj ≥ 0, j = 1, . . . , s. M = 5 A (0, m, s)-net in base b is a point set PN with N = bm points s.t. any elementary interval of volume b−M contains bm−M pts from PN, when M ≤ m. A (0, s)-sequence in base b is a sequence of points x1, x2, . . . such that {xkbm+1, . . . , x(k+1)bm} is a (0, m, s)-net for all m ≥ 0 and all k ≥ 0.

Christiane Lemieux University of Waterloo Tractability using PGFS 6 / 26

Linearly scrambled van der Corput sequence in base b: For a prime base b, it is obtained by choosing a matrix C = (cr,k)r,k≥1 with elements in Zb and an infinite number of rows and columns, and then defining the nth term of this sequence as SC

b (n) := ∞

• r=0

yn,r br+1 in which yn,r =

∞

• k=0

cr+1,k+1 · ak(n), where digits ak(n) come from n =

k≥0 ak(n)bk.

Sequence of points in I s can be constructed using (SC1

b , . . . , SCs b ),

where C1, . . . , Cs are generating matrices. Original Faure sequences (1982): Cj = Pj−1

b

, where Pb is the (upper triangular) Pascal matrix Pb in Zb ⇒ (0, s)-sequence if b ≥ s. Generalized by Tezuka (1994) to Cj = AjPj−1

b

, where each Aj is an (NLT) matrix: also a (0, s)-sequence if b ≥ s.

Christiane Lemieux University of Waterloo Tractability using PGFS 7 / 26

Insight from Korobov Lattices

Korobov lattices are extensible in the dimension, but coordinates start repeating if s ≥ N: PN = i − 1 N (1, a, a2 mod N, . . . , as−1 mod N) mod 1, i = 1, . . . , N

• .

Can take N prime and a primitive element modulo N to get cycle of maximal period N − 1. With s ≥ N means projections over indices with lag of N − 1 are bad, but N large can mitigate this issue since corresponding projections of f are often not important in that case. Bank example where s is not fixed but on average equal to 2000; use Korobov lattice with N = 1021 and a = 76 MC Kor Sobol’ QN,s(f , PN) 543 544 544 HW 2.21 0.93 0.81

Christiane Lemieux University of Waterloo Tractability using PGFS 8 / 26

2 – Periodized Generalized Faure Sequences (PGFS)

Fix the base b, and for any s ≥ 1, use sequence defined over Zb based

• n generating matrices

Cj = AjPj−1

b

, j = 1, . . . , s (can be ≥ b). Scrambling matrices Aj fixed as follows: choose a period p ≈ b/2 and then let Aj = fj mod pI (diagonal) where the multipliers fj ∈ Zb are ordered according to the quality of the one-dimensional van der Corput sequence they generate (similar idea used to define generalized Faure sequences in LF09).

Christiane Lemieux University of Waterloo Tractability using PGFS 9 / 26

(49th,50th) coordinates of 1000 first points of Faure (left) and PGFS with b = 97

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

Christiane Lemieux University of Waterloo Tractability using PGFS 10 / 26

Of course, if s ≥ b, t > 0 for a PGFS, but we have tb = 0, where:

Definition

For sequence {x1, x2, . . .} over Zb, tk=0 if for each u = {j1, . . . , jr} satisfying 1 ≤ j1 < · · · < jr ≤ s, 1 ≤ |u| = r ≤ k, and r(u) := jr − j1 + 1 ≤ k, the corresponding projection {(xn,j1, xn,j2, . . . , xn,jr ), n ≥ 1} is a (0, r)-sequence over Zb. Result: PN the first N points of a digital sequence over Zb with tb = 0; Pu

N the projection of PN over u ⊆ {1, . . . , s} with r(u) ≤ b, then

D∗(Pu

N) ≤ 1

N b + 1 2b (b logb(bN))|u|.

Christiane Lemieux University of Waterloo Tractability using PGFS 11 / 26

3 – Tractability Results

Look at worst-case error over Hilbert space Hs with norm · Hs: e(PN; Hs) = sup{|Is(f ) − QN,s(f ; PN)| : f ∈ Hs, f Hs ≤ 1}, and compare it with the initial error, defined as e(0; Hs) = sup{|Is(f )| : f ∈ Hs, f Hs ≤ 1}. We then define n(ǫ, Hs) as the smallest n for which there exists Pn such that e(Pn; Hs) ≤ ǫ · e(0; Hs), where ǫ is in (0, 1). Integration over Hs is said to be QMC-tractable if there exist non-negative numbers C, p, and q such that n(ǫ, Hs) ≤ Cǫ−psq for all ǫ ∈ (0, 1) and all s ≥ 1. (1) If q = 0 in (1), then integration over Hs is QMC-strongly tractable, and the infimum of the numbers p satisfying (1) with q = 0 is called the ǫ−exponent of QMC-strong tractability.

Christiane Lemieux University of Waterloo Tractability using PGFS 12 / 26

Can adapt the results and proofs from SWW04 to our settings. As in SWW04, we assume Hs is a reproducing kernel Hilbert space with a reproducing kernel of the form Ks(x, y) =

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

γs,u

• j∈u

ηj(xj, yj), (2) and the weights γs,u are arbitrary non-negative numbers. Hence any f ∈ Hs satisfies f (·) = f (x), Ks(x, ·). Case A: anchored Sobolev space H(Ks,A): take ηj,A(x, y) = min(|xj − aj|, |y − aj|) if (x − aj)(y − aj) > 0 and 0

• therwise. The point (a1, . . . , as) is called the anchor.

Case B: unanchored Sobolev space H(Ks,B): take ηj,B(x, y) = 1

2B2(|x − y|) +

• x − 1

2

y − 1

2

• , where B2(·) is the

Bernoulli polynomial of degree 2, i.e., B2(x) = x2 − x + 1

6

Christiane Lemieux University of Waterloo Tractability using PGFS 13 / 26

Choice of weights

Tractability over weighted spaces usually occurs by choosing weights

• f the form γs,u =

j∈u γj for some 0 ≤ γj ≤ 1, j = 1, . . . , s, or when

there exists an integer r such that γs,u = 0 for all u with |u| > r (finite order). Here, we propose a special case of the latter, which makes use of the notion of range r(u) defined earlier.

Definition

A set of weights {γs,u}u⊆{1,...,s} is said to be of finite-range if there exists an integer R ∈ {0, . . . , s − 1} (called the range) such that γs,u = 0 if r(u) = maxj{j ∈ u} − minj{j ∈ u} + 1 > R.

Christiane Lemieux University of Waterloo Tractability using PGFS 14 / 26

Theorem: Let H(Ks) be the anchored Sobolev space H(Ks,A) with an arbitrary anchor a, or the unanchored Sobolev space H(Ks,B), and assume we have finite-range weights {γs,u}u⊆{1,...,s} of range R ≥ 1. Let PN be the first N points of a digital sequence over Zb such that tR = 0, where R ≤ b. Then e2(PN; H(Ks)) ≤ 1 N2

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

r(u)≤R

γs,u b + 1 2b 2 (2b logb bN)2|u|. (3)

Christiane Lemieux University of Waterloo Tractability using PGFS 15 / 26

Theorem: Let {γs,u}s≥1,u⊆{1,...,s} be weights of finite-range R ≥ 1. Let PN be the first N points of a digital sequence over Zb with b ≥ R and such that tR = 0 for all s ≥ 1. (a) Consider the anchored Sobolev space H(Ks,A) with arbitrary anchor a and weights γs,u. Then e(PN; H(Ks,A)) e(0; H(Ks)) ≤ C(b) 1 N (logb bN)b, where C(b) = (4 √ 3b)b. Furthermore, for any arbitrary δ > 0 there exists a constant Cδ independent of s and N such that e(PN; H(Ks,A)) e(0; H(Ks,A)) ≤ CδN−1+δ. Hence we have QMC-strong tractability with ǫ−exponent 1.

Christiane Lemieux University of Waterloo Tractability using PGFS 16 / 26

(b) Unanchored Sobolev space H(Ks,B) with weights γs,u; (i) If there exists c∗

B such that γs,u ≤ c∗ B for all u and s ≥ 1, then

e(PN; H(Ks,B)) e(0; H(Ks,B)) ≤ C1(b)s1/2 N (logb bN)b, where C1(b) = c∗

B(

√ 8b)b (so QMC-tractable with q = 1/2). (ii) If we further assume that M = sups=1,2,...

• u:0≤r(u)≤R γs,u < ∞ then

e(PN; H(Ks,B)) e(0; H(Ks,B)) ≤ C2(b) 1 N (logb bN)b, where C2(b) = √ M(2b)b. Furthermore, for any arbitrary δ > 0 there exists a constant CB,δ independent of s and N such that e(PN; H(Ks,B)) e(0; H(Ks,B)) ≤ CB,δN−1+δ. Hence we have QMC-strong tractability with ǫ−exponent 1.

Christiane Lemieux University of Waterloo Tractability using PGFS 17 / 26

4 – Numerical results

Want to demonstrate that PGFS can compete with good/popular choices Will compare with Sobol’, extensible rank-1 lattice of CKN06, extensible Korobov lattice (HHLL01 and GL07) Will use 3 different problems: test-function with finite-range weights, financial option, queueing example with unbounded s Use RQMC (shift for lattices, digital shift for Sobol’ and PGFS); can then look, e.g., at 1 m

m

• l=1

|QN,s(f , ˜ PN(l)) − Is(f )| (aver. absol. err.)

• r

1 m

m

• l=1

(QN,s(f , ˜ PN(l)) − ¯ QN,s)2 (estim. var.)

Christiane Lemieux University of Waterloo Tractability using PGFS 18 / 26

(i) – Test-function with finite-range weights

Based on the test function from Sobol’ and Asotsky g(x) =

s

• j=1

(1 + c(xj − 0.5)). Here however, we slightly modify this function and use instead gk,s(x) = 1 s − k + 1

s−k+1

• l=1

g(xl, . . . , xl+k−1)

• range k

for different values of s and k, so weights are of range k.

Christiane Lemieux University of Waterloo Tractability using PGFS 19 / 26

Example with g20,96 – PGFS with b = 97, p = 42

10

3

10

4

10

5

10

−4

10

−3

10

−2

10

−1

number of points aver absol error PGFS extkor rank1 sob Christiane Lemieux University of Waterloo Tractability using PGFS 20 / 26

Example with g20,250 – PGFS with b = 97, p = 42

10

3

10

4

10

5

10

−4

10

−3

10

−2

10

−1

number of points aver absol error PGFS extkor rank1 sob Christiane Lemieux University of Waterloo Tractability using PGFS 21 / 26

(ii) – Asian call option

Corresponding function f (x) is such that option price C(T, s, r, σ) satisfies C(T, s, r, σ) = E  e−rT max  1 s

s

• j=1

S(tj) − K, 0     =

• I s f (x)dx,

where S(tj) is the price of the underlying asset at time tj = jT/s. Assume S(tj)|S(tj−1) ∼ LN((r − σ2/2)/s, σ2/s), where T is the expiration time of the option, r is the risk-free rate, and σ is the volatility of the asset. Hence the function f (x) in this case can be written as f (x) = e−rTmax  0, 1 s

s

• j=1

exp((r − σ2/2)(jT/s) +

j

• l=1

Φ−1(xl)σ

• T

s   , where Φ(·) is the CDF of a standard Normal rv.

Christiane Lemieux University of Waterloo Tractability using PGFS 22 / 26

Asian option with s = 64 – PGFS with b = 241 and p = 122

10

3

10

4

10

5

10

−3

10

−2

10

−1

10 number of points aver absol error PGFS extkor rank1 sob naivePGFS Christiane Lemieux University of Waterloo Tractability using PGFS 23 / 26

Asian option with s = 256

10

3

10

4

10

5

10

−3

10

−2

10

−1

10 aver absol error PGFS extkor rank1 sob naivePGFS

Christiane Lemieux University of Waterloo Tractability using PGFS 24 / 26

(iii)- Bank example with unbounded s

Same as motivating example (i.e., average s is 2000): added MC for comparison, PGFS with b = 727 and p = 396

10

3

10

4

10

5

10

−2

10

−1

10 10

1

10

2

number of points var PGFS extkor rank1 sob MC

Christiane Lemieux University of Waterloo Tractability using PGFS 25 / 26

Some References

LF09: C. Lemieux, H. Faure, New perspectives on (0, s)-sequences, in: P. L’Ecuyer, A. Owen (Eds.), Monte Carlo and Quasi-Monte Carlo Methods 2008, Springer-Verlag, 2009, pp. 113–130. SWW04: I. H. Sloan, X. Wang, H. Wo´ zniakowski, Finite-order weights imply tractability of multivariate integration, Journal of Complexity 20 (2004) 46–74. CKN06: R. Cools, F. Y. Kuo, D. Nuyens, Constructing embedded lattice rules for multivariate integration, SIAM Journal on Scientific Computing 28 (6) (2006) 2162–2188. HHLL01: F.J. Hickernell, H.S. Hong, P. L’Ecuyer and C. Lemieux, Extensible Lattice Sequences for Quasi-Monte Carlo Quadrature, SIAM Journal on Scientific Computing, 22 (2001), 1117-1138. GL07: H. S. Gill, C. Lemieux, A search for extensible Korobov rules, Journal of Complexity 23 (2007) 603–613. Support from National Science and Engineering Research Council (NSERC) of Canada is acknowledged.

Christiane Lemieux University of Waterloo Tractability using PGFS 26 / 26