A Tool for Custom Construction of QMC and RQMC Point Sets
Pierre L'Ecuyer, Pierre Marion, Maxime Godin, Florian Puchhammer
AA Tool for Custom Construction of QMC andRQMC Point Sets
Pierre L’Ecuyer, Pierre Marion, Maxime Godin, and Florian Puchhammer
Abstract
We present LatNet Builder, a software tool to find good parameters forlattice rules, polynomial lattice rules, and digital nets in base 2, for quasi-MonteCarlo (QMC) and randomized quasi-Monte Carlo (RQMC) sampling over the s -dimensional unit hypercube. The selection criteria are figures of merit that givedifferent weights to different subsets of coordinates. They are upper bounds on theworst-case error (for QMC) or variance (for RQMC) for integrands rescaled to havea norm of at most one in certain Hilbert spaces of functions. Various Hilbert spaces,figures of merit, types of constructions, and search methods are covered by the tool.We provide simple illustrations of what it can do. QMC methods approximate an integral of the form µ = (cid:90) · · · (cid:90) f ( u , . . . , u s ) du · · · du s = (cid:90) ( , ) s f ( u ) d u = E [ f ( U )] (1) Pierre L’EcuyerD´epartement d’Informatique et de Recherche Op´erationnelle, Universit´e de Montr´eal, Canada,e-mail: [email protected] MarionSorbonne Universit´e, CNRS, Laboratoire de Probabilit´es, Statistique et Mod´elisation, LPSM,F-75005 Paris, France, e-mail: [email protected] GodinInstitut des Actuaires, France, e-mail: [email protected] PuchhammerUniversit´e de Montr´eal, Canada, and Basque Center for Applied Mathematics, Spain, e-mail:[email protected] 1 a r X i v : . [ s t a t . C O ] D ec Pierre L’Ecuyer, Pierre Marion, Maxime Godin, and Florian Puchhammer where f : ( , ) s → R , by the average¯ µ n = n n − ∑ i = f ( u i ) (2)where P n = { u , . . . , u n − } ⊂ [ , ) s is a set of n deterministic points that cover theunit hypercube more evenly than typical independent random points. That is, thediscrepancy between their empirical distribution and the uniform distribution over [ , ) s is smaller than for independent random points and converges to 0 faster than O ( n − / ) when n → ∞ . This discrepancy can be defined in many ways. It usuallyrepresents the worst-case integration error for a given class of integrands f . Typically,this class is a reproducing-kernel Hilbert space (RKHS) H of functions, such that | E n | : = | ¯ µ n − µ | ≤ D ( P n ) V ( f ) (3)for all f ∈ H , where V ( f ) is the norm of f − µ in H (we call it the variation of f )and D ( · ) is the discrepancy measure associated with this Hilbert space [8, 18, 39].For a fixed f ∈ H with V ( f ) >
0, the error bound in (3) converges at the same rate as D ( P n ) . A traditional version of (3), whose derivation does not involve Hilbert spaces,is the classical Koksma-Hlawka inequality, in which V ( f ) is the Hardy-Krausevariation and D ( P n ) is the star discrepancy, which converges as O (( log n ) s − n − ) forwell-selected point sets [39]. Another important choice for H is a Sobolev space offunctions whose mixed partial derivatives of order up to α are square-integrable. It isknown that for this space, one can construct point sets whose discrepancy convergesas O (( log n ) ( s − ) / n − α ) , and that this is the best possible rate [4, 8, 14, 15, 16, 17].The main classes of QMC methods are lattice points and digital nets.For RQMC, the n QMC points are randomized to provide a set of random points { U , . . . , U n − } ⊂ ( , ) s for which (i) each U i individually has the uniform distribu-tion over [ , ) s , and (ii) the points keep their highly-uniform distribution collectively.Randomizations that provably preserve the low discrepancy generally depend on thetype of QMC construction: some are used for lattice points and others for digitalnets. In some cases, the randomization may even improve the convergence rate ofthe mean square discrepancy. The RQMC estimatorˆ µ n , rqmc = n n − ∑ i = f ( U i ) , (4)which is now random, is unbiased for µ and one wishes to minimize its variance. Formore details on RQMC, see for example [21, 26, 29, 30, 37, 44, 45, 46].The aim of this paper is to introduce LatNet Builder , a software tool designedto construct specific point sets for QMC and RQMC, in any number of dimensions,for an arbitrary number of points, arbitrary weights on the subsets of coordinates,arbitrary smoothness of the integrands, a variety of construction and randomizationmethods, and several choices of discrepancies. The point sets can also be extensiblein the number of points and number of dimensions. By “constructing the points” here
Tool for Custom Construction of QMC and RQMC Point Sets 3 we mean defining the set P n by selecting the parameter values for a general structure,by trying to miminize a figure of merit (FOM) that may represent a discrepancy D ( P n ) or be an upper bound on it. Once this is done, other software can be used to randomizeand generate the points for their utilization in applications; see [27, 42] for example.LatNet Builder is available in open source at https://github.com/umontreal-simul/latnetbuilder. It is a descendant of Lattice Builder [32], whose scope was limited toordinary lattice rules. Another related tool is Nuyens’ fast CBC constructions [40].The rest of this paper is organized as follows. In Section 2, we recall the types ofQMC point sets covered by our software, namely ordinary lattice points, polynomiallattice points, digital nets, and their higher-order and interlaced versions, as well as themain randomization methods to turn these point sets into RQMC points. In Section 3,we give the general form of weighted RKHS used in this paper and the correspondinggeneralized Koksma-Hlawka inequality. We also recall the most common types ofweights. In Section 4, we give examples of discrepancies that are used as FOMsto select the parameters of point set constructions. In Section 5, we summarize thesearch methods implemented in our software. In Section 6, we compare FOM valuesobtained by various point set constructions and search methods. We also compareRQMC variance for simple integrands f . Section 7 gives a conclusion. LatNet Builder handles ordinary rank-1 lattice points as well as digital nets, whichinclude polynomial lattice rules and high-order and interlaced constructions.For a rank-1 lattice rule , the point set is P n = { u i = i v mod 1 , i = , . . . , n − } where n v = a = ( a , . . . , a s ) ∈ Z sn ≡ { , · · · , n − } s . It is a Korobov rule if a =( , a , a mod n , . . . , a s − mod n ) for some integer a ∈ Z n . The parameter to selecthere is the vector a , for any given n . The usual way to turn a lattice rule into anRQMC point set is by a random shift: generate a single random point U uniformlyin ( , ) s , and add it to each point of P n , modulo 1, coordinate-wise. This satisfiesthe RQMC conditions. For more details on lattice rules and their randomly-shiftedversions, see [19, 20, 26, 29, 31, 49].The Digital nets in base 2 handled by LatNet Builder are defined as follows. Thenumber of points is n = k for some integer k . We select an integer w ≥ k and s generating matrices C , . . . , C s of dimensions w × k and of rank k , with elementsin Z ≡ { , } . The points u i , i = , . . . , k −
1, are defined as follows: for i = a i , + a i , + · · · + a i , k − k − , we take u i , j , ... u i , j , w = C j a i , ... a i , k − mod 2 , u i , j = w ∑ (cid:96) = u i , j ,(cid:96) − (cid:96) , Pierre L’Ecuyer, Pierre Marion, Maxime Godin, and Florian Puchhammer and u i = ( u i , , . . . , u i , s ) . There are more general definitions in [8, 39]. The parametersto optimize are the elements of the matrices C j . Each one-dimensional projectiontruncated to its first k digits is Z n / n = { , / n , . . . , ( n − ) / n } . Ordinary digital netsconstructed by LatNet Builder have w = k , so the points have only k digits.The most popular digital net constructions are still the Sobol’ points [51], in base b =
2, with k × k generating matrices that are upper triangular and invertible. Thesematrices are constructed by a specific method, but the bits of the first few columnsabove the diagonal can be selected arbitrarily, and their choice has an impact onthe quality of the net. General-purpose choices have been proposed in [24, 34], e.g.,based on the uniformity of two-dimensional projections. LatNet Builder allows oneto construct the matrices based on a much more flexible class of criteria.A polynomial lattice rule (PLR) in base 2 is defined as follows. We denote by Z [ z ] the ring of polynomials with coefficients in Z , by L the set of formal seriesof the form ∑ ∞ (cid:96) = (cid:96) x (cid:96) z − (cid:96) with each x (cid:96) ∈ Z and (cid:96) ∈ Z , and we define ϕ : L → R by ϕ (cid:32) ∞ ∑ (cid:96) = (cid:96) x (cid:96) z − (cid:96) (cid:33) = w ∑ (cid:96) = max ( (cid:96) , ) x (cid:96) − (cid:96) . (5)We select a polynomial modulus Q = Q ( z ) ∈ Z [ z ] of degree k , and a generating vec-tor a ( z ) = ( a ( z ) , . . . , a s ( z )) ∈ Z [ z ] s , whose coordinates are polynomials of degreesless than k having no common factor with Q ( z ) . The point set of cardinality n = k is P n = (cid:26)(cid:18) ϕ (cid:18) h ( z ) a ( z ) Q ( z ) (cid:19) , . . . , ϕ (cid:18) h ( z ) a s ( z ) Q ( z ) (cid:19)(cid:19) : h ( z ) ∈ Z [ z ] and deg ( h ( z )) < k (cid:27) . (6)Here, we want to optimize the vector a ( z ) . This point set turns out to be a digital net inbase 2 whose generating matrices C j contain the first w digits of the binary expansionof the a j ( z ) / Q ( z ) . These are Hankel matrices: each row is the previous one shiftedto the left by one position, with the last entry determined by the recurrence withcharacteristic polynomial Q ( z ) , applied to the entries of the previous row. They canhave an arbitrary number of rows, but for ordinary PLRs, they are usually truncatedto w = k rows. See [8, 33, 35, 39, 38] for further details on PLRs.A high-order polynomial lattice rule (HOPLR) of order α with n = k points isobtained by constructing an ordinary PLR with polynomial modulus ˜ Q ( z ) of degree α k having 2 α k points in s dimensions, and using only the first n = k points. See[1, 2, 7]. This type of construction can achieve a higher order of convergence for theerror (almost O ( n − α ) ) than an ordinary PLR for integrands f in a Sobolev spaceof smoothness order α (i.e., when all mixed partial derivatives of up to order α aresquare integrable). One drawback is that because of the high degree of ˜ Q , the cost of afull CBC construction (see Section 5) is much higher since there are 2 α k possibilitiesto examine each time we select a new coordinate of the generating vector.Dick [3] also proposed an interlacing construction , for digital nets in general(which includes PLRs), that can provide the higher-order convergence rate of al-most O ( n − α ) for the integration error, for integrands with smoothness order α . Foran interlacing factor d ∈ N , the method first constructs a digital net in sd dimen- Tool for Custom Construction of QMC and RQMC Point Sets 5 sions, with generating matrices C , . . . , C sd . Then the generating matrices of the s -dimensional interlaced net are C ( d ) , . . . , C ( d ) s , where the rows of C ( d ) j are the firstrows of C ( j − ) d + , . . . , C jd in this order, then the second rows of these matrices in thesame order, and so on.The simplest way to define a RQMC point set from a digital net in base 2 isto add a digital random shift modulo 2 to all the points. To do this, we generate asingle point U = ( U , . . . , U s ) uniformly in ( , ) s , and perform a bitwise exclusive-or(XOR) between the binary digits of U and the corresponding digits of each point u i .A more involved randomization method for digital nets is the nested uniformscramble (NUS) of Owen [44, 45]. In base 2, for each coordinate, we do the following.With probability 1/2, flip the first bit of all the points. Then, for the points whose firstbit is 1, with probability 1/2, flip all the second bits. Do the same for the points whosefirst bit is 0, independently. Then do this recursively for all the bits. After all flippingis done for the first (cid:96) bits, partition the points in 2 (cid:96) batches according to the values oftheir first (cid:96) bits, and for each batch, flip bit (cid:96) + ( (cid:96) − ) s random bits to flip thefirst (cid:96) bits of all coordinates. One can equivalently do this only for the first k bits, andgenerate the other bits randomly and independently across points [37].A less expensive scramble, which gives less independence than NUS but more thana digital random shift, is a (left) linear matrix scramble (LMS) followed by a digitalrandom shift (LMS+shift) [21, 23, 37, 47]. The LMS replaces C j by L j C j mod 2,where L j is a random non-singular lower-triangular w × w binary matrix.Owen [45] proved that under sufficient smoothness conditions on f , the RQMCvariance with NUS on digital nets with fixed s and bounded t converges as O ( n − ( log n ) s − ) . A variance bound of the same order was shown for LMS+shift in[21, 53]. Note that these results were proved under the assumption that w = ∞ . The FOMs used by LatNet Builder are based on generalized (weighted) Koksma-Hlawka inequalities of the form (3) where V p ( f ) = ∑ /0 (cid:54) = u ⊆{ , ,..., s } γ − p u V p ( f u ) (7)and D q ( P n ) = ∑ /0 (cid:54) = u ⊆{ , ,..., s } γ q u D q u ( P n ) , (8)where 1 / p + / q = γ u ∈ R is a weight assigned to the subset u , V ( f u ) is thevariation of f u , D u ( P n ) is the discrepancy of the projection of P n over the subset u of coordinates, and f = ∑ u ⊆{ , ,..., s } f u is the functional ANOVA decomposition of f [10, 46]. LatNet Builder allows any q ∈ [ , ∞ ] . Taking q = ∞ with p = q and taking the max instead of the sum in (8), while p = ∞ with Pierre L’Ecuyer, Pierre Marion, Maxime Godin, and Florian Puchhammer q = p and taking the max instead of the sum in (7). The mostcommon choice is p = q = D u ( P n ) , depending on thepoint set constructions. Some of these measures correspond to the worst-case errorin some function space, assuming that the points of P n are not randomized. Otherscorrespond to the mean-square error (or variance), assuming that the points arerandomized in some particular way. This is typically done by defining a RKHSwith a kernel that is invariant with respect to the given randomization (i.e., digitalshift-invariant, scramble-invariant, etc.), and taking the worst-case error in that space.The role of the weights is to better recognize the importance of the subsets u forwhich f u contributes the most to the error or variance. That is, if V ( f u ) is unusuallylarge, we want to divide it by a larger weight γ u to control its contribution to V ( f ) ,but then we have to multiply D u ( P n ) in (8) by the same weight. The final effect isthat the FOM will penalize more the discrepancy for that particular projection.In principle, the weights γ u can be arbitrary. But for large s , defining arbitraryindividual weights for the 2 s − s − product weights , for which a weight γ j isassigned to coordinate j for j = , . . . , s , and γ u = ∏ j ∈ u γ j ; order-dependent weights ,for which γ u = Γ | u | where Γ , . . . , Γ s are selected constants and | u | is the cardinality ofof u ; and the product-and-order-dependent (POD) weights , which are a combinationof the two, with γ u = Γ | u | ∏ j ∈ u γ j . These are all available in LatNet Builder.LatNet Builder can construct point sets that are extensible in the number ofdimensions and also in the number of points, which means that we can constructpoint sets that perform well in the first s dimensions for s = s min , . . . , s max , and/orif we take the first n points for n = n , n , . . . , n m , simultaneously. Typically, onewould take n j = k min + j − for j = , . . . , m , so n m = k max = k min + m − [21]. Theglobal FOM in this case will be a weighted sum or maximum of the FOMs over theconsidered dimensions s and/or cardinalities n j . The CBC construction approachdescribed in Section 5 already gives a way to implement the extension in s . Forthe extension in n (or k ), LatNet Builder implements criteria and heuristic searchmethods that account for a global FOM. Most FOMs considered in LatNet Builder have the general form (8) where typically,when the points have the appropriate special structure of a lattice, polynomial lattice,or digital net, and with an adapted FOM, we have D q u ( P n ) = n n − ∑ i = ∏ j ∈ u φ ( u i , j ) (9) Tool for Custom Construction of QMC and RQMC Point Sets 7 for some function φ : [ , ) → R . With product weights γ u = ∏ j ∈ u γ j , the FOM thenbecomes D q ( P n ) = − + n n − ∑ i = s ∏ j = ( + γ j φ ( u i , j )) , which can be computed with O ( ns ) evaluations of φ .As an illustration, for a randomly-shifted lattice rule, the variance is:Var [ ˆ µ n , rqmc ] = ∑ (cid:54) = h ∈ L ∗ s | ˆ f ( h ) | , (10)where L ∗ s ⊂ Z s is the dual lattice [29]. It is also known that for periodic continuousfunctions having square-integrable mixed partial derivatives up to order α / α ≥
2, one has | ˆ f ( h ) | = O (( max ( , h ) · · · max ( , h s )) − α ) . Thismotivates the well-known FOM [31, 39, 49]: P α : = ∑ (cid:54) = h ∈ L ∗ s ( max ( , h ) · · · max ( , h s )) − α = − + n n − ∑ i = s ∏ j = (cid:32) − ( − π ) α / α ! B α ( u i , j ) (cid:33) (11) = n n − ∑ i = ∑ /0 (cid:54) = u ⊆{ ,..., s } (cid:32) − ( − π ) α / α ! (cid:33) | u | ∏ j ∈ u B α ( u i , j ) where B α / is the Bernoulli polynomial of degree α / B ( u ) = u − / B ( u ) = u − u + /
6, etc.), and the equality in (11) holds only when α is an even integer.Moreover, there are rank-1 lattices point sets P n for which P α converges almost as O ( n − α ) [9, 48, 49], i.e., as O ( n − α + ε ) for any ε >
0. Adding projection-dependentweights γ u leads to the weighted P γ , α , defined by (8) and (9) with q = φ ( u i , j ) = − ( − π ) α / B α ( u i , j ) / α ! , and D u ( P n ) = P α , u ( P n ) is the P α for the projection of P n on the coordinates in u .There is a similar variance expression for digital nets in base 2 with a randomdigital shift, with the Fourier coefficients ˆ f ( h ) replaced the the Walsh coefficients,and the dual lattice replaced by the dual net [8, Definition 4.76] or the dual latticein the case of PLRs [25, 35]. Thus, FOMs that correspond to variance boundscan be obtained by finding easily computable bounds on the Walsh coefficients.By assuming a rate of decrease of O ( − α | h | ) of the Walsh coefficients ˜ f ( h ) with h = ( h , . . . , h s ) ∈ N s and | h | = | h | + · · · + | h s | , and using a RKHS with shift- andscramble-invariant kernel, [53] and [6] obtain a FOM of the form (8) and (9) with φ ( x ) = φ α ( x ) = µ ( α ) − I [ x > ] · ( + (cid:98) log ( x ) (cid:99) )( α − ) ( µ ( α ) + ) , Pierre L’Ecuyer, Pierre Marion, Maxime Godin, and Florian Puchhammer where I is the indicator function and µ ( α ) = (cid:0) − − α (cid:1) − for any real number α >
1. This gives µ ( ) = µ ( ) = /
3, . . . , For α =
2, this gives φ ( x ) = ( − I [ x > ] · · (cid:98) log ( x ) (cid:99) ) , which corresponds to the FOM suggested in [35, Section 6.3] for PLRs. In [22, 53], φ ( x ) is written in terms of η = α − O ( n − α ( log n ) s − ) for any α >
1. This FOM can be seen as a counterpart of P α andwe call it ˜ P α . Its value ˜ P α , u ( P n ) on the projection of P n on the coordinates in u canbe used for D u ( P n ) , with q = q = φ ( x ) = / − I [ x > ] · (cid:98) log ( x ) (cid:99)− . They show that there are digital nets for which this FOM (and therefore the squareerror) converges almost as O ( n − ) . In their Chapter 13, they find that the scramble-invariant version gives the same φ .Goda [12] examines interlaced polynomial lattice rules (IPLR), also for a Sobolevspace of order α , with an interlacing factor d ≥
1. He provides two upper bounds onthe worst-case error in a deterministic setting. These bounds can be used as FOMs.The first is valid for all positive integer values of α and δ , whereas the second holdsonly for d ≤ α , but is tighter when it applies. These two bounds have the form (8)and (9) with q = γ u replaced by ˜ γ u , and φ ( x i , j ) = − + d ∏ (cid:96) = ( + φ α , d ,(cid:96) ( x i , ( j − ) d + (cid:96) )) , where for the first bound, ˜ γ u = γ u α ( d − ) | u | / , φ α , d ,(cid:96) ( x ) = − ( min ( α , d ) − ) (cid:98) log x (cid:99) ( min ( α , d ) − ) ( α + ) / ( min ( α , d ) − − ) for all x ∈ [ , ) , where 2 (cid:98) log (cid:99) =
0, while for the second bound, ˜ γ u = γ u and φ α , d ,(cid:96) ( x ) = d − ( − ( d − ) (cid:98) log x (cid:99) ( d − )) (cid:96) ( d − − ) . One can achieve a convergence rate of almost O ( n − min ( α , d ) ) for these FOMs (andtherefore for the worst-case error).Goda and Dick [13] proposed another FOM, also for a Sobolev space of order α , for interlaced randomly-scrambled PLRs of high order. They showed that thisscheme can achieve the best possible convergence rate of O ( n − ( ( α , d )+ )+ δ ) forthe variance. The FOM, which we denote by I α , has the same form, but with q = Tool for Custom Construction of QMC and RQMC Point Sets 9 φ ( x ) = φ α , d ( x ) = − ( α , d ) (cid:98) log x (cid:99) ( ( α , d )+ − ) − − ( α , d ) , and γ u replaced by ˜ γ u = γ u D | u | α , d where D α , d = ( d − α , )+( d − ) α .Another set of FOMs are obtained from upper bounds on the star discrepancy of D ∗ ( P n ) or its projections on subsets of coordinates, when P n is a digital ( t , k , s ) -net.One such bound is D ∗ ( P n ) ≤ − ( − / n ) s + R where R = − + n n − ∑ i = s ∏ j = (cid:34) n − ∑ k = −(cid:98) log k (cid:99)− wal k ( u i , j ) (cid:35) (12)wal k is the k th Walsh function in one dimension, and we assume that the generatingmatrices C j are k × k . See [8, Theorems 5.34 and 5.36], where a more general versionwith projection-dependent weights is also given. For PLRs in base b =
2, this criterionis equal to R (cid:48) , γ given in [8, Chapter 10]: R (cid:48) , γ = − ∑ /0 (cid:54) = u ⊆{ , ,..., s } γ u + n n − ∑ i = ∑ /0 (cid:54) = u ⊆{ , ,..., s } γ u ∏ j ∈ u φ k ( u i , j ) (13)where φ k ( u ) = −(cid:98) log ( u ) (cid:99) / u ≥ − k and φ k ( u ) = + k / t -value of thedigital net: D ∗ ( P n ) ≤ t n s − ∑ j = (cid:18) k − tj (cid:19) . If we use this upper bound for each projection on the subset u of coordinates, we getthe FOM (8) with q = D u ( P n ) = t u n | u |− ∑ j = (cid:18) k − t u j (cid:19) where t u = t u ( P n ) is the t -value of the projection of P n on the coordinates in u . LatNetBuilder implements this with arbitrary weights, using algorithms described in [36].Dick [4] obtains worst-case error bounds that converge at rate almost O ( n − α ) forinterlaced digital nets, based on the t -values of the projections. For given construction type, FOM, and weights, finding the best choice of parametersmay require to try all possibilities, but their number is usually much too large. LatNetBuilder implements the following search methods.
In an exhaustive search , all choices of parameters are tried, so we are guaranteedto find the best one. This is possible only when there are not too many possibilities.A random search samples uniformly and independently a fixed number r ofparameter choices, and the best one is retained.In a full component-by-component (CBC) construction, the parameters are selectedfor one coordinate at a time, by taking into account the choices for the previouscoordinates only [5, 8, 50]. The parameters for coordinate j (e.g., the j th coordinateof the generating vector in the case of lattices), are selected by minimizing theFOM for the first j coordinates, in j dimensions, by examining all possibilitiesof parameters for this coordinate, without changing the parameter choices for theprevious coordinates. This is done for the s coordinates in succession. This greedyapproach can reduce by a huge factor (exponential in the dimension) the total numberof cases that are examined in comparison with the exhaustive search. What is veryinteresting is that for most types of QMC constructions and FOMs implemented inLatNet Builder, the convergence rate for the worst-case error or variance obtainedwith this restricted approach is provably the same as for the exhaustive search [8, 14].A random CBC construction can be used when the number of choices for eachcoordinate is too large: one examines only a fixed number of random choices foreach j .For lattice-type point sets, with certain FOMs, a fast CBC construction can beimplemented by using a fast Fourier transform (FFT), so the full CBC constructioncan be performed much faster [8, 43, 41]. LatNet Builder supports this.For lattice-type constructions, one can also further restrict the search to Korobov -type generating vectors. The first coordinate is set to 1 and only the second coordinateneeds to be selected. This can be done either by an exhaustive search or by just takinga random sample for the second coordinate ( random Korobov ).For digital nets, a mixed CBC method is also available: it uses full CBC for thefirst d − ≤ d ≤ s . Here we give a few simple examples of what LatNet Builder can do. The simulationexperiments, including the generation and randomization of the points, were doneusing SSJ [28].
One might be interested in estimating the probability distribution of FOM valuesobtained when selecting parameters at random for a given type of construction,perhaps under some constraints, and as a function of n . Here we estimate thisdistribution by its empirical counterpart with an independent sample of size 1000 Tool for Custom Construction of QMC and RQMC Point Sets 11 (with replacement), and we report the 0.1, 0.5 and 0.9 quantiles of this empiricaldistribution, for n going from 2 to 2 . We do this for PLRs, Sobol’ points, and digitalnets with arbitrary invertible and projection-regular generating matrices (randomnets), with ˜ P taken as the FOM, in s = γ u = . | u | for all u . Wealso report the value obtained by a (full) fast CBC search for a PLR. The results aredisplayed in the first panel of Figure 1. We see that the FOM distribution has a smallermean and much less variance for the Sobol’ points than for the other contructions.Even the median obtained for Sobol’ beats (slightly) the FOM obtained by a fullCBC construction with PLRs. The quantiles for random PLRs and random nets areapproximately the same.The second panel of the figure shows the results of a similar sampling for PLRswith I as a FOM, also in 6 dimensions. Note that this FOM applies only to PLRs.Here, the FOM values are more dispersed and the fast CBC gives a significantlybetter value than the best FOM obtained by random sampling. Also the search forthe point set parameters is much quicker with the fast CBC construction than withrandom sampling of size 1000.Fig. 1: The 0.1, 0.5, and 0.9 quantiles of the FOM distribution as functions of n forvarious constructions, in log-log scale. Number of points10 − − − − F O M v a l u e Convergence of P quantiles, s = PLR, 10%PLR, medianPLR, 90%PLR CBCSobol, 10%Sobol, medianSobol, 90%Random, 10%Random, medianRandom, 90% Number of points10 − − − F O M v a l u e Convergence of I quantiles, s = PLR, best among 1000PLR, 10%PLR, medianPLR, 90%PLR CBC
We now give small examples showing how searching for custom parameter valueswith LatNet Builder can make a difference in the RQMC variance compared withpre-tabulated parameter values available in software or over the Internet. We dothis for Sobol’ nets, and our comparison is with the precomputed direction numbersobtained by Joe and Kuo [24], which are arguably the best proposed values so far.These parameters were obtained by optimizing a FOM based on the t -values overtwo-dimensional projections, using a CBC construction. With LatNet Builder, wecan account for any selected projections in our FOM. For instance, if we think allthe projections in two and three dimensions are important, we can select a FOM thataccounts for all these projections. To illustrate this, we made a CBC constructionof n = Sobol’ points in s =
15 dimensions, using the sum or the maximum of t -values in the two- and three-dimensional projections. Figure 2 shows the distributionof t -values obtained with the sum, the max, and the points from [24]. Compared withthe latter, we are able to reduce the worse t -value over 3-dim projections from 8 to 5when using the max, and to reduce the average t -value when using the sum. However,when using the max, we get a few poor two-dim projections, because we comparethe t -values on the same scale for two and three dimensions. We should probablymultiply the t -value by a scaling factor that decreases with the dimension.Fig. 2: Distributions of t -values for 2-dim and 3-dim projections, for three Sobol’point sets: (1) Joe-Kuo taken from [24], (2)
Max and (3)
Sum are found by LatNetBuilder as explained in the text. For each case, we report the number of projectionshaving any given t -value, as well as the average t -value (dashed vertical lines). N u m b e r o f p r o j e c t i o n s Histogram of 2D projections per t-value
MaxSumJoe-Kuo 0 1 2 3 4 5 6 7 8t-value050100150200 N u m b e r o f p r o j e c t i o n s Histogram of 3D projections per t-value
MaxSumJoe-Kuo
In our next illustration, we compare the RQMC variances for Sobol’ points withdirection numbers taken from [24] and direction numbers found by LatNet Builderusing a custom FOM for our function. We want to integrate f ( u ) = ∏ j = ( ψ ( u j ) − µ ) + ∏ j = ( ψ ( u j ) − µ ) , Tool for Custom Construction of QMC and RQMC Point Sets 13 where ψ ( u ) = (cid:0) ( u − . ) + . (cid:1) − and µ = E [ ψ ( U )] ≈ . U ∼ U ( , ) .This function is the sum of two five-dimensional ANOVA terms for a more generalfunction taken from [11]. A good FOM for this function should focus mainly onthese two five-dim projections, namely u = { , , , , } and u = { , , , , } , andnot on the two-dim projections as in [24]. So we made a search with the ˜ P criterionwith weights γ u = n = Sobol’ points in 10 dimensions. Then we estimated the varianceof the sample RQMC average over these n points with the two choices of directionnumbers (those of [24] and ours), using m =
200 independent replications of anRQMC scheme that used only a random digital shift. The empirical variance withour custom points was smaller by a factor of more than 18.
Here we consider a family of test functions similar to those in [52]: f s , c ( u ) = s ∏ j = ( + c j · ( u j − / )) for u ∈ ( , ) s , where c = ( c , . . . , c s ) ∈ ( , ) s . The ANOVA components are, for all u ⊂ { , . . . , s } , ( f s , c ) u ( u ) = ∏ j ∈ u c j · ( u j − / )) , For an experiment, we take arbitrarly s = c = ( . , . , . ) . We use LatNetBuilder to search for good PLRs with a fast CBC construction, with product weights γ j = c j , with the FOMs P , I , and I (whose interlacing factors d are 1, 2, and3, respectively). For each n = k , k = , . . . ,
18, we estimate the RQMC variancewith m independent replications of the randomization scheme, with m = m =
100 for NUS. For the interlaced points, the randomization isperformed before the interlacing, as in [13]. Figure 3 shows the variance as a functionof n , in log-log scale. We see that the two randomization schemes give approximatelythe same variance. However, the time to generate and randomize the points is muchlarger for NUS than for LMS+shift: around 10 times longer for 2 points and 50times longer for 2 points. As expected, the variance reduction and the convergencerate are larger when the interlacing factor increases, although the curves are morenoisy. Fig. 3: Variance as a function of n in log-log scale, for PLRs with two randomizationschemes and three interlacing factors d , found with LatNet Builder.We also report the average time to generate and randomize the points. Number of points10 − − − − − − − V a r i a n c e LMS+shift - d=1NUS - d=1LMS+shift - d=2NUS - d=2LMS+shift - d=3NUS - d=3
Nb. of points Constr. time128 0.02 s512 0.04 s2048 0.12 s8192 0.50 s32768 2.25 s131072 8.32 s
LatNet Builder is both a tool for researchers to study the properties of highly uniformpoint sets and associated figures of merit, and for practitioners who want to find goodparameters for a specific task. It is relatively easy to incorporate new FOMs into thesoftware, especially if they are in the kernel form (9).Many questions remain open regarding the roles of the construction, the searchmethod, the randomization, and (perhaps more importantly) the choice of the weights.It is our hope that the software presented here will spur interest into these issues.
Acknowledgements
This work has been supported by a NSERC Discovery Grant and an IVADOGrant to P. L’Ecuyer, and by a stipend from Corps des Mines to P. Marion. F. Puchhammer wassupported by Spanish and Basque governments fundings through BCAM (ERDF, ESF, SEV-2017-0718, PID2019-108111RB-I00, PID2019-104927GB-C22, BERC 2018e2021, EXP. 2019/00432,KK-2020/00049), and the computing infrastructure of i2BASQUE and IZO-SGI SGIker (UPV).
References
1. Baldeaux, J., Dick, J., Greslehner, J., Pillichshammer, F.: Construction algorithms for higherorder polynomial lattice rules. Journal of Complexity , 281–299 (2011)2. Baldeaux, J., Dick, J., Leobacher, G., Nuyens, D., Pillichshammer, F.: Efficient calculationof the worst-case error and (fast) component-by-component construction of higher orderpolynomial lattice rules. Numerical Algorithms (3), 403–431 (2012)3. Dick, J.: Explicit constructions of quasi-Monte Carlo rules for the numerical integration ofhigh-dimensional periodic functions. SIAM Journal on Numerical Analysis (5), 2141–2176(2007)4. Dick, J.: Walsh spaces containing smooth functions and quasi-Monte Carlo rules of arbitraryhigh order. SIAM Journal on Numerical Analysis (3), 1519–1553 (2008) Tool for Custom Construction of QMC and RQMC Point Sets 155. Dick, J., Kuo, F.Y., Pillichshammer, F., Sloan, I.H.: Construction algorithms for polynomiallattice rules for multivariate integration. Mathematics of Computation (248), 1895–1921(2005)6. Dick, J., Pillichshammer, F.: Multivariate integration in weighted Hilbert spaces based on Walshfunctions and weighted Sobolev spaces. Journal of Complexity (2), 149–195 (2005)7. Dick, J., Pillichshammer, F.: Strong tractability of multivariate integration of arbitrary highorder using digitally shifted polynomial lattice rules. Journal of Complexity , 436–453(2007)8. Dick, J., Pillichshammer, F.: Digital Nets and Sequences: Discrepancy Theory and Quasi-MonteCarlo Integration. Cambridge University Press, Cambridge, U.K. (2010)9. Dick, J., Sloan, I.H., Wang, X., Wo´zniakowski, H.: Good lattice rules in weighted Korobovspaces with general weights. Numerische Mathematik , 63–97 (2006)10. Efron, B., Stein, C.: The jackknife estimator of variance. Annals of Statistics , 586–596 (1981)11. Genz, A.: Testing multidimensional integration routines. In: Proceedings of the InternationalConference on Tools, Methods and Languages for Scientific and Engineering Computation, pp.81–94. Elsevier North-Holland (1984)12. Goda, T.: Good interlaced polynomial lattice rules for numerical integration in weighted Walshspaces. Journal of Computational and Applied Mathematics , 279–294 (2015)13. Goda, T., Dick, J.: Construction of interlaced scrambled polynomial lattice rules of arbitraryhigh order. Foundation of Computational Mathematics , 1245–1278 (2019)14. Goda, T., Suzuki, K.: Recent advances in higher order quasi-Monte Carlo methods, pp. 69–102.De Gruyter (2019)15. Goda, T., Suzuki, K., Yoshiki, T.: An explicit construction of optimal order quasi-Monte Carlorules for smooth integrands. SIAM Journal on Numerical Analysis , 2664––2683 (2016)16. Goda, T., Suzuki, K., Yoshiki, T.: Optimal order quasi-Monte Carlo integration in weightedSobolev spaces of arbitrary smoothness. IMA Journal on Numerical Analysis , 505–518(2017)17. Goda, T., Suzuki, K., Yoshiki, T.: Optimal order quadrature error bounds for infinite-dimensional higher-order digital sequences. Foundation of Computational Mathematics ,433–458 (2018)18. Hickernell, F.J.: A generalized discrepancy and quadrature error bound. Mathematics ofComputation (221), 299–322 (1998)19. Hickernell, F.J.: Lattice rules: How well do they measure up? In: P. Hellekalek, G. Larcher(eds.) Random and Quasi-Random Point Sets, Lecture Notes in Statistics , vol. 138, pp. 109–166.Springer-Verlag, New York (1998)20. Hickernell, F.J.: Obtaining O ( N − + ε ) convergence for lattice quadrature rules. In: K.T. Fang,F.J. Hickernell, H. Niederreiter (eds.) Monte Carlo and Quasi-Monte Carlo Methods 2000, pp.274–289. Springer-Verlag, Berlin (2002)21. Hickernell, F.J., Hong, H.S., L’Ecuyer, P., Lemieux, C.: Extensible lattice sequences for quasi-Monte Carlo quadrature. SIAM Journal on Scientific Computing (3), 1117–1138 (2001)22. Hickernell, F.J., Yue, R.X.: The mean square discrepancy of scrambled ( t , s ) -sequences. SIAMJournal on Numerical Analysis (4), 1089–1112 (2001)23. Hong, H.S., Hickernell, F.H.: Algorithm 823: Implementing scrambled digital sequences. ACMTransactions on Mathematical Software , 95–109 (2003)24. Joe, S., Kuo, F.Y.: Constructing Sobol sequences with better two-dimensional projections.SIAM Journal on Scientific Computing (5), 2635–2654 (2008)25. L’Ecuyer, P.: Polynomial integration lattices. In: H. Niederreiter (ed.) Monte Carlo and Quasi-Monte Carlo Methods 2002, pp. 73–98. Springer-Verlag, Berlin (2004)26. L’Ecuyer, P.: Quasi-Monte Carlo methods with applications in finance. Finance and Stochastics (3), 307–349 (2009)27. L’Ecuyer, P.: SSJ: Stochastic simulation in Java (2016). http://simul.iro.umontreal.ca/ssj/28. L’Ecuyer, P., Buist, E.: Simulation in Java with SSJ. In: Proceedings of the 2005 WinterSimulation Conference, pp. 611–620. IEEE Press, Piscataway, NJ (2005)29. L’Ecuyer, P., Lemieux, C.: Variance reduction via lattice rules. Management Science (9),1214–1235 (2000)6 Pierre L’Ecuyer, Pierre Marion, Maxime Godin, and Florian Puchhammer30. L’Ecuyer, P., Lemieux, C.: Recent advances in randomized quasi-Monte Carlo methods. In:M. Dror, P. L’Ecuyer, F. Szidarovszky (eds.) Modeling Uncertainty: An Examination of Stochas-tic Theory, Methods, and Applications, pp. 419–474. Kluwer Academic, Boston (2002)31. L’Ecuyer, P., Munger, D.: On figures of merit for randomly-shifted lattice rules. In:H. Wo´zniakowski, L. Plaskota (eds.) Monte Carlo and Quasi-Monte Carlo Methods 2010,pp. 133–159. Springer-Verlag, Berlin (2012)32. L’Ecuyer, P., Munger, D.: Algorithm 958: Lattice builder: A general software tool for con-structing rank-1 lattice rules. ACM Transactions on Mathematical Software (2), Article 15(2016)33. L’Ecuyer, P., Touzin, R.: On the Deng-Lin random number generators and related methods.Statistics and Computing ∼ clemieux/randqmc.html35. Lemieux, C., L’Ecuyer, P.: Randomized polynomial lattice rules for multivariate integrationand simulation. SIAM Journal on Scientific Computing (5), 1768–1789 (2003)36. Marion, P., Godin, M., L’Ecuyer, P.: An algorithm to compute the t -value of a digital net and ofits projections. Journal of Computational and Applied Mathematics (June), 112669 (2020)37. Matousˇek, J.: On the L -discrepancy for anchored boxes. J. of Complexity , 527–556 (1998)38. Niederreiter, H.: Low-discrepancy point sets obtained by digital constructions over finite fields.Czechoslovak Math. Journal , 143–166 (1992)39. Niederreiter, H.: Random Number Generation and Quasi-Monte Carlo Methods, SIAM CBMS-NSF Reg. Conf. Series in Applied Mathematics , vol. 63. SIAM (1992)40. Nuyens, D.: Fast component-by-component constructions. http://people.cs.kuleuven.be/ ∼ dirk.nuyens/fast-cbc/ (2012). URL http://people.cs.kuleuven.be/ ∼ dirk.nuyens/fast-cbc/41. Nuyens, D.: The construction of good lattice rules and polynomial lattice rules. In: P. Kritzer,H. Niederreiter, F. Pillichshammer, A. Winterhof (eds.) Uniform Distribution and Quasi-MonteCarlo Methods: Discrepancy, Integration and Applications, pp. 223–255. De Gruyter (2014)42. Nuyens, D.: The magic point shop (2020). https://people.cs.kuleuven.be/ ∼ dirk.nuyens/qmc-generators/43. Nuyens, D., Cools, R.: Fast algorithms for component-by-component construction of rank-1lattice rules in shift-invariant reproducing kernel Hilbert spaces. Mathematics of Computation , 903–920 (2006)44. Owen, A.B.: Monte Carlo variance of scrambled equidistribution quadrature. SIAM Journal onNumerical Analysis (5), 1884–1910 (1997)45. Owen, A.B.: Scrambled net variance for integrals of smooth functions. Annals of Statistics (4), 1541–1562 (1997)46. Owen, A.B.: Latin supercube sampling for very high-dimensional simulations. ACM Transac-tions on Modeling and Computer Simulation (1), 71–102 (1998)47. Owen, A.B.: Variance with alternative scramblings of digital nets. ACM Transactions onModeling and Computer Simulation (4), 363–378 (2003)48. Sinescu, V., L’Ecuyer, P.: Variance bounds and existence results for randomly shifted latticerules. Journal of Computations and Applied Mathematics , 3296–3307 (2012)49. Sloan, I.H., Joe, S.: Lattice Methods for Multiple Integration. Clarendon Press, Oxford (1994)50. Sloan, I.H., Rezstov, A.: Component-by-component construction of good lattice rules. Mathe-matics of Computation , 262–273 (2002)51. Sobol’, I.M.: The distribution of points in a cube and the approximate evaluation of integrals.U.S.S.R. Comput. Math. and Math. Phys. (4), 86–112 (1967)52. Sobol, I.M., Asotsky, D.I.: One more experiment on estimating high-dimensional integrals byquasi-Monte Carlo methods. Mathematics and Computers in Simulation (3-6), 255–263(2003)53. Yue, R.X., Hickernell, F.J.: The discrepancy and gain coefficients of scrambled digital nets.Journal of Complexity18