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.

“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=36 and one with d=64—both using the parameter α0 set to 0 and also to 12. 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

Answer

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 ϕd, where ϕd is the unique positive root of xd+1=x+1.

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

For d=2, ϕ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, ϕ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:tn=αα0+nαα(mod1),n=1,2,3,...
whereαα=(1ϕd,1ϕ2d,1ϕ3d,...1ϕdd),

Of course, the reason this is called a recurrence sequence is because the above definition is equivalent to
R:tn+1=tn+αα(mod1)

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

Specifically for d=3, ϕ3=1.2207440846… and so for αα0=(1/2,1/2,1/2),
αα=(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:=2x1. 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!

Attribution
Source : Link , Question Author : Paul B. Slater , Answer Author : Martin Roberts

Leave a Comment