CC BY
9
Discrete & Continuous Models
#& Applied Computational Science 2021, 29 (3) 260-270
ISSN 2658-7149 (online), 2658-4670 (print) http://journals-rudn-ru/miph
Research article
UDC 519.872:519.217
PACS 07.05.Tp, 02.60.Pn, 02.70.Bf
DOI: 10.22363/2658-4670-2021-29-3-260-270
Shifted Sobol points and multigrid Monte Carlo
simulation
Aleksandr A. Belov1'2, Maxim A. Tintul1
1 M. V. Lomonosov Moscow State University 1, bld. 2, Leninskie Gory, Moscow, 119991, Russian Federation 2 Peoples' Friendship University of Russia (RUDN University) 6, Miklukho-Maklaya St., Moscow, 117198, Russian Federation
(received: August 19, 2021; accepted: September 9, 2021)
Multidimensional integrals arise in many problems of physics. For example, moments of the distribution function in the problems of transport of various particles (photons, neutrons, etc.) are 6-dimensional integrals. When calculating the coefficients of electrical conductivity and thermal conductivity, scattering integrals arise, the dimension of which is equal to 12. There are also problems with a significantly large number of variables. The Monte Carlo method is the most effective method for calculating integrals of such a high multiplicity. However, the efficiency of this method strongly depends on the choice of a sequence that simulates a set of random numbers. A large number of pseudo-random number generators are described in the literature. Their quality is checked using a battery of formal tests. However, the simplest visual analysis shows that passing such tests does not guarantee good uniformity of these sequences. The magic Sobol points are the most effective for calculating multidimensional integrals. In this paper, an improvement of these sequences is proposed: the shifted magic Sobol points that provide better uniformity of points distribution in a multidimensional cube. This significantly increases the cubature accuracy. A significant difficulty of the Monte Carlo method is a posteriori confirmation of the actual accuracy. In this paper, we propose a multigrid algorithm that allows one to find the grid value of the integral simultaneously with a statistically reliable accuracy estimate. Previously, such estimates were unknown. Calculations of representative test integrals with a high actual dimension up to 16 are carried out. The multidimensional Weierstrass function, which has no derivative at any point, is chosen as the integrand function. These calculations convincingly show the advantages of the proposed methods.
Key words and phrases: multidimensional integral, Monte Carlo method, Sobol points, multigrid calculation, a posteriori error estimates
1. Introduction
Integrals of multivariate functions occur in many areas of physics. Here are some examples. The transfer of neutrons, photons and other particles in
© Belov A. A., Tintul M.A., 2021
This work is licensed under a Creative Commons Attribution 4.0 International License http://creativecommons.org/licenses/by/4.0/
the medium is described by the equation for the distribution function; this function depends on three coordinates of the medium and three components of the particle velocity vector, that is, the number of variables is six. To determine the coefficients of thermal conductivity or electrical conductivity of a medium, it is necessary to calculate the collision integrals; they include components of the velocity vectors before the moment of collision and after the moment of collision. The total number of variables in such an integral is twelve. Problems also arise with a significantly larger number of variables.
In the simplest formulation, the calculation of the integral in the unit p-dimensional cube Vis considered. x — (xx,x2, ■■■, xp) is p-dimensional vector. Our aim is to calculate the following integral:
I — J ^ (x)dx — J ••• J , X2, •••, xp)dx~idx2 ••• dxp.
V 0 0
The accuracy of numerical grid methods drops rapidly with the increase of dimension p. In order to obtain acceptable accuracy, more and more points have to be taken, which makes the calculations exorbitantly laborious and very time consuming. Due to this fact, the local Monte Carlo method is used for high dimensions (p > 3). It involves the use of random numbers, which are mathematical abstraction. In practice, however, one has to use sequences that only imitate random numbers. Performance of the method strongly depends on the choice of such a sequence.
Calculations of the representative test integrals show that to obtain good accuracy the most important is the uniformity of the points' distribution and not its randomness. The most effective are Sobol sequences with the so-called "magic" numbers of points N — 2n , n — 0,1,....
In this work, the following results are obtained. Firstly, shifted Sobol points are proposed. It is a modification that improves uniformity of the point distribution and increases the accuracy of cubatures. Secondly, a multi-grid strategy that gives a posteriori estimate of the accuracy is constructed. The advantages of the proposed algorithms are illustrated with representative test examples.
2. Pseudorandom points
For the local Monte Carlo method, N random points Xj are selected in the cube V; in this case, the number N can be arbitrary, in contrast to cubature formulae on regular grids. The cubature formula
1 N
is similar to the formula for mean Riemann sum. However, the estimate of its error AN turns out to be radically different
AN — I-In -1/2, Df — J f2 (x)dx—
v
j f(x)dx
(2)
Here Df is variance. The estimate of the error is not majorant, but probabilistic: magnitude of the error is distributed according to the Gaussian law with the standard specified in the formula. The error does not exceed the standard deviation with a probability of 0.68.
The error estimate (2) does not depend on the dimension p. Random points are inferior in accuracy to regular grids at p = 1 or p = 2. Already at p = 4, the dependence of the error on N for random points and regular grids is the same. With further increase in dimension, random points turns out to be more advantageous; advantage increases rapidly as dimension p grows.
Formulae (2) assume that random points Xj have uniform distribution density in the cube Vand are not correlated. However, no rigorous mathematical methods for constructing such points have been found. A number of mathematical algorithms have been proposed; the resulting points are called pseudorandom. An extensive literature is devoted to the construction of pseudorandom points, for example, [1]-[13]. The following generators are most common in the literature:
— Mersenne twister and SIMD-oriented fast Mersenne twister;
— Multiplicative congruential generator;
— 64-bit multiplicative lagged Fibonacci generator;
— combined multiple recursive generator;
— generator Philo4x32;
— generator Threefry4x64;
— Marsaglia's SHR3 shift-register generator;
— modified Subtract-with-Borrow generator;
— modified Lehmer sequence.
These generators are implemented in many commercial packages (for example, Matlab).
The quality of each sequence of pseudorandom numbers is checked using some sets of tests based on the theory of probability [14]-[17]. But no set of tests can be complete and comprehensive. Therefore, such checks are limited. Even the simplest visual tests show that widespread sequences do not provide a sufficiently good uniformity of filling the unit square [18], [19]. The question of the influence of such unevenness on the actual error of cubatures remains insufficiently clear.
3. Sobol points
To construct the Sobol sequence, a set of the so-called direction numbers should be selected. There is some ambiguity in the selection of initial direction numbers. In early works [1], direction number tables were constructed for dimensions p < 13 and numbers n < 20 (total number of points N < 220). Later, direction numbers for higher p and N were constructed [20]. However, the direction numbers were also changed. The program is currently available at 21. The open access option contains p < 50 and n<31 (^«2-109). The commercial version of the program has p < 216 — 1.
It is important to note that the Sobol sequences are constructed separately for each p. It is impossible to obtain a sequence of fewer dimensions from p-dimensional Sobol sequence. This also applies to magic segments of the Sobol sequences.
The Sobol cubature formula has the same form as (2). But the estimate of its error is not entirely clear. The distribution of points only for magic N approaches uniform in properties. For intermediate N, it is obtained by discarding some of the points and loses the property of uniformity. Therefore, only magic N should be used for cubatures.
Various attempts have been made to generalize the Sobol sequences. However, the search for optimal variants of such generalizations invariably led again to the Sobol sequences. Therefore, such generalizations need to be treated with caution.
4. Shifted Sobol points
The arrangement of the Sobol points is somewhat asymmetrical. For example, if number of points N — 2n is taken, then the arithmetic mean of all points projections on any axis will not be 0.5, but 0.5 (1 — 1/N). Obviously, this asymmetry is not favourable for obtaining good cubature accuracy.
In the figure 1, black circles show two-dimensional Sobol points for the first magic numbers. For n — 0, the only point lies in the corner of the unit square. Calculation of the cubature over this point gives a formula of the first order accuracy. However, if this point is shifted by 0.5 along each coordinate, then the cubature over the shifted point (light circle) has the second order of accuracy. For the case n — 1, two points are located one in the corner of the square and one in the center, which will also give the first order of accuracy. But if these two points are shifted by 0.25 along each coordinate, then the cubature error obviously decreases. Therefore, a general shift principle for any number of dimensions can be proposed:
If N — 2n, then add to all coordinates of all points (2N).
It is advisable to apply this shift only for magic Sobol numbers. In this case, the shifts are different for different N.
5. Multigrid calculation
Test calculations show that the actual error decreases as O (N-1). This suggests that it is possible to approximate the integral (and hence its error) as a function of N. However, this approximation cannot be smooth, such as Richardson's interpolation approximation for grid methods. In this case, the points are obtained by statistical methods, therefore, their processing must be carried out using the root-mean-square approximation. To do this, the type of approximation must be chosen and some weights to the points need to be assigned.
As a working hypothesis, the law of decreasing error AN ~ N-1 was assumed. But since the nature of the error becomes clearly statistical with increasing p, the standard deviation of these errors was assumed to be proportional to N-1/2. This is the weight used for approximation.
The following multigrid procedure is proposed. The calculation with magic N — 2n, n — 10,11,... is performed. As a result, a sequence of values of the integral [IN} is obtained. Now this sequence can be approximated by the method of least squares
IN « a + bN-1. (3)
n=0 (N=1)
n=l (N=2)
n=2 (N=4)
n=3 (N=8)
A
o
«
«
..........O..........
r\
. r\
o
r\ .
A
r\
o
Figure 1. Sobol magic points for p = 2: points - unbiased, circles - shifted; the n values are indicated near the squares
Here a is the refined value of the integral. At the same time, the standard deviation aa for the value a is calculated. This standard deviation is a statistical estimate of the accuracy for the found value of the integral.
Note that the beginning of the sequence [IN} corresponding to n = 0,1,..., 9 is not taken into account in approximation (3), since these grids are not detailed enough, and the rate of decrease of the error does not yet correspond to 0(N-1).
6. Test integral
It is expedient to carry out numerical experiments on multidimensional integrals over the unit cube, the exact values of which are known. Then the error of the numerical calculation can be directly determined and its behaviour can be studied. Further, requirements that are appropriate for the integrand are discussed.
In multidimensional problems, the concept of the effective dimension of a function is used. For example, consider two functions:
f(x) = n fj (*i) (4)
3=1
and
( p
f(x) = fl ( XJ
\j=1
)
where all f^(xj) are essentially different from constants. In the first function, all variables are equally important, and the effective dimension of the function is p. The second function depends on only one combination of variables, so its effective dimension is 1. The higher the effective dimension of the function, the more difficult the problem. Therefore, the most difficult functions are of the first type.
Suppose that for a product function each f^ differs substantially from zero only on a segment of length ft of its unit edge. Then the product of one-dimensional functions will differ significantly from zero in the volume ftp. If ft is small, then as p increases, the volume ftp decreases rapidly; for example, for ft = 0.1 and p = 10 the value ftp = 10-10. In this case, to obtain acceptable accuracy, any Monte Carlo method will require the number of nodes N ft-p. It can be seen that in order for the number of points to be reasonable, ft should be taken close to one.
Taking these considerations into account, a test of the form (4) have been chosen. It is not easy, despite its seeming simplicity. All f^ are assumed to be the same and equal to the Weierstrass functions
where a is an arbitrary odd number that is not equal to one, and b is a positive number less than one. It is known that under the conditions ab > 1, a > 1, the Weierstrass function is continuous, but has no derivative at any point. This test is extremely difficult. The Weierstrass function is shown in the figure 2.
CO
(5)
1 -
> 0-
A
v.
-1 -
0,0
0.1
0,2
0,3
0,4
0,5
X
Figure 2. Weierstrass function with a = 3,6 = 0.5
Taking into account the symmetry of the Weierstrass function, the integration is carried out over a cube with sides xrj G [0,0.5]. For convenience, the Weierstrass function is normalized, the normalization condition is
J f(x)dx = 1. (6)
y
7. Calculation results
The integral of the multidimensional Weierstrass function (4), (5) was calculated using three qualitatively different approaches: regular cubature on trapezoidal formulae, the classical Monte Carlo method using the Mersenne twister and shifted Sobol magic points.
These three approaches are compared in terms of the error magnitude with a fairly modest number of points N = 220. The logarithms of the errors depending on the dimension are shown in the figure 3. Let us analyze the curves.
P
Figure 3. Logarithm of the relative error in calculating the integral of the Weierstrass function for N = 220: light triangle is AN, circle is aa for the shifted Sobol points, black inverted triangle corresponds to Mersenne twister, black square is for trapezoidal method
Mersenne twister. Beginning with dimension p = 11, the curve corresponding to the Mersenne twister lies below all. Despite the good accuracy, there are no means to confirm it. An attempt to apply the root-mean-square approximation (3) to the Mersenne twister was unsuccessful: the values of aa turn out to be either larger or smaller than the actual error depending on the dimension p, and the difference can be significant.
1/2
The standard deviation (Df /N) ' can serve as an error estimate of the Mersenne twister, but the calculation of the variance for some integrals can be problematic. In addition, the performance of the Monte Carlo method is highly dependent on the choice of a sequence that simulates random numbers, so the standard and actual error can vary greatly.
In general, the value of lg IAMKI lies in the range from -3.5 to 0 and slowly increases with increasing dimension p.
Trapezoidal formula. Its error is determined by the formula
1
12k2
IAnI C _, „ max
d2 fd
dx2
(7)
where k = N1/p is the number of nodes along each coordinate. Thus, its error is O (N-2/p); accuracy should decrease rapidly with increasing . The corresponding curve (black square marker in Fig. 3) illustrates good accuracy lg IAMTI & -5.2 at p = 2; this is much more accurate than the classical Monte Carlo method. However, with increasing p, the error rapidly increases, and already at p > 4 it exceeds the error of the Monte Carlo method. At even higher dimensions, the trapezoidal method quickly becomes uncompetitive.
Sobol sequence. Despite the fact that for high dimensions p the Mersenne twister shows the best result, the shifted Sobol points have a reasonable estimate of the accuracy. It is the standard deviation aa. Thus, even in complex problems, the actual accuracy can be estimated a posteriori using aa, the number of points can be increased and the calculation can be repeated. This is especially important for multidimensional integrals with an unknown exact answer.
8. Conclusion
The magic Sobol points are the most effective for calculating multidimensional integrals. In this paper, an improvement of these sequences is proposed. They are called the shifted Sobol magic points, which provide a more uniform distribution of points in a multidimensional cube. This significantly increases the accuracy of cubatures.
A significant difficulty with Monte Carlo methods is the a posteriori confirmation of the actual accuracy. In this paper, a multigrid algorithm is proposed that allows to find the grid value of the integral simultaneously with a statistically reliable estimate of its accuracy. Previously, such estimates were unknown.
Calculations of representative test integrals with high actual dimension p (up to p = 16) are carried out. Smooth integrands were considered, as well as the multidimensional Weierstrass function having no derivative at any point. These calculations convincingly show the advantages of the proposed methods.
Acknowledgments
This work was supported by grant MK-3630.2021.1.1.
References
[1] I. M. Sobol, Numerical Monte-Carlo methods [Chislennyye metody Monte-Karlo]. Moscow: Nauka, 1973, In Russian.
[2] D. E. Knuth, The art of computer programming, 3rd ed. Reading, Massachusetts: Addison-Wesley, 1997, vol. 2.
[3] G. S. Fishman, Monte Carlo: concepts, algorithms and applications. Berlin: Springer, 1996. DOI: 10.1007/978-1-4757-2553-7.
[4] M. Matsumoto and T. Nishimura, "Mersenne twister: a 623-dimensionally equidistributed uniform pseudo-random number generator," ACM Transactions on Modeling and Computer Simulation (TOMACS), vol. 8, no. 1, pp. 3-30, 1998. DOI: 10.1145/272991.272995.
[5] T. Nishimura, "Tables of 64-bit Mersenne twisters," ACM Transactions on Modeling and Computer Simulation, vol. 10, no. 4, pp. 348-357, 2000. DOI: 10.1145/369534.369540.
[6] "Mersenne Twister Home Page." (2021), [Online]. Available: http:// www.math.sci.hiroshima-u.ac.jp/m-mat/MT/emt.html.
[7] S. K. Park and K. W. Miller, "Random number generators: good ones are hard to find," Communications of the ACM, vol. 31, no. 10, pp. 11921201, 1998. DOI: 10.1145/63039.63042.
[8] M. Mascagni and A. Srinivasan, "Parameterizing parallel multiplicative Lagged-Fibonacci generators," Parallel Computing, vol. 30, pp. 899-916, 2004. DOI: 10.1016/j.parco.2004.06.001.
[9] P. L'Ecuyer, "Good parameter sets for combined multiple recursive random number generators," Operations Research, vol. 47, no. 1, pp. 159164, 1999. DOI: 10.1287/opre.47.1.159.
[10] J. K. Salmon, M. A. Moraes, R. O. Dror, and D. E. Shaw, "Parallel random numbers: as easy as 1, 2, 3," 2011 International Conference for High Performance Computing, Networking, Storage and Analysis (SC), pp. 1-12, 2011. DOI: 10.1145/2063384.2063405.
[11] G. Marsaglia and W. W. Tsang, "The ziggurat method for generating random variables," Journal of Statistical Software, vol. 5, pp. 1-7, 2000. DOI: 10.18637/jss.v005.i08.
[12] G. Marsaglia and A. Zaman, "A new class of random number generators," Annals of Applied Probability, vol. 1, no. 3, pp. 462-480, 1991. DOI: 10.1214/aoap/1177005878.
[13] B. A. Wichmann and I. D. Hill, "An efficient and portable pseudorandom number generator," Applied Statistics, vol. 31, no. 2, pp. 188190, 1982. DOI: 10.2307/2347988.
[14] E. A. Tsvetkov, "Empirical tests for statistical properties of some pseudorandom number generators," Mathematical Models and Computer Simulations, vol. 3, pp. 697-705, 2011. DOI: 10.1134/S20700-4821106010X.
[15] "The Marsaglia Random Number CDROM including the Diehard Battery of Tests of Randomness." (2021), [Online]. Available: http : // ftpmirror.your.org/pub/misc/diehard/.
[16] P. L'Ecuyer and R. Simard, "TestUOl: A C library for empirical testing of random number generators," ACM Transactions on Mathematical Software (TOMS), vol. 33, no. 4, pp. 1-40, 2007. DOI: 10.1145/1268776. 1268777.
[17] L. E. Bassham, A. L. Rukhin, J. Soto, J. R. Nechvatal, M. E. Smid, S. D. Leigh, M. Levenson, M. Vangel, N. A. Heckert, and D. L. Banks, "A statistical test suite for random and pseudorandom number generators for cryptographic applications," National Institute of Standards and Technology, NIST Special Publication, Gaithersburg, MD, Tech. Rep., 2010.
[18] A. A. Belov, N. N. Kalitkin, and M. A. Tintul, "Visual verification of pseudo-random number generators [Vizual'naya verifikatsiya gener-atorov psevdosluchaynykh chisel]," Keldysh IAM Preprints, Moscow, Tech. Rep. 137, 2019, In Russian. DOI: 10.20948/prepr-2019-137.
[19] A. A. Belov, N. N. Kalitkin, and M. A. Tintul, "Unreliability of pseudorandom number generators," Computational Mathematics and Mathematical Physics, vol. 60, no. 11, pp. 1747-1753, 2020. DOI: 10. 1134/S0965542520110044.
[20] I. M. Sobol, "Uniformly distributed sequences with additional uniformity properties," USSR Computational Mathematics and Mathematical Physics, vol. 16, no. 5, pp. 236-242, 1976. DOI: 10.1016/0041-5553(76)90154-3.
For citation:
A. A. Belov, M. A. Tintul, Shifted Sobol points and multigrid Monte Carlo simulation, Discrete and Continuous Models and Applied Computational Science 29 (3) (2021) 260-270. DOI: 10.22363/2658-4670-2021-29-3-260-270.
Information about the authors:
Belov, Aleksandr A. — Candidate of Physical and Mathematical Sciences, Assistant professor of Department of Applied Probability and Informatics of Peoples' Friendship University of Russia (RUDN University); Researcher of Faculty of Physics, M.V. Lomonosov Moscow State University (e-mail: aa.belov@physics.msu.ru, phone: +7(495)9393310, ORCID: https://orcid.org/0000-0002-0918-9263, ResearcherID: Q-5064-2016, Scopus Author ID: 57191950560)
Tintul, Maxim A. — Master's degree student of Faculty of Physics, M. V. Lomonosov Moscow State University (e-mail: maksim.tintul@mail.ru, phone: +7(495)9393310, ORCID: https://orcid.org/0000-0002-5466-1221)
УДК 519.872:519.217
PACS 07.05.Tp, 02.60.Pn, 02.70.Bf
DOI: 10.22363/2658-4670-2021-29-3-260-270
Сдвинутые точки Соболя и многосеточный расчёт методом Монте-Карло
А. А. Белов1'2, М. А. Тинтул1
1 Московский государственный университет им. М. В. Ломоносова Ленинские горы, д. 1, стр. 2, Москва, 119991, Россия 2 Российский университет дружбы народов ул. Миклухо-Маклая, д. 6, Москва, 117198, Россия
Многомерные интегралы возникают во многих задачах физики. Например, моменты функции распределения в задачах переноса различных частиц (фотонов, нейтронов и др.) являются 6-мерными интегралами. При расчёте коэффициентов электропроводности и теплопроводности возникают интегралы рассеяния, размерность которых равна 12. Возникают задачи и с существенно большим числом переменных. Для вычисления интегралов столь высокой кратности наиболее эффективен метод Монте-Карло. Однако работоспособность этого метода сильно зависит от выбора последовательности, имитирующей набор случайных чисел. В литературе описано большое количество генераторов псевдослучайных чисел. Их качество проверяется с помощью батарей формальных тестов. Однако простейший визуальный анализ показывает, что прохождение таких тестов не гарантирует хорошей равномерности этих последовательностей. Для вычисления многомерных интегралов наиболее эффективны магические точки Соболя. В данной работе предложено усовершенствование этих последовательностей — смещённые магические точки Соболя, обеспечивающие большую равномерность распределения точек в многомерном кубе. Это ощутимо повышает точность кубатур. Существенной трудностью методов Монте-Карло является апостериорное подтверждение фактической точности. В данной работе предложен многосеточный алгоритм, позволяющий найти сеточное значение интеграла одновременно со статистически достоверной оценкой его точности. Ранее такие оценки были неизвестны. Проведены расчёты представительных тестовых интегралов с высокой фактической размерностью до 16. В качестве подынтегральной функции выбрана многомерная функция Вейерштрасса, не имеющая производной ни в одной точке. Эти расчёты убедительно показывают преимущества предложенных методов.
Ключевые слова: многомерный интеграл, метод Монте-Карло, точки Соболя, многосеточный расчет, апостериорные оценки точности
