np.random.choice without replacement or weights is less efficient than random.sample?

I was interested in choosing k elements out of n without replacement or weighted probabilities, and with k<<n.

I've tested using the following ipython notebook: source / view.

The outcome was that using random.sample(xrange(n),k) was ~50-fold faster than using np.random.permutation(n)[:k], which is the way it is currently coded.

It seems one could basically copy paste the code from np.random.permutation. Should I do that in PR gh-2727? The thing is, this seems like an obvious and simple improvement, but it would leave the random number generator in a different state so it is probably better to do it now rather then later, since someone might want to be able to rely on reproducible results over versions.

Actually, sorry forget that, its not quite that simple of course... you would probably need a point where to switch methods or only optimize for few rather then many elements taken.

Or at list have a flag argument that allows the user to decide if k<<n

@yoavram

@yoavram, honestly I think an extra flag is likely confusing (and if you want to change it the k<<n special case might be more important since there is an existing function to handle the other case). If it bothers you, maybe you can implement this? Maybe even with a simple speed test for some simple switching rule (or likely ie. the R function already has good code to inspire from)... I am not quite sure how its handled but if its soon maybe you can convince to get it into 1.7. to avoid changing the RNG state later, which might need warnings, etc.

@seberg

Going back to this, scipy/scipy#3650 pretty much ironed out the details on how to approach this:

if n < 3k, use np.random.permutation
if n > 3k, use random.sample
@seberg would you like me to make a PR? If so, should I write a test and where?

@yoavram

@yoavram I am happy to see improvements, but take care about the random number stream. There was some thoughts on how to change versions, etc. but I am not quite sure about the state of it right now.
The thing is we must not change the stream of random numbers unless if someone used seeding, etc. I think there was some other discussion about the efficiency of choice as well if you search the issues.

Anyway, tests go just in the tests/test_random.py file in the random subfolder.

NB: the current implementation using permutation is also horribly memory-inefficient for k << n (see #5299).

And right, historically the main reason that we haven't fixed such things is that we haven't wanted to change the implementation of RNG methods because we want results to be deterministically reproducible across versions when someone uses a seed.

The solution to that problem has been discussed extensively and agreement reached on how to handle it, via adding a "version" concept to RandomState -- there's an implementation in #5911.

That PR is good on the actual versioning logic, but is stalled out b/c it doesn't make sense to merge the logic for handling improved RNG methods until we actually have at least one improved RNG method implementation :-). #5911 attempted to do this by improving random.choice actually, but it turns out that handling the k << n case is somewhat tricky because it requires a set data structure. (If you look at the implementation of random.sample, it requires a set in this case. Their algorithm is also quite inefficient when k ~ n, because it's pure rejection sampling. Floyd's algorithm is beautiful and efficient for arbitrary k, but again requires a set data structure. Though Floyd's algorithm does also require a permutation as a post-processing step.)

It's certainly possible to include an efficient set data structure for ints into np.random, which is what we'd need here -- pandas uses a Cython wrapper around khash for such things. Importing khash may be a bit more of a job than you were hoping for, but it is certainly doable :-)

OK, I'm a bit lost. Can we make a list of action bullets of what needs to be done? And also, what's wrong with python's native set?

Todo list:

Implement a better version of choice inside np.random. Ideally one whose memory usage is O(k) and with some benchmarks and tuning to make sure it's efficient across a variety of k and n values.
Implement a framework for managing the "version" of a RandomState object. The basic problem here is that if someone has existing code that does r = np.random.RandomState(seed=0); r.choice(...) then this should behave the same in all versions of numpy -- which means that we'll have to keep the old bad implementation of choice around and fall back to it in some cases. We have a good design for how to handle this, and most of the necessary implementation can be found in Allow multiple versions of RandomState: init' work. #5911. Reading the discussion on that PR should give you a sense of what's going on and what still needs to be done (if anything).

The problem with Python's native set is just that it works with Python objects instead of raw ints, so it's relatively slow and memory hungry compared to a C-level set implementation. On further thought, I guess this might be fine for a first version (though it would be good to see the benchmarks), b/c it's only a constant factor slowdown and -- crucially -- we can always swap out the set implementation without affecting the random stream.

Wow, you weren't kidding about that being complicated :-). Having stared at the code for a bit, I still have no insight into how it works or what the complexity is, but, I am highly suspicious of all the use of floating point math. It seems incredibly unlikely that the Vitter algorithm is correct or consistent on processors that cannot do exact real arithmetic (i.e. all of them). One of the first steps in the first function is to cast N to floating point, but one of our motivating use cases (#5299) is N = 2**64, which is in the range where not all integers are even representable as floating point.

The scikit-learn code doesn't seem to have any big surprises, but agreed that it's a good thing to look at for comparison.

OK, I'm a bit lost. Can we make a list of action bullets of what needs to be done?

Several np.random improvements (e.g. ziggurat for normal, stick breaking for dirichlet, Floyd's for random choice) are blocked by

#2727 (comment):

It's critical that these methods be deterministic (given the seed) in order to allow people to, e.g., reproduce scientific results across numpy releases.

and are waiting for a change like #5911 to be merged. Although I don't need this determinism across releases, there is not consensus to break it in numpy. I'd like "We have a good design for how to handle this, and most of the necessary implementation can be found in #5911" to be true but I don't know what is its actual status.

I'd like "We have a good design for how to handle this, and most of the necessary implementation can be found in #5911" to be true but I don't know what is its actual status.

AFAIK that's its actual status :-). Are you concerned about something about the design? The mailing list thread discussing this (and apparently reaching consensus that it's a good way forward) is here: https://mail.scipy.org/pipermail/numpy-discussion/2015-May/072952.html

The above comment is a little confusing, because AFAIK the main blocker on #5911 getting merged is actually having a good use case (because we don't like merging code that isn't actually used). If choice turns out to be complex and we want to get the changes in #5911 in, then maybe someone should implement a stick-breaking dirichlet to push this through? That's a relatively trivial algorithm. (Disclaimer: based on my vague memories of implementing it some time ago, YMMV :-).)

Sure, N=2**63 then for typical 64-bit platforms. The point is that it is greater than 2**53 the point at which double stops being able to represent integers exactly.

ENH: Add tail shuffle

ENH: Improvce choice without replacement

@mattip

@mattip Fixed in Generator, can't be fixed in RandomState

In [18]: g = np.random.Generator()

In [19]: n=np.iinfo(np.int).max

In [20]: n
Out[20]: 2147483647

In [21]: k = 100

In [22]: %timeit g.choice(n,k,replace=False)
31.8 µs ± 2.32 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

In [23]: %timeit random.sample(range(n),k)
191 µs ± 91.9 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

OK, closing. For performance, use the new APIs

Where is the new API documented?

This is good news, thanks all.

https://www.numpy.org/devdocs/reference/random/index.html

Navigation Menu

Repository navigation

np.random.choice without replacement or weights is less efficient than random.sample? #2764

14 remaining items

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

np.random.choice without replacement or weights is less efficient than random.sample? #2764

Description

Activity

seberg commented on Nov 23, 2012

seberg commented on Nov 23, 2012

yoavram commented on Nov 26, 2012

seberg commented on Nov 27, 2012

yoavram commented on Oct 2, 2015

seberg commented on Oct 2, 2015

njsmith commented on Oct 2, 2015

yoavram commented on Oct 2, 2015

njsmith commented on Oct 3, 2015

njsmith commented on Oct 22, 2015

argriffing commented on Oct 22, 2015

njsmith commented on Oct 22, 2015

rkern commented on Oct 23, 2015

14 remaining items

bashtage commented on Jun 1, 2019

mattip commented on Jun 1, 2019

ryanpeach commented on Jun 5, 2019

ryanpeach commented on Jun 5, 2019

bashtage commented on Jun 5, 2019

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions