Why is the diag function so slow? [in R 3.2.0 or earlier]

Asked 4/5, 2015 at 17:13 Answered 22/6, 2015 at 19:34

I was looking at the benchmarks in this answer, and wanted to compare them with diag (used in a different answer). Unfortunately, it seems that diag takes ages:

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)

Comments: I tested these with identical. I took "cond" from one of the answers to this homework question. Results are similar with a matrix of integers, 1:26 instead of letters.

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

It is just a matrix-subsetting operation, so I don't know why there's so much overhead. Looking inside the function, I see a few checks and then c(m)[v], where v is the same vector used in the "vec" benchmark. Timing these two...

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

...it seems I have found my culprit. So, the new variation on my question is: Why is there a seemingly unnecessary and very time-consuming c in diag?

Afghanistan answered 4/5, 2015 at 17:13 Comment(18)

Regarding overhead, I was looking at this similar question: https://mcmap.net/q/162076/-why-is-mean-so-slow and I was thinking "Wow, matrix algebra isn't coded in Primitives!" – Afghanistan 4/5, 2015 at 17:14

c() forces the input into a vector. Could be to deal with varying input types by coercing to a vector so as to process all types the same way. Or as a quick and dirty means of checking the input type (a data frame input throws an error because of c). – Adina 4/5, 2015 at 17:42

@AlexA. Thanks. It only reaches this point if is.matrix is true. It's a pretty small function if you want to have a look -- just type diag. (Oh and vanilla data.frames are not matrices, it turns out is.matrix(data.frame(1)) # FALSE.) – Afghanistan 4/5, 2015 at 17:45

Yeah, I was having a look at it. Missed the is.matrix part. And as it turns out that c isn't the reason why data frame input throws an error; c works on data frames. I know that data frames aren't matrices. I think a data frame input makes it all the way to .Internal. – Adina 4/5, 2015 at 17:47

Possibly they use c because it removes all attributes. You could ask on the r-devel mailing list. – Renn 4/5, 2015 at 18:0

@AlexA. This is for a different use of diag. Try .Internal(diag(1, 2, 2)) to see what it does. – Renn 4/5, 2015 at 18:1

I agree with @Roland; it might be better asked to the developers. Then once you have an answer come let us know. ;) – Adina 4/5, 2015 at 18:2

@Roland: Yeah, oops. C code is irrelevant in this case; diagonal extraction is done in R and diagonal matrix construction in C. My bad. – Adina 4/5, 2015 at 18:3

Ok, will do, thanks! Here's a copy of the thread: r.789695.n4.nabble.com/… – Afghanistan 4/5, 2015 at 18:15

BTW: m[seq.int(1,nc^2,nc+1)] is the fastest on my machine. – Nyx 4/5, 2015 at 22:22

Not sure if the following is actually realistic, but assume that there is an uncommon "[" method for a matrix-like object with a "class" attribute: xx = structure(1:9, dim = c(3, 3), class = "mymatrix"); "[.mymatrix" = function(x, i) "no '[' method". Then, maybe, the attribute loss from c could be helpful: v = (1:ncol(xx) - 1L) * ncol(xx) + 1:ncol(xx); xx[v] vs c(xx)[v] – Biracial 5/5, 2015 at 14:4

@Biracial For xx to get that far, it would have to have class=c("matrix","mymatrix") (since there's an if(is.matrix(xx)) there). I don't know how to do it in base R, but data.table has a function that works here without making a copy: setattr(xx,"class","matrix")[v]. A couple R core folks replied in that thread, but what Luke Tierney's saying is a little over my head. (I've never developed a package or made a class with dispatch methods.) – Afghanistan 5/5, 2015 at 14:26

@Afghanistan : is.matrix(xx) returns TRUE since length(dim(xx)) == 2. I'm saying, that in such a case diag outputs the correct result, although I highly doubt that such a situation will ever occur in an object that "inherits" from "matrix" and diag is not that crucial a function to have been built to forestall anything crazy. And, also, c can still have uncommon class methods to even ruin c(xx)[v] e.g. as in something like c.mymatrix = function(x) stop("no 'c' method"); c(xx)[v]. I'm interested, though, as to why c remains in diag. Interesting thread! – Biracial 5/5, 2015 at 17:6

Should this question be closed since the newest update of R included a change to diag? Should someone answer giving a summary of events? That was interesting. – Detraction 22/6, 2015 at 14:35

@Steve_Corrin Ah, cool, didn't realize it was out so soon. A summary would be great; feel free to add one. I'd hate to close it as "not reproducible", though it might arguably be justified. – Afghanistan 22/6, 2015 at 14:57

If I have time tonight or tomorrow I'll write up a neat little summary. However, it may be wise to close it to avoid confusion since the question is no longer reproducible. – Detraction 22/6, 2015 at 15:4

Yeah, I could change it to "Why was..." but closure would still be justified. I'll ask in the R chat – Afghanistan 22/6, 2015 at 15:7

@Steve_Corrin Well, first comment in chat suggested an answer is a good idea. I'll edit the question to indicate that it applied to R <= 3.2.0 and post my own (very minimal) answer some time tomorrow if you don't have time to put one up before then. – Afghanistan 22/6, 2015 at 15:43

Summary

As of R version 3.2.1 (World-Famous Astronaut) diag() has received an update. The discussion moved to r-devel where it was noted that c() strips non-name attributes and may have been why it was placed there. While some people worried that removing c() would cause unknown issues on matrix-like objects, Peter Dalgaard found that, "The only case where the c() inside diag() has an effect is where M[i,j] != M[(i-1)*m+j] AND c(M) will stringize M in column-major order, so that M[i,j] == c(M)[(i-1)*m+j]."

Luke Tierney tested @Frank 's removal of c(), finding it did not effect anything on CRAN or BIOC and so was implemented to replace c(x)[...] with x[...] on line 27. This leads to relatively large speedups in diag(). Below is a speed test showing the improvement with R 3.2.1's version of diag().

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

    microbenchmark(diagOld(m),diag(m))
    Unit: microseconds
           expr        min          lq        mean      median         uq        max neval
     diagOld(m) 451189.242 526622.2775 545116.5668 531905.5635 540008.704 682223.733   100
        diag(m)    222.563    646.8675    644.7444    714.4575    740.701   1015.459   100

Detraction answered 22/6, 2015 at 19:34 Comment(2)

Thanks. Did you copy the source for the old function to diagOld and run the benchmark in R 3.2.1, then? – Afghanistan 22/6, 2015 at 19:42

Yes I was going to post the code, but did not see any reason to. Made an edit to reflect that. – Detraction 22/6, 2015 at 19:43

Hot tags

Godot Unity Godot Help Programming Godot 4.X GUI GDScript 3D 2D Physics CSharp Godot 3.X VR XR Projects C++

Recommended topics

Hot tags