我有一个值向量(
假设我们有
val
)和一个表示组成员资格的向量(group
):vec <- 1:9
group <- rep(1:3, c(2,4,3))
假设我们有
K
个组和总共N
个值,因此两个向量的长度都为N
。目标是高效地构建一个稀疏的“块对角”矩阵,其中第一列持有第1组的值,第二列持有第2组的值,依此类推。然而,值不应该在行中“重叠”,也就是说,每行应该只有一个值,参见下面的解决方案。我需要用非常大的K
和N
数千次地执行这个操作。因此,以下基于循环的解决方案效率不够高:K <- length(unique(group))
N <- length(group)
M <- matrix(0, N, K)
for(k in 1:K){
M[group == k, k] <- vec[group == k]
}
Matrix::Matrix(M, sparse = T)
9 x 3 sparse Matrix of class "dgCMatrix"
[1,] 1 . .
[2,] 2 . .
[3,] . 3 .
[4,] . 4 .
[5,] . 5 .
[6,] . 6 .
[7,] . . 7
[8,] . . 8
[9,] . . 9
出于内存原因,最好直接构建稀疏矩阵,而不是基于密集的N乘K矩阵进行中间步骤。
编辑
对于上面给出的小例子,循环解决方案非常高效:
Unit: microseconds
expr min lq mean median uq max neval cld
ben 734.280 771.7000 826.8372 787.5230 805.2710 3185.158 100 b
CJR 711.187 745.1855 813.9948 766.9960 781.7495 4832.476 100 b
original 199.714 221.9520 235.4320 227.9395 236.7065 379.757 100 a
然而,当处理高维度的实例时(N = 10,000,K = 1,000),CJR 的解决方案在速度方面是最优的:
Unit: milliseconds
expr min lq mean median uq max neval cld
ben 128.529311 133.308972 140.032070 135.921289 139.272589 289.668852 100 b
CJR 1.841474 2.055513 2.261732 2.201557 2.395925 6.330544 100 a
original 93.387806 118.348398 171.380301 125.884493 244.421699 365.871433 100 c
Matrix::.bdiag()
如此糟糕,但我想这就是文档中警告的内容。可能在非常大的尺寸下,它避免分配密集矩阵的事实至少会使它比原始版本更好... - Ben Bolker