# The New Scramble for Faure Sequence Based on Irrational Numbers

Authors: Ali Mogharrabi O.; Behrooz Fathi V.; M. H. Behzadi; R. Farnoosh
Journal: Advances in Mathematical Physics
 (2021)

---

## Abstract

This article intends to review quasirandom sequences, especially the Faure sequence to introduce a new version of scrambled of this sequence based on irrational numbers, as follows to prove the success of this version of the random number sequence generator and use it in future calculations. We introduce this scramble of the Faure sequence and show the performance of this sequence in employed numerical codes to obtain successful test integrals. Here, we define a scrambling matrix so that its elements are irrational numbers. In addition, a new form of radical inverse function has been defined, which by combining it with our new matrix, we will have a sequence that not only has a better close uniform distribution than the previous sequences but also is a more accurate and efficient tool in estimating test integrals.


---

## Body

## 1. Introduction

It is well known that Monte Carlo calculations are based on the generation of random numbers on interval (0,1). Therefore, the generation of random numbers that have more uniformity on (0,1) guarantees better approximations in these calculations. In recent years, some researchers have employed quasirandom sequences instead of random numbers to aim producing extra uniformity of the randomly generated numbers on (0,1). Due to the breadth and complexity of some problems that are mostly unsolvable by classical mathematical methods or solving them with classical methods is associated with more time and computational cost, the stochastic solving of such cases with numerical methods and using the Monte Carlo method plays a key role. The quasirandom sequences are common in Monte Carlo calculations such as Faure, Halton, Niederreiter, and Sobol sequences, but due to the lack of complete success of these sequences in Monte Carlo computation, we use scrambled versions of them, all of which are designed to increase the uniformity of randomly (quasirandom) generated numbers on (0,1), so that we can estimate the obtained solution to the desired unknown solution of the problem.To resolve this problem, researchers are competing on the use of scrambled quasirandom generators based on their version of random number generation to provide more accurate results in Monte Carlo calculations.Today, Monte Carlo and quasi-Monte Carlo methods are widely used to solve the computations of physical and mathematical problems. Quasi-Monte Carlo (QMC) methods play an alternative role for Monte Carlo methods. The advantage of these methods is that they use numbers to provide extera uniformity on unit hypercube. This feature has led to the use of these methods to estimate high-dimensional integrals (Niederreiter, 1992; Spanier and Maize, 1994) [1].So far, several quasirandom sequences (or low discrepancy sequences) have been introduced for the QMC method. Such as the Faure sequence, the Halton sequence, and the Sobol sequence. Despite the fact that among these three sequences, the Faure sequence has better features in terms of discrepancy bound, but in practice, it is less used. Because, the convergence rate of this class of sequences is not so good compared with the other sequences [2]. In addition, due to the correlation between the different dimensions of the Faure sequence, the distribution of the sequence points is not very favorable, and we see poor two-dimensional projections (Figure 1). To overcome this problem, many scrambling methods have been proposed for the Faure sequence (see [3]). In almost all of these scrambles, there are attempts to define a new matrix by shifting numbers in the generating matrix or the placement of the elements, and there has been less talk about the properties of irrational numbers. The matrix that we introduce in this paper has been selected from several proposed matrices. Because it has good two-dimensional projections and it is also at a very high level in terms of integral estimation.Figure 1
512 points from the original Faure sequence in several bases.
(a)(b)(c)(d)In the next section, the structure of the original Faure sequence is given. We then briefly list the scramblers that have already been introduced in Section3. In Section 4, we have brought our proposed matrix. Sections 5 and 6 give the evaluation criteria for the quality of the sequence generated by our proposed matrix and compare it with previous sequences, and in Section 7, the conclusion is stated.


## 2. The Faure Sequence

Supposep≥2 is a prime number, and suppose N is the number of points we want to generate, and n=a0,a1,⋯,am−1T is a vector of integers whose components are from the expansion of the number n on the base p, where 0≤aj≤p and m=logpN. We define the radical inverse function, ϕpn, as
(1)ϕpn=a0p+a1p2+⋯+am−1pm.For the Faure sequence, we define a different generator matrix for each dimension. IfP be the Pascal matrix, then for a s-dimensional Faure sequence the generator matrix of the jth dimension is Cj=Pj−1, 1≤j≤s, where the member on the row c and the column r is defined as follows:
(2)Pj−1=c−1r−1j−1c−rmodp,c≥1,r≥1.Thus, letxnjbe represents the number n in the dimension j in the Faure sequence, then
(3)xnj=ϕpPj−1n,and so the s-dimensional Faure sequence is ϕpP0n,ϕpP1n,⋯,ϕpPs−1n.


## 3. Scrambling the Faure Sequence

Since the introduction of the Faure sequence, several methods were proposed to scramble it. In this section, we give an overview of some of such scrambles.


### 3.1. The Generalized Faure Sequence

Tezuka [2] proposed the generalized Faure sequence, GFaure, with the jth dimension generator matrix Cj=AjPj−1 and the Aj for j=1,2,⋯,s are arbitrary nonsingular lower triangular matrices over Fp. A special case for Aj is that all members are one for all dimensions [4].


### 3.2. Random Linear (Digit) Scrambling

After reviewing different versions of the Owen’s method, Matoušek introduced a scramble matrix and a transfer vector for various dimensions [5]. The sequences obtained by Matoušek have the following general form:
(4)xn=ϕpA1P0n+g1,ϕpA2P1n+g2,⋯,ϕpAsPs−1n+gs.For the random linear scrambling, the matricesAj and the vectors gj for j=1,⋯,s are of the form
(5)Aj=h1,1000h2,1h2,200h3,1h3,2h3,30⋮⋮⋮⋱,gj=g1g2g3⋮,where the gj’s and the hi,j with i≥j are chosen randomly and independently from 0,1,⋯,b1, the hj,j’s are chosen randomly and independently from 1,2,⋯,p−1.Random linear digit method is the basis of other scrambles that will follow. Even the GFaure method is a subset of this method in which the members of the shift vectors are all zero.


### 3.3. I-Binomial Scrambling

A subset of the family of random linear scrambling methods is called left I-binomial scrambling [6]. Here, the Aj is defined as
(6)Aj=h10000h2h1000h3h2h100h4h3h2h10⋱⋱⋱⋱⋱,where h1 is chosen randomly and independently from 1,2,⋯,p−1 and also hi’s i>1 are chosen randomly and independently from 0,1,⋯,p−1.


### 3.4. Striped Matrix Scrambling

The scrambling matrixAj for Striped Matrix Scrambling method, has the following form:
(7)Aj=h10000h1h2000h1h2h300h1h2h3h40⋮⋮⋮⋮⋮,where the hi’s are chosen randomly and independently from 1,2,⋯,p−1. The different types of this scramble are examined in [7], only for problems in the first dimension. Of course, it does not say what changes should be made to the matrix for higher dimensions.


### 3.5. Chi’s Optimal Scramble

When we use the I-binomial method to scramble the Faure sequence, the value of (number)h1 causes all the expansion digits of each number to be replaced (permuted). Now, if we leave out the first digit, the value of (number) h2 causes all the remaining digits to be replaced (permuted).So by cleverly selecting these two members, we can achieve better Faure sequences.In [8], Chi has shown that the best choice for these two values can be obtained based on the primitive roots of the p. Finally, the Chi’s optimal scramble matrix is as follows:
(8)Aj=h1j−1000h2h1j−1000h2h1j−1000h2h1j−1⋮⋮⋮⋱,


### 3.6. Inverse Scrambling

Fathi and Eskandari [9] adapted Chi’s optimal matrix and introduced two kind of matrix Aj as follows:
(9)A1j=dj−10000dj−10000dj−10⋮⋮⋮⋱,A2j=dj−10000dj0000dj+10⋮⋮⋮⋱,where d is chosen from 1,2,⋯,p−1 for j=1,⋯,s.Based on nonlinear congruential method, they proposed another scrambling method for the Faure sequence for which thejth coordinate of the nth point has the general form
(10)xnj=ϕpΦ−1AjΨPj−1n+gj,where Φx and Ψx are bijections that map a digit vector x to another digit vector.


## 3.1. The Generalized Faure Sequence

Tezuka [2] proposed the generalized Faure sequence, GFaure, with the jth dimension generator matrix Cj=AjPj−1 and the Aj for j=1,2,⋯,s are arbitrary nonsingular lower triangular matrices over Fp. A special case for Aj is that all members are one for all dimensions [4].


## 3.2. Random Linear (Digit) Scrambling

After reviewing different versions of the Owen’s method, Matoušek introduced a scramble matrix and a transfer vector for various dimensions [5]. The sequences obtained by Matoušek have the following general form:
(4)xn=ϕpA1P0n+g1,ϕpA2P1n+g2,⋯,ϕpAsPs−1n+gs.For the random linear scrambling, the matricesAj and the vectors gj for j=1,⋯,s are of the form
(5)Aj=h1,1000h2,1h2,200h3,1h3,2h3,30⋮⋮⋮⋱,gj=g1g2g3⋮,where the gj’s and the hi,j with i≥j are chosen randomly and independently from 0,1,⋯,b1, the hj,j’s are chosen randomly and independently from 1,2,⋯,p−1.Random linear digit method is the basis of other scrambles that will follow. Even the GFaure method is a subset of this method in which the members of the shift vectors are all zero.


## 3.3. I-Binomial Scrambling

A subset of the family of random linear scrambling methods is called left I-binomial scrambling [6]. Here, the Aj is defined as
(6)Aj=h10000h2h1000h3h2h100h4h3h2h10⋱⋱⋱⋱⋱,where h1 is chosen randomly and independently from 1,2,⋯,p−1 and also hi’s i>1 are chosen randomly and independently from 0,1,⋯,p−1.


## 3.4. Striped Matrix Scrambling

The scrambling matrixAj for Striped Matrix Scrambling method, has the following form:
(7)Aj=h10000h1h2000h1h2h300h1h2h3h40⋮⋮⋮⋮⋮,where the hi’s are chosen randomly and independently from 1,2,⋯,p−1. The different types of this scramble are examined in [7], only for problems in the first dimension. Of course, it does not say what changes should be made to the matrix for higher dimensions.


## 3.5. Chi’s Optimal Scramble

When we use the I-binomial method to scramble the Faure sequence, the value of (number)h1 causes all the expansion digits of each number to be replaced (permuted). Now, if we leave out the first digit, the value of (number) h2 causes all the remaining digits to be replaced (permuted).So by cleverly selecting these two members, we can achieve better Faure sequences.In [8], Chi has shown that the best choice for these two values can be obtained based on the primitive roots of the p. Finally, the Chi’s optimal scramble matrix is as follows:
(8)Aj=h1j−1000h2h1j−1000h2h1j−1000h2h1j−1⋮⋮⋮⋱,


## 3.6. Inverse Scrambling

Fathi and Eskandari [9] adapted Chi’s optimal matrix and introduced two kind of matrix Aj as follows:
(9)A1j=dj−10000dj−10000dj−10⋮⋮⋮⋱,A2j=dj−10000dj0000dj+10⋮⋮⋮⋱,where d is chosen from 1,2,⋯,p−1 for j=1,⋯,s.Based on nonlinear congruential method, they proposed another scrambling method for the Faure sequence for which thejth coordinate of the nth point has the general form
(10)xnj=ϕpΦ−1AjΨPj−1n+gj,where Φx and Ψx are bijections that map a digit vector x to another digit vector.


## 4. Scrambling Matrix with Irrational Members

In this section, corresponding to the method of random linear digits, we introduce a scrambling matrix that its members are a function of square root of basep (that is, they are irrational numbers).After testing many functions, we found the following function that has the most performance:(11)Aj=p−1000p−1p−100p−1p−1p−10p−1p−1p−1p−1.In (1), we introduced the common form of radical inverse function ϕ for base p. Now, we define a new form of this function as follows:
(12)ϕp′n=a0pm+a1pm−1+⋯+am−1p.We call this asreverse radical inverse function, and we denote the sequences that are made in this way, with the suffix “rev.” Therefore, by combining the matrix Aj and function ϕ′, the general form of number n in dimension j will be as follows:
(13)xnj=ϕ′pAjPj−1n.So, for example, we denote the 40-dimensional Faure sequence generated on the base 41 by the scrambled matrixAj and the function ϕ′ with Aj.41rev.In the following sections, we have examined (studied) the quality of this sequence along with its performance compared to other sequences.


## 5. Investigation of the Uniformity of Generated Sequences

### 5.1. Two-Dimensional Projections

The first step in evaluating the performance of a sequence is to see how the points in the 2D projections are distributed. From Figure2, in these designs, the Faure sequence points are located within parallel lines, which shows that there is a linear correlation between points in successive dimensions. Also in these designs, we see a lot of empty spaces. Therefore, the distribution of good points has not been done. In Figure 3, we draw the same designs using the by Aj.41rev.Figure 2
1000 points from the original Faure sequence.
(a)(b)Figure 3
1000 points from our scrambled Faure sequence.
(a)(b)These figures show that two-dimensional projections ofAj.41rev sequence are better than two-dimensional projections of the original Faure sequence.


### 5.2. Discrepancy

One way to measure the quality of a sequence is to calculate its discrepancy [10]. Warnock shows that
(14)TN2=1N2∑k=1N∑m=1N∏i=1s1−maxxki,xmi−21−sN∑k=1N∏i=1s1−xki2+3−s.where xki is the ith component of the point xk.Figure4 compares TN between the original Faure sequence and the some scrambled Faure sequences that introduced in Section 3. From Figure 4, we see that the discrepancy p diagram for our new sequence is at all points below the other sequences. This is a good indication of the high quality of our sequence points.Figure 4
L2-discrepancy for various randomized scramblings of a 40-dimensional Faure sequence.


## 5.1. Two-Dimensional Projections

The first step in evaluating the performance of a sequence is to see how the points in the 2D projections are distributed. From Figure2, in these designs, the Faure sequence points are located within parallel lines, which shows that there is a linear correlation between points in successive dimensions. Also in these designs, we see a lot of empty spaces. Therefore, the distribution of good points has not been done. In Figure 3, we draw the same designs using the by Aj.41rev.Figure 2
1000 points from the original Faure sequence.
(a)(b)Figure 3
1000 points from our scrambled Faure sequence.
(a)(b)These figures show that two-dimensional projections ofAj.41rev sequence are better than two-dimensional projections of the original Faure sequence.


## 5.2. Discrepancy

One way to measure the quality of a sequence is to calculate its discrepancy [10]. Warnock shows that
(14)TN2=1N2∑k=1N∑m=1N∏i=1s1−maxxki,xmi−21−sN∑k=1N∏i=1s1−xki2+3−s.where xki is the ith component of the point xk.Figure4 compares TN between the original Faure sequence and the some scrambled Faure sequences that introduced in Section 3. From Figure 4, we see that the discrepancy p diagram for our new sequence is at all points below the other sequences. This is a good indication of the high quality of our sequence points.Figure 4
L2-discrepancy for various randomized scramblings of a 40-dimensional Faure sequence.


## 6. Numerical Integration

Another way to compare the quality of sequences is to use them to solve high-dimensional integration problems with numerical methods. Consider the following test integrals:(15)I1f=∫01⋯∫01∏i=1sπ2sinπxidx1⋯dxs=1,(16)I2f=∫01⋯∫01∏i=1s4xi−2+ai1+aidx1⋯dxs=1,where the ai are parameters. There are four choices of parameters as follows:
(1)
ai=0 for 1≤i≤s(2)
ai=1 for 1≤i≤s(3)
ai=i for 1≤i≤s(4)
ai=i2 for 1≤i≤sNote that the most difficult case is whena=0. Because in this case, the importance of all variables is the same and the superposition dimension is approximately the same as the truncation dimension. It is important to know that the larger ai, the less important the variables are, and therefore, the effective dimension becomes smaller. The last three choices of the parameters will be considered here [11].In numerical solution of problems with qMC methods, an accepted procedure is to omit the starting points of the sequence. For example, Fox [12] has suggested that we consider the starting point of the sequence as n=QS4−1. Although, this may lead to better results, note that with this selection, a large amount of points must be omitted. For example, for s=40, we have to skip the initial 2825760 points with this formula, which is practically impossible. We found in our research that if we select the p as the starting point, it will significantly improve the results. We will probably get the best result when we start from the p2+1th point.Therefore, we skip the first 41 points and startn=42 in our calculations. For comparison purposes, we present numerical results for original Faure (Faure), our sequence (Aj-41rev), and two types of scrambled Faure sequences, I-binomial (IB) and random linear digits (RLD).Now, we compare the numerical results of different scrambled Faure sequences presented in this paper. The estimated values for the test functions are given in Figures5–8. These figures show that our proposed scramble has a very acceptable convergence compared to other scrambles.Figure 5
Estimates of the integralI1f by using various Faure sequences.
(a)(b)(c)(d)Figure 6
Estimates of the integralI2f with ai=1 by using various Faure sequences.
(a)(b)(c)(d)Figure 7
Estimates of the integralI2f with ai=i by using various Faure sequences.
(a)(b)(c)(d)Figure 8
Estimates of the integralI2f with ai=i2 by using various Faure sequences.
(a)(b)(c)(d)An observation is that estimated values by the matrix Aj-41rev very close to the actual value. This can be seen in Tables1–4. The estimation error obtained with this scramble in dimension 40, for the function (16) with parameters ai=1, ai=i, and ai=i2 are at most 11.68%, 0.32%, and 0.64%, respectively. However, there are some exceptions. For example, in the first function, when the number of dimensions increases, the accuracy of the estimation decreases. So that, for the dimension 40, the maximum estimated relative error value is 43.7%.Table 1
Estimates ofI1f by using Faure sequences.
GeneratorNs=5s=10s=20s=30s=40Faure5001.21752.230113.4022.64640.0001IB5000.00170.00160.00140.00060.0002RLD5000.00130.00130.00110.00060.0001Aj-41rev5001.03921.04391.52860.80951.1542Faure50000.98031.01283.15850.78480.0314IB50000.03400.03970.04430.04020.0306RLD50000.04690.04750.05480.07370.0341Aj-41rev50001.00151.01541.02191.14010.6935Faure100000.95840.94451.99190.65830.0464IB100000.05690.05640.05930.07740.0496RLD100000.05390.05410.05660.07580.0418Aj-41rev100000.99920.99591.01820.87800.5630Faure200000.99431.05221.49090.54400.0575IB200000.06100.06070.06280.07910.0561RLD200000.05830.05850.05920.06700.0637Aj-41rev200000.99991.00211.06140.93010.7249Faure500000.99651.01270.99720.51590.0674IB500000.06170.06180.06210.06830.0628RLD500000.06040.06050.05820.06570.0587Aj-41rev500000.99840.99831.02410.99010.9173Faure700000.99641.00580.96800.50290.0669IB700000.08090.08060.08770.20002.3467RLD700000.07110.07130.06710.07020.0583Aj-41rev700000.99850.99541.02601.08410.9623Faure1000001.00081.02481.04590.54660.1665IB1000000.46100.46160.46670.46941.7654RLD1000000.27420.27460.27720.25290.1366Aj-41rev1000000.99791.00071.02151.12511.3969Table 2
Estimates ofI2f with ai=1 by using Faure sequences.
GeneratorNs=5s=10s=20s=30s=40Faure5001.07591.23701.79940.69130.1727IB5001.50151.50001.55721.34451.1510RLD5001.51231.48201.54561.46611.2191Aj-41rev5000.98790.90140.93470.89300.9629Faure50000.99630.99081.03420.74060.4475IB50001.48361.47851.47871.42121.2841RLD50001.48071.47701.48401.44551.2824Aj-41rev50000.99760.98400.99640.98981.0902Faure100000.99350.98641.00620.78600.5359IB100001.47691.47501.46971.43701.3296RLD100001.47791.47591.47841.44381.3261Aj-41rev100000.99900.98650.99250.99681.0920Faure200000.99450.99540.98610.81400.5897IB200001.47531.47551.47061.42441.3306RLD200001.47631.47521.46401.42881.3271Aj-41rev200001.00001.00021.01511.03361.0965Faure500000.99920.99690.96090.83970.6360IB500001.47501.47441.46971.43971.3342RLD500001.47551.47511.47251.43981.3372Aj-41rev500001.00151.00391.01921.02221.0509Faure700000.99830.99350.96130.84280.6458IB700001.46541.46501.45821.42491.3238RLD700001.47101.47051.46531.43621.3491Aj-41rev700001.00201.00401.02371.04181.1168Faure1000000.99971.00030.97980.88560.7312IB1000001.27851.27831.27271.24761.1649RLD1000001.38481.38441.38211.35081.2698Aj-41rev1000001.00271.00501.02331.03711.1048Table 3
Estimates ofI2f with ai=i by using Faure sequences.
GeneratorNs=5s=10s=20s=30s=40Faure5001.02431.02001.00770.99450.9817IB5001.50101.50221.50311.50301.5023RLD5001.50551.50291.50431.50511.5045Aj-41rev5000.99510.97460.97680.97610.9770Faure50000.99720.99440.99080.98770.9851IB50001.48371.48351.48331.48331.4832RLD50001.48101.48091.48071.48061.4804Aj-41rev50000.99900.99570.99630.99680.9974Faure100000.99770.99640.99460.99290.9915IB100001.47701.47691.47681.47681.4766RLD100001.47811.47801.47791.47781.4777Aj-41rev100000.99980.99710.99730.99770.9984Faure200000.99790.99770.99650.99560.9949IB200001.47521.47531.47531.47521.4752RLD200001.47641.47631.47631.47621.4762Aj-41rev200001.00021.00021.00071.00121.0017Faure500000.99990.99960.99910.99870.9983IB500001.47501.47501.47501.47501.4749RLD500001.47551.47551.47551.47551.4755Aj-41rev500001.00091.00121.00191.00241.0027Faure700000.99930.99890.99860.99830.9980IB700001.46541.46541.46541.46541.4653RLD700001.47101.47101.47091.47091.4709Aj-41rev700001.00111.00131.00211.00261.0029Faure1000000.99970.99960.99930.99900.9989IB1000001.27851.27851.27851.27851.2784RLD1000001.38481.38481.38481.38481.3848Aj-41rev1000001.00131.00161.00241.00281.0032Table 4
Estimates ofI2f with ai=i2 by using Faure sequences.
GeneratorNs=5s=10s=20s=30s=40Faure5001.01251.01151.01101.01081.0106IB5001.50031.50051.50061.50061.5006RLD5001.50171.50151.50151.50151.5015Aj-41rev5000.99660.99360.99360.99360.9936Faure50000.99880.99860.99850.99840.9984IB50001.48381.48381.48381.48381.4838RLD50001.48101.48101.48101.48101.4810Aj-41rev50000.99930.99890.99890.99890.9989Faure100000.99910.99900.99890.99890.9988IB100001.47701.47701.47701.47701.4770RLD100001.47811.47811.47811.47811.4781Aj-41rev100001.00000.99960.99960.99960.9996Faure200000.99920.99910.99910.99910.9991IB200001.47531.47531.47531.47531.4753RLD200001.47641.47641.47641.47641.4764Aj-41rev200001.00031.00021.00031.00031.0003Faure500001.00001.00001.00001.00001.0000IB500001.47501.47501.47501.47501.4750RLD500001.47551.47551.47551.47551.4755Aj-41rev500001.00061.00061.00071.00071.0007Faure700000.99980.99970.99970.99970.9997IB700001.46551.46551.46551.46551.4655RLD700001.47101.47101.47101.47101.4710Aj-41rev700001.00061.00061.00071.00071.0007Faure1000000.99980.99980.99980.99980.9998IB1000001.27851.27851.27851.27851.2785RLD1000001.38481.38481.38481.38481.3848Aj-41rev1000001.00071.00071.00071.00081.0008


## 7. Conclusion

We studied the original Faure sequence and some of its recent years introduced scrambles. Then, we introduced a new scrambling matrix based on irrational numbers that its elements are function of square root of basep. In Figures 5–8, we have shown that this modified scrambled Faure sequence provides better results than the previous versions of its scrambles. Also, we presented that this modified scrambled Faure sequence has greatly improved the distribution of points. The 2D designs confirm this claim, good two-dimensional projections in successive dimensions. As mentioned in the previous section, using this scramble leads to very small estimated relative errors that can often be ignored. In our next research, we will use deterministic scrambling matrices based on another irrational numbers and primitive roots. The results proved the improvements of accuracy using our new scramble.



---

*Source: 6696895-2021-04-28.xml*