tag:blogger.com,1999:blog-7056990295646173627.post5476309221087153525..comments2022-12-03T16:26:38.645+00:00Comments on Learning Clojure: Gibbs Sampler : Clojure vs Julia, Java, C, MATLAB, Scala and PythonJohn Lawrence Aspdenhttp://www.blogger.com/profile/02587130870181071109noreply@blogger.comBlogger12125tag:blogger.com,1999:blog-7056990295646173627.post-43042262443516654132014-06-07T07:24:14.577+01:002014-06-07T07:24:14.577+01:00Bit late to this. Once the clojure loops are writt...Bit late to this. Once the clojure loops are written imperatively using pre-allocated double arrays, the clojure speed is about the same as Java and in my tests takes about 45% longer than the code posted by Edmund (who pre-allocates the array). The following library has some good Java PRNGs http://www.iro.umontreal.ca/~simardr/ssj/indexe.html.. I timed the following code taking 45% longer than the Julia:<br />(import (umontreal.iro.lecuyer.rng MT19937 MRG32k3a))<br />(import (umontreal.iro.lecuyer.randvar NormalGen GammaAcceptanceRejectionGen ))<br />(def SSJrngEngine (MT19937. (MRG32k3a.)))<br />(defn samples-loop [^long N ^long thin]<br /> (loop [x 0.0 y 0.0 j thin i N acc (make-array Double/TYPE N 2)]<br /> (let [x (*(GammaAcceptanceRejectionGen/nextDouble SSJrngEngine 3.0 1.0 )(+ (* y y) 4.0))<br /> y (+ (/(NormalGen/nextDouble SSJrngEngine 0.0 1.0)(Math/sqrt (+ 2.0 (* 2.0 x)))) (/ 1.0 (+ 1.0 x))) ]<br /> (if (> j 0) (recur x y (dec j) i acc)<br /> (if (> i 0) (do <br /> (let [^doubles dacc (aget ^objects acc (- N i))]<br /> (aset dacc 0 x )<br /> (aset dacc 1 y )<br /> )<br /> (recur x y thin (dec i) acc)))))))<br /><br /><br />Since Gamma PRNGs call Gaussian and Uniform PRNGs, this story largely reflects the speed of pure Java vs C uniform and normal PRNGs (calling the PRNGs dominate the numerical loops). The expensive parts of the Julia Uniform and Gaussian PRNGs are coded in c: https://github.com/JuliaLang/Rmath/blob/master/src/randmtzig.c. Not surprising that pure java version take about 45% longer which is roughly the original Java vs C comparison. There are also algorithmic differences between e.g. Julia's C Gaussian PRNG and the various Java ones. This one is well regarded in finance: http://home.online.no/~pjacklam/notes/invnorm/#Java<br /><br />A multivariate Gaussian comparison would be interesting given that Clojure matrices are a recent development but Julia has a matrix focus.Anonymousnoreply@blogger.comtag:blogger.com,1999:blog-7056990295646173627.post-79235991549700048222014-03-31T13:06:19.176+01:002014-03-31T13:06:19.176+01:00You probably timed Julia's start up time too. ...You probably timed Julia's start up time too. It's pretty funny how Julia is hyped for performance but the start-up is so slow.Anonymousnoreply@blogger.comtag:blogger.com,1999:blog-7056990295646173627.post-4230990325652769622013-09-16T12:41:21.419+01:002013-09-16T12:41:21.419+01:00There is also an additional quirk if you use leini...There is also an additional quirk if you use leiningen: https://github.com/technomancy/leiningen/wiki/Faster#tiered-compilation By default, leiningen use some JVM options that hurt performance.Dmitry Groshevhttps://www.blogger.com/profile/01265527640482779266noreply@blogger.comtag:blogger.com,1999:blog-7056990295646173627.post-41937748591712877402013-09-16T12:39:07.578+01:002013-09-16T12:39:07.578+01:00You can squeeze a little bit more performance with...You can squeeze a little bit more performance with following:<br />-(set! *unchecked-math* true)<br />-ensure that "thin" and "N" are casted to longs instead of being Objects. In your last example they are Objects;<br />-introduce type Pair instead of constructing a vector every time.<br />This gives an additional 5-10% in my tests.Dmitry Groshevhttps://www.blogger.com/profile/01265527640482779266noreply@blogger.comtag:blogger.com,1999:blog-7056990295646173627.post-89988982777421567982013-08-21T09:04:42.925+01:002013-08-21T09:04:42.925+01:00Indeed, it might sink your boat !Indeed, it might sink your boat !Edmundhttps://www.blogger.com/profile/12286644495692896674noreply@blogger.comtag:blogger.com,1999:blog-7056990295646173627.post-16881349825811496952013-08-21T09:04:15.787+01:002013-08-21T09:04:15.787+01:00Its so heavy to carry you'd never make good yo...Its so heavy to carry you'd never make good your getaway :)Edmundhttps://www.blogger.com/profile/12286644495692896674noreply@blogger.comtag:blogger.com,1999:blog-7056990295646173627.post-1067033655101559122013-08-21T00:49:01.447+01:002013-08-21T00:49:01.447+01:00Probably should, but there's not a great deal ...Probably should, but there's not a great deal of noise in the signal from that. It would also be a good idea to force garbage collection and give the code enough runs through that the JIT has time to optimize it, both of which are noticeable effects. And there are bound to be other factors. <br /><br />But to be honest, I don't really care about anything but the order-of-magnitude! Looks like I can get two of those by stealing Ed's laptop.John Lawrence Aspdenhttps://www.blogger.com/profile/02587130870181071109noreply@blogger.comtag:blogger.com,1999:blog-7056990295646173627.post-86148157487816540172013-08-21T00:45:20.044+01:002013-08-21T00:45:20.044+01:00Is your laptop really 90x faster than my netbook? ...Is your laptop really 90x faster than my netbook? What the hell? I thought Moore's Law had stopped?John Lawrence Aspdenhttps://www.blogger.com/profile/02587130870181071109noreply@blogger.comtag:blogger.com,1999:blog-7056990295646173627.post-88958257370262798252013-08-20T18:36:05.287+01:002013-08-20T18:36:05.287+01:00For comparisons like this why not seed the RNG wit...For comparisons like this why not seed the RNG with the same constant every time.Anonymoushttps://www.blogger.com/profile/12455815708286717299noreply@blogger.comtag:blogger.com,1999:blog-7056990295646173627.post-61149552072469347092013-08-20T17:36:58.129+01:002013-08-20T17:36:58.129+01:00I get to be the baddie - love it.
On my laptop ...I get to be the baddie - love it. <br /><br />On my laptop (time (last (samples-loop 50000))) takes ~1750msec.<br /><br />Lifting and updating from http://dmbates.blogspot.co.uk/2012/05/simple-gibbs-example-in-julia.html here's the equivalent code Julia. Ie no funny tricks with parallelism<br /><br /> ------------------<br />using Distributions<br />function JGibbs(N::Int, thin::Int)<br /> mat = Array(Float64, (N, 2))<br /> x = 0.<br /> y = 0.<br /> for i = 1:N<br /> for j = 1:thin<br /> x = rand(Gamma(3)) * (y*y + 4)<br /> y = 1/(x + 1) + randn()/sqrt(2(x + 1))<br /> end<br /> mat[i,:] = [x,y]<br /> end<br /> mat<br />end<br /><br />@elapsed JGibbs(50000, 200)<br />------------------------------<br /><br />That takes 680msec, which is hard to argue with.<br /><br />Brother JacksonEdmundhttps://www.blogger.com/profile/12286644495692896674noreply@blogger.comtag:blogger.com,1999:blog-7056990295646173627.post-38240293674181409312013-08-20T17:01:31.183+01:002013-08-20T17:01:31.183+01:00Instead of parallelcolt you can look to Clatrix th...Instead of parallelcolt you can look to Clatrix that wraps the JBlas...Alex Otthttps://www.blogger.com/profile/13001951608173211050noreply@blogger.comtag:blogger.com,1999:blog-7056990295646173627.post-74128517759914536342013-08-20T08:44:07.281+01:002013-08-20T08:44:07.281+01:00Brilliant as always, John. Kudos.Brilliant as always, John. Kudos.Baishampayanhttps://www.blogger.com/profile/06191868722015131112noreply@blogger.com