Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

### Introduction

Besides outstanding support for dense matrices, the Armadillo library also provides a great way to
manipulate sparse matrices in C++. However, the performance characteristics of dealing with sparse
matrices may be surprising if one is only familiar with dense matrices. This is a collection of
observations on getting best performance with sparse matrices in Armadillo.

number of substantial optimisations for sparse matrix operations, in some cases speeding things up
by as much as two orders of magnitude.

### General considerations: sparsity, row vs column access

Perhaps the most important thing to note is that the efficiency of sparse algorithms can depend
strongly on the level of sparsity in the data. If your matrices and vectors are very sparse (most
elements equal to zero), you will often see better performance even if the nominal sizes of those
matrices remain the same. This isn’t specific to C++ or Armadillo; it applies to sparse algorithms
in general, including the code used in the Matrix package for R. By contrast, algorithms for working
with dense matrices usually aren’t affected by sparsity.

Similarly, the pattern of accesssing elements, whether by rows or by columns, can have a significant
impact on performance. This is due to caching, which modern CPUs use to speed up memory access:
accessing elements that are already in the cache is much faster than retrieving them from main
memory. If one iterates or loops over the elements of a matrix in Armadillo, one should try to
iterate over columns first, then rows, to maximise the benefits of caching. This applies to both
dense and sparse matrices. (Technically, at least for dense matrices, whether to iterate over rows
or columns first depends on how the data is stored in memory. Both R and Armadillo store matrices in
column-major order, meaning that
elements in the same column are contiguous in memory. Sparse matrices are more complex but the
advice to iterate by columns is basically the same; see below.)

### Matrix multiplication

We start with a simple concrete example: multiplying two matrices together. In R, this can be done
using the %*% operator which (via the Matrix package) is able to handle any combination of sparse
and dense inputs. However, let us assume we want to do the multiplication in Armadillo: for example
if the inputs are from other C++ functions, or if we want more precise control of the output.

In Armadillo, the * operator multiplies two matrices together, and this also works for any
combination of sparse and dense inputs. However, the speed of the operation can vary tremendously,
depending on which of those inputs is sparse. To see this, let us define a few simple functions:

The outputs of these functions are all the same, but they take different types of inputs: either two
sparse matrices, or a sparse and a dense matrix, or a dense and a sparse matrix (the order
matters). Let us call them on some randomly generated data:

   user  system elapsed
0.230   0.091   0.322

   user  system elapsed
0.407   0.036   0.442

   user  system elapsed
1.081   0.100   1.181

   user  system elapsed
0.826   0.087   0.913

 TRUE


While all the times are within an order of magnitude of each other, multiplying a dense and a sparse
matrix takes about twice as long as multiplying two sparse matrices together, and multiplying a
sparse and dense matrix takes about three times as long. (The sparse-dense multiply is actually one
area where Armadillo 8.500 makes massive gains over previous versions. This operation used to take
much longer due to using an unoptimised multiplication algorithm.)

Let us see if we can help the performance of the mixed-type functions by creating a temporary sparse
copy of the dense input. This forces Armadillo to use the sparse-sparse version of matrix multiply,
which as seen above is much more efficient. For example, here is the result of tweaking the
dense-sparse multiply. Creating the sparse copy does take some extra time and memory, but not enough
to affect the result.

   user  system elapsed
0.401   0.047   0.448

 TRUE


This shows that mixing sparse and dense inputs can have significant effects on the efficiency of
your code. To avoid unexpected slowdowns, consider sticking to either sparse or dense classes for
all your objects. If one decides to mix them, it is worth remembering to test and profile the code.

### Row vs column access

Consider another simple computation: multiply the elements of a matrix by their row number, then sum
them up. (The multiply by row number is to make it not completely trivial.) That is, we want to
obtain:

Armadillo lets us extract individual rows and columns from a matrix, using the .row() and .col()
member functions. We can use row extraction to do this computation:

 1248361320

user  system elapsed
1.708   0.000   1.707


For a large matrix, this takes a not-insignificant amount of time, even on a fast machine. To speed
it up, we will make use of the fact that Armadillo uses the compressed sparse column (CSC) format
for storing sparse matrices. This means that the data for a matrix is stored as three vectors: the
nonzero elements; the row indices of these elements (ordered by column); and a vector of offsets for
the first row index in each column. Since the vector of row indices is ordered by column, and we
have the starting offsets for each column, it turns out that extracting a column slice is actually
very fast. We only need to find the offset for that column, and then pull out the elements and row
indices up to the next column offset. On the other hand, extracting a row is much more work; we have
to search through the indices to find those matching the desired row.

We can put this knowledge to use on our problem. Row slicing a matrix is the same as transposing it
and then slicing by columns, so let us define a new function to return a column slice. (Transposing
the matrix takes only a tiny fraction of a second.)

 1248361320

user  system elapsed
0.766   0.000   0.766


The time taken has come down by quite a substantial margin. This reflects the ease of obtaining
column slices for sparse matrices.

### The R–C++ interface

We can take the previous example further. Each time R calls a C++ function that takes a sparse
matrix as input, it makes a copy of the data. Similarly, when the C++ function returns, its sparse
outputs are copied into R objects. When the function itself is very simple—as it is here—all
this copying and memory shuffling can be a significant proportion of the time taken.

Rather than calling sapply in R to iterate over rows, let us do the same entirely in C++:

 1248361320

user  system elapsed
0.933   0.000   0.935


This is again a large improvement. But what if we do the same with column slicing?

 1248361320

user  system elapsed
0.005   0.000   0.006


Now the time is less than a tenth of a second, which is faster than the original code by roughly
three orders of magnitude. The moral to the story is: rather than constantly switching between C++
and R, you should try to stay in one environment for as long as possible. If your code involves a
loop with a C++ call inside (including functional maps like lapply and friends), consider writing
the loop entirely in C++ and combine the results into a single object to return to R.

(It should be noted that this interface tax is less onerous for built-in Rcpp classes such as
NumericVector or NumericMatrix, which do not require making copies of the data. Sparse data
types are different, and in particular Armadillo’s sparse classes do not provide constructors that
can directly use auxiliary memory.)

## Iterators vs elementwise access

Rather than taking explicit slices of the data, let us try using good old-fashioned loops over the
matrix elements. This is easily coded up in Armadillo, and it turns out to be quite efficient,
relatively speaking. It is not as fast as using column slicing, but much better than row slicing.

 1248361320

user  system elapsed
0.176   0.000   0.176


However, we can still do better. In Armadillo, the iterators for sparse matrix classes iterate only
over the nonzero elements. Let us compare this to our naive double loop:

 1248361320

user  system elapsed
0.002   0.000   0.002


This is the best time achieved so far, to the extent that system.time might have difficulty
capturing it. The timings are now so low that we should use the microbenchmark package to get
accurate measurements:

Unit: milliseconds
expr       min        lq      mean    median        uq       max neval
col   4.78921   4.88444   5.05229   4.99184   5.18450   5.50579    20
elem 172.84830 177.20431 179.87007 179.06447 182.08075 188.11256    20
iter   1.02268   1.05447   1.12611   1.12627   1.16482   1.30800    20


Thus, using iterators represents a greater than three-order-of-magnitude speedup over the original
row-slicing code. 