Extensible Lattice Sequences for Quasi-Monte Carlo Quadrature.
Integration lattices are one of the main types of low discrepancy sets used in quasi-Monte Carlo methods.  However, they have the disadvantage of being of fixed size.  This article describes the construction of an infinite sequence of points, the first bm of which forms a lattice for any nonnegative integer m.  Thus, if the quadrature error using an initial lattice is too large, the lattice can be extended without discarding the original points.  Generating vectors for extensible lattices are found by minimizing a loss function based on some measure of discrepancy or nonuniformity of the lattice.  The spectral test used for finding pseudorandom number generators is one important example of such a discrepancy.  The performance of the extensible lattices proposed here is compared to that of other methods for some practical quadrature problems.
Introduction
. Multidimensional integrals appear in a wide variety of applications
in finance [4, 48, 50], physics and engineering [29, 42, 49, 58], and statistics
[8, 12, 13]. The integration domain may often be assumed, after some appropriate
transformation, to be the unit cube, in which case the integral takes the form:
Z
for some known integrand, f .
Adaptive methods, such as [2], have been developed for approximating multidimensional
integrals, but their performance deteriorates as the dimension increases.
For finance problems the dimension can be in the hundreds or even thousands. An
alternative to adaptive quadrature is Monte Carlo methods, where the integral is
approximated by the sample mean of the integrand evaluated on a set, P , of N independent
random points drawn from a uniform distribution on [0;
The quadrature error for Monte Carlo methods is typically O(N \Gamma1=2 ). One reason for
this relatively low accuracy is that the points in P are chosen independently of each
other. Thus, some parts of the integration domain contain clumps of points while
other parts are empty of points.
To obtain greater accuracy one may replace the random set P by a carefully
chosen deterministic set that is more uniformly distributed on [0; 1) s . As is explained
in Section 3, one may define a discrepancy that measures how much the empirical
distribution function of P differs from the continuous uniform distribution. Then
y Department of Mathematics, Hong Kong Baptist University, Kowloon Tong, Hong Kong SAR,
China. This research was supported by a HKBU FRG grant 96-97/II-67.
z fred@hkbu.edu.hk, http://www.math.hkbu.edu.hk/~fred
x D'epartement d'Informatique et de Recherche Op'erationnelle, Universit'e de Montr'eal C.P. 6128,
Succ. Centre-Ville, Montr'eal (Qu'ebec), Canada, H3C 3J7
F. J. Hickernell, H. S. Hong, P. L' '
Ecuyer and C. Lemieux.
one chooses P in quadrature rule (1.1) to have as small a discrepancy as possible.
The quadrature methods based on low discrepancy sets are called quasi-Monte Carlo
methods. They are discussed in several review articles [3, 16, 40, 59] and monographs
[27, 44, 53].
Two important families of low discrepancy sets are:
i. integration lattices [44, Chap. 5] and [53], and
ii. digital nets and sequences [44, Chap. 4].
These two families are introduced in Section 2. One advantage of the second family
is that any number of consecutive points from a good digital sequence has low dis-
crepancy. If one needs more points, one may use additional terms from the digital
sequence without discarding the original ones. On the other hand, until now, the
number of points in an integration lattice has had to be specified in advance. So far,
there has been no systematic way of adding points to an integration lattice while still
retaining its lattice structure.
The purpose of this article is to provide a method for constructing infinite lattice
sequences, thereby eliminating the need to know N , the number of points, in advance.
Although the emphasis is on rank-1 lattices, the method may be applied to integration
lattices of arbitrary rank. Given an infinite lattice sequence one may approximate a
multidimensional integrand with a quadrature rule of the form (1.1) for moderate
number of points N 0 . If the error estimate is unacceptably high, then one may choose
an additional points from the lattice sequence to obtain a quadrature rule
with N 1 points, and so on.
The following section describes the new method for obtaining infinite lattice se-
quences. Section 3 briefly reviews some results on discrepancy and quadrature error
analysis for quasi-Monte Carlo methods. These are used to find the generating vectors
for the new lattice sequences in Section 4. The issue of error estimation is addressed
in Section 5. Two practical examples are explored in Section 6, where the new lattice
sequences are compared with existing quadrature methods. The last section contains
some concluding remarks.
2. Integration Lattices and Digital Sequences. This section begins by introducing
integration lattices. Next, digital sequences and (t; s)-sequences are de-
scribed. Finally, the idea underlying digital sequences is used to produce infinite
lattice sequences.
2.1. Integration Lattices. Rank-1 lattices, also known as good lattice point
(glp) sets, were introduced by Korobov [31] and have been widely studied since then
(see [27, 44, 53] and the references therein). The formula for a shifted rank-1 lattice
set is simply
where N is the number of points, h is an s-dimensional generating vector of integers
(a good lattice point) that depends on N , \Delta is an s-dimensional shift vector in [0; 1) s ,
and fxg denotes the fractional part of a vector x, i.e.
later generalized glp sets by introducing more than one generating vector. A shifted
integration lattice with points based on generating vectors h
is:
Integration lattices and their use for quadrature are discussed in the monograph [53].
Extensible Lattice Sequences 3
For a given N there is the problem of choosing good generating vectors. Although
theoretical constructions exist for in higher dimensions one typically finds
generating vectors by minimizing a discrepancy or measure of non-uniformity of the
lattice. Several examples of discrepancies are given in Section 3.
2.2. Digital Nets and Sequences. Digital nets and sequences are another
method of constructing low discrepancy sets (see [32] and [44, Chapter 4]). Let b
denote a positive integer greater than one. For any non-negative integer one may
extract the digits of its base b representation, finitely
many of which are nonzero:
The i th term of a digital net or sequence is given by
z
z (i)
s: (2.2d)
If the generating matrices C are m \Theta m, then this construction yields a
digital net fz points. If the generating matrices are
1 \Theta 1, i.e. each C j has entries c jkl defined for k; then one has a digital
sequence g.
The prototype digital sequence is the one-dimensional Van der Corput sequence,
g. This is defined by taking equal to the identity
In essence, the Van der Corput sequence, takes the b-ary representation of an integer
and reflects it about the decimal point.
2.3. (t; m; s)-Nets and (t; s)-Sequences. Similarly to integration lattices one
has the problem of how to choose the generating matrices C j in (2.2). Usually this
is done to optimize the quality factor of the net or sequence. For any non-negative
s-vector k, and for a base b consider the following set of disjoint boxes, whose union
is the unit cube:
Each such box in B k has volume b \Gammak 1 \Gamma\Delta\Delta\Delta\Gammak s . A (t; m; s)-net in base b is a set of
points in [0; 1) s , such that every box in B k contains b m\Gammak 1 \Gamma\Delta\Delta\Delta\Gammak s of these
points for any k satisfying t. Thus, any function that is piecewise
4 F. J. Hickernell, H. S. Hong, P. L' '
Ecuyer and C. Lemieux.
constant on the boxes in B k will be integrated exactly according to quadrature rule
The integer parameter t is called the quality parameter of the net and it takes values
between 0 and m. A smaller value of t means a better net.
A (t; s)-sequence is an infinite sequence of points in [0; 1) s such that the b m points
numbered lb m to (l +1)b non-negative integer
l. By using a (t; s)-sequence to do quasi-Monte Carlo calculations, one need not know
the number of points required in advance. If the first b m1 points do not give sufficient
accuracy, then one may add the next b m2 \Gamma b m1 points in the sequence to get a net
with b m2 points, without throwing away the first b m1 points.
There is a connection between digital nets and sequences as defined above and
(t; m; s)-nets and (t; s)-sequences [32]. Let c (m)
jiffl denote the vector containing the first
m elements of the i th row of the the generating matrix C j . Given a non-negative
integer s-vector k, let C(m; k) be the following set of the first k j rows of the j th
generating matrix for s:
c (m)
Furthermore, let
is a linearly independent set of vectors
for all k with
The following theorem (see [32]) gives the condition for which a digital net is a (t; m; s)-
net and a digital sequence is a (t; s)-sequence.
Theorem 2.1. For prime bases b the digital net defined above in (2.2) is a
s)-net. If, in addition, T (s) is finite, then the digital sequence defined in
(2.2) is a (T (s); s)-sequence.
Finding good generating matrices for digital nets and sequences is an active area
of research. Virtually all generators found so far have been based on number theoretic
arguments. Early sequences include those of Sobol' [56], Faure [9] and Niederreiter
[43]. Algorithms for these sequences can be found in the ACM Transactions on Mathematical
Software collection. The FINDER software developed at Columbia University
by Traub and Papageorgiou implements generalized Sobol' and generalized Faure se-
quences. New constructions with smaller t values are given by Niederreiter and Xing
[47].
2.4. Infinite Lattice Sequences. The idea underlying digital sequences may
be extended to integration lattices to obtain infinite lattice sequences. The i th term
of a rank-1 lattice, which is fih=Ng, depends inherently on the number of points, N .
Thus, the formula for a lattice must be rewritten in a way that does not involve N
explicitly. A way to do this was first suggested in [24].
Suppose that the number of points, N , is some integer power of a base b - 2,
that is, . This is the same assumption as for a digital or (t; m; s)-net. The
first N values of the Van der Corput sequence, defined in (2.3) areb
Extensible Lattice Sequences 5
although in a different order. Therefore, the term
appears in the definition of the rank-1 lattice may be replaced by OE b (i), a term that
does not depend on N .
The s-dimensional generating vector h in (2.1) typically depends on N also. It
may be expressed in b-ary form as:
are digits. For k ? m the digits h jk do not affect the
definition of the rank-1 lattice set with points since they only contribute
integers to the product OE b (i)h. Therefore, each component of h may be written (in
principle) as an infinite string of digits:
This single "infinite" generating vector may serve for all possible values of m.
The preceding paragraphs provide the basis for defining an infinite rank-1 lattice
sequence. Altering the original definition in (2.1) leads to the following:
Definition 2.2. An infinite rank-1 lattice sequence in base b with generating
vector h of the form (2.6) and shift \Delta is defined as:
The first b m terms of the infinite rank-1 lattice sequence (2.7) are a rank-1 lattice.
Moreover, just as certain subsets of a (t; s)-sequence are (t; m; s)-nets, so subsets of
an infinite rank-1 lattice sequence are shifted rank-1 lattices.
Theorem 2.3. Suppose that P is the set consisting of the l +1 st run of b m terms
of the infinite lattice rank-1 sequence defined in (2.7):
Then, P is a rank-1 lattice with shift OE b (l)b \Gammam h + \Delta, that is,
Proof. The proof follows directly from the definition of the Van der Corput
sequence. For all note that
Substituting the right hand side into the definition of P completes the proof.
The definition of an infinite rank-1 lattice sequence may be extended to integration
lattices of arbitrary rank.
Definition 2.4. An infinite lattice sequence (of arbitrary rank) with bases
generating vectors h of the form (2.6) is defined as:
A practical complication for an integration lattice of rank greater than 1 is that there
are multiple indices, i k , each of which may or may not tend to infinity, and each at its
own rate. Because of this complication we will focus on rank-1 lattices in the sections
that follow. Theorem 2.3 also has a natural extension to infinite lattice sequences of
arbitrary rank, and its proof is similar.
6 F. J. Hickernell, H. S. Hong, P. L' '
Ecuyer and C. Lemieux.
3. Discrepancy. Unlike (t; s)-nets for which there exist explicit constructions
of the generating matrices C j , there are no such explicit constructions of generating
vectors h for rank-1 lattices for arbitrary s. Tables of generating vectors for lattices
that do exist (see [8, 14, 27]) are usually obtained by minimizing some measure of
non-uniformity, or discrepancy, of the lattice. This section describes several useful
discrepancy measures.
the quadrature error for a rule of the form (1.1) for an
arbitrary set P . Worst case error analysis of the quadrature error leads to a Koksma-
Hlawka-type inequality of the form [23]:
where D(P ) is the discrepancy or measure of nonuniformity of the point set defining
the quadrature rule, and V (f) is the variation or fluctuation of the integrand, f . The
precise definitions of the discrepancy and the variation depend on the particular space
of integrands.
In the traditional Koksma-Hlawka inequality (see [26] and [44, Theorem 2.11]),
the variation is the variation in the sense of Hardy and Krause, and the discrepancy
is the L1 -star discrepancy:
Here F unif is the uniform distribution on the unit cube, FP is the empirical distribution
function for the sample P , and j\Deltaj denotes the number of points in a set. The notation
denotes the L p -norm or the ' p -norm, depending on the context. Error bounds of
the form (3.1) involving the L p -star discrepancy have been derived by [57, 62]. Error
bounds involving generalizations of the star discrepancy appear in [21, 22, 23, 55].
When the integrands belong to a reproducing kernel Hilbert space, the error
bound (3.1) may be easily obtained [21, 23]. The discrepancy may be written in
terms of the reproducing kernel
Z
K(x; y) dx dy \GammaN
Z
For example, the L 2 -star discrepancy, whose formula was originally derived in [60], is a
special case of the above formula with K(x;
An advantage of
considering reproducing kernel Hilbert spaces of integrands is that the computational
complexity of the discrepancy is relatively small (at worst O(N 2 ) operations). By
contrast the L1 -star discrepancy requires O(N s ) operations to evaluate.
The discrepancy of type (3.3) can also be interpreted as an average-case quadrature
error [25, 41, 61]. Suppose that the integrand is a random function lying in the
sample space A, and suppose that the integrand has zero mean and covariance kernel,
K(x; y), that is,
Extensible Lattice Sequences 7
Then the root mean square quadrature error over A is the discrepancy as defined in
If P is a simple random sample, then the mean square discrepancy is [25]:
Z
This formula serves as a benchmark for other (presumably superior) low discrepancy
sets. Since the mean square discrepancy is O(N \Gamma1 ), the discrepancy itself is typically
O(N \Gamma1=2 ) for a simple random sample. The variance of a function, f , may be defined
as
Z
f dx
The mean value of the variance over the space of average-case integrands can be shown
to be [25]:
Z
Z
which is just the term in braces in (3.5).
It may seem odd at first that the discrepancy can serve both as an average-case
and worst-case quadrature error. The explanation is that the space of integrands, A,
in the average-case analysis is much larger than the space of integrands, W , in the
worst-case analysis. See [21, 25] and the references therein for the proofs of the above
results as well as further details.
There are some known asymptotic results for the discrepancies of (t; m; s)-nets.
The L1 -star discrepancy of any (t; m; s)-net is O(N \Gamma1 [log N Theorem 4.10].
Moreover, the typical (in the sense of an average taken over all possible nets) L 2 -
star discrepancy of (0; m; s)-nets is O(N \Gamma1 [log N the best possible for
any set. For discrepancies of the form (3.3) with sufficiently smooth kernels, typical
(0; m; s)-nets have O(N \Gamma3=2 [log N
Lattice rules, the topic of this article, are known to be particularly effective for integrating
periodic functions. Suppose that the integrand has an absolutely convergent
Fourier series with Fourier coefficients -
Z
Here k 0 x denotes the dot product of the s-dimensional wavenumber vector k with x.
The quadrature error for a particular integrand with an absolutely convergent
Fourier series is simply the sum of the quadrature errors of each term:
"N
8 F. J. Hickernell, H. S. Hong, P. L' '
Ecuyer and C. Lemieux.
where the term corresponding to does not enter because constants are integrated
exactly. One can multiply and divide by arbitrary weights w(k) inside this sum. Then
by applying H-older's inequality one has the following error bound of the form (3.1)
[23]:
"N
Here 1 f\Deltag denotes the indicator function. In order to insure that the discrepancy is
finite we assume that the weights increase sufficiently fast as k tends to infinity:
If P is the node set of an integration lattice, then it is known that trigonometric
polynomials are integrated exactly for all nonzero wavenumbers not in the dual lattice,
The dual lattice is the set of all k satisfying k 0 . Thus, for node
sets of lattices the definition of discrepancy above may be simplified to
Certain explicit choices of w(k) have appeared in the literature. For example, one
may choose
where the over-bar notation is defined as
are arbitrary positive weights, and ff is a measure of the assumed smoothness
of the integrand. If fi is the node set of a lattice, then
is a traditional figure of merit for
lattices [53]. Furthermore, for the discrepancy is DF;1 (P
ae(L) is the Zaremba figure of merit for lattice rules [44, Def. 5.31]. The more
general case of P not a lattice is considered in [20], and the case of unequal weights
discussed in [22].
Extensible Lattice Sequences 9
If the weight function w(k) takes the form (3.9) for positive integer ff, then for
the infinite sum defining the discrepancy in (3.8b) may be written as a
finite sum:
s
Y
for general sets P , where B 2ff denotes the Bernoulli polynomial of degree 2ff [1, Chap.
When P is the node set of an integration lattice, the double sum can be simplified
to a single sum:
s
Y
Another choice for w(k) is a weighted ' r -norm of the vector k to some power:
again for arbitrary positive weights these weights are unity, P is
the node set of a lattice, and
kkk \Gammaff
ae
min
kkk r
oe \Gammaff
For discrepancy is equivalent to the spectral test, commonly employed
to measure the quality of linear congruential pseudo-random number generators [30,
33]. The spectral test has been used to select lattices for quasi-Monte Carlo quadrature
in [7, 34, 35, 36]. The case which one might call an ' 1 -spectral test is also
interesting. We will return to these two cases in the next section.
4. Good Generating Vectors for Lattice Sequences. As mentioned at the
beginning of the previous section, finding good generating vectors for lattices typically
requires optimizing some discrepancy measure. In this subsection we propose some
loss functions and optimization algorithms for choosing good generating vectors for
extensible rank-1 lattice sequences.
In principle one would like to have an s \Theta 1 array of digits h jk , according to (2.6).
However, in practice it is only necessary to have an s max \Theta mmax array of digits h jk ,
is the maximum number of points and s max is the maximum dimension
to be considered. In finance calculations, for example, the necessity of timely
forecasts may constrain one to a budget of
private communication).
For simplicity we consider generating vectors h that are of the form originally
proposed by Korobov, that is,
This means that only the digits need to be chosen, for which there
are b choices. The generating vector is tested for dimensions up to s max ,
but in fact it can be extended to any dimension.
F. J. Hickernell, H. S. Hong, P. L' '
Ecuyer and C. Lemieux.
The number j defining the generating vector is chosen by minimizing a loss func-
tion, G, of the form
~
Here the function ~
is related to one of the measures of discrepancy introduced
in Section 3, and the maximum over some range of values of m and s insures that the
resulting generating vector is good for a range of numbers of points and dimensions.
However, since the discrepancy itself depends significantly on the number of points
and the dimension, it must be appropriately scaled to arrive at the function ~
G. The
details of this scaling are given below.
4.1. Generating Vectors Based on Minimizing P ff . The discrepancy defined
in (3.11), which is a generalization of the P ff figure of merit for lattice rules, has
the advantage of requiring only O(sN) operations to evaluate for lattices. To remove
some of the dimension dependence of this discrepancy it is divided by the square root
of right hand side of (3.6). The root mean square of this scaled discrepancy for a random
sample is then N \Gamma1=2 , independent of s. The formula for the scaled discrepancy
of the node set of a lattice with
s
Y
The specific choice of the value of fi j here is not crucial, but seems to give good results.
For one-dimensional lattices, i.e. evenly spaced points on the interval [0; 1), this
discrepancy is N \Gamma1 , and one would expect that as the dimension increases this scaled
discrepancy would tend to (or at least do no worse than) N \Gamma1=2 . Therefore, to remove
this remaining dimension dependence the above scaled discrepancy is divided by the
For a fixed s, D asy (m,s) is asymptotically O(b \Gammam m (s\Gamma1)=2
as N tends to infinity. This is the asymptotic order for (0; m; s)-nets [25], and what
we hope to achieve for lattice sequences. Furthermore, N
and s. In summary, the resulting loss function is
~
\Theta
s
Y

The optimal values of j found by minimizing G 1 (j) for different
ranges of m and s are given in Table 4.1. The algorithm for optimizing G 1 (j) may
be described as an intelligent exhaustive search. One need not compute G 1 (j) for all
possible values of j. Suppose at any stage of the optimization j   is the best known
value of j, and one finds that ~
j. Since ~
depends only on the first digits of j, one can immediately eliminate from
consideration all j that have the same first m \Gamma 1 digits as ~
j. This same search
strategy is also used for the other loss functions described below.
Extensible Lattice Sequences 11
4.2. Generating Vectors Based on the Spectral Test. The use of the spectral
test to analyze the lattice structure of linear congruential generators is described
in [30] and tables of good integration lattices are given in [34]. The difference here is
that nearly all smaller lattices imbedded in the largest lattice considered must have
low discrepancy. In [34], only the full lattice was examined.
The length of the shortest non-zero vector in the dual lattice L ? is
which is related to the discrepancy (3.12) with 2. This length has the absolute
upper bound
d
where the constants fl s and ae s depend only on s (see [34] and the references therein).
The bound for s - 8 is the least upper bound for a general s-dimensional lattice with
real-valued coordinates, and with b \Gammam points per unit of volume. The bound for s ? 8
is not the least upper bound, but it is still reasonably tight, as our numerical results
will show. We define the normalized ' 2 -spectral test discrepancy as
~
which is larger than 1 and is the inverse of the quantity S t defined in [34]. (The
different notation here is to be consistent with the rest of this article). The loss
function to be minimized is of the form (4.2) with ~
We note that 1=d 2 (j; m; s) can be interpreted as the (Euclidean) distance between
the successive hyperplanes that contain all the points of the primal lattice L, for the
family of hyperplanes for which this distance is the largest. The problem of computing
a shortest vector in (4.4) can be formulated as a quadratic optimization problem
with s integer decision variables, because k can be written as a linear combination
of the s vectors of a basis of the dual lattice, with integer coefficients. The decision
variables are these coefficients. (See [30] for details.) We solved this problem by using
the branch-and-bound algorithm of Fincke and Pohst [10], with a few heuristic modifications
to improve the speed. The worst-case time complexity of this algorithm is
exponential in d 2 (j; m; s), and polynomial in s for d 2 (j; m; s) fixed [10]. In practice, it
(typically) works nicely even when d 2 (j; m; s) is large. For example, one can compute
2, an arbitrary j and in less than 1 second on a
Pentium-II computer.
4.3. Generating Vectors Based on the ' 1 -spectral Test. With the ' 1 norm,
the length of the shortest non-zero vector in L ? is
which is related to the discrepancy (3.12) with 1. One has the upper bound
which was established by Marsaglia [37] by applying the general convex body theorem
of Minkowski. This suggests the normalized ' 1 -spectral test quantity:
~
12 F. J. Hickernell, H. S. Hong, P. L' '
Ecuyer and C. Lemieux.
Here, we want to minimize the loss function (4.2) with ~
G 3 .
One can interpret d 1 (j; m; s) (or d 1 (j; m; s) \Gamma 1 in certain cases; see [30]) as
the minimal number of hyperplanes that cover all the points of P . We computed
via the algorithm of Dieter [6], which works fine for s up to about 10
(independently of m), but becomes very slow for larger s (the time is exponential in
s).
The test quantity in (4.8) has an interpretation similar that the quality
parameter t for (t; m; s)-nets. Define
s
Since ~
(s!)]=s.
Thus, the rank-1 lattice defined by j integrates exactly all trigonometric polynomials
of wavenumber k when
log(jk
s
If one considers T (j; m; s) as the quality parameter of the lattice, then this condition
is similar to that in (2.5). There, t determines the resolution at which piecewise
constant functions that are exactly integrated by a net. Here, T (j; m; s) determines
the resolution at which trigonometric polynomials are integrated exactly by a lattice.
The discrepancy for the node set of this lattice as defined in (3.12) with
If one can construct an infinite sequence of digits, j, for which
then the above discrepancy decays like N \Gammaff=s . Again, the parameter ff indicates the
assumed smoothness of the integrands.
4.4. Tables of Coefficients. We made computer searches to find the best j's,
based on minimizing the worst-case loss function
~
selected values of m . The bounds m
define a range of number of points in the lattice, and maximal number of dimensions,
that we are interested in. Selecting different parameters j for different ranges of
values of m and s is a convenient compromise between the extreme cases of choosing
a different j for each pair (m; s), and choosing the same j for all pairs (m; s). In
practice one typically has a general idea of how many points one is going to take. By
selecting an j specialized for that range, one can obtain a lattice with a better figure
of merit for this particular range.

Table

4.1 gives the optimal j's and the corresponding figures of merit (4.9) for
Because of computational efficiency constraints,
for the searches, we limited ourselves to m 1 - 20 for
Then, for the best j that we found, we verified the performance when m 0 was reduced
Extensible Lattice Sequences 13
or when m 1 or s 1 was increased, and retained the smallest m 0 and largest m 1 and s 1
for which G i was unchanged. The table also gives the value of G i;min (j; m
defined by replacing max by min in (4.9). This best-case figure tells us the range of
values taken by ~
over the given region. We also made exhaustive searches
with all the entries where in the table, and obtained the same
values of j and G cases.

Table
Values of j defining generating vectors of the form (4.1) for rank-1 lattice sequences with base
minimizing (4.9).
5. Estimating Quadrature Error. The advantage of an extensible lattice sequence
is that N need not be specified in advance. Therefore, in practice one would
estimate the quadrature error for an initial lattice, and continue to increase the lattice
size until the quadrature error meets the desired tolerance. Although the discrepancy,
D(P ), is a good measure for the quality of a set P , it cannot be used directly for error
estimation. The worst case error bounds (3.1) are often quite conservative, and there
is no easy way to estimate the variation of the integrand. The average case error
given, (3.4), is sensitive to how one defines the kernel - multiplying the kernel by a
factor of c 2 changes the average case error by a factor of c.
Quadrature error estimates for lattice rules have been investigated by [5, 28,
53]. Two different kinds of quadrature rules and error estimates have been proposed.
Both involve estimating the error for terms of Q(f
, where the P l are node sets of lattices and
The method proposed in [53] takes P to be composed of 2 s shifted copies of a
lattice. Rather than considering copy rules, we take P l to be the node sets of
shifted lattices of size b m1 imbedded in the extensible lattice described in Section 2.4:
14 F. J. Hickernell, H. S. Hong, P. L' '
Ecuyer and C. Lemieux.
Another quadrature rule and error estimate takes the P l to be independent random
shifts of a lattice. This can be done by taking:
where the \Delta l are independent, uniformly distributed random vectors in [0; 1] s .
Note that for both (5.1) and (5.2) the set P can be extended in size as necessary
by increasing m 1 . The theory behind the error estimates for cases (5.1) and (5.2) are
given in the following theorem.
Theorem 5.1. Suppose that a quadrature rule Q(\Delta; P ) based on some arbitrary
P , as given in (1.1), is the average value of the quadrature rules based on the sets
First, consider the case where the integrands are random functions from a sample
space A as described in the paragraph preceding (3.4). Then it follows that
where D is the discrepancy based on the covariance kernel for A.
Secondly, consider the case of a fixed integrand, f , but where the P l are random
shifts of a set P 0 , that is,
for independent, uniformly distributed
Proof. Assuming that the quadrature rules satisfy (5.3), it follows that the mean
square deviation of the Q(f ; P l ) from Q(f ; P ) may be written asM
If the integrand is a random function from a sample space A with covariance kernel
K, then average case error analysis in (3.4) plus (5.7) leads to the following equations:
Extensible Lattice Sequences 15
The equations above may be rearranged to give (5.4).
The quadrature error for rule Q(f ; P ) may also be written as:
Substituting the sum of the by equation (5.7) gives:
If the P l are random shifts as in (5.5), then the expected value of the term [I(f) \Gamma
vanishes for all k 6= l and taking the expected value of (5.8)
yields (5.6).
Some remarks are in order to explain the assumptions and conclusions of the
above theorem. These are given below.
Assumption (5.3), and thus conclusion (5.4), holds for both the cases (5.1) and
(5.2) above. In fact, this part of the theorem holds for any imbedded rule or any rule
where the P l all contain the same number of points, and their union is P or multiple
copies of P . For example, (5.4) would apply to the case where P is a (t; m; s)-net
made up of a union of subnets P l .
Assumption (5.5), and therefore conclusion (5.6) holds for (5.2). There is a difficulty
if one tries to derive a result like (5.6) for an imbedded rule of the form (5.1),
where \Delta is a random shift. The argument leading to (5.6) assumes that the points in
different P l are uncorrelated, which is not true for (5.1). However, if the extensible
lattice is a good one, it is expected that the terms
in (5.8) are on average negative. Under this assumption one may then conclude that
the right hand side of (5.6) is a conservative (too large) upper bound on the expected
square quadrature error.
The factor in (5.4) above involving the discrepancies of P and the P l does not
depend strongly on the particular choice of discrepancy, but only on the asymptotic
rate of decay. If, for example, D(P ) - CN \Gammaff for some unknown C, but known ff,
where N is the number of points in P , then
Although conclusions (5.6) and (5.9) are derived under different assumptions,
they both suggest error estimates of the form
F. J. Hickernell, H. S. Hong, P. L' '
Ecuyer and C. Lemieux.
indicates the rate of decay of the discrepancy for (5.9).
The factor c ? 1 depends on how conservative one wishes to be. The Chebyshev
inequality implies that the above inequality will hold "at least" 100(1 \Gamma c \Gamma2 )% of the
time. Error estimate (5.10) leads to the stopping criteria:
where ffl is the absolute error tolerance. Note that if the stopping criteria is not met,
one would normally increase the size of the P l by increasing m 1 , rather than increasing
the number of the P l by increasing M .
For some high dimensional problems the discrepancy of a lattice (or other low
discrepancy set) may decay as slowly as the Monte Carlo rate of O(N \Gamma1=2
N (see [18, 41]). Therefore, even when using (5.9), it may be advisable to make a
conservative choice of ff = 1=2. This choice makes the approach of error estimate
(5.9), based on random integrands, equivalent to that of (5.6), based on randomly
shifted P l .
Sloan and Joe [53, Section 10.3] suggest an error estimate of the form
s
s
where P is formed from 2 s copies of a rank-1 lattice, and each Q(f ; P l ) is an imbedded
rule based on half of the points in P . Although this case does not exactly fit Theorem
5.1, the arguments in the proof can be modified to obtain a result similar to to (5.4):
s
s
s
This would suggest that the error estimation formula of Sloan and Joe is reasonable
when
2D(P ) on average. The disadvantage of 2 s -copy rules is that they
require at least 2 s points, which may be unmanageable for large s.
To summarize, both imbedded lattice rules, (5.1), and independent random shifts
of lattices, (5.2), have similar error estimates, (5.10), and stopping criteria, (5.11).
The advantage of the independent random shifts approach is that, the theory holds
for any integrand, not the average over a space of integrands. One advantage of the
imbedded rules approach is that one need only generate a single extensible lattice.
Furthermore, the set P for the imbedded rule is the node set of a lattice (if M is a
power of b), which is never the case for the independent random shifts approach. Thus,
the accuracy of the imbedded rule approach is likely to be better. In the examples in
the next section, the imbedded lattice rules based on (5.1) are used.
6. Examples of Multidimensional Quadrature. Two example problems are
chosen to demonstrate the performance of the new rank-1 lattice sequences proposed in
Section 2.4. The first example is the computation of multivariate normal probabilities
and the second is the evaluation of a multidimensional integral arising in physics
problems.
Extensible Lattice Sequences 17
6.1. Multivariate Normal Probabilities. Consider the following multivariate
normal
Z b1
Z bs
as
where, a and b are known s-dimensional vectors, and \Sigma is a given s \Theta s positive definite
covariance matrix. Some a j and/or b j may be infinite.
Unfortunately, the original form is not well-suited for numerical quadrature.
Therefore, Alan Genz [12] proposed a transformation of variables that results in
an integral over an s \Gamma 1-dimensional unit cube. See [12, 13] for the details of the
transformation, and see [13, 15] for comparisons of different methods for calculating
multivariate normal probabilities.
The particular test problem is one considered by [13, 24] and may be described
as follows:
i.i.d. uniformly on [0;
generated randomly according to [13, 38]: (6.1d)
Numerical comparisons were made using three types of algorithms:
i. the adaptive algorithm DCHURE [2],
ii. an older Korobov rank-1 lattice rule with a different generating vector for
each N and s - this algorithm is a part of NAG and is used in [5, 12], and
iii. the new rank-1 lattice sequences proposed in Section 2.4 with generating
vectors given in Table 4.1
For the second and third algorithms we applied the periodizing transformation x 0
1j to the integrand over the unit cube. This appears to increase the accuracy
of the lattice rule methods. The computations were carried out in FORTRAN on a
Unix work station in double precision. The absolute error tolerance was chosen to be
, and this was compared with the actual error E. Since the true value of the
integral is unknown for this test problem, the value given by the Korobov algorithm
with a tolerance of used as the "exact" value for computing the error.
For the new rank-1 lattice sequences the stopping criterion (5.11) was used with M
between 4 and 7 and
For each dimension 50 random test problems were generated and solved by the
various quadrature methods. The scaled absolute errors E=ffl and the computation
times in seconds are given in the box and whisker plots of Fig. 6.1. The boxes contain
the middle half of the values and the whiskers give the range of most values except
the outliers (denoted by  ).
Ideally, the scaled error should nearly always be less than one, otherwise the
error estimate is not conservative enough. On the other hand if the scaled error is too
small, then the error estimate is too conservative. Fig. 6.1 shows that the adaptive
rule performs well for smaller dimensions, but underestimates the error and is quite
slow in higher dimensions. The lattice rules do well even in higher dimensions, and
the new rank-1 lattice sequences appear to be faster than the older Korobov-type
rule. This is likely due to the fact that the lattice sequences proposed here can re-use
the old points when N must be increased.
F. J. Hickernell, H. S. Hong, P. L' '
Ecuyer and C. Lemieux.
6.2. A Multidimensional Integral from Physics. Keister [29] considered
the following multidimensional integral that has applications in physics:
Z
R s
Z
cos@
A dy; (6.2)
where \Phi denotes the standard Gaussian distribution function. Keister gave an exact
formula for the answer and compared the quadrature methods of McNamee and
Stenger [39] and Genz and Patterson [11, 51] for evaluating this integral. Later, Papageorgiou
and Traub [49] applied the generalized Faure sequence from FINDER to
this problem.
The results of numerical experiments for the above integral for dimension 25 are
shown in Figure 6.2. The exact value of the integral is reported in [49]. To be
consistent with the numerical results reported in [29, 49], we did not perform error
estimation, but just computed the actual error for each kind of numerical method as
a function of N , the number of points. Because \Phi there is a technical
difficulty with using an unshifted lattice rule, so when performing the numerical experiments
the lattice sequences were given random shifts (modulo 1). Box and whisker
plots show how well the new rank-1 lattice sequences perform for 50 random shifts.
According to Figure 6.2 the generalized Faure sequence (in base 29) and the
lattice sequence perform much better than the other two rules. In some cases the
lattice sequences perform better than the generalized Faure sequence.
7. Conclusion. Lattice rules are simpler to code than digital nets. Given the
construction in Section 2.4, it is now possible to have extensible lattice sequences in
the same way that one has (t; s)-sequences. Good generating vectors for these lattice
sequences may be found by using the spectral test or minimizing other discrepancy
measures, as shown in Section 4. The performance of these lattice rules is in many
cases comparable to other multidimensional quadrature rules and in some cases superior

Acknowledgments

. Thanks to Alan Genz for making available his software for
computing multivariate normal distributions. Also thanks to Joe Traub and Anargy-
ros Papageorgiou for making available the FINDER software.



--R

Handbook of Mathematical Functions with Formulas
An adaptive algorithm for the approximate calculation of multiple integrals

Valuation of mortgage backed securities using Brownian bridges to reduce effective dimension
Randomization of number theoretic methods for multiple integration
How to calculate shortest vectors in a lattice


Discr'epance de suites associ'ees 'a un syst'eme de num'eration (en dimension s)
Improved methods for calculating vectors of short length in a lattice
A Lagrange extrpolation algorithm for sequences of approximations to multiple inte- grals


Parameters for integrating periodic functions of several variables
Simulation of multivariate normal rectangle probabilities and their derivatives theoretical and computational results
On the assessment of random and quasi-random point sets
Random and Quasi-Random Point Sets
A comparison of random and quasirandom points for multidimensional quadrature





Computing multivariate normal probabilities using rank-1 lattice sequences

Funktionen von beschr-ankter Variation in der Theorie der Gleichverteilung
Applications of Number Theory to Numerical Analysis
Randomization of lattice rules for numerical multiple integration
Multidimensional quadrature algorithms
The Art of Computer Programming
The approximate computation of multiple integrals
On the distribution of digital sequences




Random numbers fall mainly in the planes

Construction of fully symmetric numerical integration formu- las
Generating quas-random paths for stochastic processes


Points and sequences with small discrepancy

Shiue eds

Quasirandom points and global function fields


Faster valuation of financial derivatives
The optimum addition of points to quadrature formulae
Lattice methods for multiple integration
Lattice Methods for Multiple Integration
Lattice methods for multiple integration: Theory
An intractability result for multiple integration
The distribution of points in a cube and the approximate evaluation of integrals
"Nauka"


Computational investigations of low discrepancy point sets
Average case complexity of multivariate integration
Some applications of multidimensional integration by parts
--TR

--CTR
Fred J. Hickernell , Harald Niederreiter, The existence of good extensible rank-1 lattices, Journal of Complexity, v.19 n.3, p.286-300, June
Hee Sun Hong , Fred J. Hickernell , Gang Wei, The distribution of the discrepancy of scrambled digital (t, m, s)-nets, Mathematics and Computers in Simulation, v.62 n.3-6, p.335-345, 3 March
Pierre L'Ecuyer, Quasi-monte carlo methods in practice: quasi-monte carlo methods for simulation, Proceedings of the 35th conference on Winter simulation: driving innovation, December 07-10, 2003, New Orleans, Louisiana
Fred J. Hickernell, My dream quadrature rule, Journal of Complexity, v.19 n.3, p.420-427, June
Hee Sun Hong , Fred J. Hickernell, Algorithm 823: Implementing scrambled digital sequences, ACM Transactions on Mathematical Software (TOMS), v.29 n.2, p.95-109, June
good lattice point sets;multidimensional;discrepancy;spectral test