Random Matrices

Random matrices in Skylark

Generating random samples in Skylark

We currently base our approach on Threefry4x64, a counter-based random number generator (cbrng) from [http://www.deshawresearch.com/resources_random123.html Random123] library]. This generator, given a counter and a key (input) produces a result (output) by internally going through a number of 20 rounds of transformations:

result = cbrng(counter, key)

Each of the counter, key and result is implemented as an array of four uint64_t numbers.
counter[i], key[i], result[i] (alternatively counter.v[i], key.v[i], result.v[i]) for i=0,1,2,3 give access to these unsigned integers. counter and key are not modified during the internal transformations to produce result; no state is preserved in Threefry4x64.

MicroURNG<Threefry4x64> is a C++0x Uniform Random Number generator (urng) that can be constructed by a counter and a key. There follow the details of urng mechanics, basically its () operator. The important point is that there is an internal counter n telling us how many times a counter-based random number generator has been fired and by convention this n is stored in the high 32 bits of counter[3]. Consequently we should pass a counter with clean its high-32 bits in counter[3]; failing to do this throws a runtime exception.

• The 1st time () is invoked: result = cbrng(counter, key) and result[3] is returned.
• The next 3 times result[2], result[1] and result[0] unint64_t integers are successively returned and no cbrng is called.
• The fifth time () is called a fresh result is generated having as counter the one passed in the constructor with everything intact except for its high 32 bits of counter[3]: the first time this 32-bit part was n=0, now it has changed to n=1. n is incremented next and result[3] (in this fresh result) is returned.
• The 6th, 7th, 8th time of () invocation result[2], result[1] and result[0] of the fresh result will be returned in succession (again no new cbrng calls)...
• ...and so on.

This description says that for a specific urng(counter, key) - with the "clean high 32 bits of counter[3]" restriction - we can call 2^32^ times a cbrng and get 4 uint64_t's each time; in other words we can get up to 4 x 2^32^ such random numbers from a single urng for a given counter and key. This is definitely a small number for our purpose of populating potentially huge matrices.

Instead we could internally represent lazily computed arrays of such urngs (urng_array). The urngs would all share the same key (which could be chosen to be Skylark's context seed). Also the counter passed to the constructor of each urng would have its counter[0] set to any number in [0, 2^64^-1] (with "clean" counter[1], counter[2] and - especially counter[3]). A urng_array would have a base and a size that communicates the fact that the counter[0] of its urngs can be a number in [base, base+size) although it is exposed as zero-based for the [] operator. '''We have chosen not to provide such urng_arrays but compute each urng as needed. However the urng_array metaphor helps in the discussion that follows.'''

Finally, we want to generate number numbers drawn from a a distribution and not urngs. Ideally we would like to have ''lazily'' computed arrays of such samples for filling our matrices that could also serve as part of the data-specific section of a sketch object.

A random_samples_array fills this role; it is templated on the types of the samples and of the distribution and could be thought of as backed by a urng_array:

random_samples_array[i] <- distribution(urng_array[i])

There is a subtle issue here: urng <- urng_array[i] is essentially instantatiated by

• counter[0] = base + i, counter[1] = 0, counter[2] = 0, counter[3] = 0, and
• key = seed and it can be the case that the distribution will internally call urng's () operator more than once.

This can be due to:

• underlying algorithm "mechanics", e.g. Box-Muller for generating normal distribution samples works by calling () two times at odd-numbered stream locations (while reusing them at even-numbered steps and not calling () then) see [http://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform link]
• implementation details such as rejection semantics to ensure correctness as described [http://stackoverflow.com/questions/14689914/uniformly-distributed-random-number-generation here] or as a way for dealing with boundary values, e.g. see for(;;){...} code snippet in implementing operator () for a key building-block distribution, boost::random::uniform_01<RealType>[http://svn.boost.org/svn/boost/trunk/boost/random/uniform_01.hpp link]. The number of invocations is not known beforehand here!)

So if it were the case that urng() simply did counter[0]++ then e.g. for the generation of 2 successive normal distribution numbers let's say samples[i] and samples[i+1] what would happen is:

• urng_array[i] (for samples[i]) would make use of counter = base+i and counter = base+i+1 in the cbrng's
• urng_array[i+1] (for samples[i+1]) would make use of counter = base+i+1 and counter = base+i+2 and this means the introduction of unacceptable ''dependencies'' in the stream.

However in our implementation (where for a moment we see the counter buffer as a uint32_t counter[8] object):

• urng_array[i] would make use of counter = (0,0,0,0,0,0,0,base+i) and counter = (1,0,0,0,0,0,0,base+i) in the cbrng's
• urng_array[i+1] would make use of counter = (0,0,0,0,0,0,0,base+i+1) and counter = (1,0,0,0,0,0,0,base+i+2) (also same key, i.e. the seed) so ''no dependencies'' are introduced.

It follows that the ''multiplicity'' of the internal invocations of urng's () operator by the distribution implementation is captured in counter[7] and can be up to 2^32 - it is only 2 in Box-Muller case! This also helps in ensuring a specific size for the random_samples_array at construction time: since the number of counter increments (by one) does not necessarily coincide with the size we request, if we adopted a purely linear view of our samples then when asking the context for let's say 2 arrays, we would not know for sure the start of the second array in the global stream before generating all the elements of the first array and keeping track of the accumulated "multiplicities".

Distributions used in Skylark

Here is a table of the distributions used in each sketch. The sketch name and the type of the distribution are given (all in namespace boost::random except for rademacher_distribution_t - in skylark::utility); ValueType denotes the type of the random matrix entries.

Sketch Name	Distribution Types
CT	cauchy_distribution<ValueType>
CWT	uniform_int_distribution<int>, rademacher_distribution_t<ValueType>
FJLT	uniform_int_distribution<int>, rademacher_distribution_t<ValueType>
JLT	normal_distribution<ValueType>
MMT	uniform_int_distribution<int>, cauchy_distribution<ValueType>
Gaussian RFT	uniform_real_distribution<double>, normal_distribution<ValueType>
Laplacian RFT	uniform_real_distribution<double>, cauchy_distribution<ValueType>
WZT	uniform_int_distribution<int>, exponential_distribution<ValueType>

So this is the table of the six distribution types we need:

Distribution Types
uniform_int_distribution<int>
uniform_real_distribution<double>
rademacher_distribution_t<ValueType>
cauchy_distribution<ValueType>
exponential_distribution<ValueType>
normal_distribution<ValueType>

We could implement them based on the raw four uint64_t numbers we get in result[] each time we call a counter-based generator (for given counter and key arguments):

result = cbrng(counter, key)

Assuming for a moment we don't get the boundary values min=0, max=max_uint64 in result[0], result[1] (in that case we could utilize some of the other result[i] "parts"), here is pseudocode for producing sample samples[i] for each distribution type (key=seed (fixed) and counter[] parts not mentioned assumed zeroed - can be further simplified since min=0).

uniform_int_distribution<int>(a, b)

counter[0] = base + i
result = cbrng(counter, key)
bucket_size = (max - min + 1) / (b - a + 1)
samples[i] = a + (result[0] - min) / bucket_size

uniform_real_distribution<double>(a, b)

counter[0] = base + i
result = cbrng(counter, key)
samples[i] = a + (b - a) * (result[0] - min) / (max - min) 

rademacher_distribution_t<ValueType>

counter[0] = base + i
result = cbrng(counter, key)
value = (result[0] - min) / (max - min)
if value > 0.5 then samples[i] = 1 else samples[i] = -1

cauchy_distribution<ValueType>(median, sigma)

counter[0] = base + i
result = cbrng(counter, key)
samples[i] = median + sigma * tan (pi*(result[0] - 0.5))

exponential_distribution<ValueType>(lambda)

counter[0] = base + i
result = cbrng(counter, key)
samples[i] = -1.0/lambda * log(1.0 - result[0])

normal_distribution<ValueType>(mean, sigma)

counter[0] = base + i
result = cbrng(counter, key)
samples[i] = mean + sigma * sqrt(-2*log(1.0-result[1])) * cos(2*pi*result[0])