Quasi-Random Sequences

Quasi-random sequences, also called low discrepancy sequences, are sequences of vectors that progressively cover a multi-dimensional space with points that are uniformly distributed. They are used in quasi-Monte Carlo simulations, in the evaluation of multi-dimensional integrals, and in global optimization.

There are several types of quasi-random sequences. They all share the property that they are less random than a sequence of vectors generated from random number generators. However, this is often an advantage for tasks such as approximating integrals in higher dimensions, and in global optimization. This is because low discrepancy sequences tend to sample space "more uniformly" than random numbers. Algorithms that use such sequences may have superior convergence over algorithms that use completely random vectors.

The QuasiRandom class implements static methods that return low discrepancy sequences.

The best known algorithm for generating quasi-random sequences is due to the French mathematician Henri Fauré. Fauré sequences are popular in mathematical finance simulations.

The FaureSequence method returns a Fauré sequence. The method takes two arguments. The first is the dimension of the sequence. This is the number of elements of the vectors in the sequence. The maximum value for the dimension is 40. The second argument is the length of the sequence. The method returns an IEnumerable<T> of Vector<T> of real numbers. The sequence always returns the same vector object, which is a mutable read-only vector. The elements of the vector change in each iteration, but they cannot be modified.

In the example below, we use a Fauré sequence approximate an 8-dimensional integral:

double sum = 0.0;
int length = 10000;
foreach (var v in QuasiRandom.FaureSequence(8, length))
{
    // Evaluate the integrand:
    double functionValue = 1.0;
    for (int i = 0; i < v.Length; i++)
        functionValue *= Math.Abs(4.0 * v[i] - 2.0);
    sum += functionValue;
}
Console.WriteLine("Estimate: {0,8:F4}", sum / length);

Dim sum = 0.0
Dim length = 10000
For Each v In QuasiRandom.FaureSequence(8, length)
    ' Evaluate the integrand
    Dim functionValue = 1.0
    For i = 0 To v.Length - 1
        functionValue *= Math.Abs(4.0 * v(i) - 2.0)
    Next
    sum += functionValue
Next
Console.WriteLine("Estimate: { 0,8:F4 }", sum / length)

No code example is currently available or this language may not be supported.

let mutable sum = 0.0
let length = 10000
for v in QuasiRandom.FaureSequence(8, length) do
    // Evaluate the integrand:
    let mutable functionValue = 1.0
    for i in 0 .. (v.Length - 1) do
        functionValue <- functionValue * abs(4.0 * v.[i] - 2.0)
    sum <- sum + functionValue
printfn "Estimate: %8.4f" (sum / float length)
//

A Halton sequence in n dimensions actually consists of n separate sequences of numbers. Every sequence is generated from a different prime number. It has been used successfully for numerical integration in multiple dimensions.

A Halton sequence works very well for low-dimensional spaces. For higher-dimensional spaces (20 or more dimensions) the initial elements can be very poorly distributed over the space. Moreover, some correlation patterns have been observed for larger dimensions involving larger prime numbers, making the Halton sequence relatively unsuitable for higher dimensions.

The HaltonSequence method returns a Halton sequence. The method takes two arguments. The first is the dimension of the sequence. This is the number of elements of the vectors in the sequence. The maximum value for the dimension is 40. The second argument is the length of the sequence. The method returns an IEnumerable<T> of Vector<T> of real numbers. The sequence always returns the same vector object, which is a mutable read-only vector. The elements of the vector change in each iteration, but they cannot be modified.

In the example below, we use a Halton sequence approximate the same 8-dimensional integral as above:

double sum = 0.0;
int length = 10000;
foreach (var v in QuasiRandom.FaureSequence(8, length))
{
    // Evaluate the integrand:
    double functionValue = 1.0;
    for (int i = 0; i < v.Length; i++)
        functionValue *= Math.Abs(4.0 * v[i] - 2.0);
    sum += functionValue;
}
Console.WriteLine("Estimate: {0,8:F4}", sum / length);

Dim sum = 0.0
Dim length = 10000
For Each v In QuasiRandom.FaureSequence(8, length)
    ' Evaluate the integrand
    Dim functionValue = 1.0
    For i = 0 To v.Length - 1
        functionValue *= Math.Abs(4.0 * v(i) - 2.0)
    Next
    sum += functionValue
Next
Console.WriteLine("Estimate: { 0,8:F4 }", sum / length)

No code example is currently available or this language may not be supported.

let mutable sum = 0.0
let length = 10000
for v in QuasiRandom.FaureSequence(8, length) do
    // Evaluate the integrand:
    let mutable functionValue = 1.0
    for i in 0 .. (v.Length - 1) do
        functionValue <- functionValue * abs(4.0 * v.[i] - 2.0)
    sum <- sum + functionValue
printfn "Estimate: %8.4f" (sum / float length)
//

The SobolSequence method returns a Sobol sequence. The method takes two or three arguments. The first is the dimension of the sequence. This is the number of elements of the vectors in the sequence. The maximum value for the dimension is 1111. The second argument is the length of the sequence. The method returns an IEnumerable<T> of Vector<T> of real numbers. The sequence always returns the same vector object, which is a mutable read-only vector. The elements of the vector change in each iteration, but they cannot be modified.

In the example below, we use a Sobol sequence approximate the same 8-dimensional integral as above:

double sum = 0.0;
int length = 10000;
foreach (var v in QuasiRandom.SobolSequence(8, length))
{
    // Evaluate the integrand:
    double functionValue = 1.0;
    for (int i = 0; i < v.Length; i++)
        functionValue *= Math.Abs(4.0 * v[i] - 2.0);
    sum += functionValue;
}
Console.WriteLine("Estimate: {0,8:F4}", sum / length);

Dim sum = 0.0
Dim length = 10000
For Each v In QuasiRandom.SobolSequence(8, length)
    ' Evaluate the integrand
    Dim functionValue = 1.0
    For i = 0 To v.Length - 1
        functionValue *= Math.Abs(4.0 * v(i) - 2.0)
    Next
    sum += functionValue
Next
Console.WriteLine("Estimate: { 0,8:F4 }", sum / length)

No code example is currently available or this language may not be supported.

let mutable sum = 0.0
let length = 10000
for v in QuasiRandom.SobolSequence(8, length) do
    // Evaluate the integrand:
    let mutable functionValue = 1.0
    for i in 0 .. (v.Length - 1) do
        functionValue <- functionValue * abs(4.0 * v.[i] - 2.0)
    sum <- sum + functionValue
printfn "Estimate: %8.4f" (sum / float length)
//

It is very cheap to skip elements at the start of the sequence. This makes it simple to run calculations on vectors in a Sobol sequence in parallel. There is some startup cost when generating a Sobol sequence, both in computation and memory. The SobolSequenceGenerator class allows you to generate multiple sequences of the same dimension while incurring the startup cost only once. This makes parallel computations more efficient. The class has one constructor that takes two arguments: the first is the dimension of the sequence. The second is the total number of elements that is requested.

The code below illustrates a typical parallel computation with Sobol sequences. It computes the same integral as above, but spreads the calculation over all available cores. The method starts by constructing a suitable SobolSequenceGenerator object. It then uses the Partitioner class to partition the total range into tuples of start and end points and runs the calculation in parallel over the partitions:

object lockObject = new object();
int length = 10000;
double sum = 0.0;
var sobol = new SobolSequenceGenerator(8, length);
Parallel.ForEach(Partitioner.Create(0, length), r =>
{
    double localSum = 0.0;
    foreach (var v in sobol.Generate(r.Item2 - r.Item1, r.Item1))
    {
        // Evaluate the integrand:
        double functionValue = 1.0;
        for (int i = 0; i < v.Length; i++)
            functionValue *= Math.Abs(4.0 * v[i] - 2.0);
        localSum += functionValue;
    }
    lock (lockObject) { sum += localSum; }
});
Console.WriteLine("Estimate: {0,8:F4}", sum / length);

Dim lockObject = New Object()
Dim length = 10000
Dim sum = 0.0
Dim sobol = New SobolSequenceGenerator(8, length)
Parallel.ForEach(Partitioner.Create(0, length),
    Sub(r)
        Dim localSum = 0.0
        For Each v In sobol.Generate(r.Item2 - r.Item1, r.Item1)
            ' Evaluate the integrand
            Dim functionValue = 1.0
            For i = 0 To v.Length - 1
                functionValue *= Math.Abs(4.0 * v(i) - 2.0)
            Next
            localSum += functionValue
        Next
        SyncLock lockObject
            sum += localSum
        End SyncLock
    End Sub)
Console.WriteLine("Estimate: { 0,8:F4 }", sum / length)

No code example is currently available or this language may not be supported.

let lockObject = obj()
let length = 10000
let mutable sum = 0.0
let sobol = SobolSequenceGenerator(8, length)
let sumPartition (r : int*int) =
    let mutable localSum = 0.0
    for v in QuasiRandom.SobolSequence(8, length) do
        // Evaluate the integrand:
        let mutable functionValue = 1.0
        for i in 0 .. (v.Length - 1) do
            functionValue <- functionValue * abs(4.0 * v.[i] - 2.0)
        localSum <- localSum + functionValue
    lock lockObject (fun () -> 
        sum <- sum + localSum)
Parallel.ForEach(Partitioner.Create(0, length), sumPartition) |> ignore

printfn "Estimate: %8.4f" (sum / float length)
//

Henri Fauré, Discrepance de suites associees a un systeme de numeration (en dimension s), Acta Arithmetica, Volume XLI, 1982, pages 337-351.

J.H. Halton, On the efficiency of certain quasi-random sequences of points in evaluating multi-dimensional integrals , Numerische Mathematik, Volume 2, 1960, pages 84-90.

J.H. Halton and G.B. Smith, Algorithm 247: Radical-Inverse Quasi-Random Point Sequence, Communications of the ACM, Volume 7, 1964, pages 701-702.