进入矩阵乘法的 R 值
R values that go into matrix multiplication
保存进入矩阵乘法（不带 0）的唯一值的最快方法是什么？
例如，如果我有一个 data.table 对象
library(data.table)
A = data.table(j3=c(3,0,3),j5=c(0,5,5),j7=c(0,7,0),j8=c(8,0,8))
我想看看哪些唯一值进入 A*transpose(A)（或 as.matrix(A) %*% as.matrix(t(A))）。现在，我可以使用 for 循环来做到这一点：
B=t(A)
L = list()
models = c('A1','A2','A3')
for(i in 1:nrow(A)){
for(j in 1:ncol(B)){
u = union(unlist(A[i,]),B[,j])
u = u[u!=0] # remove 0
L[[paste(models[i],models[j])]]= u
}
}
但是，有没有更快更高效的 RAM 方法？输出不一定是列表对象，在我的例子中，它可以是 data.table (data.frame) 还有。此外，值的顺序并不重要。例如，3 5 8 与 5 3 8、8 5 3 等
一样好
感谢任何帮助。
编辑：所以 as.matrix(A) %*% as.matrix(t(A)) 是：
[,1] [,2] [,3]
[1,] 73 0 73
[2,] 0 74 25
[3,] 73 25 98
第一个元素计算为 3*3+0*0+0*0+8*8 = 73，第二个元素为 3*0+0*5+0*7+8*0 = 0，依此类推。我需要进行此计算但不包含 0 的唯一数字。
因此输出（保存在列表 L 中）为：
> L
$`A1 A1`
[1] 3 8
$`A1 A2`
[1] 3 8 5 7
$`A1 A3`
[1] 3 8 5
$`A2 A1`
[1] 5 7 3 8
$`A2 A2`
[1] 5 7
$`A2 A3`
[1] 5 7 3 8
$`A3 A1`
[1] 3 5 8
$`A3 A2`
[1] 3 5 8 7
$`A3 A3`
[1] 3 5 8
再次声明，输出不必是列表对象。如果可行的话，我更喜欢 data.table。是否可以将我的方法重写为 Rcpp 函数？
感谢您在编辑中发布附加信息。从您发布的内容来看，对于矩阵或数据的所有行对 table A，您希望这两行中的唯一非零值。
为了有效地做到这一点，我建议确保 A 是一个矩阵。数据帧或数据 tables 中的行索引比在矩阵中慢很多。 （列索引可以更快，但我怀疑是否值得转置 table 来获得它。）
一旦有了矩阵，A[i, ] 就是一个向量，其中包含行 i 中的值，这是一个非常快速的计算。您需要 c(A[i, ], A[j, ]) 中的唯一非零值。 unique 函数将产生这个，但不会遗漏零。我建议尝试一下。根据每一行的内容，可以想象在计算唯一条目之前先将零从行中删除可能比计算所有唯一值然后删除 0 更快或更慢。
你说你想做几百行，但每一行都很长。我猜你将无法在嵌套循环上有太大改进：时间将花在每个条目上，而不是循环上。但是，您可以使用 apply() 函数尝试矢量化，例如
result <- vector("list", nrows)
for (i in 1:nrows)
result[[i]] <- apply(A, 1, function(row) setdiff(unique(c(row, A[i,])), 0))
这将给出列表列表；如果你想检查条目 i，j，你可以使用 result[[c(i,j)]].
潜在的优化
跟进@user2554330 的回答，请注意，如果 A 是一个 m-by-n 矩阵，则 AAT = A %*% t(A)（相当于 tcrossprod(A) ) 是一个 m×m 对称矩阵。 AAT[i, j] 和 AAT[j, i] 是使用 A 的相同条目计算的，因此您只需要检查 m*(m+1)/2 行 A，而不是 m*m.
您可以做得更好，方法是在 配对前 查找并缓存每一行的唯一元素。以这种方式进行预处理避免了冗余计算，并且在 m << n.
时应该会显着提高性能
限制
问题的另一个方面是 unique 如何在幕后工作。 unique 有一个参数 nmax，您可以使用它来指定唯一元素的预期最大数量。来自 ?duplicated:
Except for factors, logical and raw vectors the default nmax = NA is equivalent to nmax = length(x). Since a hash table of size 8*nmax bytes is allocated, setting nmax suitably can save large amounts of memory. For factors it is automatically set to the smaller of length(x) and the number of levels plus one (for NA). If nmax is set too small there is liable to be an error: nmax = 1 is silently ignored.
Long vectors are supported for the default method of duplicated, but may only be usable if nmax is supplied.
这些评论也适用于 unique。由于您有一个 300-by-4e+07 矩阵，您将评估（通过预处理）：
unique(<4e+07-length vector>), 300 次,unique(<up to 8e+07-length vector>)，299*300/2次。
如果您对可能允许您设置 nmax 的矩阵一无所知，这会消耗大量内存。如果您无法访问许多 CPU，则可能需要很长时间。
所以我同意要求您考虑为什么您需要这样做的评论以及您的潜在问题是否有更好的解决方案。
两个答案
FWIW，这里有两种解决您的一般问题的方法，它们实际上利用了对称性。 f 和 g 没有和有预处理。 [[.utri 允许您从 return 值中提取元素，这是一个 m*(m+1)/2 长度的列表，就好像它是一个 m-by-m 矩阵一样。 as.matrix.utri 构建完整的对称 m-by-m 列表矩阵。
f <- function(A, nmax = NA) {
a <- seq_len(nrow(A))
J <- cbind(sequence(a), rep.int(a, a))
FUN <- function(i) {
if (i[1L] == i[2L]) {
x <- A[i[1L], ]
} else {
x <- c(A[i[1L], ], A[i[2L], ])
}
unique.default(x[x != 0], nmax = nmax)
}
res <- apply(J, 1L, FUN, simplify = FALSE)
class(res) <- "utri"
res
}
g <- function(A, nmax = NA) {
l <- lapply(asplit(A, 1L), function(x) unique.default(x[x != 0], nmax = nmax))
a <- seq_along(l)
J <- cbind(sequence(a), rep.int(a, a))
FUN <- function(i) {
if (i[1L] == i[2L]) {
l[[i[1L]]]
} else {
unique.default(c(l[[i[1L]]], l[[i[2L]]]))
}
}
res <- apply(J, 1L, FUN, simplify = FALSE)
class(res) <- "utri"
res
}
`[[.utri` <- function(x, i, j) {
stopifnot(length(i) == 1L, length(j) == 1L)
class(x) <- NULL
if (i <= j) {
x[[i + (j * (j - 1L)) %/% 2L]]
} else {
x[[j + (i * (i - 1L)) %/% 2L]]
}
}
as.matrix.utri <- function(x) {
p <- length(x)
n <- as.integer(round(0.5 * (-1 + sqrt(1 + 8 * p))))
i <- rep.int(seq_len(n), n)
j <- rep.int(seq_len(n), rep.int(n, n))
r <- i > j
ir <- i[r]
i[r] <- j[r]
j[r] <- ir
res <- x[i + (j * (j - 1L)) %/% 2L]
dim(res) <- c(n, n)
res
}
下面是对 4×4 整数矩阵的简单测试：
mkA <- function(m, n) {
A <- sample(0:(n - 1L), size = as.double(m) * n, replace = TRUE,
prob = rep.int(c(n - 1, 1), c(1L, n - 1L)))
dim(A) <- c(m, n)
A
}
set.seed(1L)
A <- mkA(4L, 4L)
A
## [,1] [,2] [,3] [,4]
## [1,] 0 0 2 3
## [2,] 0 1 0 0
## [3,] 2 1 0 3
## [4,] 1 2 0 0
identical(f(A), gA <- g(A))
## [1] TRUE
gA[[1L, 1L]] # used for 'tcrossprod(A)[1L, 1L]'
## [1] 2 3
gA[[1L, 2L]] # used for 'tcrossprod(A)[1L, 2L]'
## [1] 2 3 1
gA[[2L, 1L]] # used for 'tcrossprod(A)[2L, 1L]'
## [1] 2 3 1
gA # under the hood, an 'm*(m+1)/2'-length list
## [[1]]
## [1] 2 3
##
## [[2]]
## [1] 2 3 1
##
## [[3]]
## [1] 1
##
## [[4]]
## [1] 2 3 1
##
## [[5]]
## [1] 1 2 3
##
## [[6]]
## [1] 2 1 3
##
## [[7]]
## [1] 2 3 1
##
## [[8]]
## [1] 1 2
##
## [[9]]
## [1] 2 1 3
##
## [[10]]
## [1] 1 2
##
## attr(,"class")
## [1] "utri"
mgA <- as.matrix(gA) # the full, symmetric, 'm'-by-'m' list matrix
mgA
## [,1] [,2] [,3] [,4]
## [1,] integer,2 integer,3 integer,3 integer,3
## [2,] integer,3 1 integer,3 integer,2
## [3,] integer,3 integer,3 integer,3 integer,3
## [4,] integer,3 integer,2 integer,3 integer,2
mgA[1L, ] # used for first row of 'tcrossprod(A)'
## [[1]]
## [1] 2 3
##
## [[2]]
## [1] 2 3 1
##
## [[3]]
## [1] 2 3 1
##
## [[4]]
## [1] 2 3 1
## If you need names
dimnames(mgA) <- rep.int(list(sprintf("A%d", seq_len(nrow(mgA)))), 2L)
mgA["A1", ]
## $A1
## [1] 2 3
##
## $A2
## [1] 2 3 1
##
## $A3
## [1] 2 3 1
##
## $A4
## [1] 2 3 1
## If you need an 'm'-by-'m' 'data.table' result
DT <- data.table::as.data.table(mgA)
DT
## A1 A2 A3 A4
## 1: 2,3 2,3,1 2,3,1 2,3,1
## 2: 2,3,1 1 1,2,3 1,2
## 3: 2,3,1 1,2,3 2,1,3 2,1,3
## 4: 2,3,1 1,2 2,1,3 1,2
这里有两个关于两个大整数矩阵的基准测试，表明预处理可以提供很大帮助：
set.seed(1L)
A <- mkA(100L, 1e+04L)
microbenchmark::microbenchmark(f(A), g(A), times = 10L, setup = gc(FALSE))
## Unit: milliseconds
## expr min lq mean median uq max neval
## f(A) 2352.0572 2383.3100 2435.7954 2403.8968 2431.6214 2619.553 10
## g(A) 843.0206 852.5757 858.7262 858.2746 863.8239 881.450 10
A <- mkA(100L, 1e+06L)
microbenchmark::microbenchmark(f(A), g(A), times = 10L, setup = gc(FALSE))
## Unit: seconds
## expr min lq mean median uq max neval
## f(A) 290.93327 295.54319 302.57001 301.17810 307.50226 318.14203 10
## g(A) 72.85608 73.83614 76.67941 76.57313 77.78056 83.73388 10
也许我们可以试试这个
f <- function(A, models) {
AA <- replace(A, A == 0, NA)
setNames(
c(t(outer(
1:nrow(A),
1:nrow(A),
Vectorize(function(x, y) unique(na.omit(c(t(AA[c(x, y)])))))
))),
t(outer(models, models, paste))
)
}
这给出了
$`A1 A1`
[1] 3 8
$`A1 A2`
[1] 3 8 5 7
$`A1 A3`
[1] 3 8 5
$`A2 A1`
[1] 5 7 3 8
$`A2 A2`
[1] 5 7
$`A2 A3`
[1] 5 7 3 8
$`A3 A1`
[1] 3 5 8
$`A3 A2`
[1] 3 5 8 7
$`A3 A3`
[1] 3 5 8
如果你在意速度，可以试试
lst <- asplit(replace(A, A == 0, NA), 1)
mat <- matrix(list(), nrow = nrow(A), ncol = nrow(A))
mat[lower.tri(mat)] <- combn(lst, 2, function(...) unique(na.omit(unlist(...))), simplify = FALSE)
mat[upper.tri(mat)] <- t(mat)[upper.tri(mat)]
diag(mat) <- Map(function(x) unname(x)[!is.na(x)], lst)
L <- c(t(mat))
保存进入矩阵乘法（不带 0）的唯一值的最快方法是什么？
例如，如果我有一个 data.table 对象
library(data.table)
A = data.table(j3=c(3,0,3),j5=c(0,5,5),j7=c(0,7,0),j8=c(8,0,8))
我想看看哪些唯一值进入 A*transpose(A)（或 as.matrix(A) %*% as.matrix(t(A))）。现在，我可以使用 for 循环来做到这一点：
B=t(A)
L = list()
models = c('A1','A2','A3')
for(i in 1:nrow(A)){
for(j in 1:ncol(B)){
u = union(unlist(A[i,]),B[,j])
u = u[u!=0] # remove 0
L[[paste(models[i],models[j])]]= u
}
}
但是，有没有更快更高效的 RAM 方法？输出不一定是列表对象，在我的例子中，它可以是 data.table (data.frame) 还有。此外，值的顺序并不重要。例如，3 5 8 与 5 3 8、8 5 3 等
感谢任何帮助。
编辑：所以 as.matrix(A) %*% as.matrix(t(A)) 是：
[,1] [,2] [,3]
[1,] 73 0 73
[2,] 0 74 25
[3,] 73 25 98
第一个元素计算为 3*3+0*0+0*0+8*8 = 73，第二个元素为 3*0+0*5+0*7+8*0 = 0，依此类推。我需要进行此计算但不包含 0 的唯一数字。
因此输出（保存在列表 L 中）为：
> L
$`A1 A1`
[1] 3 8
$`A1 A2`
[1] 3 8 5 7
$`A1 A3`
[1] 3 8 5
$`A2 A1`
[1] 5 7 3 8
$`A2 A2`
[1] 5 7
$`A2 A3`
[1] 5 7 3 8
$`A3 A1`
[1] 3 5 8
$`A3 A2`
[1] 3 5 8 7
$`A3 A3`
[1] 3 5 8
再次声明，输出不必是列表对象。如果可行的话，我更喜欢 data.table。是否可以将我的方法重写为 Rcpp 函数？
感谢您在编辑中发布附加信息。从您发布的内容来看，对于矩阵或数据的所有行对 table A，您希望这两行中的唯一非零值。
为了有效地做到这一点，我建议确保 A 是一个矩阵。数据帧或数据 tables 中的行索引比在矩阵中慢很多。 （列索引可以更快，但我怀疑是否值得转置 table 来获得它。）
一旦有了矩阵，A[i, ] 就是一个向量，其中包含行 i 中的值，这是一个非常快速的计算。您需要 c(A[i, ], A[j, ]) 中的唯一非零值。 unique 函数将产生这个，但不会遗漏零。我建议尝试一下。根据每一行的内容，可以想象在计算唯一条目之前先将零从行中删除可能比计算所有唯一值然后删除 0 更快或更慢。
你说你想做几百行，但每一行都很长。我猜你将无法在嵌套循环上有太大改进：时间将花在每个条目上，而不是循环上。但是，您可以使用 apply() 函数尝试矢量化，例如
result <- vector("list", nrows)
for (i in 1:nrows)
result[[i]] <- apply(A, 1, function(row) setdiff(unique(c(row, A[i,])), 0))
这将给出列表列表；如果你想检查条目 i，j，你可以使用 result[[c(i,j)]].
潜在的优化
跟进@user2554330 的回答，请注意，如果 A 是一个 m-by-n 矩阵，则 AAT = A %*% t(A)（相当于 tcrossprod(A) ) 是一个 m×m 对称矩阵。 AAT[i, j] 和 AAT[j, i] 是使用 A 的相同条目计算的，因此您只需要检查 m*(m+1)/2 行 A，而不是 m*m.
您可以做得更好，方法是在 配对前 查找并缓存每一行的唯一元素。以这种方式进行预处理避免了冗余计算，并且在 m << n.
限制
问题的另一个方面是 unique 如何在幕后工作。 unique 有一个参数 nmax，您可以使用它来指定唯一元素的预期最大数量。来自 ?duplicated:
Except for factors, logical and raw vectors the default
nmax = NAis equivalent tonmax = length(x). Since a hash table of size8*nmaxbytes is allocated, settingnmaxsuitably can save large amounts of memory. For factors it is automatically set to the smaller oflength(x)and the number of levels plus one (forNA). Ifnmaxis set too small there is liable to be an error:nmax = 1is silently ignored.
Long vectors are supported for the default method of
duplicated, but may only be usable ifnmaxis supplied.
这些评论也适用于 unique。由于您有一个 300-by-4e+07 矩阵，您将评估（通过预处理）：
unique(<4e+07-length vector>), 300 次,unique(<up to 8e+07-length vector>)，299*300/2次。
如果您对可能允许您设置 nmax 的矩阵一无所知，这会消耗大量内存。如果您无法访问许多 CPU，则可能需要很长时间。
所以我同意要求您考虑为什么您需要这样做的评论以及您的潜在问题是否有更好的解决方案。
两个答案
FWIW，这里有两种解决您的一般问题的方法，它们实际上利用了对称性。 f 和 g 没有和有预处理。 [[.utri 允许您从 return 值中提取元素，这是一个 m*(m+1)/2 长度的列表，就好像它是一个 m-by-m 矩阵一样。 as.matrix.utri 构建完整的对称 m-by-m 列表矩阵。
f <- function(A, nmax = NA) {
a <- seq_len(nrow(A))
J <- cbind(sequence(a), rep.int(a, a))
FUN <- function(i) {
if (i[1L] == i[2L]) {
x <- A[i[1L], ]
} else {
x <- c(A[i[1L], ], A[i[2L], ])
}
unique.default(x[x != 0], nmax = nmax)
}
res <- apply(J, 1L, FUN, simplify = FALSE)
class(res) <- "utri"
res
}
g <- function(A, nmax = NA) {
l <- lapply(asplit(A, 1L), function(x) unique.default(x[x != 0], nmax = nmax))
a <- seq_along(l)
J <- cbind(sequence(a), rep.int(a, a))
FUN <- function(i) {
if (i[1L] == i[2L]) {
l[[i[1L]]]
} else {
unique.default(c(l[[i[1L]]], l[[i[2L]]]))
}
}
res <- apply(J, 1L, FUN, simplify = FALSE)
class(res) <- "utri"
res
}
`[[.utri` <- function(x, i, j) {
stopifnot(length(i) == 1L, length(j) == 1L)
class(x) <- NULL
if (i <= j) {
x[[i + (j * (j - 1L)) %/% 2L]]
} else {
x[[j + (i * (i - 1L)) %/% 2L]]
}
}
as.matrix.utri <- function(x) {
p <- length(x)
n <- as.integer(round(0.5 * (-1 + sqrt(1 + 8 * p))))
i <- rep.int(seq_len(n), n)
j <- rep.int(seq_len(n), rep.int(n, n))
r <- i > j
ir <- i[r]
i[r] <- j[r]
j[r] <- ir
res <- x[i + (j * (j - 1L)) %/% 2L]
dim(res) <- c(n, n)
res
}
下面是对 4×4 整数矩阵的简单测试：
mkA <- function(m, n) {
A <- sample(0:(n - 1L), size = as.double(m) * n, replace = TRUE,
prob = rep.int(c(n - 1, 1), c(1L, n - 1L)))
dim(A) <- c(m, n)
A
}
set.seed(1L)
A <- mkA(4L, 4L)
A
## [,1] [,2] [,3] [,4]
## [1,] 0 0 2 3
## [2,] 0 1 0 0
## [3,] 2 1 0 3
## [4,] 1 2 0 0
identical(f(A), gA <- g(A))
## [1] TRUE
gA[[1L, 1L]] # used for 'tcrossprod(A)[1L, 1L]'
## [1] 2 3
gA[[1L, 2L]] # used for 'tcrossprod(A)[1L, 2L]'
## [1] 2 3 1
gA[[2L, 1L]] # used for 'tcrossprod(A)[2L, 1L]'
## [1] 2 3 1
gA # under the hood, an 'm*(m+1)/2'-length list
## [[1]]
## [1] 2 3
##
## [[2]]
## [1] 2 3 1
##
## [[3]]
## [1] 1
##
## [[4]]
## [1] 2 3 1
##
## [[5]]
## [1] 1 2 3
##
## [[6]]
## [1] 2 1 3
##
## [[7]]
## [1] 2 3 1
##
## [[8]]
## [1] 1 2
##
## [[9]]
## [1] 2 1 3
##
## [[10]]
## [1] 1 2
##
## attr(,"class")
## [1] "utri"
mgA <- as.matrix(gA) # the full, symmetric, 'm'-by-'m' list matrix
mgA
## [,1] [,2] [,3] [,4]
## [1,] integer,2 integer,3 integer,3 integer,3
## [2,] integer,3 1 integer,3 integer,2
## [3,] integer,3 integer,3 integer,3 integer,3
## [4,] integer,3 integer,2 integer,3 integer,2
mgA[1L, ] # used for first row of 'tcrossprod(A)'
## [[1]]
## [1] 2 3
##
## [[2]]
## [1] 2 3 1
##
## [[3]]
## [1] 2 3 1
##
## [[4]]
## [1] 2 3 1
## If you need names
dimnames(mgA) <- rep.int(list(sprintf("A%d", seq_len(nrow(mgA)))), 2L)
mgA["A1", ]
## $A1
## [1] 2 3
##
## $A2
## [1] 2 3 1
##
## $A3
## [1] 2 3 1
##
## $A4
## [1] 2 3 1
## If you need an 'm'-by-'m' 'data.table' result
DT <- data.table::as.data.table(mgA)
DT
## A1 A2 A3 A4
## 1: 2,3 2,3,1 2,3,1 2,3,1
## 2: 2,3,1 1 1,2,3 1,2
## 3: 2,3,1 1,2,3 2,1,3 2,1,3
## 4: 2,3,1 1,2 2,1,3 1,2
这里有两个关于两个大整数矩阵的基准测试，表明预处理可以提供很大帮助：
set.seed(1L)
A <- mkA(100L, 1e+04L)
microbenchmark::microbenchmark(f(A), g(A), times = 10L, setup = gc(FALSE))
## Unit: milliseconds
## expr min lq mean median uq max neval
## f(A) 2352.0572 2383.3100 2435.7954 2403.8968 2431.6214 2619.553 10
## g(A) 843.0206 852.5757 858.7262 858.2746 863.8239 881.450 10
A <- mkA(100L, 1e+06L)
microbenchmark::microbenchmark(f(A), g(A), times = 10L, setup = gc(FALSE))
## Unit: seconds
## expr min lq mean median uq max neval
## f(A) 290.93327 295.54319 302.57001 301.17810 307.50226 318.14203 10
## g(A) 72.85608 73.83614 76.67941 76.57313 77.78056 83.73388 10
也许我们可以试试这个
f <- function(A, models) {
AA <- replace(A, A == 0, NA)
setNames(
c(t(outer(
1:nrow(A),
1:nrow(A),
Vectorize(function(x, y) unique(na.omit(c(t(AA[c(x, y)])))))
))),
t(outer(models, models, paste))
)
}
这给出了
$`A1 A1`
[1] 3 8
$`A1 A2`
[1] 3 8 5 7
$`A1 A3`
[1] 3 8 5
$`A2 A1`
[1] 5 7 3 8
$`A2 A2`
[1] 5 7
$`A2 A3`
[1] 5 7 3 8
$`A3 A1`
[1] 3 5 8
$`A3 A2`
[1] 3 5 8 7
$`A3 A3`
[1] 3 5 8
如果你在意速度，可以试试
lst <- asplit(replace(A, A == 0, NA), 1)
mat <- matrix(list(), nrow = nrow(A), ncol = nrow(A))
mat[lower.tri(mat)] <- combn(lst, 2, function(...) unique(na.omit(unlist(...))), simplify = FALSE)
mat[upper.tri(mat)] <- t(mat)[upper.tri(mat)]
diag(mat) <- Map(function(x) unname(x)[!is.na(x)], lst)
L <- c(t(mat))