Suggestion for sobol generator

[Richard Gould-6]

Hi

I've been looking at using the Sobol number generator class in QuantLib for running parallel monte-carlo simulations across many CPUs. I've modified the constructor to take an extra argument to the constructor to skip ahead in the sequence and thought you might want to consider adding this an extension to the interface. The code I used is as follows and will scale as O[log_2(N)] where N is the number of events to skip at the start of the sequence:

        // initialize the Sobol integer/double vectors
        // first draw
        unsigned long N = skip+1;
        unsigned int ops = ((unsigned int)(log(N)/M_LN2))+1;

        // Convert to Gray code
        unsigned long G = N ^ (N>>1);
        for (k=0; k<dimensionality_; k++)
        {
            integerSequence_[k] = 0;
            for(index=0; index<ops; index++)
            {
                if(G>>index & 1)
                    integerSequence_[k] ^= directionIntegers_[k][index];
            }
        }

        sequenceCounter_ = skip;

Regards
Richard Gould

NOTICE
The information contained in this email is confidential. If you are not the intended recipient, you must not disclose or use the information in this email in any way. If you received it in error, please tell us immediately by return email and delete the document. We do not guarantee the integrity of any e-mails or attached files and are not responsible for any changes made to them by any other person.

[Luigi Ballabio]

On 2/6/06, Richard Gould <[hidden email]> wrote:
> I've been looking at using the Sobol number generator class in QuantLib for running parallel monte-carlo simulations across many CPUs. I've modified the constructor to take an extra argument to the constructor to skip ahead in the sequence and thought you might want to consider adding this an extension to the interface.

Richard,
    apologies for the delay---it's been a busy month, besides one
marked by flu. I finally got around to look at your code. I'll check
it in the next few days and add it to the repository (unfortunately,
it might not get into the upcoming 0.3.12 release, as the latter is
being finalized.)

Thanks,
    Luigi

[Luigi Ballabio]

On 2/6/06, Richard Gould <[hidden email]> wrote:
> I've been looking at using the Sobol number generator class in QuantLib for
> running parallel monte-carlo simulations across many CPUs. I've modified the
> constructor to take an extra argument to the constructor to skip ahead in
> the sequence.

Richard,
    I finally got around to merge your code into the repository. It
compiles cleanly, but I had problems testing it. I expected that if I
skipped N samples and extracted 10 more, such 10 samples should be
equal to those obtained by initializing a generator without skipping,
extracting and discarding N samples, and extracting 10 more. This
should also hold in the degenerate case when N=0. However, this didn't
seem to be the case; the samples I obtained were somewhat different.
Can you check that your implementation passes such test and send me a
patch if it does?

Thanks,
    Luigi

[Richard Gould-6]

Hi Luigi

I wrote the attached harness to perform exactly this test (should have sent it with my code). I've tested skipping the first 1000 in the sequence up to dimension 50 with no problems.

Regards
Richard


-----Original Message-----
From: Luigi Ballabio [mailto:[hidden email]]
Sent: Thursday, 4 May 2006 3:51 PM
To: Richard Gould
Cc: [hidden email]
Subject: Re: [Quantlib-dev] Suggestion for sobol generator


On 2/6/06, Richard Gould <[hidden email]> wrote:
> I've been looking at using the Sobol number generator class in QuantLib for
> running parallel monte-carlo simulations across many CPUs. I've modified the
> constructor to take an extra argument to the constructor to skip ahead in
> the sequence.

Richard,
    I finally got around to merge your code into the repository. It
compiles cleanly, but I had problems testing it. I expected that if I
skipped N samples and extracted 10 more, such 10 samples should be
equal to those obtained by initializing a generator without skipping,
extracting and discarding N samples, and extracting 10 more. This
should also hold in the degenerate case when N=0. However, this didn't
seem to be the case; the samples I obtained were somewhat different.
Can you check that your implementation passes such test and send me a
patch if it does?

Thanks,
    Luigi


NOTICE
The information contained in this email is confidential. If you are not the intended recipient, you must not disclose or use the information in this email in any way. If you received it in error, please tell us immediately by return email and delete the document.  We do not guarantee the integrity of any e-mails or attached files and are not responsible for any changes made to them by any other person.