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 IEnumerableT
of VectorT
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))
{
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)
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
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 IEnumerableT
of VectorT
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))
{
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)
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
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 IEnumerableT
of VectorT
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))
{
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)
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
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))
{
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)
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
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.