I present a new low discrepancy quasirandom sequence that offers substantial improvement over current state-of-the art
sequences eg Niederreiter, Sobol,…

Figure 1. Comparison of the various low discrepancy quasirandom sequences. Note that the newly proposed $R$-sequence produces more evenly spaced points than any of the other methods. Furthermore, all other current methods require careful selection of basis parameters, and if not chosen carefully can lead to degeneracy (eg top right).

############################################  

Hot off the press.

This post is currently trending on front page of HackerNews.

Please upvote it! 😉

https://news.ycombinator.com

#############################################                 

Topics Covered

  • Low discrepancy sequences in one-dimension
  • Low discrepancy methods in two dimensions
  • Multiscale low discrepancy
  • Quasirandom sequences on the surface of sphere
  • Aperiodic Tiling of the plane
  • Low discrepancy sequences in computer graphics

Introduction: Random versus Quasirandom

In figure 1, it can be seen that the simple uniform random sampling of a point inside a unit square exhibits clustering of points and that there are also regions that contain no points at all (‘white noise’) . A low discrepancy sequence  quasirandom sequence is a method of constructing an (infinite) sequence points in a deterministic manner that reduces the likelihood of clustering (discrepancy) whilst still ensuring that the entire space is uniformly covered (‘blue noise’).

Quasirandom Sequences in 1-dimension

The methods of creating fully deterministic low discrepancy quasirandom sequences in one dimension are extremely well studied, and generally solved.  In this post, I focus almost solely on open (infinite) sequences, first in one dimension and then extending to higher dimensions. The fundamental advantage of open sequences (that is, extensible in $n$ is that if the resultant errors based on a finite number of terms is too large, the sequence can be extended without discarding all the previously calculated points. There are many methods to construct open sequences. One way to categorise the different types is by the method of constructing their basis (hyper)-parameters:

  • irrational fractions: Kronecker, Richtmyer, Ramshaw
  • (co)prime numbers: Van der Corput, Halton, Faure
  • Irreducible Polynomials : Niederreiter
  • Primitive polynomials: Sobol’

For brevity, this post generally focuses on comparing this new additive recurrence $R$-sequence, which falls into the first category as it is a recurrence method based on irrational numbers (often called Kronecker or Weyl  sequences) which are rank 1 lattices, and the Halton sequence, which is based on the canonical one-dimensional van der Corput sequence.

The canonical Kronecker Recurrence sequence is defined as:

$$R_1(\alpha): \;\; t_n = [s_0 + n \alpha], \quad n=1,2,3,…$$

where $\alpha $ is any irrational number. Note that the notation $ [x]$ indicates the fractional part of $x$. In computing, this function is more commonly expressed in the following way

$$R_1(\alpha): \;\; t_n = s_0 + n \alpha \; (\textrm{mod} \; 1); \quad n=1,2,3,…$$

For $s_0 = 0$, the first few terms of the sequence, $R(\phi)$ are:

$$t_n = 0.618, \; 0.236, \; 0.854, \; 0.472, \; 0.090, \; 0.708, \; 0.327, \; 0.944, \; 0.562, \; 0.180,\; 798,\; 416, \; 0.034, \; 0.652, \; 0.271, \; 0.888,…$$

It is important to note that the value of $s_0$ does not affect the overall characteristics of the sequence, and in nearly all cases is set to zero. However, especially in computing the option of $s \neq 0$ offers an additional degree of freedom that is often useful. If $s \neq 0$, the sequence is often called a ‘shifted lattice sequence’.

The value of \( \alpha \) that gives the lowest possible discrepancy is achieved if  \( \alpha = 1/\phi \), where \( \phi \) is the golden ratio. That is,

$$ \phi \equiv \frac{\sqrt{5}+1}{2} \simeq 1.61803398875… ; $$

It is interesting to note that there are an infinite number of other values of  $\alpha$ that also achieve optimal discrepancy, and they are all related via the Moebius transformation

$$ \alpha’ = \frac{p\alpha+q}{r\alpha+s} \quad \textrm{for all integers} \; p,q,r,s \quad \textrm{such that} |ps-qr|=1 $$

We now compare this recurrence method to the well known van der Corput reverse-radix sequences [van der Corput, 1935] . The van der Corput sequences are actually a family of sequences, each defined by a unique hyper-parameter, b. The first few terms of the sequence for b=2 are:

$$t_n^{[2]} =  \frac{1}{2}, \frac{1}{4},\frac{3}{4}, \frac{1}{8}, \frac{5}{8},  \frac{3}{8}, \frac{7}{8}, \frac{1}{16},\frac{9}{16},\frac{5}{16},\frac{13}{16}, \frac{3}{16}, \frac{11}{16}, \frac{7}{16}, \frac{15}{16},…$$

The following section compares the general characteristics and effectiveness of each of these sequences.

Consider the task of evaluating the definite integral

$$ A = \int_0^1 f(x) \textrm{d}x $$

We may approximate this by:

$$ A  \simeq A_n = \frac{1}{n} \sum_{i=1}^{n} f(x_i),  \quad x_i \in [0,1] $$

  • If the \( \{x_i\} \) are equal to  \(i/n\), this is the rectangle rule;
  • If the \( \{x_i\} \) are chosen randomly, this is the Monte Carlo method; and
  • If the \( \{x_i\} \) are elements of a low discrepancy sequence, this is the quasi-Monte Carlo method.

The following graph shows the typical error curves \( s_n = |A-A_n| \) for approximating a definite integral associated with the function, \( f(x) = \textrm{exp}(\frac{-x^2}{2}), \; x \in [0,1]  \) with: (i) quasi-random points based on the additive recurrence, where  \( \alpha = 1/\phi \), (blue); (ii) quasi-random points based on the van der Corput sequence, (orange); (iii) randomly selected points, (green); (iv) Sobol sequence (red).

 It shows that for \(n=10^6\) points, the random sampling approach results in an error of \( \simeq 10^{-4}\), the van der Corput sequence results in an error of \( \simeq 10^{-6}\), whilst the \(R(\phi)\)-sequence results in an error of \( \simeq 10^{-7}\), which is \(\sim\)10x better than the van der Corput error and \(\sim\) 1000x better than (uniform) random sampling.

Figure 2. Comparison of 1-dimensional Numerical integration using various Quasirandom Monte Carlo methods. Note smaller is better.

Several things from this figure are worth noting:

  • it is consistent with the knowledge that the errors based on uniform random sampling asymptotically decrease by $ 1/\sqrt{n} $, whereas error curves based on both quasi-random sequences tend to  $1/n $.
  • The results for the $R_1(\phi)$-sequence  (blue) and Sobol (red) are very the best.
  • It shows that the van der Corput sequence offers good, but incredibly consistent results for integration problems!
  • It shows that, although it has larger variance, the  $R_1(\phi)$-sequence produces even lower errors than the van der Corput sequence, with the worst errors being better than the typical errors of the van der Corput sequence.

The new $R_1$ sequence, which is the Kronecker sequence using the Golden ratio,  is one of the best choices for one-dimensional Quasirandom Monte Carlo (QMC) integration methods

It should also be noted that although \(\alpha = \phi\) theoretically offers the provably optimal option, \(\sqrt{2}\) and is very close to optimal, and almost any other irrational value of \(\alpha\) provides excellent error curves for one-dimensional integration. This is why, \(\alpha = \sqrt{p} \) for any prime is very commonly used.  Furthermore, from a computing perspective, selecting a random value in the interval  $ \alpha \in [0,1] $ is with almost certainty going to be (within machine precision) an irrational number, and therefore a solid choice for a low discrepancy sequence.

For visual clarity, the above figure does not show the results the Niederreiter sequence as the results  are virtually indistinguishable to that of the the Sobol and $R$ sequences.  The Neiderreiter and Sobol sequences (along with their optimized parameter selection) that were used in this post were calculated via Mathematica using what is documented as “closed, proprietary and fully optimized generators provided in Intel’s MKL library”.

Quasirandom sequences in two dimensions

Most current methods to construct higher dimension low discrepancy simply combine (in a component-wise manner), \(d\) one-dimensional sequences together. For brevity, this post mainly focuses on describing the Halton sequence [Halton, 1960] and the \(d-\)dimensional Kronecker sequence.

The Halton sequence is constructed simply by using \(d\) different one-dimensional van de Corput sequence each with a base that is relatively prime to all the others. That is, pairwise co-prime. By far, the most frequent selection, due to its obvious simplicity and sensibility, is to select the first \(d\) primes. The distribution of the first 625 points defined by the (2,3)-Halton sequence is show in figure 1.

Although many two-dimensional Halton sequences are excellent sources for low discrepancy sequences, it is also well known that many are highly problematic and do not exhibit low discrepancies. For example, figure 3 shows that the (11,13)-Halton sequence produces highly visible lines. Much effort has gone into methods of selecting which pairs of ( \(p_1, p_2) \) are exemplars and which ones are problematic. This issue is even more problematic in higher dimensions.

Recurrence methods generally suffer even greater challenges when generalising to higher dimensions. That is, although using \( \alpha = \sqrt{p}  \) produces excellent one-dimensional sequences, it is very challenging to even find pairs of prime numbers  to be used as the basis for the two dimensional case that are not problematic!

As a way around this, some have suggested using other well-known irrational numbers, such as the \( \phi,\pi,e, … \). These produce moderately acceptable solutions but not are generally not used as they are usually not as good as a well-chosen Halton sequence.

A great deal of effort has focused on these issues of degeneracy. Proposed solutions include skipping/burning, leaping/thinning. And for finite sequences scrambling is another technique that is frequently used to overcome this problem. Scrambling can not be used to create an open (infinite) low discrepancy sequence.

Figure 3. The (11,13)-Halton sequence is clearly not a low discrepancy sequence (left). Nor is the (11,13)-prime based additive recurrence sequence (middle). Some 2-dimensional additive recurrence sequences that incorporate well known irrational numbers are reasonably good (right).

Similarly, despite the generally better performance of the Sobol sequence, its complexity and more importantly the requirement of very careful choices of its hyperparameters makes it not as inviting.

  • Thus reiterating, in $d$-dimensions:
  • the typical Kronecker sequences require the selection of $d$ linearly independent irrational numbers;
  • the Halton sequence requires $d$ pairwise coprime integers; and
  • the Sobol sequence requires selecting $d$ direction numbers.


The new $R_d$ sequence is the only $d$-dimensional low discrepancy quasirandom sequence that does not require any selection of basis parameters.

Generalizing the Golden Ratio

tl;dr In this section, I show how to construct a  new class of \(d-\)dimensional open (infinite) low discrepancy sequence that do not require choosing any basis parameters, and which has outstanding low discrepancy properties.

There are many possible ways to generalize the Fibonacci sequence and/or the Golden ratio. The following proposed method of generalizing the golden ratio is not new [Krcadinac, 2005]. Also the characteristic polynomial is related to many fields of algebra, including Perron Numbers, and Pisot-Vijayaraghavan numbers. However, what is new is the explicit connection between this generalized form and the construction of higher-dimensional low discrepancy sequences.

We define the generalized version of the golden ratio,\( \phi_d\) as the unique positive root $ x^{d+1}=x+1 $. That is,

For \(d=1\),  \( \phi_1 = 1.61803398874989484820458683436563… \), which is the canonical golden ratio.

For \(d=2\), \( \phi_2 = 1.32471795724474602596090885447809… \), 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].

For \(d=3\), \( \phi_3 = 1.220744084605759475361685349108831… \),

For $ d>3$, although the roots of this equation do not have a closed algebraic form, we can easily obtain a numerical approximation either through standard means such as Newton-Rhapson, or by noting that for the following sequence, \( R_d(\phi_d) \):

$$ t_0=t_1 = … = t_{d} = 1;$$

$$t_{n+d+1} \;=\;t_{n+1}+t_n, \quad \textrm{for} \; n=1,2,3,..$$

We have the following very elegant property:

$$ \phi_d  =\lim_{n\to \infty} \;\;\frac{t_{n+1}}{t_n} .$$

This sequence, sometimes called the generalized or delayed fibonacci sequence, has been studied quite widely, [Kak 2004, Wilson 1993] and the sequence for \(d=2\) is often called the Padovan sequence [Stewart, 1996, OEIS A000931], whilst the \(d=3\) sequence is listed in [OEIS A079398].

As mentioned before, the key contribution of this post is to describe the explicit connection between this generalized sequence and the construction of \(d-\)dimensional low-discrepancy sequences.

Major Result: The following parameter-free \(d-\)dimensional open (infinite) sequence \(R_d(\phi_d)\), has superior low discrepancy characteristics when compared to other existing methods.

$$ \mathbf{t}_n = [n \pmb{\alpha}],  \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}), $$

$$ \textrm{and} \; \phi_d\ \textrm{is the unique positive root of }  x^{d+1}=x+1. $$

For two dimensions, this generalized sequence for \(n=150\), is shown in figure 1. The points are clearly more evenly distributed  for the \(R_2\)-sequence compared to the (2, 3)-Halton sequence, the Kronecker sequence based on \( (\sqrt{3},\sqrt{7}) \), the Niederreiter and Sobol sequences. (Due to the complexity of the Niederreiter and Sobol sequences they were calculated via Mathematica using proprietary code supplied by Intel.)

This type of sequence, where the basis vector $ \pmb{\alpha} $ is a function of a single real value, is often called a Korobov sequence [Korobov 1959]

However, note that despite this similarity, the construction of the $R_2$-sequence is not a lattice, as the first component of  $ \pmb{\alpha} $ is not unity (or a rational number), but rather an irrational number. Furthermore, in contrast to lattices, the $R_d(\phi_d)$ sequence is intrinsically extensible in $n$.

See figure 1 again for a comparison between various 2-dimensional low discrepancy quasirandom sequences.

Code and Demonstrations

Note that Python arrays and loops start at zero!

Template python code.


import numpy as np # Use Newton-Raphson-Method to calculate g=phi_d # or you could just hard-code it. # gamma(1) = 1.61803398874989484820458683436563 # gamma(2) = 1.32471795724474602596090885447809 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

# Number of dimensions. d=2 # number of required points n=50 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)) # this number can be any real number. # Default setting typically 0 or 0.5 seed = 0.5 for i in range(n): z[i] = (seed + alpha*(i+1)) %1 print(z)

I have written the code like this to be consistent with the maths conventions used throughout this post. However, for reasons of programming convention and/or efficiency there are two possible modifications worth considering.

Firstly, as $R_2$ is an additive recurrence sequence, an alternative formulation of $z$ which does not require floating point multiplication is

 z[i+1] = (z[i] +alpha ) %1 

Secondly, as for many modern languages that allow vectorization, much of this code coudl be vectorised. For example.
for i in range(n): z[i] = seed + alpha*(i+1)

z = z%1

Here are some demonstrations with included code, by other people based on this sequence:

Minimum Distance


The new $R_2$-sequence  is the only 2-dimensional low discrepancy quasirandom sequence where the minimum packing distance falls only by $1/\sqrt{n}$.

Although the standard technical analysis of quantifying discrepancy is by analysing the \(d^*\)-discrepancy, we first mention a couple of other geometric (and possibly more intuitive ways!) of how this new sequence is preferable to other standard methods.

If we denote the distance between points \(i\) and \(j\) to be \(d_{ij}\), and \(d_0 = \textrm{inf} \; d_{ij} \), then the following graph shows how \(d_0(n)\) varies for the \(R\)-sequence, the (2,3)-Halton, Sobol, Niederrieter and the random sequences.  This can be seen in the following figure 6.

Like the previous figure, each minimum distance measure is normalized by a factor of \( 1/\sqrt{n} \).

One can see that after \(n=300\) points, it is almost certain that for the random distribution (green) there will be two points that are extremely close to each other. It can also be seen that although the (2,3)-Halton sequence (orange) is much better than the random sampling, it also unfortunately asymptotically goes to zero.

The reason why the normalized $d_0$ goes to zero for the Sobol sequence is because Sobol himself showed that the Sobol sequence falls at a rate of $1/n$ — which is good, but obviously much worse than for $R_2$ which only falls by $1/\sqrt{n}$.

The \(R(\phi_2)\) sequence, (blue) even for all values of \( n \< 10^6\) points, the minimum distance between two points consistently falls between \(0.549/\sqrt{n} \) and \(0.868/\sqrt{n} \). Note that the optimal diameter of 0.868 corresponds to a packing fraction of 59.2%. Compare this to other circle packings. Also note that the Bridson Poisson disc sampling which is not extensible in \(n\), and on typical recommended setting,  still only produces a packing fraction of 49.4%.

Note that this concept of \(d_0\) intimately connects the \( \phi_d \) low discrepancy sequences with badly approximable numbers/vectors in \(d\)-dimensions [Hensley, 2001]. Although a little is known about badly approximable numbers in two dimensions, the construction of \(phi_d\) might offer a new window of research for higher badly approximable numbers in higher dimensions.

Figure 4. Minimum pairwise distance for various low discrepancy sequences. Note that the \(R_2(\phi) \) sequence (blue) is consistently the best option, it is also the only sequence that where the normalized distance does not tend to zero as \( n \rightarrow \infty \). The Halton sequence (orange) is next best, with the Sobol (green) and Niederreiter (red) sequences not as good but still much better than random (purple). Note larger is better as it corresponds to a larger packing distance.

Voronoi Diagrams

Another method of visualizing the evenness of a point distribution is to create a Voronoi diagram based on the first \(n\) points of  a 2-dimensional sequence, and then color each region based on its area.

In the figure below, shows the color-based Voronoi diagrams for (i) the \( R(\phi_2)\)-sequence; (ii) the (2,3)-Halton sequence, (iii) prime-based recurrence; and (iv) the simple random sampling; . The same color scale is used for all figures. Once again, it is clear that the \(R(\phi_2)\)-sequence offers a far more even distribution than the Halton or simple random sampling.

This is the same as above, but colored according the the number of vertices of each voronoi cell. Not only is it clear that the \(R\)-sequence offers a far more even distribution than the Halton or simple random sampling, but what is more striking is that for key values of $n$ it only consists of hexagons!

If we consider the generalised Fibonacci sequence, $A_1=A_2=A_3=1; \quad A_{n+3} = A_{n+1}+A{n}$. That is, $A_n:$,

\begin{array}{r} 1& 1& 1& 2& 2& 3& 4& 5& 7\\ 9& \textbf{12}& 16& 21& 28& 37& \textbf{49}& 65& 86\\ 114& \textbf{151}& 200& 265& 351& 465& \textbf{616}& 816& 1081 \\ 1432& \textbf{1897}& 2513& 3329& 4410& 5842& \textbf{7739}& 10252& 13581\\ 17991& \textbf{23833}& 31572& 41824& 55405& 73396& \textbf{97229}& 128801& 170625\\ 226030& \textbf{299426}& 396655& 525456& 696081& 922111& \textbf{1221537}& 1618192& 2143648 \\

\end{array}

All values where $n= A_{9k-2}$ or $n= A_{9k+2}$ will consist only of hexagons.


For particular values of $n$, the Voronoi mesh of the $R_2$ sequence consists only of hexagons.

Figure 5. Visualization of the shape of the Voronoi diagrams based on the number of sides of each Voronoi polygon for (i) the \( R(\phi)\)-sequence; (ii) (2,3) prime based recurrence; (iii) the (2,3)-Halton sequence, (iv) Niederretier; (v) Sobol; and (iv) the simple random sampling. The colors represent the number of sides for each Voronoi polygon. Once again, it is clear that the \(R(\phi)\)-sequence offers a far more even distribution than any of the other low discrepancy sequences.

Quasiperiodic Delaunay Tiling of the plane

Delaunay Triangulation (sometimes spelt ‘Delone Triangulation’), which is the dual of the Voronoi graph, offers another way of viewing these distributions. More importantly, though, Delaunay triangulation, offers a new method of creating quasiperiodic tilings of the plane.  Delaunay triangulation of the R(\(\phi_2)\)-sequence offer a far more regular pattern than that of the Halton or Random sequences,.

More specifically, the Delaunay triangulations of point distributions where $n$ is equal to any of the generalised Fibonacci sequence $:A_N=$ 1,1,1,2,3,4,5,7,9,12,16,21,28,37,… then the Delaunay triangulation consists of  of only 3 identically paired triangles, that is, parallelograms (rhomboids)!  (excepting those triangles that have a vertex in common with the convex hull.)

Furthermore,


For values of $n=A_k$, the Delaunay triangulation of  the $R_2$-sequence form quasiperiodic tilings that each consist only of three base triangles (red, yellow, blue) which are always paired to form a well-defined quasiperiodic tiling of the plane via three parallelograms (rhomboids).

Figure 6. Visualization of the Delaunay Triangulation based on (i) the \( R(\phi_2)\)-sequence; (ii) the (2,3)-Halton sequence, (iii) prime-based recurrence; and (iv) the simple random sampling. The colors represent the area of each triangle. The same scale is used in all four diagrams. Once again, it is clear that the \(R(\phi_2)\)-sequence offers a much more even distribution than any of the other low discrepancy sequences.

Note that $R_2$ is based on $\phi_2=1.32471795724474602596090885447809$ is the smallest Pisot number, (and $\phi = 1.61803…$ is the largest Pisot number. The association of quasiperiodic tilings with quadratic and cubic Pisot numbers is not new [Elkharrat and also Masakova] , but I believe that this is the first time a quasiperiodic tiling has been constructed based on $\phi_2$ =1.324719….

This animation below shows how the Delaunay mesh of the $R_2$ sequence changes as points are successively added. Note that whenever the number of points is equal to a term in the generalised Fibonacci sequence, then the entire Delaunay mesh consists only of red, blue and yellow parallelograms (rhomboids) arranged in a 6-fold quasiperiodic manner.

(Note that the prime-based recurrence sequence also appears quasiperiodic in the weak sense that it is an ordered pattern that is not repeating. However, its pattern over a range of $n$ is not so consistent and also critically dependent on the selection of basis parameters. For this reason, we focus our interest of quasiperiodic tilings solely on the $R_2$ sequence.)

It consists of only three triangles: red, yellow, blue. Note that for this R(\(\phi_2)\)  sequence, all the parallelograms of each color are exactly the same size and shape.

The ratio of the areas of these individual triangles is extremely elegant. Namely,

$$ \textrm{Area(red) : Area(yellow) : Area(blue) }= 1 :  \phi_2 : \phi_2^2 $$

And so is the relative frequency of the triangles, which is:

$$ f(\textrm{red}) : f(\textrm{yellow}) : f(\textrm{blue}) = 1 :  \phi_2 : 1 $$

From this, it follows that the total relative area covered i nspace by these three triangles is:

$$ f(\textrm{red}) : f(\textrm{yellow}) : f(\textrm{blue}) = 1 :  \phi_2^2 : \phi_2^2$$

Figure 7.

Although the locations of red parallelograms show considerable regularity, one can clearly see that the blue and yellow parallelograms are spaced in a quasiperiodic manner.


The $R_2$ sequence is the only low discrepancy quasirandom sequence that can be used to produce quasiperiodic tilings through its Delaunay triangulation.

Covering

In the following figure, the first \(n=400\) points for each two-dimensional low discrepancy sequence are shown. Furthermore, each of the 20×20=400 cells are colored  green only if that cell contains exactly 1 point.  That is, more green squares indicates a more even distribution of the 400 points across the 400 cells. The percentage of green cells for each of these figures is: \(R_2\) (64%), Halton (49%), Kronecker (39%), Neiderreiter (38%), Sobol 56% and Random (38%).

figure 8. A comparison of various discretised 2-dimensional quasirandom low discrepancy sequences. More green is better because it corresponds to a more even distribution.

Multiscale Low Discrepancy

Some of the low discrepancy sequences exhibit what called be termed ‘multiscale low discrepancy’.

So far we have assumed that when we need to distribute $n$ points as evenly as possible that all the points are identical and indsitrguishable. However, in many situations there are different types of points.

Consider the problem of evenly distributing the $n$ points so that not only are all the points evenly separated, but also all points of similar type are evenly distributed.

Specfically, suppose that if there are $n_k$ points of type $k$, (where $n_1+n_2+n_3 +… +n_k= n$) , then a multiset low discrepancy distribution is where each of the $n_k$ points are evenly distributed

In this case, we find that the $R$-sequence and the Halton sequence can be easily adapted to become multiset low discrepancy sequences, simply allocating consecutively allocating the points of each type.

The figure below shows how $n=150$ points have been distributed such that 75 are blue, 40 green 25 green and 10 red. For the additive recurrence sequence this is trivially achieved by simply setting the first 75 terms to correspond to blue, the next 40 to orange, the next 25 to green and the last 10 to red. This technique works  almost works for the Halton and Kronecker sequences but fares very poorly for Niederreiter and Sobol sequences. Furthermore, there are no known techniques for consistently generating multiscale point distributions for Niederreiter or Sobol sequences.

This shows that multi-set point distributions such as those in the eyes of chickens, can now be described and constructed directly via low discrepancy sequences.


The $R_2$ sequence is the only low discrepancy quasirandom sequence that admits a simple construction of multiscale low discrepancy.

Figure 9. Multiscale Low discrepancy Sequences. For the $R$ sequence, not only are all the points evenly distributed, but also points of each particular color are evenly distributed.

Quasirandom Points on a sphere

It is very common in the fields of computer graphics and physics to need to distribute points on the surface of a 3-sphere as evenly as possible. Using open (infinite) quasi-random sequences, this task merely mapping (lifting) the quasi-random points that are distributed evenly in the unit square onto the surface of a sphere via a Lambert equal-area projection. The standard Lambert projection transformation that maps a point \( (u,v) \in U[0,1] \to (x,y,z)\in S^2\), is :

$$ (x,y,z) = (\cos \lambda \cos \phi, \cos \lambda \sin \phi, \sin \lambda), $$

$$ \textrm{where} \quad  \cos \lambda = (2u-1)-\pi/2; \quad \phi = 2\pi v. $$

As this \(\phi_2-\)squence is fully open, it allows you to map an infinite sequence of points onto the surface of a sphere – one point at a time.

This is in contrast to other existing methods such as the Fibonacci spiral lattice which requires knowing in advance, the number of points.

Once again, through visual inspection, we can clearly see that for \(n=1200\), the new \(R(\phi_2)\)-sequence is far more even than the Halton mapping or Kronecker sampling which in turn is far more even than random sampling.

Figure 10.

Computer Graphics: Dithering

Most state-of-the-art dithering techniques (such as those based on Floyd-Steinberg) are based on error-diffusion, which are not well suited for parallel processing and/or direct GPU optimisation. In these instances, point-based dithering with static dither masks (ie fully independent of the target image) offer excellent performance characteristics. The most famous, and still widely used dither masks are based on Bayer matrices, however, newer ones attempt to more directly mimic blue noise characteristics.

The non-trivial challenge with creating dither masks based on low discrepancy sequences and/or blue noise, is that these are low discrepancy sequences map an integer \(Z\) to a two-dimensional point in the interval \( [0,1)^2 \). In contrast, a dither mask requires a function that maps two-dimensional integer coordinates of the rastered mask to a real valued intensity/threshold in the interval \( [0,1) \). A reasonable option for the \(R-\)sequence is to use the following approach.

For each pixel (x,y) in the mask, we set its intensity to be \(I(x,y) \), where:

$$ I(x,y) = \alpha_1 x + \alpha_2 y \; ( \textrm{mod} 1); $$

$$ \textrm{where} \quad \pmb{\alpha} =(\alpha_1,\alpha_2) = (\frac{1}{\phi_2}, \frac{1}{\phi_2^2}), $$

$$ \textrm{and} \; \phi_2\ \textrm{is the unique positive root of }  x^3\;=x+1. $$

Obviously, as this function is effectively linear, there is obvious room for improvement as it suffers from banding lines in the direction of \(   v =(\alpha_1,\alpha_2 ) \). Notwithstanding this concern, the results show that it is competitive, if not better than other current methods.

The results of this dither matrix are shown below for the classic 256×256 “Lena” test image, as well as a checkered test pattern. The results using the standard Bayer dither masks is also shown, as well as one based on blue noise and presented at Siggraph 2016. [Georgive & Fajardo, 2016].

One can see that those dithered with Blue noise, suppress detail such as boundaries between color changes and produce a somewhat ‘grainy’ effect, whereas the Bayer dither exhibits noticeable white dissonance in the light gray areas. The $R$-sequence dither is the only one that has both smooth gradients as well as clarity in highlights of Lena’s face and shoulder.

Therefore although the blue noise masks is the only mask where fourier spectra shows a clear black central disk  – an attribute which in theory is highly desirable – in practice the overall effect is only fair.  The fourier spectra of the $R$-mask is consistent with the fact that the Delaunay triangulation of the canonical R-sequence consists of aperiodic tiling of three parallelograms.

Figure 11. $R_2$-sequence (left); standard 256-Bayer (middle) ; and blue noise dither mask (right). Looking at the test patterns, it can be seen that the R-sequence dither has less dissonant gradients than the Bayer dither (eg light gray with white dots). Further, both the test pattern and Lena reveal that the R-dither produces a less grainy image the blue noise dither. This grainy effect means that the blue noise mask also masks much of the highlights and shadows on Lena’s face and shoulder.

Higher Dimensions

Similar to the prior section, but for five (5) dimensions, the graph below shows the (global) minimum distance between any two points, for the \( R(\phi_5)\)-sequence, the (2,3,5,7,11)-Halton and the random sequences. This time, the normalised minimum distance measure is normalized by a factor of \( 1/ \sqrt[5]{n} \). One can see that, as a consequence of  ‘the curse of dimensionality’, the random distribution is better than all the low discrepancy sequences — except the $R_2$-sequence. 

The \(R(\phi_5)\)-sequence , even for \( n \simeq 10^6\) points, the minimum distance between two points is still consistently around \( 0.8/\sqrt{n} \), and always above \( 0.631/\sqrt{n} \).

The $R_2$ sequence is the only $d-$dimensional low discrepancy sequence, where $d_0$ only falls at a arate of $1/n^(1/d)$.

Figure 12. This shows that the R-sequence (blue) is consistently better than the Halton (orange); Sobol (green); Niederreiter (red); and random (purple). Note larger is better as it corresponds to a larger packing distance.

Numerical Integration

The following graph shows the typical error curves \( s_n = |A-A_n| \) for approximating a definite integral associated with a gaussian function of half-width \(\sigma = \sqrt{d}\), \( f(x) = \textrm{exp}(\frac{-x^2}{2d}), \; x \in [0,1]  \), with: (i) \( R_{\phi} \)  (blue); (ii) Halton sequence, (orange); (iii) randomly  (green); (iv) Sobol  (red) 

It shows that for \(n=10^6\) points, the differntial between random sampling is the Halton sequence is somewhat less now. However, as was seen in the 1-dimensional case, the $R$-sequence and the Sobol are consistently better than the Halton sequence. It also suggests that the Sobol sequence is marginally better than the $R$-sequence.

Figure 13. Quasirandom Monte Carlo methods for 8-dimensional integration. Note smaller is better. The new R-sequence and Sobol peform consistently better than the Halton sequence.

Other related topics

  • Quasirandom Walks / Propp rotor
  • Combinatorial discrepancy
    • probabilistic taguchi methods / factor analysis for large d, [draft]
    • low discrepancy subsets
    • Quasirandom permutations
  • Hyperuniformity,
  • Machine learning,
    • quasirandom stochastic descent
    • extreme learning machines
  • quasicrsytals


My name is Martin Roberts. I have a PhD  in theoretical physics. I love maths and computing. I’m open to new opportunities as a data scientist – consulting, contract or full-time – so let’s have a chat on how we can work together! My contact details can be found here.

Other Posts you may like