Vectorize matrix operation in R

Question

I have a R x C matrix filled to the k-th row and empty below this row. What i need to do is to fill the remaining rows. In order to do this, i have a function that takes 2 entire rows as arguments, process these rows and output 2 fresh rows (these outputs will fill the empty rows of the matrix, in batches of 2). I have a fixed matrix containing all 'pairs' of rows to be processed, but my for loop is not helping performance:

# the processRows function:
processRows = function(r1, r2)
{
  # just change a little bit the two rows and return it in a compact way    
  nr1 = r1 * 0.1
  nr2 = -r2 * 0.1

  matrix (c(nr1, nr2), ncol = 2)
 }

# M is the matrix
# nrow(M) and k are even, so nLeft is even

M = matrix(1:48, ncol = 3)
# half to fill (can be more or less, but k is always even)
k = nrow(M)/2

# simulate empty rows to be filled
M[-(1:k), ] = 0

cat('before fill')
print(M)

# number of empty rows to fill
nLeft  = nrow(M) - k
nextRow = k + 1

# each row in idxList represents a 'pair' of rows to be processed 
# any pairwise combination of non-empty rows could happen
# make it reproducible

set.seed(1)
idxList = matrix (sample(1:k, k), ncol = 2, byrow = TRUE)

for ( i in 1 : (nLeft / 2))
{
   row1 = M[idxList[i, 1],]
   row2 = M[idxList[i, 2],]

   # the two columns in 'results' will become 2 rows in M
   results = processRows(row1, row2)

   # fill the matrix
   M[nextRow, ] = results[, 1]
   nextRow = nextRow + 1
   M[nextRow, ] = results[, 2]
   nextRow = nextRow + 1
 }

 cat('after fill')
 print(M)

You'll get much better answers if you make your question [reproducible](http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example). If `processRows()` is really involved or not entirely relevant to the question, simplify that bit and give us something we can work with. — Chase, Jul 09 '12 at 01:55
also look at [this recent question](http://stackoverflow.com/questions/11377677/r-vectorised-conditional-replace/11379749) for some free speed suggestions via `complier` package or something more involved via `Rcpp` package — Chase, Jul 09 '12 at 02:02
It is relatively trivial to make two long vectors of all pairs of rows. The real question is what is `processRows()` doing? The assignment operations you are doing in the loop are relatively trivial (you'll likely get virtually no gain from using the `compiler` package in an instance like this). **IF** `processRows()` can take all the data in one pass, you may get a significant performance improvement, but that is likely the bottleneck (and what you have told us _nothing_ about). You may get a slight gain by `M[nextRow:(nextRow+1), ] <- t(results)` to replace what you do in two assignments — Joshua, Jul 09 '12 at 02:13
Now the code is _reproducible_. I removed processRows() because it's not relevant, it takes two vector as args, perform a fast calculation and returns a matrix with 2 columns, each column will become a row in M...the problem is the for loop. — Fernando, Jul 09 '12 at 03:18
Will the new rows always be computed from old rows only, or could there be new rows computed from preceeding new rows? The latter would cause severe rpoblems for parallelization. — MvG, Jul 09 '12 at 14:09
The list of rows to use is fixed, only from the 'old' matrix...i will change the question description to make it more clearer, tks. — Fernando, Jul 09 '12 at 18:34

score 3 · Answer 1 · answered Jul 09 '12 at 04:29

Okay, here is your code first. We run this so that we have a copy of the "true" matrix, the one we hope to reproduce, faster.

#### Original Code (aka Gold Standard) ####
M = matrix(1:48, ncol = 3)
k = nrow(M)/2
M[-(1:k), ] = 0
nLeft  = nrow(M) - k
nextRow = k + 1
idxList = matrix(1:k, ncol = 2)
for ( i in 1 : (nLeft / 2))
{
    row1 = M[idxList[i, 1],]
    row2 = M[idxList[i, 2],]
    results = matrix(c(2*row1, 3*row2), ncol = 2)
    M[nextRow, ] = results[, 1]
    nextRow = nextRow + 1
    M[nextRow, ] = results[, 2]
    nextRow = nextRow + 1
}

Now here is the vectorized code. The basic idea is if you have 4 rows you are processing. Rather than passing them as vectors one at a time, do it at once. That is:

(1:3) * 2
(1:3) * 2
(1:3) * 2
(1:3) * 2

is the same (but slower) as:

c(1:3, 1:3, 1:3, 1:3) * 2

So first, we will use your same setup code, then create the rows to be processed as two long vectors (where all 4 original rows are just strung together as in my simple example above). Then, we take those results, and transform them into matrices with the appropriate dimensions. The last trick is to assign the results back in in just two steps. You can assign to multiple rows of a matrix at once, so we use seq() to get odd and even numbers so assign the first and second column of the results to, respectively.

#### Vectorized Code (testing) ####
M2 = matrix(1:48, ncol = 3)
k2 = nrow(M2)/2
M2[-(1:k2), ] = 0
nLeft2  = nrow(M2) - k2
nextRow2 = k2 + 1
idxList2 = matrix(1:k2, ncol = 2)

## create two long vectors of all rows to be processed
row12 <- as.vector(t(M2[idxList2[, 1],]))
row22 <- as.vector(t(M2[idxList2[, 2],]))

## get all results
results2 = matrix(c(2*row12, 3*row22), ncol = 2)

## add results back
M2[seq(nextRow2, nextRow2 + nLeft2-1, by = 2), ] <- matrix(results2[,1], nLeft2/2, byrow=TRUE)
M2[seq(nextRow2+1, nextRow2 + nLeft2, by = 2), ] <- matrix(results2[,2], nLeft2/2, byrow=TRUE)

## check that vectorized code matches your examples
all.equal(M, M2)

Which on my machine gives:

> all.equal(M, M2)
[1] TRUE

Wow, thanks! Very tricky, i'm trying to understand it. There's one thing i noticed im my code: im creating the **idxList** matrix using byrow = TRUE...so i'm nothing getting the same results. The vector used to create **idxList** is not ordered, and may contain duplicates...is this a problem?! — Fernando, Jul 09 '12 at 05:30
@Fernando you're welcome. It is late and I am brain dead at the moment, but if you're still unsure about it, comment and tomorrow, I will write a more detailed explanation with more examples to show how/why it works (particularly the matrix -> vector -> matrix bits). — Joshua, Jul 09 '12 at 05:30
Now i'm getting the same results, the problem is solved. I also edited the question, providing more details.Thanks everybody! — Fernando, Jul 09 '12 at 21:11

Vectorize matrix operation in R

1 Answers1