hardware-counters

Modelling methodology and results

How the profiling data is turned into a performance-prediction model, and what the models actually achieve. Companion to methodology.md, which covers data collection.

The dataset

191 configurations spanning 8 applications × core counts 1–128 × 3–4 problem sizes × 1–4 MPI/OpenMP decompositions.

Each configuration was executed five times with rotating hardware-counter sets (the EPYC 7742 exposes only 5 counters at once). Treating those five runs as five independent samples would give 815 rows and would be a mistake in two ways. It is pseudo-replication — the five runs share a target that differs only by run-to-run noise — but more damagingly, each run only measures its own five counters, so a per-run row has most of its feature vector missing. Measured availability on those 815 rows: l2_miss_per_instr, l1_access_per_instr and prefetch_l2_frac were 0% present and always imputed; ipc 50%; arith_intensity 20%. The average row was ~58% imputed.

Merging to one row per configuration — union of counters across sets, median runtime — gives 191 rows with a 97.4% populated feature matrix, and roughly halves the tail error against the per-run table (p90 2.17 → 1.34), the largest single accuracy difference in the whole project.

Target. The median wall time across counter sets. Set A is always the first run in a job and carries measurable CrayPat warm-up overhead (6.2% slower than B/D/E at identical configurations while recording only 0.6% more cycles), so the median rejects it.

Why leave-one-application-out validation

A random split would place STREAM at 32 cores in training and STREAM at 64 in test — near-duplicates on a smooth scaling curve. A model could score well by interpolating between neighbours without learning anything transferable.

Leave-one-application-out (LOAO) trains on 7 applications and predicts the 8th. The test application is genuinely unseen: different algorithm, different bottleneck, different counter signature.

Error metric. Runtime spans 0.1 s to 215 s, so absolute error would be dominated by the slowest applications. Errors are reported as a multiplicative factor max(pred/actual, actual/pred): 1.00× is perfect, 2.00× is wrong by a factor of two either way. Because this metric is minimised by the median, the constant baseline predicts the median of the target. A baseline predicting the mean of log10(eta) would carry a systematic 15% bias on this left-skewed target (mean −0.169 vs median −0.099), and every comparison against it would be inflated.

Problem framing: the efficiency factor

Regressing log10(runtime) directly on counters fails badly — worse than predicting a constant, with linear models emitting physically absurd values (6 ms for a run that took 6.1 s).

The reason is a measurement detail: CrayPat reports counters for rank 0 only. Rank-0 cycle count correlates strongly with wall time (log–log r = 0.973) but runtime is not simply cycles ÷ peak clock. The models therefore decompose

    runtime = cycles / (f_peak · η)

where η ∈ (0, 1] absorbs stalls, imbalance and communication. Cycles are measured; η is learned.

One caveat on the motivation. The implied clock does vary from 1.98 GHz (GROMACS) to 0.23 GHz (STREAM at 128 cores), but only at STREAM’s largest problem size: at 10M and 40M elements STREAM’s implied clock is a flat 1.85–1.93 GHz at every core count, and only at 80M does it fall (1.91 → 1.08 GHz across 1→128 cores). The minimum of 0.24 GHz occurs at one point, not across the board. A larger driver of low η is simply short runtime — median η is 0.47 for runs under 0.3 s versus 0.83 above 3 s — which is a fixed additive overhead (MPI_Init, I/O, instrumentation) that a pure multiplicative ratio cannot express. Fitting runtime = cycles/(f·η) + b reduces p90 from 1.76 to 1.37, so an additive term is worth adding in future work.

Do the hardware counters earn their keep?

This is the central question for the project and the most robust positive result in it, so it is reported first. Some features are free configuration facts (core, rank and thread counts) requiring no profiling at all, and if those were sufficient the entire five-run-per-configuration campaign would have been wasted. Features were grouped by measurement cost and added cumulatively to a random forest under the same leave-one-application-out protocol.

The answer is yes, but only the full set, and it is worth being precise about what “only the full set” means, because there are two different baselines and they give different-looking but consistent answers.

Against config-only, every counter group helps, and all three survive Holm-Bonferroni within the family of three:

Against the no-feature constant, which is the baseline that actually matters because it is what a practitioner gets for free, only the full set survives:

HL is the Hodges-Lehmann pseudomedian of the paired differences and r_rb the matched-pairs rank-biserial effect size; p-values are Holm-Bonferroni adjusted within each family.

Stated plainly: the full five-counter-set profiling campaign is the only feature configuration that beats a zero-parameter constant, and it does so by 1.5 per cent in median error factor with a large effect size (r_rb = -0.35, 130 of 191 configurations improved). Everything cheaper is statistically indistinguishable from, or significantly worse than, guessing the training median. Free configuration metadata is not enough, one counter is not enough, and one counter set is not enough. This is what justifies the five-set rotation.

Note also that Ridge degrades sharply as features are added (1.146 → 1.293) while the random forest improves, so the relationship between counters and efficiency is not linear.

Figure: ../results/figures/fig_ablation.png.

What the models actually buy you

Having established that the counters carry the signal, the next question is which model extracts it, and the answer is that it barely matters. This section follows the ablation rather than preceding it because the comparison below is close to null, and leading with a near-null horse race would misdirect the reader about where the evidence in this project actually sits.

The honest summary: on typical accuracy a zero-parameter constant is competitive with everything. The random forest beats it significantly but by one percentage point, the gradient-boosted quantile regressor is indistinguishable from the random forest, and the three confidence intervals overlap heavily. The MLP row shows why the paired estimator is the right one: its marginal median (1.076) is higher than the constant’s (1.065), but on the paired differences it wins 119 of 191 and its Hodges-Lehmann pseudomedian is -0.0166, so it is significantly better. The next section explains.

Where models genuinely earn their place is the tail: p90 falls from 1.53 to 1.27 to 1.39. The defensible claim is “models reduce worst-case error by ~17%, while a constant suffices for typical cases”, which is a statement about when profiling is worth doing and is arguably more useful.

Statistical methodology and the precision of the numbers

Every comparison in this document is reported by one routine, analysis/stats_rigour.py, whose output for every key comparison in the project is ../results/model/stats_corrected.txt. This section sets out what it does and why, and how many digits any of these numbers are entitled to.

Direction from the paired differences

A paired Wilcoxon signed-rank test is a test about the distribution of the paired differences, and median(a) - median(b) is not median(a - b), so the direction of a comparison cannot be read from two marginal medians. The two quantities coincide only under conditions this data does not satisfy. Error factors are bounded below by 1 and unbounded above, so their distributions are strongly right-skewed, and a model can shift the great majority of pairs downwards while a handful of large-error pairs hold its own marginal median up. Under those circumstances the paired test and the marginal comparison genuinely disagree, and the marginal comparison is the one that is wrong, because it is not the quantity the p-value refers to. A model that wins 118 of 191 paired comparisons at p = 0.0004 is better, whatever its marginal median says, and a model at p = 0.31 is not “better”, whatever its marginal median says.

The direction is therefore taken from the sign of the Hodges-Lehmann pseudomedian of the paired differences, which is the location estimate the signed-rank test actually targets. The plain median of the differences is reported alongside; the two agree everywhere here. Because the skew runs one way, a marginal comparison is systematically pessimistic about the models on this data.

Multiple comparisons, intervals and effect sizes

The project makes 28 model-versus-baseline comparisons. Unadjusted, at a nominal alpha of 0.05, twenty of them alone would put the family-wise error rate at roughly 1 - 0.95^20, or 64 per cent, so every declared family gets Holm-Bonferroni step-down correction, with monotonicity enforced. Holm rather than plain Bonferroni because it is uniformly more powerful and remains valid under arbitrary dependence, which matters here because the comparisons share a baseline and are therefore strongly dependent.

The correction costs little. Of 28 comparisons, 20 are significant at raw alpha and 19 survive Holm. The single casualty is Cirrus fixed-size LOAO, GBM quantile(0.5) versus constant, raw p = 0.0279 and Holm-adjusted p = 0.0558. So few are lost because the surviving p-values are mostly 1e-4 and below, which the correction cannot touch: the substantive findings are not artefacts of multiple testing.

Every arm’s median error factor carries a percentile bootstrap 95 per cent confidence interval (10,000 resamples; 2,000 for the Hodges-Lehmann statistic, which is O(n²) per resample), and every comparison carries a matched-pairs rank-biserial correlation, the effect size that belongs with the signed-rank test, plus explicit win and loss counts. This matters because a p-value of 0.0002 can sit beside a median improvement of about 0.011 in error factor, which is of very little engineering consequence, and the two facts need to be visible together.

How many digits a random forest median is worth

The same nominal experiment, “random forest, leave-one-application-out, ARCHER2, 191 merged configurations, median error factor”, is run by three scripts and comes out as 1.054497 (retrain2.py), 1.054812 (advanced_nn.py) and 1.055408 (ablation2.py); ../results/model/reconcile_rf.txt reproduces all three exactly. The data, the merge, the 21 feature values, the eight folds and the estimator are byte-identical across the three scripts, because two of them import the third’s loader, and thread non-determinism is excluded because the forest is bit-reproducible across n_jobs settings. What remains is two incidental implementation differences, neither of which is a difference in method.

The dominant term is feature column order. ablation2.py rebuilds the feature matrix in cumulative measurement-cost order rather than definition order, which moves log_t_analytic from last position to fourth. The columns and their values are identical; only the positions differ. But max_features="sqrt" draws candidate feature subsets from a PRNG stream indexed by column position, so permuting the columns builds genuinely different trees even at random_state=0. The second term is the preprocessing chain: a forest splits on thresholds and is therefore invariant to any strictly increasing per-feature map in exact arithmetic, and indeed the sign-log and plain-standardisation chains agree bit for bit, but advanced_nn.py’s RobustSmoothScaler applies 3·tanh(z/3), which saturates. It maps |z| > 12 into a band narrower than 1e-5 around ±3, so values that were distinct collapse onto the same float64 and the split that separated them disappears.

The decisive observation is that neither term matters, because both are smaller than the noise. Sweeping random_state over ten values on the canonical chain gives a mean of 1.0572 with a standard deviation of 0.0049 and a spread of 0.0158. The implementation spread across the whole eight-cell reproduction grid is 0.0030. The seed spread is five times larger, and all three figures lie comfortably inside it. A third source, scikit-learn version, moves the same code from 1.054497 (1.7.2, ARCHER2) to 1.055414 (1.9.0, Cirrus) through a change in tree-builder tie-breaking.

Canonical value. Random forest, leave-one-application-out, ARCHER2, 191 merged configurations, full 21-feature set: median error factor 1.05. Seed-averaged point estimate 1.0572 over ten seeds (sd 0.0049). Bootstrap 95 per cent CI on the median [1.044, 1.072]. The canonical pipeline is retrain2.py’s (median imputation, sign-log, standardisation) applied to build_features output in definition order.

The reporting rule that follows, used throughout, is to quote two decimal places with an interval, because the further digits encode seed, column order and library version rather than anything about performance modelling. No conclusion depends on the choice, since every variant beats the constant with 123 to 126 wins of 191 at raw p of about 2e-4. Any reported forest median is seed-averaged over at least ten seeds, and the scikit-learn version matters enough to pin. The full reproduction grid is in rigour.md.

Learning curve: is 191 configurations enough?

The question of whether more data would help has a different answer depending on what “more data” means, and separating the two answers turns out to be the most directly actionable finding in the project.

The outer loop is the unchanged leave-one-application-out protocol. Inside each fold the training set is subsampled to 25, 50, 75 or 100 per cent before fitting, with the test fold never subsampled, so every point on the curve is scored on the same 191 configurations and the points are directly comparable. Twenty random draws per point, summarised as a median with a percentile bootstrap 95 per cent confidence interval. Two subsampling schemes are run because they answer different questions: random subsamples training rows and therefore measures the density of coverage (“would more runs of the codes I already have help?”), while by-application subsamples training applications and keeps all their rows, measuring coverage of the application space (“would more codes help?”). The second is the question that actually governs leave-one-application-out generalisation.

On the density axis the dataset is saturated. Under random row subsampling the curve flattens: the step from 75 to 100 per cent of the rows buys 0.0037 for the random forest and is slightly negative for the gradient-boosted model, and the 75 and 100 per cent confidence intervals overlap. Collecting more configurations of the same eight codes would buy essentially nothing.

On the coverage axis it is badly under-sampled. Under application subsampling the curve is still clearly descending at the right-hand edge. The final step, from five training codes to seven, is worth 0.0223 for the random forest, six times the corresponding density step, and the 75 and 100 per cent confidence intervals are disjoint. There is no sign of a plateau. The 25 per cent point (two training codes, mean 47 training rows) sits at 1.1333, which is worse than the constant baseline at that size (1.0881), and that is a direct measure of how much the models depend on application diversity rather than on sample size.

Design implication, stated plainly. The profiling budget was spent widening the configuration sweep within each application: three or four problem sizes and up to four MPI/OpenMP decompositions per code, five counter sets per configuration. The learning curve says that budget should have gone on more applications with a coarser sweep each. Density is exhausted; coverage is not, by a factor of six on the final measured step.

Two consequences follow for the rest of this document. The first is that the instability documented above, a random forest median that moves in the third decimal place between seeds, is not mysterious: with eight applications, LOAO has eight folds and every fold trains on seven codes, which sits on the steep part of the coverage curve, so the quantity being estimated is dominated by which seven codes happen to be in the training set. The second is that the cross-platform comparison is partly reframed. Cirrus has six applications where ARCHER2 has eight, so every Cirrus and transfer number in the project sits further left on the coverage curve than the ARCHER2 numbers do. Part of the gap between within-ARCHER2 (1.05) and within-Cirrus (1.12 on the matched data) is a training-set coverage effect rather than a hardware effect, so a conclusion about hardware is here partly a conclusion about sample composition. See crossplatform.md.

Figure: ../results/figures/fig_learning_curve.png. Data: ../results/model/learning_curve.{txt,csv}.

Neural networks

Neural approaches were pursued through three rounds of increasing sophistication. All converge on the same answer.

Round 1 — plain MLP, 40 rows. Significantly worse than Random Forest (11/40 wins, p = 0.019) and indistinguishable from a baseline.

Round 2 — retuned, 147 rows. Diagnosing round 1 found five candidate problems; an ablation showed only some helped. Recovering the per-counter-set samples and lowering an over-aggressive regularisation strength (α 10 → 1, the single largest gain) brought the MLP to 1.069, level with trees. Two proposed fixes were rejected on evidence: bagging was marginally worse, and missingness-indicator features actively harmed the model — each application turned out to have a unique indicator signature, so the columns leaked application identity and wrecked the STREAM fold (1.096 → 1.552). Aggregate metrics hid this; only the per-application breakdown caught it.

Round 3 — published state of the art, 191 merged configurations. Three techniques the tabular deep-learning literature identifies as closing the gap to trees:

• PLR numerical embeddings (Gorishniy et al., NeurIPS 2022): each scalar becomes concat[sin(v), cos(v)] with v = 2πcx and trainable frequencies c ~ N(0, σ²), then Linear+ReLU. Intended to give MLPs the high-frequency capacity they structurally lack — the mechanism Grinsztajn et al. identify as the reason trees win on non-smooth tabular targets.
• Robust scaling + smooth clipping (Holzmüller et al., NeurIPS 2024): median/IQR scaling then 3·tanh(z/3), replacing standardisation.
• RealMLP (same paper, via pytabkit): a pre-tuned MLP with meta-learned defaults, PBLD embeddings and a diagonal weight layer.

Every neural approach in this round loses to the zero-parameter constant, and PLR, the most-recommended technique, is the worst and significantly so. Nested selection barely helped (1.185 → 1.165), confirming this is not a tuning failure. The one exception in the wider set of neural results is the plain MLP (32,16) of the main comparison, which is significantly better than the constant on the paired test (Holm p = 0.0070) despite a higher marginal median; a small effect, 0.017 in error factor, that the marginal medians alone would not show.

The explanation is the data regime, and it is exactly what the literature predicts. With 191 rows the dataset sits well below the ~1000 where these techniques begin to pay, and PLR adds parameters (24 frequencies × 21 features plus embedding layers). More fundamentally, 62% of η values lie in [0.75, 0.90] — the target is nearly constant, so there is little signal for extra capacity to capture, which is precisely why a constant is hard to beat.

The one positive: RealMLP achieves the best p90 of any model tested (1.354) despite a worse median, reinforcing that models earn their keep on the tail.

This is reported as a negative result with an identified mechanism rather than a deficiency of method. Reviving neural networks here would require η to vary more than it does on a single machine — i.e. data from different hardware, not more runs on the same node.

Feature importance under structural missingness

Feature importance on a counter-multiplexed dataset needs more care than it usually gets, and this section shows why with the feature most exposed to the problem: core_energy_frac, the PP0/PACKAGE energy ratio, which RAPL records on every run. On the 191 merged configurations it carries 0.008 impurity importance, rank 15 of 21; on a smaller, differently assembled table of the same measurements it reaches 0.509 and tops the ranking. analysis/importance_artefact.py runs the whole analysis; outputs are analysis/out/importance_artefact.txt and analysis/out/fig_importance_artefact.png.

The ranking on 191 merged configurations (figure: ../results/figures/fig_importance_artefact.png):

The two measures agree in direction (Spearman rho = +0.48) but disagree on detail, which is the expected outcome for 191 rows and is a reason not to over-read any single ranking. They agree on the point that matters: the leading term is the analytic time scale, the rest is instruction mix and cache behaviour, and energy is in the bottom half either way.

How the same feature reaches 0.509

On data/runs.csv, the project’s first table (40 rows, five applications at eight core counts, one problem size each, 17 features; train_model.py), core_energy_frac scores 0.509. Three things compound to put it there, and each is isolated by changing one thing at a time.

A. No run-length covariate. eta is dominated by how long a run lasts, because short runs carry a fixed MPI_Init, I/O and instrumentation overhead (median eta 0.47 below 0.3 s versus 0.83 above 3 s). That table has no log_t_analytic term, so the forest has to split on a correlate of run length instead. Adding log_t_analytic and log_instr_per_rank to the same 40 rows and changing nothing else drops core_energy_frac from 0.526 to 0.148, removing 72% of the importance, and log_t_analytic takes the top slot at 0.582.

B. One problem size per application, at n = 40. In runs.csv each application appears at exactly one problem size, so application identity, problem size and counter signature are perfectly confounded. The correlation between core_energy_frac and log10(eta) is +0.44 on those 40 rows but only +0.088 on the same five applications once several problem sizes each are present. Subsampling the merged data to n = 40 while keeping several sizes per application gives 0.028 (rank 10), so small n alone accounts for only part of it; small n plus the single-size design is what puts the feature at the top.

C. Structural missingness on the 815-row per-counter-set table. This is the effect that generalises. The EPYC 7742 exposes five counters at once, so the feature space is widened by rotating counter sets, and a per-set row therefore has most of its feature vector missing: only 44.1% of feature cells are observed. RAPL sits on separate registers and is appended to every set for free, so energy was among the best-observed features while most others appeared in one set in four.

Median imputation collapses a mostly-missing column to a near-constant, a constant cannot reduce impurity, and the importance mass therefore migrates to whatever happens to be recorded. Measured on that table, impurity importance correlates with per-feature observation frequency (Spearman rho = +0.57, p = 0.007); on the merged table the same correlation is absent (rho = +0.08, p = 0.73). The four fully-observed features take 63% of the total importance there, against 38% on the merged data, and the three features never observed there take 0.000 against 0.078 merged, so real signal is suppressed as well as misattributed.

The control settles it. Re-imposing the per-set missingness pattern on the merged 191-configuration data, averaged over 20 random masks, lifts core_energy_frac from 0.007 (rank 15) to 0.018 (rank 7) and reproduces the inflated ranking at Spearman rho = +0.74, against +0.11 for the unmasked merged data. Nothing about the underlying signal changed; only the observation pattern did.

Why this is worth reporting

Rotating counter sets is the standard way to widen a feature space past a small hardware register budget, and it produces structurally missing columns by construction. Any counter that sits on a separate register, and is therefore recorded in every set, will absorb the displaced importance and look like a discovery. RAPL energy is the usual such counter. We could find no counter-based performance-modelling paper that flags this, despite multiplexing being near-universal in the field.

The safeguards are cheap:

• merge to one row per configuration before interpreting importance at all;
• include an explicit scale covariate, so run length is not attributed to whatever happens to correlate with it;
• vary problem size within application, so importance cannot proxy application identity;
• confirm with held-out permutation importance rather than training-set impurity.

The finding is consistent with the rest of the project: the cross-platform intersection ablation (docs/crossplatform.md, finding 2) shows that dropping both energy features costs 0.4% of accuracy, which is what an importance of 0.008 predicts.

Threats to validity

• 191 configurations, but only 8 applications. Differences of ~0.01 in median error are within noise, and the seed-to-seed spread of a single random forest median is 0.0158, so most inter-model gaps in this document are not resolvable. Only the counter-group ablation reaches convincing significance. The learning curve shows the binding constraint is the application count, not the configuration count: the configuration axis is saturated and the application axis is not.
• Single hardware platform for the results in this document. Every measurement here is from ARCHER2’s EPYC 7742. What is demonstrated is transfer to an unseen application, not to unseen hardware, despite the project title. crossplatform.md addresses the gap and finds transfer works but weakly.
• Significance families are retrospective. The families used for Holm-Bonferroni were declared after the comparisons were run, not pre-registered, so the correction controls the error rate within each family as declared and not across the project as a whole. The 28 comparisons are grouped into six families; a single project-wide family would be more conservative still, though with p-values of 1e-9 and below in several places it would change few verdicts.
• Rank-0 counters only. Load imbalance between ranks is invisible to the features and appears only implicitly in η.
• Two applications remain confounded. GROMACS and OpenFOAM have a single problem size each, so for them application identity and problem size cannot be separated.
• Structural missingness. pat_report omits counters reading zero, and counter multiplexing means most counters are absent from most per-set rows, so gaps are not random. After merging only ~2.6% of cells are missing and the choice of imputation barely matters (median 1.073, zero-fill 1.089, KNN 1.105), but on unmerged data it is severe enough to make feature importance meaningless (see Feature importance under structural missingness). Accuracy is much less affected than interpretation is, which is itself the point: a pipeline can look healthy on error metrics while its importances are meaningless.
• Hyperparameter optimism. Non-nested selection was measured to be optimistic by ~0.058 — larger than most margins reported between models. Results above use nested selection where a search is involved.

Reproducing

module load cray-python/3.10.10
python3 -m venv --system-site-packages $B/venv     # numpy/pandas from cray-python
source $B/venv/bin/activate
pip install scikit-learn scipy matplotlib torch pytabkit

python scripts/analysis/parse_expand.py    # CrayPat dirs -> runs_expanded.csv
python scripts/analysis/retrain2.py        # main results
python scripts/analysis/ablation2.py       # do the counters earn their keep?
python scripts/analysis/advanced_nn2.py    # PLR / RealMLP / ensembling
python scripts/analysis/importance_artefact.py  # feature importance under multiplexing
python scripts/analysis/stats_rigour.py    # the reporting routine, every comparison
python scripts/analysis/reconcile_rf.py    # how many digits a forest median is worth
python scripts/analysis/learning_curve.py  # density versus coverage

Install scikit-learn into a venv with --system-site-packages, not with pip install --user: a user install pulls a newer numpy that shadows and breaks cray-python’s pandas.