基于Rcpp实现hist计数功能的性能优化求助

Hey there! Let's fix that performance bottleneck in your Rcpp histogram count function. Since your current get_freq is slowing down with larger x vectors (3000+ elements), here are two optimized approaches that’ll give you a big speedup while matching R’s hist(x, breaks, plot=FALSE)$counts behavior:

1. Sorted + Binary Search (Best for Uneven Breaks) #

If your breaks aren’t evenly spaced, sorting the input vector once and then using binary search to count elements per bin is way more efficient than checking every element against every break. This cuts the time complexity from O(n*m) (n = length of x, m = number of breaks) to O(n log n), which makes a huge difference as n grows.

Here’s the code:

#include <Rcpp.h>
#include <algorithm>
using namespace Rcpp;

// [[Rcpp::export]]
IntegerVector fast_hist_counts(NumericVector x, NumericVector breaks) {
  // Clone and sort the input vector (preserves original x)
  NumericVector x_sorted = clone(x);
  std::sort(x_sorted.begin(), x_sorted.end());
  
  int bin_count = breaks.size() - 1;
  IntegerVector counts(bin_count, 0);
  
  auto current_start = x_sorted.begin();
  for (int i = 0; i < bin_count; ++i) {
    // Find the first element greater than the upper break of the current bin
    auto current_end = std::upper_bound(current_start, x_sorted.end(), breaks[i+1]);
    // The number of elements in the bin is the distance between iterators
    counts[i] = current_end - current_start;
    // Move the start to the end of the current bin for the next iteration
    current_start = current_end;
  }
  
  return counts;
}

This matches R’s default hist behavior (right-closed bins, right=TRUE, and includes the lowest boundary with include.lowest=TRUE).

2. Direct Bin Index Calculation (Best for Evenly Spaced Breaks) #

If your breaks are evenly spaced, you can skip sorting entirely and calculate each element’s bin index directly. This is O(n) time—blazingly fast even for very large x.

Here’s the code, with checks to ensure breaks are evenly spaced and match R’s boundary handling:

#include <Rcpp.h>
#include <cmath>
using namespace Rcpp;

// [[Rcpp::export]]
IntegerVector fast_hist_equal_breaks(NumericVector x, NumericVector breaks) {
  int bin_count = breaks.size() - 1;
  IntegerVector counts(bin_count, 0);
  
  double min_break = breaks[0];
  double bin_width = breaks[1] - breaks[0];
  
  // Verify breaks are evenly spaced (optional but prevents bugs)
  bool is_equal_spaced = true;
  for (int i = 1; i < bin_count; ++i) {
    if (std::abs(breaks[i+1] - breaks[i] - bin_width) > 1e-8) {
      is_equal_spaced = false;
      break;
    }
  }
  if (!is_equal_spaced) stop("Breaks must be evenly spaced for this function!");
  
  for (double val : x) {
    if (val >= breaks.back()) {
      // Place values equal to the highest break in the last bin (matches R)
      counts[bin_count - 1]++;
    } else if (val < min_break) {
      // Place values below the lowest break in the first bin (matches R)
      counts[0]++;
    } else {
      // Calculate bin index directly
      int idx = static_cast<int>(std::floor((val - min_break) / bin_width));
      // Ensure index stays within valid range
      idx = std::max(0, std::min(idx, bin_count - 1));
      counts[idx]++;
    }
  }
  
  return counts;
}

How to Test Performance #

Use the microbenchmark package to compare your original get_freq with these new functions against R’s built-in hist:

library(microbenchmark)
set.seed(123)
x <- rnorm(5000) # Larger test vector to see the difference
breaks <- seq(-3, 3, by=0.25)

microbenchmark(
  original = get_freq(x, breaks),
  sorted_bin = fast_hist_counts(x, breaks),
  equal_bin = fast_hist_equal_breaks(x, breaks),
  r_hist = hist(x, breaks=breaks, plot=FALSE)$counts,
  times = 100
)

You’ll see that both new Rcpp functions outperform your original get_freq—and even beat R’s built-in hist in most cases!

内容的提问来源于stack exchange，提问作者Adrian