2
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ What does "low-discrepancy" mean? $\endgroup$ – David G. Stork Apr 12 '17 at 20:23
  • $\begingroup$ "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) $\endgroup$ – Paul B. Slater Apr 12 '17 at 20:34
  • $\begingroup$ RandomReal[] supports the Sobol' and Niederreiter sequences via "MKL"; see the docs for details. $\endgroup$ – J. M. will be back soon Apr 12 '17 at 23:46
  • $\begingroup$ Thanks, J.M. Can you provide a more specific pointer to the docs in question? $\endgroup$ – Paul B. Slater Apr 13 '17 at 1:52
  • $\begingroup$ See this. (Also, entering "Sobol" or "Niederreiter" in the docs' search box would have returned this part of the docs among the results.) $\endgroup$ – J. M. will be back soon Apr 14 '17 at 18:50
2
$\begingroup$

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!

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.