[Rd] Why is the diag function so slow (for extraction)?

[franknarf]

(I  asked this question on StackOverflow
<http://stackoverflow.com/q/30035939/1191259>   a short time ago; sorry if
you're seeing it again. Feel free to answer there as well if you like. The
code formatting and such on that site can be nice.)

I benchmarked matrix and vector subsetting to extract the diagonal of a
square matrix against the diag() function, and the latter lost by a wide
margin:

nc  <- 1e4
set.seed(1)
m <- matrix(sample(letters,nc^2,replace=TRUE), ncol = nc)

microbenchmark(
  diag = diag(m),
  cond = m[row(m)==col(m)],
  vec  = m[(1:nc-1L)*nc+1:nc],
  mat  = m[cbind(1:nc,1:nc)],
times=10)

# results
Unit: microseconds
 expr         min          lq         mean       median          uq        
max neval
 diag  604343.469  629819.260  710371.3320  706842.3890  793144.019 
837115.504    10
 cond 3862039.512 3985784.025 4175724.0390 4186317.5260 4312493.742
4617117.706    10
  vec     317.088     329.017     432.9099     350.1005     629.460    
651.376    10
  mat     272.147     292.953     441.7045     345.9400     637.506    
706.860    10

Looking in the diag() function, I suspect the fault is mostly due to its
core operation of c(m)[v], where v is the vector in the "vec" benchmark
above. This code seems to confirm it:

v <- (1:nc-1L)*nc+1:nc
microbenchmark(diaglike=c(m)[v],vec=m[v])
# Unit: microseconds
#      expr        min          lq        mean     median          uq       
max neval
#  diaglike 579224.436 664853.7450 720372.8105 712649.706 767281.5070
931976.707   100
#       vec    334.843    339.8365    568.7808    646.799    663.5825  
1445.067   100

But I'm still wondering why diag() uses c()...? With it being so slow, I'd
be inclined to write a qdiag() without the c() and just use that the next
time I need matrix algebra. Any insight would be appreciated; thanks!




--
View this message in context: http://r.789695.n4.nabble.com/Why-is-the-diag-function-so-slow-for-extraction-tp4706780.html
Sent from the R devel mailing list archive at Nabble.com.