.. _ProcSobolSequence:

Procedure: Generating a Sobol's Sequence
========================================

Description and Background
--------------------------

The numbers in a Sobol's sequence are generated directly as binary
fractions of length :math:`w` bits. These numbers are created from a
special set of :math:`w` binary fractions :math:`V_i` called direction
numbers. In Sobol's original method, the j-th number :math:`X_j` is
generated by the operation XOR among those direction numbers :math:`V_i`
such that the i-th bit of :math:`j` is not zero, see comments below. See
:ref:`DiscSobol<DiscSobol>` for additional discussion.

Detailed description of the procedure
-------------------------------------

The :math:`i`-th number in Sobol's sequence is

.. math::
   x_i=\bigoplus_{j=1}^\infty b_j(i)v_j,

where :math:`b_j(i)` is the :math:`j`-th digit of the binary expansion of
:math:`i`; :math:`v_j` is the :math:`j`-th direction number and :math:`\oplus` is
the binary (bitwise) XOR operation. The number :math:`b_j(i)` is a bit and
thus it only takes values :math:`0` or :math:`1`, i.e. for an integer :math:`i`
we have its binary expansion

.. math::
   i=\sum_{j=1}^\infty b_j(i)2^{j-1}.

As for any finite number :math:`i`, its binary expansion has a finite
number of non-zero digits, the Sobol's number is computed with a XOR sum
of a finite number of direction numbers and thus no convergence
properties are needed in the above equations.

An important element in Sobol's sequence are the direction numbers. The
direction numbers are binary (dyadic) fractions :math:`v_j\in(0,1)` and the
resulting Sobol's number is also a binary fraction :math:`x_i\in(0,1)`. The
direction numbers are represented as

.. math::
   v_j=\frac{m_j}{2^j}

where :math:`m_j` takes integer values, :math:`1<m_j<2^j`, i.e. the direction
number :math:`v_j` is a shift of the binary digits in :math:`m_j`. The
direction numbers are computed by the recursive relation

.. math::
   m_j=\left(\bigoplus_{k=1}^d2^ka_km_{j-k} \right)\bigoplus m_{j-d},

where the initial values :math:`m_1,\ldots,m_d` are odd integers and the
numbers :math:`a_0,a_1,\ldots,a_d` are given by the coefficients of
:math:`p_d(x)`, a primitive polynomial of degree :math:`d` in :math:`{Z}_2,`
which is written as

.. math::
   p_d(x)=\sum_{k=0}^da_kx^{d-k}.

Recall that a primitive polynomial is a polynomial that cannot be
represented as the product of two polynomials of lower degree. Primitive
polynomials are available in tables in the literature, see Press et al.
(1994). For a primitive polynomial, the numbers :math:`a_0` and :math:`a_d`
take the value one by definition, and :math:`a_0` is not used in the
computations.

Running example
---------------

In what follows, the subindex :math:`_{\textbf{2}}` refers to binary
numbers, everywhere else usual representation (base ten) is assumed.

Take the primitive polynomial :math:`p_3(x)=x^3+x+1`, and thus :math:`d=3` and
:math:`(a_0,a_1,a_2,a_3)=(1,0,1,1)`. The recurrence relation for the
direction numbers is

.. math::
   \begin{array}{rcl}m_j&=2a_1m_{j-1}\oplus 4a_2m_{j-2}\oplus
   8a_3m_{j-3}\oplus m_{j-3}\\&=4m_{j-2}\oplus 8m_{j-3}\oplus m_{j-3} \\
   \end{array}

Using a different primitive polynomial will give a different recurrence
relation. To iterate the recurrence relation, initial odd values
:math:`m_1,m_2,m_3` are needed. Take :math:`m_1=1,m_2=3,m_3=7.`

Now we compute a few extra numbers :math:`m_j`. We have

.. math::
   \begin{array}{rcl}m_4&=4(3)\oplus 8(1)\oplus 1=12\oplus 8\oplus
   1 \\ &=1100_{\textbf{2}}\oplus 1000_{\textbf{2}}\oplus
   1_{\textbf{2}}=101_{\textbf{2}}=5 \end{array}

and

.. math::
   m_5=28\oplus 24\oplus 3=111_{\textbf{2}}=7.

The direction numbers are computed by shifting the values obtained, e.g.

.. math::
   v_1 &= \frac{1}{2}=0.1_{\textbf{2}} \\
   v_2 &= \frac{3}{4}=0.11_{\textbf{2}} \\
   v_3 &= \frac{7}{8}=0.111_{\textbf{2}} \\
   v_4 &= \frac{5}{16}=0.0101_{\textbf{2}}

and

.. math::
   v_5=\frac{7}{32}=0.00111_{\textbf{2}}.

The first Sobol's number is obtained with the first direction number

.. math::
   x_1=1(v_1)=v_1=0.1_{\textbf{2}}=0.5.

As the binary expansion of two is :math:`10_{\textbf{2}}` then

.. math::
   x_2=1(v_2)\oplus 0(v_1) =v_2=0.11_{\textbf{2}}=0.75

and

.. math::
   x_3=v_1\oplus v_2=0.01_{\textbf{2}}=0.25.

The following table summarises the construction performed, with the
Sobol's sequence on the two rightmost columns, the next to last in binary
and the last in base ten.

================ =========== =============== =========== =========== =========== ===============
:math:`i,j (10)` :math:`i,j` :math:`m_j(10)` :math:`m_j` :math:`v_j` :math:`x_i` :math:`x_i(10)`
================ =========== =============== =========== =========== =========== ===============
1                1           1               1           .1          .1          .5
2                10          3               11          .11         .11         .75
3                11          7               111         .111        .01         .25
4                100         5               101         .0101       .111        .875
5                101         7               111         .00111      .011        .375
6                110                                     .001        .125
7                111                                     .101        .625
8                1000                                    .0101       .3125
9                1001                                    .1101       .8125
10               1010                                    .1001       .5625
================ =========== =============== =========== =========== =========== ===============

Additional Comments, References and Links
-----------------------------------------

It is important to note that for the first Sobol's number, only the first
direction number was needed; then for the following two Sobol's numbers
the second direction number was included; for the following four Sobol',
the third direction number was included. Without iterating again the
recursive relation, :math:`31` Sobol's numbers can be constructed using the
first five direction numbers in the above table.

In short, to construct the first :math:`2^k-1` Sobol's numbers, we need
:math:`k` direction numbers. If more Sobol's numbers are needed, then the
recursive equation must be iterated to obtain direction numbers as
required.

Note that by selecting odd initial values :math:`m_1,\ldots,m_d`, all the
subsequent :math:`m_{d+1},m_{d+2},\ldots` are guaranteed to be odd numbers
and thus the :math:`i`-th bit of the direction number :math:`v_i` is one. This
has the important consequence of allowing Sobol's numbers to lie in
consecutive finer binary meshes. In other words, a latin hypercube is
constructed with the first :math:`2^k-1` Sobol's points.

The construction of multivariate Sobol's sequences is achieved by using
different primitive polynomials for each dimension. For a table with
different primitive polynomials see Press et al. (1994). Sobol's gave a
list of recommended primitive polynomials, to avoid high correlations
between different dimensions.

An alternative version of Sobol's sequence was due to Antonov and
Saleev, who proposed taking instead

.. math::
   x_i=\bigoplus_{j=1}^\infty g_j(i)v_j

where :math:`g_j(i)` is the :math:`j`-th digit of the Gray code representation
of :math:`i`. This different Sobol's proposal is faster than the original,
as it simplifies the computation to :math:`x_{i+1}=x_i\oplus v_c,` where
:math:`b_c` is the rightmost zero bit in the representation of :math:`i`.

References:

Antonov, Saleev (1979). USSR Comput. Math. Math. Phys. 19, 252-256.

Bratley, Fox (1988), ACM Trans. Math. Soft. 14(1), 88-100.

Press et al. (1994). Numerical Recipes in C, Cambridge.