Really slow beta distribution

[marcospassos]

Description

Computing the beta distribution is really slow when alpha >= 1.5 and beta <= 1. For example, computing the beta distribution for alpha = 1001, and beta = 1, takes up to 8 seconds.

Any suggestion on how to make it faster?

Reproduction

import iterator from '@stdlib/random-iter-beta';

iterator(1001, 1, {iter: 10000, seed: 1234});

Version

0.0.5

Environments

Node.js

Node.js / npm Version

v14.17.1

Platform

OSX

Checklist

• Read and understood the Code of Conduct.
• Searched for existing issues and pull requests.

[github-actions]

🎉 Welcome! 🎉

And thank you for opening your first issue! We will get back to you shortly. 🏃 💨

[kgryte]

@marcospassos Just did a quick test in the stdlib REPL and I am not able to reproduce the slowdown. One question: when you refer to 8 seconds are you doing something like the following in order to generate an array of random numbers:

In [138]: var it = random.iterators.beta(1001,1, {iter:10000,seed:1234});

In [139]: var arr = [...it]
Out[139]: [...]

or does it take 8 seconds to get a single value (e.g., it.next())?

If the former, while generating 10K random numbers on my 2013 Macbook Pro took 1-2 seconds, one possible speedup is to do

import gfillBy from '@stdlib/blas-ext-base-gfill-by';
import beta from '@stdlib/random-base-beta';

let x = new Float64Array(10000);
gfillBy(x.length, x, 1, beta.factory(1001,1,{'seed':1234}));

When I did this locally, I did not notice a significant perf difference TBH, but your mileage my vary.

[kgryte]

I also just ran the benchmarks for @stdlib/random/iter/beta (both in Node.js and in Chrome) and got the following results:

Running benchmark: /stdlib/stdlib/lib/node_modules/@stdlib/random/iter/beta/benchmark/benchmark.js
TAP version 13
# @stdlib/random/iter/beta
  ---
  iterations: 1000
  elapsed: 0.090648618
  rate: 11031.607784687903
  ...
ok 1 benchmark finished
# @stdlib/random/iter/beta
  ---
  iterations: 1000
  elapsed: 0.095458372
  rate: 10475.77052749234
  ...
ok 2 benchmark finished
# @stdlib/random/iter/beta
  ---
  iterations: 1000
  elapsed: 0.076719305
  rate: 13034.528923326403
  ...
ok 3 benchmark finished
# @stdlib/random/iter/beta::iteration
  ---
  iterations: 1000000
  elapsed: 0.442267812
  rate: 2261073.4330356377
  ...
ok 4 benchmark finished
# @stdlib/random/iter/beta::iteration
  ---
  iterations: 1000000
  elapsed: 0.440454849
  rate: 2270380.272280758
  ...
ok 5 benchmark finished
# @stdlib/random/iter/beta::iteration
  ---
  iterations: 1000000
  elapsed: 0.433772116
  rate: 2305357.9589703274
  ...
ok 6 benchmark finished
#
1..6
# total 6
# pass  6
#
# ok

Based on the above, the package was able to generate 1M random numbers in less than half a second.

So, I am not sure the cause for the slowdown you are seeing. 😞

[marcospassos]

Yeah, that's almost exactly the snippet I'm running (we're using the lib for Bayesian A/B computing):

const samples = [
    ...iterator(1001, 1, {
         iter: 10000,
         seed: 124,
     }),
];

I got fast results for any other case where alpha == beta, or beta > 1. So it seems something related to the sample3 case.

[kgryte]

Update: okay. I am seeing the slowdown, as well. Ran the benchmarks with alpha = 1001 and beta = 1 and got the following results locally:

TAP version 13
# @stdlib/random/iter/beta
  ---
  iterations: 1000
  elapsed: 0.100106084
  rate: 9989.40284188921
  ...
ok 1 benchmark finished
# @stdlib/random/iter/beta
  ---
  iterations: 1000
  elapsed: 0.109215775
  rate: 9156.186457496638
  ...
ok 2 benchmark finished
# @stdlib/random/iter/beta
  ---
  iterations: 1000
  elapsed: 0.074298839
  rate: 13459.160512588896
  ...
ok 3 benchmark finished
# @stdlib/random/iter/beta::iteration
  ---
  iterations: 1000
  elapsed: 0.799214944
  rate: 1251.2278549185887
  ...
ok 4 benchmark finished
# @stdlib/random/iter/beta::iteration
  ---
  iterations: 1000
  elapsed: 0.783351437
  rate: 1276.566241876952
  ...
ok 5 benchmark finished
# @stdlib/random/iter/beta::iteration
  ---
  iterations: 1000
  elapsed: 0.805283159
  rate: 1241.7992215828767
  ...
ok 6 benchmark finished
#
1..6
# total 6
# pass  6
#
# ok

So, extrapolating a bit, the 8 seconds you are seeing makes sense.

[marcospassos]

@kgryte this a simple snippet that reproduces the issue:

console.time('test');

[...iterator(1001, 1, {
    iter: 10000,
    seed: 124,
})];

console.timeEnd('test');

  console.time
    test: 4840 ms

This is the best stat library for JS we've seen so far. Congrats on the great work! Unfurntanatelly, this performance is not only unacceptable for our use case but also a security issue. Do you see any workaround?

[kgryte]

@marcospassos You are correct that the sample3 case is the issue. In particular, L53 and the reason why becomes evident from the preceding lines:

x = pow( u, 1.0/alpha );
y = pow( v, 1.0/beta );
xy = x + y;

For u = 0.5 (the average value), alpha = 1001, and beta = 1, we have

In [148]: base.pow(0.5,1/1001)
Out[148]: 0.9993077849656112

In [149]: base.pow(0.5,1/1)
Out[149]: 0.5

In which case, we are going to end up tossing away many values until finding something that meets the exit criterion.

As for workarounds, one option is that you can provide a different PRNG to the iterator function. For example, if you don't care about seeding, you can provide the built-in Math.random.

it = iterBeta(1001,1, {iter:10000,seed:1234,prng:Math.random});

Otherwise, if you need seeding, you can try using an LCG generator which should be faster, although the randomness quality is not as good. There is also @stdlib/random/base/minstd-shuffle, which is also fairly fast and has better PRNG qualities.

Otherwise, a few other options:

1. You can try preallocating your array, possibly even via an array pool in order to avoid allocating so much memory each time.
2. You can try using the strided array API I mentioned above, in order to avoid creating and using an iterator.
3. You can consider decreasing the number of iterations (i.e., size of the allocated array).

Otherwise, @Planeshifter you have any other suggestions?

[marcospassos]

Well, seeding is essential for us, as well as good randomness. I tried pre-allocating the array, but it resulted in almost no performance improvement. For bayesian experiments, 10.000 simulations are like an industry standard. So I don’t see any viable alternative, at least for us :(. Hope @Planeshifter can help out here

[kgryte]

@stdlib/random/base/minstd-shuffle could be a good option. While LCGs have a shorter period and the output can have "structure", they are good enough for many use cases.

[marcospassos]

But do you believe the problem is related to the PRNG? Cause I think it's more connected to the distribution algorithm. Libraries like SciPy or jStats doesn't have the same issue (maybe a good source of inspiration for the solution)

[kgryte]

@marcospassos One thing to note is that, looking at the PDF of the beta distribution for alpha = 1001 and beta = 1, you have a very skewed distribution.

jStat always uses gamma random variates. This may or may not produce good results, depending on the regime.

We use a similar algorithm as NumPy for generating random variates for the beta distribution. However, we are investigating, and are looking into including a special case.

Will report back on our findings.

[marcospassos]

To give you more context, such a skewed distribution is an actual use case for uninformative prior.