R：部分修改矩阵对角线的有效方法

Question

我有一个尺寸范围从 100x100 到 10,000x10,000 的方阵。该矩阵表示函数的参数值。我经历了一个循环，尝试各种参数组合。每次迭代

W

都有不同的值，因此我必须更新矩阵中具有值

W

的位置。这恰好是对角线的偶数条目，所以

[2,2]

、

[4,4]

等。我很好奇是否有比我当前的方法更有效的方法：

W<-treeDepth*newVar
for (iW in evenDiag) {
  matrixSource[iW,iW]<-W
}

到目前为止，我只测试了 134x134 大小的矩阵，但根据

profvis

的分析，循环似乎是最快的方法。

当我尝试时

diag(matrixSource)[evenDiag]<-W

似乎花费了相似的时间，但平均每 5 次调用左右后它就开始调用 < GC>，但实际调用的时间 < GC> 似乎是随机的。

我认为 < GC> 是垃圾收集，但无论它是什么，它都会花费很长时间，并且当我有上面的循环版本时，它很少被调用。

我认为有比逐一循环更好的方法是错误的吗？用

Rcpp

编写循环会使其更快吗？ Hadley 在“Advanced R”中并没有说使用 Rcpp 将值写入矩阵会更快，因此可能不会。如果确实如此，我如何将这一小行更改为 Rcpp 而不使用 C++ 函数或任何复杂的东西（我不懂 C++）。

根据我的研究，不可能只写这样的东西

matrixSource[evenDiag,evenDiag]<-W

但如果是的话，R 将擅长矢量化处理。

最好的方法是什么？

如果有帮助，上下文是需要将矩阵输入

negLogLik<- -mvtnorm::dmvnorm(flattenedData,sigma=matrixSource,log=T,checkSymmetry = FALSE)

在该函数内部，它被输入到 chol() （有时并行调用 < GC>）

因此，如果有一种方法可以修改该函数以仅处理不完全的部分矩阵 创建这样也许我只需要分配 W 的值一次，那也很好。

我还需要一种方法来将主对角线正上方和正下方的对角线分配为相同的值 Z。有没有好的方法可以做到这一点？

谢谢你<3

Answer 1

@MichaelChirico 的建议似乎是最有效的。它以矩阵维度的时间常数运行，与

diag<-

不同。
下面使用的两种方法是

# create an index
i <- rep(c(FALSE, TRUE), ncol(A) %/% 2)
j <- which(i)

# can also be diag(A)[i]
diag(A)[j] <- W

# cannot be cbind(i, i)
A[cbind(j, j)] <- 0

cbind

显然效率更高。

library(microbenchmark)
library(ggplot2)

testFun <- function(N) {
  out <- lapply(seq.int(N)[-1L], \(n) {
    A <- matrix(1:n^2, n, n)
    W <- rep(0L, ncol(A) %/% 2L)
    i <- rep(c(FALSE, TRUE), ncol(A) %/% 2)
    j <- which(i)
    mb <- microbenchmark(
      diag = {diag(A)[i] <- W},
      cbind = {A[cbind(j, j)] <- W}
    )
    mb$dim <- n
    mb
  })
  out <- do.call(rbind, out)
  aggregate(time ~ ., out, median)
}

testData <- testFun(1000)
ggplot(testData, aes(dim, time, colour = expr)) +
  geom_line() +
  geom_point() +
  theme_bw()

^{创建于 2024-01-14，使用 reprex v2.0.2}