This thread is there to enable further discussion on the topic raised in
There are two different points in Owen’s paper:
- First, start the sequence at 0 to preserve the net property
- Second, scramble things if you want to control the error in a quasi-Monte Carlo algorithm
In OT, the first point is the responsibility of the LowDiscrepancySequence class and its derived classes, and the second point is the responsibility of the LowDiscrepancyExperiment class
In my view only the low discrepancy sequences have to be fixed, the LowDiscrepancyExperiment class is correct.
Is scrambling (that is, performing random permutations on digits) implemented in OT? Apparently, since all permutations depend on one another, it is computationally expensive.
There is no scrambling in OT. Pamphile’s suggestion is excellent on the maths, but adding the 0 in e.g. the current Sobol’ sequence does not seem to be a good idea, because it will make most algorithms fail: scrambling is the correct solution for this. It is rather an improvement than a bug.
I adapted a C++ implementation at:
Below is a list of implementations, based on an analysis done for Scilab in 2013:
(Apparently, my Scilab account was hacked in 2019!)
Perhaps there are more up-to-date implementations as of 2020, e.g. Boost or the GSL ?
Regards,
Michaël
Ouppsss sorry for the mess between scrambling and randomization. I may have missed something but I understood that Owen advocated the use of both the initial point at (0,…,0) and scrambling. Am I wrong?
I agree on the fact that it is more an improvement than a bug. Nevertheless, it would be good to allow for the Sobol sequence to start from 0 when needed, which is not currently possible.
You are correct. Having the zero is important as well as randomization. So it is good to scramble, but even better to scramble and have the zero. Having the zero without scrambling performs worse than than scrambling without zero.
From @ArtOwen on Github:
Randomizing by adding a vector \boldsymbol{U} \sim U[0,1]^d and taking the result modulo 1 (ie wraparound) is known as a Cranley-Patterson rotation. They proposed it for lattice rules in the 1970s. It does give you randomized QMC points. But if you use a digital net like a Sobol’ sequence then shifting like this does not produce a new net. It also has an inferior convergence rate increasing variance by a factor of about n for smooth integrands compared to scrambling. Intuitively the problem is that it does not introduce enough randomness.
In the d=1 case, where we don’t need QMC, it is easy to see what is going on. The unrandomized Sobol’ sequence would have x_i = (i-1)/n for 1 \leq i \leq n. Shifting would give points x_i = (i-u)/n for some u in [0,1] after relabeling. Scrambling would give x_i ~ U[(i-1)/n,i/n) independently (i.e., stratified sampling) bringing an error cancellation that shifting does not bring.
Second, scramble things if you want to control the error in a quasi-Monte Carlo algorithm
I think that you may have overlooked one point here. Scrambling a (t, m, s)-net low discrepancy sequence provides a way to estimate the variability in the estimator due to sampling, indeed. But this is only one side: it also improves the convergence rate. See the abstract in [1]:
This paper shows that for smooth integrands over s dimensions, the variance (note by M.B.: of a scrambled (t, m, s)-net) is of order n^{-3}(\log(n)^{s - 1}), compared to n^{-1} for ordinary Monte Carlo. Thus the integration errors are of order n^{-3/2} (\log n)^{(s - 1)/2} in probability. This compares favorably with the rate n^{-1} (\log n)^{s - 1} for unrandomized (t, m, s)-nets.
To see this, please consider the “Sum of 5 exp” example from The first zero point should be moved back into the low discrepancy sequences · Issue #2686 · openturns/openturns · GitHub . This example, from [2], is designed to show the good properties of randomizing and keeping the first 0. We see that the best result is obtain with the scrambled Sobol’ sequence with 0. This is why @tupui said:
So it is good to scramble, but even better to scramble and have the zero.
To get this image, I used the scrambled Sobol’ sequence implemented by @tupui in Scipy. The number in the parentheses in each legend is the slope in the log-log space, estimated using least squares. The picture clearly shows the three different rates of convergence:
- the Monte Carlo rate, close to -0.5,
- the unscrambled Sobol’ sequence rate (from OpenTURNS), close to -1,
- the scrambled Sobol’ sequence with the 0, close -1.5.
If we use the scrambled Sobol’ sequence without the initial 0, then the rate is close to -1. Hence, removing the initial 0 deteriorates the order of convergence by a factor \sqrt{n}. Notice that there is some variability in the previous picture. We can average several repetitions, as done in [2:1], to get a cleaner picture, but this does not change the overall trend.
Notice not only the qualitative, but also the massive quantitative change: the error is lower than 10^{-8} with a sample size close to 10^{6}. For this particular function, we would need a sample size as large as 10^8 using the unscrambled Sobol’ sequence.
Best regards,
Michaël
Owen, A. B. (1997). Scrambled net variance for integrals of smooth functions. The Annals of Statistics, 25(4), 1541-1562. ↩︎
Owen, A. B. (2020, August). On dropping the first Sobol’ point. In International conference on Monte Carlo and quasi-Monte Carlo methods in scientific computing (pp. 71-86). Cham: Springer International Publishing. ↩︎ ↩︎