- Email: [email protected]
Discrepancy behaviour in the non-asymptotic regime
Applied Numerical Mathematics 50 (2004) 227–238 www.elsevier.com/locate/apnum
Discrepancy behaviour in the non-asymptotic regime Ch. Schlier Physikalisches Institut, Universität Freiburg, 79104 Freiburg, Germany Available online 29 January 2004
Abstract We have computed true star-discrepancies D ∗ (N) for Halton and Niederreiter point sequences in dimensions s = 2, 3, 4, and 6 with N up to 250 000, 10 000, 2000 and 300, respectively. For comparison, we also calculated some L2 discrepancies T ∗ . The mean behaviour of D ∗ (N) can well be approximated by power laws with an exponent between about −0.7 and −0.9, which slowly decreases with N. This behaviour is far from the generally assumed asymptotic one ∼ C(s)(log N)s /N. The factor between the true and the asymptotic behaviour increases strongly with s and reaches many orders of magnitude for large s. The ratios of D ∗ (N) for different low discrepancy sequences are not proportional to the presumed asymptotic pre-factors C(s). Especially for the range of bases investigated, Niederreiter sequences have about the same D ∗ as Halton ones in a range of N, which we conjecture to be at least of the order of 1010 . 2003 IMACS. Published by Elsevier B.V. All rights reserved. Keywords: Quasi-Monte Carlo methods; Low-discrepancy sequences; Numerical integration
1. Introduction ∗ ) defined by [12] The star-discrepancy of a s-dimensional point sequence (sometimes also called D∞ N 1 ∗ χJ (xn ) − V (J ), (1) D (N) = sup s N J ⊆[0,1) n=1
is a measure for the uniformity of distribution of the points of the sequence xn , n = 1, 2, . . . in the unit hypercube [0, 1)s . Here, J denotes any subinterval of the s-dimensional unit cube which has faces parallel to the coordinate axes and one corner at the origin, and V (J ) its volume. The characteristic function χJ (n) of the points is 1 if the point is in J , else 0. I.e., D ∗ measures the largest difference E-mail address: [email protected] (Ch. Schlier). 0168-9274/$30.00 2003 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.apnum.2003.12.004
228
Ch. Schlier / Applied Numerical Mathematics 50 (2004) 227–238
between the relative number of points in J and the volume of J . Customarily, a sequence for which an estimate of its discrepancy of the type ∗ (N) = C(s) · (log N)s /N + O (log N)s−1 /N , (2) D ∗ (N) Das is proved is called a low-discrepancy sequence (LDS) [12,8], and it is generally conjectured [12] that the order O(log N)s /N is the best possible for the discrepancy of an infinite sequence. In the sequel, for any low-discrepancy sequence we use the smallest constant C(s), for which such an inequality is proved. The constants C(s) for many LDSs can be found in [12,8], and play an ubiquitous role in the discussion of quasi-Monte Carlo methods. The main reason for this role is the Koksma–Hlawka inequality for the worst case error ε max of quasi-Monte Carlo integration εmax V (f ) · D ∗ (N).
(3)
Here, V (f ) is the variation of the integrand [10,12], a measure of its smoothness, which, however, is unknown for most functions. The pre-factors C(s) in Eq. (2) are generally thought to determine the practical quality of a sequence. However, Eq. (2) describes an asymptotic limit, valid for N → ∞, and the next term of order (log N)s−1 /N is not known. Astonishingly, while everybody mentions Eqs. (2) and (3), nobody knows where, in practice, the asymptotic regime begins. Is it at N = 105 or 1025 or 10100 ? And, strangely enough, there are few, if any, numerical computations of D ∗ (N) for smaller N , i.e., those which are usually used in integration problems. Recently we have performed extensive tests of quasi-Monte Carlo integration with integrands whose variation we could compute [17]. These tests confirmed not only the well-known fact that practical integration errors are many orders of magnitude smaller than Eq. (3) would suggest, but showed also that the relative magnitudes of the errors for different point sequences do not at all behave like the prefactors C(s) predict, at least in a range up to 1010 trials. In Fig. 1 we show an example. The reduced integration error for two test functions from [1], i.e., the error divided by the variation of the integrand, is plotted, which can be directly compared with D ∗ . This led us to ask whether this situation might be a consequence of the non-asymptotic regime of D ∗ , or if everything is due to the difference between practical and worst case errors. The answer is obvious from the figure and discussed below. In the literature one finds formulas for the computation of the root mean square discrepancy T ∗ (also called D2∗ ), defined by replacing the absolute norm in Eq. (1) by the quadratic norm [12]. These algorithms have an expense ∼ s 2 N 2 , so we could compute T ∗ to rather large N , e.g., for dimension 6 up to N = 106 . For all sequences and dimensions which we investigated the functions T ∗ (N) observe an approximate power law in the non-asymptotic region. Unfortunately, there is no known way to draw useful conclusions from T ∗ to D ∗ , which is needed in Eq. (3). So, even if this behaviour of T ∗ could be a hint that D ∗ might behave similarly, this is not proved theoretically. Since obviously D ∗ (1) < 1, and intuitively one would expect D ∗ to decrease with N , it is trivial that D ∗ (N) for small N must behave differently from the asymptotic behaviour. But, again, where does the asymptotic regime begin? A handwaiving argument says that it must be much farther out than the maximum of the asymptotic D ∗ (N) in Eq. (2), which is at es . But how much? A more restrictive argument [18] maintains that D ∗ (N) should be at least < D ∗ (2) before the asymptotic region can begin. This, then, would mean that the limiting N is about 3000 for s = 3, but already 1.6 × 109 for s = 6, increasing to computationally unfeasible values for larger s. So we started a project to compute D ∗ for several Halton, Niederreiter, Niederreiter–Xing and Sobol sequences, which we describe in the next sections. The maximum N for which we could compute exact
Ch. Schlier / Applied Numerical Mathematics 50 (2004) 227–238
229
Fig. 1. Asymptotic discrepancy (thick smooth curve) and true discrepancy (dots with thin fitted line) of the Halton sequence in 6 dimensions (primes 2 to 13) together with reduced errors of integration tests (dashed lines with dashed fits) performed with this sequence. Plotted is the root mean square error of repeated integrations of test functions 1 (lower) and 4 (upper) from [1] divided by the integrand’s variation (data from [17]). The heavy dots show the square root of the variances of the test functions. These dots are the origins of the expectation value of the reduced errors of pseudo-Monte Carlo integration. Note that the asymptotic discrepancy overestimates the true one appreciably, but the practical reduced errors are still much lower than their Koksma–Hlawka limit.
values of D ∗ was limited by our computer resources to about 250 000, 10 000, 2500, and 300 for s = 2, 3, 4, and 6, respectively. (We had to progress in steps > 1 for the larger values of N .) Niederreiter and Halton sequences were computed with different bases respectively base sets. For these s, and some larger ones we have also computed approximate D ∗ . The convergence of our simple approximation to D ∗ for large s and N is, however, very bad. In view of the many orders of magnitude difference from the asymptotic values they are, however, quite useful even with their errors. Our main result can be summed up in two statements: (1) D ∗ (N) behaves approximately according to a power law with very slowly decreasing negative power for those N which determine practical integration results, and (2) the relative magnitude of D ∗ (N) for different low-discrepancy sequences does not behave as the pre-factors C(s) of the asymptotic behaviour predict. Especially, at practical values of N and for the range of bases investigated, D ∗ (N) does not justify a preference of Niederreiter sequences over Halton ones. As a consequence for the practitioner, we propose that he should feel free to choose the simpler and faster available Halton sequences and also other LDSs considered “non-optimal” in the literature (e.g., Niederreiter sequences to base 2 even if this is not the “optimal” base) without fear of large penalties in performance. As a consequence, repeated integrations with different LDSs can be recommended as a method to increase the overall confidence in the results.
230
Ch. Schlier / Applied Numerical Mathematics 50 (2004) 227–238
2. Calculation methods We have used one approximate and two exact algorithms to compute D ∗ . Our first idea was to replace the supremum over all hyper-rectangles in Eq. (1) by the maximum over a large number, M, of randomly ∗ (N, M) chosen ones. The computational effort is s · N · M. Unfortunately, the convergence of this Dappr with increasing M is not good. Some convergence plots can be found in [17]. Actually, one can argue that the number of hyper-rectangles needed here must be larger than that which one must visit in the exact programs described below. On the other hand, this algorithm is very simple, and can, therefore, easily be checked for correctness. So we could use it with small N (say < 20) and large M (up to 1010 ) to check, at least punctually, the correctness of the implementation of the two exact algorithms described below, whose structures are much more complex. Also, if one accepts an error of one order of magnitude, which is not much compared to the many orders difference from the asymptotic behaviour, one can use even these approximate values for a crude estimate of the true ones at larger N , cf. Fig. 2. This figure includes ∗ (N) and T ∗ (N), and shows that they are approximately parallel. This was found also in the plots of Dappr other cases, and means that these two types of discrepancy differ only by a slowly changing factor. For an exact computation of D ∗ , one must consider all hyper-rectangles at whose corners the number of included points may possibly change. For Halton and Niederreiter (and other) low-discrepancy sequences these lie on a grid, whose size can be computed. In the simple Halton case it is the product of (in Mathematica notation) pk Ceiling[Log[pk , N]] , where the pk are the primes used in the construction
Fig. 2. True discrepancies in 3 dimensions for sequences (see text) G2 (!) and H1 (e). Thin straight lines: linear fits of log10 (D ∗ ) vs. log10 (N) of G2 points in the 2nd, respectively, 4th decade of N . Thick lines: approximate D ∗ for G2 (long dash) and H1 (short dash). Thin lines: T ∗ for G2 (long dash) and H1 (short dash).
Ch. Schlier / Applied Numerical Mathematics 50 (2004) 227–238
231
of the sequence. This is of the order of s! · N s . For each such possible hyper-rectangle one must check which points of the sequence are contained in it, which leads to an effort of s! · N s+1 for this algorithm. We have used it in our older calculations [17], but could only handle N = 1000 and 100 for s = 2 and 3, respectively. A better algorithm has been published for s = 2 and 3 in a difficult to find paper [4]. It can easily be extended to other dimensions, and we have done so for s = 3 to 12. The idea is to investigate only those hyper-rectangles at whose corners the number of included sequence points does change. The algorithm can best be described by an excerpt from our Mathematica implementation, here in 3 dimensions, the extension to other dimensions being obvious. Note that Mathematica does not allow 0 as index, so all indexes are increased by one compared to mathematical notation. Note also that Union[] includes sorting: ptab =Union[{{0,0,0}}, Take[pointlist,n]]; xtab = Union[Map[ex31, ptab], {1}]; (∗ ex31 extracts the first component of three ∗ ) mx1 = mx2 = 0; For[l = 0, l n, l++, yta = Map[ex323, Take[ptab, l+1]]//Sort; (∗ ex323 extracts components 2 and 3 of three ∗ ) ytab = Union[Map[ex21, yta],{1}]; (∗ ex21 extracts the 1st of two components ∗ ) mx10 = mx20 = 0; For[k = 0, k l, k++, zta = Map[ex22, Take[yta, k+1]]//Sort; (∗ ex22 extracts the 2nd of two components ∗ ) ztab = Union[zta,{1}]; mx100 = mx200 = 0; For[ j = 0, j k, j++, taba = Abs[j/n - xtab[[l+1]]∗ytab[[k+1]]∗ztab[[j+1]] ]; tabb = Abs[xtab[[l+2]]∗ytab[[k+2]] ∗ ztab[[j+2]]- j/n ]; mx100 = Max[Abs[taba],mx100]; mx200 = Max[Abs[tabb],mx200]; ]; mx10 = Max[mx10, mx100]; mx20 = Max[mx20, mx200]; ]; mx1= Max[mx1,mx10]; mx2 = Max[mx2,mx20]; ]; mx = Max[mx1,mx2]; The input is a point list of length > n of the points of the—in this case 3-dimensional—sequence, and the output is mx = D ∗ (n) as rational number. The sequel of projections on hyper-rectangles of ever lower dimension restricts the grid points which must be visited to those where the number of included points of the sequence actually changes, and the kind of indexing makes this number known. This leads to an asymptotic number of N s /s! hyper-rectangles, which must be investigated. The limit on single values of N for our equipment (four 450 MHz Xeons, computations up to one week accepted) was about 35 000, 1500, 350, and 90 for s = 2, 3, 4, and 6, respectively. An implementation of this algorithm in Fortran 90 (now with double precision “real” numbers instead of integers) speeds it up by about three orders of magnitude. An example program for s = 5 can be found
232
Ch. Schlier / Applied Numerical Mathematics 50 (2004) 227–238
on the author’s web page http://www.phya8.uni-freiburg.de/abt/papers/papers.html. With this program we could compute single values of D ∗ (N) up to N = 250 000, 10 000, 2500 and 300 for dimensions s = 2, 3, 4, and 6, respectively. Another factor of 10 to 100 might be gained with one or several faster processors. But note that 1001/6 is only 2.15! So, the calculation of D ∗ in higher dimensions for practically interesting N is still unfeasible today.
3. Results We have computed exact discrepancies in 2, 3, 4, and 6 dimensions for the following sequences: Halton points [5] produced from consecutive primes starting from the first (2), the third (5), or the 7th (17), abbreviated as H1, H3, and H7, points from Niederreiter’s (t, s)-sequences [11], Niederreiter points for short, with bases 2, 3, 5, and 7, abbreviated G2, G3, G5, and G7. In addition we computed some sequences with much larger bases, Halton points H20 and H30 (same notation as above) and Niederreiter points G29 and G47, and finally some Niederreiter–Xing and Sobol sequences. Halton points were produced by our simple algorithm [2], Niederreiter ones by the TOMS algorithm 728 [1], Niederreiter– Xing points using [3,15], and Sobol points by the subroutine in Ref. [16]. As we discussed above, the maximum feasible N depends on the dimension. To avoid clutter, we plot only the Halton results for H1, and H7, those for H3 lying generally between them. For Niederreiter sequences we restrict our plots to only G2 and the so called “optimal” sequence, i.e., that with the smallest C(s) for the given dimension. These are G2, G3, G3, and G7 for the four dimensions involved. For s = 2 we included also G7. Note, however, that on today’s binary computers any Niederreiter sequence other than G2 is much more slowly produced. We do not plot Niederreiter–Xing or Sobol discrepancies, since they are on the average not much different from G2. The computed values are depicted in Figs. 3–6 as points. For small N we stepped N by 1, while for higher N larger steps had to be used. D ∗ (N) is not a monotonously decreasing function but fluctuates about some average trend. Local maxima and minima follow each other in short runs, and are, therefore, not too interesting for applications. These fluctuations are smallest for the sequences with small bases, as one would perhaps expect from pictures of two-dimensional projections of the sequence points, which do not show large “holes” in this case. Their relative amplitude seems to be almost constant, but we could not really ascertain this. Support for this statement comes, however, from a short test (in dimension 3, with N = 1 . . 1500) performed with a series of ten sequences of G2 type scrambled according to the prescription of [6], which is a simplified version of Owen’s scrambling [14]. As expected, the average discrepancy of this set fluctuates much less than that of a single one. In contrast, the relative standard deviation of the Di∗ (N) (i.e., the standard deviation divided by the mean) decreases very slowly with N . In the following we discuss several observations, which the reader can verify from the figures. Extrapolation from the computed data to higher values of N lies, however, at hand. It seems justified for several reasons: (1) From the construction of a low-discrepancy sequence one certainly does not expect a sudden break in its mean behaviour. ∗ and for the (2) The same general behaviour as for the true D ∗ is also observed for the approximate Dappr ∗ L2 -discrepancy T , in which cases we could compute data for much larger N .
Ch. Schlier / Applied Numerical Mathematics 50 (2004) 227–238
233
Fig. 3. True discrepancies in 2 dimensions for sequences (see text) G2 (!), G7 (1), H1 (e), and H7 (P). Points of H7 and later also G7 lie distinctly above the others. The insert shows the same points plus the asymptotic discrepancies for the sequences (from top to bottom: H7, G7, H1, G2).
(3) This behaviour is also observed for the error of test integrations of functions with known variation, which we computed up to N = 1010 trials [17] (cf. Fig. 1). Indeed, it was these test integrations which led use to the investigation of D ∗ . We also think that some extrapolation to larger dimensions is allowed, which is, again, supported by points (2) and (3) above. The most striking observation if one has Eq. (2) in mind, is the practically linear decrease of D ∗ (N) in the log–log plot after a short initial region. I.e., apart from an early bend, which insures that D ∗ (N) 1, the smoothed D ∗ (N) obeys an approximate power law. The fitted powers are between about −0.7 and −0.9, and a detailed analysis shows (for s = 2 and 3, where we have enough data), that the local exponent (fitted separately for each decade) decreases slowly with increasing N . E.g.—and this is typical—for s = 3 and sequence G2 (H1) the exponent decreases from −0.70 (−0.73) in the second decade to −0.73 (−0.74) in the third one and −0.83 (−0.82) in the fourth one. An example can be seen in Fig. 2, where two such fitted lines are included. This slowly increasing steepness resembles, however faintly, ∗ , Eq. (2), for large N , where also power fits over one decade are not too bad with the behaviour of Das ∗ ∗ and Das are similar. powers approaching −1 for N → ∞. But that is also the only trait in which Dtrue The second observation is the small ratio between the discrepancies of different sequences. Indeed, D ∗ for H7 is larger than that for H1 (Figs. 3–6), and Figs. 4 and 6 show a modest further increase when we proceed to much larger bases, but in all investigated dimensions and in the N -range covered the Halton
234
Ch. Schlier / Applied Numerical Mathematics 50 (2004) 227–238
Fig. 4. True discrepancies in 3 dimensions for sequences (see text) G2 (!), G3 (1), H1 (e), G47 (2), H7 (P), and H30 (Q). The latter three ones lie distinctly above the others. Included below are graphs of T ∗ for H1 (long dash) and H7 (short dash). The insert shows the points for G2, G3, H1, and H7 plus the asymptotic discrepancies for these sequences (from top to bottom: ∗ still diverges from the extrapolated true one. H7, H1, G2, G3). Note that the asymptotic Das
sequences H1 have same discrepancy within 20% as the Niederreiter ones G2, even though the factors C(s) for them have a ratio of 2.5, 3.4, 2.3, and 3.8 for s = 2, 3, 4, and 6, respectively. Similarly, from the values of C(s) for different bases the ratios of G2 to the “optimal” Niederreiter sequence for dimensions 3, 4, and 6, i.e., G3, G3, G7, are 2.0, 6.3, and 86.0, respectively, while the true average ratios turn out to be only about 1.13, 1.03, and 1.07 with a rms fluctuation of about 0.15. This almost negligible difference between the average discrepancies of these sequences (in the covered ranges of N ) is striking if we ∗ (N) together with the true points, which is done in the inserted figures. plot the asymptotic curves Das To some extent the basic observation that for practical values of N different sequences behave not as expected from their value C(s) seems to have been made by others, but we found only a single remark in this respect [8]. A third observation is that the ratio between the true and asymptotic curves (i.e., their distance in the ∗ is still increasing within the range of log–log plot) for values of N well behind the maximum of Das our computations. This is true even for s = 2, where our data go more than 104 times farther than the maximum. That leads to an interesting speculation about where the “asymptotic region” begins. A weak but ∗ and the fitted mean true simple definition would be to demand a parallel behaviour of the asymptotic Das ∗ D in the log–log plot. If we take our fits in the last decade of N for which we have data, and extrapolate
Ch. Schlier / Applied Numerical Mathematics 50 (2004) 227–238
235
Fig. 5. True discrepancies in 4 dimensions for sequences (see text) G2 (!), G3 (1), H1 (e), and H7 (P). Only H7 lies distinctly above the others. The insert shows the same points plus the asymptotic discrepancies for the sequences (from top to bottom: H7, H1, G2, G3). Note that the asymptotic curve of G2 lies clearly above the “optimal” G3, which is not the case for the true ∗ still diverges from the extrapolated true one. points. Note further that the asymptotic Das
them to larger values of N , this definition leads to values of N from 105 for s = 2 to 108 for s = 6. For this estimate, we must only assume that the true average D ∗ (N) does not bend upwards in the log–log plot. However, since the data show definite downward bends, we assume that these so defined asymptotic regions begin still much farther out. From our integration tests we would rather propose values of N of 1020 or more for the parameter range covered there, which included dimensions 12 and 24. Finally we would like to call the reader’s attention to the observation that the disparate behaviour of true and asymptotic discrepancies grows strongly with increasing dimension. I.e., while the plotted true discrepancies span an almost constant range for all dimensions treated, the factors between the lowest and the highest plotted asymptotic curves (Gopt respectively H7) increase from 33 for s = 2 to 147 000 for s = 6. From integration errors [17], where, again, the difference of errors between the sequences for increasing dimension remains limited, we conjecture that this behaviour extrapolates to higher dimensions also for true D ∗ .
4. Conclusion We have numerically computed discrepancies of several low discrepancy sequences of the Halton, Niederreiter, Niederreiter–Xing, and Sobol type. Since the number of operations for this computation goes as N s /s!, only a limited range of N can be handled, which becomes very small for large s. With
236
Ch. Schlier / Applied Numerical Mathematics 50 (2004) 227–238
Fig. 6. True discrepancies in 6 dimensions for sequences (see text) G2 (!), G7 (1), G29 (2), H1 (e), H7 (P), and H20 (Q). H7, G29, and H20 lie distinctly above the others. Further down, T ∗ for G2 (long dash) and H1 (short dash) is plotted. The insert shows the points of G2, G7, H1, and H7 plus the asymptotic discrepancies for the sequences (from top to bottom: H7, H1, G2, G7). Note that the asymptotic curve of G2 lies appreciably above the “optimal” G7, which is not the case for the true ∗ still diverges from the extrapolated true one. points. Note further that the asymptotic Das
our resources s = 6 was the largest useful dimension. But we feel justified to extrapolate our results to larger N as well as larger s, because quasi-Monte Carlo integration tests, which are feasible for such larger N and s, show a behaviour of the error quite analogous to that of the true discrepancies, even if the practical error is still much smaller than the discrepancy multiplied with the variation of the integrand. Our essential result is that the presumed asymptotic behaviour of the sequences gives no really useful hint for the true behaviour in regions, which are typical for today’s applications. Especially we find that (1) True discrepancies, when smoothed over local maxima and minima, drop monotonously according to an approximate power law, slightly bending to larger negative powers with N . (2) True discrepancies for different low discrepancy sequences (produced from not too large bases) fall into a narrow range of values, which has no similarity with that suggested by the values of C(s) for these sequences. A detailed selection of low discrepancy sequences on the basis of their asymptotic pre-factor C(s) is, therefore, not justified for practical purposes. Other criteria, such as speed of production or simplicity of programming, can and should then be taken into account. Especially, the larger values of C(s) are no argument against Halton points, which are easy to program and fast to compute. (3) True discrepancies are far below the asymptotic ones, and their ratio, i.e., the logarithmic distance ∗ ∗ and Das increases strongly with dimension. E.g., for s = 2 this distance is 12 to 1 12 between Dtrue
Ch. Schlier / Applied Numerical Mathematics 50 (2004) 227–238
237
orders of magnitude for sequences G2 . . . H7, while it is 2 to 6 orders for s = 6, and we conjecture that it increases further with increasing dimension. Our results give a partial explanation for the behaviour of the errors of integration experiments. E.g., the authors of a textbook [7] wonder why the errors of an example integrated with Niederreiter and Halton sequences are so similar in size, because they expect them to be very different on the basis of C(s). The same was found in other integration tests, also by the author [17] with N up to 1010 trials. It must be repeated, however, that the use of true discrepancies instead of the asymptotic ones in the Koksma– Hlawka inequality removes only some part of the large ratio between the practical integration error and the worst case one given by Eq. (2). So, for typical integrands (i.e., not those specifically constructed to produce the worst case error of a given sequence, cf. [12]), Eq. (2) is not helpful at all. Reliable error estimates should always be deduced from repeated integrations, either with different sequences or with “scrambled” ones [14,6,9]. Similar error estimates are also obtained from repeated integrations with different pieces from a single LDS, which the author intuitively used in recent years [2,3,17]. More precise numerical results for the slowly varying power of the local mean of D ∗ (N) and for its fluctuation are clearly desired, but need much larger computer resources than we had available. In addition, it remains as a task for future theoretical research to find non-asymptotic analytical approximations to the properties of low discrepancy sequences.
Acknowledgement This work was supported by the Deutsche Forschungsgemeinschaft (SFB 238). I am indebted to Peter Bundschuh for sending me a reprint of Ref. [4], which made many of the reported calculations feasible.
References [1] P. Bratley, B.L. Fox, H. Niederreiter, Algorithm 748: Programs to generate Niederreiter’s low-discrepancy sequences, ACM Trans. Math. Software 20 (1994) 494–495. [2] M. Berblinger, Ch. Schlier, Monte Carlo integration with quasi-random numbers: Some experience, Comput. Phys. Comm. 66 (1991) 157–166. [3] M. Berblinger, Ch. Schlier, Monte Carlo integration with quasi-random numbers: Experience with discontinuous integrands, Comput. Phys. Comm. 99 (1997) 151–162. [4] P. Bundschuh, Y. Zhu, A method for exact calculation of the discrepancy of low-dimensional finite point sets, Abh. Math. Sem. Univ. Hamburg 63 (1993) 115–133; also: P. Bundschuh, Y. Zhu, The formulas of exact calculation of the discrepancy of low-dimensional finite point sets, Chinese Sci. Bull. 38 (1993) 1318–1319. [5] J.H. Halton, On the efficiency of certain quasi-random sequences of points in evaluating multidimensional integrals, Numer. Math. 2 (1960) 84–90, and 196. [6] H.S. Hong, F.J. Hickernell, Algorithm 823: Implementing scrambled digital sequences, ACM Trans. Math. Software 29 (2003) 95–109. [7] A.R. Krommer, C.W. Ueberhuber, Numerical Integration on Advanced Computer Systems, Lecture Notes in Comput. Sci., vol. 848, Springer, Berlin, 1994, p. 134. [8] G. Larcher, Digital point sets: Analysis and applications, in: P. Hellekalek, G. Larcher (Eds.), Random and Quasi-Random Point Sets, Lecture Notes in Statist., vol. 138, Springer, Berlin, 1998, pp. 167–222.
238
Ch. Schlier / Applied Numerical Mathematics 50 (2004) 227–238
[9] P. L’Ecuyer, C. Lemieux, Recent advances in randomized quasi-Monte Carlo methods, in: M. Dror, P. L’Ecuyer, F. Szidarovszki (Eds.), Modelling Uncertainty: An Examination of Stochastic Theory, Methods, and Applications, Kluwer Academic, Amsterdam, 2001. [10] W.J. Morokoff, R.E. Caflisch, Quasi-random sequences and their discrepancies, SIAM J. Sci. Comput. 15 (1994) 1251– 1279. [11] H. Niederreiter, Low-discrepancy and low-dispersion sequences, J. Number Theor. 30 (1988) 51–70. [12] H. Niederreiter, Random Number Generation and Quasi-Monte Carlo Methods, SIAM, Philadelphia, PA, 1992. [13] H. Niederreiter, C. Xing, The algebraic-geometry approach to low-discrepancy sequences, in: H. Niederreiter, P. Hellekalek, G. Larcher, P. Zinterhof (Eds.), Monte Carlo and Quasi-Monte Carlo Methods, 1996, Lecture Notes in Statist., vol. 127, Springer, Berlin, 1998, pp. 139–160. [14] A.B. Owen, Randomly permuted (t, m, s)-nets and (t, s) sequences, in: H. Niederreiter, P.J.-S. Shiue (Eds.), Monte Carlo and Quasi-Monte Carlo Methods in Scientific Computing, Lecture Notes in Statist., vol. 106, Springer, Berlin, 1995, pp. 299–315. [15] G. Pirsic, A software implementation of Niederreiter–Xing sequences, in: K.-T. Fang, F.J. Hickernell, H. Niederreiter (Eds.), Monte Carlo and Quasi-Monte Carlo Methods, 2000, Springer, Berlin, 2002, pp. 434–445. [16] W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical Recipes in Fortran, second ed., Cambridge University Press, Cambridge, 1992, Chapter 7.7. [17] Ch. Schlier, A practitioner’s view on QMC integration, Talk given at MCQMC 2001 in Salzburg, Austria, submitted for publication. Available from http://www.phya8.uni-freiburg.de/abt/papers/papers.html. [18] E. Veach, Rubust Monte Carlo Methods for Light Transport Simulation, Ph.D. Thesis, Stanford University, 1997, http://graphics.stanford.edu/papers/veach_thesis.
Discrepancy behaviour in the non-asymptotic regime Ch. Schlier Physikalisches Institut, Universität Freiburg, 79104 Freiburg, Germany Available online 29 January 2004
Abstract We have computed true star-discrepancies D ∗ (N) for Halton and Niederreiter point sequences in dimensions s = 2, 3, 4, and 6 with N up to 250 000, 10 000, 2000 and 300, respectively. For comparison, we also calculated some L2 discrepancies T ∗ . The mean behaviour of D ∗ (N) can well be approximated by power laws with an exponent between about −0.7 and −0.9, which slowly decreases with N. This behaviour is far from the generally assumed asymptotic one ∼ C(s)(log N)s /N. The factor between the true and the asymptotic behaviour increases strongly with s and reaches many orders of magnitude for large s. The ratios of D ∗ (N) for different low discrepancy sequences are not proportional to the presumed asymptotic pre-factors C(s). Especially for the range of bases investigated, Niederreiter sequences have about the same D ∗ as Halton ones in a range of N, which we conjecture to be at least of the order of 1010 . 2003 IMACS. Published by Elsevier B.V. All rights reserved. Keywords: Quasi-Monte Carlo methods; Low-discrepancy sequences; Numerical integration
1. Introduction ∗ ) defined by [12] The star-discrepancy of a s-dimensional point sequence (sometimes also called D∞ N 1 ∗ χJ (xn ) − V (J ), (1) D (N) = sup s N J ⊆[0,1) n=1
is a measure for the uniformity of distribution of the points of the sequence xn , n = 1, 2, . . . in the unit hypercube [0, 1)s . Here, J denotes any subinterval of the s-dimensional unit cube which has faces parallel to the coordinate axes and one corner at the origin, and V (J ) its volume. The characteristic function χJ (n) of the points is 1 if the point is in J , else 0. I.e., D ∗ measures the largest difference E-mail address: [email protected] (Ch. Schlier). 0168-9274/$30.00 2003 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.apnum.2003.12.004
228
Ch. Schlier / Applied Numerical Mathematics 50 (2004) 227–238
between the relative number of points in J and the volume of J . Customarily, a sequence for which an estimate of its discrepancy of the type ∗ (N) = C(s) · (log N)s /N + O (log N)s−1 /N , (2) D ∗ (N) Das is proved is called a low-discrepancy sequence (LDS) [12,8], and it is generally conjectured [12] that the order O(log N)s /N is the best possible for the discrepancy of an infinite sequence. In the sequel, for any low-discrepancy sequence we use the smallest constant C(s), for which such an inequality is proved. The constants C(s) for many LDSs can be found in [12,8], and play an ubiquitous role in the discussion of quasi-Monte Carlo methods. The main reason for this role is the Koksma–Hlawka inequality for the worst case error ε max of quasi-Monte Carlo integration εmax V (f ) · D ∗ (N).
(3)
Here, V (f ) is the variation of the integrand [10,12], a measure of its smoothness, which, however, is unknown for most functions. The pre-factors C(s) in Eq. (2) are generally thought to determine the practical quality of a sequence. However, Eq. (2) describes an asymptotic limit, valid for N → ∞, and the next term of order (log N)s−1 /N is not known. Astonishingly, while everybody mentions Eqs. (2) and (3), nobody knows where, in practice, the asymptotic regime begins. Is it at N = 105 or 1025 or 10100 ? And, strangely enough, there are few, if any, numerical computations of D ∗ (N) for smaller N , i.e., those which are usually used in integration problems. Recently we have performed extensive tests of quasi-Monte Carlo integration with integrands whose variation we could compute [17]. These tests confirmed not only the well-known fact that practical integration errors are many orders of magnitude smaller than Eq. (3) would suggest, but showed also that the relative magnitudes of the errors for different point sequences do not at all behave like the prefactors C(s) predict, at least in a range up to 1010 trials. In Fig. 1 we show an example. The reduced integration error for two test functions from [1], i.e., the error divided by the variation of the integrand, is plotted, which can be directly compared with D ∗ . This led us to ask whether this situation might be a consequence of the non-asymptotic regime of D ∗ , or if everything is due to the difference between practical and worst case errors. The answer is obvious from the figure and discussed below. In the literature one finds formulas for the computation of the root mean square discrepancy T ∗ (also called D2∗ ), defined by replacing the absolute norm in Eq. (1) by the quadratic norm [12]. These algorithms have an expense ∼ s 2 N 2 , so we could compute T ∗ to rather large N , e.g., for dimension 6 up to N = 106 . For all sequences and dimensions which we investigated the functions T ∗ (N) observe an approximate power law in the non-asymptotic region. Unfortunately, there is no known way to draw useful conclusions from T ∗ to D ∗ , which is needed in Eq. (3). So, even if this behaviour of T ∗ could be a hint that D ∗ might behave similarly, this is not proved theoretically. Since obviously D ∗ (1) < 1, and intuitively one would expect D ∗ to decrease with N , it is trivial that D ∗ (N) for small N must behave differently from the asymptotic behaviour. But, again, where does the asymptotic regime begin? A handwaiving argument says that it must be much farther out than the maximum of the asymptotic D ∗ (N) in Eq. (2), which is at es . But how much? A more restrictive argument [18] maintains that D ∗ (N) should be at least < D ∗ (2) before the asymptotic region can begin. This, then, would mean that the limiting N is about 3000 for s = 3, but already 1.6 × 109 for s = 6, increasing to computationally unfeasible values for larger s. So we started a project to compute D ∗ for several Halton, Niederreiter, Niederreiter–Xing and Sobol sequences, which we describe in the next sections. The maximum N for which we could compute exact
Ch. Schlier / Applied Numerical Mathematics 50 (2004) 227–238
229
Fig. 1. Asymptotic discrepancy (thick smooth curve) and true discrepancy (dots with thin fitted line) of the Halton sequence in 6 dimensions (primes 2 to 13) together with reduced errors of integration tests (dashed lines with dashed fits) performed with this sequence. Plotted is the root mean square error of repeated integrations of test functions 1 (lower) and 4 (upper) from [1] divided by the integrand’s variation (data from [17]). The heavy dots show the square root of the variances of the test functions. These dots are the origins of the expectation value of the reduced errors of pseudo-Monte Carlo integration. Note that the asymptotic discrepancy overestimates the true one appreciably, but the practical reduced errors are still much lower than their Koksma–Hlawka limit.
values of D ∗ was limited by our computer resources to about 250 000, 10 000, 2500, and 300 for s = 2, 3, 4, and 6, respectively. (We had to progress in steps > 1 for the larger values of N .) Niederreiter and Halton sequences were computed with different bases respectively base sets. For these s, and some larger ones we have also computed approximate D ∗ . The convergence of our simple approximation to D ∗ for large s and N is, however, very bad. In view of the many orders of magnitude difference from the asymptotic values they are, however, quite useful even with their errors. Our main result can be summed up in two statements: (1) D ∗ (N) behaves approximately according to a power law with very slowly decreasing negative power for those N which determine practical integration results, and (2) the relative magnitude of D ∗ (N) for different low-discrepancy sequences does not behave as the pre-factors C(s) of the asymptotic behaviour predict. Especially, at practical values of N and for the range of bases investigated, D ∗ (N) does not justify a preference of Niederreiter sequences over Halton ones. As a consequence for the practitioner, we propose that he should feel free to choose the simpler and faster available Halton sequences and also other LDSs considered “non-optimal” in the literature (e.g., Niederreiter sequences to base 2 even if this is not the “optimal” base) without fear of large penalties in performance. As a consequence, repeated integrations with different LDSs can be recommended as a method to increase the overall confidence in the results.
230
Ch. Schlier / Applied Numerical Mathematics 50 (2004) 227–238
2. Calculation methods We have used one approximate and two exact algorithms to compute D ∗ . Our first idea was to replace the supremum over all hyper-rectangles in Eq. (1) by the maximum over a large number, M, of randomly ∗ (N, M) chosen ones. The computational effort is s · N · M. Unfortunately, the convergence of this Dappr with increasing M is not good. Some convergence plots can be found in [17]. Actually, one can argue that the number of hyper-rectangles needed here must be larger than that which one must visit in the exact programs described below. On the other hand, this algorithm is very simple, and can, therefore, easily be checked for correctness. So we could use it with small N (say < 20) and large M (up to 1010 ) to check, at least punctually, the correctness of the implementation of the two exact algorithms described below, whose structures are much more complex. Also, if one accepts an error of one order of magnitude, which is not much compared to the many orders difference from the asymptotic behaviour, one can use even these approximate values for a crude estimate of the true ones at larger N , cf. Fig. 2. This figure includes ∗ (N) and T ∗ (N), and shows that they are approximately parallel. This was found also in the plots of Dappr other cases, and means that these two types of discrepancy differ only by a slowly changing factor. For an exact computation of D ∗ , one must consider all hyper-rectangles at whose corners the number of included points may possibly change. For Halton and Niederreiter (and other) low-discrepancy sequences these lie on a grid, whose size can be computed. In the simple Halton case it is the product of (in Mathematica notation) pk Ceiling[Log[pk , N]] , where the pk are the primes used in the construction
Fig. 2. True discrepancies in 3 dimensions for sequences (see text) G2 (!) and H1 (e). Thin straight lines: linear fits of log10 (D ∗ ) vs. log10 (N) of G2 points in the 2nd, respectively, 4th decade of N . Thick lines: approximate D ∗ for G2 (long dash) and H1 (short dash). Thin lines: T ∗ for G2 (long dash) and H1 (short dash).
Ch. Schlier / Applied Numerical Mathematics 50 (2004) 227–238
231
of the sequence. This is of the order of s! · N s . For each such possible hyper-rectangle one must check which points of the sequence are contained in it, which leads to an effort of s! · N s+1 for this algorithm. We have used it in our older calculations [17], but could only handle N = 1000 and 100 for s = 2 and 3, respectively. A better algorithm has been published for s = 2 and 3 in a difficult to find paper [4]. It can easily be extended to other dimensions, and we have done so for s = 3 to 12. The idea is to investigate only those hyper-rectangles at whose corners the number of included sequence points does change. The algorithm can best be described by an excerpt from our Mathematica implementation, here in 3 dimensions, the extension to other dimensions being obvious. Note that Mathematica does not allow 0 as index, so all indexes are increased by one compared to mathematical notation. Note also that Union[] includes sorting: ptab =Union[{{0,0,0}}, Take[pointlist,n]]; xtab = Union[Map[ex31, ptab], {1}]; (∗ ex31 extracts the first component of three ∗ ) mx1 = mx2 = 0; For[l = 0, l n, l++, yta = Map[ex323, Take[ptab, l+1]]//Sort; (∗ ex323 extracts components 2 and 3 of three ∗ ) ytab = Union[Map[ex21, yta],{1}]; (∗ ex21 extracts the 1st of two components ∗ ) mx10 = mx20 = 0; For[k = 0, k l, k++, zta = Map[ex22, Take[yta, k+1]]//Sort; (∗ ex22 extracts the 2nd of two components ∗ ) ztab = Union[zta,{1}]; mx100 = mx200 = 0; For[ j = 0, j k, j++, taba = Abs[j/n - xtab[[l+1]]∗ytab[[k+1]]∗ztab[[j+1]] ]; tabb = Abs[xtab[[l+2]]∗ytab[[k+2]] ∗ ztab[[j+2]]- j/n ]; mx100 = Max[Abs[taba],mx100]; mx200 = Max[Abs[tabb],mx200]; ]; mx10 = Max[mx10, mx100]; mx20 = Max[mx20, mx200]; ]; mx1= Max[mx1,mx10]; mx2 = Max[mx2,mx20]; ]; mx = Max[mx1,mx2]; The input is a point list of length > n of the points of the—in this case 3-dimensional—sequence, and the output is mx = D ∗ (n) as rational number. The sequel of projections on hyper-rectangles of ever lower dimension restricts the grid points which must be visited to those where the number of included points of the sequence actually changes, and the kind of indexing makes this number known. This leads to an asymptotic number of N s /s! hyper-rectangles, which must be investigated. The limit on single values of N for our equipment (four 450 MHz Xeons, computations up to one week accepted) was about 35 000, 1500, 350, and 90 for s = 2, 3, 4, and 6, respectively. An implementation of this algorithm in Fortran 90 (now with double precision “real” numbers instead of integers) speeds it up by about three orders of magnitude. An example program for s = 5 can be found
232
Ch. Schlier / Applied Numerical Mathematics 50 (2004) 227–238
on the author’s web page http://www.phya8.uni-freiburg.de/abt/papers/papers.html. With this program we could compute single values of D ∗ (N) up to N = 250 000, 10 000, 2500 and 300 for dimensions s = 2, 3, 4, and 6, respectively. Another factor of 10 to 100 might be gained with one or several faster processors. But note that 1001/6 is only 2.15! So, the calculation of D ∗ in higher dimensions for practically interesting N is still unfeasible today.
3. Results We have computed exact discrepancies in 2, 3, 4, and 6 dimensions for the following sequences: Halton points [5] produced from consecutive primes starting from the first (2), the third (5), or the 7th (17), abbreviated as H1, H3, and H7, points from Niederreiter’s (t, s)-sequences [11], Niederreiter points for short, with bases 2, 3, 5, and 7, abbreviated G2, G3, G5, and G7. In addition we computed some sequences with much larger bases, Halton points H20 and H30 (same notation as above) and Niederreiter points G29 and G47, and finally some Niederreiter–Xing and Sobol sequences. Halton points were produced by our simple algorithm [2], Niederreiter ones by the TOMS algorithm 728 [1], Niederreiter– Xing points using [3,15], and Sobol points by the subroutine in Ref. [16]. As we discussed above, the maximum feasible N depends on the dimension. To avoid clutter, we plot only the Halton results for H1, and H7, those for H3 lying generally between them. For Niederreiter sequences we restrict our plots to only G2 and the so called “optimal” sequence, i.e., that with the smallest C(s) for the given dimension. These are G2, G3, G3, and G7 for the four dimensions involved. For s = 2 we included also G7. Note, however, that on today’s binary computers any Niederreiter sequence other than G2 is much more slowly produced. We do not plot Niederreiter–Xing or Sobol discrepancies, since they are on the average not much different from G2. The computed values are depicted in Figs. 3–6 as points. For small N we stepped N by 1, while for higher N larger steps had to be used. D ∗ (N) is not a monotonously decreasing function but fluctuates about some average trend. Local maxima and minima follow each other in short runs, and are, therefore, not too interesting for applications. These fluctuations are smallest for the sequences with small bases, as one would perhaps expect from pictures of two-dimensional projections of the sequence points, which do not show large “holes” in this case. Their relative amplitude seems to be almost constant, but we could not really ascertain this. Support for this statement comes, however, from a short test (in dimension 3, with N = 1 . . 1500) performed with a series of ten sequences of G2 type scrambled according to the prescription of [6], which is a simplified version of Owen’s scrambling [14]. As expected, the average discrepancy of this set fluctuates much less than that of a single one. In contrast, the relative standard deviation of the Di∗ (N) (i.e., the standard deviation divided by the mean) decreases very slowly with N . In the following we discuss several observations, which the reader can verify from the figures. Extrapolation from the computed data to higher values of N lies, however, at hand. It seems justified for several reasons: (1) From the construction of a low-discrepancy sequence one certainly does not expect a sudden break in its mean behaviour. ∗ and for the (2) The same general behaviour as for the true D ∗ is also observed for the approximate Dappr ∗ L2 -discrepancy T , in which cases we could compute data for much larger N .
Ch. Schlier / Applied Numerical Mathematics 50 (2004) 227–238
233
Fig. 3. True discrepancies in 2 dimensions for sequences (see text) G2 (!), G7 (1), H1 (e), and H7 (P). Points of H7 and later also G7 lie distinctly above the others. The insert shows the same points plus the asymptotic discrepancies for the sequences (from top to bottom: H7, G7, H1, G2).
(3) This behaviour is also observed for the error of test integrations of functions with known variation, which we computed up to N = 1010 trials [17] (cf. Fig. 1). Indeed, it was these test integrations which led use to the investigation of D ∗ . We also think that some extrapolation to larger dimensions is allowed, which is, again, supported by points (2) and (3) above. The most striking observation if one has Eq. (2) in mind, is the practically linear decrease of D ∗ (N) in the log–log plot after a short initial region. I.e., apart from an early bend, which insures that D ∗ (N) 1, the smoothed D ∗ (N) obeys an approximate power law. The fitted powers are between about −0.7 and −0.9, and a detailed analysis shows (for s = 2 and 3, where we have enough data), that the local exponent (fitted separately for each decade) decreases slowly with increasing N . E.g.—and this is typical—for s = 3 and sequence G2 (H1) the exponent decreases from −0.70 (−0.73) in the second decade to −0.73 (−0.74) in the third one and −0.83 (−0.82) in the fourth one. An example can be seen in Fig. 2, where two such fitted lines are included. This slowly increasing steepness resembles, however faintly, ∗ , Eq. (2), for large N , where also power fits over one decade are not too bad with the behaviour of Das ∗ ∗ and Das are similar. powers approaching −1 for N → ∞. But that is also the only trait in which Dtrue The second observation is the small ratio between the discrepancies of different sequences. Indeed, D ∗ for H7 is larger than that for H1 (Figs. 3–6), and Figs. 4 and 6 show a modest further increase when we proceed to much larger bases, but in all investigated dimensions and in the N -range covered the Halton
234
Ch. Schlier / Applied Numerical Mathematics 50 (2004) 227–238
Fig. 4. True discrepancies in 3 dimensions for sequences (see text) G2 (!), G3 (1), H1 (e), G47 (2), H7 (P), and H30 (Q). The latter three ones lie distinctly above the others. Included below are graphs of T ∗ for H1 (long dash) and H7 (short dash). The insert shows the points for G2, G3, H1, and H7 plus the asymptotic discrepancies for these sequences (from top to bottom: ∗ still diverges from the extrapolated true one. H7, H1, G2, G3). Note that the asymptotic Das
sequences H1 have same discrepancy within 20% as the Niederreiter ones G2, even though the factors C(s) for them have a ratio of 2.5, 3.4, 2.3, and 3.8 for s = 2, 3, 4, and 6, respectively. Similarly, from the values of C(s) for different bases the ratios of G2 to the “optimal” Niederreiter sequence for dimensions 3, 4, and 6, i.e., G3, G3, G7, are 2.0, 6.3, and 86.0, respectively, while the true average ratios turn out to be only about 1.13, 1.03, and 1.07 with a rms fluctuation of about 0.15. This almost negligible difference between the average discrepancies of these sequences (in the covered ranges of N ) is striking if we ∗ (N) together with the true points, which is done in the inserted figures. plot the asymptotic curves Das To some extent the basic observation that for practical values of N different sequences behave not as expected from their value C(s) seems to have been made by others, but we found only a single remark in this respect [8]. A third observation is that the ratio between the true and asymptotic curves (i.e., their distance in the ∗ is still increasing within the range of log–log plot) for values of N well behind the maximum of Das our computations. This is true even for s = 2, where our data go more than 104 times farther than the maximum. That leads to an interesting speculation about where the “asymptotic region” begins. A weak but ∗ and the fitted mean true simple definition would be to demand a parallel behaviour of the asymptotic Das ∗ D in the log–log plot. If we take our fits in the last decade of N for which we have data, and extrapolate
Ch. Schlier / Applied Numerical Mathematics 50 (2004) 227–238
235
Fig. 5. True discrepancies in 4 dimensions for sequences (see text) G2 (!), G3 (1), H1 (e), and H7 (P). Only H7 lies distinctly above the others. The insert shows the same points plus the asymptotic discrepancies for the sequences (from top to bottom: H7, H1, G2, G3). Note that the asymptotic curve of G2 lies clearly above the “optimal” G3, which is not the case for the true ∗ still diverges from the extrapolated true one. points. Note further that the asymptotic Das
them to larger values of N , this definition leads to values of N from 105 for s = 2 to 108 for s = 6. For this estimate, we must only assume that the true average D ∗ (N) does not bend upwards in the log–log plot. However, since the data show definite downward bends, we assume that these so defined asymptotic regions begin still much farther out. From our integration tests we would rather propose values of N of 1020 or more for the parameter range covered there, which included dimensions 12 and 24. Finally we would like to call the reader’s attention to the observation that the disparate behaviour of true and asymptotic discrepancies grows strongly with increasing dimension. I.e., while the plotted true discrepancies span an almost constant range for all dimensions treated, the factors between the lowest and the highest plotted asymptotic curves (Gopt respectively H7) increase from 33 for s = 2 to 147 000 for s = 6. From integration errors [17], where, again, the difference of errors between the sequences for increasing dimension remains limited, we conjecture that this behaviour extrapolates to higher dimensions also for true D ∗ .
4. Conclusion We have numerically computed discrepancies of several low discrepancy sequences of the Halton, Niederreiter, Niederreiter–Xing, and Sobol type. Since the number of operations for this computation goes as N s /s!, only a limited range of N can be handled, which becomes very small for large s. With
236
Ch. Schlier / Applied Numerical Mathematics 50 (2004) 227–238
Fig. 6. True discrepancies in 6 dimensions for sequences (see text) G2 (!), G7 (1), G29 (2), H1 (e), H7 (P), and H20 (Q). H7, G29, and H20 lie distinctly above the others. Further down, T ∗ for G2 (long dash) and H1 (short dash) is plotted. The insert shows the points of G2, G7, H1, and H7 plus the asymptotic discrepancies for the sequences (from top to bottom: H7, H1, G2, G7). Note that the asymptotic curve of G2 lies appreciably above the “optimal” G7, which is not the case for the true ∗ still diverges from the extrapolated true one. points. Note further that the asymptotic Das
our resources s = 6 was the largest useful dimension. But we feel justified to extrapolate our results to larger N as well as larger s, because quasi-Monte Carlo integration tests, which are feasible for such larger N and s, show a behaviour of the error quite analogous to that of the true discrepancies, even if the practical error is still much smaller than the discrepancy multiplied with the variation of the integrand. Our essential result is that the presumed asymptotic behaviour of the sequences gives no really useful hint for the true behaviour in regions, which are typical for today’s applications. Especially we find that (1) True discrepancies, when smoothed over local maxima and minima, drop monotonously according to an approximate power law, slightly bending to larger negative powers with N . (2) True discrepancies for different low discrepancy sequences (produced from not too large bases) fall into a narrow range of values, which has no similarity with that suggested by the values of C(s) for these sequences. A detailed selection of low discrepancy sequences on the basis of their asymptotic pre-factor C(s) is, therefore, not justified for practical purposes. Other criteria, such as speed of production or simplicity of programming, can and should then be taken into account. Especially, the larger values of C(s) are no argument against Halton points, which are easy to program and fast to compute. (3) True discrepancies are far below the asymptotic ones, and their ratio, i.e., the logarithmic distance ∗ ∗ and Das increases strongly with dimension. E.g., for s = 2 this distance is 12 to 1 12 between Dtrue
Ch. Schlier / Applied Numerical Mathematics 50 (2004) 227–238
237
orders of magnitude for sequences G2 . . . H7, while it is 2 to 6 orders for s = 6, and we conjecture that it increases further with increasing dimension. Our results give a partial explanation for the behaviour of the errors of integration experiments. E.g., the authors of a textbook [7] wonder why the errors of an example integrated with Niederreiter and Halton sequences are so similar in size, because they expect them to be very different on the basis of C(s). The same was found in other integration tests, also by the author [17] with N up to 1010 trials. It must be repeated, however, that the use of true discrepancies instead of the asymptotic ones in the Koksma– Hlawka inequality removes only some part of the large ratio between the practical integration error and the worst case one given by Eq. (2). So, for typical integrands (i.e., not those specifically constructed to produce the worst case error of a given sequence, cf. [12]), Eq. (2) is not helpful at all. Reliable error estimates should always be deduced from repeated integrations, either with different sequences or with “scrambled” ones [14,6,9]. Similar error estimates are also obtained from repeated integrations with different pieces from a single LDS, which the author intuitively used in recent years [2,3,17]. More precise numerical results for the slowly varying power of the local mean of D ∗ (N) and for its fluctuation are clearly desired, but need much larger computer resources than we had available. In addition, it remains as a task for future theoretical research to find non-asymptotic analytical approximations to the properties of low discrepancy sequences.
Acknowledgement This work was supported by the Deutsche Forschungsgemeinschaft (SFB 238). I am indebted to Peter Bundschuh for sending me a reprint of Ref. [4], which made many of the reported calculations feasible.
References [1] P. Bratley, B.L. Fox, H. Niederreiter, Algorithm 748: Programs to generate Niederreiter’s low-discrepancy sequences, ACM Trans. Math. Software 20 (1994) 494–495. [2] M. Berblinger, Ch. Schlier, Monte Carlo integration with quasi-random numbers: Some experience, Comput. Phys. Comm. 66 (1991) 157–166. [3] M. Berblinger, Ch. Schlier, Monte Carlo integration with quasi-random numbers: Experience with discontinuous integrands, Comput. Phys. Comm. 99 (1997) 151–162. [4] P. Bundschuh, Y. Zhu, A method for exact calculation of the discrepancy of low-dimensional finite point sets, Abh. Math. Sem. Univ. Hamburg 63 (1993) 115–133; also: P. Bundschuh, Y. Zhu, The formulas of exact calculation of the discrepancy of low-dimensional finite point sets, Chinese Sci. Bull. 38 (1993) 1318–1319. [5] J.H. Halton, On the efficiency of certain quasi-random sequences of points in evaluating multidimensional integrals, Numer. Math. 2 (1960) 84–90, and 196. [6] H.S. Hong, F.J. Hickernell, Algorithm 823: Implementing scrambled digital sequences, ACM Trans. Math. Software 29 (2003) 95–109. [7] A.R. Krommer, C.W. Ueberhuber, Numerical Integration on Advanced Computer Systems, Lecture Notes in Comput. Sci., vol. 848, Springer, Berlin, 1994, p. 134. [8] G. Larcher, Digital point sets: Analysis and applications, in: P. Hellekalek, G. Larcher (Eds.), Random and Quasi-Random Point Sets, Lecture Notes in Statist., vol. 138, Springer, Berlin, 1998, pp. 167–222.
238
Ch. Schlier / Applied Numerical Mathematics 50 (2004) 227–238
[9] P. L’Ecuyer, C. Lemieux, Recent advances in randomized quasi-Monte Carlo methods, in: M. Dror, P. L’Ecuyer, F. Szidarovszki (Eds.), Modelling Uncertainty: An Examination of Stochastic Theory, Methods, and Applications, Kluwer Academic, Amsterdam, 2001. [10] W.J. Morokoff, R.E. Caflisch, Quasi-random sequences and their discrepancies, SIAM J. Sci. Comput. 15 (1994) 1251– 1279. [11] H. Niederreiter, Low-discrepancy and low-dispersion sequences, J. Number Theor. 30 (1988) 51–70. [12] H. Niederreiter, Random Number Generation and Quasi-Monte Carlo Methods, SIAM, Philadelphia, PA, 1992. [13] H. Niederreiter, C. Xing, The algebraic-geometry approach to low-discrepancy sequences, in: H. Niederreiter, P. Hellekalek, G. Larcher, P. Zinterhof (Eds.), Monte Carlo and Quasi-Monte Carlo Methods, 1996, Lecture Notes in Statist., vol. 127, Springer, Berlin, 1998, pp. 139–160. [14] A.B. Owen, Randomly permuted (t, m, s)-nets and (t, s) sequences, in: H. Niederreiter, P.J.-S. Shiue (Eds.), Monte Carlo and Quasi-Monte Carlo Methods in Scientific Computing, Lecture Notes in Statist., vol. 106, Springer, Berlin, 1995, pp. 299–315. [15] G. Pirsic, A software implementation of Niederreiter–Xing sequences, in: K.-T. Fang, F.J. Hickernell, H. Niederreiter (Eds.), Monte Carlo and Quasi-Monte Carlo Methods, 2000, Springer, Berlin, 2002, pp. 434–445. [16] W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical Recipes in Fortran, second ed., Cambridge University Press, Cambridge, 1992, Chapter 7.7. [17] Ch. Schlier, A practitioner’s view on QMC integration, Talk given at MCQMC 2001 in Salzburg, Austria, submitted for publication. Available from http://www.phya8.uni-freiburg.de/abt/papers/papers.html. [18] E. Veach, Rubust Monte Carlo Methods for Light Transport Simulation, Ph.D. Thesis, Stanford University, 1997, http://graphics.stanford.edu/papers/veach_thesis.