Keywords: Quasi-Monte Carlo, scrambling, digital net, gain coefficient, Sobol’ sequence

AMS subject classifications: 65C05, 65D30

1 Introduction↩︎

We consider numerical approximation of \(d\)-dimensional integrals of the form \[\label{eq:target95integral} I(f) := \int_{[0,1)^d} f(\boldsymbol{x})\,d\boldsymbol{x},\tag{1}\] for Riemann-integrable functions \(f:[0,1)^d\to\mathbb{R}\). A common practical estimator is the equally weighted quadrature rule \[\label{eq:qmc-estimator} I(f;P_n) := \frac{1}{n}\sum_{i=0}^{n-1} f(\boldsymbol{x}_i),\tag{2}\] where \(P_n=\{\boldsymbol{x}_0,\dots,\boldsymbol{x}_{n-1}\}\subset[0,1)^d\) is a collection of sample points. When the \(\boldsymbol{x}_i\) are independent uniform samples, 2 reduces to the classical Monte Carlo (MC) estimator which is unbiased and has variance \[\label{eq:MCvar} \operatorname{Var}(I(f;P_n)) = \frac{\operatorname{Var}{f}}{n}, \qquad \text{where } \operatorname{Var}f := I(f^2) - I(f)^2,\tag{3}\] for \(f\in L^2([0,1)^d)\). Despite its robustness, the root mean square error of the MC estimator converges at the relatively slow rate of \(O(n^{-1/2})\), motivating the use of more structured point sets.

Quasi-Monte Carlo (QMC) methods replace independent sampling by carefully constructed low-discrepancy point sets such as lattice rules [1], [2], Halton sequences [3], and digital nets [4]; [5]. For integrands with bounded variation, these constructions often achieve much faster error decay than MC as \(O(n^{-1+\epsilon})\). A key drawback of deterministic QMC, however, is the lack of a practical and reliable error estimate. Randomized QMC (RQMC) resolves this issue by introducing a randomization that preserves the low-discrepancy structure of the points while restoring the benefits of a stochastic method: the estimator becomes unbiased, and its variance can be estimated from independent replications, thus providing practical error bars.

Among randomized QMC schemes, Owen’s scrambling [6], [7] has become a standard method for digital nets. In Owen’s scrambling each digit of an input point’s base-\(b\) expansion is permuted by an independently chosen random permutation (possibly dependent on higher-order digits). For a set of \(n\) points \(P_n\), Owen proved a bound of the form \[\label{eq:gain-intro2} \operatorname{Var}\bigl(I(f;\widetilde{P}_n)\bigr)\le \frac{\Gamma(P_n)}{n}\operatorname{Var}(f),\tag{4}\] where \(\Gamma(P_n)\) is the maximal gain coefficient associated with the point set. Moreover, for sufficiently smooth integrands the variance decays substantially faster as \(O(n^{-3 + \epsilon})\) (see [8], [9]), explaining the large practical gains of scrambling over plain MC. While recent work has also focused on alternative estimators, such as the median-of-means, to achieve robust and almost optimal order of convergence without a knowledge of integrands [10]–[13], the focus of this paper is on the variance properties of the classical mean-based estimator.

In practice, scrambled Sobol’ sequences [14]–[16] are arguably the most widely used in scientific computing and industry due to their excellent empirical performance. A theoretical drawback, however, is that their maximal gain coefficient \(\Gamma(P_n)\) can grow exponentially with dimension \(d\) [17], [18]. In very high dimensions, this can lead to a situation where scrambled QMC performs substantially worse than plain Monte Carlo, undermining its primary motivation.

Some alternative RQMC constructions yield smaller maximal gain coefficients. For example, scrambled Faure sequences in prime power base \(b\) (with \(d \leq b\)) satisfy \(\Gamma \leq e \approx 2.718\) [7], but they are generally not considered efficient. One reason is that, numerically, the benefits of scrambling appear only when \(n \gtrsim b^{d}\). Other drawbacks are that large bases are required in high dimensions, and arithmetic in base \(b\) is slower than in base \(2\). Scrambled Halton sequences also have small gain coefficients, growing only as \(O(\log d)\) [19], but they typically do not achieve the higher-order variance decay observed for scrambled digital nets.

This paper introduces coarse scrambling, a randomization technique for digital sequences designed to operate on blocks of digits rather than individual digits. The motivation for this block-wise approach stems directly from our primary goal: to preserve the powerful \((0,\boldsymbol{e},d)\)-sequence property of Sobol’ and Niederreiter sequences. Since this property is defined over mixed bases of the form \((b^{e_1},\dots,b^{e_d})\), scrambling entire blocks of digits is the natural and necessary way to maintain this structure—a feat not achieved by standard digit-wise scrambling. By preserving this crucial property, our construction achieves two critical objectives: it controls the maximal gain coefficient to grow only logarithmically with dimension, \(O(\log d)\), while simultaneously retaining the \(O(n^{-3+\epsilon})\) convergence rate for smooth functions.

Our main contributions are as follows:

1. A Unifying Framework and Structure-Preserving Randomization: We introduce the concept of an equidistributed sequence in a mixed base, a novel theoretical framework that, for the first time, provides a unified treatment of both Halton-type sequences and digital \((0,\boldsymbol{e},d)\)-sequences like Sobol’s. Within this framework, we formally define coarse scrambling and prove it preserves the powerful \((0,\boldsymbol{e},d)\)-structure, a key property for high-performance QMC integration.

2. Theoretical Performance Guarantees: We derive an analytical expression for the gain coefficients of coarsely scrambled sequences by extending recent results for scrambled Halton sequences [19]. Our analysis reveals two key properties: (i) for sufficiently smooth integrands, the variance converges at the canonical \(O(n^{-3+\epsilon})\) rate, matching the performance of Owen’s scrambling, and (ii) the maximal gain coefficient grows only logarithmically with dimension, \(O(\log d)\), ensuring robustness in high-dimensional settings.

3. Empirical Validation of the Theoretical Trade-offs: We conduct numerical experiments to validate our theoretical findings and illustrate the practical trade-off between coarse and usual scrambling. The results highlight that the superior accuracy of usual scrambling is driven by its excellent performance on low-dimensional projections. Conversely, for functions with low effective dimensionality, where the total variance is concentrated in the first few coordinates, the performance of coarse scrambling is comparable.

The remainder of the paper is organized as follows. Section 2 recalls basic definitions and notation for digital nets, scrambling, and gain coefficients. In Section 3 we develop the theoretical analysis of gain coefficients for coarse scrambling. Section 4 compares the usual (digitwise) scrambling and coarse theoretically. Section 5 reports numerical experiments.

2 Preliminaries↩︎

We use the following notation. Let \(\mathbb{N}\) be the set of positive integers, \(\mathbb{N}_0 := \mathbb{N}\cup \{0\}\), and \(\mathbb{Z}_b:= \{0,1,\dots, b-1\}\). Let \(\chi(A)\) denote the characteristic function of the event \(A\). Let \(1{:}d := \{1, \dots, d\}\). For \(u \subseteq 1{:}d\), we denote the cardinality of \(u\) by \(|u|\).

2.1 Digital Constructions and Equidistribution↩︎

Digital constructions in prime power base \(b\) were introduced by Niederreiter. We adopt the definition given in [5]. Let \(\mathbb{F}_b\) be a \(b\)-element field and fix a bijection \(\phi \colon \mathbb{Z}_b\to \mathbb{F}_b\) with \(\phi(0) = 0\). If \(b\) is a prime, we usually identify \(\mathbb{F}_b\) with \(\mathbb{Z}_b\) with addition, subtraction, and multiplication defined modulo \(b\).

To assess the quality of point sets and sequences, the notions of \((t,m,d)\)-nets and \((t,d)\)-sequences are widely used to establish their equidistribution property. To introduce these objects, we first define elementary intervals.

2.2 Generalized Niederreiter Sequences↩︎

A well-known and important class of digital sequences is given by the Sobol’ and (generalized) Niederreiter sequences. The Sobol’ sequence, defined in [14], can be regarded as a generalized Niederreiter sequence in base \(2\) whose base polynomials are chosen to be primitive (see [20] for details). The Niederreiter sequence was introduced in [21] and later fully generalized by Tezuka in [22].

Generalized Niederreiter sequences are important examples of \((t,d)\)-sequences, as stated in the following theorem [5].

As the theorem indicates, the quality of the generalized Niederreiter sequence depends on the degrees of the base polynomials. To achieve a small \(t\)-value, one must choose base polynomials with low degrees. Table 1 summarizes the number of available irreducible and primitive polynomials over \(\mathbb{F}_2\) for a given degree. For example, when constructing a Sobol’ sequence (where primitive polynomials are typically used), if one wishes to keep the degrees of the base polynomials at most \(7\), the maximum possible dimension is \(d = 1 + \sum_{n=1}^{7} P_2(n) = 37\), where the first dimension is specially chosen with \(p_1(x) = x\).

2.3 Stronger Equidistribution and Further Generalization↩︎

A smaller \(t\)-value is preferable, since it indicates that the point set satisfies an equidistribution property with respect to finer partitions of \([0,1)^d\). However, the ideal case \(t=0\) is unattainable in high dimensions. For instance, if \(b\) is prime, a \((0,b+1)\)-sequence in base \(b\) does not exist. To overcome this limitation, generalized notions of \((t,\boldsymbol{e},m,d)\)-nets and \((t,\boldsymbol{e},d)\)-sequences were introduced [23]. In these notions, the class of elementary intervals under consideration is restricted, which in turn allows for smaller \(t\)-values.

Under this generalized framework, generalized Niederreiter sequences exhibit an ideal equidistribution property, as shown by Tezuka [23].

In this paper, we further generalize the concept of \((0,\boldsymbol{e},d)\)-sequences to arbitrary bases.

This framework includes well-known sequences as special cases. For example, the Halton sequence is an equidistributed sequence in base \((p_1, \dots, p_d)\), where \(p_j\) denotes the \(j\)th prime number. Furthermore, the notion of a \((0,\boldsymbol{e},d)\)-sequence in base \(b\) is identical to that of an equidistributed sequence in base \((b^{e_1}, \dots, b^{e_d})\). Therefore, Tezuka’s theorem implies that the generalized Niederreiter sequences are equidistributed sequences in base \((b^{e_1}, \dots, b^{e_d})\), where \(e_j = \deg(p_j)\).

2.4 Scrambling↩︎

Scrambling is a procedure that applies a random permutation to the digits of a point’s coordinates, aiming to improve the uniformity of the point set. We begin with a general definition of scrambling.

For a scramble to be effective in randomized quasi-Monte Carlo, it should satisfy certain uniformity properties. These uniformity properties are crucial for variance analysis.

Note that the two-fold uniformity implies the one-fold uniformity. Several methods exist to implement a scramble with these properties. Owen [6] introduced fully nested scrambling, where the permutations are chosen independently and uniformly from the \(b!\) possible permutations on \(\mathbb{Z}_b\). This method is a conceptual ideal but can be costly both in storing the permutations and applying them. As a fast and practical alternative, the affine matrix scramble has been proposed [24]–[26]. As in Section 2.1, we fix a bijection \(\phi \colon \mathbb{Z}_b\to \mathbb{F}_b\) with \(\phi(0)=0\). Further, for a real number \(x \in [0,1)\) with the base-\(b\) expansion \(x = \xi_1 b^{-1} + \xi_2 b^{-2} + \cdots\), where \(\xi_k \in \{0, 1, \dots, b-1\}\), we denote the corresponding vector of its digits as \[\label{eq:real2vec} \vec{x} := (\phi(\xi_{1}),\phi(\xi_{2}), \dots)^\top.\tag{6}\]

Both fully nested scrambling and affine matrix scrambling possess two-fold uniformity. A key consequence of such scrambles is that a single scrambled point \(\sigma(\boldsymbol{x})\) is uniformly distributed on \([0,1)^d\). Furthermore, scrambling preserves the equidistribution property of the entire sequence. The proof proceeds similarly to the case of \((t,d)\)-sequences [6] and is therefore omitted.

For sequences like the Sobol’ and Niederreiter sequences, which are \((0,\boldsymbol{e},d)\)-sequences in base \(b\), their strong equidistribution property is defined with respect to the mixed base \((b^{e_1},\dots,b^{e_d})\). To preserve this specific structure, it is natural to apply scrambling in this mixed base. We refer to this as a coarse scramble, a central concept of this paper, in contrast to the usual fine or digitwise scramble in base \((b,b,\dots,b)\). This can be achieved by full nested scramble or affine matrix scrambling over a field \(\mathbb{F}_{b^{e_j}}\) with mixed bases. Here we suggest the block affine matrix scramble, a method that generalizes the affine matrix scramble to operate on blocks of digits and directly corresponds to performing operations in the mixed base.

This generalized scramble retains the essential properties of the original. The following result is an analogue of the standard case [24], [26], and its proof is omitted.

2.5 Nested ANOVA decomposition and gain coefficients↩︎

To analyze the variance of an RQMC estimator, we decompose it into contributions from the integrand function and the point set. The function’s structure is captured by the nested ANOVA (Analysis of Variance) decomposition, which we introduce first.

The nested ANOVA decomposition breaks down a function \(f: [0,1)^d \to \mathbb{R}\) into a sum of orthogonal components. This decomposition is typically defined using an orthonormal basis, such as the Haar basis [7] or Walsh basis [5]; the equivalence of our formulation to such definitions is established in Appendix 6. The decomposition is constructed with respect to the mixed base \(\boldsymbol{b}\); while this dependency is often omitted from the following notation for brevity, all resulting ANOVA components implicitly depend on \(\boldsymbol{b}\). The decomposition is given by \[\label{eq:ANOVA} f(\boldsymbol{x}) = f_{\emptyset} + \sum_{\emptyset \neq u \subseteq 1{:}d} \sum_{\boldsymbol{k}\in \mathbb{N}_0^{|u|}} \beta_{u,\boldsymbol{k}}(\boldsymbol{x}_u),\tag{7}\] where \(f_{\emptyset} = \int_{[0,1)^d} f(\boldsymbol{x}) \,d\boldsymbol{x}\), and the component \(\beta_{u,\boldsymbol{k}}(\boldsymbol{x}_u)\) is defined via auxiliary functions \(\tilde{\beta}_{u,v,\boldsymbol{k}}(\boldsymbol{x}_u)\) as \[\begin{align} \beta_{u,\boldsymbol{k}}(\boldsymbol{x}_u) &= \sum_{v \subseteq u} (-1)^{|u|-|v|} \tilde{\beta}_{u,v,\boldsymbol{k}}(\boldsymbol{x}_u), \tag{8} \\ \tilde{\beta}_{u,v,\boldsymbol{k}}(\boldsymbol{x}_u) &= m_{u,v,\boldsymbol{k}} \int_{\operatorname{BOX}_{u,v,\boldsymbol{k}}(\boldsymbol{x}_u)} d\boldsymbol{y}_u \int_{[0,1)^{d-|u|}} f(\boldsymbol{y}_u,\boldsymbol{z}_{-u}) d\boldsymbol{z}_{-u},\tag{9} \end{align}\] where \(\operatorname{BOX}_{u,v,\boldsymbol{k}}(\boldsymbol{x}_u)\) is the elementary interval containing \(\boldsymbol{x}_u\) defined by \[\label{eq:BOX} \operatorname{BOX}_{u,v,\boldsymbol{k}}(\boldsymbol{x}_u) = \prod_{j \in v} \left[\frac{\lfloor x_j b_j^{k_j+1} \rfloor}{b_j^{k_j+1}}, \frac{\lfloor x_j b_j^{k_j+1} \rfloor+1}{b_j^{k_j+1}} \right) \times \prod_{j \in u\setminus v} \left[\frac{\lfloor x_j b_j^{k_j} \rfloor}{b_j^{k_j}}, \frac{\lfloor x_j b_j^{k_j} \rfloor+1}{b_j^{k_j}} \right),\tag{10}\] and \((m_{u,v,\boldsymbol{k}})^{-1}\) is the volume of \(\operatorname{BOX}_{u,v,\boldsymbol{k}}(\boldsymbol{x}_u)\), namely \[m_{u,v,\boldsymbol{k}} := \prod_{j \in v} b_j^{k_j+1} \prod_{j \in u \setminus v} b_j^{k_j}.\] By definition, \(\tilde{\beta}_{u,v,\boldsymbol{k}}\) is constant on each elementary interval \(\operatorname{BOX}_{u,v,\boldsymbol{k}}(\boldsymbol{x}_u)\), taking the value of the average of \(f\) over it. These components are related by the following identity: \[\label{eq:beta2} \sum_{\substack{\boldsymbol{k}' \in \mathbb{N}_0^{|u|} \\ k'_j \le k_j \, (\forall j \in u)}}\beta_{u,\boldsymbol{k}'}(\boldsymbol{x}_u) = \tilde{\beta}_{u,u,\boldsymbol{k}}(\boldsymbol{x}_u).\tag{11}\] This is also shown in Appendix 6. The components \(\beta_{u,\boldsymbol{k}}\) are mutually orthogonal. By setting \[\sigma_{u,\boldsymbol{k}}^2 := \operatorname{Var}\beta_{u,\boldsymbol{k}},\] it follows from orthogonality (Parseval’s identity) that the total variance of the function is \[\operatorname{Var}f = \sum_{\emptyset \neq u \subseteq 1{:}d} \sum_{\boldsymbol{k}\in \mathbb{N}_0^{|u|}} \sigma_{u,\boldsymbol{k}}^2.\] From this, the variance of the standard Monte-Carlo (MC) estimator with \(n\) independent random points is given by \[\label{eq:MC-variance} \frac{\operatorname{Var}f}{n} = \dfrac{1}{n}\sum_{\emptyset \neq u \subseteq 1{:}d} \sum_{\boldsymbol{k}\in \mathbb{N}_0^{|u|}} \sigma_{u,\boldsymbol{k}}^2.\tag{12}\]

Next, we introduce a tool to characterize the structure of the point set. The gain coefficients quantify how a specific \(n\)-point set \(P_n\) is expected to perform for each ANOVA component when used in an RQMC estimate. The definition for base \(b\) was given in [7], and its extension to mixed bases was introduced in [19].

An explicit formula for the gain coefficients is known [19], \[\label{eq:Gtilda95formula} \widetilde{G}_{u,\boldsymbol{k}}(P_n) = \sum_{v \subseteq u} H_{u,v} C_{u,v,\boldsymbol{k}}(P_n),\tag{13}\] where \[\begin{align} H_{u,v} &= (-1)^{|u|-|v|} m_{u,v,\boldsymbol{0}} = \prod_{j \in v} b_j \prod_{j \in u \setminus v} (-1), \tag{14}\\ C_{u,v,\boldsymbol{k}}(P_n) &= \frac{1}{n^2} \sum_{i=0}^{n-1}\sum_{i'=0}^{n-1} \chi(\text{\boldsymbol{x}_i and \boldsymbol{x}_{i'} are in the same elementary (\boldsymbol{k}+ \boldsymbol{1}_v)-interval}). \tag{15} \end{align}\]

With these two tools—the ANOVA decomposition for the function and the gain coefficients for the point set—we can now state the variance of the RQMC estimator [7], [19].

By comparing the RQMC variance ?? with the MC variance 12 , it becomes clear that the gain coefficients \(G_{u,\boldsymbol{k}}(P_n)\) measure how much the underlying point set \(P_n\) reduces the variance compared to ordinary MC sampling for each frequency component \((u, \boldsymbol{k})\) of the integrand. A smaller gain coefficient indicates a greater variance reduction for that component.

3 Gain coefficients of scrambled points↩︎

3.1 General results↩︎

The theoretical properties of gain coefficients for scrambled equidistributed sequences can be established by generalizing the recent analysis of Owen and Pan for Halton sequences [19]. The key insight enabling this generalization is that their foundational formula for the counting term \(C_{u,v,\boldsymbol{k}}(P_n)\) (re-derived here as Lemma 3.1) does not depend on the specific structure of the Halton sequence, but only on the defining property of an equidistributed sequence. This allows us to extend their main results to our broader framework.

This observation enables us to generalize results in [19], as follows.

As a corollary, the variance of the RQMC error using a scrambled equidistributed sequence decays as \(o(1/n)\), which is strictly better than plain Monte Carlo. The proof is identical to that of [19], with \(\boldsymbol{x}_i\) replaced by \(\tilde{\boldsymbol{x}}_i\).

3.2 Gain coefficients of scrambled \((0,\boldsymbol{e},d)\)-sequences↩︎

We now specialize these general results to scrambled \((0,\boldsymbol{e},d)\)-sequences, in particular Sobol’ and Niederreiter sequences. Since \((0,\boldsymbol{e},d)\)-sequences in base \(b\) are equidistributed in base \((b^{e_1}, \dots, b^{e_d})\), we will state the results for the latter.

First we consider non-asymptotic results for \(n=b^m\).

Indeed, the exact maximum gain coincides with the upper bound given in Corollary 2.

From Corollary 2 and Eq. ?? , we have the bound of the variance of the scrambled points estimator.

Lastly, we show a slow growth of the maximal gain coefficients of the coarsely scrambled Sobol’ or Niederreiter sequences.

This logarithmic growth of the maximal gain coefficient stands in stark contrast to the exponential growth observed for standard Owen scrambling of Sobol’ sequences, highlighting the primary theoretical advantage of coarse scrambling in high dimensions for the worst case.

4 Comparison of coarse and usual scrambling↩︎

In this section, we compare the coarse scrambling and the usual scrambling in base \(b\). We denote the two bases of interest as the coarse base, \(\boldsymbol{b}^{({\mathrm{C}})} := (b^{e_1}, \dots, b^{e_d}),\) and the usual base, \(\boldsymbol{b}^{({\mathrm{U}})} := (b, \dots, b).\) Hereafter, superscripts or subscripts C and U will refer to quantities related to the coarse and usual bases, respectively. We first give a relation between the nested ANOVA decomposition in coarse and usual bases.

Hereafter, let \(\mathcal{S}\) be a \((0,\boldsymbol{e},d)\)-sequence in base \(b\). We denote the variance of the estimator using points scrambled in the coarse and usual bases by \(\sigma^2_C(n)\) and \(\sigma^2_U(n)\), respectively: \[\begin{align} \sigma^2_{\mathrm{C}}(n) &:= \operatorname{Var}\left(\frac{1}{n}\sum_{i=0}^{n-1}f(\tilde{\boldsymbol{x}}_i^{(C)})\right), \tag{17} \\ \sigma^2_{\mathrm{U}}(n) &:= \operatorname{Var}\left(\frac{1}{n}\sum_{i=0}^{n-1}f(\tilde{\boldsymbol{x}}_i^{(U)})\right). \tag{18} \end{align}\] Lemma 3 implies that \[\label{eq:sigma-seq-basechange} \sigma^2_{\mathrm{C}}(n) = \sum_{u \subseteq 1{:}d} \sum_{\boldsymbol{k}\in \mathbb{N}_0^{|u|}} G_{u,\boldsymbol{k}}^{({\mathrm{C}})}(n) (\sigma_{u,\boldsymbol{k}}^{({\mathrm{C}})})^2 = \sum_{u \subseteq 1{:}d} \sum_{\boldsymbol{l}\in \mathbb{N}_0^{|u|}} G_{u,(\lfloor l_j/e_j \rfloor)_{j \in u}}^{({\mathrm{C}})}(n) (\sigma_{u,\boldsymbol{l}}^{({\mathrm{U}})})^2.\tag{19}\] Furthermore, for \(n=b^m\), Lemma 3 combining with Corollary 2 and the inequality \[e_j(\lfloor l_j/e_j \rfloor + 1) \le e_j(l_j/e_j + 1) = l_j + e_j\] implies that \[\begin{align} \sigma^2_{\mathrm{C}}(b^m) &\le \sum_{\emptyset \neq u \subseteq 1{:}d} \sum_{\substack{\boldsymbol{k}\in \mathbb{N}_0^{|u|} \\ \sum_{j \in u} e_j(k_j+1) > m}} G_{u,\boldsymbol{k}}^{({\mathrm{C}})}(b^m) (\sigma_{u,\boldsymbol{k}}^{({\mathrm{C}})})^2 \notag\\ &= \sum_{\emptyset \neq u \subseteq 1{:}d} \sum_{\substack{\boldsymbol{l}\in \mathbb{N}_0^{|u|} \\ \sum_{j \in u} e_j(\lfloor l_j/e_j \rfloor + 1) > m}} G_{u,(\lfloor l_j/e_j \rfloor)_{j \in u}}^{({\mathrm{C}})}(b^m) (\sigma_{u,\boldsymbol{l}}^{({\mathrm{U}})})^2 \tag{20}\\ &\le \sum_{\emptyset \neq u \subseteq 1{:}d} \Gamma_u \sum_{\substack{\boldsymbol{k}\in \mathbb{N}_0^{|u|} \\ \sum_{j \in u} l_j > m - \sum_{j \in u} e_j}} (\sigma_{u,\boldsymbol{l}}^{({\mathrm{U}})})^2, \tag{21} \end{align}\] where \(\Gamma_u\) is defined as in ?? . Let us compare 21 with the general bound for the usual scramble of the (digital) \((t,m,d)\)-net (see [17], [18]) as \[\label{eq:general95bound95dense} \sigma^2_{\mathrm{U}}(b^m) \le \sum_{\emptyset \neq u \subseteq 1{:}d} b^{t_u}\left(\dfrac{b}{b-1}\right)^{|u|-1} \sum_{\substack{\boldsymbol{l}\in \mathbb{N}_0^{|u|} \\ \sum_{j \in u} l_j > m - t_u}} (\sigma_{u,\boldsymbol{l}}^{({\mathrm{U}})})^2,\tag{22}\] where \(t_u\) denotes the \(t\)-value of the projection of the digital sequence \(\mathcal{S}\) onto the coordinates in \(u\). In our case, by Theorem 1 we can take \(t_u = \sum_{j \in u}(e_j-1)\). It follows that the bound 21 has more terms, but each with a smaller constant factor. We can formalize this observation in the following “meta” theorem. This can be applied to [9] and we obtain a corresponding corollary.

Let us proceed more detailed comparison. Comparing 21 and 22 , The variance \(\sigma^2_{\mathrm{C}}(b^m)\) of the estimator based on coarse scrambling may be larger for subsets \(u\) with small \(t_u\) and small cardinality \(|u|\). In some cases, it can be shown that the precise \(t\)-value satisfies \(t_u = \sum_{j \in u} (e_j-1)\) [29], although this does not hold for all subsets.

The most important case is the one-dimensional projection, i.e., \(u = \{j\}\). For such one-dimensional projections, the Sobol’ sequence (or, more generally, Niederreiter sequences with non-singular upper-triangular generating matrices) forms a \((0,1)\)-sequence in base \(b\), i.e., \(t_{\{j\}} = 0\). Thus the variance for the \(\{j\}\)-component for the usual scramble is given by \[\label{eq:1dim-dense} \sum_{l \ge m} (\sigma_{\{j\},l}^{({\mathrm{U}})})^2.\tag{23}\] On the other hand, for the coarse scramble, Corollary 2 and 20 imply that the variance for the \(\{j\}\)-component is bounded by \[\label{eq:1dim-coarse} \sum_{(\lfloor l/e_j \rfloor + 1)e_j > m} G_{\{j\},(\lfloor l_j/e_j \rfloor)}^{({\mathrm{C}})}(b^m) (\sigma_{\{j\},l}^{({\mathrm{U}})})^2 \le \sum_{l \ge \lfloor m/e_j \rfloor e_j} (\sigma_{\{j\},l}^{({\mathrm{U}})})^2\tag{24}\] Thus, since \(l \ge m\) implies \(l \ge \lfloor m/e_j \rfloor e_j\), the exact variance for the usual scramble in 23 is no greater than the bound for the coarse scramble in 24 . This means that, for one-dimensional projections, the coarse scramble can be significantly worse than the usual one. To avoid this drawback, we may choose \(m\) to be a multiple of \(e_j\), in which case the two bounds coincide. This avoidance is observed in our numerical experiment.

In practical applications, such as those in finance, integrands often have what is broadly termed a low effective dimension [30]. This notion is typically specified in two main ways. The first is having a low superposition dimension, where the variance is dominated by interactions among a small number of variables (small \(|u|\)). The second is having a low truncation dimension, where the variance is concentrated in the first few coordinates. As our analysis shows, the trade-off between coarse and usual scrambling depends on the integrand’s effective dimension. For functions with a low superposition dimension, where low-order interactions (especially 1D projections) dominate, usual scrambling is superior due to its perfect \(t=0\) property on 1D components. For functions with a low truncation dimension in high-dimensional settings, where the influence of higher-order interactions is non-negligible, the superior \(O(\log d)\) scaling of the gain coefficient \(\Gamma_u\) for coarse scrambling can become the decisive advantage, making it a compelling alternative.

5 Numerical Experiments↩︎

To empirically validate our theoretical findings, we conduct a series of numerical integration experiments comparing the performance of coarse scrambling with usual scrambling on Sobol’ sequences.

5.1 Experimental Setup↩︎

For both scrambling strategies, we estimate the integral of test functions using a RQMC estimator. The number of sample points is set to \(n=2^m\), where \(m\) ranges from 1 to 16. To assess the estimator’s stability, we compute the root mean square error (RMSE) over \(R=200\) independent replications for each value of \(n\). For sufficiently smooth integrands, theory predicts a convergence rate of approximately \(O(n^{-1.5})\) for the RMSE. The usual and coarse scrambling methods are implemented using affine matrix scramble and block affine matrix scramble, respectively.

We use two distinct test functions designed to probe different aspects of the scrambling methods:

1. A 37-dimensional linear function. As its variance is entirely contained in one-dimensional projections (\(|u|=1\)), which has low superposition dimension.

2. A 100-dimensional weighted smooth function. This function is designed to have a low truncation dimension.

The results, presented as log-log plots of the RMSE against the number of points \(n\), are shown in Figure 1.

5.2 Results and Discussion↩︎

5.2.1 Linear Function↩︎

Our first test integrand is the 37-dimensional linear function: \[f(\boldsymbol{x}) = \sum_{j=1}^{37} x_j.\] The dimension \(d=37\) was chosen because for the Sobol’ sequence, the corresponding generator polynomial degrees result in block sizes \(e_j\) up to 7 (see Table 1).

As shown in Figure 1 (a), the results align perfectly with the theoretical analysis in Section 4. The usual scrambling method consistently outperforms coarse scrambling, exhibiting a smooth convergence rate close to \(O(n^{-1.5})\). In contrast, the coarse scrambling method shows a stepped convergence. The error convergence follows a baseline rate of approximately \(O(n^{-1})\), superimposed with significant drops when \(m\) is a multiple of 7 (i.e., at \(n=2^7\) and \(n=2^{14}\)). This confirms our theoretical prediction: for one-dimensional projections, the error bound for coarse scrambling only matches that of usual scrambling when the number of points is chosen appropriately relative to the block sizes. Since many dimensions have a block size of \(e_j=7\), choosing \(m\) as a multiple of 7 synchronizes the error reduction across these dimensions, leading to the observed drops.

5.2.2 Weighted Smooth Function↩︎

Next, we consider the 100-dimensional weighted function: \[g(\boldsymbol{x}) = \prod_{j=1}^{100} \left(1+\frac{1}{j^2}(x_j e^{x_j}-1)\right).\] This integrand is dominated by its first few dimensions. For \(d=100\), roughly half of the dimensions in the Sobol’ sequence correspond to a block size of \(e_j=9\).

The results in Figure 1 (b) show that both scrambling methods perform comparably, with convergence rates faster than \(O(n^{-1})\) but slightly slower than the ideal \(O(n^{-1.5})\). In this scenario, the potential drawbacks of coarse scrambling are negligible. The function’s variance is concentrated in the low-index dimensions, which have small block sizes (\(e_j=1, 2, \dots\)). The performance degradation from the larger block sizes associated with high-index dimensions is insignificant because these dimensions contribute very little to the total variance. This suggests that for functions with low truncation dimension, coarse scrambling provides a comparable level of accuracy to usual scrambling.

6 Equivalence of ANOVA Formulations↩︎

This appendix serves two purposes. First, we demonstrate that our formulation of the ANOVA decomposition components is equivalent to the standard formulation via the Walsh basis. Second, we prove the identity in 11 . For simplicity, we consider the case of a common base, i.e., \(\boldsymbol{b}= (b,b,\dots,b)\), as the generalization to the mixed base case is straightforward.

Let \(\mathrm{wal}_{\boldsymbol{h}}(\boldsymbol{x})\) be the \(\boldsymbol{h}\)-th Walsh function and \(\widehat{f}(\boldsymbol{h})\) be the corresponding Walsh coefficient of \(f\), as defined in [5]. The ANOVA components in [5], which we denote by \(\beta'_{\boldsymbol{l}}(\boldsymbol{x})\), are defined for \(\boldsymbol{l}\in \mathbb{N}_0^d\) as \[\label{eq:beta-dash} \beta'_{\boldsymbol{l}}(\boldsymbol{x}) = \sum_{\boldsymbol{h}\in L_{\boldsymbol{l}}} \widehat{f}(\boldsymbol{h}) \mathrm{wal}_{\boldsymbol{h}}(\boldsymbol{x}),\tag{25}\] where the index set \(L_{\boldsymbol{l}}\) is given by \[L_{\boldsymbol{l}} = \{ \boldsymbol{h}= (h_1,\dots,h_d) \in \mathbb{N}_0^d \mid b^{l_j-1} \le h_j < b^{l_j} \text{ if } l_j>0, \text{ and } h_j=0 \text{ if } l_j=0 \}.\] Our goal is to show that this definition coincides with ours. Specifically, for an index \(\boldsymbol{l}\in \mathbb{N}_0^d\), let \(u := \{j \in 1{:}d \mid l_j > 0 \}\) and let \(\boldsymbol{k}= (l_j-1)_{j \in u}\). We will show that \(\beta'_{\boldsymbol{l}}(\boldsymbol{x}) = \beta_{u,\boldsymbol{k}}(\boldsymbol{x}_u)\), where \(\beta_{u,\boldsymbol{k}}\) is defined as in 11 .

First, note that for any \(j \notin u\), we have \(l_j=0\) and thus \(h_j=0\) for all \(\boldsymbol{h}\in L_{\boldsymbol{l}}\). Since \(\mathrm{wal}_0(x_j) = 1\) and \(\mathrm{wal}_{\boldsymbol{h}}(\boldsymbol{x}) = \prod_{j=1}^d \mathrm{wal}_{h_j}(x_j)\), the term \(\mathrm{wal}_{\boldsymbol{h}}(\boldsymbol{x})\) does not depend on \(x_j\) for \(j \notin u\). Consequently, \(\beta'_{\boldsymbol{l}}(\boldsymbol{x})\) is a function of only \(\boldsymbol{x}_u\).

The set \(L_{\boldsymbol{l}}\) can be expressed using the inclusion-exclusion principle. Let \[\tilde{L}_{\boldsymbol{l},v} := \{ \boldsymbol{h}\in \mathbb{N}_0^d \mid h_j < b^{l_j-\chi(j \in u \setminus v)} \text{ for } 1 \le j \le d\}\] for any \(v \subseteq u\). The indicator function for \(\boldsymbol{h}\in L_{\boldsymbol{l}}\) can then be written as \[\chi(\boldsymbol{h}\in L_{\boldsymbol{l}}) = \sum_{v \subseteq u} (-1)^{|u|-|v|} \chi(\boldsymbol{h}\in \tilde{L}_{\boldsymbol{l},v}).\] Substituting this into 25 yields \[\label{eq:betadash-IE} \beta'_{\boldsymbol{l}}(\boldsymbol{x}) = \sum_{v \subseteq u}(-1)^{|u|-|v|} \left( \sum_{\boldsymbol{h}\in \tilde{L}_{\boldsymbol{l},v}} \widehat{f}(\boldsymbol{h}) \mathrm{wal}_{\boldsymbol{h}}(\boldsymbol{x}) \right).\tag{26}\] The inner sum over \(\boldsymbol{h}\in \tilde{L}_{v}\) is a partial Walsh series of \(f\). As shown in the proof of [5], this sum has an integral representation: \[\begin{align} \sum_{\boldsymbol{h}\in \tilde{L}_{\boldsymbol{l},v}} \widehat{f}(\boldsymbol{h}) \mathrm{wal}_{\boldsymbol{h}}(\boldsymbol{x}) &= \int_{[0,1)^d} f(\boldsymbol{y}) \prod_{j=1}^d b^{l_j-\chi(j \in u \setminus v)}\chi\big(x_j \ominus y_j \in [0,b^{-(l_j-\chi(j \in u \setminus v))})\big) \,d\boldsymbol{y}\notag \\ &= \tilde{\beta}_{u,v,\boldsymbol{k}}(\boldsymbol{x}_u), \label{eq:beta-betadash} \end{align}\tag{27}\] where \(\ominus\) is the \(b\)-adic bitwise subtraction. The integral expression in the first line corresponds precisely to the definition of \(\tilde{\beta}_{u,v,\boldsymbol{k}}(\boldsymbol{x}_u)\) in our framework (see 9 ).

By substituting 27 into 26 , we obtain \[\beta'_{\boldsymbol{l}}(\boldsymbol{x}) = \sum_{v \subseteq u}(-1)^{|u|-|v|} \tilde{\beta}_{u,v,\boldsymbol{k}}(\boldsymbol{x}_u).\] The right-hand side is, by definition in 8 , equal to \(\beta_{u,\boldsymbol{k}}(\boldsymbol{x}_u)\). This establishes the equivalence \(\beta'_{\boldsymbol{l}}(\boldsymbol{x}) = \beta_{u,\boldsymbol{k}}(\boldsymbol{x}_u)\) and completes the proof.

Finally, we prove the identity in 11 . For any \(u \subseteq 1{:}d\) and \(\boldsymbol{k}\in \mathbb{N}_0^{|u|}\), we consider the sum of our ANOVA components \(\beta_{u,\boldsymbol{k}'}\) over all \(\boldsymbol{k}' \le \boldsymbol{k}\) (component-wise). Using the established equivalence \(\beta_{u,\boldsymbol{k}'}(\boldsymbol{x}_u) = \beta'_{\boldsymbol{l}(u,\boldsymbol{k}')}(\boldsymbol{x})\), where \(\boldsymbol{l}(u,\boldsymbol{k}')\) is defined such that \(l_j = k'_j+1\) for \(j \in u\) and \(l_j=0\) otherwise, we have \[\begin{align} \sum_{\substack{\boldsymbol{k}' \in \mathbb{N}_0^{|u|} \\ \boldsymbol{k}' \le \boldsymbol{k}}} \beta_{u,\boldsymbol{k}'}(\boldsymbol{x}_u) &= \sum_{\substack{\boldsymbol{k}' \in \mathbb{N}_0^{|u|} \\ \boldsymbol{k}' \le \boldsymbol{k}}} \beta'_{\boldsymbol{l}(u,\boldsymbol{k}')}(\boldsymbol{x}) \\ &= \sum_{\substack{\boldsymbol{k}' \in \mathbb{N}_0^{|u|} \\ \boldsymbol{k}' \le \boldsymbol{k}}} \sum_{\boldsymbol{h}\in L_{\boldsymbol{l}(u,\boldsymbol{k}')}} \widehat{f}(\boldsymbol{h}) \mathrm{wal}_{\boldsymbol{h}}(\boldsymbol{x}). \end{align}\] The union of the disjoint sets \(L_{\boldsymbol{l}(u,\boldsymbol{k}')}\) over all \(\boldsymbol{k}' \le \boldsymbol{k}\) covers exactly the set of all Walsh indices \(\boldsymbol{h}\) where \(h_j < b^{k_j+1}\) for all \(j \in u\) and \(h_j=0\) for \(j \notin u\). This corresponds precisely to the index set we previously denoted as \(\tilde{L}_{\boldsymbol{l}(u,\boldsymbol{k}),u}\). Therefore, we can combine the summations: \[\begin{align} \sum_{\substack{\boldsymbol{k}' \in \mathbb{N}_0^{|u|} \\ \boldsymbol{k}' \le \boldsymbol{k}}} \sum_{\boldsymbol{h}\in L_{\boldsymbol{l}(u,\boldsymbol{k}')}} \widehat{f}(\boldsymbol{h}) \mathrm{wal}_{\boldsymbol{h}}(\boldsymbol{x}) = \sum_{\boldsymbol{h}\in \tilde{L}_{\boldsymbol{l}(u,\boldsymbol{k}),u}} \widehat{f}(\boldsymbol{h}) \mathrm{wal}_{\boldsymbol{h}}(\boldsymbol{x}) = \tilde{\beta}_{u,u,\boldsymbol{k}}(\boldsymbol{x}_u). \end{align}\] The final equality follows directly from 27 . This confirms 11 .

References↩︎

---

1. Faculty of Science, Yamagata University, 1-4-12 Kojirakawa-machi, Yamagata, 990-8560, Japan (kosuke-suzuki@sci.kj.yamagata-u.ac.jp)↩︎

2. The work was supported by JSPS KAKENHI Grant Number 24K06857.↩︎