I thought I’d share a comment I made on Perry de Valpine‘s post on the NIMBLE blog,
Perry was posting about a paper that tried to compare efficiency of Stan, NIMBLE, and JAGS for some linear modeling problems. You can find the link in Perry’s post if you really care and then the paper links to their GitHub. Perry was worried they were misrepresenting NIMBLE and users would get the wrong idea. We decided to take the approach of established companies and simply ignore this kind of criticism. But this time I couldn’t resist as it’s not really about us.
Here’s my comment (lightly edited):
Comparing systems on an equal footing is a well-nigh impossible task, which is why we shy away from doing it in Stan. The polite thing to do with these kinds of comparisons is to send your code to the devs of each system for tuning. That’s what we did for our autodiff eval paper.
I don’t like this paper’s evaluation any more than you do! I’d like to see an evaluation with (a) arithmetic on an equal footing, (b) the kinds of priors we actually use in our applied work, and (c) something higher dimensional than p = 100 (as in p = 10,000 or even p = 100,000 like in the genomics regressions I’m working on now). Then the evaluation I care about is time to ESS = 100 as measured by our conservative cross-chain ESS evaluations that also allow for antithetical sampling (Stan can produce samples whose ESS is higher than the number of iterations; may estimators just truncate at the number of iterations because they don’t understand ESS and its relation to square error through the MCMC CLT). The problem with this kind of eval is that we want to represent actual practice but also minimize warmup to put systems in as favorable a light as possible. In simple GLMs like these, Stan usually only needs 100 or maybe 200 iterations of warmup compared to harder models. So if you use our default of 1000 warmup iterations and then run sampling until you hit ESS = 100, then you’ve wasted a lot of time in too much iteration. But in practice, you don’t know if you can get away with less warmup (though you can use something like iterative deepening algorithms to probe deeper, they’re not built in yet).
One way to get around this sensitivity to warmup is to evaluate ESS/second after warmup (or what you might call “burnin” if you’re still following BUGS’s terminology). But given that we rarely need more than ESS = 100 and want to run at least 4 parallel chains to debug convergence, that’s not many iterations and you start getting a lot of variance in the number of iterations it takes to get there. And things are even more sensitive to getting adaptation right. And also, I don’t think ESS/second after warmup is the metric practitioners care about unless they’re trying to evaluate tail statistics, at which point they should be seriously considering control variates rather than more sampling.
In other words, this is a really hard problem.
I then read Perry’s follow up and couldn’t help myself. I actually looked at their Stan code. Then I had a follow up comment.
I just read their source code. It’s not exactly Stan best practices for statistics or computation. For instance, in mixture_model.stan, there are redundant data computations per iteration (line 25), redundant distributions (also line 25), inefficient function calls (line 31), and a conjugate parameterization inducing extra work like sqrt (line 23). Then in AFT_non_informative.stan, they use very weakly informative priors (so it’s misnamed), there are missing constraints on constrained variables (lines 17, 32), redundant computation of subexpressions and of constants (lines 26, 27), missing algebraic reductions (also lines 26 and 27), redundant initialization and setting (lines 22/23 and 26/27), redundant computations (line 32).
The worst case for efficiency is in their coding of linear models where they use a loop rather than a matrix multiply (LinearModel_conjugate.stan, lines 30–32, which also violates every naming convention in the world by mixing underscore separators and camel case). This code also divides by 2 everywhere when it should be multiplying by 0.5 and it also has a bunch of problems like the other code of missing constraints (this one’s critical—`sigma` needs to be constrained to be greater than 0).
Then when we look at LinearModel_non_informative_hc.stan, things get even worse. It combines the problems of LinearModel_conjugate with two really bad ones for performance: not vectorizing the normal distribution and needlessly truncating the Cauchy distribution. These would add up to at least a factor of 2 and probably much more.
And of course, none of these exploited within-chain parallelization or GPU. And no use of sufficient stats in the conjugate cases like Linear_model_conjugate.stan.
If it takes longer to get the code efficient (more time reading docs, figuring out bottlenecks, installing dependencies, etc) that should also be counted in overall efficiency.
The extreme example is coding the model in assembly. Obviously almost no one does that for a reason. But differences in that up-front cost can still be significant when comparing different software.
Indeed, the whole reason “frequentist” methods for estimation are so popular amongst practicioners is that it is very easy to just plug in your data to a default model and hope the output approximates what you would get from a real model.
Hm, does anyone know if brms could have generated (better versions of) the discussed models?
This is the stupidest thing I’ve read all day, and that even includes some of your other comments.
“Here are my benchmarks of Java against R at multiplying a thousand numbers. R is like a gajillion times slower.”
“If you vectorize instead of looping in R they’re comparable”
“It took you 5 seconds to enunciate that second, which counts against R”
Really? I would have expected it is obvious to anyone who has ever had to do anything with limited resources. It is just more difficult to measure than the other metrics. I can’t imagine not taking that into account for a project though.
If Anoneuoid’s comment is the stupidest thing you’ve read all day then you must be reading some pretty stimulating material. Of course you want to factor in all the time it takes to run a model, which includes the build / coding time. Should one not care about the elapsed time from starting a project to finishing it? That would be news to me
As others below have pointed out while important in a theoretical sense in a practical sense it is a very difficult thing to assess given the great disparity between operators / coders among other reasons. But that doesn’t mean it’s not worth thinking about. We have to think about the hard things sometimes!
Of course it’s worth thinking about. That doesn’t mean you write bad code and then include the inefficiencies in your benchmarking numbers. The concept of an “overall efficiency” benchmark is nonsense. For starters, you’re bringing in your starting levels of familiarity across languages–what if your reader has the opposite biases? What if your reader will only ever write one program in the alien language? Then if it’s hard to learn it’s not worth it. But then if your reader will write all their programs in the alien language, it might be worth taking their time for some performance improvements down the road.
What you do is you write “this syntax has a well-deserved reputation for difficulty; it took the authors 10,000 hours to write the code for this comparison, yadayadayada.” Or you get an subject matter expert to write the code in the unfamiliar language.
People are already writing “bad code and then includ[ing] the inefficiencies in [their] benchmarking numbers.” The OP is an example of this. They do it because that is much faster than writing good code.
> If it takes longer to get the code efficient, that should also be counted in overall efficiency.
The problem is that if you don’t control for inexperience of the operator, the comparison is going to be screwed.
“We took 3 undergrads who had never used a statistical package and who scored close to 115 on a standardized IQ test, assigned each to a different software package, allowed them 6 hours with a tutorial recommended by the publisher of the software, and then timed them as they solved a set of problems” vs. “we optimized all code to current best practices for the platform”. If you do neither, your comparison might be quite biased!
Yea, you’d different estimates depending on expertise. Eg, the OP we saw whoever wrote the model had far less experience with Stan than Bob. So their results do not really apply to him.
I don’t see the point of “controlling” for it though. I would want an estimate for my approximate scenario, not the hypothetical average.
The thing about sufficient statistics caught my interest, are there any case studies on modeling with sufficient stats in Stan?
> The worst case for efficiency is in their coding of linear models where they use a loop rather than a matrix multiply (LinearModel_conjugate.stan, lines 30–32, which also violates every naming convention in the world by mixing underscore separators and camel case). This code also divides by 2 everywhere when it should be multiplying by 0.5 and it also has a bunch of problems like the other code of missing constraints (this one’s critical—`sigma` needs to be constrained to be greater than 0).
Hmm. I didn’t think that loops would be so bad in stan since it’s compiled code, and that multiplying by 0.5 would be substantially more efficient than dividing by 2. How much of a difference would this stuff make?
(Also, I think the filename is okay if e.g. the underscore and camel case denote different separations. e.g. if they had another file called BinomialModel_conjugate etc, then this will improve clarity imo.)
Even sorting horses by speed is hard. Speeds vary by distance and by conditions, etc. etc.
Anoneuoid: I agree that time taken to do the coding is another factor to take into account when evaluating a tool. A better way to measure would be time taken to code vs. run time. Then if you’re building a front-end server app, you care more about run time and if you have a model you are going to run exactly once, you might care more about code time.
But then as someone pointed out, that depends on who’s doing the coding. Stan’s very simple and I have 10 years experience. I also know the compilation chain and compilation targets in C++ and how C++ is compiled and optimized. So I can write fairly optimal Stan code faster than most users are going to be able to remember which order the blocks go in and if it’s “parameter” or “parameters” and “transformed data” or “transformed_data”.
Where I find myself optimizing Stan code beyond our best practices (e.g., going to GPU, multi-core, specialized functions, etc.) is when I can’t get the model to run in the time I have. Say I need to get something that would take a month to fit down to a day.
Mendel: Yup—the operator matters. You can learn to write faster code in a given language, but that might take hours or days or months or even careers depending on how tight the code has to be. That cost gets amortized over everything you then do in the language, but that’s hard to predict. Something better might come along next decade and all that time you (I mean I) spent learning Prolog programming tricks wouldn’t even come back into play until I started on template meta programming in C++ for Stan.
Anoneuooid: Assembly is a tough target. It’s challenging to write C that’s faster than modern C++ and that’s already a lot lower level.
Niko: I’m not sure if those models could even be expressed in brms as written as they don’t use the kinds of priors we recommend. Also, our GLM tools tend to do things like QR decompositions to allow more generic priors and to precondition the posterior (if you remove correlation and scale issues in predictors, it’s a huge win for making the posterior more like a standard normal).
Anonymous: As I wrote above, I think Anoneuoid made an important point. And as they said, it’s also harder to measure the programmer time as it’s not a fixed number with narrow variance like the time to run code on a fixed piece of hardware with fixed compiler settings.
psych defector: The discussion is in its own section in the user’s guide chapter on efficiency: https://mc-stan.org/docs/2_28/stan-users-guide/optimization.html
The key to this kind of tuning is that Stan’s target density only needs to be equal to the posterior log density up to a constant. It doesn’t need to be the joint log densitydensity. That means you can unfold conjugacy and sufficient statistics (closely related through exponential families) in the Stan code. We’re starting to build-in some of these conjugacy relations into compound distributions.
@Zhou: That chapter on efficiency also talks about vectorization in its own chapter and explains what’s going on. If you have a loop over data in Stan, it’s very fast as it’s just compiled down to primitive (double or int) operations in C++. But if you have functions of parameters, you have to autodiff through all the intermediates, which means Stan builds up an expression graph which essentially turns into interpreted code in the backward pass. Compressing this graph through partial evaluation is huge and our vectorized functions do that automatically. But it can be faster to vectorize even double-based operations. Consider evaluating a normal log density normal(y[1:N] | mu, sigma). If we just loop over 1:N and sum up log normal(y[n] | mu, sigma), we wind up calculating log(sigma) in every loop. That log will dominate the whole calculation (the rest is just subtraction, multiplication, and division). So if you calculate it once and then re-use it in the loop, it’s much faster. It’d be even faster if you just multiplied by N on the outside but we haven’t even optimized our normal density that far last time I checked.
Multiplying by 0.5 rather than dividing by 2 is like shaving the heads of your bolts in a race car (i.e., not much). Using repeated vector.transpose() * vector operations rather than matrix multiplication is a huge loss, especially for big matrices. Not only can we reduce the expression graph size to a single element where the adjoint-Jacobin is a simple matrix operation, we can use blocked matrix multiply algorithms that tune themselves to your cache size and exploit SSE and AVX pipelining in the CPU (that’s all in the Eigen code that Stan calls).
It’s a tough battle vs. the world to try to innovate on settle conventions like naming. Of course, it helps to declare your convention so you can follow it. We’ve wound up in Stan with a lot of files named foo-bar_baz.hpp for no reason I can fathom. We’d have solved that problem if we declared a convention and auto-checked the file names as we do with things like code formatting and dependencies.
I worked on one of these comparison papers a while back (although (a) we only had one MCMC platform in the mix; (b) we didn’t focus too much on timing). A bunch of researchers – with expertise distributed across the three platforms we were evaluating (R, AD Model Builder, BUGS) got together and spent a week developing the case studies and coding. Something I really wanted to do, but which would have been *way* too painful to actually implement, was to measure (1) coding times for naive and expert users; (2) run times for naive and expert users. The need/desirability/feasibility of having a dedicated expert to do statistical computation, rather than relying on naive subject-area researchers doing their best on their own, has been rehashed many times on this blog and elsewhere, but I think it is really interesting at least to know what the approximate worst- and best-case timings are (as one example, naive Python code tends to be a whole lot faster than naive R coding, because R punishes you for not vectorizing …)
Bolker, Benjamin M., Beth Gardner, Mark Maunder, Casper W. Berg, Mollie Brooks, Liza Comita, Elizabeth Crone, et al. “Strategies for Fitting Nonlinear Ecological Models in R, AD Model Builder, and BUGS.” Methods in Ecology and Evolution 4, no. 6 (June 2013): 501–12. https://doi.org/10.1111/2041-210X.12044.
Thanks, Ben. That looks really I interesting.
P.S. Pythons loops aren’t exactly speed demons. Last time I timed, Python’s loops were about an order of magnitude faster than R’s, which is still orders of magnitude slower than C++.
Fair enough. But “just” an order of magnitude is a pretty big subjective difference (“the same code took 6 seconds in Python and a minute in R … 1 minute in Python and 10 minutes in R …”)
Don’t get your hopes too much about the R/ADMB/BUGS comparison – the target audience is quantitative ecologist types who are much less computationally savvy then you …
10 minutes in R, 1 minute in Python and 138 milliseconds in Julia.
:-)
So, you’re still a Julia stan these days? (sorry couldn’t resist haha) I’ve thought about checking it out some more, but I honestly won’t have the bandwidth to learn a new programming language for a while..
Chris:
Julia’s our sister language! Stan and Julia were both created from the same Department of Energy grant back in 2011. Julia founder Alan Edelman is my friend; we used to go running and swimming together back when we both taught at the University of California. Also, there’s Stan.jl as well as native Julia implementations of Nuts, Pathfinder, etc.
Neat! It looks pretty slick. I might be able to carve a little time to trial it next summer. Do you use Julia Andrew?
Chris, it’s kinda a joke on the Julia discourse that if your Julia code isn’t fast, you just post it to the discourse asking why and suddenly people will come out of the woodwork and drop the execution time by 3 orders of magnitude. It’s a bit of a fetish for Julia programmers.
I’ve been working on a project in Julia and it’s just a pleasure to use. It helps that I was a Common Lisp fan from way back, since Secretly Julia is a Lisp dialect (with a Matlabish syntax on top).
If you have experience in languages with:
* some kind of ad-hoc polymorphism
* higher order functions
as well as experience with R data structures/tidyverse or numpy/pandas, it should be pretty fast to learn. It’s a pretty straightforward procedural scripting language with analogous libraries to the rest of the data science ecosystem, the one big twist being that it selects which version of a function to execute based on the types of the arguments.
If you do a lot of scientific work (less so traditional software engineering) I’d say it’s worth it. I’m normally very skeptical of trendy new technology, but the fragmentation problem that Julia solves is very clear. Taking python as an example, if you want to get meaningful numerical work done, you’ve gotta use numpy. On the other hand, numpy doesn’t have any kind of derivative support. To make differentiable functions, you’ll have to subclass numpy arrays to add support for dual number types and the algebra they induce. But the base numpy functions won’t know what to do with your arrays, so you’ll have to replace calls to “np.dot(arr1, arr2)” with “dnp.dot(arr1, arr2)”. So you end up with this fragmentation where any time someone needs slightly different functionality, they make their own version of numpy with clones of all the numpy functions. Think pytorch and tensorflows “tensors” (amazingly pretentious name) and JAX’s differentiable arrays. Then, to move between them, you find all instances of “pt.op(tensor1, tensor2…)” and change it to “tf.op(tensor1, tensor2…)”. Because Julia will automatically choose which version of a function to dispatch based on the datatypes involved, you don’t have to do anything beyond changing the datatypes at the source. All the libraries just talk to or override each other’s interfaces seamlessly, so there’s less pointless copy-paste nonsense.
And you can fit statistical models by calling Stan using stan.jl