计算叉积、矩阵相乘并求和的有效方法

Question

首先，让我们生成一些虚拟数据以供参考：

X <- matrix(runif(27*27, 0, 1), nrow=27, ncol=27)
Y <- matrix(runif(27*27, 0, 1), nrow=27, ncol=27)

我有 2 个矩阵，X 和 Y。首先，我将使用 R 的命令

计算 X 前两列向量的叉积矩阵

cp <- tcrossprod(X[,1], X[,2])

结果 cp 现在与矩阵 Y 相乘，所有乘积相加：

res <- sum(cp * Y, na.rm=T)

现在我正在寻找一种快速有效的方法来在 R 中针对矩阵 X 的所有列向量组合执行此计算。结果应保存在与 X 和 Y 相同维度的第三个矩阵中，矩阵 Z , 在 Z[i,j] 对于 X.

的第 i 列和第 j 列

我已经用两个堆叠的 for 循环完成了这项工作：

Z <- matrix(nrow=27, ncol=27)
for (i in 1:ncol(X)) {
 for (j in 1:ncol(X)) {
  cp     <- tcrossprod(X[,i], X[,j])
  Z[i,j] <- sum(cp * Y)
 }
}

然而，它并没有我想要的那么快。

因此，如果您能帮我找到一个比我的堆叠 for 循环解决方案更快的解决方案，我将不胜感激。

非常感谢！

PS: 我在一个列表中存储了 13 个矩阵 X。应对所有这些矩阵执行计算。但是，我假设一旦我们找到一种使用 1 个矩阵进行计算的有效方法，我可以将这种方法与 lapply 一起使用以对完整列表进行所有操作？！

Answer 1

我们可以使用outer来申请每个列的组合

fun <- function(x, y) sum(tcrossprod(X[,x], X[,y]) *Y)
outer(seq_len(ncol(X)), seq_len(ncol(X)), Vectorize(fun))

或嵌套sapply

t(sapply(seq_len(ncol(X)), function(x) 
         sapply(seq_len(ncol(X)), function(y)  sum(tcrossprod(X[,x], X[,y]) *Y))))

或 apply

t(apply(X, 2, function(x) apply(X, 2, function(y) sum(tcrossprod(x, y) *Y))))

这与您的 Z 具有两个 for 循环的结果相同。我不确定使用上述任何方法是否有任何性能提升，因为我们在这里没有做任何完全不同的事情。

Answer 2

每个元素Z[i,j]都可以写成bilinear form。剩下的就是：将矩阵 Z.
的所有类似计算放在一起你可以这样做：

Z <- t(X) %*% Y %*% X  ### or
Z <- crossprod(X, Y) %*% X

要将此计算与您的代码进行比较：

set.seed(42)
n <- 27
X <- matrix(runif(n*n, 0, 1), nrow=n, ncol=n)
Y <- matrix(runif(n*n, 0, 1), nrow=n, ncol=n)

Z <- matrix(nrow=n, ncol=n)
for (i in 1:ncol(X)) {
  for (j in 1:ncol(X)) {
    cp     <- tcrossprod(X[,i], X[,j])
    Z[i,j] <- sum(cp * Y)
  }
}

Z2 <- t(X) %*% Y %*% X
Z3 <- crossprod(X, Y) %*% X
sum(abs(Z2-Z))
sum(abs(Z3-Z))

如果L是你的13个矩阵X的列表，你可以这样做：

lapply(L, function(X) crossprod(X, Y) %*% X)

这是基准测试：

Z1 <- function(X) {
  Z <- matrix(nrow=27, ncol=27)
  for (i in 1:ncol(X)) {
    for (j in 1:ncol(X)) {
      cp     <- tcrossprod(X[,i], X[,j])
      Z[i,j] <- sum(cp * Y)
    }
  }
  return(Z)
}

library("microbenchmark")

microbenchmark(Z1=Z1(X), Z2=t(X) %*% Y %*% X, Z3=crossprod(X, Y) %*% X)
#> microbenchmark(Z1=Z1(X), Z2=t(X) %*% Y %*% X, Z3=crossprod(X, Y) %*% X)
#Unit: microseconds
# expr      min        lq       mean    median       uq      max neval cld
#   Z1 3563.167 3671.6355 4391.00888 3721.3380 3874.617 9423.808   100   b
#   Z2   26.558   27.3420   34.31214   35.5865   39.815   56.426   100  a 
#   Z3   24.779   25.1675   27.43546   26.0965   28.034   47.268   100  a

Ronak 的解决方案并不比原始代码快，即它们 loop-hiding:

fun <- function(x, y) sum(tcrossprod(X[,x], X[,y]) *Y)

microbenchmark(Z1=Z1(X), 
               R1=outer(seq_len(ncol(X)), seq_len(ncol(X)), Vectorize(fun)), 
               R2=t(sapply(seq_len(ncol(X)), function(x) 
                 sapply(seq_len(ncol(X)), function(y)  sum(tcrossprod(X[,x], X[,y]) *Y)))),
               R3=t(apply(X, 2, function(x) apply(X, 2, function(y) sum(tcrossprod(x, y) *Y)))),
               unit="relative")
# Unit: relative
# expr      min       lq     mean   median       uq       max neval cld
#   Z1 1.000000 1.000000 1.000000 1.000000 1.000000  1.000000   100  a 
#   R1 1.207583 1.213846 1.195597 1.216147 1.223139  1.060187   100  ab
#   R2 1.225521 1.230332 1.487811 1.230852 1.299253 13.140022   100   b
#   R3 1.156546 1.158774 1.217766 1.160142 2.012623  1.098679   100  ab

计算叉积、矩阵相乘并求和的有效方法

Efficient way to calculate cross products, multiply matrix and sum it up

r

sum

multiple-columns

matrix-multiplication

cross-product