The Strider Package

Replacing index loops with standard library algorithms

As a prelude to using the strided_iterator class, I start out with an example of vector convolution because it is one of the original examples in the Rcpp documentation, which I have copied verbatim here. I have placed some definitions into a header file whose contents can be found in the RMarkdown source file used to generate this output. Otherwise, all compiler settings were left at their default values.

// [[Rcpp::depends(BH)]]
// [[Rcpp::depends(strider)]]
// [[Rcpp::plugins(cpp11)]]

#include <vignette.h>

// [[Rcpp::export]]
NumericVector
convolve_cpp(const NumericVector& a,
             const NumericVector& b) {
  // Declare loop counters, and vector sizes
  int i, j,
      na = a.size(), nb = b.size(),
      nab = na + nb - 1;
  // Create vector filled with 0
  NumericVector ab(nab);
  // Crux of the algorithm
  for(i = 0; i < na; i++) {
    for(j = 0; j < nb; j++) {
          ab[i + j] += a[i] * b[j];
    }
  }
  // Return result
  return ab;
}

The next code block shows my translation using the standard library transform algorithm. This version eschews integer indexes for incrementing iterators along the input and output vectors.

// [[Rcpp::depends(BH)]]
// [[Rcpp::depends(strider)]]
// [[Rcpp::plugins(cpp11)]]

#include <vignette.h>

// [[Rcpp::export]]
NumericVector
stl_algo_convolve(const NumericVector& a,
                  const NumericVector& b) {
  // Declare vector sizes
  int
    na = a.size(), nb = b.size(), 
    nab = na + nb - 1;
  // Create vector filled with 0
  NumericVector ab(nab);
  // Crux of the algorithm
  transform(begin(a), end(a),
            begin(ab), begin(ab),
            [&b](const double t, double& u) {
              transform(begin(b), end(b), &u, &u,
                        [t](const double v, const double w) {
                          return w + t * v; });
              return u; });
  // Return result
  return ab;
}

While the code is certainly more busy, replacing explicit loops with standard library algorithms is considered by some to be a best-practice that is less susceptible to indexing errors and buffer overruns.

C++11 introduced range-based for-loops that many of the advantages of the standard library algorithms in terms of avoiding indexing errors and are simpler and somewhat more natural to reason about.

// [[Rcpp::depends(BH)]]
// [[Rcpp::depends(strider)]]
// [[Rcpp::plugins(cpp11)]]

#include <vignette.h>

// [[Rcpp::export]]
NumericVector
range_loop_convolve(const NumericVector& a,
                    const NumericVector& b) {
  // Declare vector sizes
  int
    na = a.size(), nb = b.size(), 
    nab = na + nb - 1;
  // Create vector filled with 0
  NumericVector ab(nab);
  // Crux of the algorithm
  auto t = &ab[0];
  for (const auto u : a)
  {
    auto v = t++;
    for (const auto w : b)
      *v++ += u * w;
  }
  // Return result
  return ab;
}

The main difference here is we have to manually increment the output iterator. That could be solved using a zip iterator adapter. (Convenient tuple syntax for this case is available in C++17.) Nonetheless, this example demonstrates how expressive range-based for-loops can be while avoiding explicit indices.

There is no particular performance advantage to avoiding indices, so the choice is largely a matter of taste. Modern compilers will produce similar if not the exact same instructions regardless of the type of loop. Benchmarking comparisons can however point out problems with a particular implementation. It is useful to recall that benchmark differences within \(\pm 10\%\) are likely meaningless and results will vary in different environments. Relative numbers below are therefore printed using digits = 1 to demphasize small differences. With two vectors of length \(10^4\), I get the following performances for the different loop implementations.

Striding over matrix margins

As a minimal example, consider computing the column sums of a matrix. I provide three example implementations using indices, make_strided and make_strided_range.

// [[Rcpp::depends(BH)]]
// [[Rcpp::depends(strider)]]
// [[Rcpp::plugins(cpp11)]]

#include <vignette.h>

// [[Rcpp::export]]
NumericVector indx_col_sum(const NumericMatrix& x)
{
  auto
    nr = x.nrow(),
    nc = x.ncol();
  NumericVector res(nc, 0.0);
  for (int c = 0; c != nc; ++c)
    for (int r = 0; r != nr; ++r)
      res[c] += x(r, c);
  return res;
}

// [[Rcpp::export]]
NumericVector range_col_sum(const NumericMatrix& x)
{
  auto
    nr = x.nrow(),
    nc = x.ncol();
  NumericVector res(nc);
  auto t = res.begin();
  for (const auto& u : make_strided_range(begin(x), nr, nc))
  {
    for (auto v = &u; v != &u + nr; ++v) *t += *v;
    ++t;
  }
  return res;
}

// [[Rcpp::export]]
NumericVector strided_col_sum(const NumericMatrix& x)
{
  auto nr = x.nrow();
  NumericVector res(x.ncol());
  transform(make_strided(begin(x), nr), make_strided(end(x)), begin(res),
            [nr](const double& v){ return accumulate(&v, &v + nr, 0.0); });
  return res;
}

In this case, owing to the column-major memory layout of R matrices, the STL transform algorithm combined with make_strided is the most compact in terms of code. Furthermore, the strided iterator is incremented in the outer loop meaning the scanning of the matrix is in the optimal order. The performance of these three implementations are equivalent, as is R’s built in colSumms implemented in C.

Row sums

Summing over rows is more interesting because R matrices are stored in column-major format. That means that looping over a row will skip from column to column by-passing all the rows in between. This non-local memory access is a drag on performance. Four implementations of the row-summing algorithm are given below.

// [[Rcpp::depends(BH)]]
// [[Rcpp::depends(strider)]]
// [[Rcpp::plugins(cpp11)]]

#include <vignette.h>

// [[Rcpp::export]]
NumericVector indx_row_sum(const NumericMatrix& x)
{
  auto
    nr = x.nrow(),
    nc = x.ncol();
  NumericVector res(nr, 0.0);
  for (int c = 0; c != nc; ++c)
    for (int r = 0; r != nr; ++r)
      res[r] += x(r, c);
  return res;
}

// [[Rcpp::export]]
NumericVector range_row_sum(const NumericMatrix& x)
{
  auto
    nr = x.nrow(),
    nc = x.ncol();
  NumericVector res(nc);
  for (const auto& u : make_strided_range(begin(x), nr, nc))
  {
    auto t = begin(res);
    for (auto v = &u; v != &u + nr; ++v) *t++ += *v;
  }
  return res;
}

// [[Rcpp::export]]
NumericVector strided_row_sum(const NumericMatrix& x)
{
  auto
    nr = x.nrow(),
    nc = x.ncol();
  NumericVector res(nr);
  transform(begin(x), begin(x) + nr, begin(res), [nr, nc](const double& v) {
    return accumulate(make_strided(&v, nr), make_strided(&v, nr, nc), 0.0); });
  return res;
}

// [[Rcpp::export]]
NumericVector strided_row_sum2(const NumericMatrix& x)
{
  auto nr = x.nrow();
  NumericVector res(nr, 0.0);
  for_each(make_strided(begin(x), nr), make_strided(end(x)), [&res, nr](const double& y) {
    transform(&y, &y + nr, begin(res), begin(res), plus<double>()); });
  return res;
}

Again, the transform with accumulate implementation is compact, but here the strided iterator is in the inner loop causing poor memory locality. The alternative implementation uses for_each with make_strided to skip to the top of each column and then uses transform to accumulate the values. The final range-based version is equivalent to the for_each with transform version.

The benchmark results show that indeed placing the strided iterator in the inner loop gives poor results. The index-based, range-based and for_each-transform version are all equivalent in efficiency. Interestingly, R’s built in rowSums is several times slower, most likely indicating that the column index is incremented in the inner loop.

2D convolution

A much better demonstration of the strided approach is to compute the 2-dimensional convolution. This requires four nested loops and much more expensive indexing calculations. Furthermore, because C++ does not provide a bivariate index operators, accessing 2D data structures requires defining a new method and the details of this method will strongly determine overhead. In the strided approach, we simply walk the output matrix by iteration, so no indexing is required.

// [[Rcpp::depends(BH)]]
// [[Rcpp::depends(strider)]]
// [[Rcpp::plugins(cpp11)]]

#include <vignette.h>

// [[Rcpp::export]]
NumericMatrix
convolve2_cpp(const NumericMatrix& a,
              const NumericMatrix& b)
{
  int
    nra = a.nrow(), nca = a.ncol(),
    nrb = b.nrow(), ncb = b.ncol();
  NumericMatrix ab(nra + nrb - 1, nca + ncb - 1);
  for(int i = 0; i != nca; ++i)
    for (int j = 0; j != ncb; ++j)
      for (int k = 0; k != nra; ++k)
        for (int l = 0; l != nrb; ++l)
          ab(k + l, i + j) += a(k, i) * b(l, j);
  return ab;
}

// [[Rcpp::export]]
NumericMatrix
range_loop_convolve2(const NumericMatrix& a,
                     const NumericMatrix& b)
{
  int
    nra = a.nrow(), nca = a.ncol(),
    nrb = b.nrow(), ncb = b.ncol(),
    nrab = nra + nrb - 1,
    ncab = nca + ncb - 1;
  NumericMatrix ab(nrab, ncab);
  auto iter1 = make_strided(begin(ab), nrab);
  for (const auto& t : make_strided_range(begin(a), nra, nca)) {
    auto iter2 = make_strided(&*iter1++, nrab);
    for (const auto& u : make_strided_range(begin(b), nrb, ncb)) {
      auto iter3 = &*iter2++;
      for (const auto& v : make_strided_range(&t, 1, nra)) {
        auto iter4 = iter3++;
        for (const auto& w : make_strided_range(&u, 1, nrb)) {
          *iter4++ += v * w; }}}}
  return ab;
}

// [[Rcpp::export]]
NumericMatrix
stl_algo_convolve2(const NumericMatrix& a,
                   const NumericMatrix& b)
{
  int
    nra = a.nrow(), nca = a.ncol(),
    nrb = b.nrow(), ncb = b.ncol(),
    nrab = nra + nrb - 1,
    ncab = nca + ncb - 1;
  NumericMatrix ab(nrab, ncab);
  transform(make_strided(begin(a), nra), make_strided(end(a)),
            make_strided(begin(ab), nrab), make_strided(begin(ab), nrab),
            [&](const double& t, double& u) {
    transform(make_strided(begin(b), nrb), make_strided(end(b)),
              make_strided(&u, nrab), make_strided(&u, nrab),
              [&](const double& v, double& w) {
      transform(&t, &t + nra, &w, &w,
                [&](const double x, double& y) {
        transform(&v, &v + nrb, &y, &y,
                  [&](const double z, const double zz) {
          return zz + x * z; });
        return y; });
      return w; });
    return u; });
  return ab;
}

As a quick test, I check whether the convolution is invariant to the delta function as required.

a = matrix(c(1, 2, 3,
             4, 5, 6), 2, 3, byrow = TRUE)
b = matrix(c(0, 0, 0,
             0, 0, 0,
             0, 1, 0,
             0, 0, 0), 4, 3, byrow = TRUE)
stl_algo_convolve2(a, b)

##      [,1] [,2] [,3] [,4] [,5]
## [1,]    0    0    0    0    0
## [2,]    0    0    0    0    0
## [3,]    0    1    2    3    0
## [4,]    0    4    5    6    0
## [5,]    0    0    0    0    0

The output is correct. Here are benchmark results for \(50\times50\) random matrices.

The transform-based implementation appears to have a slight advantage, however these results may change on different systems. It is possible the index-offset calculations are making the indexed version slower. The importance of offset calculations should increase with the number of dimensions. It should be noted that these are not the fasted algorithms for performing convolution. If the kernel is separable, each dimension can be convolved individually. For large matrices, Fourier methods are faster still.