R 中的 3D 矩阵乘法
3D Matrix Multiplication in R
我有一个简单的问题。我想在不使用 for 循环的情况下将 3D 数组乘以 R 中的另一个 3D 数组。
举例:
假设我有一个 1x3 矩阵 A:
[A1, A2, A3]
我有一个 3x3 矩阵 B:
[B1, B2, B3 \
B4, B5, B6 \
B7, B8, B9]
我的主要操作是 A %*% B 生成 1x3 矩阵。
但现在我想重复该过程 10,000 次，每次使用与上面相同尺寸的不同 A 和 B。我可以使用 for-loop
for (i in 1:10000) {
A[i] %*% B[i]
}
然后我可以存储 10,000 个值。
但是有没有什么方法可以在不使用 for 循环的情况下实现同样的事情。我正在考虑可能是 3D 数组乘法。但我不确定如何在 R 中执行此操作。
Matrix A: 1 x 3 x 10000
[A1, A2, A3]
Matrix B: 3 x 3 x 10000
[B1, B2, B3
B4, B5, B6
B7, B8, B9]
此外，矢量化会有帮助吗？
你们能帮忙吗？谢谢！
如果你的A和B是list，你可以使用mapply():
> nn <- 1e1
> set.seed(1)
> A <- replicate(nn,matrix(rnorm(3),nrow=1),simplify=FALSE)
> B <- replicate(nn,matrix(rnorm(9),nrow=3),simplify=FALSE)
> head(mapply("%*%",A,B,SIMPLIFY=FALSE),3)
[[1]]
[,1] [,2] [,3]
[1,] -1.193976 0.1275999 -0.6831007
[[2]]
[,1] [,2] [,3]
[1,] 1.371143 1.860379 -1.639078
[[3]]
[,1] [,2] [,3]
[1,] 0.8250047 -0.6967286 1.949236
有多种方法可以通过数组乘法实现这一点。您付出的代价是将矩阵重新格式化为包含许多零的更大的张量。根据定义，这些是稀疏的，因此主要成本是转换的开销。当您有 10,000 个数组要相乘时，它实际上优于循环。
让 n 乘以 (A,B) 对的数量和 k=3 的维度。
最圆滑的解决方案似乎是将 A 的 n 行（n by k 矩阵）重组为 n*k by n*k block-diagonal k 乘 k 块的矩阵。块 i，i=1..n，在其顶行包含 A 的行 i，否则为零。将此（在右侧）乘以 B（排列为 k*n 乘以 k 矩阵，该矩阵由 "stack" 的 n 维块组成 k by k) 计算所有单个产品，将它们存放在结果的第 1、k+1、2k+1、...行，在那里可以挑选它们。
f3 <- function(a, b) {
require(RcppArmadillo) # sparseMatrix package
n <- dim(b)[3]
k <- dim(b)[2]
i0 <- (1:n-1)*k+1
i <- rep(i0, each=k)
j <- 1:(k*n)
aa <- sparseMatrix(i, j, x=c(t(a)), dims=c(n*k, n*k))
bb <- matrix(aperm(b, c(1,3,2)), nrow=n*k)
t((aa %*% bb)[i0, ])
}
如您所见，数组操作是基本的：创建稀疏矩阵、转置数组（使用 aperm 和 t）以及相乘。它 returns 它的结果是 k by n 数组（如果你愿意，你可以转置），每列一个结果向量。
作为测试，这里有一个使用相同数组数据结构的 brute-force 循环。
f1 <- function(a, b) sapply(1:nrow(a), function(i) a[i,] %*% b[,,i])
我们可以将这些解决方案应用于相同的输入并比较结果：
#
# Create random matrices for testing.
#
k <- 3
n <- 1e6 # Number of (a,B) pairs
a <- matrix(runif(k*n), ncol=k)
b <- array(runif(k^2*n), dim=c(k,k,n))
system.time(c1 <- f1(a,b)) # 4+ seconds
system.time(c3 <- f3(a,b)) # 2/3 second
mean((c1-c3)^2) # Want around 10^-32 or less
结果并不完全相等，但它们的均方差小于 10^-32，表明它们在浮点舍入误差上可以被认为是相同的。
array-oriented 过程 f3 最初比循环过程 f1 慢，但在 n 为 10,000 时赶上了。之后它的速度大约是原来的两倍或更好（在这台机器上；YMMV）。两种算法都应在 n 内线性扩展（时间表明它们确实如此，至少扩展到 n=10,000,000）。
for-loop比你想象的更高效
你的 n (A,B) 对相乘问题并不等同于通常意义上的张量乘法，尽管 whuber 提出了一种非常巧妙的方法，通过堆叠 B 将其转化为矩阵乘法作为稀疏矩阵中的块。
你说你想避免for-loop，但for-loop方法在高效编程时实际上非常有竞争力，我建议你重新考虑它。
我将使用与 whuber 相同的符号，A 的维度为 n x k，B 的维度为 k x k x n，例如：
n <- 1e4
k <- 3
A <- array(rnorm(k*n),c(n,k))
B <- array(rnorm(k*k*n),c(k,k,n))
一个简单高效的 for-loop 解决方案是这样的
justAForLoop <- function(A,B) {
n <- nrow(A)
for (i in 1:n) A[i,] <- A[i,] %*% B[,,i]
A
}
生成一个 n x k 矩阵的结果。
我修改了whuber的f3函数加载Matrix包，否则sparseMatrix函数不可用。我的 f3 版本比原来的版本稍微快一点，因为我在返回结果之前消除了最后一个矩阵转置。
通过此修改，它 returns 与 justAForLoop.
相同的数值结果
f3 <- function(a, b) {
require(Matrix)
n <- dim(b)[3]
k <- dim(b)[2]
i0 <- (1:n-1)*k+1
i <- rep(i0, each=k)
j <- 1:(k*n)
aa <- sparseMatrix(i, j, x=c(t(a)), dims=c(n*k, n*k))
bb <- matrix(aperm(b, c(1,3,2)), nrow=n*k)
(aa %*% bb)[i0, ]
}
现在我在新的 R 会话中重新运行 whuber 的模拟：
> k <- 3
> n <- 1e6
> a <- matrix(runif(k*n), ncol=k)
> b <- array(runif(k^2*n), dim=c(k,k,n))
>
> system.time(c1 <- f1(a,b))
user system elapsed
3.40 0.09 3.50
> system.time(c3 <- f3(a,b))
Loading required package: Matrix
user system elapsed
1.06 0.24 1.30
> system.time(c4 <- justAForLoop(a,b))
user system elapsed
1.27 0.00 1.26
for-loop 方法实际上是最快的，只差一点点。它比依赖于 sapply 的 f1 快得多。 （我的机器是 Windows 10 PC，32Gb RAM 运行ning R 3.6.0）。
如果我第二次 运行 所有三种方法，那么 f3 将成为最快的，因为这次 Matrix 包已经在搜索路径中并且不必重新加载：
> system.time(c1 <- f1(a,b))
user system elapsed
3.23 0.04 3.26
> system.time(c3 <- f3(a,b))
user system elapsed
0.33 0.20 0.53
> system.time(c4 <- justAForLoop(a,b))
user system elapsed
1.28 0.01 1.30
但是 f3 使用的 RAM 比 for-loop 多。在我的电脑上，我可以 运行 justAForLoop 使用 n=1e8 成功，而 f1 和 f3 都 运行 内存不足并且失败。
总结
直接 for-loop 方法比 sapply 更有效。
对于 n=10,000 矩阵乘法的问题，运行宁 for-loop 简单高效，耗时 <0.02 秒。相比之下，仅加载具有稀疏矩阵函数的包大约需要 2/3 秒。
对于 n 1-1000 万之间的情况，whuber 的稀疏矩阵解决方案开始表现出色，尤其是在已加载 Matrix 包的情况下。
for-loop 使用三种方法中最少的 RAM。对于 n 在我的 32Gb RAM PC 上的 1 亿，只有 for-loop 方法有效。
我有一个简单的问题。我想在不使用 for 循环的情况下将 3D 数组乘以 R 中的另一个 3D 数组。
举例:
假设我有一个 1x3 矩阵 A:
[A1, A2, A3]
我有一个 3x3 矩阵 B:
[B1, B2, B3 \
B4, B5, B6 \
B7, B8, B9]
我的主要操作是 A %*% B 生成 1x3 矩阵。
但现在我想重复该过程 10,000 次，每次使用与上面相同尺寸的不同 A 和 B。我可以使用 for-loop
for (i in 1:10000) {
A[i] %*% B[i]
}
然后我可以存储 10,000 个值。
但是有没有什么方法可以在不使用 for 循环的情况下实现同样的事情。我正在考虑可能是 3D 数组乘法。但我不确定如何在 R 中执行此操作。
Matrix A: 1 x 3 x 10000
[A1, A2, A3]
Matrix B: 3 x 3 x 10000
[B1, B2, B3
B4, B5, B6
B7, B8, B9]
此外，矢量化会有帮助吗？
你们能帮忙吗？谢谢！
如果你的A和B是list，你可以使用mapply():
> nn <- 1e1
> set.seed(1)
> A <- replicate(nn,matrix(rnorm(3),nrow=1),simplify=FALSE)
> B <- replicate(nn,matrix(rnorm(9),nrow=3),simplify=FALSE)
> head(mapply("%*%",A,B,SIMPLIFY=FALSE),3)
[[1]]
[,1] [,2] [,3]
[1,] -1.193976 0.1275999 -0.6831007
[[2]]
[,1] [,2] [,3]
[1,] 1.371143 1.860379 -1.639078
[[3]]
[,1] [,2] [,3]
[1,] 0.8250047 -0.6967286 1.949236
有多种方法可以通过数组乘法实现这一点。您付出的代价是将矩阵重新格式化为包含许多零的更大的张量。根据定义，这些是稀疏的，因此主要成本是转换的开销。当您有 10,000 个数组要相乘时，它实际上优于循环。
让 n 乘以 (A,B) 对的数量和 k=3 的维度。
最圆滑的解决方案似乎是将 A 的 n 行（n by k 矩阵）重组为 n*k by n*k block-diagonal k 乘 k 块的矩阵。块 i，i=1..n，在其顶行包含 A 的行 i，否则为零。将此（在右侧）乘以 B（排列为 k*n 乘以 k 矩阵，该矩阵由 "stack" 的 n 维块组成 k by k) 计算所有单个产品，将它们存放在结果的第 1、k+1、2k+1、...行，在那里可以挑选它们。
f3 <- function(a, b) {
require(RcppArmadillo) # sparseMatrix package
n <- dim(b)[3]
k <- dim(b)[2]
i0 <- (1:n-1)*k+1
i <- rep(i0, each=k)
j <- 1:(k*n)
aa <- sparseMatrix(i, j, x=c(t(a)), dims=c(n*k, n*k))
bb <- matrix(aperm(b, c(1,3,2)), nrow=n*k)
t((aa %*% bb)[i0, ])
}
如您所见，数组操作是基本的：创建稀疏矩阵、转置数组（使用 aperm 和 t）以及相乘。它 returns 它的结果是 k by n 数组（如果你愿意，你可以转置），每列一个结果向量。
作为测试，这里有一个使用相同数组数据结构的 brute-force 循环。
f1 <- function(a, b) sapply(1:nrow(a), function(i) a[i,] %*% b[,,i])
我们可以将这些解决方案应用于相同的输入并比较结果：
#
# Create random matrices for testing.
#
k <- 3
n <- 1e6 # Number of (a,B) pairs
a <- matrix(runif(k*n), ncol=k)
b <- array(runif(k^2*n), dim=c(k,k,n))
system.time(c1 <- f1(a,b)) # 4+ seconds
system.time(c3 <- f3(a,b)) # 2/3 second
mean((c1-c3)^2) # Want around 10^-32 or less
结果并不完全相等，但它们的均方差小于 10^-32，表明它们在浮点舍入误差上可以被认为是相同的。
array-oriented 过程 f3 最初比循环过程 f1 慢，但在 n 为 10,000 时赶上了。之后它的速度大约是原来的两倍或更好（在这台机器上；YMMV）。两种算法都应在 n 内线性扩展（时间表明它们确实如此，至少扩展到 n=10,000,000）。
for-loop比你想象的更高效
你的 n (A,B) 对相乘问题并不等同于通常意义上的张量乘法，尽管 whuber 提出了一种非常巧妙的方法，通过堆叠 B 将其转化为矩阵乘法作为稀疏矩阵中的块。
你说你想避免for-loop，但for-loop方法在高效编程时实际上非常有竞争力，我建议你重新考虑它。
我将使用与 whuber 相同的符号，A 的维度为 n x k，B 的维度为 k x k x n，例如：
n <- 1e4
k <- 3
A <- array(rnorm(k*n),c(n,k))
B <- array(rnorm(k*k*n),c(k,k,n))
一个简单高效的 for-loop 解决方案是这样的
justAForLoop <- function(A,B) {
n <- nrow(A)
for (i in 1:n) A[i,] <- A[i,] %*% B[,,i]
A
}
生成一个 n x k 矩阵的结果。
我修改了whuber的f3函数加载Matrix包，否则sparseMatrix函数不可用。我的 f3 版本比原来的版本稍微快一点，因为我在返回结果之前消除了最后一个矩阵转置。
通过此修改，它 returns 与 justAForLoop.
f3 <- function(a, b) {
require(Matrix)
n <- dim(b)[3]
k <- dim(b)[2]
i0 <- (1:n-1)*k+1
i <- rep(i0, each=k)
j <- 1:(k*n)
aa <- sparseMatrix(i, j, x=c(t(a)), dims=c(n*k, n*k))
bb <- matrix(aperm(b, c(1,3,2)), nrow=n*k)
(aa %*% bb)[i0, ]
}
现在我在新的 R 会话中重新运行 whuber 的模拟：
> k <- 3
> n <- 1e6
> a <- matrix(runif(k*n), ncol=k)
> b <- array(runif(k^2*n), dim=c(k,k,n))
>
> system.time(c1 <- f1(a,b))
user system elapsed
3.40 0.09 3.50
> system.time(c3 <- f3(a,b))
Loading required package: Matrix
user system elapsed
1.06 0.24 1.30
> system.time(c4 <- justAForLoop(a,b))
user system elapsed
1.27 0.00 1.26
for-loop 方法实际上是最快的，只差一点点。它比依赖于 sapply 的 f1 快得多。 （我的机器是 Windows 10 PC，32Gb RAM 运行ning R 3.6.0）。
如果我第二次 运行 所有三种方法，那么 f3 将成为最快的，因为这次 Matrix 包已经在搜索路径中并且不必重新加载：
> system.time(c1 <- f1(a,b))
user system elapsed
3.23 0.04 3.26
> system.time(c3 <- f3(a,b))
user system elapsed
0.33 0.20 0.53
> system.time(c4 <- justAForLoop(a,b))
user system elapsed
1.28 0.01 1.30
但是 f3 使用的 RAM 比 for-loop 多。在我的电脑上，我可以 运行 justAForLoop 使用 n=1e8 成功，而 f1 和 f3 都 运行 内存不足并且失败。
总结
直接 for-loop 方法比 sapply 更有效。
对于 n=10,000 矩阵乘法的问题，运行宁 for-loop 简单高效，耗时 <0.02 秒。相比之下，仅加载具有稀疏矩阵函数的包大约需要 2/3 秒。
对于 n 1-1000 万之间的情况，whuber 的稀疏矩阵解决方案开始表现出色，尤其是在已加载 Matrix 包的情况下。
for-loop 使用三种方法中最少的 RAM。对于 n 在我的 32Gb RAM PC 上的 1 亿，只有 for-loop 方法有效。