Loïs Paulin
email: lois.paulin@ens-lyon.fr

MatBuilder: Mastering Sampling Uniformity Over Projections

Keywords: CCS Concepts:, Computing methodologies → Rendering, • Mathematics of computing → Permutations and combinations, Quadrature, Number-theoretic computations Path tracing, quasi-Monte Carlo integration, low discrepancy sequences, generator matrices, integer linear programming

Many applications ranging from quasi-Monte Carlo integration over optimal control to neural networks benefit from high-dimensional, highly uniform samples. In the case of computer graphics, and more particularly in rendering, despite the need for uniformity, several sub-problems expose a low-dimensional structure. In this context, mastering sampling uniformity over projections while preserving high-dimensional uniformity has been intrinsically challenging. This difficulty may explain the relatively small number of mathematical constructions for such samplers. We propose a novel approach by showing that uniformity constraints can be expressed as an integer linear program that results in a sampler with the desired properties. As it turns out, complex constraints are easy to describe by means of stratification and sequence properties of digital nets. Formalized using generator matrix determinants, our new MatBuilder software solves the set of constraints by iterating the linear integer program solver in a greedy fashion to compute a problem-specific set of generator matrices that can be used as a drop-in replacement in the popular digital net samplers. The samplers created by MatBuilder achieve the uniformity of classic low discrepancy sequences. More importantly, we demonstrate the benefit of the unprecedented versatility of our constraint approach with respect to low-dimensional problem structure for several applications.

INTRODUCTION

Generating samples covering a domain as uniformly as possible is core to many scientific domains including computer graphics. It is well known that simple Monte Carlo integration converges faster when using well-spaced samples than using random sampling, which is particularly useful in physically based rendering [START_REF] Pharr | Physically Based Rendering: From Theory to Implementation[END_REF]. Stratified sampling of a population improves the precision of statistical surveys [START_REF] Kish | Survey Sampling[END_REF]]. Generating highly uniform samples is also important in Deep Learning [START_REF] Keller | Artificial Neural Networks generated by Low Discrepancy Sequences[END_REF] and stratification can be used to prevent class imbalance during learning [START_REF] Yuan | A Regularized Ensemble Framework of Deep Learning for Cancer Detection from Multi-class, Imbalanced Training Data[END_REF]. Qualifying and quantifying this uniformity is a challenge in itself and has been the scope of decades of research, involving tools such as discrepancy [START_REF] Hickernell | A generalized discrepancy and quadrature error bound[END_REF]] and optimal transport [START_REF] Paulin | Sliced optimal transport sampling[END_REF].

Applications may require uniformity, not only in the possibly high-dimensional space of samples, but also for particular lowerdimensional projections. This is notably the case for rendering [START_REF] Perrier | Sequences with Low-Discrepancy Blue-Noise 2-D Projections[END_REF], where integration is performed over a path space consisting of the union of domains of increasing dimensions as light bounces -uniformity is then desired for each integration subdomain. While generating a specific set of samples satisfying complex combinations of constraints may be performed via energy optimization, certain applications require the ability to generate sequences of points instead of a single point set, that is, the ability of adding more samples without touching the existing ones, while preserving uniformity constraints. For instance, rendering applications may progressively show the rendered image being refined or adaptively terminate rendering based on image quality. In addition, the efficiency of sample generation is crucial and it may not be acceptable for point set optimization time to exceed the Monte Carlo simulation itself.

Such constraints are dictated by the user and are application specific. While low-discrepancy sequences have been designed and used as general-purpose samplers -this is the case of the popular family of Sobol' sequences -they may not be adapted to the enduser's requirements.

In this paper, we introduce the mathematical tools and a description language that allow one to define sample point sets and sequences that satisfy complex uniformity constraints. We achieve this feat with two key ideas. First, we we show that these constraints can be written in the form of constraints on matrix determinants. This formalization allows us to represent the constraints as an integer linear program [START_REF] Wolsey | Integer Programming[END_REF]]. Second, we show that higher-order finite fields offer many more degrees of freedom to satisfy constraints, overcoming restrictions of classic algebraic constructions. Solving the integer linear program lets us create application-specific samplers. We demonstrate the superior performance of the resulting samplers for rendering, texture exploration, and optimal control applications. The source code of the proposed MatBuilder system has been released to facilitate further research [START_REF] Paulin | MatBuilder: Mastering Sampling Uniformity over Projections[END_REF]].

In numerical integration, quasi-Monte Carlo methods can outperform Monte Carlo methods by using correlated samples that are more uniformly distributed than independent uniformly distributed random samples. Quasi-Monte Carlo methods have profoundly impacted computer graphics and are now standard in rendering technologies [START_REF] Pharr | Physically Based Rendering: From Theory to Implementation[END_REF]]. We will review the work most relevant to our new construction.

Low-discrepancy sequences. Deterministic number-theoretic constructions of highly uniform point sets and sequences form the foundation of quasi-Monte Carlo methods [START_REF] Halton | Radical-inverse quasi-random point sequence [G5[END_REF][START_REF] Niederreiter | Random Number Generation and quasi-Monte Carlo Methods[END_REF][START_REF] Ilya | On the distribution of points in a cube and the approximate evaluation of integrals[END_REF]. Their uniformity is measured by discrepancy [START_REF] Hickernell | A generalized discrepancy and quadrature error bound[END_REF]] that is defined as the supremum of the absolute difference between the relative number of samples that fall into an arbitrary convex subset of the unit cube and the volume of the subset. Point sequences that achieve a discrepancy of order O (log 𝑁 /𝑁 ) are called low-discrepancy sequences. The classic theory and constructions of low discrepancy sequences are extensively covered in reference text books by [START_REF] Niederreiter | Random Number Generation and quasi-Monte Carlo Methods[END_REF], [START_REF] Lemieux | Monte Carlo and Quasi Monte Carlo Sampling[END_REF], and [START_REF] Dick | Digital Nets and Sequences: Discrepancy Theory and Quasi-Monte Carlo Integration[END_REF].

Our contribution draws from the algorithmic framework for the generation of low-discrepancy sequences and hence can be a drop-in replacement for the improvement of existing codes.

Projective samplers. Already Sobol ' [1967] recognized that there are combinatorial lower bounds on the achievable uniformity of low-dimensional projections of his algebraic high-dimensional construction. Later Sobol ' et al. [2011] introduced the concept of constraints to improve the uniformity of low-dimensional projections at the cost of relaxed high-dimensional uniformity. Similarly, improvements to the low dimensional projections of the Sobol' sequence have been identified by computer search [START_REF] Joe | Constructing Sobol' sequences with better twodimensional projections[END_REF]. More recently, computer search has been key to improving the quality parameter 𝑡 of a generalization of Sobol' construction [Faure andLemieux 2016, 2019].

Recognizing the particular structure of physically-based rendering integration domains, for example obtained by repeated scattering in path tracing algorithms, computer graphics has been at the forefront of improving the uniformity of point sets in low dimensional projections using both computer search [START_REF] Abdalla | Low-Discrepancy Blue Noise Sampling[END_REF][START_REF] Marques | Spectral Analysis of Quadrature Rules and Fourier Truncation-Based Methods Applied to Shading Integrals[END_REF][START_REF] Paulin | Sliced optimal transport sampling[END_REF][START_REF] Perrier | Sequences with Low-Discrepancy Blue-Noise 2-D Projections[END_REF]] and mathematical constructions [START_REF] Paulin | Cascaded Sobol' Sampling[END_REF].

As compared to the aforementioned specific approaches, we come up with a general formulation of constraints and a system to solve them in an efficient way.

Orthogonal Arrays. Orthogonal Arrays (OA) are widely used in engineering and the statistical design of scientific experiments [START_REF] Bose | Orthogonal arrays of strength two and three[END_REF][START_REF] Bush | Orthogonal arrays of index unity[END_REF][START_REF] Plackett | The design of optimum multifactorial experiments[END_REF]. Given a set of values in {0, . . . , 𝑏 -1}, an OA in dimension 𝑠 of strength 𝑡 ≤ 𝑠 is set of 𝑛 points seen as an array 𝐴 of size 𝑛 × 𝑠 where each submatrix of size 𝑛 × 𝑡 (i.e. 𝑡 dimensional projection) contains all 𝑏 𝑡 possible distinct rows [START_REF] Owen | Monte Carlo Theory, Methods and Examples[END_REF]].

As introduced to computer graphics [START_REF] Jarosz | Orthogonal Array Sampling for Monte Carlo Rendering[END_REF], OAs generalize the notion of Latin hypercube and (multi-)jittered sampling, focusing on stratification properties on some 𝑡-dimensional projections of the samples. So-called strong OAs fill the gap between orthogonal arrays and low discrepancy nets [START_REF] He | Strong orthogonal arrays and associated Latin hypercubes for computer experiments[END_REF].

The constraints implied by orthogonal arrays are a special case of our general approach to designing generator matrices.

Randomization. Low-discrepancy sequences may be randomized to allow for variance estimation and unbiased Monte Carlo integration [START_REF] Owen | Monte Carlo extension of Quasi-Monte Carlo[END_REF]]. The underlying principles are recursive random scrambling and random shifting on the unit torus. These techniques may also be used for the optimization of sampling points according to various criteria: besides minimizing discrepancy, the mutual minimum distance of the points may be maximized and the spectral properties may be tweaked to approximate blue noise characteristics in low dimensions [START_REF] Abdalla | Low-Discrepancy Blue Noise Sampling[END_REF]Ahmed andWonka 2020, 2021;[START_REF] Heitz | A Low-Discrepancy Sampler that Distributes Monte Carlo Errors as a Blue Noise in Screen Space[END_REF].

While many improvements of scrambled algebraic samplers are intrinsic to our constraint-based approach, our samplers are compatible with the aforementioned randomization techniques.

PRELIMINARIES

Our new method is based on the classic algorithms (reviewed in Section 3.1) to generate digital nets and sequences [START_REF] Dick | Digital Nets and Sequences: Discrepancy Theory and Quasi-Monte Carlo Integration[END_REF][START_REF] Lemieux | Monte Carlo and Quasi Monte Carlo Sampling[END_REF][START_REF] Niederreiter | Random Number Generation and quasi-Monte Carlo Methods[END_REF] in the 𝑠dimensional unit cube [0, 1) 𝑠 , where a set of generator matrices 𝐶 0 , . . . , 𝐶 𝑠 -1 ∈ F 𝑚×𝑚 p defines an 𝑠-dimensional point set of size p 𝑚 . Rather than constructing such matrices using number theory, we aim to design the generator matrices 𝐶 0 , . . . , 𝐶 𝑠 -1 by specifying constraints that we introduce in Section 4 as design principles derived from the algebraic properties of the classic constructions (see Section 3.2).

Digital Nets and Sequences from Generator Matrices

To compute the 𝑖-th component 𝑥 𝑎,𝑖 of the 𝑎-th point x 𝑎 , the index (or its ordinal number) 𝑎 ∈ N 0 is represented in base p as the vector (𝑎 0 , . . . , 𝑎 𝑚-1 ) of 𝑚 digits, where 𝑎 0 is the least-significant digit. Given the 𝑖-th generator matrix 𝐶 𝑖 , the vector

𝑏 𝑚-1 . . . 𝑏 0 := 𝐶 𝑖 • 𝑎 0 . . . 𝑎 𝑚-1 (1)
is mapped to the unit interval by

𝑥 𝑎,𝑖 := 1 p 𝑚 𝑚-1 ∑ 𝑗=0 𝑏 𝑗 • p 𝑗 ∈ [0, 1) .
(2)

The matrix multiplication is performed in the finite (Galois) field F p (sometimes denoted GF(p)). Using a finite field with a base p that is a prime power rather than a prime would involve mapping the digits 𝑎 𝑘 to the finite field and mapping the digits 𝑏 𝑘 back into the integers [Niederreiter 1992, Section 4.3]. Applying only prime bases in this article, these mappings are not required and operations are simply performed modulo p.

The above algorithm has become extremely popular, as in base p = 2, the vector operations on F 𝑚 2 can be efficiently implemented as bitvector 'xor' and 'and' operations for addition and multiplication, ACM Trans. Graph.,Vol. 41,No. 4 . . .

0 2 𝑚 -1 (a) 𝐶 1 𝐶 0 𝑀 𝐾 (b)
Fig. 1. Generator matrix structure: (𝑎) The first row of the generator matrix 𝐶 𝑖 acts on the most significant digit of 𝑏 𝑚-1 selecting one two intervals of size 2 𝑚-1 in base 2, the first two rows of 𝐶 𝑖 determine one intervals of size 2 𝑚-2 and so on. (𝑏 ) When coupling the matrices to design a 2D sampler, the first rows of 𝐶 0 select intervals along the 𝑥-axis, while the first rows of the 𝐶 1 select intervals along the 𝑦-dimension. If 𝑀 𝐾 is of full rank, samples generated by 𝐶 0 and 𝐶 1 will be stratified in the cells induced by selected rows of the two matrices.

respectively. For higher bases, efficient implementations process multiple digits simultaneously using lookup tables [START_REF] Keller | Advanced (Quasi) Monte Carlo Methods for Image Synthesis[END_REF]].

Stratification of Digital Nets and Sequences

Digital nets and sequences generated by matrices imply a stratification of the unit cube as illustrated in Figure 3. These strata are specified by Definition 3.1 (see [Niederreiter 1992, p.48]). For a fixed dimension 𝑠 ≥ 1 and an integer base 𝑏 ≥ 2 the subinterval

𝐸 = 𝑠 𝑗=1 𝑎 𝑗 • 𝑏 -𝑑 𝑗 , (𝑎 𝑗 + 1) • 𝑏 -𝑑 𝑗 ⊆ [0, 1) 𝑠 with 0 ≤ 𝑎 𝑗 < 𝑏 𝑑 𝑗 , 𝑎 𝑗 , 𝑑 𝑗 ∈ N 0 , is called an elementary interval in base p.
The stratification and hence uniformity constraints given by all the elementary intervals can be used to characterize the uniformity of point sets and sequences: Definition 3.2 (see [Niederreiter 1992, Def. 4.1]). For integers 0 ≤ 𝑡 ≤ 𝑚, a (𝑡, 𝑚, 𝑠)-net in base 𝑏 is a point set of 𝑏 𝑚 points in [0, 1) 𝑠 such that there are exactly 𝑏 𝑡 points in each 𝑏-adic elementary interval 𝐸 with volume 𝑏 𝑡 -𝑚 .

Definition 3.3 (see [Niederreiter 1992, Def. 4.2]). For an integer 𝑡 ≥ 0, a sequence x 0 , x 1 , . . . of points in [0, 1) 𝑠 is a (𝑡, 𝑠)-sequence in base 𝑏 if, for all integers 𝑘 ≥ 0 and 𝑚 > 𝑡, the point set x 𝑘𝑏 𝑚 , . . . , x (𝑘+1)𝑏 𝑚 -1 is a (𝑡, 𝑚, 𝑠)-net in base 𝑏.

Revisiting the point set in Figure 3, we see a realization of a (0, 4, 2)-net in base p = 2 that is an 𝑠 = 2-dimensional point set with 2 4 = 16 points, where 𝑡 = 0, i.e. exactly one point is in each possible elementary interval. A smaller quality parameter 𝑡 implies a lower discrepancy and hence better uniformity of a point set or sequence.

Historically, Sobol's construction [START_REF] Ilya | On the distribution of points in a cube and the approximate evaluation of integrals[END_REF] has led to the above definitions and in fact, the Sobol' low discrepancy sequence is a (𝑡, 𝑠)-sequence in base 2. Its generator matrices 𝐶 𝑖 are constructed from the 𝑖-th primitive polynomial over F 2 . [START_REF] Joe | Constructing Sobol' sequences with better twodimensional projections[END_REF] have published Sobol' generator matrices with optimized low-dimensional projections up to dimension 𝑠 = 21201 that are widely used in computer graphics and numerical simulation in general. [START_REF] Grünschloß | t, m, s)-Nets and Maximized Minimum Distance[END_REF] have performed an exhaustive search on the binary matrix 𝐶 1 such that starting from the identity matrix 𝐶 0 = 𝐼 , the point sets generated by (𝐶 0 , 𝐶 1 ) are (0, 𝑚, 2)-nets maximizing the minimal distance between the samples. More recently, [START_REF] Paulin | Cascaded Sobol' Sampling[END_REF] have proposed a cascaded construction of the 𝐶 𝑖 matrices in base 2 in a way that any pair of consecutive matrices (𝐶 𝑖 , 𝐶 𝑖+1 ) defines a (0, 𝑚, 2)-net in base 2, with optimized discrepancy on the remaining two-dimensional projections. Ahmed and Wonka [2021] further explore all possible dyadic nets in base 2 but only for 2 dimensions. Still in base 2, L'Ecuyer and Munger [2016] propose a generic tool, LatNetBuilder, for constructing matrices in higher dimensions optimizing a given figure of merit (e.g. full-space or projective discrepancy), using either an exhaustive or a randomized exploration.

There remain severe limitations to a low base. For example, achieving the best quality parameter 𝑡 = 0 is impossible for 𝑠 > p. While there exist constructions with 𝑡 = 0 in higher bases, a higher base results in a higher discrepancy and hence reduced uniformity. In general, it is considered hard to construct algebraic low-discrepancy sequences. Since we aim to improve the uniformity of low-dimensional projections, too, the problem becomes even harder. Therefore, we propose a novel, constraint-based approach to digital nets and sequences.

GENERATOR MATRIX CONSTRAINTS

In this section we show how classic uniformity constraints translate into determinant inequalities that can be used as design constraints.

How Generator Matrices Act

The algorithm of Equation ( 1) is illustrated in Figure 1: the 𝑗-th row of a generator matrix 𝐶 𝑖 determines the 𝑗-th digit 𝑏 𝑚-𝑗 . The first 𝑘 rows of 𝐶 𝑖 determine the 𝑘 most significant digits of the vector 𝑏. They in turn determine in which of the p 𝑘 uniformly sized intervals 𝑙 • p 𝑚-𝑘 , (𝑙 + 1) • p 𝑚-𝑘 , 0 ≤ 𝑙 < p 𝑘 the sample will be placed.

If the generator matrix is of full rank, that is its determinant is non-zero, the matrix multiplication in Equation (1) will be bijective and thus corresponds to permuting the set {0, . . . , p 𝑚 }. Hence a set of generator matrices with non-zero determinants generates a point set by Equation ( 2) that is perfectly stratified in each single dimension 𝑖 and thus a Latin hypercube sample. The most prominent example of this property is the Sobol' sequence.

As a generator matrix with non-zero determinant implies perfect stratification, we conclude that non-zero determinants may be used to express uniformity constraints.

Design Constraint: Stratification and Net Properties

The intuition developed in Section 4.1 may be used to define uniformity constraints across dimensions. Let 𝐾 = (𝑘 0 , . . . , 𝑘 𝑠 -1 ) ∈ N 𝑠 be a set of non-negative integers such that 𝑠 -1 𝑖=0 𝑘 𝑖 = 𝑚. Then, for each 𝑘 𝑖 the first rows of 𝐶 𝑖 map the index 𝑎 to one of the intervals

𝑙 • p 𝑚-𝑘 𝑖 , (𝑙 + 1) • p 𝑚-𝑘 𝑖 , 0 ≤ 𝑙 < p 𝑘 𝑖 .
Compositing all 𝑘 𝑖 first rows of the corresponding generator matrix 𝐶 𝑖 results in an 𝑚 × 𝑚 matrix 𝑀 𝐾 that maps the index 𝑎 to one of p 𝑚 elementary intervals of size p 𝑚-𝑘 0 × • • • × p 𝑚-𝑘 𝑠 -1 as illustrated for two dimensions in Figure 1.

By analogy with Section 4.1, a composite matrix 𝑀 𝐾 with a nonzero determinant is a bijection between the index 𝑎 and the elementary intervals specified by 𝐾. In case the dimension 𝑠 divides the number 𝑚 of digits, the generator matrices 𝐶 0 , . . . , 𝐶 𝑠 -1 generate a stratified point set if the determinant det 𝑀 (𝑚/𝑠,...,𝑚/𝑠 ) ≠ 0 , because each stratum of volume 𝑝 𝑚/𝑠 × • • • × 𝑝 𝑚/𝑠 contains exactly one point. In case 𝑠 does not divide 𝑚, the property can be generalized to elementary intervals with a ratio of at most p between the longest and shortest dimensions. This leads to Design Constraint 1 (Generalized Stratification). A point set of size p 𝑚 is said to be stratified if

∀𝐾 = (𝑘 0 , . . . , 𝑘 𝑠 -1 ) ∈ 𝑚 𝑠 , 𝑚 𝑠 𝑠 with 𝑠 -1 ∑ 𝑖=0 𝑘 𝑖 = 𝑚 , all cells of size p 𝑚-𝑘 0 × • • • × p 𝑚-𝑘 𝑠 -1 contain exactly 1 point. Thus 𝐶 0 , . . . , 𝐶 𝑠 -1 generate a stratified point set if and only if det (𝑀 𝐾 ) ≠ 0 . (3)
Classic stratification along the canonical axes is a special case of the more general concept of stratification by elementary intervals. The constraint resulting from the definition of (0, 𝑚, 𝑠)-nets [Niederreiter 1992] is:

Design Constraint 2 ((0,m,s)-net). An 𝑠-dimensional point set of size p 𝑚 is said to be a (0, 𝑚, 𝑠)-net if ∀𝐾 = (𝑘 0 , . . . , 𝑘 𝑠 -1 ) ∈ N 𝑠 with 𝑠 -1 ∑ 𝑖=0 𝑘 𝑖 = 𝑚 , all cells of size p 𝑚-𝑘 0 × • • • × p 𝑚-𝑘 𝑠 -1 contain exactly one point. Thus 𝐶 0 , . . . , 𝐶 𝑠 -1 generate a (0, 𝑚, 𝑠)-net if and only if det (𝑀 𝐾 ) ≠ 0 . (4)
We conclude that stratification and (0, 𝑚, 𝑠)-net constraints of digital nets and sequences can be verified by computing determinants of matrices that are composited from sets of first rows of generator matrices {𝐶 0 , . . . , 𝐶 𝑠 -1 } in F p .

Design Constraint: Sequences

Many applications use adaptive sampling, where the number of samples is controlled by a termination criterion, e.g. a threshold on the variance in Monte Carlo integration. Progressive random sampling is straightforward, as it just requires to draw more random samples. In the case of our digital construction, points have to respect uniformity constraints that are only satisfied for sample counts of powers of p, see Section 3.2.

In theory a sequence property shall be valid for an infinite number of samples. In practice the number of samples is bounded due to finite floating point precision and finite computation time. Given the base p, this allows us to select a suitable 𝑚 determining an upper bound of p 𝑚 samples that covers the range of samples representable in floating point arithmetic. As a consequence, one can select a point set of maximum size p 𝑚 and consider it to be sequential if, for 1 ≤ 𝑖 ≤ 𝑚, all subsets of the first p 𝑖 points satisfy the Design Constraint 3 (Seqential property). 𝐶 0 , . . . , 𝐶 𝑠 -1 have a sequence property if the property is valid for all 𝑗 × 𝑗 submatrices 𝐶 𝑗 0 , . . . , 𝐶 𝑗 𝑠 -1 of 𝐶 0 , . . . , 𝐶 𝑠 -1 anchored at the most significant digit, 1 ≤ 𝑗 ≤ 𝑚.

Design Constraint: The Prime Base p

The base p in which one decomposes the index 𝑎 in Equation ( 1) is an important choice to consider in the design of generator matrices. It is a known result [START_REF] Niederreiter | Random Number Generation and quasi-Monte Carlo Methods[END_REF]] that in base p it is neither possible to construct a (0,m,s)-net for 𝑠 > p + 1 nor a (0,s)-sequence for 𝑠 > p. The base may also have an impact on the satisfiability of the entire set of other constraints. It determines the number of values any specific matrix entry can take. In turn, each constraint det(•) ≠ 0 may reduce the range of possible matrix entries.

As an example, given a matrix 𝐶 0 , we search for a matrix 𝐶 1 such that both generate a (0, 2)-sequence in base p = 2. By induction, one can show that there is only one possible value for each entry of 𝐶 1 and thus there exists exactly one 𝐶 1 that results in a (0, 2)-sequence given 𝐶 0 . In particular, if 𝐶 0 is the identity matrix, 𝐶 1 corresponds to the Pascal triangle, and (𝐶 0 , 𝐶 1 ) forms the first two dimensions of the Sobol sequence. This restriction due to the small base can be severely limiting. For instance, we show in Figure 2 the samples obtained by generalizing the construction of Cascaded Sobol' Sampling [START_REF] Paulin | Cascaded Sobol' Sampling[END_REF] where all consecutive pairs of dimensions form (0, 𝑚, 2)nets, to instead impose the constraint that all consecutive pairs of dimensions form (0, 2)-sequences, using base p = 2. Although consecutive pairs of this construction are (0, 2)-sequences, half of the 2-dimensional projections have extremely poor uniformity that can be detrimental in many applications such as rendering. This is due to the fact that since satisfying the (0, 2)-sequence constraint has only one solution in base p = 2, all matrices with even indices equal 𝐶 0 and others equal 𝐶 1 . Selecting a higher base comes at the price of a larger distance between sample counts at which design constraints derived in the previous sections can be verified. We further discuss changing the base in Section 6.2.

SOLVING THE CONSTRAINTS

Given a problem dimension 𝑠, a base p, and set of uniformity constraints as discussed in the previous section, the goal is to determine 𝑠 generator matrices 𝐶 𝑖 with entries in {0, . . . , p -1} that satisfy the constraints. In this section, we prove that these constraints are equivalent to an integer linear program [START_REF] Wolsey | Integer Programming[END_REF]], whose solution -if existent -yields a set of generator matrices with the desired uniformity properties.

Our approach is motivated by Definition 3.3 stating that (𝑡, 𝑠)sequences are sequences of (𝑡, 𝑚, 𝑠)-nets for 𝑚 > 𝑡. We hence implement matrix construction as a greedy algorithm that iterates over the number of digits 𝑚 starting with 1 × 1 matrices until the target size is reached. To grow the matrices, we only need to consider triangular matrices, see the proof in Appendix A: Theorem 5.1. If a set of constraints can be satisfied by square matrices, it will equivalently be satisfiable by triangular matrices up to an Owen scrambling step.

All the aforementioned uniformity constraints have the common form of det(𝑀) ≠ 0 .

(5) Using the recursive determinant formula yields a polynomial of degree 𝑚 with the 𝑚 2 matrix entries as variables in F p for each constraint. Already solving for polynomials on integers is a NPhard problem with few efficient approximations [START_REF] Belotti | Mixed-integer nonlinear optimization[END_REF]]. However, we can transform the Inequality (5) into a set of linear constraints to be solved in Z instead of F p : For p a prime number, the first linear constraint amounts to

0 < det(𝑀) -𝑘p < p (6)
with 𝑘 ∈ Z an additional variable representing the modulo arithmetic of F p . As we no longer operate in F p , the matrix elements need to be constrained to {0, . . . ,

p -1} ⊂ Z that is ∀0 ≤ 𝑗 ≤ 𝑚 -1, 0 ≤ 𝑀 𝑗,𝑚-1 < p. ( 7
)
Note that since we want to iteratively grow the matrices, the 𝑚 -1 × 𝑚 -1 triangular matrix is already known. Hence only the newly added last column of matrix entries 𝑀 𝑗,𝑚-1 adds to the constraints.

The key observation is that these 𝑚+1 constraints form an integer linear program that can be used to determine the next larger set of generator matrices -much in the reverse spirit of the recursive determinant formula. The resulting sequence of linear integer programming problems remains NP-hard but is much easier to solve. We formalize this in Theorem 5.2. For given 𝑚 × 𝑚 generator matrices 𝐶 0 , . . . , 𝐶 𝑠 -1 , adding a column of variables 𝑐 𝑖 𝑙,𝑚+1 to the right and adding a row of 𝑚 zeros at the bottom yields a set of (𝑚 + 1) × (𝑚 + 1) matrices 𝐶 ′ 0 , . . . , 𝐶 ′ 𝑠 -1 , i.e.

𝐶 ′ 𝑖 = 𝑐 𝑖 1,𝑚+1 𝐶 𝑖 𝑐 𝑖 2,𝑚+1 . . . 0 • • • 0 𝑐 𝑖 𝑚+1,𝑚+1
, satisfying a set of properties that is an integer linear program.

Proof. Let 𝐿 be the set of constraints we wish to satisfy. Each constraint is in the form of det(𝑀) ≠ 0 for a given matrix 𝑀. According to our design principles 𝑀 is an (𝑚 + 1) × (𝑚 + 1) matrix composed of lines from multiple 𝐶 ′ 𝑗 . Thus all cells are known except for the last column. Since the determinant is linear by columns, it is linear for all variables 𝑐 𝑖 𝑙,𝑚+1 in the new column, leading to the fact that Equation ( 6) and ( 7) correspond to a pair of integer linear constraints. Packing all determinants together, 𝐿 defines an integer linear program. □

When the constraints are transformed into an integer linear program they can be prioritized by assigning a weight 𝑤 𝑗 . This is implemented by extending the hard constraint of Equation ( 6) to

𝑥 𝑗 ≤ det(𝑀) -𝑘p < p, where 𝑥 𝑖 ∈ {0, 1} , (8)
and adding 𝑤 𝑗 𝑥 𝑗 to the objective function to be maximized. Unless det(𝑀) = 0, in which case the constraint cannot be satisfied and hence 𝑥 𝑗 = 0 must follow, 𝑥 𝑗 will be set to 1 and 𝑤 𝑗 𝑥 𝑗 will count towards the objective function.

u0 u4 u2 Det = 1 Det = 1 Det = -1 Det = 1 Det = -1 𝑀 𝐾 = 𝐶 0 = 𝐶 1
Fig. 3. Modifiers u<m'> illustrated for the relaxed (0, 𝑚, 𝑠 )-net constraint in 2 dimensions, base 2, and for 2 4 samples: u0 corresponds to a stratification test enforcing one sample per 1 4 × 1 4 strata, u4 amounts to a complete (0, 𝑚, 2)-net check, and u2 leads to intermediate constraints. Such properties can be assessed by checking the determinant of the 𝑀 𝐾 matrices: 2 rows from 𝐶 0 and 𝐶 1 (u0), 1 row from 𝐶 0 , three from 𝐶 1 and conversly (u2), or the full 𝐶 0 and 𝐶 1 matrices (u4). The strata depicted in the unit squares on the right are the elementary intervals underlying the displayed (0, 4, 2)-net.

In summary, we implement matrix construction as a greedy algorithm that iterates over 𝑚 and digit by digit determines the columns of all 𝐶 𝑖 by solving the set of linear constraints derived in Equations ( 6) and ( 7) while maximizing a target objective function that is the weighted sum of satisfied weak constraints (see Figure 4). As solving this problem remains NP-hard, we rely on efficient heuristics as implemented in the CPLEX [IBM 2022] library in order to find instances that satisfy the constraints.

With our greedy approach it may happen, especially for restrictive sets of constraints, that choices made for a previous column may render finding the next column unsatisfiable. To tackle this issue we add a backtracking feature that when faced with an unsatisfiable set of constraints drops the most recently added column and tries again with different choices. When this strategy fails again, the algorithm starts from scratch with empty matrices. This behavior can be observed in the timings of the mixed profile in Figure 5.

If even after backtracking the linear integer programs remain infeasible, our algorithm will signal that the chosen set of constraints does not admit a solution. In this case, the user needs to consider a higher prime base or turning some of the hard constraints into weak ones. We also provide a timeout parameter that limits the maximum time spent on maximizing the number of satisfied weak constraints.

Uniformity beyond Generator Matrices. Up to this point both constraints and solver are deterministic. In order to optimize with respect to uniformity metrics that cannot be expressed by generator matrices alone, we extend the objective function by a stochastic 𝑖 are determined by solving an ILP, which consists of an objective function to maximize subject to a set of constraints. The range constraints enforce that 𝑐 𝑖 𝑙,𝑚+1 ∈ {0, . . . , p -1} and the hard uniformity constraints enforce a non-zero determinant to guarantee the design constraints of stratification, net, and sequence properties as introduced in Section 4. Remember that matrices 𝑀 are constructed from first rows of the 𝐶 ′ 𝑖 matrices and hence include some of the 𝑐 𝑖 𝑙,𝑚+1 . Indicated by 𝑥 𝑗 = 1, a satisfied weak constraint adds its weight 𝑤 𝑗 to the objective function. Otherwise, a zero determinant comes along with 𝑥 𝑗 = 0. The stochastic regularization term is detailed in Equation ( 9).

𝑤 𝑗 • 𝑠 -1 𝑖=0 𝑚+1 𝑙=1 𝑐 𝑖 𝑙,𝑚+1 -𝑟 𝑖 𝑙 𝑠 • 𝑚 • (p -1) (9)
Besides maximizing the objective function by satisfying a maximum number of the 𝑉 weak constraints, we generate random numbers 𝑟 𝑖 𝑙 ∈ {0, . . . , p -1} and subtract the normalized ℓ 1 -distance between them and the matrix elements of the new columns. The normalization accounts for the maximum distance, which is the product of the largest possible number p -1 and the number 𝑠 •𝑚 of variables to be determined by Theorem 5.2. In order not to override the weak constraints, the normalized difference is scaled by the smallest weight of the weak constraints. As a result, the solution of the ILP will be as close as possible to the random choice 𝑟 𝑖 𝑙 , which allows us to sample a large diversity of solutions satisfying all strong constraints, while optimizing the weak ones.

Given such a set of solutions, any application specific selection may be exercised. Options are uniformity metrics that cannot be expressed as generator matrix properties such as maximized minimum distance [START_REF] Grünschloß | t, m, s)-Nets and Maximized Minimum Distance[END_REF]], spectral analysis [START_REF] Pilleboue | Variance analysis for Monte Carlo integration[END_REF], or Optimal Transport [START_REF] Paulin | Sliced optimal transport sampling[END_REF].

For the example of our rendering application in Section 7.1, we generate 64 sets of matrices and select the best candidate with respect to lowest 𝐿 2 discrepancy. Usually, the quality of the generated matrices exposes low variance as can be seen Figure 6-top.

THE MATBUILDER SOLVER

The input to the MatBuilder program is a text file that specifies the constraints imposed on the desired point set as introduced in the next section. The solver returns the 𝑠 generator matrices 𝐶 𝑖 ∈ F 𝑚×𝑚 p that are used to generate up to p 𝑚 samples in [0, 1) 𝑠 according to Equation (2). The performance of the algorithm is explored in Section 6.2.

Constraint Specification Language

The header of the constraint text file specifies the dimension 𝑠, the finite field F p , and the number of uniformly distributed points p 𝑚 to generate: s= < dimension s of point set > p= < base p > m= < matrix dimension m > Hard constraints on selected dimensions are added using the keywords net <i_1 > <i_2 > ... <i_s ' > # (0 ,m ,s ') -net constraint where the dimension 𝑠 ′ amounts to the number of dimensions enumerated after the keyword and stratified <i_1 > <i_2 > ... <i_s '> By default, the net and stratified constraints are enforced for submatrices, i.e. all 𝑚 ′ < 𝑚, which allows for progressive sampling with the selected constraints guaranteed. As for example, it does not make sense to check the stratification for less than p 2 points in dimension 𝑠 = 2, using the from and to keywords, constraints can be restricted to sample indices between p 𝑚 1 and p 𝑚 2 : from <m_1 > to <m_2 > net <i_1 > <i_2 > ... <i_s ' > from <m ' > stratified <i_1 > <i_2 > ... <i_s ' >

We have added a modifier to the net keyword, denoted u<m'> for 0 ≤ 𝑚 ′ < 𝑚 to define relaxed versions of the net property: we only keep constraints that rely on an 𝑀 𝐾 matrix (see design constraint 2) with a maximum difference of 𝑚 ′ between the number of rows taken from the 𝐶 𝑖 . This way u0 corresponds to the Property A of [START_REF] Ilya | On the distribution of points in a cube and the approximate evaluation of integrals[END_REF] (contiguous blocks of p 𝑠 samples are stratified), u1 corresponds to design constraint 1, u<m> would correspond to the original net constraint (see Figure 3). This modifier allows for a control of both the positively correlated uniformity strength and computational cost of the problem. The efficiency of different values 𝑚 ′ with respect to discrepancy is shown in Figure 6.

When hard constraints cannot be satisfied, the weak keyword allows one to attach a weight resulting in a weak constraint: weak <w > from <m_1 > to <m_2 > net <i_1 > <i_2 > ... <i_s '> weak <w > from <m '> stratified <i_1 > <i_2 > ... <i_s '>

Numerical Evidence

To evaluate the performances of the MatBuilder program, we consider several generic scenarios: first, a low discrepancy sequence in base 3 with up to 3 10 samples in dimension 8, with as low discrepancy as possible for the full dimension (net):

# Generic -full -space -LDS s= 8 p= 3 m= 10 weak 1 net 0 1 2 3 4 5 6 7

Solutions to this profile will be as close as possible to (0, 8)sequences up to 3 10 samples (note that exact (0, 𝑠)-sequences do not exist for 𝑠 > p in base p [START_REF] Niederreiter | Random Number Generation and quasi-Monte Carlo Methods[END_REF]).

Second, we consider a perfect projective profile in dimension 6, where each consecutive pair of dimensions is a (0, 2)-sequences up to 3 10 in base 3, and all remaining pairs must be as close as possible to a generalized stratification (see Figure 7 for a comparison to the Sobol' sequence, and the net obtained by LatNetBuilder [L'Ecuyer and Munger 2016], when optimizing the same pairs of projection):

# Generic -proj -LDS s= 6 p= 3 m= 10 net 0 1 net 1 2 net 2 3 net 3 4 net 4 5 weak 1 net 0 2 weak 1 net 0 3 weak 1 net 0 4 weak 1 net 0 5 weak 1 net 1 3 weak 1 net 1 4 weak 1 net 1 5 weak 1 net 2 4 weak 1 net 2 5 weak 1 net 3 5 Third, we take a look at an Orthogonal Array profile in dimension 9: we search for a sequential sampler in base 3 with OA strength 3 (each triplet of dimension are stratified starting from 3 3 samples), and a final weak stratification constraint in 9 dimensions: # Generic -OA s= 9 p= 3 m= 10 from 3 stratified 0 1 2 from 3 stratified 1 2 3 from 3 stratified 2 3 4 from 3 stratified 3 4 5 from 3 stratified 4 5 6 from 3 stratified 5 6 7 from 3 stratified 6 7 8 from 5 weak 1 stratified 0 1 2 3 4 5 6 7 8

Finally, a more complex scenario combines stratified and net properties in dimension 10 for p = 3: # Mixed s= 10 p= 3 m= 10 net 0 1 net 1 2 net 2 3 net 3 4 net 4 5 from 3 stratified 0 1 2 from 4 stratified 1 2 3 ACM Trans. Graph., Vol. 41, No. 4, Article 84. Publication date: July 2022. from 3 stratified 2 3 4 from 4 stratified 3 4 5 from 4 to 6 stratified 0 1 2 3 4 from 4 to 6 stratified 1 2 3 4 5

In Figure 5, we evaluate the performance of the solver for all profiles as 𝑚 increases ( 2x Intel(R) Xeon(R) CPU E5-2650 v2 @ 2.60GHz, 32 cores, 64Gb memory, the solver has been given a timeout of 120 seconds per solution). In less than 13 minutes, the solver outputs the generator matrices for all profiles. Note that for the mixed problem, the solver has to backtrack many times for 𝑚 ≤ 4 before being able to complete the matrices. For the generic-OA and the mixed profiles, the largest the integer linear programs to solve has 321 rows, 371 columns, and 1426 nonzero entries (resp. 117 × 228 with 713 nonzero entries), which is a relatively small problem. The generic-full-space-LDS profile induces larger problems as 𝑚 increases (see Figure 5). The problem size is not strictly monotonous in 𝑚, but increases with the matrix size. For this profile, 12871 out of 19448 weak constraints have been satisfied. For the generic-OA, all weak constraints (9/9) have been satisfied. Note that once matrices have been produced our sampler performs with the same efficiency as Halton or Sobol' sequences as they rely on the same kind of matrix vector multiplication.

To evaluate the quality of the MatBuilder outputs, we compare the samples generated by the generic-LDS profile to a random uniform sampler and a Sobol' sampler [START_REF] Ilya | On the distribution of points in a cube and the approximate evaluation of integrals[END_REF] in base 2, using the 𝐿 2 discrepancy. We have considered 64 realizations of the samplers, i.e. 64 runs of the solver for the generic-full-space-LDS profile). As illustrated in Figure 6-(top), we generate samples with a discrepancy comparable to the classic Sobol' sequence with Owen scrambling, while preserving the (0, 2)-sequence property on each pair of dimensions (contrary to Sobol'). Note that our solver generates randomized solutions with very low variance with respect to the 𝐿 2 discrepancy. In Figure 6-(bottom), we evaluate the impact of the u<m'> modifier on the net property: we have considered increasing modifier values {u0, u1, u2, u3, u4, u5, u6 ,u7} and evaluate the associated discrepancy. When relaxing the constraints using a u<m'> modifier, problems are easier to solve (see the timings in Figure 5), but it comes with a trade-off between the efficiency and the quality of the produced generator matrices.

In the supplementary material (Section 1), we further evaluate the impact of timeout constraints and overall quality of the matrices. We also experiment with various prime bases for the Generic-full-space-LDS profile showing both the efficiency of the solver and the uniformity quality of the samples (Section 2).

APPLICATIONS

We demonstrate the advantages of samplers generated from application specific constraints to master uniformity over projections for path tracing, texture exploration, and optimal control.

Path Tracing

In rendering, the path space exposes a clear nested structure: each high-dimensional light transport path is built from low-dimensional events: sampling a pixel footprint (2D), sampling the lens aperture (2D), sampling the time (1D), sampling a point on a light source For the mixed profile, the plots show several backtracks for small 𝑚 before being able to expand the matrices (u4 and u7 have the same timings). At the bottom, the size of the reduced integer linear programs are shown during execution (number of rows and number of nonzero entries) for the generic-full-space-LDS profile (most difficult one in our tests). For this profile, up to 60 calls to the internal integer linear program simplification procedure have been performed.

(3D, usually 1D to select a light source and another 2D to sample a point on this light source), sampling a direction scattered from a material (3D, usually 1D to select a layer and another 2D to sample a direction from this layer), and so on. Previous work [START_REF] Kollig | Efficient Multidimensional Sampling[END_REF][START_REF] Schlick | An Adaptive Sampling Technique for Multidimensional Integration by Ray-Tracing[END_REF]] has composited high-dimensional path space samples from low-dimensional samplers connected by some form of randomization.

We design a path-tracing profile PT-Profile to leverage perfect low dimensional stratification and good stratification over increasing dimensions to replicate the path space structure: Generic full space LDS u1 Generic full space LDS u2 Generic full space LDS u3 Generic full space LDS u4 Generic full space LDS u5 Generic full space LDS u6 Generic full space LDS u7 Fig. 6. Discrepancy test for the generic-full-space-LDS profile: At the top, we compare the generalized 𝐿 2 discrepancy the solver outputs to the Sobol' sampler in base 2 as well as random sampling, jittered sampling, a Rank-1 lattice [START_REF] Keller | Stratification by rank-1 lattices[END_REF]] with generators obtained from [L'Ecuyer and Munger 2016], and Orthogonal Arrays (CMJND) [START_REF] Jarosz | Orthogonal Array Sampling for Monte Carlo Rendering[END_REF]. At the bottom, we evaluate the impact of the u<m'> modifiers on the discrepancy: the higher u, the better the discrepancy as compared to Sobol', of course at the cost of higher timings, see Figure 5).

weak 1 net u4 1 2 3 weak 1 net u4 2 3 4 weak 1 net u4 3 4 5 weak 1 net u2 0 1 2 3 weak 1 net u2 1 2 3 4 weak 1 net u2 2 3 4 5 weak 1 net u2 0 1 2 3 4 weak 1 net u2 1 2 3 4 5 weak 1 net u4 0 1 2 3 4 5

We compare our profile to a commonly used randomization of Sobol'01 (ZeroTwo [START_REF] Pharr | Physically Based Rendering: From Theory to Implementation[END_REF]) with perfect 2D dyadic stratification, to the Sobol' sequence with Owen scrambling, and Cascaded Sobol' points [START_REF] Paulin | Cascaded Sobol' Sampling[END_REF]. PT-Profile works well for smooth components such as depth of field, motion blur, and direct lighting as showcased in Figure 8. A base 𝑝 = 3 sampler behaves differently when selecting 2 light sources or 3, as an odd number of samples cannot be split equally into two. These differences are further illustrated in the first three rows of the same Figure 8 with 2, 3, and 5 light sources. Two light sources have benefit from common base 2 samplers, whereas 5 light sources are a more difficult case for all samplers. The last row depicts a more typical scene setup with 3 light sources.

The performance of the samplers generated by MatBuilder is quite remarkable: It widely matches and outperforms the classic constructions, although it is automatically generated by constraints only.

Correlating the base with the number of light sources leads to an interesting improvement in quality as shown in figure 8. This stems from the fact that the partitioning of points that happens when selecting a light source happens along the border of strata from the net property and thus guarantees the high uniformity of samples used on each different light source. When there is no way to reach such a correlation, the remedy is to use one area sampler per light source in combination with one sampler constraint to the paths [START_REF] Kollig | Efficient Multidimensional Sampling[END_REF].

The supplementary material includes many more renderings and more comparisons.

Exploring Parametric Materials

In this scenario, the objective is to sample the space of a parametric material with explicit correlation between some parameters. This is an initial step in (semi-)procedural texture exploration or texture space analysis [START_REF] Guehl | Semi-Procedural Textures Using Point Process Texture Basis Functions[END_REF]Lasram et al. 2012a,b]. As an example, we consider a procedural texture with seven dimensions to define brick pattern, window brace, and lintel. The structure of the parameters leads to the following profile: s= 7 p= 2 m= 5 # ( Brick -Amount -X , Brick -Amount -Y ) net 0 1 # ( Window -Brace -Amount -X , Window -Brace -Width -X , # Window -Brace -Amount -Y , Window -Brace -Width -Y) weak 1 net 3 4 5 6 # Overall uniformity , including Brick -Lintel -Width weak 2 net 0 1 2 3 4 5 6

The constraints ensure the (0, 𝑚, 2)-net stratification of the first pair of dimensions, a weak 4-dimensional net property of the last four dimensions, and an overall net property across all seven dimensions. Figure 9 illustrates the first 16 realizations as compared to a uniform random sampling of each parameter, exhibiting a more uniform exploration of the parameter space. For this profile, matrices are generated in 3.7 seconds by our MatBuilder software.

Optimal Control

Reinforcement learning can be used to control a robot, a virtual character, or a dynamic system. This amounts to discretizing the states of the system as well as possible actions applied to it, defining a reward for achieving particular states, and iteratively solving Bellman equation through value iterations. This results in a value function and optimal control policy at each state, which can then be interpolated during the simulation in real-time to control the system. We apply the technique described by [START_REF] Coros | Robust task-based control policies for physics-based characters[END_REF] to control a simple inverted pendulum on a chart to reach an upward orientation, a specific position, and a near zero velocity and angular velocity. The state space is 4-dimensional (position 𝑥, velocity 𝑥, angle 𝜃 , dim 0 1 2 3 4

1 2 3 4 5 (a) Sobol', F 2 (b) LatNetBuilder, F 2 (c) Ours, F 3 • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • Ours • Sobol' • LatNetBuilder • Random • log nb
• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • Ours • Sobol'
• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •
Ours Sobol' LatNetBuilder Random 34 35 36 37 38 39 310

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • Ours • Sobol'
• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •
Ours Sobol' LatNetBuilder Random 34 35 36 37 38 39 310

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •
Ours Sobol' LatNetBuilder Random 34 35 36 37 38 39 310

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • Ours • Sobol'
• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •
Ours Sobol' LatNetBuilder Random 34 35 36 37 38 39 310

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •
Ours Sobol' LatNetBuilder Random 34 35 36 37 38 39 310

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •
Ours Sobol' LatNetBuilder Random 34 35 36 37 38 39 310

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • Ours • Sobol'
• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •
Ours Sobol' LatNetBuilder Random 34 35 36 37 38 39 310

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •
Ours Sobol' LatNetBuilder Random 34 35 36 37 38 39 310

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •
Ours Sobol' LatNetBuilder Random 34 35 36 37 38 39 310 and angular velocity 𝜃 ), and the action space is 1-dimensional (a horizontal force 𝑢). The action space is uniformly discretized. We compare state-of-the-art samplings of the state space with samples created by MatBuilder with the following profile that seeks 4-d uniformity:

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • Ours • Sobol'
s= 4 p= 3 m= 17 # particular stratifications were skipped due to # lack of solutions from 3 to 5 stratified 0 1 2 3 from 7 to 9 stratified 0 1 2 3 from 11 to 13 stratified 0 1 2 3 from 15 stratified 0 1 2 3 weak 1 net 0 1 2 3

The state-of-the-art samplers we experimented with are random sampling, a Sobol' sequence [START_REF] Joe | Constructing Sobol' sequences with better twodimensional projections[END_REF][START_REF] Ilya | On the distribution of points in a cube and the approximate evaluation of integrals[END_REF], and a rank-1 lattice produced by LatNetBuilder [L' Ecuyer and Munger 2016]. In all cases, we bias these samples so that they follow the same Gaussian distribution around the target state and get a finer discretization near the solution. In practice, we used the inverse cumulative density function of the Gaussian distribution to generate 3 12 = 531, 441 samples (or 2 19 = 524, 288 for Sobol in base 2), plus an additional sample placed at the exact desired state. To compare these discretizations, we simulate 200 trajectories for 2000 time steps each and compute the average number of time steps 𝛼 the pendulum is within an area of non-zero reward (i.e., within the area considered as success), the number of trajectories 𝛽 which succeed at least in one time step to get a non-zero reward, and the average value 𝛾 of the value function (higher values mean that more time is spent near the desired state). For all these metrics, higher values indicate better behavior. For our stratified sampling, 𝛼 = 170, 𝛽 = 200 and 𝛾 = 530 (the reward for achieving success is 1000, and the maximum value of the optimized value function is near 2500). For 200 trajectories with a single realization of a random sampling, 𝛼 = 12, 𝛽 = 141 and 𝛾 = 349. For this setup, the Sobol' sequence yields 𝛼 = 7, 𝛽 = 44, and 𝛾 = 611 while the rank-1 lattice yields 𝛼 = 1, 𝛽 = 20, 𝛾 = 189. For 200 trajectories consisting in one (different) trajectory for 200 random sampling realizations, 𝛼 = 30, 𝛽 = 121 and 𝛾 = 594. Figure 10 shows a (𝑥, 𝜃 ) slice of 3369 samples of our 3 12 samples within a small ( 𝑥, 𝜃 ) interval, as well as 3 representative trajectories. Clearly, the sampler generated by MatBuilder excels.

DISCUSSION AND PERSPECTIVES

Our algorithm is greedy with respect to the size of matrices and uniformity properties consider whole matrix rows at a time (and not just the last column). Consequently, some solution may be exponentially hard to find as they may require a specific succession of choices among a wide set of locally satisfactory ones. It is interesting to search for a solution generating lines of a matrix at a time. Furthermore, rather than using a greedy algorithm iterating integer linear programs, one may opt for an algorithm generating matrices one at a time through a non-linear integer program solver, which may work better in this aspect.

The design of the objective function of our integer linear programs implies a relatively small control over the way how weak sub-properties of net constraints are distributed. The solver may find the most optimal solution to be satisfying many of the constraints for a single group of dimensions while leaving a pair of dimensions entirely correlated. In perspective, this may be solved using a non-linear convex objective that is only be slightly slower to solve.

In our paper we only deal with prime Galois fields, since it relies on modular arithmetic that can easily be translated into linear constraints in N 0 . This technique cannot be performed for powers of prime bases as they define a specific set of arithmetic operators. To extend prime Galois fields to power bases would require a solver that is able to process addition and multiplication tables specific to the individual Galois fields over prime power bases.

In our rendering section we discuss the fact that the base of the generator matrices sometimes have a huge impact on the quality of integration. The reason for this is that our selectors split the point set at separations of strata considered by the net property and we Exploring parametric texture space: the first 16 realizations of the parameter space sampling using a random sampling strategy (first row) and using samples constraint such that the space of variations is explored more uniformly in accordance with the structure of parameters of the texture description (second row). In this example, we can see that the brick pattern is more uniformly sampled using our approach, the Brick-Lintel-Width is sequentially sampled in a low-discrepancy manner (values between adjacent instances are maximized). For the window brace, all possible variations in 4D are more captured using our profile (parametric material courtesy Adobe Substance 3D Asset). thus guarantee that the points sampling each light source have an optimal distribution. That is true for selectors that split points in a power of p subsets of equal size. We indicated a solution of how more attention could be given as to how to use that fact to improve rendering.

It is noteworthy that the results achieved by the novel constraintbased approach already are competitive or outperform classic constructions of low-discrepancy sequences while providing improved sampling uniformity over projections. We are confident that the evolving the constraint-based approach will allow for more exploration and novel approaches to low-discrepancy sequences in the future.

CONCLUSION

Highly uniform point sets can be tailored specific to the structure of an application as a set of constraints that can be transformed into an integer linear program. Solving such an integer linear program yields a construction for samples that in specific applications have been demonstrated to outperform classic constructions. We hence introduced a new and powerful way for the computer-aided design of generator matrices for digital nets and sequences.

Fig. 2 .

 2 Fig. 2. The base p constrains the space of cascaded sequences: Complying with a (0, 𝑚, 2)-net constraint in base p = 2 for all pairs of subsequent dimensions, there exist exactly two 2-dimensional projections.

Fig. 4 .

 4 Fig. 4. Anatomy of an Integer Linear Program (ILP): To grow the matrices 𝐶 𝑖 according to Theorem 5.2, the values of the 𝑐 𝑖 𝑙,𝑚+1 in the last column of their respective 𝐶 ′𝑖 are determined by solving an ILP, which consists of an objective function to maximize subject to a set of constraints. The range constraints enforce that 𝑐 𝑖 𝑙,𝑚+1 ∈ {0, . . . , p -1} and the hard uniformity constraints enforce a non-zero determinant to guarantee the design constraints of stratification, net, and sequence properties as introduced in Section 4. Remember that matrices 𝑀 are constructed from first rows of the 𝐶 ′ 𝑖 matrices and hence include some of the 𝑐 𝑖 𝑙,𝑚+1 . Indicated by 𝑥 𝑗 = 1, a satisfied weak constraint adds its weight 𝑤 𝑗 to the objective function. Otherwise, a zero determinant comes along with 𝑥 𝑗 = 0. The stochastic regularization term is detailed in Equation (9).

Fig. 5 .

 5 Fig.5. Performance: On the top we show the timings of the MatBuilder solver as a function of the matrix size 𝑚. For the mixed profile, the plots show several backtracks for small 𝑚 before being able to expand the matrices (u4 and u7 have the same timings). At the bottom, the size of the reduced integer linear programs are shown during execution (number of rows and number of nonzero entries) for the generic-full-space-LDS profile (most difficult one in our tests). For this profile, up to 60 calls to the internal integer linear program simplification procedure have been performed.

  u4 0 1 2 ACM Trans. Graph., Vol. 41, No. 4, Article 84. Publication date: July 2022.

Fig. 7 .

 7 Fig. 7. Projective LDS profile: (𝑎) Projection of the first 2048 Sobol' samples in dimension 6; (𝑏 ) projections obtained from optimized LatNetBuilder matrices using a projective profile (see Appendix B);(𝑐 ) the first 3 7 = 2187 samples for the matrices obtained with the generic-proj-LDS profile exhibiting better uniformity in the 2D projections. In (𝑑 ), we quantify our benefit through a 2D discrepancy comparison for all projections. The lower is discrepancy, the better is uniformity.

6D

Fig. 8 .

 8 Fig.8. Rendering results: In 6 dimensions, the dimensions are used to sample the lens (2D), the time (1D), and the light sources (3D). For light sources, the first dimension is used to select a light source (2, 3, or 5 light sources for rows 1, 2, and 3), the remaining 2D coordinates sample a point in the selected light quad. In 8 dimensions, we sample the pixel (2D), the direct lighting (3D) and the indirect lighting (3D). For sampling indirect bounce directions, the first dimension is used to select a material component, and the remaining 2D coordinates importance sample a direction. The first column shows the ground-truth images, the second one a per pixel labelling according to the sampler (Random, ZeroTwo, Sobol', Cascaded, and our PT-Profile) with lowest error at 243spp for our sampler in base 3 or 256spp for others. Finally the last column presents the mean squared error (MSE) for various sample counts.

Fig

  Fig.9. Exploring parametric texture space: the first 16 realizations of the parameter space sampling using a random sampling strategy (first row) and using samples constraint such that the space of variations is explored more uniformly in accordance with the structure of parameters of the texture description (second row). In this example, we can see that the brick pattern is more uniformly sampled using our approach, the Brick-Lintel-Width is sequentially sampled in a low-discrepancy manner (values between adjacent instances are maximized). For the window brace, all possible variations in 4D are more captured using our profile (parametric material courtesy Adobe Substance 3D Asset).

Fig. 10 .

 10 Fig.10. Optimal control: We perform value iterations to solve the Bellman equation for the problem of balancing an inverted pendulum on a chart, upward, centered, with low velocity (left). We discretize the 4-dimensional state space (position 𝑥, velocity 𝑥, angle 𝜃 , and angular velocity 𝜃 ) and give rewards only for vertical positions (rectangle in light gray). We display a 2-d slice (𝑥, 𝜃 ) for values of ( 𝑥, 𝜃 ) near the target state and three sample trajectories; the value function is color coded on each sample. Middle: We use a random discretization, biased towards the target (see Section 7.3). Right: We use a stratified sampling generated with MatBuilder, using the same bias. The space is explored much more uniformly.

  , Article 84. Publication date: July 2022.

			𝑎 0		𝑏 𝑚-1	0	1
			𝑎 1		𝑏 𝑚-2
						00	01	10	11
	𝐶 𝑖	•	𝑎 2	=	𝑏 𝑚-3	000 001 010 011 100 101 110 111
			. . .		. . .
			𝑎 𝑚-1		𝑏 0

ACM Trans. Graph., Vol. 41, No. 4, Article 84. Publication date: July 2022. 84:2 • Loïs Paulin, Nicolas Bonneel, David Coeurjolly, Jean-Claude Iehl, Alexander Keller, and Victor Ostromoukhov

ACM Trans. Graph., Vol. 41, No. 4, Article 84. Publication date: July 2022.

ACKNOWLEDGMENTS

This work was partially funded by ANR-16-CE33-0026 (CALiTrOp) and ANR-20-CE45-0025 (MoCaMed).

A PROOF OF THEOREM 5.1

Proof. Owen scrambling [START_REF] Owen | Monte Carlo extension of Quasi-Monte Carlo[END_REF]] is a randomization technique that does not alter the structure of the elementary intervals and hence preserves all properties of a digital net. To do so, it applies a bijective function 𝑓 𝑏 0 ,...,𝑏 𝑖 -1 : F p → F p to digit 𝑏 𝑖 of the vector 𝑏 in Equation (1). Given factors 𝑤 𝑖,𝑗 ∈ F p , bijective functions of the form 𝑓 𝑏 0 ,...,𝑏 𝑖 -1 (𝑏 𝑖 ) = 𝑖 -1 𝑗=0 𝑤 𝑖,𝑗 𝑏 𝑗 + 𝑏 𝑖 are a realization of Owen scrambling.

Let 𝐶 be the generator matrix of size 𝑚 that generated 𝑏 and denote its rows by 𝐶 𝑖 . Row 𝐶 𝑖 encodes the digit 𝑏 𝑖 . Hence one realization of Owen scrambling is the application of the function 𝑓 𝐶 0 ,...,𝐶 𝑖 -1 (𝐶 𝑖 ) = 𝑖 -1 𝑗=0 𝑤 𝑖,𝑗 𝐶 𝑗 + 𝐶 𝑖 that acts on rows the same way as 𝑓 𝑏 0 ,...,𝑏 𝑖 -1 acts on digits.

We apply 𝑚 such linear functions to each row of 𝐶. By definition, Gaussian elimination can be written as such a set of linear transformations, defining the factors 𝑤 𝑖,𝑗 . Thus for any 𝐶 we can find a realization of Owen scrambling that transforms it into an upper triangular matrix 𝐶 ′ . Since all these functions 𝑓 are bijective, their inverse defines one realization of an Owen scrambling that transforms 𝐶 ′ into 𝐶.

It follows that all nets generated by any generator matrices can also be generated by upper triangular matrices with an additional scrambling step in the sense of Owen. It is therefore sufficient to only consider triangular matrices. □

B LATNETBUILDER PROJECTIVE PROFILE

For Figure 7, we have used the following configuration for LatNet-Builder to generate optimized matrices:

latnetbuilder -o projnet --set -type net \ --construction explicit -s 2^16 \ --dimension 6 --norm -type 1 \ --exploration -method random -CBC :100000 \ --figure -of -merit projdep :t -value \ --weights order -dependent :0:1 ,1 \ --multilevel true --combiner sum