首先，让我们生成一些虚拟数据以供参考：
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，对于X的第 i 列和第 j 列，位于Z[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)
}
}
但是，它并不像我想要的那么快。
因此，如果您能帮助我找到比我的堆叠循环解决方案更快的解决方案，我将不胜感激。
提前非常感谢！
PS：我已经在一个列表中存储了 13 个矩阵 X。应为所有这些矩阵执行计算。但是，我认为一旦我们找到了使用 1 个矩阵进行计算的有效方法，我就可以将这种方式与 lapply 一起使用来对完整列表进行整个操作？！
每个元素Z[i,j]都可以写成双线性形式。剩下的就是：把矩阵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 的解决方案并不比原始代码快，即它们是循环隐藏的：
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
我们可以使用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))))
这给出了与具有两个for循环的Z相同的结果。我不确定使用上述任何方法是否有任何性能提升，因为我们在这里没有做任何完全不同的事情。