The way we wrote the code in all 5 languages does not use SIMD so the comparison is not fair.
Assuming one can do so (which we have not looked at) we would be happy to include a separate comparison between all 5 languages that use SIMD.
For now, we will add the SIMD Julia code and timing to the appendix.
We already changed the appendix to include @turbo. We would love to credit @Jorge_Vieyra but want his approval first.
Jon was impressed with this proposal, he runs it every day on his website http://extremerisk.org/ (the estimation is done in Julia) and it will significantly speed up his code. And be more elegant.
In response to @jlperla’s point that we don’t need to replace R, python, and Stata because they’re good at what they do, well, if you’re a modest person that may be correct 
But (in the spirit of the Julia manifesto), I’m greedy: I want one environment to do everything: the data cleaning, the linear regressions, and the numerical coding. I want it to be fast, and user-friendly. I want it for free, and I want it easily reproducible, because I’m doing science.
And, you know, it does seem to be all within reach! Just a couple of packages more, and we won’t need to glue different environments together.
For those that, like me, miss the Stata syntax, I’ve put together a package that implements a version of it in Julia: GitHub - jmboehm/Douglass.jl: Stata-like toolkit for data wrangling on Julia DataFrames
(because in Julia this is actually pretty straightforward!) If there’s interest in this from the community, please let me know and I’ll test and develop it a bit more. Active help in the form of tests, issues, PRs, is even more welcome.
Any idea why the numba @jit version is ~20% faster than the Julia version in the original codes?
What optimization @jit is doing that the Julia compiler isn’t doing there? (I confirmed the result here).
Maybe it’s worth a separate thread to get to the bottom of the performance issue? Some discussion on slack suggests it may be to do with our log implementation.
I believe optimized-routines/log.c at master · ARM-software/optimized-routines · GitHub is a slightly faster log than ours. If someone wanted to take a stab at porting it, I would be happy to help.
It’s OK that you don’t want to change the algorithm, and that might prevent simd, fine. Changing the algorithm to enable simd changes the game, OK.
But you’re talking about this as if ‘using simd’ is somehow ‘unfair’ in general, which does not make any sense. What would you do if you compare performance for an ordinary loop, and one language automatically vectorizes it and another language does not? I would call that ‘being a faster language’.
It would seem that the log function takes up significant time. Using @fastmath together with @inbounds results in a ~20% speedup.
using LoopVectorization
nValues = 10_000
y2 = ones(nValues)
function likelihood(o, a, b, h, y2, N)
local lik = 0.0
for i in 2:N
h = o+a*y2[i-1]+b*h
lik += log(h)+y2[i]/h
end
return(lik)
end
orig = likelihood(0.5,0.5,0.5,0.5,ones(nValues),nValues)
println("This is orig = ", orig)
function likelihood_fm(o, a, b, h, y2, N)
lik = 0.0
@fastmath @inbounds for i in 2:N
h = o+a*y2[i-1]+b*h
lik += log(h)+y2[i]/h
end
return lik
end
fm = likelihood_ae(0.5,0.5,0.5,0.5,ones(nValues),nValues)
println("This is fm = ", fm)
function likelihood_turbo(o, a, b, h, y2, N)
lik = 0.0
h_arr = Array{Float64,1}(undef,N)
h_arr[1] = h
@inbounds for i in 2:N
h_arr[i] = o+a*y2[i-1]+b*h_arr[i-1]
end
@turbo for i in 2:N
lik += log(h_arr[i])+y2[i]/h_arr[i]
end
return lik
end
turbo = likelihood_turbo(0.5,0.5,0.5,0.5,ones(nValues),nValues)
println("This is turbo = ", turbo)
julia> @benchmark likelihood($0.5,$0.5,$0.5,$0.5,$y2,$nValues)
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
Range (min … max): 82.595 μs … 206.673 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 85.549 μs ┊ GC (median): 0.00%
Time (mean ± σ): 86.780 μs ± 6.957 μs ┊ GC (mean ± σ): 0.00% ± 0.00%
█
▄▂▂█▆▂▂▂▁▅▂▂▂▂▂▁▁▁▁▁▁▂▁▂▂▂▂▂▂▂▂▂▂▁▁▂▁▂▁▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂ ▂
82.6 μs Histogram: frequency by time 129 μs <
Memory estimate: 0 bytes, allocs estimate: 0.
julia> @benchmark likelihood_fm($0.5,$0.5,$0.5,$0.5,$y2,$nValues)
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
Range (min … max): 70.755 μs … 164.299 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 71.907 μs ┊ GC (median): 0.00%
Time (mean ± σ): 72.928 μs ± 5.002 μs ┊ GC (mean ± σ): 0.00% ± 0.00%
▂ █▃ ▅▁ ▁
█▁██▅█▇▄▁▁▁██▃▆▅▁▁▁▃▁▁▁▁▃▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▃▅▅▅▅▄▄▄▃▃▄▄▄▅▃▄▄▃▄ █
70.8 μs Histogram: log(frequency) by time 95.4 μs <
Memory estimate: 0 bytes, allocs estimate: 0.
julia> @benchmark likelihood_turbo($0.5,$0.5,$0.5,$0.5,$y2,$nValues)
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
Range (min … max): 43.133 μs … 923.304 μs ┊ GC (min … max): 0.00% … 92.41%
Time (median): 43.849 μs ┊ GC (median): 0.00%
Time (mean ± σ): 46.107 μs ± 28.207 μs ┊ GC (mean ± σ): 2.43% ± 3.77%
▇█▆▂ ▁▇▇▅▂▁ ▃▁ ▂
██████▅██████▇▆▇███▆▄▁▄▄▁▇▇▁▁▃▁▄▁▃▅▃▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▄▃▅▅ █
43.1 μs Histogram: log(frequency) by time 60.6 μs <
Memory estimate: 78.17 KiB, allocs estimate: 2.
Sure, feel free to publish it. I just took yours and sped it up with the help from the #performance-helpdesk at the Julia Slack. I will make the github repo over the weekend.
@turbo doesn’t use multi-threading but rewrites some native code to have simd instructions. Numba.@jit doesn’t do it by itself, but LLVM might apply them which I suspect is happening and that’s how you get the speed up, wrt to @inbounds.
The reason I posted the code on Slack was because @turbo was rewriting it wrongly. I modified the algorithm a little bit in order to investigate what was going on. Ideally @turbo should be able to do it without any rewrite.
Yeah, I think we’ll make this a motivating example in the rewrite.
Will be interesting to see whether the best approach is to stack allocate a small buffer, and then chunk the loop based on the buffer side, or, (I’d put my money on) a single “vectorized” loop that scalarizes the dependency chain inside it, packing the result to pass to log. The latter approach may give better ILP, as it shouldn’t be so dominated by the dependence chains.
Sorry for the long response @y.lin but I think the question of whether a single language should be recommended is crucial for economics education.
All the more reason to be extremely careful in the recommendations so they aren’t misunderstood. Maybe I misread the conclusion, but it felt like you were suggesting that the answer to which single language people should standardize on was R. I object to the question itself, not the conclusion. I don’t think people should choose one language because different tasks have different tradeoffs.
If we need to do an unconditional modal language which should be chosen for an unconditional draw of an economic researcher weighted by where they work, what field they are in, etc. then my answer is Stata without a second of hesitation and not R. But that doesn’t really help people choose what languages to learn conditional on what they are trying to accomplish. The intersection of problems where I would choose between R and Julia is essentially empty (whereas something like python could compete head-to-head against both julia and R in different circumstances). I am not even sure that I would benchmark Julia vs. R because the performance of R is very high in all of the places you would want to use it If it is too slow for your needs, it is a good signal that R shouldn’t be considered in the first place.
I think that choosing between R and Stata is a different question - and much tougher than convincing matlab users to try julia given my experience with applied microeconomists. They are extremely productive with stata and I understand their hesitation to switch to something less tailored to their needs and workflow. Where they are willing to give ground is that it is better to educate undergraduates with R instead of Stata.
I disagree. Everyone in economics already learns matlab + stata (and most people in grad school or at a central bank learn dynare). They are used to having different languages for different purposes and can handle having a diverse toolbox.
Even if it was currently true that they stick with one language, it is something we need to convince people to change. There will never be one language to rule them all - whether it is Julia, R, Python, or some future language.
There is an alternative story: lots of professors have tried Python and Julia, and rationally decided it wasn’t worth the trouble to switch until they have a problem that Matlab can’t handle. Matlab comes across as a pretty solid experience in spite of the quirks in the language, and people are willing to pay for that (well… they never personally pay for it in practice). The editor and debugger is polished, has almost no latency, and is completely intuitive and even runtime performance on many problems in practice can surpass Julia because so many methods basically just come down to how well tuned the BLAS/LAPACK implementation is by default, and if you ever want to open a CSV file or run an optimizer in matlab it just works and you never struggle because everything you need is built-in for simple problems.
Of course, there are many horrible things about the matlab language/etc. and gaps in the packages they sell but productivity comes down to the quality of the package ecosystem, its documentation, latency in the development environment, etc. more than the quality of the language. I think that Julia (and python, for different reasons) can make a compelling case to people to switch from matlab, but it isn’t just that people are simply ignorant or old.
I think you need to consider the entire end-to-end workflow of those tasks and look at cleaning/plotting/exporting results/changing the specification/etc… in the whole process. Dedicated and specialized DSLs with tailored development environments - which have no latency where it matters - are tough to beat.
If you haven’t done this sort of research, then it would be tough to convey what is missing because it is about the low-latency environment and specialization of the language. Things like flexibility in renaming variables and trying out different empricial functional forms, not having to code any of the post-estimation statistics that are required for a given method, getting the necessary statistics and clustering of standard errors into latex files you can include directly in your papers, having RAs trained to help you where it is unlikely they will make coding errors, having terse scripts to load and clean data which anyone would understand.
Even R is tough to get many people to use because it is so much more verbose and lower productivity than a specialized language like Stata. The biggest blocker for julia is probably latency since the real performance comparison needs to require latency since the runtime of many operations is reaonably fast and you are constantly changing and manipulating them. For someone who knows Julia and where they aren’t spending all of their time manipulating the regression equation, staring at pretrends for an event study, and examining robust standard errors it is tough to explain how much slower it is.
But for others, these activities are only a small part of their job and Julia would be perfectly fine - even for those tasks.
100% and to defend what @jmboehm and others have done - it is very useful, especially because in many cases the regression may be just one small part of a broader problem, which Julia might be idea for. Same with the DataFrames/QueryVerse/GLM/etc. packages. They are essential.
My point is that if all you are doing is cookie-cutter linear applied microeconomics without any sort of twist, then your productivity radically lower relative to Stata (or even R). That is true even if you are a talented computer scientist who knows a dozen language. If we start suggesting to people they use Julia for those tasks (even if they are experienced programmers), they are going to be miserable. For others who already know Julia, it might be very comfortable and enable them to take the next steps with the project where Stata/R would be limited.
I don’t think that is true in this case. It isn’t driven by ignorance. I know many accomplished coders with a solid background in CS and numerical methods still choose Stata (or R) where appropriate, despite Julia knowledge.
For nonlinear stuff, Turing is great, but it isn’t at the same level of stability and productivity as using R (or Stan if being Bayesian is appropriate and the model fits into its DSL). This is also a grossly unfair comparison considering the amount of $$$ that Stata, R, and Stan developers have pumped into their environments. It is amazing what the Turing and Zygote team have delivered given their relative resources. And there are many cases where Stan or R would themselves fall apart and Turing is a better choice. It depends on the task.
It seems like people here in this conversation believe that Julia is really quite difficult and annoying to use for reading in data sets, transforming them, plotting them, and then fitting 5-20 different models on those datasets. I find this bizarre, because it’s what I literally just did for the last few weeks and it couldn’t have been more enjoyable.
Are there a few warts? Yes, for example there was some situation where by reading in all the census ACS microdata for 5 years I wound up hitting a bug/issue in memory management. My actual datasets were like 200MB after subsampling, but Julia was taking up 10GB due to memory fragmentation or something. But I worked around it in a few minutes and I’ve had to work around other issues in every language I’ve ever used, so I don’t feel like it was different from just the general issue of working around bugs in any language.
But, at its core, Julia is a JOY to use because it’s fast and enables me to do things that just couldn’t be done in R, and most likely also couldn’t be done in STATA though I admit to having zero experience with STATA. Still, I imagine it’d be hair-pulling to for example write a spatial Agent Based model or a differential equations model in STATA.
But even when I’m doing things that 100% could be done in R or STATA, it’s just clean and well thought out in Julia over for example all the BS in R/Tidyverse with its nonstandard evaluation.
TTFX has 100% not been an issue for me, even when I’m not using my custom made sysimage with all the plotting, CSV, DataFrames, etc stuff included.
I mean, to each his own I guess, but seriously, I’m going to start a new thread for head-to-head comparisons of Julia vs X for data analysis… I think we need to dispel some of these myths. The difference between Julia today and Julia in 2019 even is huge, and that’s not true for R, STATA, or Matlab so I feel like people have an outdated idea of how annoying it is using Julia as of this actual moment.
Immature in which sense? We are looking into using it in our company and we have found a few surprises, mostly related to documentation or the fact that algorithms implemented in Julia sometimes are the latest from scientific articles but the old “traditional” ones are not present. Besides that it works quite well. Regarding data pipelines I am surprised since a big chunk of the community focuses on that. Do you have a particular example? Which libraries are you missing?
Please do!! The problem with Julia is not its capabilities but the fact that people cannot find copy/paste-able examples on how to do things since some people need “recipes” to get started.
We ran into several problems where the Julia data libraries were either buggy or lacking functionality, we could not understand the documentation or where the documentation was reading the code. It is quite possible that all of these are okay today, and it is on us to understand the documentation, but since we got it working in Python, there is no need to try Julia again. Since you asked, here is the list (apologies to the community).
AWS.jl. This is not fair because Amazon has written its own library for Python (Boto), but we needed to interact with Amazon and could not make it work with AWS.jl for the services we needed.
Parquet.jl. This is our preferred data format for sharing data between programs. We had a lot of problems with CSV (both size and precision being reduced). Parquet.jl usually works, but there were issues with column types.
SQLite.jl, there was an issue with how it translates sqlite data types into Julia data types.
HTTP.jl, Downloads.jl vs. python Requests. We have to download data from a vendor that needs a password, and we could not make it work in Julia but it worked quite well in Python.
Locks, Proper maintenance and care of multi-threading locks · The Julia Language we ran into a problem with locks and multithreading and we could not make sense of the documentation to fix it.
serialize(), this is writing and reading binary files. Extremely useful when saving data frames and dictionaries and complicated data structures. The relevant Julia page. Serialization · The Julia Language says “In general, this process will not work if the reading and writing are done by different versions of Julia, or an instance of Julia with a different system image.” So one cannot use this in production.
hash and sha. We use hashes a lot to determine whether we have to reestimate, and an easy way is to take a hash of a data frame or dictionary, and if that changes, reestimate. The last time we checked, the built in hash() gives different hashes for different Julia versions, and the SHA functions don’t work with dictionaries.
PS. We put your code on our website, and will now give you credit by name and link. Feel free to get in touch if you want this done differently.
Thank you very much for your thoughtful response. There is much to agree with. Can we suggest you write this up for VoxEU where it will get in front of the economists? Feel free to write to Jon for suggestions and even to collaborate on it.
Yes. And it is unfortunate that this discussion is taking place on this programing language forum and on https://news.ycombinator.com/item?id=32447728
That was not our intention, we tried to do a fair comparison between these four languages, they each have their own advantages and we did not intend to recommend a single one. Rereading the conclusion, I don’t think we did either.
Yes. And would be a worthy topic for a VoxEU if you choose to write it.
That we find curious. The way Jon does the risk calculations for http://extremerisk.org/ used to be R and is now Julia. Both work quite well, and there’s no reason to recommend one or the other for such applications. (Matlab would work but not Python, and certainly not Stata).
don’t think so. We like to ask PhD students and faculty what they have learned and use. Stata is most common, then (often along with) matlab, R or (increasingly common) python, (julia is catching up), and a surprising number use fortran and c.
If you look at the LSE masters courses using computations in relevant fields, it is usually R. And since these feed into the PhD programs, it gives R a lot of weight.
Jon has 2 risk forecasting masters courses and both use R (because of Rstudio). His sample code is in all 4 languages Financial Risk Forecasting code.
Yes. And if the Julia advocates understood that and responded rationally to such concerns, Julia would have a brighter future than otherwise.
Just some comments on the thread since I last commented…
SIMD isn’t entirely manual, or more generally there is “automatic vectorization”. Julia’s @simd lets the compiler know it has more leeway, but the compiler can do it without you writing @simd, or it may decline to do it despite you writing @simd. Since Numba uses LLVM, it also has this. Numba doesn’t seem to have an analog to the @simd hint, but the fastmath=True argument is documented to get around one reason that vectorization may not happen, and it has a way to check the stderr to see if vectorization happened, which is cool. Many C compilers have this, so it’s worth checking for your benchmark report.
The Numba and Julia versions of likelihood do runtime dispatch in the benchmark loop when they take global variables, and Numba’s dispatch scheme is different. If it’s something in the methods themselves, maybe someone (not me) can spot the difference in the assembly;.inspect_llvm() and .inspect_asm() seem to be a Numba function’s equivalent to @code_llvm and @code_native.
But yes a separate thread is warranted; microbenchmarks alone aren’t enough to make any general statement about the relative performance of languages (point 1 of this blog points out that low % differences in small simple methods don’t reflect performance differences of larger applications that use packages), so this definitely isn’t enough either.
To be fair, many do, and they’re developing or contributing to improvements, not only to Julia’s ecosystem but also to its integration with other languages. Any enthusiast of a particular language will promote it wherever feasible, but I think most understand that in the real world, one uses the languages their colleagues use to get a project done smoothly and reliably, whether their favorite language is ready for that or not.
I agree with this, and as evidence I’ll point out that some HN users seemed to think that the article was recommending Julia as the best language.
While “R is the best overall language” sounds like a singular conclusion out of context, they immediately go on to list its negatives too; I believe the intention here was to say “R is the best language on average across several dimensions”. I.e. Julia (or other languages) might be better in particular dimensions, and those might matter the most to you depending on your priorities and needs. (A Jeep might have the best overall rating across several terrains, but that doesn’t matter to me if I’m only looking for a city commute vehicle.)
Thanks for the write-up.
A.
I noticed here, likely added after posted date(?), in case others missed that link: Web appendix for the Python, R, Julia and Matlab language comparison in 2022 on Jon Danielsson's ModelsandRisk.org
p.s.
It is possible to significantly speed up the calculations by taking advantage of SIMD. That is easily accomplished in Julia with @turbo. We have not explored how to do the same in the other 4 languages, but are open to suggestions.
[code with @turbo]Using
@turbospeeds this Julia code up 2.3 times, but that is not a fair comparison to the other languages.
So it looks like your other Julia code is 1.54x slower than C (and slower than Python+Numba), showing in the graphs, is actually not the fastest, with that Julia code 33% faster than C? I think it’s actually fair to Julia to show best case, is it really unfair to the other languages? At least if such is not available there (or much more complicated to implement)?
B.
When using Julia, we use two versions, standard and without bonds checking, @inbounds.
You have a typo “bond”, but I would clarify further:
When using Julia, we use two versions of the code, standard and code applying the @inbounds macro, i.e. disabling bounds-checking.
It felt like you were saying you used non-standard Julia [tools], but using @inbound is very much standard practice, especially in packages. Note though you can disable globally (and that’s the fair comparison to C) with:
$ julia --inline=no
I think actually you can fix your code to avoid bounds-checks without @inbounds or the global option, and maybe merging the loops back is possible (was it really, a needed, part of the speed-up trick?).
B.
On reading compressed CSV files, and you using CSV.jl. It was the fastest package of all language (when subtracting first-use latency).
I’m not sure which is fastest or easiest API, others to consider:
C.
All four languages easily support GPU programming.
A bit surprising, I thought others not as good for GPU use, and at some point at least Python not good at all.
All four languages easily support GPU programming.
I use Matlab for work, and it does indeed not easily support GPU programming, not without buying the parallel processing toolbox.