SZ Sequences: Binary-Based $(0,2^{q})$-Sequences

The term “discrepancy” refers to the error of using a point set to estimate the volume of region by counting the ratio of points falling inside. Star discrepancy, commonly notated as $D^{*}$, refers to the maximum discrepancy measured over all axis-aligned hyper-rectangles in the concerned domain, and gives a reliable measure of the performance of the point set in numerical integration. A point set is considered a low-discrepancy (LD) set if it attains $\mathcal{O}\left(\log^{s-1}(N)/N\right)$ discrepancy for $N$ points in $s$-dimensions. A low-discrepancy sequence is an infinite sequence of points that provides an LD set for certain contiguous blocks of points of certain sizes, and attains $\mathcal{O}\left(\log^{s}(N)/N\right)$ discrepancy for all contiguous blocks.

A significant construction of a one-dimensional LD sequence was presented by van der Corput (1935). If

denotes the binary encoding of the sample index $v$, then the $v$th sample is computed by mirroring the bits of $v$ at the fraction point:

This approach can be generalized to higher dimensions by generating the $k$th component/dimension of the sample using a distinct linear mapping

where $V$ is a vector representing the $v$th digit-reversed van-der-Corput-like sample in some base $b$, $X^{(k)}_{v}$ is a vector representing the same-base fractional digits of the final sample location along the $k$th axis, and $C^{(k)}$ is a matrix in Galois Field (GF) $b$, known as a generator matrix, that linearly maps reversed digits of the sample index to digits of the computed sample, hence the name “digital” to collectively describe these methods. Next we briefly outline the three most established digital construction approaches to LD sequences.

The first known higher-dimensional extension to van der Corput sequences in due to Halton (1960), and simply uses a different prime base for encoding the index $v$ along each axis. Identity matrices may be used for $C^{(k)}$ in Eq. (3), but more complex scrambling matrices have been proposed to improve the pairwise 2D projections.

A very different approach was taken by Sobol (1967), who stuck to a binary-base encoding, and constructed a different generator matrix for each dimension, using a recurrence operation to compute each column of the matrix from the preceding columns, taking the coefficients from an irreducible polynomial. These polynomials, in a sense, replace the role of different bases in Halton’s construction. We skip the details and only highlight the most relevant one here that the first two irreducible polynomials lead to the pair

that constitutes the first two, invariant, and most well-known dimensions of Sobol’s construction, where

happens to be the Pascal matrix in GF(2), where $i$/$j$ are row/column indices, indexed from 0. Please note that we hereinafter adopt the convention of showing binary matrices as arrays of (0) gray and (1) black dots for better visibility, and also occasionally using dots instead of zeros for the same purpose.

These first two dimensions of Sobol inspired a third construction by Faure (1982) that uses powers of a Pascal matrices

in a prime base $b$, and generates a $b$-dimensional LD sequence analogous to the 2D Sobol sequence, in a sense explained in the following subsection. Niederreiter (1987) extended Faure construction to power-of-prime bases $b^{q}$ using symbolic matrices and arithmetic over the GF($b^{q}$).

Thus, we have three well-established constructions approaches to high-dimensional LD sequences:

1. (1)

   Halton, using different prime bases.

2. (2)

   Sobol and variants, using binary-polynomial-based matrices.

3. (3)

   Faure and Niederreiter extension, using Pascal matrices over fields.

Among the three, two factors made Sobol the most adopted one in computer graphics: a) the binary base, leading to extremely fast computation, and b) that higher distribution quality is obtained for lower dimensions, as will be explained in the following subsection. Halton comes next, sharing the second but not the first property of Sobol. Because the prime bases are known at compile-time, the cost of integer divisions and modulus operations can be reduced (Warren, 2012), bringing the digit reversal of Halton to within tolerable limits; this has led to making Halton sequences being used in common rendering platforms like pbrt (Pharr et al., 2016).

In contrast, we are unaware of any mention of Faure sequences and their derivatives in CG, and we may attribute this primarily to the large base used for all dimensions, requiring a considerably large number of points to attain the advertised uniformity. The lookup-based arithmetic of Niederreiter’s extension exacerbate this by adding considerable computational overhead.

While our work is developed from first principles, the resulting construction may be considered as an extension to Faure’s, specifically Niederreiter’s extension, that addresses both limitations: we ensemble higher-dimensional sequences from lower-dimensional ones, and use GF(2) matrices like Sobol’s. Thus, our construction effectively introduces the third category of LD sequences to CG.

Niederreiter (1987) also established a complete theoretical framework for developing and studying a large class of LD constructions, including the radix-based ones listed above. At the heart of this theory is so-called $(t,m,s)$-Nets in base $b$, which refers to a multi-stratified distribution of $b^{m}$ points in an $s$-dimensional hypercube such that, for all possible stratification (slicing) ways of the domain into $b^{m-t}$ similar rectangular strata (cells), we find exactly $b^{t}$ points in each stratum. An illustration of the essence of a $(0,4,4)$-net in base 4, for example, is shown in Fig. 1. Then come $(t,s)$-sequences in base $b$, which are infinite sequences of $s$-dimensional points such that, for all integer $m$, the first and all subsequent blocks of $b^{m}$ points comprise a $(t,m,s)$-net in base $b$. While these definitions are independent of the construction method, the model possibly existed thanks to the existence of algebraic recipes for construction. Specifically, the definition of $(t,m,s)$-nets translates directly into a simple condition on digital: if we build a hybrid $m\times m$ matrix

by combining the first $m_{k}$ rows of the $k$th generator matrix in Eq. (3), with $m_{k}$s summing to $m$, then each such a hybrid matrix multiplied by $V$ computes the allocation of the points to strata in a specific stratification. The $(0,m,s)$-net condition, for example, then corresponds to requiring that each hybrid matrix is invertible, ensuring a one-to-one correspondence of points to strata. A $(0,s)$-sequence requires a set of generator matrices that maintains this hybrid matrix condition progressively over their top-left corners.

These definitions suggest that it is desirable to keep both the base $b$ and the parameter $t$ small, though different needs might call for different configurations. To give some examples, the van der Corput construction is a (0, 1)-sequence in base 2. Faure sequences are $(0,b)$-sequences in their respective prime bases, while Niederreiter’s extension to Faure creates $(0,b^{q})$-sequences in power-of-prime bases $b^{q}$. The first two dimensions of Sobol sequence comprise a $(0,2)$-sequence in base 2, but the $t$ value increases as we take more dimensions, which makes lower dimensions attain better quality, as we mentioned in Sec. 2.2. Halton sequences, in contrast, do not readily fit within the $(t,s)$-sequences model, unless we adapt it to accept different bases. They still exhibit a multistratified structure, though. Hence, the common ground of all the three families of constructions is that they attain optimal quality at powers of the base, or product of bases in case of Halton, and here we can see the reason of superiority of Sobol and Halton over Faure for graphics applications where lower dimensions are usually more important.

In this paper we develop a novel construction of $(0,2^{q})$-sequences using binary GF(2) matrices. Technically speaking, we are breaking two traditional rules of nets and sequences: that a composite number cannot directly be used as a base, and that one can not obtain a $(0,s)$-sequence with a base $b<s$. We achieve this result, however, by combining the two violations: instead of using a base-4 directly, for example, we use $2\times 2$ binary matrices as building blocks to emulate base-4 matrices. Please note that $b=s$ is implied in the following if no base is mentioned.

Beyond the algebraic recipes for constructing LD sets and sequences, the design space is greatly enriched by scrambling techniques that can derived new valid nets/sequences from given ones. Tezuka (1994) showed that multiplying arbitrary lower-triangular matrices from left by the generator matrices produces valid generators, since doing so preserves that ranks of all hybrid matrices. Owen (1995), in contrast, introduced a post-processing technique, “Owen scrambling”, that preserves the net/sequence structure. In contrast to scrambling that applies to sample locations, the term shuffling refers to reordering (permuting) the sample indices (Faure and Tezuka, 2002). As $(0,s)$-sequences, all these are applicable to our sequences, and we will use them in our development.

The interest in LD construction has grown in the graphics community over the past decade, mostly devoted to finding handles to impose control over these distributions as commonly needed in graphics applications. Ahmed et al. (2016) presented an original 2D low-discrepancy construction that enables imposing a blue noise spectrum, then Perreir et al. (2018) imposed a blue noise spectrum onto selected pairs of dimensions of a Sobol sequence. While both methods compromised the discrepancy of the underlying distribution, they presented a proof of concept that encouraged further research.

In a different path, Christensen et al. (2018) tried searching for dyadic (i.e., base-2) $(0,2)$-sequences sequences with good spatial and spectral properties, and Pharr (2019) presented a more efficient search, but Ahmed and Wonka then presented a theorem ([)Theorem 4.2]Ahmed2021Optimizing implying that the search space is actually confined to Owen-scrambled 2D Sobol, and soon the first group of authors, lead by Helmer (2021), developed an extremely fast implementation of Owen scrambling. We take this as an excellent example of the fast-paced learning cycles of the community to analyze and manipulate the inherently difficult-to-understand LD construction. In this paper we present an almost complete cycle from experimentation to elements of a theory.

In more recent years, the research on LD in CG reached the state of enriching the quasi-Monte-Carlo (QMC) literature with significant findings and novel constructions, including an efficient algorithm that spans the whole space of dyadic $(0,m,2)$-nets (Ahmed and Wonka, 2021), a cascade of dyadic $(0,m,2)$-nets over successive pairs of dimensions (Paulin et al., 2021), a family of self-similar dyadic $(0,2)$-sequences (Ahmed et al., 2023), a gradient-decent optimization of discrete Owen scrambling (Doignies et al., 2024), and a base-3 analog to Sobol sequences (Ostromoukhov et al., 2024). Our current work is partially inspired by this last work of Ostromoukhov et al., where we thought of using base 4 instead to gain the advantage of binary computation.

For all combinations of valid initial 2$\times$2 blocks of the three matrices, we consistently obtained 768 distinct valid sets for the first extension to 4$\times$4, and 384 distinct valid extension on each branch for subsequent extension to 6$\times$6 and 8$\times$8, which suggests an analytical construction model. The fixed number of alternatives in each extension step indicates that the degrees of freedom are possibly limited to the top and bottom blocks of the extension columns. Following this hint, we note that the corner extension block $X$ in Eq. (8) is not involved in any of the hybrid matrices, hence comprises, for each of the three matrices, a free choice among the four options in Eq. (12). This accounts for a factor of $4^{3}=64$, and we verified that fixing the extension corner blocks $X$ to identity reduces the list to 12/6 sets in the first/subsequent extension steps. This leaves us with the degrees of freedom in the top row, which we analyze next.

Let

denote the first block-row in one matrix $X$ of the three upon the $j$th extension, where each $X_{i}$ is a 2$\times$2 block matrix. By constructing a hybrid matrix comprising $j$ rows of the first identity matrix and this row of $X$, we find that

leading to the conclusion that, for each of the three matrices, all $2\times 2$ blocks of the first block row have to be invertible.

Now consider taking a hybrid of the first $j-1$ rows of $I$ with the first row of each of two other matrices $A$ and $B$, then we get

Since all four elements are invertible by virtue of Eq. (15), we can use our earlier “$RA^{-1}C+X$” reduction of Eq. (9) to get

or

Let us define a symbol

to designate the product of the inverse of the $j-1$st by the $j$th entry in the first row of block-matrix $X$. Then Eq (18), rewritten as

says that the sum of corresponding symbols in the first row of each pair $\{A,B\}$ of generator matrices must be invertible, which also implies that the symbols at corresponding slots must be distinct for the three matrices. Thus, for each slot in the first block row, the inserted block-matrix symbols must be

1. (1)

   Invertible.

2. (2)

   Distinct for each matrix.

3. (3)

   Have invertible sums.

The population of six invertible 2$\times$2 matrices comprises two disjoint sets,

that satisfy these requirements, each comprising two diagonally-mirrored triangular elements, and a counter diagonal element that switches between them via addition, hence we call them $S$ (think ) and $Z$ (think ) to make it easier to recall which is which. Hence comes the name SZ for these sequences. We will also stick hereinafter to this “S/Z” convention in discriminating left/right-hand elements and operations, respectively.

We conclude this subsection by noting that these three simple rules comprise all the requirements for satisfying the hybrid matrix condition for the first $2\times 2$ block matrices. Further, they analogously extend to other powers of 2. Thus, constructing binary matrices for generating a “$(0,2,2^{q})$-net in base $2^{q}$” reduces to finding an alphabet of $q\times q$-block matrices satisfying these rules. This name in Niederreiter’s notation just means a set of $2^{2q}$ points in $2^{q}$ dimensions that is stratified in all $2^{q-1}(2^{q}-1)$ pairs of 2D projections; cf. Jarosz et. al. (2019). Next, we show how to extend the matrices to arbitrary sizes for $(0,2^{q})$-sequences in base $2^{q}$.

Instead of defining the symbols in Eq. (19) in terms of adjacent entries, it would be more constructive to populate the matrices directly with symbols. Realizing that all first-row blocks must be invertible helps us to achieve this via this subtle factoring

The common $S$ (left-hand-side) factor represents a Tezuka scambling, hence may be dropped. Thus, without loss of generality, we may populate all blocks in the first row of the second-dimension (first of the three) matrix with identity matrix blocks. Once we implement this last factoring, the exhaustive search reduces to this one-and-only quadruple of generator matrices:

here shown to 32-bit depth. As per our exhaustive search, this quadruple represents the template of all binary-matrix-based (0, 4)-sequences, and all other possible variants may be derived by plugging back the degrees of freedom factored earlier, as we will discuss later.

Upon careful inspection of these matrices, the first important observation is that they are built exclusively from the $S$ block symbols in Eq. (21), plus empty (0) blocks. Further inspection, guided by a preliminary assumption that these sequences bear some relation to Faure’s, reveals that these are actually Pascal matrices

constructed from these four block symbols.