Unmasking Algorithmic Energy: A Two-Stage GLM Methodology

Author

Methodology Explainer

Published

June 17, 2026

1. The Challenge: Hardware Noise vs. Algorithmic Truth

When measuring the energy consumption of an evolutionary algorithm, the raw measurements (PKG or delta_PKG) are heavily contaminated by the operating environment. A single execution’s energy footprint is determined by two competing forces:

  1. The Algorithmic Intent: The theoretical cost demanded by the problem parameters (Dimension, Population Size, Max Generations).
  2. The Environmental Friction: Operating system background tasks, CPU frequency scaling, and cache misses that occur randomly during execution.

If we attempt to model energy directly against time and parameters simultaneously, we run into multicollinearity—the algorithm’s parameters inherently dictate the execution time, which confuses the statistical model.

To solve this, we employ a Two-Stage Generalized Linear Model (GLM) to surgically separate the environmental noise from the pure algorithmic cost.

1.1 Data Loading and Preparation

First, we load our 5 independent experimental blocks. By relying on the process_deltas() function from the energyR package, we ensure that static baseline power is subtracted, leaving us with purely dynamic execution metrics (the Sandwich Protocol).

Code 
library(ggplot2)
library(dplyr)
library(marginaleffects)

# Load the Sandwich protocol delta processing logic
source("energyR/R/process_deltas.R")

# Load the 5 experimental blocks
files <- c("data/cec-1.12.4-25.10-cec-mixed-v1-26-Jan-08-30-26.csv",
           "data/cec-1.12.4-25.10-cec-mixed-v2-26-Jan-12-20-59.csv",
           "data/cec-1.12.4-25.10-cec-mixed-v3-26-Jan-17-42-27.csv",
           "data/cec-1.12.4-25.10-cec-mixed-v4-27-Jan-08-09-20.csv",
           "data/cec-1.12.4-25.10-cec-mixed-v5-27-Jan-10-16-05.csv")

# Combine and process the deltas
all_data <- do.call(rbind, lapply(files, function(f) process_deltas(read.csv(f))))

# Ensure parameters are treated as discrete categorical factors
all_data$dimension <- as.factor(all_data$dimension)
all_data$population_size <- as.factor(all_data$population_size)
all_data$max_gens <- as.factor(all_data$max_gens)

1.2 The Co-linearity Trap (Why a single model fails)

Before applying the two-stage method, let’s visualize why time and parameters cannot share the same predictive space without adjustments.

Code 
ggplot(all_data, aes(x = seconds, y = PKG, color = dimension)) +
  geom_point(alpha = 0.4, size = 1.5) +
  geom_smooth(method = "lm", se = FALSE, linewidth = 1.2) +
  labs(
    title = "The Confounding Effect: Time vs. Energy",
    subtitle = "Larger dimensions inherently take more time, tangling the variables.",
    x = "Execution Time (Seconds)",
    y = "Energy Consumption (PKG)",
    color = "Problem\nDimension"
  ) +
  theme_minimal() +
  theme(legend.position = "right")

Severe correlation between execution time and problem dimension.

Because dimension and the rest of the algorithm variables fundamentally drive seconds, a single-stage model cannot determine if an energy spike is caused by the problem size or simply because the clock ticked longer.

2. Stage 1: The Time Model (Extracting the Context)

To break the co-linearity, we first model the algorithm’s execution time based exclusively on its configuration parameters and variables.

Code 
# Stage 1: How LONG SHOULD the algorithm take based on its parameters?
multi_run_time_model <- glm(seconds ~ dimension * population_size * max_gens, data = all_data)

# Extract the residuals: The difference between expected time and actual time
all_data$residual_seconds <- residuals(multi_run_time_model)

Understanding the Residuals

By extracting the residual_seconds, we have quantified the “Environmental Friction”: * Residual = 0: The operating system executed the algorithm exactly in the expected time. * Residual > 0: The server was “noisy” or busy, causing the execution to drag on longer than the pure parameters required. * Residual < 0: The server executed the task unusually fast (e.g., optimal cache hits).

Let’s visualize the shape of this environmental noise across our 5 experimental blocks.

Code 
ggplot(all_data, aes(x = residual_seconds)) +
  geom_density(fill = "steelblue", alpha = 0.5, color = "darkblue", linewidth = 1) +
  geom_vline(xintercept = 0, color = "red", linetype = "dashed", linewidth = 1) +
  annotate("text", x = 0.5, y = 0.1, label = "Slower than expected\n(OS Friction)", color = "red", hjust = 0) +
  annotate("text", x = -0.5, y = 0.1, label = "Faster than expected", color = "darkgreen", hjust = 1) +
  labs(
    title = "Isolated Environmental Noise (Residual Seconds)",
    subtitle = "This metric is mathematically entirely orthogonal to our parameters.",
    x = "Residual Time (Seconds)",
    y = "Density"
  ) +
  theme_minimal()

Distribution of Operating System noise isolated from algorithmic complexity.

What we are seeing here is a two peak distribution. We might need to understand it better.

Code 
#! label: two-regime-distribution
#! fig-cap: understanding the two "time" regimes

# 1. Compute the density explicitly
dens <- density(all_data$residual_seconds)
dens_df <- data.frame(x = dens$x, y = dens$y)

# 2. Find the local maxima (peaks)
# Mathematically, a peak occurs where the slope changes from positive to negative
peaks_indices <- which(diff(sign(diff(dens_df$y))) == -2) + 1
peaks <- dens_df[peaks_indices, ]

# Sort by density (y-axis) to grab the two most prominent peaks, 
# then re-sort them by the x-axis (left to right)
top_peaks <- peaks[order(-peaks$y), ][1:2, ]
top_peaks <- top_peaks[order(top_peaks$x), ]

# 3. Plot the data with the detected regimes
ggplot(all_data, aes(x = residual_seconds)) +
  geom_density(fill = "steelblue", alpha = 0.5, color = "darkblue", linewidth = 1) +
  
  # Peak 1: High-Power / Fast State (Left)
  geom_vline(xintercept = top_peaks$x[1], color = "darkgreen", linetype = "dashed", linewidth = 1) +
  annotate("text", x = top_peaks$x[1], y = top_peaks$y[1] + 0.05, 
           label = paste("Turbo Regime\nPeak:", round(top_peaks$x[1], 2), "s"), 
           color = "darkgreen", hjust = 0.5) +
  
  # Peak 2: Low-Power / Throttled State (Right)
  geom_vline(xintercept = top_peaks$x[2], color = "darkred", linetype = "dashed", linewidth = 1) +
  annotate("text", x = top_peaks$x[2], y = top_peaks$y[2] + 0.05, 
           label = paste("Throttled Regime\nPeak:", round(top_peaks$x[2], 2), "s"), 
           color = "darkred", hjust = 0.5) +
  
  labs(
    title = "Observing operating context",
    subtitle = "Bimodal distribution reveals two distinct processor performance regimes.",
    x = "Residual Time (Seconds)",
    y = "Density"
  ) +
  theme_minimal()

By creating a model and working with the residuals of that model, we begin to understand how the operating context is influencing the energy we are measuring: we have two different peaks, one where an experiment is taking longer than expected: the system has lowered frequency, and then the CPU can extract less power, making it take longer. This is less frequent than the opposite: taking less time than expected because the experiment has been assigned to a higher-frequency core or simply because the CPU is being maintained at a higher frequency.

3. Stage 2: The Energy Model (The Clean Footprint)

With the OS noise mathematically isolated (residual_seconds), we can now build our definitive energy model. Because the residuals are completely uncorrelated with the dimension or population_size variables, multicollinearity is reduced to zero.

Code 
# Stage 2: Modeling Energy using isolated noise + pure algorithmic parameters
multi_run_PKG_model <- glm(PKG ~ residual_seconds + I(residual_seconds^2) + 
                                 dimension * population_size * max_gens, 
                           data = all_data)

3.1 Visualizing the Hardware Penalty

First, let’s look at what the residual_seconds terms are doing. They act as a mathematical sponge, absorbing the non-linear energy penalty of a busy operating system.

Code 
# Generate predictions varying ONLY the noise, keeping base parameters constant
pred_noise <- predictions(
  multi_run_PKG_model, 
  newdata = datagrid(residual_seconds = seq(min(all_data$residual_seconds), 
                                            max(all_data$residual_seconds), 
                                            length.out = 100))
)

ggplot(pred_noise, aes(x = residual_seconds, y = estimate)) +
  geom_point(data = all_data, aes(x = residual_seconds, y = PKG), 
             color = "grey", alpha = 0.2, size = 1) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), fill = "firebrick", alpha = 0.3) +
  geom_line(color = "darkred", linewidth = 1.2) +
  geom_vline(xintercept = 0, linetype = "dotted", color = "black", linewidth = 1) +
  labs(
    title = "The Energy Cost of Operating System Friction",
    subtitle = "The quadratic curve absorbs the hardware penalty, cleaning the parameters.",
    x = "Environmental Noise (Residual Seconds)",
    y = "Predicted Energy (PKG)"
  ) +
  theme_minimal()

The isolated effect of environmental friction on total energy consumption.

Here you can also observe how the two regimes play out: high-frequency, the majority, would be below the line, lower-frequency above the line. The paradox here is that since we can only measure energy-consumed-while-we-run, since these take longer they consume more energy.

3.2 Visualizing the Pure Algorithmic Truth

Finally, because the Stage 2 model has absorbed the hardware noise, the coefficients for our parameters (dimension, population_size, max_gens) now represent the pure algorithmic truth.

We can ask the model: “What would the energy consumption be if the operating system were perfectly silent (residual_seconds = 0)?”

Code 
# Predict energy across different dimensions and populations, FORCING noise to 0
pred_algo <- predictions(
  multi_run_PKG_model,
  newdata = datagrid(
    dimension = unique(all_data$dimension),
    population_size = unique(all_data$population_size),
    max_gens = unique(all_data$max_gens)[1], # Keep max_gens constant for 2D visualization
    residual_seconds = 0                     # The crucial step: silencing the hardware noise
  )
)

ggplot(pred_algo, aes(x = dimension, y = estimate, color = population_size, group = population_size)) +
  geom_point(size = 4) +
  geom_line(linewidth = 1.2) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.2, linewidth = 0.8) +
  labs(
    title = "Unmasked Algorithmic Energy Footprint",
    subtitle = "Modeled under 'perfect' OS conditions (Residual = 0) revealing the true interaction.",
    x = "Problem Dimension",
    y = "Predicted Pure Energy Cost (PKG)",
    color = "Population\nSize"
  ) +
  theme_bw() +
  theme(legend.position = "bottom")

The pure, unconfounded energy cost of evolutionary parameters.

4. Conclusion

By utilizing this two-stage methodology alongside 5 blocks of experimental replications: 1. We avoided the multicollinearity trap that invalidates standard single-stage energy models. 2. We isolated the macro-noise of the hardware environment into an orthogonal vector (residual_seconds). 3. We secured narrow confidence intervals for complex 3-way interactions, allowing us to definitively prove how algorithm geometry and search density combine to consume physical server energy.

References

  • Gelman, A., & Hill, J. (2006). Data Analysis Using Regression and Multilevel/Hierarchical Models (1st ed.). Cambridge University Press. DOI: 10.1017/CBO9780511790942
  • McCullagh, P., & Nelder, J. A. (1989). Generalized Linear Models (2nd ed.). Chapman and Hall/CRC. DOI: 10.1201/9780203753736