Quad-Optimized Low-Discrepancy Sequences

Integrating functions is at the core of many computer graphics – and other – systems. Perhaps the most common example in our community is the case of physically-based rendering, that involves integrating the rendering equation. This is most often performed through Monte Carlo integration, involving the evaluation of the integrand at random locations, or quasi-Monte Carlo integration, evaluating the integrand for pseudo-random well-located samples. Quasi-Monte Carlo integration benefits from Koksma-Hlawka's inequality: the convergence rate of the Monte Carlo integration estimator improves with the integrand regularity, and with samples uniformity. While the regularity of the integrand itself can be somewhat controlled via other variance reduction techniques, we focus on samples uniformity that does not require knowledge of the integrand. Over the years, much work has been devoted to producing extremely uniform samples that effectively improve Monte Carlo convergence rate, e.g., by optimizing samples according to various criteria. Still, Sobol’ sequences with Joe and Kuo initialization tables [Joe and Kuo 2008; Sobol’ 1967] and Owen scrambling [Owen 1995] remain used in many production rendering systems and simulators, despite suboptimal sample uniformity. This is due to many other advantages offered by this approach that are often overlooked by other researchers. First, Sobol's method produces a sequence of arbitrary length: this allows for progressively refining a rendering by adding samples while maintaining good Monte Carlo convergence. In many cases, competing methods only offer point sets that do not allow for progressive rendering, or provide sequences only up to a finite predetermined number of samples. Second, it works in (almost) arbitrary dimensions while many approaches remain low-dimensional or even 2-d. Third, samples are fast to generate, with a compact representation (a small binary matrix per dimension and low-order polynomial binary coefficients), while many other approaches require complex optimizations or storing the precomputed points themselves. Nevertheless, other production renderers may only use the first few dimensions of Sobol’ sequences due to their “increasingly poor distributions in pairs of higher dimensions” (regarding Renderman [Christensen et al. 2018]).

Based on this observation, we seek a way to improve Sobol’ samples uniformity on consecutive 2-d and 4-d projections while preserving their advantages. Sobol’ samples uniformity is guaranteed in part by construction, but also via appropriate parameter search. We leverage both aspects. First, we benefit from the additional degrees of freedom offered by a construction in base 3 rather than base 2 [Paulin et al. 2022]. Second, we benefit from other degrees of freedom offered by a search over irreducible polynomials instead of primitive polynomials [Faure and Lemieux 2016]. The relatively small number of irreducible polynomials still allows us for a comprehensive search, and allows us for finding solutions that provide the following desirable properties, that will be detailed in this paper. We show one quadruplet of dimensions (or quad hereafter, for short) that is a (1, 4)-sequence with (0, 2)-sequences pairwise projections, proven by construction. This provides ideal uniformity as measured in terms of discrepancy. We provide 11 additional quads of dimensions with the similar properties but only numerically verified up to nearly 59k samples. All quads we present are compatible: their dimensions do not overlap, and we use them together to form higher-dimensional sequences that posess (t, s)-sequence properties guaranteed by the Sobol’ construction, in up to 48 dimensions. We show that the discrepancy of our sequence is improved over Sobol’ on these consecutive projections. In exchange, a base-3 construction prevents the use of fast xor-based binary arithmetic, and can incur minor additional costs.

In its simplest form, quasi-Monte Carlo integration consists in numerically estimating an integral as

Much work has been devoted to improving the uniformity of generated point sets for improving the integration error.

Point sets.. Point sets have been numerically optimized specifically for improving uniformity. Optimal transport has been used as a measure of uniformity to minimize [De Goes et al. 2012] but this has remained limited to 2 or 3-dimensions. A sliced optimal transport variant allows to reach about 20 dimensions [Paulin et al. 2020; Salaün et al. 2022]. By instead optimizing the variance of the function obtained by summing Gaussians centered at each sample, Ahmed et al. [Ahmed et al. 2022] reach excellent uniformity, demonstrated up to 8 dimensions. Closer to our work, the cascaded construction of Paulin et al. [2021] based on alternating Sobol’ process with bit reversal, guarantees a perfect uniformity of the point set for projections onto consecutive pairs of dimensions in terms of discrepancy. The discrepancy of the first nested 100 dimensions is numerically optimized, although this does not guarantee a low discrepancy behavior (see next for a discussion on low-discrepancy sequences). All these approaches offer excellent uniformity, which translates into excellent numerical integration convergence, through Koksma-Hlavka's theorem but also other theorems bounding integration error with different kinds of discrepancy measures [Harman 2010], optimal transport metric [Paulin et al. 2020], or Fourier spectra behavior [Pilleboue et al. 2015]. Their main drawback is that they do not allow for progressive evaluation of the integral: once a numerical estimate has been obtained using a set of samples, refining this computation to improve its accuracy is not directly possible. This requires throwing away the previous result and producing a new completely different point set with more samples. Recently, Ahmed et al. [2023] introduced a method that allows converting certain point sets to sequences up to a finite predetermined number of samples. It is unfortunately limited to two dimensions, and using this process independently on different pairs of dimensions would deteriorate high-dimensional uniformity (as illustrated in Sec. 5 with ZeroTwo and Padded samplers). It is therefore not the solution we are looking for in higher dimensions.

Low-discrepancy sequences (LDS). Matrix-based constructions allow to produce sequences of points using finite field arithmetic. Among others, the Sobol’ sequence [Sobol’ 1967] has been most widely used. It produces points via modulo 2 matrix-vector multiplications, where the per-dimension binary triangular matrix is in part determined as a preprocess by optimizing generators (the first e × e upper left part of these matrices, where typically e ∈ [1, 18] depending on the dimension) and the rest of each matrix is obtained using a recursion involving a linear combination of previous columns weighted by the coefficients of a well-chosen primitive polynomial. While all these generators and primitive polynomials were carefully optimized, the obtained uniformity has remained limited (see Fig. 1 and 4) due to several constraints: the base 2 construction offers limited degrees of freedom and the set of primitive polynomials is small. However, they satisfy the important property of (t, s)-sequences, described hereafter. A point set of n = b^m samples within the s-dimensional unit hypercube, where b is some base (b = 2 for Sobol), is said to be a (t, m, s)-net if all subintervals of volume b^t/n of the form $\Pi _{j=1}^s [a_j b^{-d_j}, (a_j+1) b^{-d_j}]$ for all j contain b^t samples, with $0\le a_j < b^{d_j}$, and $a_j, d_j \in \mathbb {N}_0$. Intuitively, t = 0 enforces high uniformity since it requires many small intervals each to contain a single point, while larger t values only enforce larger intervals to contain more points. A (t, s)-sequence is a sequence that enforces the (t, m, s)-net property for all possible values of m.

An advantage of (t, s)-sequence is that their discrepancy —and hence integration error— is bounded as a function of t [Niederreiter 1992] as:

In MatBuilder [Paulin et al. 2022], an integer linear program is designed to satisfy various user constraints on generative matrices. Contrary to Sobol’, these matrices are directly optimized, without resorting to polynomials, and they show that base-3 constructions offer significantly improved results due to the largely increased degrees of freedom. We take inspiration from this work to optimize Sobol’ generators in base b = 3. However, for a set of m × m optimized matrices, properties are only enforced up to b^m samples. They hence cannot strictly produce “sequences” (where properties ought to hold for an arbitrary number of added points). We call these samplers “progressive” as in Ahmed et al. [2023]. In the work of Faure and Lemieux [2016], a Sobol’ construction based on irreducible polynomials instead of primitive polynomials is used for higher degrees of freedom. This construction has subsequently been used to numerically optimize generators in base 2 [Faure and Lemieux 2019]¹. We also benefit from their result, seeking solutions for a subset of irreducible polynomials.

Other constructions have been proposed to produce low-discrepancy sequences, notably benefitting from the extremely high uniformity of low-dimensional Sobol’ sequences. The ZeroTwo approach presented in PBRT [Pharr et al. 2023] uses the first two dimensions of the Sobol’ sequence, and repeat them along with a random permutation of points to produce higher dimensional point sets. While pairing randomly loses the sequence property, shuffling indices instead with an Owen permutation allows to remain a sequence [Burley 2020; Helmer et al. 2021]. The first four dimensions can similarly be used [Burley 2020]. While these produce sequences with perfect projections on consecutive pairs of dimensions, their higher-dimensional behavior is not LDS anymore and loses benefits for general higher-dimensional integration problems [Helmer et al. 2021].

For completeness, this section describes the extension of Sobol’ sequences as described in the work of Faure and Lemieux [2016].

The construction is based on the following recipe. The general idea is to decompose the index of a point (say, point number i) in a chosen prime base b: ${\bf i}=(i_0, i_1,..., i_{m-1})_b = \sum _{k=0}^{m-1} i_kb^k$ and use this decomposition to produce point coordinates for each dimension. For a given dimension d, a matrix C_d with coefficients in {0,.., b − 1} is built such that the coordinate of that point for that dimension is obtained with the help of the matrix-vector multiplication j = C_di^T. The d’th coordinate of point i is then expressed as x_{i, d} = ∑_kj_kb^{− k − 1}.

The power of Sobol’ construction lies in the way matrices C_d are obtained and the properties they guarantee on the generated points. Matrices C_d are upper triangular, with a non-zero diagonal. Its upper left e_d × e_d block consists of a set of e_d · (e_d − 1)/2 degrees of freedom ; these values are optimized in order to achieve good properties on the resulting points, as in the work of Joe and Kuo [Joe and Kuo 2008]. These first e_d columns are called generating vectors. Typically, e_d increases with the dimension d (e₁ = 1 while e₁₀₀ = 9 in the work of Joe and Kuo [Joe and Kuo 2008]).

Each column {V_n}_{n = e + 1..m} of C_d is obtained as a linear combination of its e_d previous columns, plus a shifted column (see Fig. 2). The coefficients in this combination are given by the coefficients of a polynomial $p_d(x) = \sum _{i=0}^{e_d} a_i x^i$ of degree e_d:

While in the initial work of Sobol’ these polynomials were taken as primitive polynomials, Faure and Lemieux [2016] showed that similar guarantees were obtained using a larger set consisting of irreducible polynomials. Irreducible polynomials simply correspond to polynomials that cannot be divided by other polynomials other than themselves and 1 (similarly to prime numbers), while a primitive polynomial of degree m in base b carries the additional restriction that the smallest positive integer n such that it divides xⁿ − 1 is n = b^m − 1.

In particular, it is shown that for s matrices {C₀,..., C_{s − 1}} obtained using irreducible polynomials of degrees {e₀,..., e_{s − 1}}, the value for t is bounded by $t=\sum _{d=0}^{s-1} (e_d-1)$ [Faure and Lemieux 2016]. This in turn guarantees the low discrepancy of this sequence, though in practice, lower values of t may be found for specific generating vectors. This larger choice of polynomials (over primitive polynomials) allows to reduce the value of t because the degree of available polynomials grows much more gracefully. Similarly, the choice of a basis b = 3 is guided, in our work, by the added degrees of freedom hinted at by Paulin et al. [2022] as well as the slower growth of irreducible polynomials degree.

As an example, the first few primitive polynomials in base b = 3 are x + 1, x² + x + 2, x² + 2x + 2, x³ + 2x + 1, while the first few irreducible polynomials in base b = 3 are x, x + 1, x + 2, x² + 1. Similarly, the first few irreducible polynomials in base b = 2 are x, x + 1, x² + x + 1, x³ + x + 1. For primitive polynomials in base 3 and for irreducible polynomials in base 2, the fourth polynomial is already cubic while for base-3 irreducible polynomials it is only quadratic.

While in the case of b = 2, Joe and Kuo have extensively searched for good generating vectors with primitive polynomials [Joe and Kuo 2008] and Faure and Lemieux for irreducible polynomials [Faure and Lemieux 2019], such work has not been conducted for base 3 irreducible polynomials. The generating vectors we propose give desirable properties and benefit from Faure and Lemieux’ advances on Sobol’ sequences theory.

We first list the properties obtained by our sequences. The first 4 dimensions exhibit the following exceptional properties:

• (a) The first pair of dimensions (0, 1) is (0, 2)-progressive (up to m = 100)
  
• (b) The second pair (2, 3) is a (0, 2)-sequence (up to an infinite number of points)
  
• (c) Pairs (0, 2) and (0, 3) are (0, 2)-sequence (ditto)
  
• (d) Pairs (1, 2) and (1, 3) are (1, 2)-sequence (ditto)
  
• (e) The quad (0, 1, 2, 3) is a (1, 4)-sequence (ditto)
  

The other quads (4i, 4i + 1, 4i + 2, 4i + 3) exhibit the following good properties:

• (f) The first pair (4i, 4i + 1) is (0, 2)-progressive (up to m = 10, i.e., up to 3¹⁰ ≈ 59k points)
  
• (g) The second pair (4i + 2, 4i + 3) is (0, 2)-progressive (ditto)
  
• (h) The quad (4i, 4i + 1, 4i + 2, 4i + 3) is at most (2, 4)-progressive (ditto)
  

The resulting sequence is a (t, s)-sequence, for some t, for any subset of dimensions – this is granted by design due to Sobol’ properties.

We analyze our sequences in terms of their 2-d projections, in terms of higher-dimensional discrepancies, integration error on synthetic integrands and integration error for rendering applications. Extensive experiments can be found in our supplementary document for all pairs and quads of dimensions.

5.2 2-d projections

We illustrate the uniformity of our sequences in all pairs of projections in Fig. 4 as compared to Sobol, ZeroTwo (the first two dimensions of Sobol’ with random permutations), the ISN-dec LDS sequence optimized by Faure and Lemieux [2019], the first four projections of Sobol’ with random permutations as used by Burley [Burley 2020] for rendering, and the cascaded Sobol’ construction of Paulin et al. [Paulin et al. 2021]. The approach of Faure and Lemieux [2019] results in more uniform projections for distant pairs of dimensions while we improve uniformity on adjacent pairs of dimensions. Adjacent pairs of cascaded Sobol’ appear more uniform, but non-adjacent pairs are much less uniform. The ZeroTwo and Padded Sobol’ (0123) approaches result in white noise behavior for distant pairs of dimensions. Note that while we illustrate 2-d projections for 243 (resp. 256 in base 2) samples, behavior may change for other sample counts. We guarantee t = 0 for consecutive 2-d projections and t = 1 or t = 2 for consecutive 4-d projections at least for all sample counts until m = 59, 049, but t = 0 may be occasionally achieved by our sampler and others as well for a given sample count and pair of dimension. Also, t = 0 in base 2 involves a better uniformity than t = 0 in base 3, as discrepancy generally grows with b (Eq. 1). While some pairs of non-optimized projections appear very non-uniform (e.g., dimensions (7,11)), their 2-d discrepancy remains counterintuitively lower than that of white noise (see Fig. 5 in supplementary materials). If such structures need to be avoided, ZeroTwo or Padded Sobol’ remain good candidates in consecutive 2-d or 4-d projections but at the cost of lower convergence rate (see Sec. 5.3).

We evaluate the 2-d discrepancy and integration error for Gaussian and Heaviside integrands of 2-d consecutive projections as the number of samples increases. Fig. 3 (top) shows this data for two pairs, while graphs for all consecutive pairs are shown in supplementary materials. The supplementary materials also showcase an example of discrepancy above 3¹⁰ samples (up to 531k samples), and for sample counts different from powers of 3. Cascaded Sobol’ [Paulin et al. 2021] offers best uniformity on consecutive 2-d projections as (0,2)-net property is enforced by construction, though their method does not result in a sequence. Aside from cascaded Sobol’, as dimensionality increases, our approach offers better uniformity and integration error over competitors, including Faure and Lemieux [2019]. Fig. 1 (a) summarizes all 2-d discrepancy curves on the same graph for our approach and Sobol’ for comparison.

5.4 Rendering

We performed rendering experiments, in 6-d and 10-d. In the first experiments shown in Fig. 5 (top two rows), the scenes only have direct illumination. The first two dimensions are used for sampling light sources and the next two dimensions for sampling the BRDF ; a multi-sample MIS estimation is used to combine these estimators. The next two dimensions are used to sample within each pixel. The second experiment, shown in Fig. 5 (bottom row), adds one bounce of indirect lighting, involving 4 additional dimensions. Here, the first two dimensions are used for sampling within pixels, the next two dimensions for sampling light sources, the next two dimensions for sampling the BRDF (used both for MIS estimation of direct lighting like in the 6-d case, and for generating an indirect ray if the direct ray did not find a light source), the next four dimensions are used similarly at the next light bounce (sampling light sources and BRDF). In practice, dimension order did not significantly affect results. For direct light only, we improve the rendering error over other state-of-the-art samplers [Burley 2020; Faure and Lemieux 2019; Paulin et al. 2021; Pharr et al. 2023; Sobol’ 1967] on the simpler scene (middle row), while in the more complex scene (top row), Sobol’ with Owen scrambling, Faure and Lemieux’ and cascaded Sobol’ perform similarly to our method. In the case of indirect light (bottom row), we perform better than ZeroTwo [Pharr et al. 2023] that involves the first two Sobol’ dimensions, and similarly to the approach of Burley [2020] that involves the first four Sobol’ dimensions (padded Sobol’ (0123)). This underlines the importance of uniformity on consecutive quads. Neither the padded Sobol’ of Burley [2020] nor cascaded Sobol’ [Paulin et al. 2021] guarantee higher dimensional uniformity, although they work remarkably well on this example.

We provide evidence for the existence of (0,2)-progressive pairs of dimensions in base b = 3 when using irreducible polynomials in Sobol’ construction [Faure and Lemieux 2016], that can be combined to form (2,4)-progressive quads. These quads allow for reducing the discrepancy in higher-dimensional problems and effectively reduce integration errors. We provide initialization tables in supplementary materials for moderate dimensions only (48 dimensions, as compared to the 10k dimensions of Joe & Kuo [2008] or 16,510 in the work of Faure and Lemieux [2019]). This limitation is not fundamentally linked to our approach, but by computer power and fine optimization of our unoptimized search implementation. The base-3 construction requires a base-3 Owen scrambling for rendering applications as well as modulo 3 arithmetic, which reduce speed compared to efficient bitwise operations. It also implies that sample counts of highest quality correspond to powers of 3, which grows faster than more commonly used powers of 2. Finally, aside from the first four dimensions, the only other criteria we enforced concerned sample sizes of up to nearly 59k samples. While we do not guarantee the best possible behavior for more than 59,049 samples and outside of the quadruplets of dimensions we optimized, a reasonable uniformity is guaranteed due to Sobol’ construction that enforces the (t, s)-sequence property and Faure and Lemieux's result bounding t by a sum of polynomial degrees [Faure and Lemieux 2016].

Within the context of moderate-dimensional problems where high uniformity in specific pairs of dimensions is desired, such as expected from rendering applications, our sampler should provide benefits for reducing integration error. We extensively tested our sampler based on a discrepancy metric and synthetic integrands, and offered preliminary experiments that exhibit such benefit in practice. We believe our sampler can spur further research in low discrepancy sampling.