sobolset

Sobol quasirandom point set

Description

sobolset is a quasirandom point set object that produces points from the Sobol sequence. The Sobol sequence is a base-2 digital sequence that fills space in a highly uniform manner.

Creation

Syntax

p = sobolset(d)

p = sobolset(d,Name,Value)

Description

p = sobolset(d) constructs a d-dimensional point set p, which is a sobolset object with default property settings. The input argument d corresponds to the Dimensions property of p.

p = sobolset(d,Name,Value) sets properties of p using one or more name-value pair arguments. Enclose each property name in quotes. For example, sobolset(5,'Leap',2) creates a five-dimensional point set from the first point, fourth point, seventh point, tenth point, and so on.

The returned object p encapsulates properties of a Sobol quasirandom sequence. The point set is finite, with a length determined by the Skip and Leap properties and by limits on the size of the point set indices (maximum value of 2⁵³). Values of the point set are generated whenever you access p using net or parenthesis indexing. Values are not stored within p.

Properties

Dimensions — Number of dimensions

Read-only: positive integer scalar in interval [1,1111]

This property is read-only.

Number of dimensions of the points in the point set, specified as a positive integer scalar in the interval [1,1111]. For example, each point in the point set p with p.Dimensions = 5 has five values.

Use the d input argument to specify the number of dimensions when you create a point set using the sobolset function.

Use the reduceDimensions object function to reduce the number of dimensions after you create a point set.

Leap — Interval between points

0 (default) | positive integer scalar

Interval between points in the sequence, specified as a positive integer scalar. In other words, the Leap property of a point set specifies the number of points in the sequence to leap over and omit for every point taken. The default Leap value is 0, which corresponds to taking every point from the sequence.

Leaping is a technique used to improve the quality of a point set. However, you must choose the Leap values with care. Many Leap values create sequences that fail to touch on large sub-hyper-rectangles of the unit hypercube and, therefore, fail to be a uniform quasirandom point set. For more information, see [4].

Example: p = sobolset(__,'Leap',50);

Example: p.Leap = 100;

PointOrder — Point generation method

'standard' (default) | 'graycode'

Point generation method, specified as 'standard' or 'graycode'. The PointOrder property specifies the order in which the Sobol sequence points are produced. When PointOrder is set to 'standard', the points produced match the original Sobol sequence implementation. When PointOrder is set to 'graycode', the sequence is generated by an implementation that uses the Gray code of the index instead of the index itself.

You can use the 'graycode' option for faster sequence generation, but the software then changes the order of the generated points. For more information on the Gray code implementation, see [1].

Example: p = sobolset(__,'PointOrder','graycode');

Example: p.PointOrder = 'standard';

ScrambleMethod — Settings that control scrambling

0x0 structure (default) | structure with Type and Options fields

Settings that control the scrambling of the sequence, specified as a structure with these fields:

• Type — A character vector containing the name of the scramble

• Options — A cell array of parameter values for the scramble

Use the scramble object function to set scrambles. For a list of valid scramble types, see the type input argument of scramble. An error occurs if you set an invalid scramble type for a given point set.

The ScrambleMethod property also accepts an empty matrix as a value. The software then clears all scrambling and sets the property to contain a 0x0 structure.

Skip — Number of initial points in sequence to omit

0 (default) | positive integer scalar

Number of initial points in the sequence to omit from the point set, specified as a positive integer scalar.

Initial points of a sequence sometimes exhibit undesirable properties. For example, the first point is often (0,0,0,...), which can cause the sequence to be unbalanced because the counterpart of the point, (1,1,1,...), never appears. Also, initial points often exhibit correlations among different dimensions, and these correlations disappear later in the sequence.

Example: p = sobolset(__,'Skip',2e3);

Example: p.Skip = 1e3;

Type — Sequence type

Read-only: 'Sobol' (default)

This property is read-only.

Sequence type on which the quasirandom point set p is based, specified as 'Sobol'.

Object Functions

net	Generate quasirandom point set
reduceDimensions	Reduce dimensions of Sobol point set
scramble	Scramble quasirandom point set

You can also use the following MATLAB^® functions with a sobolset object. The software treats the point set object like a matrix of multidimensional points.

length	Length of largest array dimension
size	Array size

Examples

Create Sobol Point Set

Generate a three-dimensional Sobol point set, skip the first 1000 values, and then retain every 101st point.

p = sobolset(3,'Skip',1e3,'Leap',1e2)

p = 
Sobol point set in 3 dimensions (89180190640991 points)

Properties:
              Skip : 1000
              Leap : 100
    ScrambleMethod : none
        PointOrder : standard

Apply a random linear scramble combined with a random digital shift by using scramble.

p = scramble(p,'MatousekAffineOwen')

p = 
Sobol point set in 3 dimensions (89180190640991 points)

Properties:
              Skip : 1000
              Leap : 100
    ScrambleMethod : MatousekAffineOwen
        PointOrder : standard

Generate the first four points by using net.

X0 = net(p,4)

X0 = 4×3

    0.7601    0.5919    0.9529
    0.1795    0.0856    0.0491
    0.5488    0.0785    0.8483
    0.3882    0.8771    0.8755

Generate every third point, up to the eleventh point, by using parenthesis indexing.

X = p(1:3:11,:)

X = 4×3

    0.7601    0.5919    0.9529
    0.3882    0.8771    0.8755
    0.6905    0.4951    0.8464
    0.1955    0.5679    0.3192

Tips

• The Skip and Leap properties are useful for parallel applications. For example, if you have a Parallel Computing Toolbox™ license, you can partition a sequence of points across N different workers by using the function spmdIndex (Parallel Computing Toolbox). On each nth worker, set the Skip property of the point set to n – 1 and the Leap property to N – 1. The following code shows how to partition a sequence across three workers.

  Nworkers = 3;
  p = sobolset(10,'Leap',Nworkers-1);
  spmd(Nworkers)
      p.Skip = spmdIndex - 1;
  
      % Compute something using points 1,4,7...
      % or points 2,5,8... or points 3,6,9...
  end

Algorithms

Sobol Sequence Generation

Consider a default sobolset object p that contains d-dimensional points. Each p(i,:) is a point in a Sobol sequence. The jth coordinate of the ith point, p(i,j), is equal to

$$\{\begin{array}{cc}0,& i=1\\ {\gamma }_i(1)v_j(1)\oplus {\gamma }_i(2)v_j(2)\oplus \mathrm{...},& i>1.\end{array}$$

• The $${\gamma }_i(n)$$ values are 0s or 1s such that

  $$i-1={\displaystyle \sum _{n=1}{\gamma }_i(n)2^{n-1}}.$$

  In other words, the $${\gamma }_i(n)$$ values are the binary digits of the integer i – 1.

• The v_j(n) values are called direction numbers. They are uniquely defined for each coordinate j. For more details on these values, see Direction Numbers Generation.

• The ⊕ operator is the bitwise exclusive-or operator. For two numbers expressed in binary, the ⊕ operator compares the digits in each position. For a given digit position, the ⊕ operator returns a 1 if the digits in that position differ and returns a 0 if the digits in that position are the same.

  • For example, $$19\oplus 24={(10011)}_2\oplus {(11000)}_2={(01011)}_2=11.$$

  • Similarly, $$\frac{1}{2}\oplus \frac{3}{4}={(0.1)}_2\oplus {(0.11)}_2={(0.01)}_2=\frac{1}{4}.$$

For more information, see [3].

Direction Numbers Generation

The set of direction numbers v_j(n) depends on the coordinate j. Define the direction numbers in terms of m_j(n) values:

$$v_j(n):=\frac{m_j(n)}{2^n}.$$

For each j, you can generate the direction numbers by selecting the following:

• A primitive polynomial in $${\mathbb{Z}}_2$$ of some degree s_j

  $$x^{s_j}+a_j(1)x^{s_j-1}+a_j(2)x^{s_j-2}+\mathrm{...}+a_j(s_j-1)x+1.$$

  Each coefficient in the polynomial is either 0 or 1.

• s_j initial direction numbers. For each initial direction number, the corresponding m_j(n) value must be either 1 or an odd number less than 2ⁿ.

The remaining direction numbers are determined by the following recurrence relation, which uses the coefficients of the primitive polynomial, the previous direction numbers, and the ⊕ bitwise exclusive-or operator.

$$m_j(n):=2a_j(1)m_j(n-1)\oplus 2^2{a}_j(2)m_j(n-2)\oplus \mathrm{...}\oplus 2^{s_j-1}{a}_j(s_j-1)m_j(n-s_j+1)\oplus 2^{s_j}{m}_j(n-s_j)\oplus m_j(n-s_j).$$

sobolset uses the same primitive polynomials and initial direction numbers described in [3]. These parameters are provided for the first 1111 dimensions.

References

[1] Bratley, P., and B. L. Fox. “Algorithm 659 Implementing Sobol's Quasirandom Sequence Generator.” ACM Transactions on Mathematical Software. Vol. 14, No. 1, 1988, pp. 88–100.

[2] Hong, H. S., and F. J. Hickernell. “Algorithm 823: Implementing Scrambled Digital Sequences.” ACM Transactions on Mathematical Software. Vol. 29, No. 2, 2003, pp. 95–109.

[3] Joe, S., and F. Y. Kuo. “Remark on Algorithm 659: Implementing Sobol's Quasirandom Sequence Generator.” ACM Transactions on Mathematical Software. Vol. 29, No. 1, 2003, pp. 49–57.

[4] Kocis, L., and W. J. Whiten. “Computational Investigations of Low-Discrepancy Sequences.” ACM Transactions on Mathematical Software. Vol. 23, No. 2, 1997, pp. 266–294.

[5] Matousek, J. “On the L2-Discrepancy for Anchored Boxes.” Journal of Complexity. Vol. 14, No. 4, 1998, pp. 527–556.

Version History

Introduced in R2008a