# How can one generate an open-ended sequence of low-discrepancy points in 3D?

I'd like a low-discrepancy sequence of points over a 3D-hypercube $$[-1,1]^3$$, but don't want to have to commit to a fixed number $$n$$ of points beforehand, that is just see how the numerical integration estimates develop with increasing numbers of low-discrepancy points.

I'd like to avoid have to start all over again, if the results with a fixed $$n$$ are unsatisfactory. Of course, one could just employ random numbers, but then the convergence behavior would be poorer.

Since in his answer below, Martin Roberts advances a very interesting, appealing approach to the open-ended low-discrepancy problem, I’d like to indicate an (ongoing) implementation of his approach I’ve just reported in https://arxiv.org/abs/1809.09040 . In sec. XI (p. 19) and Figs. 5 and 6 there, I analyze two problems—one with sampling dimension $$d=36$$ and one with $$d=64$$—both using the parameter $$\bf{\alpha}_0$$ set to 0 and also to $$\frac{1}{2}$$. To convert the quasi-uniformly distributed points yielded by the Roberts’ algorithm to quasi-uniformly distributed normal variates, I use the code developed by Henrik Schumacher in his answer to Can I use Compile to speed up InverseCDF?

• What does "low-discrepancy" mean? – David G. Stork Apr 12 '17 at 20:23
• "A sequence of n-tuples that fills n-space more uniformly than uncorrelated random points, sometimes also called a low-discrepancy sequence. Although the ordinary uniform random numbers and quasirandom sequences both produce uniformly distributed sequences, there is a big difference between the two." (mathworld.wolfram.com/QuasirandomSequence.html) – Paul B. Slater Apr 12 '17 at 20:34
• RandomReal[] supports the Sobol' and Niederreiter sequences via "MKL"; see the docs for details. – J. M. will be back soon Apr 12 '17 at 23:46
• Thanks, J.M. Can you provide a more specific pointer to the docs in question? – Paul B. Slater Apr 13 '17 at 1:52
• See this. (Also, entering "Sobol" or "Niederreiter" in the docs' search box would have returned this part of the docs among the results.) – J. M. will be back soon Apr 14 '17 at 18:50

As the OP cross-posted this question from Math stackexchange, I have also cross-posted the answer that I wrote there.

The simplest traditional solution to the $$d$$-dimensional which provides quite good results in 3-dimensions is to use the Halton sequence based on the first three primes numbers (2,3,5). The Halton sequence is a generalization of the 1-dimensional Van der Corput sequence and merely requires that the three parameters are pairwise-coprime. Further details can be found at the Wikipedia article: "Halton Sequence".

An alternative sequence you could use is the generalization of the Weyl / Kronecker sequence. This sequence also typically uses the first three prime numbers, however, in this case they are chosen merely because the square root of these numbers is irrational.

However, I have recently written a detailed blog post, "The Unreasonable Effectiveness of Quasirandom Sequences", on how to easily create an open-ended low discrepancy sequences in arbitrary dimensions, that is:

• algebraically simpler
• faster to compute
• produces more consistent outputs
• suffers less technical issues

than existing existing low discrepancy sequences, such as the Halton and Kronecker sequences.

The solution is an additive recurrence method (modulo 1) which generalizes the 1-Dimensional problem whose solution depends on the Golden Ratio. The solution to the $$d$$-dimensional problem, depends on a special constant $$\phi_d$$, where $$\phi_d$$ is the value of smallest, positive real-value of $$x$$ such that $$x^{d+1}\;=x+1$$

For $$d=1$$$$\phi_1 = 1.618033989...$$, which is the canonical golden ratio.

For $$d=2$$, $$\phi_2 = 1.3247179572...$$, which  is often called the plastic constant, and has some beautiful properties. This value was conjectured to most likely be the optimal value for a related two-dimensional problem [Hensley, 2002]. Jacob Rus has posted a beautiful visualization of this 2-dimensional low discrepancy sequence, which can be found here.

And finally specifically relating to your question, for $$d=3$$, $$\phi_3 = 1.2207440846...$$

With this special constant in hand, the calculation of the $$n$$-th term is now extremely simple and fast to calculate:

$$R: \mathbf{t}_n = \pmb{\alpha}_0 + n \pmb{\alpha} \; (\textrm{mod} \; 1), \quad n=1,2,3,...$$ $$\textrm{where} \quad \pmb{\alpha} =(\frac{1}{\phi_d}, \frac{1}{\phi_d^2},\frac{1}{\phi_d^3},...\frac{1}{\phi_d^d}),$$

Of course, the reason this is called a recurrence sequence is because the above definition is equivalent to $$R: \mathbf{t}_{n+1} = \mathbf{t}_{n} + \pmb{\alpha} \; (\textrm{mod} \; 1)$$

In nearly all instances, the choice of $$\pmb{\alpha}_0$$ does not change the key characteristics, and so for reasons of obvious simplicity, $$\pmb{\alpha}_0 =\pmb{0}$$ is the usual choice. However, there are some arguments, relating to symmetry, that suggest that $$\pmb{\alpha}_0=\pmb{1/2}$$ is a better choice.

Specifically for $$d=3$$, $$\phi_3 = 1.2207440846...$$ and so for $$\pmb{\alpha}_0= (1/2,1/2,1/2)$$, $$\pmb{\alpha} = (0.819173,0.671044,0.549700)$$ and so the first 5 terms of the canonical 3-dimensional sequence are:

1. (0.319173, 0.171044, 0.0497005)
2. (0.138345, 0.842087, 0.599401)
3. (0.957518, 0.513131, 0.149101)
4. (0.77669, 0.184174, 0.698802)
5. (0.595863, 0.855218, 0.248502) ...

Of course, this sequence ranges between [0,1], and so to convert to a range of [-1,1], simply apply the linear transformation $$x:= 2x+1$$. The result is

1. (-0.361655, -0.657913, -0.900599)
2. (-0.72331, 0.684174, 0.198802)
3. (0.915035, 0.0262616, -0.701797)
4. (0.55338, -0.631651, 0.397604)
5. (0.191725, 0.710436, -0.502995),...

The Mathematica Code for creating this sequence is as follows:

f[n_] := N[Root[-1 - # + #^(n + 1) &, 2 - Boole[EvenQ[n]]]];

d = 3;
n = 5

gamma = 1/f[d];
alpha = Table[gamma^k , {k, Range[d]}]
ptsPhi =  Map[FractionalPart, Table[0.5 + i alpha, {i, Range[n]}], {2}]


Similar Python code is

# Use Newton-Rhapson-Method
def gamma(d):
x=1.0000
for i in range(20):
x = x-(pow(x,d+1)-x-1)/((d+1)*pow(x,d)-1)
return x

d=3
n=5

g = gamma(d)
alpha = np.zeros(d)
for j in range(d):
alpha[j] = pow(1/g,j+1) %1
z = np.zeros((n, d))
for i in range(n):
z = (0.5 + alpha*(i+1)) %1

print(z)


Hope that helps!