为什么diag函数如此缓慢？[在R 3.2.0或更早版本中]

Question

我在阅读这个回答中的基准测试结果，想要将它们与另一个答案中使用的diag进行比较。不幸的是，看起来diag需要很长时间：

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)

评论：我使用 identical 进行了测试。我从一个回答中的“cond”来自这个作业问题。用一个整数矩阵，1:26 代替 letters 结果类似。

结果：

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

这只是一个矩阵子集操作，所以我不知道为什么会有这么多开销。查看函数内部，我看到了一些检查，然后是 c(m)[v]，其中 v 是在“vec”基准测试中使用的相同向量。计时这两个操作...

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

看起来我找到了罪魁祸首。因此，我的新问题变成了：为什么在diag中会有一个看似不必要且非常耗时的c？

- Frank

---

关于开销，我看了这个类似的问题：https://dev59.com/72Ml5IYBdhLWcg3wMkgo ，并且想到“哇，矩阵代数没有用基元编码!” - Frank

可能他们使用 c 是因为它可以删除所有属性。您可以在 r-devel 邮件列表上提问。 - Roland

@AlexA。这是针对diag的不同用法。尝试使用.Internal(diag(1, 2, 2))来查看它的作用。 - Roland

好的，我会的，谢谢！这是线程的副本：http://r.789695.n4.nabble.com/Why-is-the-diag-function-so-slow-for-extraction-td4706780.html - Frank

顺便说一句，m[seq.int(1,nc^2,nc+1)] 在我的电脑上是最快的。 - cryo111

Answer 1

简介

自从R版本3.2.1（世界著名宇航员）更新以来，diag()已经得到了更新。讨论移到了r-devel，在那里指出c()剥离非名称属性并可能是为什么放在那里的原因。虽然有些人担心删除c()会对类似矩阵的对象造成未知问题，但Peter Dalgaard发现，“diag()内部的c()只有在M [i,j]！= M [(i-1）* m + j]且c(M)将以列为主的顺序字符串化M时才起作用，M[i,j] == c(M)[(i-1)*m+j]。”

Luke Tierney测试了@Frank删除c()后，在CRAN或BIOC上没有影响，因此将 x [ ... ] 替换为[...] ，实施于第27行。这导致在diag()中速度相对较快。下面是使用R 3.2.1版本的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

- Steve Bronder

---

谢谢。你把旧函数的源代码复制到diagOld中并在R 3.2.1中运行了基准测试吗？ - Frank

是的，我本来想发布代码的，但是没有看到任何理由这样做。我已经进行了编辑以反映这一点。 - Steve Bronder