Does anyone success in initialising Sobol sequence

[mmillet]

I have tried to initialize Sobol sequence with my own polynomial and direction number as it is described in the vslNote.

First of all there is a lack of infornation in the documentation as we do not know the convention used to enter polynomial and direction( it is said that they must be 32 bit array, indeed there are many possibility: for polynomial the first coeficient is alwys one so many implementation of sobol take just the other coeficient...)

Anyway I didn't succes just to use my initialisation for the first 40 dimention, whatever I use as initial direction I got the same sequence with the following code.

paramToFill[0] = fGeneratorSize;

paramToFill[1] = VSL_USER_QRNG_INITIAL_VALUES;

//here I use a mix of the VSL note and the header definition. Anyway none of the other possibilityies work

paramToFill[2] = VSL_USER_PRIMITIVE_POLYMS | VSL_USER_INIT_DIRECTION_NUMBERS | VSL_USER_DIRECTION_NUMBERS;

if(overRideFirstDim) paramToFill[2] = paramToFill[2] | VSL_QRNG_OVERRIDE_1ST_DIM_INIT;

paramToFill[3+fGeneratorSize] = maxDegree;

for (size_t compteur = 0; compteur<fGeneratorSize; ++compteur)

{

SSSobolDimension& refSobolDim = sobolDimension[compteur];

size_t upperToFill = refSobolDim.degree;

if(tryComplete)

{

upperToFill = maxDegree;

refSobolDim.CompleteInitialisation(maxDegree);

}

paramToFill[3+compteur] = refSobolDim.polynom;

for (size_t idIndex = 0; idIndex<upperToFill; ++idIndex)

{

paramToFill[3+maxDegree*compteur+idIndex] = refSobolDim.initialDirection[idIndex];

}

}

VSLStreamStatePtr * streamGenerator = (VSLStreamStatePtr *) fGeneratorVector;

int status = 0;

for (size_t j=0; j<fThreadCount; ++j )

status = vslNewStreamEx( &streamGenerator[j], VSL_BRNG_SOBOL, fGeneratorSize,paramToFill));

Did I miss something?

[mmillet]

I found some misunderstanding in my code:

now I am using

vslNewStreamEx

( &streamGenerator[j], VSL_BRNG_SOBOL, parameterSize, paramToFill);

where parameterSize is the size of paramToFill.

Anyways if I am using less than or exactly 40dimension my initialisation is not taken into account and the sequences generated are the 40 dimension sequence of Bradley&Fox. If I try 41 or more dimention then myt initilisation is not taken into account but now only one dimension is generated and when I am calling

vdRngUniform

( VSL_METHOD_SUNIFORM_STD_ACCURATE, streamGenerator[threadID-1], fGeneratorSize, fLocalBuffer.GetReference(), 0.0, 1.0 );

I get the sequence generated by using a congruence modulo fGeneratorSize on the first sequence which is very very bad in case of sobol...

[Andrey_N_Intel]

Thank you for your post and feedback to the VSL documentation.

Ideally, in order to understand why things go wrong during Sobol initialization on your side it could be nice to have working short test case which would demonstrate the issue.

To make initialization of VSL Sobol QRNG you need to pass into the library either polynomials and initial direction numbers or table of direction numbers computed on you side. Initialization of the Sobol QRNG is done using vslNewStreamEx function:

errcode = vslNewStreamEx( &stream, brng, nparams, params );

Prior call to the function nparams and params should be properly initialized. Array of unsigned integers params of size nparams is used to pack Sobol initialization values.
First element is maximal dimension of quasi-random vector you are planning to have:

params[0] = dimen;

The next field defines that QRNG parameters are passed into the library. It should be set to VSL_USER_QRNG_INITIAL_VALUES:

params[1] = VSL_USER_QRNG_INITIAL_VALUES;

Value of the next field determines what you are going to pass in the library, polynomials/initial direction numbers or table of Sobol direction numbers.
Let's set it to VSL_USER_PRIMITIVE_POLYMS | VSL_USER_INIT_DIRECTION_NUMBERS if there is a need to pass both polynomials & initial direction numbers:

params[2] = VSL_USER_INIT_DIRECTION_NUMBERS | VSL_USER_PRIMITIVE_POLYMS;

The next dimen positions are used to pack polynomials:
for ( i = 0; i < dimen; i++ ) params[3+i] = uSobolPoly;

The polynomials are packed in the following way: least significant bit contains coefficient at zero degree, the next bit contains coefficient at 1st degree and so forth.

The library needs to know about the maximal dimension of polynomials (over all polynomials used to initialize Sobol):

params[3+dimen] = maxdeg;

The rest positions are used to hold initial direction numbers:
k = 4+dim;
for ( i = 0; i < dimen-1; i++ )
{
for ( j = 0; j < maxdeg; j++ )
{
params[k++] = uSobolMInit;
}
}

where 2 dimensional table uSobolMInit of size num_of_poly*maxdeg contains initial direction numbers. If you do not want to override default initialization for 1st dimension num_of_poly equals to dimen-1 (otherwise, it equals to dimen and params[2] = params[2] | VSL_QRNG_OVERRIDE_1ST_DIM_INIT).

Thus, size nparams of array params should be enough to hold at least 3 + num_of_poly + 1 + num_of_poly * maxdeg:

nparams = 3 + num_of_poly + 1 + num_of_poly * maxdeg;

Value of parameter nparams depends on type of initialization and, thus, number of Sobol parameters to be passed in VSL.

Finally, vslNewStreamEx routine is called to initialize Sobol QRNG.

Please, let me know if you have more questions/concerns on this and, also, how this scheme works for you. Also, please fill free to attach a test case with your initialization of VSL Sobol.

Also, what version of MKL do you use?

Andrey

[mmillet]

Thank for your explanation.

As far as the convention for polynimial, ifI understand correctly x^3+x+1 is 1011 = 1*8+0*4+1*2+1=11.

For the initial direction it is not clear. Let us says that my first two dimension are 0.5 and 0.75 it can be coded in an int as 1 and 3 or (1 << 31) and (2 << 30). What is the convention?

If I understand correctly the flag VSL_QRNG_OVERRIDE_1ST_DIM_INIT is to say if I overwrite the first dimension. But in any case I can fill dimen polynomial and dimen*maxdeg initial direction. If I do not overwrite I just give 1 dimension useless data. I am right?

I am using MKL 10.0.

Thanks for you help.

[mmillet]

I finally manage to initialize sequence:

the initial direction convention is 0.5 and 0.75must be coded as 1 and 3.

when using VSL_QRNG_OVERRIDE_1ST_DIM_INIT the initialisation data are exactly as explained in the vsl note. If you do not use then for a sequence of dimention dimen you must use dimen in parameter 0 but for the positions in vector param you must replace dimen by dimen-1 everywhere.

[Andrey_N_Intel]

As far as the convention for polynimial, if I understand correctly x^3+x+1 is 1011 = 1*8+0*4+1*2+1=11.
Yes, this is correct.

For the initial direction it is not clear. Let us says that my first two dimension are 0.5 and 0.75 it can be coded in an int as 1 and 3 or (1 << 31) and (2 << 30). What is the convention?
First variant, that is 1 and 3, should be used for initial direction numbers.

If I understand correctly the flag VSL_QRNG_OVERRIDE_1ST_DIM_INIT is to say if I overwrite the first dimension.
Yes, that's right.

But in any case I can fill dimen polynomial and dimen*maxdeg initial direction. If I do not overwrite I just give 1 dimension useless data. I am right?
Let us consider the example. You want to pass initial parameters for dimen=40 and do not want to override 1st dimension.In this case you need to pass only 40-1=39 sets of polynomials/initial direction numbers.If you want to override 1st dimension you need to pass 40 sets of polynomials/initial direction numbers.Value of maxdeg should be properly defined in each case. In other words, the library expects only "useful" data. Passing 1st dimension related data into the library when flag VSL_QRNG_OVERRIDE_1ST_DIM_INIT is not set would result in incorrect generation of quasi-random sequences.

Please, let me know how generation of quasi-random vectors with your initialization parameters work. Please, feel free to address any questions/concerns on this. Also, you may want to have a look into soboluserdirnums.c file which is available in examplesvslcsource directory of your MKL installation package.

Thanks,
Andrey

[mmillet]

Thanks for the explanation. I will look in the .c files.

In fact my initialisation work only when I am asking for less than 40 dimensions simultaneously. If I am asking for let says 41 dimension then all dimension return the first dimension sobol sequences. I do not know if it is an internal limitation or if it is a remaining mistake in my code.

[Andrey_N_Intel]

Can you please provide short test case which would demonstrate the issue and/or output of the generator with parameters of the generation? Thanks, Andrey

[mmillet]

As it is difficult to provide my original code (cpp with many object and method) I have modified the code of the example provided by the MKL: I have added one polynom in the list and the corresponding dimention. Then I initialize a 40 dimensions generator and if you look in the array named r after the initialisation you will see that all the 40 sequences are differents. Then I do the same initialisation with dimension 41 and now if you look in array r you will see that only the first and second dimension are differents, all the others are the same as the first one.

Did I make a mistake somewhere?

//the code:

/* Dimension of quasi-random vectors */

#define

DIM 41

/* Number of bits used for storing unsigned integer */

#define

NBITS 32

/* Number of primitive polynomials */

#define

NPOLY 40

/* Maximum degree of primitive polynomials */

#define

MAXDEG 8

/* Number of digits to scramble */

#define

S_DIGITS 32

/* Table of initial direction numbers */

static

unsigned int uSobolMInit[NPOLY][MAXDEG] = {

1,0,0, 0, 0, 0, 0, 0,

1,1,0, 0, 0, 0, 0, 0,

1,3,7, 0, 0, 0, 0, 0,

1,1,5, 0, 0, 0, 0, 0,

1,3,1, 1, 0, 0, 0, 0,

1,1,3, 7, 0, 0, 0, 0,

1,3,3, 9, 9, 0, 0, 0,

1,3,7,13, 3, 0, 0, 0,

1,1,5,11,27, 0, 0, 0,

1,3,5, 1,15, 0, 0, 0,

1,1,7, 3,29, 0, 0, 0,

1,3,7, 7,21, 0, 0, 0,

1,1,1, 9,23,37, 0, 0,

1,3,3, 5,19,33, 0, 0,

1,1,3,13,11, 7, 0, 0,

1,1,7,13,25, 5, 0, 0,

1,3,5,11, 7,11, 0, 0,

1,1,1, 3,13,39, 0, 0,

1,3,1,15,17,63, 13, 0,

1,1,5, 5, 1,27, 33, 0,

1,3,3, 3,25,17,115, 0,

1,1,3,15,29,15, 41, 0,

1,3,1, 7, 3,23, 79, 0,

1,3,7, 9,31,29, 17, 0,

1,1,5,13,11, 3, 29, 0,

1,3,1, 9, 5,21,119, 0,

1,1,3, 1,23,13, 75, 0,

1,3,3,11,27,31, 73, 0,

1,1,7, 7,19,25,105, 0,

1,3,5, 5,21, 9, 7, 0,

1,1,1,15, 5,49, 59, 0,

1,1,1, 1, 1,33, 65, 0,

1,3,5,15,17,19, 21, 0,

1,1,7,11,13,29, 3, 0,

1,3,7, 5, 7,11,113, 0,

1,1,5, 3,15,19, 61, 0,

1,3,1, 1, 9,27, 89, 7,

1,1,3, 7,31,15, 45,23,

1,3,3, 9, 9,25,107,39,

1,3,1,11,11,11,77,249

};

/* Array of primitive polynom ials */

static

unsigned int uSobolPoly[NPOLY] = {

3, 7, 11, 13, 19, 25, 37, 59, 47,

61, 55, 41, 67, 97, 91, 109, 103, 115, 131,

193, 137, 145, 143, 241, 157, 185, 167, 229, 171,

213, 191, 253, 203, 211, 239, 247, 285, 369, 299,

301

};

/* Array of degrees of primitive polynomials */

static

int iSobolPolyDeg[NPOLY] = {

1, 2, 3, 3, 4, 4, 5, 5, 5,

5, 5, 5, 6, 6, 6, 6, 6, 6, 7,

7, 7, 7, 7, 7, 7, 7, 7, 7, 7,

7, 7, 7, 7, 7, 7, 7, 8, 8, 8,

8

};

/* Generator matrix, table of direction numbers */

static

unsigned int uSobolV[DIM][NBITS];

/* Generator matrix, table of modified direction numbers */

static

unsigned int uSobolSV[DIM][NBITS];

/* Routine to initialize matrix of Sobol */

static

void InitSobolMatrix( int dimen );

/* Routine to initialize Sobol matrix for Owen-like scrambled sequence */

static

void InitSobolMatrixOs( int dimen );

/* Routine to generate sequence of low triangular matrices for

Owen scrambling */

static

void GenerateLSM( int dimen, unsigned int lsm[ DIM][S_DIGITS] );

/* Dimension of vectors to generate */

/* Size of params array */

#define

NU DIM*NBITS + 4

#define

MaxDim 41

#define

NN 10

/* Size of buffer for quasi-random sequence */

#define

N 1000*MaxDim

/* Buffer for quasi-random sequence */

double

r[N];

double

rt[N];

int

testGeneration()//main()

{

VSLStreamStatePtr stream; MKL_INT brng; int dim; int errcode; double a, b; int i, j, k; MKL_INT nparams; int errnum; /***** Initialize *****/ brng = VSL_BRNG_SOBOL; a = 0.0; b = 1.0; /* dim=40 work fine */ dim = 40; size_t sizeToUse = dim-1;// must be dim if we overwrite 1st dim unsigned int maxDegree = iSobolPolyDeg[sizeToUse-1]; unsigned int parameterSize = 4+sizeToUse*(1+maxDegree); unsigned int *params = new unsigned int[parameterSize]; params[0] = dim; params[1] = VSL_USER_QRNG_INITIAL_VALUES; params[2] = VSL_USER_PRIMITIVE_POLYMS | VSL_USER_INIT_DIRECTION_NUMBERS; params[3+sizeToUse] = maxDegree; k = 3; for ( i = 0; i < sizeToUse; i++ )

{

params[3+i] = uSobolPoly[i]; for ( j = 0; j < maxDegree; j++ ) params[4+sizeToUse+maxDegree*i+j] = uSobolMInit[i][j< /FONT>];

}

/***** Initialize ****/ errcode = vslNewStreamEx( &stream, brng, parameterSize, params ); //CheckVslError( errcode ); /***** Call RNG *****/ errcode = vdRngUniform( VSL_METHOD_SUNIFORM_STD_ACCURATE, stream, 1000*dim, r, a, b ); //CheckVslError( errcode ); delete [] params; /***** Printing results *****/ printf("Sobol sequence with user-defined direction numbers: "); printf("-------------------------------------------------- "); printf("Parameters: "); printf(" a=%.4f ",a); printf(" b=%.4f ",b); k = 12; printf("Results: %d component of quasi-random sequence (first %d of %d): ", k, NN, N/dim); for(i=0;i<NN;i++) { printf("%.4f ",r[k+i*dim]);

}

printf(" "); /***** Deinitialize *****/ errcode = vslDeleteStream( &stream ); //CheckVslError( errcode ); /* dim=41 does not work */ dim = 41; sizeToUse = dim-1;// must be dim if we overwrite 1st dim maxDegree = iSobolPolyDeg[sizeToUse-1]; parameterSize = 4+sizeToUse*(1+maxDegree); params = new unsigned int[parameterSize]; params[0] = dim; params[1] = VSL_USER_QRNG_INITIAL_VALUES; params[2] = VSL_USER_PRIMITIVE_POLYMS | VSL_USER_INIT_DIRECTION_NUMBERS; params[3+sizeToUse] = maxDegree; k = 3; for ( i = 0; i < sizeToUse; i++ )

{

params[3+i] = uSobolPoly[i]; for ( j = 0; j < maxDegree; j++ ) params[4+sizeToUse+maxDegree*i+j] = uSobolMInit[i][j];

}

/***** Initialize ****/ errcode = vslNewStreamEx( &stream, brng, parameterSize, params ); //CheckVslError( errcode ); /***** Call RNG *****/ errcode = vdRngUniform( VSL_METHOD_SUNIFORM_STD_ACCURATE, stream, 1000*dim, r, a, b ); //CheckVslError( errcode ); delete [] params; /***** Printing results *****/ printf("Sobol sequence with user-defined direction numbers: "); printf("-------------------------------------------------- "); printf("Parameters: "); printf(" a=%.4f ",a); printf(" b=%.4f ",b); k = 12; printf("Results: %d component of quasi-random sequence (first %d of %d): ", k, NN, N/dim); for(i=0; i<NN;i++) { printf("%.4f ",r[k+i*dim]);

}

printf(" "); /***** Deinitialize *****/ errcode = vslDeleteStream( &stream ); //CheckVslError( errcode ); return 0;

}

[mmillet]

Just a remark. I am using the debugger to look in the array. If you look the result (the 12th dimension) then you must see that it is not the same when I am usig 40 and when I am using 41.

[Andrey_N_Intel]

Thank you for the test case. It appears there isan issue in the QRNG initialization procedure. We'll investigate and fix it. Meanwhile, feel free to report about the issue to Premier Support. As a workaround you may want to compute table of the Sobol direction numbers using your polynomials/initial direction numbers and passit into the generator as demonstrated in the example soboluserdirnums.c(which is availablein MKL installation package).Thanks, Andrey

[mmillet]

OK I will report the problem tomorrow. I have checked that the workarround allow to define more than 40 simultaneous dimensions.

Thanks for your help.

[s_kucherenko]

"I have tried to initialize Sobol sequence .."

Initialization of the Sobol' sequence or developingthe direction numbers is not easy and may take many months-men work. There are known well tested high-dimensional Sobol' sequence generators.

Currently, the Sobol' sequence generator with highest maximum dimension isSobolSeq6144. It wasdeveloped byBRODA ( www.broda.co.uk ) in collaboration with I. M. Sobol'. Not only this generator has very high dimensionality and employs the super fast generation algorithm but generated Sobol' sequences satisfy homogeneity Property A in all dimensions.