Investigate SIMD support

[pitdicker]

Disclaimer: I don't really know what I am talking about here.

The first step, as far as Rand is concerned, is to generate a small array of values from an RNG. Some RNG's may lend itself well for this, like ChaCha, and possibly HC-128, although that one does not use SIMD. Other options are packing together multiple small RNGs, like Xorshift+ or Xoroshiro+. The result type of such an RNG should be some array, so maybe we can make use of the BlockRngCore trait?

Making use of the generated values can be a more interesting problem. This involves the distributions code.

Converting to floating point: the Standard uniform distribution should not cause any problems, and neither should the range code.

I suppose the other distributions are going to give trouble, at least anything that needs branches or a look-up table (i.e. the current ziggurat code) is not going to work well. And what should happen when an algorithm needs to do rejections sampling, and one of a couple of SIMD values get rejected? Return NaN? For most distributions there should exist an algorithm suitable for SIMD. The Box-Muller transform for example looks promising for the Normal distribution, although it would need a library to provide functions like sin and cos.

Integer ranges: I don't expect the widening multiply method to work for integer range reduction, and division or modulus also do not seem to be available. And I suppose we have no way to indicate a value is rejected because it falls outside the acceptable zone and can cause bias.

Now I don't know how far we should go, and if Rand should provide any implementations at all. But I think we should make sure the API does not prevent Rand from being used in this area.

[TheIronBorn]

I’ve actually been playing with this a lot. Let me see about putting something worth scrutiny on github.

…

Sent from my mobile device

On Apr 5, 2018, at 12:58, Paul Dicker ***@***.***> wrote: Disclaimer: I don't really know what I am talking about here. The first step, as far as Rand is concerned, is to generate a small array of values from an RNG. Some RNG's may lend itself well for this, like ChaCha, and possibly HC-128, although that one does not use SIMD. Other options are packing together multiple small RNGs, like Xorshift+ or Xoroshiro+. The result type of such an RNG should be some array, so maybe we can make use of the BlockRngCore trait? Making use of the generated values can be a more interesting problem. This involves the distributions code. Converting to floating point: the Standard uniform distribution should not cause any problems, and neither should the range code. I suppose the other distributions are going to give trouble, at least anything that needs branches or a look-up table (i.e. the current ziggurat code) is not going to work well. And what should happen when an algorithm needs to do rejections sampling, and one of a couple of SIMD values get rejected? Return NaN? For most distributions there should exist an algorithm suitable for SIMD. The Box-Muller transform for example looks promising for the Normal distribution, although it would need a library to provide functions like sin and cos. Integer ranges: I don't expect the widening multiply method to work for integer range reduction, and division or modulus also do not seem to be available. And I suppose we have no way to indicate a value is rejected because it falls outside the acceptable zone and can cause bias. Now I don't know how far we should go, and if Rand should provide any implementations at all. But I think we should make sure the API does not prevent Rand from being used in this area. — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub, or mute the thread.

[TheIronBorn]

I think this might give us some direction https://github.com/lemire/SIMDxorshift, specifically the avx_randombound_epu32 function: https://github.com/lemire/SIMDxorshift/blob/1766a98f4be5e1d2e78d9e2c3623800fbf92931f/src/simdxorshift128plus.c#L90

[TheIronBorn]

Also of note: for multi-rngs we could use several sets of constants at once. More or less alternating between Xorshift64(10, 13, 10), Xorshift64(8, 9, 22), Xorshift64(2, 7, 3), and Xorshift64(23, 3, 24) for example. Might help raise statistical quality.

[pitdicker]

Great if you would want to help out a bit here!

I believe @vks has also been exploring this, for example in https://github.com/vks/xoroshiro/blob/master/src/rng/xoroshiro128simd.rs

I think this might give us some direction lemire/SIMDxorshift

👍 I can't read those function calls yet, but he is doing two widening multplies, even and odd, and keeping the high parts? And without checking against low they will be biased, but this is not always a problem, especially if the range is small.

Also of note: for multi-rngs we could use several sets of constants at once.

We could, and that seems like a good idea. But just to make sure, that doesn't mean all four (/eight?) can use the same seed.

@TheIronBorn Would you be interested in exploring what it would take to have a Distribution take the output of such an RNG, and transforming it wile staying in the SIMD area?

Converting it to a float, or a floating point range are the simplest methods.

[TheIronBorn]

We could, and that seems like a good idea. But just to make sure, that doesn't mean all four (/eight?) can use the same seed.

Are you asking if such an implementation would require each lane to have the same value upon seeding?

[TheIronBorn]

Float gen is pretty simple. Here's a quick gist https://gist.github.com/TheIronBorn/1ba4da07d6e8e05ab20cb780c383dbec

[TheIronBorn]

Added range sampling

[dhardy]

Speeding up block RNGs with SIMD sounds reasonable; anything else sounds application-specific IMO.

Distribution::sample(..) -> T so in theory the trait can support SIMD distributions, given native consumers of the SIMD types. I guess it should be possible to make adapter implementations, i.e. automatically implement Distribution<(f64, f64)> for any type supporting Distribution<f64x2>.

For anything else, e.g. attempting to produce Normal samples in blocks for individual consumption, I guess something like the MutDistribution or CachedDistribution (i.e. equivalent to the old Sample trait) is needed to allow caching — but this significantly increases complexity and I'm not sure it would even be faster in many cases.

One specific case might be worth a specific solution: StandardBlockRng<T>. It may be possible to make a wrapper around any BlockRngCore which immediately applies the Standard distribution to convert its results buffer to an array of type T, then draws values from this buffer. In theory this might allow some optimisation between the RNG (generation) and the distribution (conversion) steps, although probably whatever loop the RNG step requires to fill the buffer would prevent this. This has the added advantages of a common design with & without SIMD and that it can simply implement Iterator<Item=T>.

[pitdicker]

Thank you for code @TheIronBorn, something to look into. Can I also ask you to create a branch in your fork of rand with the changes? Or did you have another way set up to test things?

Speeding up block RNGs with SIMD sounds reasonable; anything else sounds application-specific IMO.

When I was exploring papers around RNGs and related stuff GPU processing was the current trend, and exploring the possibilities of SIMD also seems to come up regularly. As it is a topic we did not give any attention until know, at least good to explore a bit...

But I have no idea yet how useful it will really be. Especially for some distributions it may be hard for SIMD to provide any benefit over the traditional methods.

One thing that seems like an obvious problem is that Distribution takes an Rng, but Rng exposes no methods to give something like next_u32x4. I think a lot can be solved by adding another method and associated type to RngCore. This type would indicate the native value produced by this RNG. It is something I thought that would be useful for months, but never spoke about...

So our current Xorshift implementation natively generates u32s, Xoroshiro128+ u64s, ChaCha [u32; 16], etc. With an addition like this is is possible to implement RNGs for specific types not covered by the trait. For example one specialized for u16 (maybe for embedded?), for SIMD with a type like next_u32x4, or one designed for floats with f64. I created an incomplete branch with changes. @dhardy Please don't shoot it down instantly...

[pitdicker]

For anything else, e.g. attempting to produce Normal samples in blocks for individual consumption, I guess something like the MutDistribution or CachedDistribution (i.e. equivalent to the old Sample trait) is needed to allow caching

Why would we need caching? To return one value at a time instead of a SIMD type? Or did you write this because the smaller number of sampling algorithms may need some mutable state?

[dhardy]

Caching would just be to return one value at a time, yes; the trouble with that approach is that there is significant overhead storing/retrieving the cached values.

@pitdicker your branch looks interesting, though you could probably just do RngCore: Iterator instead. Also what would the OsRng implementation look like?

The BlockRngCore::Results return values are probably going to be too big for many uses and there is no provided way of breaking the block up into smaller segments; this can be worked around however by choosing an RNG with Item type of the correct size.

There is a significant reason I don't want to add an associated type to RngCore: it is no longer object-safe without specifying the type (the same problem we hit in #307). It might however make sense to use a completely separate trait, possibly just Iterator, and then just implement SIMD distributions as simple functions usable via Iterator::map — this is not well integrated however.

But I hold to my previous point: other than a few cases where SIMD can fit into the existing API (e.g. ChaCha core) I see little point trying to develop the infrastructure without at least one (preferably several) application able to consume SIMD values.

[TheIronBorn]

Here's my branch https://github.com/TheIronBorn/rand/tree/simd/src/simd

[TheIronBorn]

Blending sets of constants for Xorshift seems to only optimize well for machines later than somewhere around Haswell: https://godbolt.org/g/G2x1E2. (Adapted from O'Neill's Xorshift library) I don't have such a machine so I can't confirm actual speed.

Removing the multiply helps.

[pitdicker]

There is a significant reason I don't want to add an associated type to RngCore: it is no longer object-safe without specifying the type

I realized just after posting. So it was a bad idea from me, feel free to shoot it down 😄.

Here's my branch TheIronBorn/rand:src/simd@simd

Nice, I have been playing with it. And SFC seems a pretty good choice here.

I think I have found a way to support SIMD without extra traits like SimdRng and SimdDistribution. The idea is that RngCore::fill_bytes does a simple copy of its result (like u32x4) to the byte slice if they have exactly the same length. And the integer distributions for Standard do the opposite: they initialize an empty u32x4, and pass it to fill_bytes as a byte slice. The optimizer is able to work things out.

I don't think the BlockRng wrapper is going to work well in combination with a simple RNG adapted for SIMD. I expect the bookkeeping to take as much time as just generating a new batch of values. So we may just as well always generate a new batch, and use as much as needed.

I have a proof of concept in https://github.com/pitdicker/rand/commits/simd_support. The code in sfc.rs is very bad and ugly. I am not sure using BlockRngCore in there is useful. But the trick with using fill_bytes seems to work.

But I hold to my previous point: other than a few cases where SIMD can fit into the existing API (e.g. ChaCha core) I see little point trying to develop the infrastructure without at least one (preferably several) application able to consume SIMD values.

I agree we shouldn't turn everything upside down only to support SIMD, especially when we are not sure how useful it will be.

Does implementing Standard and Range for SIMD types fall for you in the category 'develop infrastructure' or 'extension of existing functionality'? Without tests it takes about 125 extra loc: pitdicker/rand@c03deb6...pitdicker:6918e11. The first commit is a necessary change in rand_core to get the inlining right for fill_bytes, and I think that is something that should be in rand_core 0.1. The rest is definitely not interesting until after rand 0.5.

[dhardy]

Nice, I didn't expect fill_bytes to work so well!

Implementing existing traits/distributions is fine — well, we have to look at how much code gets added, but probably fine.