Models of benchmark performance for software refactoring?

Modeling
1

I’m writing to ask if people have recommendations for models of profiling/benchmarking.

The data I have at hand is 25K paired observations of new code vs. old code timing (or old code vs. old code for a control). There is also a new code vs. new code data set to measure false positives. Each experimental condition comes with covariates about things like whether the test was single, double, or quadruple precision, whether it involves complex numbers or not, whether it is a matrix, array, or scalar operation, etc. etc. There are 12 paired replicates of each experimental condition. Other benchmarking setups may vary, but the general organization will be the same.

The covariates I have are either binary (e.g., does it involve a complex number, is the operation stripped over matrices) or finite categorical (which test suite it is in, which position it was among the 12 replicates, etc.). Some can be treated as continuous, like problem dimension.

There are two obvious choices here. We can model each system independently and then compare in the posterior, or we can model the differences in the paired comparisons. Modeling differences is a lot easier because the log of the ratio of timings is a natural scale that removes the actual timing baseline and makes the units ratios of performance. Modeling the actual time that all of these processes take is possible and I have the data broken down for it, but it’s going to be a lot harder, I think. So the model I’m thinking about is something like this for n pairwise observations of the A code timing, t^\text{A}_n and B code timing, t^\text{B}_n.

\log(t^\text{A}_n / t^\text{B}_n) \sim \textrm{Student-t}(\nu, \alpha + \beta \cdot x_n + ..., \sigma)

The \ldots is for the varying effects I don’t know how to easily write down in a sketchy way.

Edit: (oops, fat-fingered that out too early). Here’s some of the literature I’ve found so far:

I would love to get feedback on what I should actually be reading, especially if there’s someone who understands Bayesian statistics and has worked in this area. I’m OK reading the frequentist stuff, too.

2

To make this concrete, here’s my current Stan model.

data {
  int<lower=0> N;
  vector<lower=0, upper=1>[N] has_scalar;
  vector<lower=0, upper=1>[N] has_matrix;
  vector<lower=0, upper=1>[N] has_strided_slice;
  vector[N] log2_size;
  int<lower=0> K_suite;
  array[N] int<lower=1, upper=K_suite> suite;
  int<lower=0> K_value_type;
  array[N] int<lower=1, upper=K_value_type> value_type;
  int<lower=0> K_round;
  array[N] int<lower=1, upper=K_round> round;
  int<lower=0> K_run;
  array[N] int<lower=1, upper=K_run> run;
  int<lower=0> K_operands;
  array[N] int<lower=1, upper=K_operands> operands;
  
  vector[N] cpu_time_log_ratio;
}
parameters {
  real alpha;
  real beta_scalar;
  real beta_matrix;
  real beta_strided; 
  real beta_size;
  sum_to_zero_vector[K_suite] gamma_suite;
  sum_to_zero_vector[K_value_type] gamma_value_type;
  sum_to_zero_vector[K_round] gamma_round;
  sum_to_zero_vector[K_run] gamma_run;
  sum_to_zero_vector[K_operands] gamma_operands;
  real<lower=0> sigma;
  real<lower=0> nu;
}
model {
  alpha ~ normal(0, 5);
  beta_scalar ~ normal(0, 2);
  beta_matrix ~ normal(0, 2);
  beta_strided ~ normal(0, 2);
  beta_size ~ normal(0, 2);
  gamma_suite ~ normal(0, 2);
  gamma_value_type ~ normal(0, 2);
  gamma_round ~ normal(0, 2);
  gamma_run ~ normal(0, 2);
  gamma_operands ~ normal(0, 2);
  sigma ~ lognormal(0, 1);
  nu ~ lognormal(log(2), 2);
 
  cpu_time_log_ratio
    ~ student_t(nu,
		alpha
		+ beta_scalar * has_scalar
		+ beta_matrix * has_matrix
		+ beta_strided * has_strided_slice
		+ beta_size * log2_size
		+ gamma_suite[suite]
		+ gamma_value_type[value_type]
		+ gamma_round[round]
		+ gamma_run[run]
		+ gamma_operands[operands],
		sigma);
}

This model takes about 30 minutes to fit robustly on my new Mac Studio (M3 Ultra) at home and about 45 minutes to fit on my older Mac Studio at work (M2 Max). The new sum-to-zero structures are just amazing and I’d highly recommend trying them out for identifiability. Ditto for the new Mac memory architecture.

If you have suggestions for improvements or alternatives, please let me know. I like the approach of Chen et al. for Julia in that it models an optimal time and then system things that can cause slowdown directly.

Hopefully I’ll be able to share some real data soon.

3

This is quite the coincidence; I just got some weird benchmark results back, and decided to pop by here on a whim while I mulled them over.

Your proposed model is decent, though there’s more accurate ones out there. You can think of a processor as a collection of micro operations (“fetch a sequence of bytes,” “decode what those bytes mean,” and so on). A high-level operation code (“add the two numbers in register A and B, store in C”) gets broken up into a series of micro operations, carried out in a specific order. Each of these micro operations takes a minimum number of clock ticks to complete, and that number cannot be less than one, but sometimes they can “stall out” due to external factors (they need to fetch data from memory, and the cache doesn’t have a copy) and take an additional number of ticks. This is a multi-stage Poisson process, and if you add a tonne of them together you get more of an offset Gamma distribution than a log Gaussian/Student. The latter is a decent approximation of the former, though.

Julia’s benchmarking process seems the best of the bunch (using the minimum was wise), but all suffer a bit from spending too much time around frequentism. That framework encourages you to lump together observations to minimize noise. But if you read classic papers on algorithms, they’re full of estimates that theirs performs equivalent to (say) 4n^2 + \frac{3}{4}n + 8. The resulting execution time model would be

t(n) \thicksim \sum_{j=1}^{3} m_j n^{j-1} + \Gamma( \alpha(n), \beta(n) ), \\ E[ \Gamma( \alpha(n), \beta(n) ) ] = \sum_{j=1}^{3} c_j n^{j-1},

for m_j > 0 and c_j \propto m_j. But if that model’s accurate, then the runtime we observe at n strongly predicts the runtime at n-1. So what’s the point in ever repeating a test? Instead think of all the parameters that could influence your algorithm’s performance as axes of a hypercube, and benchmark by taking point samples from that space. Noise still gets eliminated via the structure imposed by the model, but now you cover a much broader range of potential input parameters faster and get finer resolution over an algorithm’s behaviour.

While you’ve already finished your benchmark sampling, the above does suggest that choice of \log( t_n^A / t_n^B ) could be improved. It can work perfectly fine, but only if c_A n^d and c_B n^d are accurate estimators of performance for both algorithms with a fixed value of d. If you didn’t pick large enough values of n, or d isn’t the same for both algorithms, that starts to fall apart. We tend to pretend the term with the highest degree is the only relevant one, yet when reaching for sort algorithms we prefer the O( n \log n ) Quicksort over the O( n ) Radix sort.

You can see me figuring all this out on-the-fly in the last dozen pages of this pre-print:

Haysn Hornbeck, Fast Cubic Spline Interpolation , January 25, 2020. ArXiv, 2001.09253v1.

I missed the additive portion back then, alas, and I probably should have used the mean instead of the mode. If I do publish what I’m benchmarking currently, though, I’ll get a chance to fix that.

4

Ack! I really botched that equation! It should have been something more like

Pr\left( m_j ~|~ t(n) \right) = \Gamma\left( t(n) - \left(\sum_{j=1}^3 m_j n^{j-1} \right) ~|~ \alpha(n), \beta(n) \right), \\ E\left[ \Gamma( \alpha(n), \beta(n) ) \right] \propto \sum_{j=1}^3 m_j n^{j-1} .

That alas will be a huge pain for cmdstan’s NUTS implementation, as there’s a high chance of negative values being fed into the Gamma distribution and initialization failing. Solutions include using Pathfinder to initialize, as it copes better, or scanning the data for an estimate of m_j and plugging that into the offset/multiplierarguments for each. A third is to return to the approximation

\log t(n) \thicksim \mathcal{N}\left( \log \left( \sum_{j=1}^3 m_j n^{j-1} \right) , \sigma \right),

or swap in Student’s t. While I’m back making corrections, I forgot to mention that

E \left[ ~ \left(\Gamma( \alpha(n), \beta(n) ) - E[\Gamma( \alpha(n), \beta(n) )]\right)^2 ~\right] \propto \sum_{j=1}^3 m_j n^{j-1}

under the model I outlined above.

5

Thanks, @hjhornbeck!

I think of what I was proposing here as the obvious baseline homework answer you’d give if you got assigned this in a Bayesian stats class. That’s why I was hoping I could get some suggestions for better literature to read on how to improve the model.

I’m a computer scientist and spent 10+ years in industry writing production code for speech recognition, search, and NLP before starting to work on Stan. So I’m pretty familiar with how both CPUs and compilers work (though they’re both now so complicated that intuition largely fails me beyond trying to maintain memory locality and reduce branch points). The latency/throughput tradeoff in building speech rec systems really drove home the issues you bring up about there being baseline timings then shocks of some kind (for me, in 2000, it was the garbage collection shocks of the early Javascript implementations, so we wound up building up and tearing down the virtual machine every query to avoid long latency at the cost of halving our throughput).

My guess is that there’s not much signal we’re going to be able to get out of the bigger shocks. So my plan now is to not model them. In the data I have, if I just look at the standard deviation of all the log timing ratios, there are values that are dozens of standard deviations away from the mean (all to the high side, of course). That is, my goal isn’t to get a good statistical model of how long something is going to take, but rather to judge whether a putative performance improvement is worth merging. For the latter problem, I’d really have to figure out those mixtures of potential delays.

I believe one of the papers I mentioned (or maybe one I found elsewhere) was all about trying to do experimental design to minimize compute needed for these benchmarks. This particular example is from profiling a system like Eigen, so it’s not going to get too out of control combinatorially, compared to say, testing an operating system or compiler. But I’ve so far skimmed through only half a dozen papers and that’s all I know about modeling optimization beyond what I’ve learned in the field over the last 30 years. So any useful literature on modeling this problem would help.

I couldn’t figure out what the c, d, n are in in your discussion of \log(t^\text{A}_n / t^\text{B}_n). In the notation, n is just the index of the experiment, so I’m assuming you meant the total number of trials by n (the usual convention outside of Stan code). In numerical algorithm analysis, the term with the highest degree we write is usually the error because we’ve gotten the earlier order terms included, and when those are additive, you usually drop them and concentrate on the lower order terms being more informative. But these are usually least meaningful because \epsilon < 1. So I wasn’t quite sure what the analogy with quicrksort was getting at, either.

6

My suggestion would be to consider heteroskedasticity here. If any of the covariates you mention can plausible increase the uncertainty of the difference in speedup, that could be valuable to estimate: new version could, for example, mostly speed up complex numbers but sometimes slow them down?

7

The model already has covariates for things like is-complex, so the fits will diagnose which things got more expensive and which got cheaper. The total prediction is the sum of the relevant varying effects and the fixed effects (like dimensionality). Was there something else you had in mind to deal with the variation in test makeup?

8

What I meant is also making the sigma vary with the covariates. Like for example:

vector[N] log_sigma = log_sigma_base
+ delta_scalar * has_scalar
+ delta_matrix * has_matrix
+ delta_strided * has_strided_slice;

cpu_time_log_ratio ~ student_t(nu, mu, exp(log_sigma));

This – if it fits (which can be tricky with heteroskedastic student t) – would have the benefit that delta are higher when a specific feature has a more variable effect. Say, matrices got sped up by 2 on average, so their effect on mu is log(2), but in practice the speed up is from x0.5 to x4.0, which would show up as high delta. This is quite useful – it can indicate very uneven gains, multitude of unhandled cases, etc.

9

Thanks—I should have figured out that’s what you meant when you said “heteroskedastic.” This can be tricky to fit even without the tails of the Student-t. Any time you allow more than one explanation for the data, identifiability becomes a problem. And now larger values can be (a) higher mean, (b) more dispersion, (c) more variance. I can give it a try, though, and see if I can fit it. So far, the model’s very well identified with 20K+ paired observations.

10

… I’m pretty familiar with how both CPUs and compilers work (though they’re both now so complicated that intuition largely fails me beyond trying to maintain memory locality and reduce branch points).

I am too, so I was fairly astonished to find myself writing this six years ago:

This provides consistent evidence for an unexpected conclusion, that floating-point division does not carry
the same penalty for the Intel i7-7700k that it once did for older processors. Of the newer algorithms,
the fastest is the one which calculates the division last, without loading it from memory, and converting
a division into a multiplication in Numerical Recipes’ algorithm resulted in a loss in performance. More
surprisingly, [one algorithm variant’s] two additional memory access result in a performance improvement over
[the original algorithm, with division as multiplication].

Division was faster than multiplication, one of the most heavily optimized instructions? Adding extra memory accesses lead to faster execution?! My years of experience would have told me both were flatly impossible, and yet for certain Intel processors the numbers said otherwise. Sometimes it’s handy to peel back some of your assumptions and reassess, even if that sometimes retreads old ground.

My guess is that there’s not much signal we’re going to be able to get out of the bigger shocks. So my plan now is to not model them. In the data I have, if I just look at the standard deviation of all the log timing ratios, there are values that are dozens of standard deviations away from the mean (all to the high side, of course).

That’s a fair call. My instinct is to instead incorporate them into the model in some way, by spamming nuisance parameters, but that’s not always worth the effort. The fact that all the outliers are on the high side suggests to me that a symmetric distribution like Student’s t would be a poor choice. You could switch to a skew normal, or apply the same transform used to create that skew onto Student’s t.

So any useful literature on modeling this problem would help.

Agreed, and six years ago I didn’t think there was anything useful, for the reasons I outlined previously. Raiding my old citations, though, led me to this paper:

J. J. Yi et al., “A Statistically Rigorous Approach for Improving Simulation Methodology,” The Ninth International Symposium on High-Performance Computer Architecture, 2003. HPCA-9 2003. Proceedings., February 2003, 281–91, https://doi.org/10.1109/HPCA.2003.1183546.

It’s a close match to the testing methodology I thought I invented, and the date in my citation manager tells me I knew of it back then. That in turn was heavily based on a paper published in 1946 about controlling industrial processes, and the output of other researchers expanding those ideas. In hindsight I really should have cited the 2003 paper, even though I wander off in a very different direction when it comes to the analysis. Whoops.

I couldn’t figure out what the 𝑐,𝑑,𝑛 are in in your discussion of \log⁡(𝑡^A_𝑛 / 𝑡^B_𝑛).

c_j was just one of several unknown constants, d the degree of the polynomial approximation, and n a parameter that effects the algorithm’s performance. Think the n in Quicksort’s O( n \log n ) average runtime on random data.

In numerical algorithm analysis, the term with the highest degree we write is usually the error because we’ve gotten the earlier order terms included, and when those are additive, you usually drop them and concentrate on the lower order terms being more informative.

Ah yes, that explains it. Maclaurin series often deal with x < 1, and in that context higher powers lead to smaller values and thus are safer to ignore. I was instead thinking in terms of algorithm performance, where n >> 1 and thus the opposite is true. But my use of n (as a parameter effecting performance) didn’t match your use (as the total number of experiments conducted), which only added to the confusion.