1. Understanding the Algorithms
1.1 Exact vs Greedy: When to Use Each
corrselect offers two algorithmic approaches for
corrPrune():
Exact Mode (Graph-Theoretic)
Algorithm: Eppstein–Löffler–Strash (ELS) or Bron–Kerbosch Complexity: O(2^p) - exponential in number of predictors Guarantee: Finds the largest maximal independent set
data(mtcars)
# Exact mode: guaranteed optimal
exact_result <- corrPrune(mtcars, threshold = 0.7, mode = "exact")
cat("Exact mode kept:", ncol(exact_result), "variables\n")
#> Exact mode kept: 5 variablesUse exact mode when:
- p ≤ 20 (feasible runtime)
- You need guaranteed optimal solution
- Reproducibility is critical
- You’re writing a paper (justify optimality)
Greedy Mode (Heuristic)
Algorithm: Deterministic iterative removal Complexity: O(p² × k) where k = iterations Guarantee: Near-optimal, deterministic
# Greedy mode: fast approximation
greedy_result <- corrPrune(mtcars, threshold = 0.7, mode = "greedy")
cat("Greedy mode kept:", ncol(greedy_result), "variables\n")
#> Greedy mode kept: 5 variablesUse greedy mode when:
- p > 20 (exact becomes slow)
- Speed is priority
- Near-optimal is acceptable
- High-dimensional data (p > 100)
1.2 Complexity Analysis with Benchmarks
Let’s measure runtime scaling:
# Generate datasets with increasing p
library(microbenchmark)
benchmark_corrPrune <- function(p_values) {
results <- data.frame(
p = integer(),
exact_time_ms = numeric(),
greedy_time_ms = numeric()
)
for (p in p_values) {
# Generate correlated data
set.seed(123)
cor_mat <- 0.5^abs(outer(1:p, 1:p, "-"))
data <- as.data.frame(MASS::mvrnorm(n = 100, mu = rep(0, p), Sigma = cor_mat))
# Exact mode (skip if p too large)
exact_time_ms <- if (p <= 500) {
median(microbenchmark(
corrPrune(data, threshold = 0.7, mode = "exact"),
times = 3, # Few iterations for speed
unit = "ms"
)$time) / 1e6 # Convert nanoseconds to milliseconds
} else {
NA
}
# Greedy mode
greedy_time_ms <- median(microbenchmark(
corrPrune(data, threshold = 0.7, mode = "greedy"),
times = 3, # Few iterations for speed
unit = "ms"
)$time) / 1e6 # Convert nanoseconds to milliseconds
results <- rbind(results, data.frame(
p = p,
exact_time_ms = round(exact_time_ms, 1),
greedy_time_ms = round(greedy_time_ms, 1)
))
}
results
}
# Benchmark (extended range to show comprehensive scaling behavior)
p_values <- c(10, 20, 50, 100, 200, 300, 500, 1000)
benchmark <- benchmark_corrPrune(p_values)
print(benchmark)
#> p exact_time_ms greedy_time_ms
#> 1 10 0.7 0.2
#> 2 20 0.8 0.4
#> 3 50 1.0 0.7
#> 4 100 1.9 1.4
#> 5 200 5.1 3.1
#> 6 300 12.1 5.7
#> 7 500 32.0 12.5
#> 8 1000 NA 43.2# Visualize scaling
# Separate data for exact mode (only where available) and greedy mode (all points)
exact_valid <- !is.na(benchmark$exact_time_ms) & benchmark$exact_time_ms > 0
greedy_valid <- !is.na(benchmark$greedy_time_ms) & benchmark$greedy_time_ms > 0
# Replace any zeros with small positive value for log scale
exact_times <- benchmark$exact_time_ms
greedy_times <- benchmark$greedy_time_ms
exact_times[exact_times <= 0 & !is.na(exact_times)] <- 0.01
greedy_times[greedy_times <= 0 & !is.na(greedy_times)] <- 0.01
# Determine y-axis limits from valid data
all_valid_times <- c(exact_times[exact_valid], greedy_times[greedy_valid])
ylim <- c(min(all_valid_times) * 0.5, max(all_valid_times) * 2)
# Plot greedy mode (all points)
plot(benchmark$p[greedy_valid], greedy_times[greedy_valid],
type = "b", col = rgb(0.2, 0.5, 0.8, 1), pch = 19, lwd = 2,
xlab = "Number of Predictors (p)",
ylab = "Time (milliseconds, log scale)",
main = "Exact vs Greedy Scaling",
ylim = ylim,
xlim = range(benchmark$p),
log = "y")
# Add exact mode (only where available)
lines(benchmark$p[exact_valid], exact_times[exact_valid],
type = "b", col = rgb(0.8, 0.2, 0.2, 1), pch = 19, lwd = 2)
# Mark exact mode limit
abline(v = 500, lty = 2, col = "gray50")
text(500, ylim[2] * 0.5, "Exact mode limit", pos = 4, col = "gray30")
legend("topleft",
legend = c("Exact", "Greedy"),
col = c(rgb(0.8, 0.2, 0.2, 1), rgb(0.2, 0.5, 0.8, 1)),
pch = 19,
lwd = 2,
bty = "o",
bg = "white")
Key insight: The scaling behavior depends heavily on correlation structure. For this moderately correlated dataset (ρ = 0.5^|i-j|), exact mode remains practical up to p ≈ 500, while greedy mode scales efficiently beyond p = 1000. The exact mode provides guaranteed optimality at the cost of computational complexity, while greedy mode offers substantial speed advantages for large p with near-optimal results in most practical scenarios.
1.3 Deterministic Tie-Breaking
When multiple variables have identical correlation profiles, corrselect uses deterministic tie-breaking:
# Create data with identical correlations
set.seed(123)
x1 <- rnorm(100)
x2 <- x1 + rnorm(100, sd = 0.1) # Almost identical to x1
x3 <- x1 + rnorm(100, sd = 0.1) # Also almost identical to x1
x4 <- rnorm(100) # Independent
data_ties <- data.frame(x1, x2, x3, x4)
# Run multiple times - always same result
result1 <- corrPrune(data_ties, threshold = 0.95)
result2 <- corrPrune(data_ties, threshold = 0.95)
cat("Run 1 selected:", names(result1), "\n")
#> Run 1 selected: x2 x4
cat("Run 2 selected:", names(result2), "\n")
#> Run 2 selected: x2 x4
cat("Identical:", identical(names(result1), names(result2)), "\n")
#> Identical: TRUETie-breaking rules: 1. Prefer variables with lower mean absolute correlation with others 2. If still tied, prefer lexicographically first (alphabetical)
This ensures reproducibility across runs, machines, and R versions.
