Adrian Saldanha

Sobol numbers: understanding and computing them

Nov 09, 2018 ยท

1 Introduction

The Sobol numbers are a sequence of \(k\)-tuples which are pseudo-random. They are very useful for quasi-monte carlo techniques as they provide a source of quasi-random numbers. There are many scattered pages around the internet already describing the Sobol numbers in detail. I will first write a practical guide to computing them, then explain the mathematics in some detail.

2 Computing the Sobol numbers

The Sobol numbers are a sequence of tuples of some length \(D\) in the real numbers.

To compute the \(d\)th dimension of the \(n\)th random number, we use the binary representation for \(n\) and take the xor of some other set of numbers, called our direction numbers, of which there is one set of such numbers per dimension. To compute the nth Sobol number, we take the xor of the direction numbers whose corresponding bit in n is set. This can be explained mathematically as follows.

Define an operator \(\bigoplus\) like \(\sum\), except using \(\oplus\) instead of \(+\):

$$\bigoplus_{k=s}^n x_k = x_s \oplus x_{s+1} \oplus ... \oplus x_n$$

Using the binary representation for \(n = n_0\cdot2^0 + n_1\cdot2^1 + \cdots + n_{\lfloor\log_2 n\rfloor}\cdot2^{\lfloor \log_2 n\rfloor}\), the \(n\)th Sobol number is:

$$s_n = \bigoplus_{i=0}^{\lfloor\log_2 n\rfloor} n_id_i$$

Direction numbers

For 32-bit integers there are 32 direction numbers for each dimension we want to compute.

To start, we are given the first \(s\) direction numbers. \(s\) is the degree of a polynomial with the \(s\)th coefficient multiplied by the direction number. To extend our \(s\) direction numbers to the bit size of our integers, we can get the next direction number by taking the xor of the \(s\) previous direction numbers, each multiplied by its coefficient in the polyomial for this dimension, multiplied by \(2^i\), where \(i\) is the index into the coefficient sequence.

If our initial direction numbers are some sequence \(d=d_1\ldots d_s\), for \(i \ge s\) the next direction number is:

$$d_{i+1} = d_{i-s} \oplus \bigoplus_{k=1}^{s} d_{i-k}\cdot a_{k} \cdot 2^k$$

The quality of the sequence depends on the direction numbers chosen for each dimension \(d\).

Luckily high-quality direciton numbers are available to us. To simplify things, I will provide all the direction numbers required to compute the first 21201 dimensions. These numbers are derived from the set provided over here: Sobol sequence generator. This takes thrity-two multiplies and xors to compute the nth number. While not terribly computationally expensive, there is a huge improvement which can be made.

Increasing efficiency

We can do better than thirty-two multiplies and xors if we know the previous sobol number. We can do it in two to six hardware instructions. We only need one xor as well as the find least significant zero bit instruction, which is often just a few hardware instructions. The mnemonic for this instruction is often ffz. See below for implementaiton of ffz in terms of the other bit scan operations like clz.

There is a re-ordering of the k-bit integers, \(g : \mathbb{Z}_{2^{k}}\to \mathbb{Z}_{2^{k}}\) called the Gray code. To simplify, we will include its definition later.

All the \(g(i)\) are unique, which means instead of \(n\) when computing \(s_n\) we can use \(g(n)\) to compute \(s_{g(n)}\) and get the Sobol numbers in a different order. This is just as inefficient as before, though. The main property we need of \(g(i)\) is that they differ in only one bit from the previous Gray code. Therefore, we simply need to xor against this bit's direction number:

$$ s_{i+1} = s_{i} \oplus d_{f(i)} $$

\(f(i)\) is the bit that will change from \(g(i)\) to \(g(i+1)\), which happens to be the least significant zero bit in the binary representation of \(i\).

We don't actually require it if computing numbers sequentially, but the function \(g(l)\) is called the Gray code of \(l\), and to compute it:

$$ g(l) = l \oplus {l \over 2} $$

The gray code can be used to compute a skipahead, so it can be extremely useful for a massively parallel implementation to skip ahead to the nth point, and then use the formula above.

3 Properties of the direction numbers required

The primitive polynomial is a polynomial in \(\mathbb{Z}_2\)--coefficients are zero or one--and the polynomial is primitive, meaning it cannot be factored in \(\mathbb{Z}_2\). For \(N\)-bit integers it defines the sequence \(a=a_1 \ldots a_{s}\), where \(s\) is the dimension of the polynomial:

$$ p(x) = x^s + a_1 \cdot x^{s-1} + \cdots + a_{s-1} \cdot x + a_s $$

where \(a_s = 1\)

4 Original motivation

Sobol initially was trying to compute some integral of \(f : I^n \to R\) over the hypercube \(I^n\):

$$ \int_{I^n} f d\mu $$

where \(\mu\) is the standard metric on \(\mathbb{R}^n\).

This can be done by using a monte-carlo estimator with a uniform distribution, using a set of size \(n\) with points \(P \in I^n\): $$ {1\over N} \sum_{x\in I^n} f(x) $$

Appendix A: Implementation of ffz in terms of other common operators

Sometimes ffz is not available directly in hardware, so we may want to use other common operations to compute it.

If we have a bitset \(011010101111\), we want to find the index from the right (least-sigificant bit) of the first zero, in this case it is \(4\).

On almost all architectures, we have the find most significant set bit operation. There is a table here with various options here. We will refer to this operation as clz().

If we have this operation, we can negate the integer bit-wise, replacing all zeros with ones and vice versa. Then we can clear all the ones except for the least sigificant one, ie. the rightmost.

This is (in c pseudocode):

 function  negate(uint32_t x)     { return ~x; }
 function  clear_ones(uint32_t x) { return  x & -x; }
 function  ffs(uint32_t x)        { return 31 - clz(clear_ones(x)); }
 function  ffz(uint32_t x)        { return 31 - 
                                      clz(clear_ones(negate(x))); }

Appendix B: Direction numbers

One can use the formulas above to compute the direction numbers for each dimension given the primitive polynomials.

However, it may be easier to simply load them from the disk. Here are all the direction numbers computed using the primitive polynomials available here.