randn performance - implement Ziggurat algorithm

[ViralBShah]

This is the julia implementation of randn, which is exactly the same as the C version. It runs about 6 times slower than the C version - even if I remove all the global variable stuff altogether.

jl_randn_next = -42.0
function randn()
    global jl_randn_next

    if (jl_randn_next != -42.0)
        s = jl_randn_next
        jl_randn_next = -42.0
        return s
    end

    s = 1.0
    vre = 0.0
    vim = 0.0
    while (s >= 1.0)
        ure = rand()
        uim = rand()
        vre = 2.0*ure - 1.0
        vim = 2.0*uim - 1.0
        s = vre*vre + vim*vim
    end

    s = sqrt(-2.0*log(s)/s)
    jl_randn_next = s * vre
    return s * vim
end

Here is the performance analysis for a million calls. I even tried using the libm log, but only sqrt made a difference. I think Jeff has to tell us whats going on. I thought we were not that much slower than C.

randn (fully in C): .06 sec
randn (fuly in Julia, with libm sqrt): 0.38 sec - 6x slower
randn (fuly in Julia, with libfdm sqrt): 0.59 sec - 10x slower

[JeffBezanson]

I see a 3x difference (I'll focus on the fast sqrt case, since that's a simple fix).

The sources of overhead are the function call, and boxing the return value. I tried to separate these.

times for 1e6 iterations:

C randn: 0.05454802513122559

original julia randn: 0.15975308418273926

return nothing instead of value: 0.12742090225219727

no function call (1e6 loop inside randn): 0.10554385185241699

making jl_randn_next a local variable: 0.07610011100769043

I think there is still a gap between llvm and gcc -O3 performance. In the core arithmetic part of the function, the disassembly is clean, we are generating good code. I guess I can gradually work on call overhead. And I have a planned optimization to eliminate the return value boxing when the function is called in a well-typed context. This can probably only get us to the 0.1055 range, or 2x slower.
But, long story short, there aren't any glaring performance bugs in the generated code.

[StefanKarpinski]

So in other words, while it would be nice for us to eventually match C here, for now we should just call the C version and be done with it.

[ViralBShah]

Once we move to dsfmt, this will require a bit more work. 2x slower than C is a great place to be at.

My take is that if its only 2x slower, we can leave it in julia.

-viral

On Jul 8, 2011, at 12:24 PM, StefanKarpinski wrote:

So in other words, while it would be nice for us to eventually match C here, for now we should just call the C version and be done with it.

Reply to this email directly or view it on GitHub:
#96 (comment)

[StefanKarpinski]

It's great that we could be at 2x the speed of C if Jeff did a bunch of fantastic optimization work overnight. That's not going to happen immediately, however, and in the meantime, we can be at 1x the speed of C just by using the C version. What's the point of not doing that?

[ViralBShah]

Ok, I think that the call overhead is dominating in many places. I have replaced sqrt and all rand stuff with faster versions, but don't see speedups in the randmatstat benchmark. Perhaps we will see speedups after your optimizations to call overhead.

-viral

On Jul 8, 2011, at 11:42 AM, JeffBezanson wrote:

I see a 3x difference (I'll focus on the fast sqrt case, since that's a simple fix).

The sources of overhead are the function call, and boxing the return value. I tried to separate these.

times for 1e6 iterations:

C randn: 0.05454802513122559

original julia randn: 0.15975308418273926

return nothing instead of value: 0.12742090225219727

no function call (1e6 loop inside randn): 0.10554385185241699

making jl_randn_next a local variable: 0.07610011100769043

I think there is still a gap between llvm and gcc -O3 performance. In the core arithmetic part of the function, the disassembly is clean, we are generating good code. I guess I can gradually work on call overhead. And I have a planned optimization to eliminate the return value boxing when the function is called in a well-typed context. This can probably only get us to the 0.1055 range, or 2x slower.
But, long story short, there aren't any glaring performance bugs in the generated code.

Reply to this email directly or view it on GitHub:
#96 (comment)

[StefanKarpinski]

What's the problem with just using the C version?

[ViralBShah]

Already done before you wrote. See commit 411bc0b. And for bonus, it is much faster than before. Originally, 1 million randn() took 0.06 seconds, and now it is 0.04 seconds for me.

-viral

On Jul 8, 2011, at 6:59 PM, StefanKarpinski wrote:

It's great that we could be at 2x the speed of C if Jeff did a bunch of fantastic optimization work overnight. That's not going to happen immediately, however, and in the meantime, we can be at 1x the speed of C just by using the C version. What's the point of not doing that?

Reply to this email directly or view it on GitHub:
#96 (comment)

[StefanKarpinski]

Excellent :-)

[ViralBShah]

No problem, so long as the pressure for higher performance for julia code is always in our face. Otherwise we are no different from R or Matlab, where you have to write a lot of stuff in C for performance.

-viral

On Jul 8, 2011, at 7:20 PM, StefanKarpinski wrote:

What's the problem with just using the C version?

Reply to this email directly or view it on GitHub:
#96 (comment)

[StefanKarpinski]

Well, I think we have to strike a balance. For a lot of code being 2x slower than C is perfectly fine. But for really core things like BLAS, FFT, and random number generation, we need to be really fast. And one of the huge strengths of the language is how simple it is to provide a correct but relatively slow general implementation of some operation as a fallback, and then plug in highly optimized special cases transparently so that the user never has to think about it. It lets us present a consistent abstraction without sacrificing performance. Every other language has to pick one or the other.

[StefanKarpinski]

Can we close this now? Seems like the practical effect is taken care of, which is what's important for v1.

[JeffBezanson]

The goal is not so much to eliminate C as to reduce the need. Yes, it might happen that a user needs to write something in C for performance. But it would be silly if the thing they needed to rewrite in C was our own random number generation library.
Remember, even C programmers call a BLAS written in assembly language.

[StefanKarpinski]

I'm confused by this answer. I'm not arguing that we shouldn't have performance as close to C as possible. We definitely should. However, I just don't get what's wrong with using high-performance C code that we're never going to beat like libMT, FFTW, etc. I feel like we're talking past each other here though.

[ViralBShah]

I guess eventually, we need to implement the Ziggurat algorithm, rather than using the Box Mueller transformation for generating normal random numbers:

http://en.wikipedia.org/wiki/Ziggurat_algorithm

Some good analysis on ziggurat performance:

http://seehuhn.de/pages/ziggurat

[ViralBShah]

A textbook description and implementation of the Ziggurat algorithm (and others) by Cleve Moler. Matlab also uses this algorithm for randn().

http://www.mathworks.com/moler/random.pdf
http://www.mathworks.com/moler/ncm/randntx.m

http://www.mathworks.com/company/newsletters/news_notes/clevescorner/spring01_cleve.html

[ViralBShah]

Immediate randn performance is addressed with the Box-Mueller implementation in jl_random.c. However, this is 4x slower than rand(), and needs to be faster for a lot of applications. I have updated the subject of this issue to reflect the implementation of Ziggurat, changed milestone to v2.0, and assigned it to myself for now.

[StefanKarpinski]

Why do we need to use Ziggurat over Box-Muller? Wikipedia of course has detailed descriptions of each, but does no comparison of them and why you might prefer one over the other. From Cleve's writeup, I'm gathering that Ziggurat is less expensive in that you're less likely to reject values when generating random numbers. Is that the reason?

[StefanKarpinski]

Ah, I see the performance analysis now. Makes sense. Seems like the best performance overall would be a high-quality Ziggurat implementation in C.

[ViralBShah]

Yes, Ziggurat just requires two random numbers and some cheap computation, whereas Box-Mueller requires two random numbers, and a sqrt and log, which are expensive. Also, Ziggurat rejects much fewer values.

I did a bunch of digging around - and it seems that it is widely accepted that ziggurat is the way to go.

-viral

On Jul 9, 2011, at 8:51 PM, StefanKarpinski wrote:

Why do we need to use Ziggurat over Box-Muller? Wikipedia of course has detailed descriptions of each, but does no comparison of them and why you might prefer one over the other. From Cleve's writeup, I'm gathering that Ziggurat is less expensive in that you're less likely to reject values when generating random numbers. Is that the reason?

Reply to this email directly or view it on GitHub:
#96 (comment)

[ViralBShah]

Implemented in commit 7ff80c9

[stevengj]

Too bad we didn't keep Box-Muller as a fallback. The advantage of Box-Muller is that it is precision-independent, so it would be usable by user-defined types (like DecimalFloats JuliaMath/DecFP.jl#137 or DoubleFloats JuliaMath/DoubleFloats.jl@14d9131) once they implement rand().

[StefanKarpinski]

Could always be added back as a generic fallback. It’s there in the git history somewhere and probably not that hard to recreate.