Scrambled geometric net integration over general product spaces
SScrambled geometric net integration over generalproduct spaces
Kinjal BasuStanford University Art B. OwenStanford UniversityMarch 2015
Abstract
Quasi-Monte Carlo (QMC) sampling has been developed for integra-tion over [0 , s where it has superior accuracy to Monte Carlo (MC) forintegrands of bounded variation. Scrambled net quadrature gives allowsreplication based error estimation for QMC with at least the same accu-racy and for smooth enough integrands even better accuracy than plainQMC. Integration over triangles, spheres, disks and Cartesian products ofsuch spaces is more difficult for QMC because the induced integrand on aunit cube may fail to have the desired regularity. In this paper, we presenta construction of point sets for numerical integration over Cartesian prod-ucts of s spaces of dimension d , with triangles ( d = 2) being of specialinterest. The point sets are transformations of randomized ( t, m, s )-netsusing recursive geometric partitions. The resulting integral estimates areunbiased and their variance is o (1 /n ) for any integrand in L of the prod-uct space. Under smoothness assumptions on the integrand, our random-ized QMC algorithm has variance O ( n − − /d (log n ) s − ), for integrationover s -fold Cartesian products of d -dimensional domains, compared to O ( n − ) for ordinary Monte Carlo. Keywords : Discrepancy, Multiresolution, Rendering, Scrambled net, Quasi-Monte Carlo.
Quasi-Monte Carlo (QMC) sampling is designed for problems of integration overthe unit cube [0 , s . Sampling over more complicated regions, such as the tri-angle or the sphere is more challenging. Measure preserving mappings from theunit cube to those spaces work very well for plain Monte Carlo. Unfortunately,the composition of the integrand with such a mapping may fail to have even themild smoothness properties that QMC exploits.In this paper, we consider quasi-Monte Carlo integration over product spacesof the form X s where X is a bounded set of dimension d . We are especiallyinterested in cases with d = 2 such as triangles, spherical triangles, spheres,1 a r X i v : . [ c s . NA ] M a r emispheres and disks. Integration over such sets is important in graphicalrendering (Arvo et al., 2001). For instance, when X is a triangle, an integralof the form (cid:82) X f ( x , x ) d x d x describes the potential for light to leave onetriangle and reach another. The function f incorporates the shapes and relativepositions of these triangles as well as whatever lies between them.Recent work by Basu and Owen (2014) develops two QMC methods for usein the triangle. One is a lattice like construction that was the first constructionto attain discrepancy O (log( n ) /n ) in that space. The other is a generalizationof the van der Corput sequence that makes a recursive partition of the triangle.In this paper, we generalize that van der Corput construction from the unittriangle to some other sets. We also replace the van der Corput sequence bydigital nets in dimension s , to obtain QMC points in X s . The attraction ofdigital nets is that they can be randomized in order to estimate our quadratureerror through independent replication of the estimate. Those randomizationshave the further advantage of reducing the error by about O ( n − / ) comparedto unrandomized QMC, when the integrand is smooth enough. For a survey ofrandomized QMC (RQMC) in general, see L’Ecuyer and Lemieux (2002). Foran outline of QMC for computer graphics, see Keller (2013).We study QMC and RQMC estimates of µ = 1 vol ( X ) s (cid:90) X s f ( x ) d x . Our estimates are equal weight rulesˆ µ = 1 n n (cid:88) i =1 f ( x i ) , where x i = φ ( u i ) (1)for random points u i ∈ [0 , s . The transformation φ maps [0 ,
1) into X and isapplied componentwise. We assume throughout that vol ( X ) = 1 whenever weare integrating over X , which simplifies several expressions.For any f ∈ L ( X s ) we find that Var(ˆ µ ) = o (1 /n ), so it is asymptoticallysuperior to plain Monte Carlo. We also find that for each finite n , scramblednets have a variance bounded by a finite multiple of the Monte Carlo variance,uniformly over all f ∈ L ( X s ).Our main result is that under smoothness conditions on f and a sphericityconstraint on the partitioning of X we are able to show that the estimate (1)attains Var(ˆ µ ) = O (cid:18) (log n ) s − n /d (cid:19) (2)when u i are certain scrambled digital nets. This variance rate is obtained viaa functional ANOVA decomposition of the integrand. The case d = 1 in (2)corresponds to the rate for scrambled net RQMC from Owen (1997b). Theprimary technical challenge in lifting that result from d = 1 to d > f ◦ φ is a well-behaved integrand.The statements above are for integration over X s but our proof of the mainresult is for integration over (cid:81) sj =1 X ( j ) where X ( j ) ⊂ R d are potentially different2ets of dimension d . Some of our results allow different dimensions d j for thespaces X ( j ) .The rest of the paper is organized as follows. In Section 2, we give back-ground material on digital nets and their scrambling. In Section 3, we presentrecursive geometric splits of a region X ⊂ R d and geometric van der Corput se-quences based on them. Section 4 generalizes those constructions to Cartesianproducts of s (cid:62) sd dimensions.We conclude this section by citing some related work on QMC over tensorproduct spaces. Tractability results have been obtained for integration over the s -fold product of the hypersphere S d = { x ∈ R d +1 | x T x = 1 } by Kuo andSloan (2005). Basu (2014) obtained such results for the s -fold product of thesimplex T d = { x ∈ [0 , d | (cid:80) j x j (cid:54) } . Those results are non-constructive. For s -fold tensor products of S there is a component-by-component constructionby Hesse et al. (2007). Both QMC and ordinary Monte Carlo (MC) correspond to the case with d = 1and X = [0 , , s takes x i ∼ U [0 , s . Thelaw of large numbers gives ˆ µ → µ with probability one when f ∈ L . If also f ∈ L then the root mean square error is σ/ √ n where σ is the variance of f ( x ) for x ∼ U [0 , s .QMC sampling improves upon MC by taking x i more uniformly distributedin [0 , s than random points usually are. Uniformity is measured via discrep-ancy. The local discrepancy of x , . . . , x n ∈ [0 , s at point a ∈ [0 , s is δ ( a ) = δ ( a ; x , . . . , x n ) = 1 n n (cid:88) i =1 x i ∈ [0 , a ) − vol ([0 , a )) . The star discrepancy of those points is D ∗ n ( x , . . . , x n ) = D ∗ n = sup a ∈ [0 , d | δ ( a ) | . The Koksma-Hlawka inequality is | ˆ µ − µ | (cid:54) D ∗ n V HK ( f ) , where V HK is the s -dimensional total variation of f in the sense of Hardy andKrause. For a detailed account of V HK see Owen (2005). Numerous construc-tions are known for which D ∗ n = O ((log n ) s − /n ) (Niederreiter, 1992) and soQMC is asymptotically much more accurate than MC when V HK ( f ) < ∞ .3 .1 Digital nets and sequences Of special interest here are QMC constructions known as digital nets (Niederre-iter, 1987; Dick and Pillichshammer, 2010). We describe them through a seriesof definitions. Throughout these definitions b (cid:62) s (cid:62) Z b = { , , . . . , b − } . Definition 1.
For k j ∈ N and c j ∈ Z b for j = 1 , . . . , s , the set s (cid:89) j =1 (cid:104) c j b k j , c j + 1 b k j (cid:17) is a b -adic box of dimension s . Definition 2.
For integers m (cid:62) t (cid:62)
0, the points x , . . . , x b m ∈ [0 , s area ( t, m, s )-net in base b if every b -adic box of dimension s with volume b t − m contains precisely b t of the x i .The nets have good equidistribution (low discrepancy) because boxes [0 , a ]can be efficiently approximated by unions of b -adic boxes. Digital nets canattain a discrepancy of O ((log( n )) s − /n ). Definition 3.
For integer t (cid:62)
0, the infinite sequence x , x , · · · ∈ [0 , s is a( t, s )-sequence in base b if the subsequence x rb m , . . . , x ( r +1) b m is a ( t, m, s )-netin base b for all integers r (cid:62) m (cid:62) t .The ( t, s )-sequences (called digital sequences) are extensible versions of ( t, m, s )-nets. They attain a discrepancy of O ((log( n )) s /n ). It improves to O ((log( n )) s − /n )along the subsequence n = λb m for integers m (cid:62) (cid:54) λ < b . Here we consider scrambling of digital nets and give several theorems for [0 , s that we generalize to X s . Let a ∈ [0 ,
1) have base b expansion a = (cid:80) ∞ k =1 a k b − k where a k ∈ Z b . If a has two base b expansions, we take the one with a tail of0s, not a tail of b − a k yielding x k ∈ Z b and deliver x = (cid:80) ∞ k =1 x k b − k . There are many different ways to choosethe permutations (Owen, 2003). Here we present the nested uniform scramblefrom Owen (1995).In a nested uniform scramble, x = π • ( a ) where π • is a uniform randompermutation (all b ! permutations equally probable). Then x = π • a ( a ), x = π • a a ( a ) and x k +1 = π • a a ...a k ( a k +1 ) where all of these permutations areindependent and uniform. Notice that the permutation applied to digit a k +1 depends on the previous digits. A nested uniform scramble of a = ( a , . . . , a s ) ∈ [0 , s applies independent nested uniform scrambles to all s components of a , sothat x j,k +1 = π j • a j a j ,...,a jk ( a j,k +1 ). A nested uniform scramble of a , . . . , a n ∈ [0 , s applies the same set of permutations to the digits of all n of those points.Propositions 1 and 2 are from Owen (1995).4 roposition 1. Let a ∈ [0 , s and let x be the result of a nested uniformrandom scramble of a . Then x ∼ U [0 , s . Proposition 2.
If the sequence a , . . . , a n is a ( t, m, s ) -net in base b , and x i are a nested uniform scramble of a i , then x i are a ( t, m, s ) -net in base b withprobability 1. Similarly if a i is a ( t, s ) -sequence in base b , then x i is a ( t, s ) -sequence in base b with probability 1. In scrambled net quadrature we estimate µ = (cid:82) [0 , s f ( x ) d x byˆ µ = 1 n n (cid:88) i =1 f ( x i ) , (3)where x i are a nested uniform scramble of a digital net a i .It follows from Proposition 1 that E (ˆ µ ) = µ for f ∈ L [0 , s . When f ∈ L [0 , s we can use independent random replications of the scramblednets to estimate the variance of ˆ µ . If V HK ( f ) < ∞ then we obtain Var(ˆ µ ) = O (log( n ) s − /n ) = O ( n − (cid:15) ) for any (cid:15) > Theorem 1.
Let f : [0 , s → R with continuous ∂ s ∂x ··· ∂x s f . Suppose that x i are a nested uniform scramble of the first n = λb m points of a ( t, s ) -sequence inbase b , for λ ∈ { , , . . . , b − } . Then for ˆ µ given by (3) , Var(ˆ µ ) = O (cid:16) log( n ) s − n (cid:17) = O ( n − (cid:15) ) as n → ∞ for any (cid:15) > .Proof. Owen (1997b) has this under a Lipschitz condition. Owen (2008) removesthat condition and corrects a Lemma from the first paper.Smoothness is not necessary for scrambled nets to attain a better rate thanMonte Carlo. Bounded variation is not even necessary:
Theorem 2.
Let x , . . . , x n be a nested uniform scramble of a ( t, m, s ) -net inbase b . Let f ∈ L ([0 , s ) . Then for ˆ µ given by (3) , Var(ˆ µ ) = o (cid:16) n (cid:17) as n → ∞ .Proof. This follows from Owen (1998). The case t = 0 is in Owen (1997a).The factor log( n ) s − is not necessarily small compared to n for reasonablesizes of n and large s . Informally speaking those powers cannot take effect forscrambled nets until after they are too small to make the result much worsethan plain Monte Carlo: 5 heorem 3. Let x , . . . , x n be a nested uniform scramble of a ( t, m, s ) -net inbase b . Let f ∈ L ([0 , s ) with Var( f ( x )) = σ when x ∼ U [0 , s . Then for ˆ µ given by (3) , Var(ˆ µ ) (cid:54) b t (cid:16) b + 1 b − (cid:17) s − σ n . If t = 0 , then Var(ˆ µ ) (cid:54) eσ /n . = 2 . σ /n .Proof. The first result is in Owen (1998), the second is in Owen (1997a).
The van der Corput sequence is constructed as follows. We begin with aninteger i (cid:62)
0. We write it as i = (cid:80) ∞ k =1 a k ( i ) b k − for digits a k ( i ) ∈ Z b . Definethe radical inverse function φ b ( i ) = (cid:80) ∞ k =1 a k ( i ) b − k ∈ [0 , b is x i = φ b ( i −
1) for i (cid:62)
1. It is a (0 , b .The original sequence of van der Corput (1935a,b) was for base b = 2. Any n consecutive van der Corput points have a discrepancy O (log( n ) /n ) where theimplied constant can depend on b .The lowest order base b digit of i determines which of b subintervals [ a/b, ( a +1) /b ) will contain x i . The second digit places x i into one of b sub-subintervals ofthe subinterval that the first digit placed it in, and so on. Basu and Owen (2014)used a base 4 recursive partitioning of the triangle to generate a triangularvan der Corput sequence. Discrepancy in the triangle is measured throughequidistribution over trapezoidal subsets (Brandolini et al., 2013). Triangularvan der Corput points have trapezoidal discrepancy of O ( n − / ). We begin with a notion of splitting sets. Splits are like partitions, except thatwe don’t require empty intersections among their parts.
Definition 4.
Let
X ⊂ R d have finite and positive volume. A b -fold split of X isa collection of Borel sets X a for a ∈ Z b with X = ∪ b − a =0 X a , vol ( X a ) = vol ( X ) /b for a ∈ Z b , and vol ( X a ∩ X a (cid:48) ) = 0 for 0 (cid:54) a < a (cid:48) < b .In all cases of interest to us, any overlap between X a and X a (cid:48) for a (cid:54) = a (cid:48) takes place on the boundaries of those sets. The unit interval [0 ,
1) is custom-arily partitioned into subintervals [ a/b, ( a + 1) /b ) in QMC. Handling X = [0 , X it could be burdensome to keeptrack of which subsets had which parts of their boundaries. Using splits allowsone for example to divide a closed triangle into four congruent closed triangles.Of course a partition is also a valid split. Our preferred approach uses a ran-domization under which there is probability zero of any sample point appearingon a split boundary. 6 BC ab c abc b = A BC a bc ab c abc b = A BC ab c a bc a bc a bc b = Figure 1: Splits of a triangle X for bases b = 2, 3 and 4. The subtriangles X j are labeled by the digit j ∈ Z b .Figure 1 shows a triangle X = ABC split into subtriangles. The left panelhas b = 2 subtriangles, the middle panel has b = 3 and the right panel has b = 4.For b = 2, the vertex labeled ‘A’ is connected to the midpoint of the oppositeside. The subset ABC is the one containing ‘B’. In each case, the new ‘A’ isthe mean of the old ‘B’ and ‘C’. The new ‘B’ is to the right as one looks fromthe new A towards the center of the triangle, and the new ‘C’ is on the left. Analgebraic description is more precise: Using lower case abc to describe the new ABC , the case described above (base 2 and digit 0) has abc = / /
21 0 00 1 0
ABC . Similar rules apply to the other bases. From such rules we may obtain thevertices of a split at level k by multiplying the original vertices by a sequenceof k × Definition 5.
Let
X ⊂ R d have finite and positive volume. A recursive b -foldsplit of X is a collection X of sets consisting of X and exactly one b -fold split ofevery set in the collection. The members of X are called cells.The original set X is said to be at level 0 of the recursive split. The cells X , . . . , X b − of X are at level 1. A member of a recursive split of X is atlevel k (cid:62) k splits of X . The cell X a ,a is the subset of X a corresponding to a = a and similarly, an arbitrary cell at level k (cid:62) X a ,a ,...,a k for a j ∈ Z b .We will need to enumerate all of the cells in a split X . For this we write t = (cid:80) kj =1 a j b j − ∈ Z b k and then take X ( k,t ) = X a ,a ,...,a k The cells in the splitare now X ( k,t ) for k ∈ N and t ∈ Z b k , with X (0 , = X .Figure 2 shows the first few levels of recursive splits for each of the splits fromFigure 1. The base 3 version has elements that become arbitrarily elongated as k Decomposition 3 Decomposition 4 Decomposition
Figure 2: The base b splits from Figure 1 carried out to k = 6 or 3 or 4 levels.increases. That is not desirable and our best results do not apply to such splits.The base 4 version was used by Basu and Owen (2014) to define a triangular vander Corput sequence. The base 2 version at 6 levels has a superficially similarappearance to the base 4 version at 3 levels. But only the latter has subtrianglessimilar to X . A linear transformation to make the parent ABC an equilateraltriangle yields congruent equilateral subtriangles for b = 4, while for b = 2 onegets isoceles triangles of two different shapes, some of which are elongated. The triangular sets were split in the same way at each level. Splits can be moregeneral than that, as we illustrate by splitting a disk. A convenient way to definea subset of a disk is in polar coordinates via upper and lower limits on the radiusand an interval of angles (which could wrap around 2 π ≡ b cells with aspectratios near one for b as large as several hundred. But their decompositions arenot recursive. Figure 3 shows eight levels in a recursive binary split of the diskinto cells. A cell with aspect ratio larger than one is split along the radial linethrough its centroid. Other cells are split into equal areas by an arc throughthe centroid.A spherical triangle is a subset of the sphere in R bounded by 3 great circles.By convention, one only considers spherical triangles for which all internal anglesare less than π radians. The spherical triangle can be split into b = 4 cells likethe rightmost panel in Figure 1. If one does so naively, via great circles throughmidpoints of the sides of the original triangle, the four cells need not have equalareas. Song et al. (2002) present a four-fold equal area recursive splitting for8 daptive binary splits of a disk Figure 3: A recursive binary equal area splitting of the unit disk, keeping theaspect ratio close to unity.spherical triangles, but their inner boundary arcs are in general small circles. Ifone wants a recursive splitting of great circles into great circles, then it can bedone by generalizing the b = 2 construction of the leftmost panel in Figure 1. Agreat circle with vertices ABC can be split into two by finding a point P on BCso that APB has half the area of ABC, using the first step of Arvo’s algorithmArvo (1995). Given a set X and a recursive splitting of it in base b we can construct ageometric van der Corput sequence for X . The integer i is written in base b as i = (cid:80) ∞ k =1 a k ( i ) b k − . To this i we define a sequence of sets X i : K = X a ( i ) ,a ( i ) ,...,a K ( i ) . Then x i is any point in ∩ ∞ K =1 X i : K . The volume of X i : K is b − K which convergesto 0 as K → ∞ . For most of the constructions we are interested in, each x i is auniquely determined point. For decompositions with bad aspect ratios, like thebase 3 decomposition in Figure 2, some of the set sequences converge to a linesegment. For instance if i ∈ { , , } , then a k ( i ) = 0 for k (cid:62) x i on one of the edges of the triangle.To get a unique limit x i , we use the notion of a sequence of sets convergingnicely to a point. Here is the version from Stromberg (1994).9 efinition 6. The sequence S k ∈ R d of Borel sets for k ∈ N converges nicelyto x ∈ R d as k → ∞ if there exists α < ∞ and d -dimensional cubes C k suchthat x ∈ C k , S k ⊆ C k , 0 < vol ( C k ) (cid:54) α vol ( S k ), and lim k →∞ diam( S k ) = 0.A sequence of sets that converges nicely to x cannot also converge nicely toany x (cid:48) (cid:54) = x . We generally assume the following condition. Definition 7.
A recursive split X in base b is convergent if for every infinitesequence a , a , a , · · · ∈ Z b , the cells X a ,a ,...a K converges nicely to a point as K → ∞ . That point is denoted lim K →∞ X a ,a ,...,a K .In a geometric van der Corput sequence, we take a convergent recursive splitand choose x i = lim M →∞ X a ( i − ,a ( i − ,...,a K ( i − , , , . . . , (cid:124) (cid:123)(cid:122) (cid:125) M where K is the last nonzero digit in the expansion of i − M (cid:62) x i is simply the center point of X a ( i − ,a ( i − ,...,a K ( i − . For b = 2, x i is an interior but noncentral point. Therecursive split for b = 3 is not convergent. Definition 8.
Let X be a recursive split of X ∈ R d in base b . Then X satisfiesthe sphericity condition if there exists C < ∞ such that diam( X a ,...,a k ) (cid:54) Cb − k/d holds for all cells X a ,...,a k in X .A recursive split that satisfies the sphericity condition is necessarily conver-gent. The constant C can be as low as 1 when d = 1 and the cells are intervals.The smallest possible C is usually greater than 1. We will assume without lossof generality that 1 (cid:54) C < ∞ . Definition 9.
Given a set
X ⊂ R d and a convergent recursive split X of X in base b , the X -transformation of [0 ,
1) is the function φ = φ X : [0 , → X given by φ ( x ) = lim K →∞ X x ,x ,...,x K where x has the base b representation0 .x x . . . . If x has two representations, then the one with trailing 0s is used. Let X be a bounded subset of R d with finite nonzero volume. Here we definedigital geometric nets in X s via splittings. It is convenient at this point togeneralize to an s -fold Cartesian product of potentially different spaces, eventhough they may all be copies of the same X .For s ∈ N , we represent the set { , , . . . , s } by 1: s . For j ∈ s we havebounded sets X ( j ) ⊂ R d j with vol ( X ( j ) ) = 1. For sets of indices u ⊆ s , thecomplement 1: s − u is denoted by − u . We use | u | for the cardinality of u . TheCartesian product of X ( j ) for j ∈ u is denoted X u . A vector x ∈ X s hascomponents x j ∈ X ( j ) . The vector in X u with components x j for j ∈ u isdenoted x u . 10 point in X s has (cid:80) sj =1 d j components. We write it as x = ( x , x , . . . , x s ).The components in this vector of vectors are those of the x j concatenated. Thenotation x j for d j consecutive components of x is the same as we use for the j ’th point in a quadrature rule. The usages are different enough that contextwill make it clear which is intended. Definition 10.
For j = 1 , . . . , s , let X j be a recursive split of X ( j ) in a commonbase b . Denote the cells of X j by X j, ( k,t ) for k ∈ N and t ∈ Z b k . Then a b -adiccell for these splits is a Cartesian product of the form (cid:81) sj =1 X j, ( k j ,t j ) for integers k j (cid:62) t j ∈ Z b kj . Definition 11.
Let X ( j ) ⊂ R d j have volume 1 for j ∈ s and let X j be arecursive split of X ( j ) in a common base b . For integers m (cid:62) t (cid:62)
0, the points x , . . . , x b m ∈ X s are a geometric ( t, m, s )-net in base b if every b -adic cellof volume b t − m contains precisely b t of the x i . They are a weak geometric( t, m, s )-net in base b if every b -adic cell of volume b t − m contains at least b t ofthe x i .Some of the b -adic cells can get more than b t points of a weak geometric( t, m, s )-net because the boundaries of those cells are permitted to overlap. Proposition 3.
Let a , . . . , a n be a ( t, m, s ) -net in base b . Let u , . . . , u n be anested uniform scramble of a , . . . , a n . For j ∈ s , let X j be a recursive base b split of the unit volume set X ( j ) ⊂ R d j with transformation φ j . Then z i = φ ( a i ) (componentwise) is a weak geometric ( t, m, s ) -net in base b and x i = φ ( u i ) (componentwise) is a geometric ( t, m, s ) -net in base b with probability one.Proof. In both cases the transformation applied to half open intervals placesenough points in each b -adic cell to make those points a weak geometric ( t, m, s )-net. The result for scrambled nets follows because each x i is uniformly dis-tributed and vol (cid:16) X j, ( k,t ) (cid:92) X j, ( k,t (cid:48) ) (cid:17) = 0for all j ∈ s , k ∈ Z b k and 0 (cid:54) t < t (cid:48) < b k . Let φ : [0 , → X ⊂ R d be the function that takes a point in the unit intervaland maps it to X according to the convergent recursive split X . We show herethat φ preserves the uniform distribution. This is the only section in which weneed to distinguish Lebesgue measures of different dimensions. To that end, weuse λ for Lebesgue measure in R and λ d for R d . Proposition 4.
Let
X ⊂ R d with vol ( X ) = 1 . Let X be a convergent recursivesplit of X in base b (cid:62) . Let φ be the X -transformation of [0 , and let A ⊆ X be a Borel set. Then λ ( φ − ( A )) = λ d ( A ) where φ − ( A ) = { x ∈ [0 , | φ ( x ) ∈ A } . roof. First, suppose that A = X a ,a ,...,a k for a j ∈ Z b . Then φ − ( A ) =[ t/b k , ( t + 1) /b k ) for some t ∈ Z b k , and so λ d ( A ) = 1 b k = λ (cid:16)(cid:104) tb k , t + 1 b k (cid:17)(cid:17) = λ ( φ − ( A )) . Now let A be any Borel subset of X . Given (cid:15) >
0, there exists a level k (cid:62) n cells B i at level k of X such that A ⊆ (cid:83) ni =1 B i and λ d ( (cid:83) ni =1 B i \ A ) < (cid:15) .Similarly, there exists a level k (cid:62) m cells C i at level k of X such that (cid:83) mi =1 C i ⊆ A and λ d ( A \ (cid:83) mi =1 C i ) < (cid:15) . Thus we get, λ d ( A ) = λ d (cid:18) n (cid:83) i =1 B i (cid:19) − λ d (cid:18) n (cid:83) i =1 B i \ A (cid:19) ≥ λ d (cid:18) n (cid:83) i =1 B i (cid:19) − (cid:15) = n (cid:88) i =1 λ d ( B i ) − (cid:15) = n (cid:88) i =1 λ ( φ − ( B i )) − (cid:15) = λ (cid:18) φ − (cid:18) n (cid:83) i =1 B i (cid:19)(cid:19) − (cid:15) ≥ λ ( φ − ( A )) − (cid:15), where the second equality follows because φ is bijective and therefore φ − ( A ) ⊆ φ − ( (cid:83) ni =1 B i ). Similarly, λ d ( A ) (cid:54) λ ( φ − ( A )) + (cid:15) . Since (cid:15) was arbitrary wehave the proof.Measure preservation extends to the multidimensional case. The next propo-sition combines that with uniformity under scrambling. Proposition 5.
Let X ( j ) ⊂ R d j with vol ( X ( j ) ) = 1 for j ∈ s have convergentrecursive splits X j in bases b j (cid:62) with corresponding transformations φ j . Let a ∈ [0 , s and let x j be a base b j nested uniform scramble of a j . Then φ ( x ) =( φ ( x ) , . . . , φ s ( x s )) ∼ U ( X s ) .Proof. By Proposition 1, x is uniformly distributed on [0 , s . By Proposition 4, φ preserves uniform measure. Thus φ j ( x j ) are independent U ( X ( j ) ) randomelements. L not requiring smoothness Some of the basic properties of scrambled nets go through for geometric scram-bled nets, without requiring any smoothness of the integrand. They don’t evenrequire that the same base be used to define both the transformations and thedigital net.
Theorem 4.
Let X ( j ) ⊂ R d j with vol ( X ( j ) ) = 1 for j = 1 , . . . , s . Let X j bea convergent recursive split of X ( j ) in base b j (cid:62) with transformation φ j . Let u , . . . , u n be a nested uniform scramble of a ( t, m, s ) -net in base b (cid:62) and let x i = φ ( u i ) componentwise. Then for any f ∈ L ( X s ) , Var(ˆ µ ) = o (cid:16) n (cid:17) as n → ∞ . roof. Since f ∈ L ( X s ) we have f ◦ φ ∈ L [0 , s . Then Theorem 2 applies. Theorem 5.
Under the conditions of Theorem 4,
Var(ˆ µ ) (cid:54) b t (cid:16) b + 1 b − (cid:17) s − σ n , where σ = Var( f ( x )) for x ∼ U ( X s ) . If t = 0 , then Var(ˆ µ ) (cid:54) eσ /n . =2 . σ /n .Proof. Once again, f ◦ φ ∈ L [0 , s . Therefore Theorem 3 applies. X s There is a well known analysis of variance (ANOVA) for [0 , s . Here we presentthe corresponding ANOVA for X s . Then we give a multiresolution of L ( X s )adapting the base b wavelet multiresolution in Owen (1997a) for [0 , X s For f ∈ L ( X s ) the ANOVA decomposition provides a term for each u ⊆ s .These are defined recursively via f u ( x ) = (cid:90) X − u (cid:16) f ( x ) − (cid:88) v (cid:40) u f v ( x ) (cid:17) d x − u . (4)The function f u represents the ‘effect’ of x j for j ∈ u above and beyond whatcan be explained by lower order effects of strict subsets v (cid:40) u . While f u is a function defined on X s its value only depends on x u . By convention f ∅ ( x ) = (cid:82) X s f ( x ) d x = µ for all x . We define variances σ u = (cid:82) X s f u ( x ) d x for | u | > σ ∅ = 0. The ANOVA decomposition satisfies (cid:80) | u | > σ u = σ where σ = (cid:82) X s ( f ( x ) − µ ) d x . It also satisfies f ( x ) = (cid:80) u ⊆ s f u ( x ) by thedefinition of f s , wherein a 0-fold integral of a function leaves it unchanged. We begin with a version of base b Haar wavelets adapted to
X ⊂ R d using arecursive split X of X in base b (cid:62)
2. Recall that the cells at level k of a splitare represented by one of X ( k,t ) for 0 (cid:54) t < b k . Those cells are in turn split atlevel k + 1 via X ( k,t ) = b − (cid:91) c =0 X ( k,t,c ) , where X ( k,t,c ) = X ( k +1 ,bt + c ) . X in terms of X has a function ϕ ( x ) = 1 for all x ∈ X as well as functions ψ ktc = b ( k +1) / x ∈X ( k,t,c ) − b ( k − / x ∈X ( k,t ) ≡ b ( k − / (cid:0) bN ktc ( x ) − W kt ( x ) (cid:1) , (5)where N ktc and W kt are indicator functions of the given narrow and wide cellsrespectively. The scaling in (5) makes the norm of ψ ktc independent of k : (cid:82) ψ ktc ( x ) d x = ( b − /b .For f , f ∈ L ( X ) define the inner product (cid:104) f , f (cid:105) = (cid:82) X f ( x ) f ( x ) d x .Then let f K ( x ) = (cid:104) f, ϕ (cid:105) ϕ ( x ) + K (cid:88) k =1 b k − (cid:88) t =0 b − (cid:88) c =0 (cid:104) f, ψ ktc (cid:105) ψ ktc ( x ) . For x belonging to only one cell at level K + 1, as almost all x do, f K ( x )is the average of f over that cell. By Lebesgue’s differentiation theorem, localaverages over sets that converge nicely to x satisfylim K →∞ (cid:82) S K f ( x ) d x vol ( S K ) = f ( x ) , a.e.for f ∈ L ( R d ). So if X is convergent, then lim K →∞ f K ( x ) = f ( x ) almosteverywhere. Thus we may use the representation f ( x ) = (cid:104) f, ϕ (cid:105) ϕ ( x ) + ∞ (cid:88) k =1 b k − (cid:88) t =0 b − (cid:88) c =0 (cid:104) f, ψ ktc (cid:105) ψ ktc ( x ) . (6)Equation (6) resembles a Fourier analysis with basis functions ϕ and ψ ktc . Un-like the Fourier case, the functions ψ ktc and ψ ktc (cid:48) are not orthogonal. Indeed (cid:80) c ∈ Z b ψ ktc = 0 a.e.. Non-orthogonal bases that nonetheless obey (6) are knownas tight frames.We may extend (6) to the multidimensional setting by taking tensor prod-ucts. For j ∈ s , let X ( j ) ⊂ R d have recursive split X j in base b (cid:62)
2. Let thebasis functions be ϕ j and ψ j ( ktc ) with narrow and wide cell indicators N jktc and W jkt . For u ⊆ s , let κ ∈ N | u | have elements k j (cid:62) j ∈ u . Similarly let τ have elements t j ∈ Z b kj and γ have elements c j ∈ Z b , both for j ∈ u . Then for x ∈ X s define ψ uκτγ ( x ) = (cid:89) j ∈ u ψ jk j t j c j ( x j ) (cid:89) j (cid:54)∈ u ϕ j ( x j ) . (7)Our multiresolution of L ( X s ) is f ( x ) = (cid:88) u ⊆ s (cid:88) κ | u (cid:88) τ | u,κ (cid:88) γ | u (cid:104) ψ uκτγ , f (cid:105) ψ uκτγ ( x )= µ + (cid:88) | u | > (cid:88) κ | u (cid:88) τ | u,κ (cid:88) γ | u (cid:104) ψ uκτγ , f (cid:105) ψ uκτγ ( x ) . The sum over κ is over all possible values of κ given the subset u . The othersums are similarly over their entire ranges given the other named variables.14 .3 Variance and gain coefficients Here we study the variance of averages over scrambled geometric nets. We startwith arbitrary points a i ∈ [0 , s . For now, they need not be from a digital net.They are given a nested uniform scramble, yielding points u i ∈ [0 , s . Thosepoints are then mapped to x i ∈ X s using recursive splits in base b .It follows from Proposition 5 thatVar(ˆ µ ) = E (cid:18) n n (cid:88) i =1 n (cid:88) i (cid:48) =1 (cid:88) | u | > (cid:88) κ | u (cid:88) τ | u,κ (cid:88) γ | u (cid:88) | u (cid:48) | > (cid:88) κ (cid:48) | u (cid:48) (cid:88) τ (cid:48) | u (cid:48) ,κ (cid:48) (cid:88) γ (cid:48) | u (cid:48) (cid:104) f, ψ uκτγ (cid:105)(cid:104) f, ψ u (cid:48) κ (cid:48) τ (cid:48) γ (cid:48) (cid:105) ψ uκτγ ( x i ) ψ u (cid:48) κ (cid:48) τ (cid:48) γ (cid:48) ( x i (cid:48) ) (cid:19) . This formula simplifies due to properties of the randomization. Lemma 4 fromOwen (1997a) shows that if u (cid:54) = u (cid:48) or κ (cid:54) = κ (cid:48) or τ (cid:54) = τ (cid:48) , then, E ( ψ uκτγ ( x i ) ψ u (cid:48) κ (cid:48) τ (cid:48) γ (cid:48) ( x i (cid:48) )) = 0 . (8)Consequently Var(ˆ µ ) = (cid:88) | u | > (cid:88) κ | u Var (cid:18) n n (cid:88) i =1 ν uκ ( x i ) (cid:19) , where ν uκ ( x ) = (cid:88) τ | u,κ (cid:88) γ | u (cid:104) f, ψ uκτγ (cid:105) ψ uκτγ ( x )with ν ∅ , () = µ . The function ν uκ is constant within elementary regions of theform (cid:89) j ∈ u X j, ( k j ,t j ,c j ) (cid:89) j (cid:54)∈ u X ( j ) for 0 ≤ t j < b k j and 0 (cid:54) c j < b .Define σ uκ = (cid:90) X s ν uκ ( x ) d x . The multiresolution-based ANOVA decomposition is σ = (cid:90) X s ( f ( x ) − µ ) d x = (cid:88) | u | > (cid:88) κ | u σ uκ (9)which follows from the orthogonality in (8).The equidistribution properties of a , . . . , a n determine the contribution ofeach ν uκ to Var(ˆ µ ). Write a i = ( a i , . . . , a is ) and defineΥ i,i (cid:48) ,j,k = 1 b − (cid:16) b (cid:98) b k +1 a ij (cid:99) = (cid:98) b k +1 a i (cid:48) j (cid:99) − (cid:98) b k a ij (cid:99) = (cid:98) b k a i (cid:48) j (cid:99) (cid:17) . | u | > κ ∈ N | u | defineΓ u,κ = 1 n n (cid:88) i =1 n (cid:88) i (cid:48) =1 (cid:89) j ∈ u Υ i,i (cid:48) ,j,k j . It follows from Theorem 2 of Owen (1997a) thatVar(ˆ µ ) = 1 n (cid:88) | u | > (cid:88) κ | u Γ u,κ σ uκ . We must have Γ u,κ (cid:62) µ ) (cid:62)
0. In plain Monte Carlo sampling,Var(ˆ µ ) = σ /n which corresponds to all Γ u,κ = 1 (compare (9)). The Γ u,κ arecalled ‘gain coefficients’ because they describe variance relative to plain MonteCarlo. If the points a i are carefully chosen, then many of those coefficients canbe reduced and an improvement over plain Monte Carlo can be obtained.If a , . . . , a n are a ( t, m, s )-net in base b , we can put bounds on the gaincoefficients using lemmas from Owen (1998). In particular, Γ u,κ = 0 if m − t (cid:62) | u | + | κ | , and otherwise Γ u,κ ≤ b t (cid:18) b + 1 b − (cid:19) s . Thus finally we have,Var(ˆ µ ) ≤ b t n (cid:18) b + 1 b − (cid:19) s (cid:88) | u | > (cid:88) | κ | + | u | >m − t σ uκ . (10)Equation (10) shows that the scrambled net variance depends on the rateat which σ uκ decay as | κ | + | u | increases. For smooth functions on [0 , s theydecay rapidly enough to give Var(ˆ µ ) = O (log( n ) s − /n ). To get a variance rateon X s we study the effects of smoothness on σ uκ for d -dimensional spaces X . Our main results for scrambling geometric nets require some smoothness ofthe integrand. We also use some extensions of the integrand and its ANOVAcomponents to rectangular domains.Let f be a real-valued function on X ⊆ R m . The dimension m will usuallybe d × s , for an s -fold tensor product of a d -dimensional region. For v ⊆ m ,the mixed partial derivative of f taken once with respect to x j for each j ∈ v is denoted ∂ v f . By convention ∂ ∅ f = f , as differentiating a function 0 timesleaves it unchanged. We present the Sobol’ extension through a series of definitions.16 efinition 12.
Let
X ⊆ R m for m ∈ N . The function f : X → R is said to besmooth if ∂ m f is continuous on X . Definition 13.
Let
X ⊂ R m . The rectangular hull of X is the Cartesianproduct rect( X ) = m (cid:89) j =1 (cid:2) inf { x j | x ∈ X } , sup { x j | x ∈ X } (cid:3) , which we also call a bounding box. For two points a , b ∈ R m we write rect[ a , b ]as a shorthand for rect[ { a , b } ].For later use, we note thatdiam(rect( X )) (cid:54) √ d × diam( X ) , (11)for X ⊂ R d . Definition 14.
A closed set
X ⊆ R m with non-empty interior is said to be Sobol’ extensible if there exists a point c ∈ X such that z ∈ X implies rect[ c , z ] ⊆X . The point c is called the anchor.Figure 4 shows some Sobol’ extensible regions. Figure 5 shows some setswhich are not Sobol’ extensible, because no anchor point exists for them. Setslike the right panel of Figure 5 are of interest in QMC for functions with inte-grable singularities along the diagonal. Proposition 6. If X ( j ) j ⊂ R d j is Sobol’ extensible with anchor c j for j =1 , . . . , s , then (cid:81) sj =1 X ( j ) is Sobol’ extensible with anchor c = ( c , . . . , c s ) .Proof. Suppose that x ∈ (cid:81) sj =1 X ( j ) . We write x as ( x , . . . , x s ) where each x j ∈ X ( j ) . Then rect( c , x ) ⊂ (cid:81) sj =1 rect( c j , x j ) ∈ (cid:81) sj =1 X ( j ) .Given points x , y ∈ R m and a set u ⊆ m , the hybrid point x u : y − u is thepoint z ∈ R m with z j = x j for j ∈ u and z j = y j for j (cid:54)∈ u . We will also requirehybrid points x u : y v : z w whose j ’th component is that of x or y or z for j in u or v or w respectively, where those index sets partition 1: m .A smooth function f can be written as f ( x ) = (cid:88) u ⊆ m (cid:90) [ c u , x u ] ∂ u f ( c − u : y u ) d y u (12)where (cid:82) [ c u , x u ] denotes ± (cid:82) rect[ c u , x u ] . The sign is negative if and only if c j > x j holds for an odd number of indices j ∈ u . The term for u = ∅ equals f ( c )under a natural convention. Equation (12) is a multivariable version of thefundamental theorem of calculus. For m = 1 it simplifies to f ( x ) = f ( c ) + (cid:82) xc f (cid:48) ( y ) d y . 17 c Figure 4: Sobol’ extensible regions. At left, X is the triangle with vertices(0 , , √ √ ,
0) and the anchor is c = (0 , X is a circular diskcentered its anchor c . The dashed lines depict some rectangular hulls joiningselected points to the anchor. Definition 15.
Let f be a smooth real-valued function on the Sobol’ extensibleregion X ⊂ R m . The Sobol’ extension of f is the function ˜ f : R m → R given by˜ f ( x ) = (cid:88) u ⊆ m (cid:90) [ c u , x u ] ∂ u f ( c − u : y u )1 c − u : y u ∈X d y u (13)where c is the anchor of X .The Sobol’ extension can be restricted to any domain X (cid:48) with X ⊂ X (cid:48) ⊂ R m .We usually use the Sobol’ extension for ˜ f from X to rect( X ). This extensionwas used in Sobol’ (1973) but not explained there. An account of it appearsin Owen (2005, 2006). For x ∈ X the factor 1 c − u : y u ∈X is always 1, making˜ f ( x ) = f ( x ) so that the term “extension” is appropriate.Two simple examples serve to illustrate the Sobol’ extension. If X = [0 , f ( x ) = x on 0 (cid:54) x (cid:54)
1, then the Sobol’ extension of f to [0 , ∞ ) is ˜ f ( x ) =min( x, X = [0 , and f ( x ) = x x on X , then the Sobol’ extension of f to [0 , ∞ ) is ˜ f ( x ) = min( x , × min( x , f has acontinuous mixed partial derivative ∂ m ˜ f for x in the interior of X and alsoin the interior of R m \ X where ∂ m ˜ f = 0 (Owen, 2005). As our examplesshow, ∂ u ˜ f for | u | > ∂ X . ASobol’ extensible X ⊆ R m has a boundary of m -dimensional measure 0, sowhen forming integrals of ∂ u ˜ f we may ignore those points or simply take thosepartial derivatives to be 0 there.The Sobol’ extension has a useful property that we need. It satisfies themultivariable fundamental theorem of calculus, even though some of its partialderivatives may fail to be continuous or even to exist everywhere. We can evenmove the anchor from c to an arbitrary point z .18 Figure 5: Non-Sobol’ extensible regions. At left, X is an annular regioncentered at the origin. At right, X is the unit square exclusive of an (cid:15) -widestrip centered on the diagonal. Theorem 6.
Let
X ⊂ R m be a Sobol’ extensible region and let f have a con-tinuous mixed partial ∂ m f on X . Let ˜ f be the Sobol’ extension of f and let z ∈ R m . Then ˜ f ( x ) = (cid:88) u ⊆ m (cid:90) [ z u , x u ] ∂ u ˜ f ( z − u : y u ) d y u . (14) Proof.
Define g ( x ) = (cid:88) v ⊆ m (cid:90) [ z v , x v ] ∂ v ˜ f ( z − v : t v ) d t v (15)where ∂ v ˜ f is taken to be 0 on those sets of measure zero where it might notexist when | v | >
0. We need to show that g ( x ) = ˜ f ( x ). Now˜ f ( z − v : t v ) = (cid:88) u ⊆ m (cid:90) [ c u , z u ∩− v : t u ∩ v ] ∂ u y f ( c − u : y u ) X ( c − u : y u ) d y u , (16)where for typographical convenience we have replaced 1 ·∈X by X ( · ). The sub-script in ∂ v y makes it easier to keep track of the variables with respect to whichthat derivative is taken. Substituting (16) into (15), we find that g ( x ) equals (cid:88) v ⊆ m (cid:90) [ z v , x v ] ∂ v t (cid:34) (cid:88) u ⊆ m (cid:90) [ c u , z u ∩− v : t u ∩ v ] ∂ u y f ( c − u : y u ) X ( c − u : y u ) d y u (cid:35) d t v = (cid:88) v ⊆ m (cid:90) [ z v , x v ] ∂ v t (cid:34) (cid:88) u ⊇ v (cid:90) [ c u , t u ] ∂ u y f ( c − u : y u ) X ( c − u : y u ) d y u (cid:35) d t v = (cid:88) v ⊆ m (cid:90) [ z v , x v ] (cid:88) u ⊇ v (cid:90) [ c u − v , t u − v ] ∂ u y f ( c − u : y u − v : t v ) X ( c − u : y u − v : t v ) d y u − v d t v . w = u − v and rewrite the sum, getting (cid:88) w ⊆ m (cid:88) v ⊆− w (cid:90) [ z v , x v ] (cid:90) [ c w , t w ] ∂ w + v y f ( c − w − v : y w : t v ) X ( c − w − v : y w : t v ) d y w d t v = (cid:88) w ⊆ m (cid:88) v ⊆− w (cid:90) [ z v , x v ] ∂ v y (cid:90) [ c w , t w ] ∂ w y f ( c − w − v : y w : t v ) X ( c − w − v : y w : t v ) d y w d t v . Any term above with v (cid:54) = ∅ vanishes. Therefore g ( x ) = (cid:88) w ⊆ m (cid:90) [ z ∅ , x ∅ ] (cid:90) [ c w , t w ] ∂ w y f ( c − w : y w ) X ( c − w : y w ) d y w d t ∅ = (cid:88) w ⊆ m (cid:90) [ c w , t w ] ∂ w y f ( c − w : y w ) X ( c − w : y w ) d y w = ˜ f ( x ) . Here we assume that X is a bounded closed set with non-empty interior, notnecessarily Sobol’ extensible. Sobol’ extensible spaces may fail to have a non-empty interior, but outside such odd cases, they are a subset of this class.Non-Sobol’ extensible regions like those in Figure 5 are included. To handledomains X of greater generality, we require greater smoothness of f .Let k ∈ N m be any multi-index with | k | = k + . . . + k m ≤ m . We denotethe k -th order partial derivative as D k f ( x ) = ∂ | k | ∂x k · · · ∂x k m m f ( x , . . . , x m ) . Definition 16.
A real-valued function f on X ⊂ R m is in C m ( X ) if all partialderivatives of f up to total order m are continuous on X .Whitney’s extension of a function in C m ( X ) to a function in C m (rect( X ))is given by the following lemma. Lemma 1.
Let f ∈ C m ( X ) for a bounded closed set X ⊂ R m with non-emptyinterior. Then there exists a function ˜ f ∈ C m (rect( X )) with the followingproperties:1. ˜ f ( x ) = f ( x ) for all x ∈ X ,2. D k ˜ f ( x ) = D k f ( x ) for all | k | ≤ m and x ∈ X , and3. ˜ f is analytic on rect( X ) \ X .Proof. The extension we need is the one provided by Whitney (1934). A functionin C m ( X ) in the ordinary sense is a fortiori in C m ( X ) according to Whitney’sdefinition. We use the restriction of Whitney’s function to the domain rect( X ).20e will need one more condition on X . We require the boundary of X tohave m -dimensional measure zero. Then Theorem 6 in which the fundamentaltheorem of calculus applies to ˜ f , holds also for the Whitney extension. Here we show that the ANOVA components of our smooth extensions are alsosmooth. We suppose that each X j ⊂ R d j and we let m = (cid:80) sj =1 d j . TheCartesian product X s is now a subset of R m . Lemma 2.
Let f be a smooth function on Sobol’ extensible X s ⊂ R m and for u ⊆ s let f u be the ANOVA component from (4) . Then f u is smooth on X s .Proof. We prove this by induction on | u | . Let | u | = 0, that is u = ∅ . Then f u ( x ) = (cid:82) X s f ( x ) d x which is a constant µ and is therefore smooth on X s .Let us suppose that the hypothesis holds for | u | = k − < s and we shall proveit for | u | = k .Fix any u ⊆ s such that | u | = k . By (4) we have, f u ( x ) = (cid:90) X − u f ( x ) d x − u − (cid:88) w ⊂ u f w ( x ) , using the fact that f w ( x ) does not depend on x j for j (cid:54)∈ w . Each term inthe summation is f w for | w | ≤ k − v ⊆ m . Now since f is smooth, ∂ v f ( x ) is continuous on X s and henceapplying Leibniz’s integral rule we have, ∂ v (cid:90) X − u f ( x ) d x − u = (cid:90) X − u ∂ v f ( x ) d x − u . Now the right hand side is the integral of a continuous function and is therefore acontinuous function. Thus the induction hypothesis hold for | u | = k completingthe proof. Lemma 3.
Let f ∈ C m ( X s ) for a bounded closed set X s ∈ R m and for u ⊆ s let f u be the ANOVA component in (4) . Then f u ∈ C m ( X s ) .Proof. The proof goes along the same lines as Lemma 2. We replace v in thatproof by any multi-index (cid:96) with | (cid:96) | ≤ m . Now since f is smooth D (cid:96) f ( x ) iscontinuous on X s and hence applying the Leibniz’s integral rule we have, D (cid:96) (cid:90) X − u f ( x ) d x − u = (cid:90) X − u D (cid:96) f ( x ) d x − u . Now the right hand side is the integral of a continuous function over certainvariables and is therefore a continuous function.21ow for a smooth function f defined on a product X s of Sobol’ extensiblesets, or on a product of more general spaces but with the smoothness requiredfor a Whitney extension, there exists an extension ˜ f on rect( X s ) such that˜ f ( x ) = (cid:88) u ⊆ m (cid:90) [ c u , x u ] ∂ u ˜ f ( c − u : y u ) d y u (17)for some point c ∈ rect( X s ). Here we prove that the variance of averages over scrambled geometric nets is O ( n − − /d log( n ) s − ), under smoothness and sphericity conditions. The proofis similar to the one in Owen (2008) for scrambled nets. We begin with notationfor some Cartesian products of cells. For this section we assume that d j = d isa constant dimension for all j ∈ s .Let b be the common base for recursive splits X j of X ( j ) ⊂ R d for j ∈ s .Let κ = ( k , . . . , k s ) and τ = ( t , . . . , t s ) be s -vectors with k j ∈ N and t j ∈ Z b kj .Then we write B uκτ = (cid:89) j ∈ u X j, ( k j ,t j ) (cid:89) j (cid:54)∈ u X ( j ) and (cid:101) B uκτ = rect( B uκτ ). For j = 1 , . . . , s , let S j = (( j − d + 1):( jd ) and thenfor u ⊆ s , define S u = (cid:91) j ∈ u S j . (18)Now let S u = { T ⊆ S u | T ∩ S j (cid:54) = ∅ , ∀ j ∈ u } . These are the subsets of S u thatcontain at least one element of S j for each j ∈ u . There are 2 d − S j , and so | S u | = (2 d − | u | . (19) Lemma 4.
Suppose that f is a smooth function on the Sobol’ extensible region X s ⊆ R ds , with extension ˜ f . Let each X ( j ) have a convergent recursive splitin base b whose k -level cells have diameter at most Cb − k/d for (cid:54) C < ∞ . Let u ⊆ s and let κ and τ be | u | -tuples with components k j ∈ N and t j ∈ Z b kj ,respectively for j ∈ u . Let ψ uκτγ be the multiresolution basis function (7) definedby the splits of X s . Then |(cid:104) f, ψ uκτγ (cid:105)| ≤ (cid:16) − b (cid:17) | u | b − | κ | (1+ d ) − | u | (cid:88) v ∈ S u d | v | / C | v | sup y ∈ (cid:101) B uκτ | ∂ v ˜ f u ( y ) | . (20) If f ∈ C ds ( X s ) with Whitney extension ˜ f , where each X ( j ) is a bounded closedset with non-empty interior and a boundary of measure zero, then (20) holdsregardless of whether X s is Sobol’ extensible. ∂ v ˜ f u takes the value 0 in places where itis not well defined. Alternatively one could use the essential supremum insteadof the supremum in (20). Later when we use (cid:107) · (cid:107) ∞ it will denote the essentialsupremum of its argument. Proof.
The same proof applies to both smoothness assumptions. From the def-inition we have (cid:104) f, ψ uκτγ (cid:105) = (cid:104) f u , ψ uκτγ (cid:105) = (cid:90) X − u (cid:90) X u f u ( x ) ψ uκτγ ( x ) d x u d x − u = b − ( | κ | + | u | ) / (cid:90) X u f u ( x ) (cid:89) j ∈ u b k j (cid:0) bN jk j t j c j ( x j ) − W jk j t j ( x j ) (cid:1) d x u . By either Lemma 2 or Lemma 3, f u is smooth and we let ˜ f u be its extension.We know ˜ f u ( x ) = f u ( x ) for all x ∈ X s . As the above integral is over X u , wecan write it as b − ( | κ | + | u | ) / (cid:90) X u ˜ f u ( x ) (cid:89) j ∈ u b k j (cid:0) bN k j t j c j ( x j ) − W k j t j ( x j ) (cid:1) d x u . (21)Now ˜ f u is smooth on rect( X u ) and depends only on x u . Applying (17) we canwrite, ˜ f u ( x ) = (cid:88) v ⊆ S u (cid:90) [ z v , x v ] ∂ v ˜ f u ( z − v : y v ) d y v , (22)choosing to place the anchor z at the center of (cid:101) B uκτ . Note that if v (cid:54)∈ S u , thenthere exists an index j ∈ u such that S j ∩ v = ∅ and then the integral in (22)above does not depend on x j making it orthogonal to bN jk j t j c j ( x j ) − W jk j t j ( x j ).Also the integrand in (21) is supported only for x u ∈ B uκτ . Putting thesetogether we get, b ( | κ | + | u | ) / |(cid:104) f, ψ uκτγ (cid:105)| = (cid:12)(cid:12)(cid:12)(cid:12)(cid:90) (cid:88) v ∈ S u (cid:90) [ z v , x v ] ∂ v ˜ f u ( z − v : y v ) d y v (cid:89) j ∈ u b k j (cid:0) bN jk j t j c j ( x j ) − W jk j t j ( x j ) (cid:1) d x u (cid:12)(cid:12)(cid:12)(cid:12) ≤ (cid:88) v ∈ S u sup x u ∈ B uκτ (cid:12)(cid:12)(cid:12)(cid:12)(cid:90) [ z v , x v ] ∂ v ˜ f u ( z − v : y v ) d y v (cid:12)(cid:12)(cid:12)(cid:12) × (cid:90) X u (cid:89) j ∈ u b k j (cid:12)(cid:12) bN jk j t j c j ( x j ) − W jk j t j ( x j ) (cid:12)(cid:12) d x u = (cid:16) − b (cid:17) | u | (cid:88) v ∈ S u sup x u ∈ B uκτ (cid:12)(cid:12)(cid:12)(cid:12)(cid:90) [ z v , x v ] ∂ v ˜ f u ( z − v : y v ) d y v (cid:12)(cid:12)(cid:12)(cid:12) . ∂ v ˜ f u is bounded we can write, (cid:12)(cid:12)(cid:12)(cid:12)(cid:90) [ z v , x v ] ∂ v ˜ f u ( z − v : y v ) d y v (cid:12)(cid:12)(cid:12)(cid:12) ≤ vol (rect[ z v , x v ]) sup y ∈ ˜ B uκτ | ∂ v ˜ f u ( y ) | . Because z ∈ ˜ B uκτ we have vol (rect[ z v , x v ]) = (cid:89) (cid:96) ∈ v | z (cid:96) − x (cid:96) | (cid:54) C | v | d | v | / (cid:89) j ∈ v ( b − k j /d ) | v ∩ S j | (cid:54) ( Cd / ) | v | b −| κ | /d . The last inequality follows because | v ∩ S j | (cid:62) j ∈ u and also usesequation (11) on the diameter of a bounding box. Finally, putting it all together,we get |(cid:104) f, ψ uκτγ (cid:105)| ≤ (cid:16) − b (cid:17) | u | b − | κ | ( d ) − | u | (cid:88) v ∈ S u C | v | d | v | / sup y ∈ (cid:101) B uκτ | ∂ v ˜ f u ( y ) | . The factor d | v | / in the bound can be as large as d s/ in applications, whichmay be quite large. It arises as a | v | -fold product of ratios diam(rect( · )) / diam( · )for cells. For rectangular cells that product is 1. Similarly for cells that are ‘axisparallel’ right-angle triangles, the product is again 1. Lemma 5.
Let h u ( z ) = max v ∈ S u | ∂ v ˜ f u ( z ) | for z ∈ (cid:101) B uκτ . Under the conditionsof Lemma 4, σ uκ ≤ (cid:20) (cid:101) C (cid:16) − b (cid:17) (cid:21) | u | b − | κ | /d (cid:107) h u (cid:107) ∞ where (cid:101) C = d / (2 d − C d .Proof. The supports of ψ uκτγ and ψ uκτ (cid:48) γ (cid:48) are disjoint unless τ = τ (cid:48) . Therefore ν uκ ( x ) = (cid:88) τ | u (cid:88) γ,γ (cid:48) | u (cid:104) f, ψ uκτγ (cid:105)(cid:104) f, ψ uκτγ (cid:48) (cid:105) ψ uκτγ ( x ) ψ uκτγ (cid:48) ( x ) . Now σ uκ = (cid:90) X s ν uκ ( x ) d x = (cid:88) τ | u,κ (cid:88) γ,γ (cid:48) | u (cid:104) f, ψ uκτγ (cid:105)(cid:104) f, ψ uκτγ (cid:48) (cid:105) (cid:90) X s ψ uκτγ ( x ) ψ uκτγ (cid:48) ( x ) d x = (cid:88) τ | u,κ (cid:88) γ,γ (cid:48) | u (cid:104) f, ψ uκτγ (cid:105)(cid:104) f, ψ uκτγ (cid:48) (cid:105) (cid:89) j ∈ u (1 c j = c (cid:48) j − b − ) ≤ (cid:16) − b (cid:17) | u | b −| κ | (1+ d ) −| u | (cid:88) τ | u,κ (cid:32) (cid:88) v ∈ S u d | v | / C | v | sup y ∈ (cid:101) B uκτ | ∂ v ˜ f u ( y ) | (cid:33) (cid:88) γ,γ (cid:48) | u (cid:89) j ∈ u (1 c j = c (cid:48) j − b − ) . (cid:80) γ,γ (cid:48) | u (cid:81) j ∈ u (1 c j = c (cid:48) j − b − ) = (2 − /b ) | u | . Thesupremum above is at most (cid:107) h u (cid:107) ∞ . From equation (19), we have | S u | = (2 d − | u | and also C | v | (cid:54) C d | u | for v ∈ S u . There are b | κ | indices τ in the sum given u and κ . From these considerations, σ uκ (cid:54) (cid:16) − b (cid:17) | u | b − | κ | /d (cid:107) h u (cid:107) ∞ (2 d − | u | d | u | C d | u | (cid:54) (cid:20) (cid:101) C (cid:16) − b (cid:17) (cid:21) | u | b − | κ | /d (cid:107) h u (cid:107) ∞ . Theorem 7.
Let u , . . . , u n be the points of a randomized ( t, m, s ) -net in base b .Let x i = φ ( u i ) ∈ X s for i = 1 , . . . , n where φ is the componentwise applicationof the transformation from convergent recursive splits in base b . Suppose as n → ∞ with t fixed, that all the gain coefficients of the net satisfy Γ uκ ≤ G < ∞ .Then for a smooth f on X s , Var(ˆ µ ) = O (cid:18) (log n ) s − n /d (cid:19) . If f ∈ C ds ( X s ) where each X ( j ) is a bounded closed set with non-empty interiorand a boundary of measure zero, then (20) holds regardless of whether X s isSobol’ extensible.Proof. We know from (10) thatVar(ˆ µ ) ≤ Gn (cid:88) | u | > (cid:88) | κ | > ( m − t −| u | ) + σ uκ ≤ Gn (cid:88) | u | > (cid:34) (cid:101) C (cid:18) − b (cid:19) (cid:35) | u | (cid:107) h u (cid:107) ∞ (cid:88) | κ | > ( m − t −| u | ) + b − | κ | /d ≤ (cid:101) Gn (cid:88) | u | > (cid:88) | κ | > ( m − t −| u | ) + b − | κ | /d where (cid:101) G = G (cid:34) ˜ c d (cid:18) − b (cid:19) (cid:35) | u | max | u | > (cid:107) h u (cid:107) ∞ . Since we are interested in the limit as m → ∞ , we may suppose that m > s + t .For such large m , we have (cid:88) | κ | > ( m − t −| u | ) + b − | κ | /d = ∞ (cid:88) r = m − t −| u | +1 b − r/d (cid:18) r + | u | − | u | − (cid:19) where the binomial coefficient is the number of | u | -vectors κ of nonnegative25ntegers that sum to r . Making the substitution s = r − m + t + | u | , (cid:88) | κ | > ( m − t −| u | ) + b − | κ | /d = b ( − m + t + | u | )2 /d ∞ (cid:88) s =1 b − s/d (cid:18) s + m − t − | u | − (cid:19) ≤ b ( t + | u | )2 /d n /d ( | u | − ∞ (cid:88) s =1 b − s/d ( s + m − t − | u |− = b ( t + | u | )2 /d n /d ( | u | − ∞ (cid:88) s =1 b − s/d | u |− (cid:88) j =0 (cid:18) | u | − j (cid:19) s j ( m − t − | u |− − j = b ( t + | u | )2 /d n /d | u |− (cid:88) j =0 ( m − t − | u |− − j j !( | u | − − j )! ∞ (cid:88) s =1 b − s/d s j ≤ b ( t + | u | )2 /d n /d m | u |− | u | ∞ (cid:88) s =1 b − s/d s | u |− . Note by the ratio test it is easy to see that (cid:80) ∞ s =1 b − s/d s | u |− converges. Alsoas m ≤ log b ( n ) and | u | ≤ s we get (cid:88) | κ | > ( m − t −| u | ) + b − | κ | /d = O (cid:18) (log n ) s − n /d (cid:19) . Plugging this back into the bound for the variance we get the desired result.
Our integration of smooth functions over an s -fold product of d -dimensionalspaces has root mean squared error (RMSE) of O ( n − / − /d (log( n )) ( s − / ).Plain QMC might map [0 , sd to R . If the composition of the integrand withsuch a mapping is in BVHK, then QMC attains an error rate of O ( n − log( n ) sd − ).Our mapping then has the advantage for d = 1 and 2. When the composition isnot in BVHK then QMC need not even converge to the right integral estimate.Then scrambled nets provide much needed assurance as well as error estimates.When the composed integrand is smooth, then scrambled nets applied di-rectly to [0 , sd would have an RMSE of O ( n − / log( n ) ( sd − / ). That is abetter asymptotic rate than we attain here, and it might really be descriptive offinite sample sizes even for very large sd , if the composite integrand were of loweffective dimension (Caflisch et al., 1997). If however, the composed integrandis in L but is not smooth, then scrambled nets applied in sd dimensions wouldhave an RMSE of o ( n − / ) but not necessarily better than that. Our proposalis then materially better for small d .The composed integrand that we actually use is not smooth on [0 , s . Itgenerally has discontinuities at all b -adic fractions t/b k for any of the componentsof u . For example in the four-fold split of Figure 1, an (cid:15) change in u can move26 point from the top triangle to the right hand triangle. These are howeveraxis-aligned discontinuities. Wang and Sloan (2011) call these QMC-friendlydiscontinuities. They don’t induce infinite variation.We have used nested uniform scrambles. The same results apply to otherscrambles, notably the linear scrambles of Matouˇsek (1998). Those scramblesare less space-demanding than nested uniform scrambles. A central limit theo-rem applies to averages over nested uniform scrambles (Loh, 2003), but has notbeen shown for linear scrambles. Hong et al. (2003) find that nested uniformscrambles have stochastically smaller values of a squared discrepancy measure.The splits we used allowed overlaps on sets of measure zero. We couldalso have relaxed X = ∪ b − a =0 X a to vol ( X \ ∪ b − a =0 X a ) = 0 . That could cause thedeterministic construction to fail to be a weak geometric ( t, m, s )-net but thescrambled versions would still be geometric ( t, m, s )-nets with probability one.Our main result was proved assuming that all d j = d . We can extend itto unequal d j by taking d = max j ∈ s d j . To make the extension, one can add d j − d ‘do nothing’ dimensions to X ( j ) . The splits never take place along thosedimensions, so the cells become cylinder sets and the function does not dependon the value of those components. We can make the extent of those do-nothingdimensions as small as we like to retain control of the diameter of the splits andthen apply Theorem 7. References
Arvo, J. (1995). Stratified sampling of spherical triangles. In
Proceedings ofthe 22nd annual conference on Computer graphics and interactive techniques ,pages 437–438. ACM.Arvo, J., Fajardo, M. Hanrahan, P., Jensen, H. W., Mitchell, D., Pharr, M.,and Shirley, P. (2001). State of the art in Monte Carlo ray tracing for realisticimage synthesis. In
ACM Siggraph 2001 , New York. ACM.Basu, K. (2014). Quasi-Monte Carlo tractability of high dimensional inte-gration over product of simplices. Technical report, Stanford University.arXiv:1411.0731.Basu, K. and Owen, A. B. (2014). Low discrepancy constructions in the triangle.Technical report, Stanford University. arXiv:1403.2649.Beckers, B. and Beckers, P. (2012). A general rule for disk and hemispherepartition into equal-area cells.
Computational Geometry , 45(2):275–283.Brandolini, L., Colzani, L., Gigante, G., and Travaglini, G. (2013). A Koksma–Hlawka inequality for simplices. In
Trends in Harmonic Analysis , pages 33–46.Springer.Caflisch, R. E., Morokoff, W., and Owen, A. B. (1997). Valuation of mort-gage backed securities using Brownian bridges to reduce effective dimension.
Journal of Computational Finance , 1:27–46.27ick, J. and Pillichshammer, F. (2010).
Digital sequences, discrepancy andquasi-Monte Carlo integration . Cambridge University Press, Cambridge.Hesse, K., Kuo, F. Y., and Sloan, I. H. (2007). A component-by-componentapproach to efficient numerical integration over products of spheres.
Journalof Complexity , 23(1):25–51.Hong, H., Hickernell, F. J., and Wei, G. (2003). The distribution of the dis-crepancy of scrambled digital (t,m,s)-nets.
Mathematics and Computers inSimulation , 62(3-6):335–345. 3rd { IMACS } Seminar on Monte Carlo Meth-ods.Keller, A. (2013). Quasi-Monte Carlo image synthesis in a nutshell. In Dick, J.,Kuo, F. Y., Peters, G. W., and Sloan, I. H., editors,
Monte Carlo and Quasi-Monte Carlo Methods 2012 , volume 65 of
Springer Proceedings in Mathemat-ics & Statistics , pages 213–249. Springer, Berlin.Kuo, F. Y. and Sloan, I. H. (2005). Quasi-Monte Carlo methods can be efficientfor integration over products of spheres.
Journal of Complexity , 21(2):196–210.L’Ecuyer, P. and Lemieux, C. (2002). A survey of randomized quasi-MonteCarlo methods. In Dror, M., L’Ecuyer, P., and Szidarovszki, F., editors,
Modeling Uncertainty: An Examination of Stochastic Theory, Methods, andApplications , pages 419–474. Kluwer Academic Publishers.Loh, W.-L. (2003). On the asymptotic distribution of scrambled net quadrature.
Annals of Statistics , 31(4):1282–1324.Matouˇsek, J. (1998).
Geometric Discrepancy : An Illustrated Guide . Springer-Verlag, Heidelberg.Niederreiter, H. (1987). Point sets and sequences with small discrepancy.
Monat-shefte fur mathematik , 104:273–337.Niederreiter, H. (1992).
Random Number Generation and Quasi-Monte CarloMethods . S.I.A.M., Philadelphia, PA.Owen, A. B. (1995). Randomly permuted ( t, m, s )-nets and ( t, s )-sequences. InNiederreiter, H. and Shiue, P. J.-S., editors,
Monte Carlo and Quasi-MonteCarlo Methods in Scientific Computing , pages 299–317, New York. Springer-Verlag.Owen, A. B. (1997a). Monte Carlo variance of scrambled equidistributionquadrature.
SIAM Journal of Numerical Analysis , 34(5):1884–1910.Owen, A. B. (1997b). Scrambled net variance for integrals of smooth functions.
Annals of Statistics , 25(4):1541–1562.Owen, A. B. (1998). Scrambling Sobol’ and Niederreiter-Xing points.
Journalof Complexity , 14(4):466–489. 28wen, A. B. (2003). Variance with alternative scramblings of digital nets.
ACMTransactions on Modeling and Computer Simulation , 13(4):363–378.Owen, A. B. (2005). Multidimensional variation for quasi-Monte Carlo. In Fan,J. and Li, G., editors,
International Conference on Statistics in honour ofProfessor Kai-Tai Fang’s 65th birthday .Owen, A. B. (2006).
Quasi-Monte Carlo for integrands with point singularitiesat unknown locations . Springer.Owen, A. B. (2008). Local antithetic sampling with scrambled nets.
The Annalsof Statistics , 36(5):2319–2343.Sobol’, I. M. (1973). Calculation of improper integrals using uniformly dis-tributed sequences.
Soviet Math Dokl , 14(3):734–738.Song, L., Kimerling, A. J., and Sahr, K. (2002). Developing an equal area globalgrid by small circle subdivision. In Goodchild, M. F. and Kimerling, A. J.,editors,
Discrete Global Grids . National Center for Geographic Information& Analysis, Santa Barbara, CA.Stromberg, K. R. (1994).
Probability for analysts . Chapman & Hall, New York.van der Corput, J. G. (1935a). Verteilungsfunktionen I.
Nederl. Akad. Wetensch.Proc. , 38:813–821.van der Corput, J. G. (1935b). Verteilungsfunktionen II.
Nederl. Akad. Weten-sch. Proc. , 38:1058–1066.Wang, X. and Sloan, I. H. (2011). Quasi-Monte Carlo methods in financialengineering: An equivalence principle and dimension reduction.
OperationsResearch , 59(1):80–95.Whitney, H. (1934). Analytic extensions of differentiable functions defined inclosed sets.