加快对称矩阵的计算；外用答案

【问题标题】：Speeding up calculation of symmetric matrices; use of outer加快对称矩阵的计算；外用
【发布时间】：2019-02-22 08:31:47
【问题描述】：

我需要加快生成对称矩阵的计算。目前我有这样的事情：

X <- 1:50
Y<- 1:50
M <- outer(X, Y, FUN = myfun)

其中 myfun 是一个相当复杂的矢量化但对称函数 (myfun(x, y) = myfun(y, x))。

所以我的代码不必要地浪费时间计算下三角矩阵和上三角矩阵。

如何在不使用缓慢的 for 循环的情况下避免这种重复？

【问题讨论】：

这与您的问题没有直接关系，但我相信它会有所帮助 - the R open distribution from Microsoft 使用多线程数学库 - 英特尔 MKL 显着提升了矩阵运算（超过 100 倍的加速是一些情况）与单线程 BLAS/LAPACK 库相比。我发现安装它非常值得。
M不是方阵怎么可能是对称矩阵？
一种方法是使用记忆并使用包装器来确保在调用myfun(y,x) 时它会查找myfun(x,y)。这可能有用：github.com/r-lib/memoise
@989 不一定；我只是用整数来说明问题。
@989 该函数的形式为 Vectorize(fn(x, y))。 fn 接受两个标量并返回一个标量。

标签： r matrix vectorization

【解决方案1】：

如果你的函数很慢并且时间随输入的大小而变化，你可以使用combn:

X <- 1:50
Y <- 1:50

#a slow function
myfun <- function(x, y) {
  res <- x * NA
  for (i in seq_along(x)) {
    Sys.sleep(0.01)
    res[i] <- x[i] * y[i]
    }
  res
}

system.time(M <- outer(X, Y, FUN = myfun))
#user  system elapsed 
#0.00    0.00   26.41 

system.time({
  inds <- combn(seq_len(length(X)), 2)
  M1 <- matrix(ncol = length(X), nrow = length(Y))

  M1[lower.tri(M1)] <-  myfun(X[inds[1,]], Y[inds[2,]])
  M1[upper.tri(M1)] <- t(M1)[upper.tri(M1)]
  diag(M1) <- myfun(X, Y)
})
#user  system elapsed 
#0.00    0.00   13.41

all.equal(M, M1)
#[1] TRUE

但是，最好的解决方案可能是通过 Rcpp 在 C++ 中实现。

【讨论】：

非常感谢。我不会写 C++，但将你的 R 代码的时间减半是值得的。
也许seq_along(X) 代替seq_len(length(X))？
当然。不过，这里的速度并没有什么不同。