index.Rmd (4886B)
1 --- 2 title: "Indexing matrix elements in R" 3 author: "Eamon Caddigan" 4 date: "2015-10-22" 5 output: html_document 6 layout: post 7 summary: Why I prefer zero-based numbering. 8 categories: R programming 9 --- 10 11 ```{r global_options, include=FALSE} 12 knitr::opts_chunk$set(cache=TRUE) 13 ``` 14 15 I came to science with a background in engineering, but most of my scientist friends didn’t. I often have a hard time articulating why I’m so annoyed by one-based indexing--which R and MATLAB use, but most other programming languages don’t. Here’s a recent example that might help. 16 17 For a simulation I’m running, I use the values in several of the columns of a data frame as indexes into separate vectors. Here’s an example of using indexes in one vector (instead of a column from a `data.frame`) to access the elements in another: 18 19 ```{r} 20 valueVector <- c(2, 5, 8, 3, 0) 21 indexVector <- c(2, 3, 3, 1, 4, 5, 5) 22 valueVector[indexVector] 23 ``` 24 25 R veterans will point out that you can use factors for this, and that's definitely true here. However, when the values in the small vector are changing often but the indices are relatively stable, I prefer this approach. Either strategy works. 26 27 Unfortunately, things aren't so easy when the data is in a matrix (a 2D vector) and you want to access its elements using two index vectors (i.e., one indexing the matrix’s rows, and the second indexing its columns). R’s default behavior might not be what you expect: 28 29 ```{r} 30 valueMatrix <- matrix(LETTERS[1:15], ncol = 3) 31 valueMatrix 32 33 rowIndices <- c(4, 2, 2, 4, 4, 4) 34 colIndices <- c(2, 3, 3, 2, 3, 2) 35 valueMatrix[rowIndices, colIndices] 36 ``` 37 38 Instead of returning the six values associated with the six row/column pairs, it returns a 6 × 6 matrix with the elements of interest along the diagonal. When the number of elements is this small, it’s easy to wrap this in a call to `diag()`; that approach isn’t appropriate for long index vectors, since memory requirements are squared. 39 40 Readers who’ve done data manipulation in C are probably familiar with pointer arithmetic. Zero-based numbering makes it easy to combine row and column indexes and make a one-dimensional array behave like a matrix: 41 42 ``` 43 /* When the code has to deal with matrices of varying size, you can’t allocate 44 an array right away. */ 45 double *valueMatrix; 46 47 /* Stuff happened, and now you know how big your matrix will be (nrow x ncol), 48 so you allocate memory on the “heap” to store the data. */ 49 valueMatrix = (double *)malloc(sizeof(double) * nrow * ncol); 50 51 /* This is how you’d access specific elements. */ 52 oldValueAtRowCol = valueMatrix[row + nrow*col]; 53 valueMatrix[row + nrow*col] = newValueAtRowCol; 54 55 /* Can’t forget this when you’re done; I don’t miss C. */ 56 free(valueMatrix); 57 valueMatrix = 0; 58 59 /* NB: All programs work this way, but most scripting languages take care of 60 this stuff for you. Progress! */ 61 ``` 62 63 To its credit, R makes this easy too; everything is stored as a one-dimensional vector behind the scenes, and can be accessed as such. R’s use of one-based indexing just makes it look a bit more awkward: 64 65 ```{r} 66 valueMatrix[rowIndices + nrow(valueMatrix) * (colIndices - 1)] 67 ``` 68 69 I hope this is useful for anybody who wants to access matrix elements using index vectors. More importantly, maybe some scientists-turned-programmers can gain a bit of insight into why zero-based numbering make sense to those of us who cut our teeth on C (or languages like it). 70 71 *** 72 73 ### Update 74 75 After I posted this, [Andrie de Vries](http://rfordummies.com/) sent [a tweet](https://twitter.com/RevoAndrie/status/657187336833388545) sharing an alternate syntax that's more R-like: 76 77 ```{r} 78 valueMatrix[cbind(rowIndices, colIndices)] 79 ``` 80 81 I definitely agree that it looks better, but there is a slight performance hit due to the call to `cbind()`. If you're repeatedly accessing a matrix with the same pair of indices, it might be worth it store the bound pair as a variable and reuse that. Here's the benchmark performance of each of the approaches I've discussed: 82 83 ```{r} 84 library("microbenchmark") 85 86 rowIndices <- sample(1:5, 1e4, replace=TRUE) 87 colIndices <- sample(1:3, 1e4, replace=TRUE) 88 boundIndices <- cbind(rowIndices, colIndices) 89 90 op <- microbenchmark(diag_matrix = diag(valueMatrix[rowIndices, colIndices]), 91 pointer_math = valueMatrix[rowIndices + nrow(valueMatrix) * (colIndices - 1)], 92 array_indexing = valueMatrix[cbind(rowIndices, colIndices)], 93 array_indexing_prebound = valueMatrix[boundIndices], 94 times = 100) 95 print(op) 96 ``` 97 98 Indexing with the pre-bound pair is fastest, using arithmetic on the indexes is a close second, and calling `cbind()` inside the brackets is in third place. Creating the n × n matrix and extracting its diagonal is excessively slow (and uses up a lot of RAM), so don't ever do that.