Re: [boost] Review Request: random number library

From: troyer (troyer_at_[hidden])
Date: 2000-05-29 12:11:49

I have now looked into the library in more detail with our 
local random number expert. I admire the effort that went into
creating it and would like to see it improved so that we can use
it in our production codes soon. As we have a large number of
comments I will order them by topics:

1.) Floating point generators
    a) problems with min() and max()
    b) validation
2.) Distributions
    a) problems with the distributions
    b) uniform_sphere is incorrect
    b) more distributions
3.) Performance issues
    a) uniform_smallint
    b) uniform_sphere 
    c) floating point generators
4.) Quality
    a) uniform_int
    b) mersenne_twister
    c) tests of the generators
5.) Other useful additions

Despite the length I hope that interested people will find the time 
to read this mail and to reply.
 

1.) Floating point generators

Most random number gernators that we use are floating point
generators as they can be much faster than integer generators.
Floating point generators are also supprted by the library but
I see the following problems:

1.a.) the concept UniformRandomNumberGenerator specifies

 X::min_value
 X::max_value
 v.min()
 v.max()

to be strict lower respective upper bounds. This is perfect
for integer generators but is problematic for floating point
generators that usually produce numbers in an interval
[0,1), i.e. 0 <= x < 1. This is a subtle issue, but in some
cases it is absolutey essential that the number 1 is NEVER
generated. The simples example is provided by the 
exponential_distribution class template in file random.hpp:

  result_type operator()() { 
    using std::log;
    return -1.0 / _lambda * log(1-_rng());
  }

Here log(1-_rng()) is used to make sure that the argument is
never zero and the log always defined!

What are the possible solutions? I see two options:
  
i)  add more members to the UniformRandomNumberGenerator 
    concept such as
  
     static const bool can_assume_lower_bound = /* */;
     static const bool can_assume_upper_bound = /* */;

    This allows all possible cases for the interval: 
    (min,max), [min,max), (min,max] and [min,max]

ii) always assume the interval to be [min,max) for floating point 
    generators and [min,max] for integer generators.

iii) use min()<= x < max() for floating point and integer.
     If an integer generator can assume the largest possible 
     integer value we cannot define max() and this option
     is useless.

The second solution has the disadvantage that the semantics
of max() is different for floating point and integer generators.
I still prefer it since in all my applications I need floating 
point random numbers to be in an interval a <= x < b while I 
never need them to ever assume the upper bound. If howver
anybody knows of an application where having random numbers
in an interval [a,b] is essential we have to go for a solution
such as i).

These is a related inconsisteny with the classes uniform_real 
and uniform_01. The former produces values in [min,max], the
latter in [0,1). If one goes for solution ii) the uniform_real
class should be changed to produce numbers in [min,max) and 
uniform_01 would just be a specialization of uniform_real.

1.b.) validation of floating point generators.

The PseudoRandomNumberGenerator concept defines a function 

T validation() 

that returns a precomputed 10001th element to be used for validation. 
While this is great for integer based RNGs it fails for floating 
point generators. Even assuming IEEE arithmetic as far as I know there 
is some leeway in how rounding is done and the 10001th elemt can differ
slightly from one machine to another. As solutions I propose a 
modified validation function

bool validation(T x)

that checks if the calculated number x is within the expected numerical 
errors of a precalculated hardcoded one. Other solutions, such as
providing another member function returning the accepted numerical
error are also feasible. I however find above method the most elegant
one.

2.) Distributions

2.a) Problems with the distributions

The library already provides a number of distributions, like:

bernoulli_distribution
geometric_distribution
exponential_distribution
normal_distribution
lognormal_distribution
uniform_sphere

I have not looked at lognormal_distribution and 
geometric_distribution but have checked the others

The documentation for all these distributions reads:
It transforms a uniform distribution into a ....

This is not completely correct. For the code to work
it must be a distribution in [0,1) (compare 1.a), such
as provided by the "uniform_01" generator. Without much
effort one could rewrite them to use an interval [a,b),
but as most floating point generators provide output in [0,1)
I don't think this to be a problem. However we definitively
need at least runtime checks on the range, such as, instead of:

  exponential_distribution(base_type& rng, result_type lambda)
    : _rng(rng), _lambda(lambda) { assert(lambda > 0); iterator_init(); }

we need code like:

  exponential_distribution(base_type& rng, result_type lambda)
    : _rng(rng), _lambda(lambda) 
  { assert(lambda > 0 && rng.min()==0. &&
rng.max()==1.);iterator_init(); }

However before starting any coding the conventions for floating
point generators (see 1.a) have to be defined.
             
2.b) uniform_sphere

Our local random number expert has written a paper on this topic,
and we have used his methods extensively. Thus I realized that
the implementation of uniform_sphere is incorrect as it is.
The incorrectness is simple to see: it will never produce negative 
numbers. Furthermore a cube is projected onto a sphere the
distribution is not uniform. It IS correct if the input random number
generator produces Gaussian random numbers. Either the documentation
should be changed or the normal_distribution class used inside
uniform_sphere. Further discussion of this class under 3.)

2.c) more distributions

I propose to replace uniform_sphere by two classes, uniform_on_sphere
for points on the surface on a sphere and uniform_in_sphere for points
anywhere inside a sphere. We can contribute these classes as well as
Poisson, beta and exponential distributions in a finite range some time
in the next weeks.

3.) Performance issues

Even using fast generators we often spend 10%-20% of our CPU time on
generating random numbers. Using a slow generator like the Mersenne 
twister this grows to more than 95%. Thus performance is a crucial
factor.

3.a.) uniform_smallint

uses the following function to convert to a smaller interval:

  result_type operator()()
  {
    // we must not use the low bits here, because LCGs get very bad then
    return ((_rng() - _rng.min()) / _factor) % _range + _min;
  }

This is very costly if _factor is not a power of two. My proposal
is to change it to change the division by a bit shift >>
and to use a rejectance step to avoid quantization errors.

3.b.) uniform_sphere

In my own uniform_sphere generator I had to struggle with
this generator since in contrast to the others it does
NOT return a single number but one number per dimension,
such as a std::vector object.
An implementation like

template<class UniformRandomNumberGenerator, class RealType>
std::vector<RealType>
uniform_sphere<UniformRandomNumberGenerator, RealType>::operator()()
{
  result_type result(_dim);
  //.... 
  return result;
}

might make sense but is extremely inefficient. In my applications
I often need points on the surface of a 3-sphere. The overhead of
allocating a full std::vector object is immense compared to the
small computational effort. I came up with several solutions, 
one of which I used in my code:

  i) pass a std::vector<RealType> by refernce and have the function
     copy into it:

template<class UniformRandomNumberGenerator, class RealType>
void
uniform_sphere<UniformRandomNumberGenerator, RealType>::operator()
(std::vector<RealType>& result)
{
  result.resize(_dim);
  //.... 
}

   ii) pass a container using an iterator pair:
template<class UniformRandomNumberGenerator, class RealType>
template<class Iterator>
void
uniform_sphere<UniformRandomNumberGenerator, RealType>::operator()
(Iterator begin, Iterator end)
{
   //.... 
}

   iii) this is the solution I prefer for large dimensions. 
        Provide an extra data member of type std::vectory<RealType>,
        have the operator() write into that data member and
        have it return a const reference to this data member.
        As far as I see at the moment this should allow all possible
        uses. of course one need not return a std::vector but could
        provide the output container as an additional template
        parameter.

Another big performance issue are low dimensions, e.g. 3 dimensions.
There instead of a generator based on Gaussianrandom numbers faster 
uniform random numbers combined with an acceptance/rejactance step 
can be used. Also in this case I prefer the loops in 

  for(int i = 0; i < _dim; i++) {
    RealType val = _rng();
    result[i] = val;
    sqsum += val * val;
  }
  using std::sqrt;
  // for all i: result[i] /= sqrt(sqsum)
  std::transform(result.begin(), result.end(), result.begin(),
                 std::bind2nd(std::divides<RealType>(), sqrt(sqsum)));

to be unrolled explicitly. I did that using Todd's marvelous 
TinyVector classes from the Blitz++ library. Although there are
no TinyVectors in boost a specially optimized implementation
can still be done explicitly using template meta programs. As 
this is crucial to my research work I am willing to contribute
an optimized uniform_sphere for low dimension where the dimensionality
is known at compile time.

3.c.) floating point generators

Floating point distributions are created from integer ones by
costly floating point divisions. A faster approach is to use
generators that work on floating point numbers, such as
lagged Fibonacci generators. We can contribute such floating
point generators to the library and can do so probably next week.

4.) Quality

4.a) uniform_int

This class provides random integers with more bits
than the base random number generatot by concatenating
several unform random integers. As in LCGs the low bits are
usually bad it is useless to concatenate it with another
random number. Instead of supporting this bad habit by a 
library in the class I would remove it and rather encourage
the use of generators with more bits. Anyways, at least a strong
warning should be written into the documentation. 

4.b) mersenne_twister

As these generators seem to be favorite ones at the moment I tried
using them in my codes and was surprised that the execution time
increased by a factor of 8! Before incuring such a penalty we wanted
to check the generator on our tests, such as random walk and Ising model
tests. Since the Merseene Twister failed the random walk test badly
we do not use it anymore.

4.c.) tests of the generators

We are currently running extensive tests with a number of generators,
including the Mersenne Twisters, lagged Fibonacci generators and the
ranlux generator of the CERN library. We now plan to also test all the
generators in this library and to publish the results once the tests 
are finished.

5.) Other useful additions

A disadvantage of the current design is that there is only one
construcor that works for all generators, namely the default
constructor. And there is no seed() function that is common
to all of them. Thus without knowledge of the specific generator
we can produce only one sequence. I suggest that at least a function

void seed(uint32_t s);

and the corresponding constructor should be provided. That way
one can use more than one random sequence.

We have another extension that is probably not suitable for boost, 
but is very useful if one wants to switch generators at run time. 
It is an abstract base class and a templated derived class using the
Barton-Nackmann trick to provide both an efficient inline operator()
as well as run time polymorphism. If there is interest I can convert
it easily to this random number library and contribute it.

I hope these comments by a non-computer scientist are useful
and look forward to a good and standardized random number generator 
library.

Matthias Troyer