Proceedings of the Tomsk Polytechnic Institute named after S. M. Kirov, Order of the October Revolution and Order of the Red Banner of Labor

Vol. 277　　　　　　　　　　　　　　　　　　　　　　　　　1977

ON THE PROSPECTS FOR THE PRACTICAL APPLICATION OF I. M. SOBOL PSEUDORANDOM SEQUENCES

V. V. ZAKHAROV, V. G. SMIRNOVA

(Submitted by the scientific seminar of the Laboratory of Computing Technology and Automation, Research Institute of Nuclear Physics, Electronics, and Automation)

In recent years, methods for computing multidimensional integrals based on the use of pseudorandom sequences to form integration nodes have found wide application [1—3]. In this connection, related areas of application of pseudorandom sequences may be noted.

— In most problems of modeling complex objects, criterion relations can also be represented in the form of certain multidimensional integrals [2], although not always in explicit form.

— In search methods of the type [4, 5], based on preliminary scattering of trials of the function being optimized, the search time depends substantially on the uniformity of the scattering.

— In problems of experimental design [6], the optimal strategy consists in filling the space of variables with trials as uniformly as possible.

All this compels special attention to be paid to the quality of generation of pseudorandom sequences. In the present work an attempt is made at a numerical comparison of certain deterministic sequences [1, 2, 3] according to the criterion of uniformity of the generated pseudorandom variables. Establishing the effectiveness of these methods theoretically is difficult, since the corresponding relations are asymptotic in character.

On the BESM-4 computer, Sobol sequences [3] (LP\(_\tau\)-sequences) were experimentally compared with sequences often used in practice, specified by standard random-number generators (DSCh-sequences, generator No. 4 from [7]), and also with Korobov sequences (K-sequences, [2]). The effectiveness of the sequences is established by comparing the accuracy of computing multidimensional integrals whose integration nodes were constructed from the corresponding sequences.

Integrals of the following four functions of increasing complexity were computed:

\[ f_1=\sum_{i=1}^{n} a_i x_i,\qquad f_2=\sum_{i=1}^{n} b_i x_i^2, \]

\[ f_3=\sum_{i=1}^{n}\left(a_i x_i^3+b_i x_i^2+c_i x_i\right), \]

\[ f_4=\sum_{i,j=1}^{n} a_{ij}x_i x_j, \]

whose coefficients in each problem were chosen with equal probability from \((0,1)\). Thus, the experiments were carried out on fairly smooth functions.

In the main experiment (Tables 1–4), a criterion was chosen that made it possible to find the minimal value of the quantity \(N\) that ensured a prescribed value of the error \(\varepsilon\). The choice of this somewhat unusual criterion was due to the need for an infinite sequence of nodes independent of \(N\), and used in schemes for searching for the minimum of functions of many variables [4, 5].

For each of the functions given, 5 problems were solved, differing from one another in the set of coefficients. The dimensionality \((n)\) of the computed integrals was equal to 3, 5, 7, 9, 11, 13. The tables give the smallest values of \(N\), averaged over 5 problems, for which the prescribed accuracy \(\varepsilon\) of the computation of the integrals was satisfied.

Table 1

\(f_1\)

Table 2

\(f_2\)

Table 3

\(f_3\)

Table 4

The fact that Sobol sequences are more efficient in the computation of integrals is not among the unexpected ones. Theoretically, these sequences provide convergence of order $1/N$, whereas DS sequences provide only $1/\sqrt{N}$. From the tables presented, however, it follows that this asymptotics begins to hold only for $n \gg 10$ and $\varepsilon \leq 10^{-4}$. In the region of values $n < 10$, $\varepsilon > 10^{-4}$, both types of sequences are in fact more efficient than could have been expected from their asymptotic properties. Moreover, LP$_\tau$ sequences are practically everywhere in the range of values of $n$, $N$, and $\varepsilon$ considered more efficient than the asymptotic forecast.

In experimenting with Sobol sequences, some anomalous results were also obtained. Thus, in Table 3, for $n=3$ the accuracy $\varepsilon = 10^{-4}$ was achieved, on average, only after 13 thousand trials in each of the 5 problems. Thus, computing the integral of the function

\[ f_3(x)=0.5177\,x_1^3+0.9313\,x_2^3+0.3243\,x_3^3+0.7804\,x_1^2+0.661\,x_3^2+ \]

\[ +0.1328\,x_2^2+0.6661\,x_1+0.5141\,x_2+0.2726\,x_3 \]

with error $\varepsilon=10^{-4}$ required $N>15000$, whereas for the function

\[ f_3(x)=0.3335\,x_1^3+0.0644\,x_2^3+0.2619\,x_3^3+0.0452\,x_1^2+ \]

\[ +0.4834\,x_2^2+0.5603\,x_3^2+0.3452\,x_1+0.7152\,x_2+0.2776\,x_3 \]

$N=488$ proved sufficient. Apparently, this is explained by the strong scatter in the asymptotic behavior of LP sequences.

At the same time, in general the asymptotic laws do hold. This is seen from Fig. 1, where the asymptotic functions

Fig. 1. Asymptotic properties of DS (curve $E_{1/2}$) and LP$_\tau$ (curves $E_1'$ and $E_1''$) sequences.

\[ E_{1/2}(N)=\frac{\varepsilon(N)}{\sqrt{N}},\qquad E_1(N)=\frac{\varepsilon(N)}{N},\qquad N=0(32)\ 16000. \]

for the linear function \(f_1(x)\). Under an ideal approximation of the asymptotics, the curves should approach certain straight lines \(E_{1/2}(N)=\mathrm{const}\), \(E_1=\mathrm{const}\). In fact, the asymptotics of PRN sequences (\(E_{1/2}\), the solid line in Fig. 1) reproduces complex low-frequency oscillations, while the values of the asymptotic function of Sobol sequences (\(E_1\)) have a very large scatter and are contained mainly between the dotted lines in Fig. 1. About 45% of all values of \(E_1(N)\) are concentrated in the hatched region. Thus, if the central line of this region is taken approximately as the asymptotic constant \(C_1\), it turns out that \(C_1\approx 0.8\). For PRN sequences the most probable asymptotic constant is \(C_{1/2}\approx 0.2\).

Fig. 2. Dependence of the accuracy of integration on the length of the LP\(_\tau\)-sequence

From Fig. 1 one can also see that the quality of PRN sequences is very nonuniform in different broad ranges of values of \(N\). Thus, for \(6000<N<1200\) the asymptotics \(E_{1/2}\) is well satisfied with a very small value of the asymptotic constant \(C_{1/2}\). For values \(N>12000\), the convergence due to PRN sequences is clearly lower than \(N^{-1/2}\). Such behavior makes it important to pose the problem of finding the most uniform segments of PRN sequences by means of a statistical analysis of the corresponding random-number generator.

The strong scatter of the asymptotics of LP-sequences, as well as the analysis of Tables 1–4, compels one to suppose the existence of such values of \(N\) for which the error in computing multidimensional integrals is substantially lower than the mean-asymptotic one. Computational experiments showed that such values \(N=N_c\) do indeed exist and form the sequence:

\[ N_c=32,\ 64,\ 128,\ 248,\ 256,\ 400,\ 448,\ 488,\ 648,\ 992,\ 1000,\ 1176,\ 1408,\ 1504, \]

1792, 2048, 2272, 3008, 3488, 4096, 5120, 5568, 6208, 6400, 6720, 6848,
7040, 7360, 8192, 8356.

As can be seen from Fig. 2, the accuracy of the computation of integrals for these values of \(N\) is substantially higher than in a certain very small neighborhood of them, and the asymptotics has order \(1/N^{1+\alpha}\), \(\alpha>0\).

It is very important to emphasize that the subsequence \(N_c\) retains its best properties, in the sense of uniformity, for all investigated \(n\) and for any of the functions \(f_1—f_4\). Thus, in practical computations it makes sense to choose \(N=N_c\), since it is natural to suppose that the selected LP\(_\tau\)-sequence retains the properties of increased uniformity for the entire class of continuous sufficiently smooth functions.

Fig. 3. Dependence of the accuracy of integration on the length of the best subsequences of PRN and LP\(_\tau\).

For PRN sequences, a certain improvement in uniformity was also observed for particular \(N=N_D\). However, these \(N_D\) are different for each function, which makes practical application of these properties of PRN sequences impossible.

In Fig. 3 are shown graphs of the dependence \(\varepsilon(N)\) when integrating a function of the form \(f_3\) (\(n=5\)) for LP\(_\tau\)-sequences (dash-dot line), for PRN sequences (dotted line), and for Korobov grids (solid line). To construct the curves, the best subsequence lengths (\(N_D\) and \(N_c\)) were used. In other words, the curves for PRN and LP\(_\tau\) in Fig. 3 are the lower envelopes of broken lines similar to the one shown in Fig. 2.

It is evident from Fig. 3 that Korobov’s method permits computations to be performed with the smallest error. However, in high dimensions the smallest value of \(N\) that guarantees high accuracy is quite large. For example, for \(n=10\) the minimum value is \(N=85633\). In some practical simulation problems the value of \(N\) must be chosen smaller even by an order of magnitude. Thus, the Korobov grid does not allow one to control continuously the accuracy of compu-

...calculations as a function of the value of \(N\). The need to use arbitrary values of \(N\) (even if with lower computational accuracy) also arises when using the indicated search methods [4], [5]. In this case, the resulting error is smoothed out by the search method itself. However, arbitrary values of \(N\) cannot be chosen on the basis of Korobov’s method. The method for generating pseudorandom sequences proposed by I. M. Sobol is free of these two substantial shortcomings.

REFERENCES

1. N. P. Buslenko, D. I. Golenko, I. M. Sobol, Yu. A. Shreider. The Monte Carlo Method. Moscow, Fizmatgiz, 1962.

2. N. M. Korobov. Number-Theoretic Methods in Approximate Analysis. Moscow, Fizmatgiz, 1963.

3. I. M. Sobol. Multidimensional Quadrature Formulas and Haar Functions. Moscow, “Nauka,” 1969.

4. I. B. Motsius. Multiextremal Problems in Design. Moscow, “Nauka,” 1967.

5. V. V. Zakharov. One method of random search. Izv. AN SSSR, “Technical Cybernetics,” No. 4, 1969.

6. V. V. Nalimov, N. A. Chernova. Statistical Methods for Planning Extremal Experiments. Moscow, “Nauka,” 1965.

7. V. F. Lyashchenko. Programming for M-20 Digital Computers. Moscow, “Sovetskoe radio,” 1967.