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.