On the generation of direction numbers for Sobol sequences and their application to Quasi Monte Carlo methods

Keywords Monte Carlo · Quasi-random number · Sobol sequences

Monte Carlo methods are a class of algorithms that utilize random numbers for solving problems in numerical computing such as simulation of complex systems, optimization or integration. In the context of finance, Monte Carlo simulations are often utilized for pricing financial instruments and scenario analysis, whereby Monte Carlo is used to generate possible realizations of a model by varying the input model parameters drawn from a predefined distributions. In scenario analysis Monte Carlo is often used to generate possible future trajectories for quantities of interest. After generating multiple scenarios, one can combine them to compute aggregated metrics such as (but not limited to) means and variances of the metrics at different time points, from which risk estimates can be calculated.
Both the task of generating scenarios and computing aggregated metrics can be recast as solutions to integration problems, whereby we attempt to evaluate an integrand $f(x)$ over an integration domain $\Omega$, with $x \in \Omega$ representing the model parameters, $f(\cdot)$ representing the model and $I$ representing the result of the integration:
$$I = \int_\Omega f(x) dx$$
While analytic solutions to some integrals exist and can be derived by hand, this endeavour becomes impractical for even moderate increases in integral dimensionality or integrand complexity. In such cases, users may resort to symbolic mathematics software to derive analytic solutions or to numeric approximations. Numeric techniques for evaluating integrals are termed quadrature methods - these attempt to compute a numeric estimate $\hat{I}$ from a set of evaluations of the integrand at different locations in the integration domain $\Omega$. The set of integrand evaluations $f(x_i)$ is then combined in a weighted sum, with weights determined by the quadrature rule, to give the resulting estimate:
$$\hat{I} = \sum_i \omega_i f(x_i)$$
The points $x_i \in \Omega$ at which the integrand $f(x)$ is evaluated are termed abscissa while the weights of each evaluations contribution to the final result are denoted $\omega_i$. The specification of the abscissa $x_i$ and weights $\omega_i$ act to define a given quadrature method. One of the most well known low dimensional example of a quadrature method is the trapezium rule which approximates a definite integral:
$$\begin{align} I &= \int_a^b f(x) dx \\ \hat{I} &= \frac{b-a}{N} \sum_{i=1}^{N}\frac{f(x_{i-1})+f(x_i)}{2}\end{align}$$
For uniformly spaced abscissa $x_i = a + \frac{b-a}{N} i$, the rate of at which the error $|I - \hat{I}|$ shrinks as we increase the number of abscissa $N$ is bounded from the top by:
$$|I - \hat{I}| = -\frac{(b-a)^2}{12N^2}[f`(b) - f`(a)] + O(N^{-3})$$
The estimate error in this case thus scales as $O(N^{-2})$ [12] with the number of points and depends also on the gradient of the function $f`(x)$. Generalizing to $d$ dimensions the above rate scales as $O(N^{-2/d})$ - slowing down and thus requiring more integrand evaluations to achieve the same level of accuracy. Combined with the exponentially increasing number of abscissa in higher dimensions, this limits applicability this technique to high dimensional integrals. To understand why - consider using only 10 abscissa per dimension. In five dimensions ($d=5$), this would required $10^5$ integrand evaluations, while in $d=20$ dimensions this becomes $10^{20}$ - an impractical number.
This motivates the development of stochastic quadrature techniques where, instead of using deterministic abscissa, the points at which the integrand is evaluated are selected at random. Counter-intuitively, such a random selection of points results in a generic procedure allowing a much broader class of integrals to be evaluated, with the precision of the estimate determined only by the number of points and independent of the dimensionality $O(N^{-1/2})$. Furthermore, only some minor smoothness assumptions are placed on the integrand.

Monte Carlo integration are stochastic quadrature techniques utilizing random numbers to perform the integration. Similarly to the deterministic quadratures above, the Monte Carlo integral estimates a $d$-dimensional integral. A key feature of MC is that the $N$ abscissa are selected at random throughout the domain and each evaluation of the integrand is given an equal weighting $\forall i \; \omega_i = 1/N$
$$\hat{I} = V \left( \frac{1}{N} \sum_{i=1}^{N} f(x_i) \right) = V \langle f \rangle$$
Here the `volume` of the integration domain is represented by $V = \int_\Omega dx^n$ while $\langle f \rangle$ represents the average of all integrand evaluations. The choice of abscissa to be independent and identically distributed (IID) across the domain allows us reframe the problem as an evaluation of the expectation:
$$I = \mathbb{E}[f(x)] = \int_\Omega f(x) u(x) dx$$
Where $u(x)$ denotes the uniform density on the domain. Imposing an additional assumption of a bounded integrand satisfying $Var[f(x)] = \sigma^2 \le \infty$ enables the Central Limit Theorem (CLT) that ensures MC estimator $\hat{I}$ convergence to the true value $I$ as the number of abscissa tends to infinity ($N\rightarrow \infty$), a result known as the Law of Large Numbers (LLN):
$$\mathbb{P}\left(\left| \int_\Omega f(x) dx - \frac{1}{N}\sum_i f(x_i) \right| \geq \epsilon \right) \leq \frac{\sigma^2}{N \epsilon^2}$$
This is the key result that enables one to be confident in the results produced by MC integration, providing both an estimate of the value as well as a bound on the error associated with the estimate - obtained by computing the variance $Var[f(x)]$. Finally, the rate of convergence of the estimate increase can be seen to scale as $O(N^{-1/2})$, the standard error of the MC estimate (applying the CLT). This $O(N^{-1/2})$ scaling of estimator error, independently of dimension, is key benefit of Monte Carlo. Additionally, unlike in a deterministic quadrature rule, the number of points $N$ need not be known before-hand - we can continually improve our estimate by simplify evaluating more points without re-discretizing the grid (of abscissa) from scratch - giving the MC algorithm the any-time property. Through using this technique we rely not on the unpredictability of the random numbers but on the uniformity of their distribution within the integration volume and the rate at which we can generate the random numbers. For purposes of performing this calculation, we therefore do not need true random number generators - instead we can use pseudo-random number generators (PRNG). Such generators mimic the properties of true RNGs but are fully deterministic algorithms, where an initial seed (number) controls the resulting sequence, indistinguishable from truly random sequences for our application. A widely used instance of a PRNG is the Mersenne-Twister algorithm [13] that is implemented in many software libraries used for numeric computing.
While such generators are a workhorse for modern numeric computing one may ask if it is possible to attain faster convergence using sequences with improved uniformity properties, without attempting to emulate random numbers at all? The answer to this is yes! From the early 1960s such methods termed Quasi-Monte Carlo have been developed and applied for exactly this purpose - accelerating convergence of numeric integration. In this setting, the PRNG is replaced with so-called low discrepancy sequence generators, designed to produce sequences equi-distributed number sequences. This allows for both deterministic error convergence analysis and improved rates, of order $O(log(N)^d/N)$ [11]. We discuss these sequences next.

An intuitive explanation for why Monte-Carlo works for integration is that in the limit as the number of points, $N\rightarrow \infty$, the point set fills the entire domain $\Omega$ uniformly and averages the integrand together to get an estimate of the integral. A key requirement for this work is that the points fill the space uniformly. Formally, this criterion is defined using a measure called discrepancy (denoted $D^*_N)$ which in turn can be used to define the uniformity.

Definition 1 Uniformly distributed sequence
A sequence $x_1,x_2,...,x_N$ is uniformly distributed if $D^*_N(x_1,...,x_N) \rightarrow 0$ as $N \rightarrow \infty$

To simplify discussion, we will consider the integration domain to by a $d$-dimensional unit hypercube $\Omega = [0,1)^d$, with a lower bound including the origin and the upper bound of unity. Before defining the star-discrepancy $D^*_N$ used above, we, following Owen [23], consider local discrepancy:
$$\delta ( a; x_1, ..., x_N ) = \frac{1}{N} \sum_{i=1}^{N} 1_{x_i \in [0,a)} - \prod_{j=1}^d a_j$$
The local discrepancy is the difference between the fraction of points from the sequence $x_i$ that fall into the interval $[0,a)$ ($\frac{1}{N} \sum_{i=1}^{N} 1_{x_i \in [0,a)}$) and the volume occupied by this interval ($\prod_{j=1}^d a_j$) relative to the entire domain $\Omega = [0,1)^d$. In this way, the local discrepancy is a measure of deviation of the expected number of points in that interval relative to the volume occupied by the interval. The star-discrepancy $D^*_N$ is then the maximum absolute deviation in local discrepancy over all possible intervals $[0,a)$ in the unit hypercube:

Definition 2 Star-Discrepancy
$$D^*_N(x_1,...,x_N) = \underset{a \in [0,1)^d}{sup} |\delta ( a; x_1, ..., x_N )|$$

While in practice there are a multitude of other discrepancy measures, we will not discuss them here, noting only that for a practical sequences their estimation is often an NP-hard problem [14], [15], something hinted at by the definition of the star-discrepancy, requiring the evaluation of all possible sub-volumes of the unit cube. Practical techniques for estimating thus attempt to provide bounds on the discrepancies instead of estimating them exactly.

The key reason for defining the discrepancy measures is their application to deriving theoretical bounds on the convergence of the error. A central result of this analysis is the Koksma-Hlawka (KH) inequality relating a sequences` discrepancy to the approximation error produced when using the sequence for numeric integration. This acts as the analogue of the central limit theorem (CLT) for Monte-Carlo, but has the additional constraint that the function being integrated is required to be Riemann integrable. Unlike the CLT the bound is deterministic (not probabilistic) - it holds always and applies for finite N (not like the CLT which holds in the limit as $N\rightarrow \infty$). The inequality relates the error of approximating the integral $\int_{[0,1)^d}f(x)dx$ with a quadrature approximation $\frac{1}{N}\sum_{i=1}^{N}f(x_i)$ and gives a bound on the absolute error as the product of $V(f)$, describing the variation of the integrand $f(x)$ and the discrepancy $D^*(x_1,...,x_N)$ of the sequence of numbers $x_1,...,x_N$ at which we evaluate the integrand.
$$\epsilon = \left|\frac{1}{N} \sum_{i=1}^{N} f(x_i) - \int_{[0,1)^d}f(x) dx \right| \leq V(f)D^*(x_1,...,x_N)$$
In practice, it is often overly pessimistic and difficult to apply due to the computational complexity of estimating the Hardy-Krause (HK) variation $V(f)$ of the function. This is often more computationally demanding than the computation of the integral itself, requiring the calculation of partial derivations of $f$. In the one-dimensional case, this reduces to:
$$V(f) = \int_\Omega |f`(x)|dx$$
We rarely compute this bound in practice, but it does provide the key connection between the error convergence and the discrepancy $D^*_N$ of the point set used to evaluate the integral. Thus, one can hope to do better than the MC estimate using low discrepancy sequences instead. A sequence is considered to be low discrepancy if it satisfies the following:

Definition 3 Low discrepancy sequence
A sequence $x_1,x_2,...$ of points in $[0,1)^d$ is a low discrepancy sequence if for any $N>1$
$$D^*(x_1,...,x_N) \leq C_d \frac{log(N)^d}{N}$$
where $C_d$ is a constant depending only on the dimension

Combining this definition with the HK inequality, under the assumption that the functions variation $V(f)$ is finite (and $f(x)$ is Riemann integrable), we obtain the the upper bound of the error convergence for a given number of function evaluations $N$ with the integral of dimension $d$:
$$\epsilon \sim O\left(\frac{log(N)^d}{N} \right) \sim O(N^{-1})$$
This compares favourably with the MC bound of $\epsilon \sim O(N^{-1/2})$ and scales with dimension $d$ better than the quadrature bound of $O(N^{-2/d})$ for high dimensions. However, it relies on the central assumption of finite $V(f)$: an assumption that can break when we consider transforming the unit hypercube $[0,1)^d$. An example of such a transformation is the inverse Gaussian cumulative distribution function that maps the limits of the interval $(0,1) \rightarrow (-\infty, \infty)$. This mapping produces an function with an infinite Hardy-Krause variation $V(f) = \infty$. While of concern in theory, practical applications of low discrepancy sequences still demonstrate good performance. Finally, practical applications should consider that, the asymptotic convergence rates may kick in only for sufficiently large sample counts. An example of this is reported by Owen [23] who derived a bound $n \geq b^{2s}$ an s-dimensional sequence in base-b.

There are many techniques for generating low discrepancy sequences. A common theme among them is partitioning the unit interval and generating suitably stratified points by utilizing an $n$ digit (number of positions after the decimal point) base $b$ representation of the each number $x_i$ in the sequence:
$$x_i = \sum_{j=1}^n v_{ij} b^{-j} = (0.v_{i1}v_{i2}v_{i3}...)_b$$
The process of generating the sequence then reduces to generating a set of coefficients $v_ij$ for each $x_i$ in the sequence the sequence $x_1,x_2,...$ such that the resulting sequence has the desired low discrepancy property. Given a base $b$, we can thus represent each value in the sequence by the set of coefficients $v_{i1}v_{i2}v_{i3}...$ allows us to work using integer arithmetic instead of a floating point format.

Sobol sequences are a base 2 low discrepancy sequence, making them comparatively faster to generate than sequences in a another base, with the generation rate being comparable to the rate of generation by PRNGs. To generate a Sobol sequence, we need a primitive polynomial and a set of direction numbers; these define and initialize a recurrence relation that, given the current point $x_i$, generates the next point $x_{i+1}$ for every $x_i$ in the sequence. The first of these, the primitive polynomial, has the form given by:
$$p(x) = x^s + a_1 x^{s - 1} + a_2 x^{s - 2} + \cdots + a_{s - 1} x + 1$$
Here, $s$ denotes the degree of the polynomial, while the coefficients $a_k$ can take values from the field $GF(2) = \{0,1\}$. By treating the coefficients of the polynomials are binary digits, we view each polynomial as a binary number, defining operations of $(+,\times)$ (and their inverses) on the set of all such polynomials, forming a field. The polynomial is then considered primitive if it is minimal and cannot be factored into a product of two any other polynomials, i.e. it is irreducible. Given a primitive polynomial, we can define a recurrence relation that generates each of the digits of the integer representation of $x_i = (0.v_{i1}v_{i2}v_{i3}...)_b$ [6]:
$$v_{i,j} = \left( \bigoplus_{k=1}^{s} a_k v_{i-k,j} \right) \oplus v_{i-d,j} \oplus [v_{i-d,j}/2^d] \text{ for } i > d$$
This recurrence allows us to generate the next term for the given digit, and thereby construct the sequence $x_i$. However, to initialize this recurrence relation, we need the first $v_{0,j}, ..., v_{d,j}$ values; these are the direction numbers. Owing to the parallels between Sobol sequences and linear congruential random number generators we can view these direction numbers as analogous to random seeds. Nonetheless, unlike seeds, we must be careful with the choice of direction numbers, as these can drastically impact the quality of the sequences, particularly high dimensional sequences where we need one set of direction numbers for every dimension. Thus, two central questions to the practical application of Sobol sequences are:

• How do we find good direction numbers?
  

• Given a set of direction numbers, how do we quantify what we mean by good?

As Sobol sequences consist of $N$ randomly generated points $x_1,\dots x_N$, they can be described as point sets $P_N = \{x_1,\dots,x_N\}$ [20], [19]. When the base of the set is chosen to be $b$ and $N = b^m$, each dimension $j$ can be partitioned into $b^{q_j}$ equal parts for some $q_i>0$. In other words, the unit hypercube $[0,1)^s$ is partitioned into $2^{q_1+\dots+q_s}$ rectangles of equal sizes. If each rectangle contains exactly the same number of points from $P_N$, then $P_N$ is $(q_1, q_2,\dots,q_N)$-equidistributed in base $b$. Ultimately, if the equidistribution property holds for smaller values of $q_j$, then the point set can be seen as more uniformly distributed. In order to explore the uniformity properties of $P_N$ further, a $(t,m,s)$-net is defined as follows:

Definition 4 $(t,m,s)$-net
$P_N$ is a $(t,m,s)$-net if and only if it is $(q_1, q_2,\dots,q_N)$-equidistributed whenever $q_1+\dots+q_s=m-t$.

The $t$-value is thus the smallest value for which the $t$-value of $P_N$ is smallest. Given that a greater $t$-value for a given $m$ minimises $q_1+\dots+q_s$, the $t$-value can be seen as a uniformity measure of a point set, where a higher $t$-value indicates greater uniformity.

In addition to being $(t,m,s)$-nets, Sobol sequences can be referred to as $(t,s)$-sequences of points that possess certain stratification properties, which describe how the points in the sequence fill sub-regions of a unit hypercube. More specifically, a $(t,s)$-sequence is defined as:

Definition 5 $(t,s)$-sequence
A $(t,s)$-sequence constructed in base b can be considered an infinite series of $(t,m,s)$-nets for every positive integer m, each net being a set of $b^m$ points, with each elementary interval of volume $1/b^m$ containing at most $b^t$ points.

While being a key feature of Sobol and other low discrepancy sequences, such stratification properties are not unique to them. To illustrate, consider Figure [1] which divides the unit square into a $8\times 8$ grid. The left-most plot shows that deterministic quadratures like the trapezoid rule discussed above also stratify the unit square. Comparatively, the lack of stratification of (pseudo)random sequences can be seen on the right, where some of the boxes have an excess number of points while others are empty. Another feature of the plot is the projections of the sequences onto 1-dimension at the bottom of each figure, shown as histograms. Here, we observe the benefit of the `LDS` sequence relative to the `Grid` scheme where all points in a given axis project onto the same coordinate in 1-dimension, thus reducing the `effective` number of points. On the other hand, the `LDS` points have distinct values and thus do not suffer from this problem, as can be seen by the uniform histogram.

Before considering numeric properties calculated from Sobol sequences, we can consider properties of the direction numbers used to initialize a given sequence. As the direction numbers used to initialize the recurrence relations determine how the sequence is generated, their properties can influence the quality of the resulting sequence. Initially formulated in [17], properties A and A` properties represent uniformity conditions [16] that determine how the points in the sequence stratify the unit hypercube.

Definition 6 Property A [6]
A d-dimensional Sobol sequence with $2^d$ points satisfies property A if any binary segment of the unit cube, obtained by diving in half the unit hypercube along each dimension, contains exactly one point from the sequence.

Definition 7 Property A` [6]
A d-dimensional Sobol sequence with $4^d$ points satisfies property A` if any binary segment of the unit cube, obtained by diving by four the unit hypercube along each dimension, contains exactly one point from the sequence.

Points from a two-dimensional Sobol sequence satisfying both properties A and A` are shown in Figure [2], where one can observe that there is a single point in each subinterval. Thus, properties A and A` can be related to stratified sampling schemes, used for variance reduction [23].
Property A
A d-dimensional Sobol` sequence posesses Property A if and only if:
$$det(V_d) \equiv 1 \:\: (\text{mod } 2)$$
Where the matrix $V_d \in (\{0,1\})^{d \times d}$ and is given by:
$$V_d = \left[ {\begin{array}{ccccc} 1 & 0 & 0 & \cdots & 0 \\ 1 & v_{2,2,1} & v_{3,2,1} & \cdots & v_{d,2,1} \\ 1 & v_{2,3,1} & v_{3,3,1} & \cdots & v_{d,3,1} \\ \vdots & \vdots & \vdots & \ddots & \vdots\\ 1 & v_{2,d,1} & v_{3,d,1} & \cdots & v_{d,d,1} \\ \end{array} } \right]$$
Here, the binary values $v_{i,j,k}$ denote the $k^{th}$ binary digit of the $i^{th}$ direction number for the dimension $j$, where the entire direction number is represented by $v_{i,j} = (0.v_{i,j,1}v_{i,j,2}...)_2$.
Property A`
A d-dimensional Sobol` sequence posesses Property A` if and only if:
$$det(U_d) \equiv 1 \:\: (\text{mod } 2)$$
Where the matrix $U_d \in (\{0,1\})^{2d \times 2d}$ is given by:
$$U_d = \left[ {\begin{array}{ccccccc} v_{1,1,\textcolor{red}{1}} & v_{1,1,\textcolor{blue}{2}} & v_{2,1,\textcolor{red}{1}} & v_{2,1,\textcolor{blue}{2}} & \cdots & v_{d,1,\textcolor{red}{1}} & v_{d,1,\textcolor{blue}{2}} \\ v_{1,2,\textcolor{red}{1}} & v_{1,2,\textcolor{blue}{2}} & v_{2,2,\textcolor{red}{1}} & v_{2,2,\textcolor{blue}{2}} & \cdots & v_{d,2,\textcolor{red}{1}} & v_{d,2,\textcolor{blue}{2}} \\ v_{1,3,\textcolor{red}{1}} & v_{1,3,\textcolor{blue}{2}} & v_{2,3,\textcolor{red}{1}} & v_{2,3,\textcolor{blue}{2}} & \cdots & v_{d,3,\textcolor{red}{1}} & v_{d,3,\textcolor{blue}{2}} \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & . \\ v_{1,2d,\textcolor{red}{1}} & v_{1,2d,\textcolor{blue}{2}} & v_{3,2d,\textcolor{red}{1}} & v_{3,2d,\textcolor{blue}{2}} & \cdots & v_{d,2d,\textcolor{red}{1}}& v_{d,2d,\textcolor{blue}{2}} \\ \end{array} } \right]$$
Similarly to Property A, the binary values $v_{i,j,k}$ represent the $k^{th}$ binary digit of the $(i,j)^{th}$ direction number $v_{i,j}$ given by: $v_{i,j} = (0.v_{i,j,1}v_{i,j,2}...)_2$.
Properties A and A` in adjacent dimensions
In the case when the dimensionality, $d$, of the Sobol sequence is large, we may choose to impose properties A and A` on a subset of the $\{1,...,d\}$ dimensions. One such choice can be subsets of variables forming a consecutive sequence, such that $\forall i \in \{1,...d-k\}$ the dimensions $\{i,i+1,...,i+k\}$ satisfy $A$ and/or $A`$. These `adjacency` properties generalize $A$ and $A`$ into $A_k$ and $A`_k$, with the former being special cases where $k=d$. Such a generalization is useful when attempting to construct Sobol sequences in practice in order to reduce the search algorithm time complexity, while satisfying $A$ and/or $A$ only in adjacent dimensions. For a more in-depth discussion of $A_k$ and $A`_k$, see [1].
Given that the conditions for properties A and A` can be satisfied independently, a given sequence can satisfy both/one of/neither of these properties. As mentioned above, when searching for practical sequences, we may choose to utilize any combination of $A_k$ and $A`_k$, as well as potentially different values of $k$. This affects the algorithmic complexity through the size of the determinants of direction numbers we need to compute.

Although the properties described above are useful theoretically, we must also consider practical sequence equi-distribution properties for finite numbers of samples. While univariate Sobol sequences are always well distributed, higher dimensional sequences may exhibit significant co-dependence between individual dimensions for a finite number of samples. This can be illustrated by considering two dimensional projections of multidimensional Sobol sequences (Figure [3]), where we can contrast a well distributed two dimensional projection in Figure [3]a with two highly co-dependent projections showing in Figure [3]b and [3]c. Aside from filling the space non-uniformly, the poor equidistribution of the points lead to integration errors as some regions in the integration domain are completely empty. While increasing the number of points will asymptotically fill in such regions, these asymptotes will occur for impractically large samples. Thus, it is desirable to design sequences without such pathologies in the first place. One practical measure of the quality of the two dimensional projections is the Pearson correlation coefficient which we can compute after transforming each one-dimensional sequences through the inverse cumulative Gaussian distribution. Using this, we can quickly estimate the quality of the two-dimensional projections in the multidimensional sequence. We can construct a metric to minimize using this estimate, such as the maximum absolute correlation, that can be used when searching for direction number sets with the aim of generating sequences with improved projections.
While the two-dimensional projection evaluated at a low sample counts may show high correlations, a simple approach to reducing them, without searching for an alternative low discrepancy sequence, is to simply increase the number of points sampled - in the limit as the number of points $n\rightarrow \infty$ the sequence discrepancy tends to zero. This space is filled uniformly and thus the correlations go to zero as well.
An alternative approach to reducing, but not eliminating, correlations is to scramble the individual one dimensional sequences using a method like Owen`s scrambling (see Section [3]). This preserves the discrepancy of the sequence and shifts the distribution of correlations towards lower values by permuting the elements of a given univariate sequence (Figure [7]). Being based on permutations of the unit interval, the scrambling method is thus unlikely completely wipe out all pathological two-dimensional projections, but one would expect it to improve the overall quality of the multidimensional sequence.

In this section, we describe two mathematical operations that need to performed when generating low discrepancy sequences: direction number generation and determinant calculation. Direction number generation involves mapping any integer to a direction number list, while determinant calculation is used to ascertain whether the direction numbers satisfy properties A and A`. We provide the theory and mathematical formulation for these concepts here and refer to them in section [Sequence: tkrg-a-ap5 (this work)] when searching for direction numbers.

To understand the process of direction number generation consider that they form an ordered sequence. A consequence of this is that there exists a mapping that converts the index of the sequence $\in \mathbb{N}$ into the corresponding set of direction numbers. To help explain this we begin by considering integer sequences and their representation in an arbitrary in-place system like binary or decimal. When using such `in-place` systems a number is written as a sequence of digits in some base. The `base` assigns a possible number of values that a given digit can take. For example, every binary digit can take on either 0 or 1, every decimal digit can take on any value in {0,1,2...9}.
For common systems each digit has the same number of values, so the digit denoting the number of 10`s (in base 10) can take on the same number of values as the digit denoting the number of millions ($10^6$). For a given base $a$, we denote $a_i$ as the number of combinations for the $i^{th}$ digit, with $i \geq 1$. A given number with $n$ digits in base $a$ then consists of the following digits combination values:
$$a_1 \;a_2\; a_3\; a_4\; a_5\; a_5\; ... \;a_n$$
The total number of possible combinations this digit string can represent is $\prod_{i=1}^{n} a_i$. In the example of base 2 with $n=10$ bits, we get $\prod_{i=1}^{10} 2 = 2^{10}$. The number of combinations $a_i$ relates to how a number $x$ can be represented in the base. Consider the binary system again:
$$x = x_1 2^0 + x_2 2^1 + x_1 2^2 + ... x_n 2^{n-1} = \sum_{i=1}^{n} x_i 2^{i-1}$$
Denoting the powers of $2^{i-1}$ as $b_i$ allows the product to be represented as:
$$x = \sum_{i=1}^{n} x_i b_i$$
Where $b_i = 2^{i-1}$. For the case of the binary system, this means that the $n^{th}$ term is:
$$b_n = \prod_{i=0}^{n-1} a_i$$
Where we define $a_0 = 1$ to make $b_1 = 1$. In this way, we can represent a given integer $x$ using a non-constant base system, where the base of a given digit $i$ is $a_i$. Using this, we can then represent the coefficients $x_n$ of the number $x = \sum_{i=1}^{n} x_n b_n$ by computing:
$$x_n = (x \;mod\; b_n) \;/\; a_n$$Counting with direction numbers
Consider the recurrence relation 3.6 of [1]. To be able to construct a given Sobol sequence using direction numbers $v_{i,j}$ (where $i$ denotes the position of the direction number in the dimensions sequence, and $j$ denotes the dimension), we need a total of $s_j$ initialization direction numbers: $v_{0,j}, v_{1,j}, ... v_{s_j,j}$. These direction numbers can be obtained from their integer representation $m_{i,j} = v_{i,j} 2^i$. For a given position $i$, the corresponding direction number $m_{i,j}$ must satisfy:

• $m_{i,j} \geq 1$
  

• $m_{i,j} \leq 2^i$
  

• $m_{i,j}$ is odd.

Thus, for the position $i$, there are $2^{i-1}$ possible direction numbers $m_{i,j}$. The number of coefficients $a_i$ is therefore calculated using the following equation:
$$a_i = 2^{i}$$
Hence, we can view the generation of direction numbers as decomposing a given number $x$ into a non-uniformly sized digit representation with the number of combinations per digit given by $a_i$. Using this method, we can construct functions that map a given number $x$ onto a corresponding direction number list by treating each of the $s_j$ direction numbers as a digit.

A Galois field consists of only two elements $\{0, 1\}$, and defines the following addition ($+$) and multiplication ($\times$) operators:

While this appears to be a trivial modification to standard addition and multiplication defined on the real numbers, the consequences of this for the calculation of determinants in high dimensions and on computers are profound. The first of these is the ability to use one bit per matrix element giving a $\times 64$ per-element memory reduction. This allows us to represent a $\times 8$ larger matrix in $GF(2)$ using the same memory, compared to a 64-bit floating point number. The second set of consequences relates to the algorithms for calculating determinants; we must modify these to account for the new arithmetic. In python, this can be done by casting numpy matrices to contain elements in GF(2) using the `galois` python library. There are other aspects of GF(2) matrices that are more mathematical. The first of these is that in GF(2) there are vectors of elements which are self-orthogonal, i.e. $\langle v_i,v_i \rangle = 0$ when $v_i \neq 0$. To illustrate, consider the following:
$$\begin{bmatrix} 1 \\ 1 \end{bmatrix} \cdot \begin{bmatrix} 1 \\ 1 \end{bmatrix} = 1 \times 1 + 1 \times 1 = 1 + 1 = 0$$
In addition, for a random $n \times n$ matrix of elements, the probability that the matrix is non-singular (has a determinant = 1) is given by [24]:
$$\prod_{i=1}^{i=n} 1 - 2^{i-1-n}$$
Taking the limit as $n\rightarrow \infty$, the probability of finding a non-singular matrix at random tends to $0.28878$ [24]. This contrasts with the case for random matrices with elements in $\mathbb{R}$, where the probability of finding a singular matrix tends to $0$; high dimensional random vectors become approximately orthogonal. Thus, the problem of finding a non-singular matrix is comparatively harder for GF(2) than for $\mathbb{R}$. However, it may still be possible to find a non-singular matrix within a tractable number of attempts.

Low discrepancy sequences (LDS), unlike random numbers, are deterministic and thus unable to provide a standard error on the estimate due to the non-applicability of the CLT. This is problematic as we are then unable to quantify the uncertainty around the estimate. Additionally, problems might arise when a `bad` number of points in a given sequence is used. While we can simply add a random number to each point of the sequence, this would increase the discrepancy of the resulting sequence. Hence, the challenge is to restore `randomness` (for the purposes of applying the CLT) without sacrificing the low discrepancy of the sequences (and by proxy - the convergence rate of the estimator). This can be achieved with `Randomized Quasi Monte-Carlo` where we apply a `scrambling` operation to the low discrepancy sequences. A good scrambling operation (defined by the `seed` parameter) will map the input LDS to another LDS with sufficient `randomness` to allow the computation of a standard error estimate. This is because scrambling the sequence with different seeds allows the application of CLT. Moreover, a good scrambling algorithm may also help reduce the non-uniformity of the Sobol sequences.

The simplest method for scrambling a low discrepancy is simply adding a uniform random number separately to each of the d-dimensions. This is known as a Cranley-Patternson (CP) rotation [23].
$$x_{ij} = z_{ij} + u_j\; mod\; 1$$
The rotation preserves the low discrepancy nature of a sequence (Equation \ref{eq:cranley_patterson_discrepancy}), but does not eliminate pathological non-uniform regions in the sequences; it merely shifts them. As the shifts are continuous, the CP rotation also does not preserve the stratification properties of the digital sequences achieved by subdividing the unit hypercube.
$$D_n(x_1,...,x_n) \leq 2^d D_n(a_1,...,a_n)$$

Digital scrambling is a similar operation to the CP rotation described above. However, here we work with the integer representation of the low discrepancy sequences, performing instead a summation in base-b between values in each dimension of the sequence and the uniform number:
$$x_{ij} = z_{ij} \oplus u_j$$
Representing both the low discrepancy value $z$ and the scrambling value $u$ in base-b where an arbitrary value $y$ can be written as:
$$y = \sum_{k=0}^p y_k b^{-k-1}$$
The per-digit scrambling operation can then be expressed as:
$$x_k = z_k + u_k \; mod \; b$$
For the case of Sobol sequences, the base-2 operation summation is the XOR making it fast to compute on modern computers. Similar to the CP rotation, a key downside of digital shift scrambling is that they do not randomize the sequence sufficiently, and thus do not achieve a significant improvement of estimate variance properties.

A more powerful approach to randomization is Owen`s scrambling, also known as nested uniform scrambling. Like the previous scrambling methods, Owen scrambling applies to each dimension separately. Within each dimension the nested uniform scrambling can be viewed as an extension of the digital scrambling where the bits of a number are XOR-ed with bits of a scrambling sequences. However, unlike digital scrambling where the scrambling sequence is constant, the scrambling sequence is sampled in a path-dependent way from a set of random bits generated for the scrambling. To illustrate with an example, consider the following binary sequences of length $m=4$:
$$\begin{align} \color{blue}{z_1 = 0110} \\ \color{violet}{z_2 = 0010}\end{align}$$
Generation of scrambling sequences for these requires a $2^m - 1 = 15$ random values which can be interpreted as forming a binary tree, as shown in Figure [5]. Given this binary tree, we can generate scrambling sequences for $\color{blue}{z_1}$ and $\color{violet}{z_2}$ by traversing the binary tree, taking the left branch when we encounter a $0$ in the value of $\color{blue}{z_1}$ (or $\color{violet}{z_2}$) and right when we encounter a $1$. Examples of this sort of tree traversal are shown in Figure [6].

The scrambling sequences generated by traversing the binary tree for $\color{blue}{z_1}$ and $\color{violet}{z_2}$ are:
$$\begin{align} \color{blue}{u_1 = 0101}\\ \color{violet}{u_2 = 0111}\end{align}$$
The path through the binary tree used to obtain these scrambling sequences is shown in Figure [6] which highlights the nodes of the binary tree that make up the digits of each value of $u$. The edges of the binary tree are annotated with the digits of the scrambled value ($\color{blue}{z_1}$/$\color{violet}{z_2}$) to show how this path through the tree was obtained.

Having generated the scrambling sequences we can now calculate the resulting scrambled values using digital scrambling:
$$\begin{align} \color{blue}{x_1 = z_1 \oplus u_1 = 0110 \oplus 0101 = 0011}\\ \color{violet}{x_2 = z_2 \oplus u_2 = 0010 \oplus 0111 = 0101} \end{align}$$
Formally, nested uniform scrambling can be viewed as a path-dependent permutation of the intervals. In this way, the scrambling rearranges the sub-intervals, but does not change the number of points in each sub-interval, thus preserving the uniformity properties. In this context, the source of randomness is the binary tree - defining the permutations to be performed. A key limitation of Owen`s scrambling is that for a d-dimensional Sobol sequence with a scrambling depth of $m$ in base-b, the scrambling requires $O(d\times b^m)$ values in order to store the binary trees for each dimension. Thus, scrambling 32-bit Sobol sequences requires $\approx 512$ MiB/dimension. Practical techniques for reducing memory usage are to limit the scrambling depth $m$ to $\le 32$ bits and apply a digital scramble to the remaining digits. An alternative is to approximate a binary tree with a hash function to scramble the bits.

The hash-based Owen scrambling replaces an explicit binary tree used to calculate the scrambling bit with an implicit traversal - a hash function. Here, instead of traversing down the tree we note that, for a given bit, the previous bits form a path or address of the next random bit to XOR with the bit being scrambled. Thus, to scramble a given bit $a_i$ of $z = 0.z_1z_2z_3...z_k$, we use the preceding bits: $v_i = z_1z_2...z_{i-1}$ we can treat the binary tree as a scrambling function into which we pass this address. Replacing the binary tree by a hash-function we can get the effect - the next random bit. To emulate randomness, we replace the RNG used to generate the binary tree with a random seed, also passed into the hash function. Hence, to compute the random bit $u_i$ that is then XOR`ed with the input bit to generate the scrambling, we perform:
$$\begin{aligned} x_i &= z_i \oplus u_i \\ u_i &= hash(v_i,seed) = hash(z_1z_2...z_{i-1}, seed) \end{aligned}$$
A naive implementation would iteratively apply this formula to perform a nested scrambling. However, applying a per-bit hash function is computationally expensive. An alternative approach described in [8] and based on the Laine-Karras (LK) hash function [9]. This hash function approximates the set of bit-wise hash calculations and vector XOR operations in one go, improving efficiency. Finally, the LK-hash must be applied to Sobol values in reverse order starting from the least significant bit. With this, the LK-hashing becomes a good approximation to Owen scrambling, permuting the elementary intervals and retaining the good stratification properties \\cite{building_better_lk_hash}.

While scrambling acts on the values of a points in the sequence, shuffling acts on the index (the order), reordering them without changing the values. This approach, described in [8], is aimed solely at reducing the correlations between dimensions, as opposed to scrambling which aims to improve sequence properties. While [8] described a `nested uniform shuffle`, our implementation utilizes a random shuffling to remove the correlations between individual one-dimensional Sobol sequences within the multidimensional sequence. As the one-dimensional points remain unchanged by the shuffling, the distribution of the points remains unchanged. However, as shown in Figure [8], the shuffling degrades the scaling of discrepancy with increasing point count for higher dimensions. The relation, on a log-log scale, between the centered discrepancy and number of points, $n$, is described by the formula:
$$log_{10}(CD_n) = \alpha log_{10}(n) + \beta$$
Prior to shuffling, the gradient, $\alpha$, of the tkrg-a-ap5 and new-joe-kuo-6.21201 were $-1.63$ and $-1.67$ respectively. Meanwhile after shuffling, the new values of $\alpha$ are $-1.03$ and $-1.05$, significantly closer to that of the pcg64 generator, at $-1.01$. This acts to illustrate that by shuffling the sequence we can decorrelate the individual dimensions, but at the cost of reducing the high dimensional equidistribution properties. An example of how scrambling and shuffling influence the two dimensional projections is illustrated in Figure [9] which illustrates that we can decorrelate individual dimensions by shuffling and scrambling, and thus achieve both randomization and decorrelation by combining the techniques.

To evaluate the quality of the directed numbers generated by the search procedure outlined above, we compare the performance of the generated Sobol sequences with the new-joe-kuo-6.21201 direction numbers described in [22],[4], using pseudo-random sequences generated by the pcg64 [10] as a baseline. In the following section, we present the key findings of our benchmarking, while the full set of test integrals and discrepancy measures is included in Appendix [C]. We note that low discrepancy sequences are expected to outperform pseudo-random numbers in the case when the test problem has a low `effective` dimension (distinguishing it from the `nominal dimension` of a problem - the number of terms in the integral). Thus, we plot a mean dimension of the test problem alongside the test metric where applicable to allow the reader to verify the suitability of applying low discrepancy sequences. Discussion on mean dimension alongside the relevant formulae are presented in the Appendix [A].

The approach employed by Joe and Kuo for searching for new direction numbers accounts for the fact that property A is not sufficient to prevent correlations between pairs of dimensions (we denote dimension as $d$). Joe and Kuo, therefore, take advance of the fact that a Sobol sequence can be viewed as a $(t, d)$-sequence in base 2, where every block of $2^m$ points forms a $(t, m, d)$-net. Hence, the $t$-value, given by ,$t(j, d; m) = \sum^d_{j=1}(s_j - 1)$, can measure the uniformity of point sets; smaller $t$-values lead to finer partitions, and hence more uniformly distributed points. Using the aforementioned concepts, Joe and Kuo adopted the following method:

• Search for direction numbers that satisfy property A for dimension 1,111 by assigning polynomials for each dimension up to degree 13.
  

• After 1,111 dimensions, search for direction numbers that optimize the t-values of the two-dimensional projections, reaching 21,201 dimensions.

Joe and Kuo made the following assumptions when defining their search criterion: direction numbers were already chosen up to dimension $d-1$ and that $m$ is given and fixed. They would thus aim to choose direction numbers for dimension $d$ so as to minimize the t-value.

Similarly to [1], our search strategy to generate tkrg-a-ap5 was to search for Sobol sequences satisfying the following properties:

• Property A for all 195,000 dimensions ($A_{k=195,000}$)
  

• Property A` for 5 adjacent dimensions ($A`_{k=5}$)

The reason for choosing to search only for $A_{k=195,000}$ is Lemma 2 of [1] which states that it is not possible to construct a Sobol sequence for $d\geq 4$ satisfying $A_k$ for all $k\leq d$. In order to obtain Sobol sequences which satisfy properties A ($A_{k=195,000}$) and A` ($A`_{k=5}$), we needed to find sets of direction numbers whose digits satisfy the equations \ref{eq:property_A} and \ref{eq:property_Ap}. In [1], the suggested proposal appeared to formulate this as a set of constraints, in effect transforming this search into a satisfiability problem - solvable by automated theorem provers. A potential benefit of this approach is to generate all possible direction numbers satisfying the given constraints. However, such boolean satisfiability is in general NP-hard. Instead of this the key to our approach is the efficient calculation of the determinants of high-dimensional matrices in the Galois field GF(2) (Section [Determinants in GF(2)]).
One of the factors that could affect the uniformity properties of tkrg-a-ap5 is the correlation between dimensions. To investigate the effect of the numper of sampled points in the sequence, we compare the correlation coefficient histograms in two dimenions for the first 1,000 ([10]a) and 10,000 ([10]b) dimensions of the new-joe-kuo-6.211201 and the tkrg-a-ap5 sequences. We show the heat maps of both sequences in Figure [10]. In both cases, we observe that as the number of points `n` increase, the distribution of correlation coefficients shifts towards zero. This is most apparent for the new-joe-kuo-6.211201 direction numbers in the case of 1,000 dimensions. The figures also help to illustrate that choosing an `n` arbitrarily, instead an integer power of 2, produced unbalanced point sets - increasing the number of highly correlated dimensions unnecessarily.

(a) A logarithmic-scale heatmap of correlation coefficients for 1,000 dimensions.

(b) A logarithmic-scale heatmap of correlation coefficients for 10,000 dimensions.

Figure 10: A logarithmic-scale (log10) heatmap of correlation coefficients for 1,000 (a) and 10,000 (b) dimensions is plotted for the tkrg-a-ap5 and new-joe-kuo-6.21201 sequences. The heatmap illustrates that increasing the number of points reduces the the number dimension pairs with a high correlation coefficient. For the first 1,000 dimensions we see that the new-joe-kuo-6.21201 sequence exhibits better correlation decay, while for the 10,000 dimension case - tkrg-a-ap5 has an advantage, particularly for high ($\geq 0.7$) correlation values. The plot also underscores the need to choose the number of points `n` to be a power of 2 - otherwise the the `unbalanced` point sets produces high correlation between dimensions.

Here, we examine the spurious variance of a multidimensional Sobol sequence. This quantity is a measure of the the average pairwise correlation between dimensions of a Sobol sequence (see Appendix Section [B] for full description) and is given by:
$$\rho_s = 1 + \frac{1}{d}\sum_{i=1}^{d}\sum_{j \neq i}^d \rho_{ij}$$
Plotting $\rho_S - 1$ against dimensions $d$ in Figure [11]a, we observe that the tkrg-a-ap5 direction numbers exhibit substantially less deviation from zero than the new-joe-kuo-6.21201 sequences. This may be attributed to the increased uniformity of the tkrg-a-ap5 sequence that gives rise to a more symmetric distribution of correlation coefficients, in turn resulting in the maximum deviation of the spurious variance being smaller. In contrast, the new-joe-kuo-6.21201 sequences have consecutive sets of correlation coefficients whose distribution is highly asymmetric (around zero), thus resulting in the spurious variance drifting away from zero until approximately $d\approx1,000$. This behaviour of having highly biased distributions of correlation coefficients between dimensions appears to be a feature of the new-joe-kuo-6.21201, as other sequences such as sobol-levitan-lemieux or jaeckel do not exhibit such deviations; such a structure may have arisen as a result of the employed direction number search strategy used in [4] to generate the sequence.
An effective strategy to decorrelate and make symmetric the correlation coefficient distribution is to scramble the underlying low discrepancy sequence, before recomputing the spurious variance for each scrambled sequence. The effect of this is shown in Figure [11]b where the mean, obtained by averaging 30 spurious variance traces, of both the tkrg-a-ap5 and new-joe-kuo-6.21201 are shown - with the latter no longer exhibiting the large drift away from zero, present in Figure [11]a.

A second method for assessing quality of the low discrepancy sequences is to evaluate numerical integrals, whose value can be calculated analytically, and comparing the results obtained (see Appendix [C] for a full comparison of test integrals). Here, one considers the absolute error between the numeric and analytic result against the number of points used, as well as other parameters like dimension of the integral or other terms in its definition that can affect its difficulty. Unlike with Monte-Carlo, in the case of QMC we cannot rely on the central limit theorem, and thus must know analytic results to measure convergence. Using randomization (RQMC), we can attempt to recover the ability to estimate the error using the CLT. However, we must compute multiple scrambled versions of the sequence - increasing the computation time.
Consider first the integral $I_{1}$ below. Its integrand is a product of individual dimensions $x_i$, with $d$ being the nominal integral dimension and $c_i$ being a constant term, parameterized by the index of the dimension $c_i = i^k$. We can scale the integrals difficulty, its effective dimension (see Appendix [A]), by setting $c_i$ which determines the size of the denominator relative to the numerator of each of the product terms.
$$\int_{[0,1]^d} \prod_{i=1}^d \frac{|4x_i - 2 | + c_i}{1+c_i} dx_i$$
In [21], the authors discuss three variants for $c_i$, a constant $c_i=1$, a quadratic $c_i = i^2$ and a partial power dependence $c_i = i^{1/3}$. In all cases the analytic solution is $I_{1}=1$. However, the integrals have different effective dimensions, with the quadratic scaling $c_i = i^2$ having the lowest effective dimension for a given actual dimension $d$. In Figure [12], we consider the first case, with $c_i = 1$ for $d=50$, where we show that the error of QMC does not appear to decay at a faster rate that of MC, represented by the pcg64 line. This may be attributed to the form of the product terms, where each term contributes equally to the effective dimension - albeit with a scaling coefficient of $\sim 0.06$, meaning that each nominal dimension contributes only slightly to the effective dimensionality. Of the three cases considered, convergence benefits of low discrepancy sequences for test integral $I_1$ are observed only for the case of $c_i = i^2$, where the effective dimension tends to $\sim 1.016$ as the nominal dimension $d \rightarrow \infty$.

(a) Absolute numeric integration error for $I_1$.

(b) Mean effective (superposition) dimension for $I_1$.

Figure 12: a) Shows the absolute error of the numeric estimate for $I_1$ with $c_i = 1$ for all $i$ with $d=50$. We see that for this integral there is little benefit to using Sobol sequences as the pseudo-random pcg64 appear to converge at the same rate as both of the low discrepancy sequences b) The effective dimension (multiplicative superposition sense) for the integral $I_1$ with $c_i = 1$. The effective and nominal dimension relationship is approximated linear with a gradient of $0.06$, thus every nominal dimension contributes approximately $0.06$ to the mean effective dimension.

A second test integral $I_2$, also considered in [21], is defined as:
$$\int_{[0,1]^d} \prod_{i=1}^d\frac{i + 2x_i}{i + 1} dx_i$$
Here, the decay of each nominal dimension to the mean effective dimension (in the multiplicative superposition sense) is similar to the case of $I_1$ with $c_i = i$ thus the contribution of successive nominal dimensions decays to zero with increasing dimension index. Consequently the effective dimension for $I_2$ tends to an asymptotic value of $\sim 1.084$ as $d\rightarrow \infty$, as shown in Figure [13b]. The result of this small value of effective dimension gives rise to the fast onset of the $O(n^{-1})$ convergence rate, shown in Figure [13a].

(a) Absolute numeric integration error for $I_2$.

(b) Mean effective (superposition) dimension for $I_2$.

Figure 13: a) Shows the absolute error of the numeric estimate for integral $I_2$. b) The effective dimension (multiplicative superposition) for the integral $I_2$ reaches $\sim 1.08$ in the limit as $d\rightarrow \infty$ with higher dimensions being dominated by index $i$ in both the numerator and denominator - resulting in an almost constant value of higher order terms in the product.}

Figure [14] shows plots of centred discrepancy versus number of points for the tkrg-a-ap5, new-joe-kuo-6.21201 and pseudorandom sequences at dimensions 5, 10 and 30. For all dimensions and sample sizes, the performance of tkrg-a-ap5 and new-joe-kuo-6.21201 are similar. The pseudorandom sequences are the worst performing on dimensions 5 and 10, while the convergence rate of the pseudorandom sequence at dimension 20 is similar to that of tkrg-a-ap5 and new-joe-kuo-6.21201.

(a) Centred discrepancy for $d=5$.

(b) Centred discrepancy for $d=10$.

(c) Centred discrepancy for $d=30$.

Figure 14: Centred discrepancy for dimensions a) d=5 b) d=10 c) d=30. In all cases the plots show an expected reduction of discrepancy with increasing point count with the pseudo-random pcg64 performing worse that the low discrepancy sequences for d=5 and d=10. However for d = 30 this difference appears smaller. This may be attributed to the difficulty in estimating the discrepancy as well as the number of points required to see improvements in discrepancy increases with increasing dimension, as described by Owen.

In this report, we reviewed the basic theory of low discrepancy sequences and the conditions under which such they can be applied to Monte-Carlo numeric integration. When applied to suitable problems, the Quasi-Monte Carlo methods estimator convergence rates of $O(n^{-1})$, converging faster than conventional Monte-Carlo which exhibits an $O(n^{-1/2})$ rate. We then considered Sobol sequences which are base-2 low discrepancy sequences, defined using recurrence relations from an initial set of direction numbers and primitive polynomials. We showed that when sampling, it is critical to keep the number of points in a sequence a power of 2, required to maintain `balanced` sequences, in maintaining the improved convergence rate and minimizing bad two-dimensional projections quantified by computing correlation coefficients. We then presented properties A and A`, relating them to properties of the direction numbers and outlined a strategy to search for such direction numbers by performing determinant calculations in GF(2) - allowing us to generate the tkrg-a-ap5 sequence. Following this, we covered Randomized QMC introducing three different scrambling techniques followed by presenting a comprehensive set of benchmarks for evaluating LDS sequences in both low and high dimension.

While the `nominal` dimension of $\Omega$ (the integration domain) can give an indication of the difficulty of the integration problem, some integration problems may have a `effective` dimension that is significantly lower this nominal value. In this case, the majority of the contribution lies in a lower dimensional manifold, sampling from which we can converge to the target solution at a rate substantially faster than predicted by the nominal dimension. This class of problems is the one most suitable for QMC integration. Calculating an `effective` dimension in practice is challenging, but analytic equations exist for integrands that have a certain functional forms. Two forms, as discussed in [7] are given below. Example of function with these forms are shown in [3].
In the subsequent sections, we use the following terms:

• $s$: the dimensionality of the function;
  

• $g_j$: a subfunction applied the $j$th dimension;
  

• $x_j$: the $j$th element of a given random sample.

Integrands in the additive form can be represented as a sum of simpler terms:
$$f(x) = \sum^s_{j=1}g_j(x_j)$$
The mean effective dimension of the additive function in the superposition sense is one:
$$d_{sup} = 1$$
The mean effective dimension of the additive function in the truncation sense is given by:
$$d_{trunc} = \frac{1}{\sigma^2}\sum^s_{j=1}j\gamma^2_j$$

Integrands in the multiplicative form can be represented as a product of simpler terms as follows:
$$f(x) = \prod^s_{j=1}g_j(x^j)$$
The mean effective dimension of the multiplicative function in the superposition sense is:
$$d_{sup} = \frac{\sum^s_{j=1}\gamma^2_j/(\gamma^2_j + \mu^2_j)}{1 - \prod^s_{j=1}\mu^2_j/(\gamma^2_j + \mu^2_j)}$$
The mean effective dimension in the truncation sense is not defined.

In this section, we describe the `spurious variance` test used to evaluate quality of the low discrepancy sequences. In it we consider independent and identically distributed random variables $X_i \sim Normal(0,1)$. Performing a summation over $d$ sum variables, we obtain a new random variable $Y$:
$$Y = \sum_{i=1}^d X_i$$
In the case when the random variables $X_i$ satisfy the IID assumption, the resulting variance of $Y$ reduces to:
$$Var(Y) = \sum_{i=1}^d Var(X_i) = d$$
While this equality will not hold exactly for finite sample sizes, it should hold approximately and can thus be used to check whether the samples generated are approximately independent and identically distributed. These assumptions may not hold in practice when the underlying number generators exhibit pairwise correlations between dimensions of the random numbers. Such correlations will give rise to an additional unwanted or `spurious` variance term - causing the sum of the random variable samples to deviate from the theoretically predicted value. Thus, by calculating this spurious variance term, we can assess the quality of multidimensional low discrepancy sequence due to the one-to-one mapping between random variables $X_i$ and the dimensions $i$ of the low discrepancy sequence.

Generalizing the above two variable sum to $d$ variables, and applying the bilinearity of covariance:
$$Cov\left[\sum_{i}a_iX_i, \sum_j b_j Z_j\right] = \sum_i \sum_j a_i b_j Cov[X_i,Z_j]$$
where $a_i$ and $b_j$ are scaling factors, while $X_i$ and $Z_i$ are random variables drawn from a Gaussian distribution.
We can obtain the following expression for the variance of the sum $Y$:
$$Var(Y) = d + \sum_{i=1}^{d}\sum_{j \neq i}^d \rho_{ij}$$
Here, the second term is the pairwise sum of correlation coefficients running over all unequal indices $j\neq i$, where again we have applied the reductions $\sigma_i = 1$ for all $i\in [1,\cdots, d]$. In the case when the individual random variables $X_i$ are independent and identically distributed the correlations, $ho_{ij} = 0$ for all $i \neq j$. However, when the two dimension projections of the sequences exhibit non-zero correlations, the resulting correlations directly impact the variance $Var(Y)$, as $Cov(x_i, x_j)\neq0$. Thus, by computing variance of $Y$ and measuring its deviation from the ideal case when $\rho_{ij} = 0$ for $i \neq j$, we can detect problematic two dimensional projections in a d-dimensional sequence.

One way of comparing low discrepancy sequences is by estimating the error of various integrals and comparing these estimates to the analytic solutions of the same integrals. Given that the analytic solution of the test integral is known, the approximation errors and its rate of convergence (or lack thereof) can be used to assess the quality of the low discrepancy sequence, or other MC/QMC approach. In general the rate of convergence will depend both the sample generation method and the integral being estimated. By plotting the estimate error with increasing number of samples on a log-log plot and extracting the gradient one can then estimate the convergence rate. This is the topic of this appendix

• Section [Classes of test integrals]: presents a table of test integrals with additive and multiplicative form (see Appendix [A]).
  

• Section [Non-Decomposable test integrals]: presents several special integrals that have neither additive or multiplicative form.
  

• Section [Experiments]: presents experimental results on some of the additive and multiplicative alongside a discussion.

The first set of integrals we examine are those that can be categorised as either additive or multiplicative. These categories refer to the means of decomposing the integrals: either as a summation or a product of several sub-integrals (one sub-integral per dimension). We provide the mathematical formulations of each category, as well as possible instantiations of the two classes in the following two tables. these integrals have analytic solutions against which the estimation error can be calculated. The sources of the test integrals are given in the right-most column of the second table. Additionally, $d$ is defined as the number of dimensions, $g_j$ is the sub-integrand for dimension $j$ and $c(d)$ is the scaling factor.