# 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[-1,1]^3$$, but don’t want to have to commit to a fixed number $$nn$$ 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 $$nn$$ are unsatisfactory. Of course, one could just employ random numbers, but then the convergence behavior would be poorer.

“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)

This question has also just been put on the mathematica.stack.exchange
(https://mathematica.stackexchange.com/questions/143457/how-can-one-generate-an-open-ended-sequence-of-low-discrepancy-points-in-3d)

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=36d=36$$ and one with $$d=64d=64$$—both using the parameter $$α0\bf{\alpha}_0$$ set to 0 and also to $$12\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 https://mathematica.stackexchange.com/questions/181099/can-i-use-compile-to-speed-up-inversecdf

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 $$dd$$-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 $$dd$$-dimensional problem, depends on a special constant $$ϕd\phi_d$$, where $$ϕd\phi_d$$ is the unique positive root of $$xd+1=x+1x^{d+1}=x+1$$.

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

For $$d=2d=2$$, $$ϕ2=1.3247179572… \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=3d=3$$, $$ϕ3=1.2207440846… \phi_3 = 1.2207440846...$$

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

$$R:tn=αα0+nαα(mod1),n=1,2,3,... R: \mathbf{t}_n = \pmb{\alpha}_0 + n \pmb{\alpha} \; (\textrm{mod} \; 1), \quad n=1,2,3,...$$
$$whereαα=(1ϕd,1ϕ2d,1ϕ3d,...1ϕdd), \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:tn+1=tn+αα(mod1) R: \mathbf{t}_{n+1} = \mathbf{t}_{n} + \pmb{\alpha} \; (\textrm{mod} \; 1)$$

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

Specifically for $$d=3d=3$$, $$ϕ3=1.2207440846…\phi_3 = 1.2207440846...$$ and so for $$αα0=(1/2,1/2,1/2)\pmb{\alpha}_0= (1/2,1/2,1/2)$$,
$$αα=(0.819173,0.671044,0.549700)\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 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_] := x /. FindRoot[x^(1 + n) == x + 1, {x, 1}];

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!