logo标签列表

稀疏矩阵的向量表示

rmatrixsparse-matrix
4
我有一个值向量（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
3个回答
5

Matrix::.bdiag()允许您直接从一组矩阵构造一个块对角线（稀疏）矩阵：

mm <- lapply(split(vec, group), matrix)
Matrix::.bdiag(mm)

.bdiag(mm) 大致相当于 do.call(Matrix::bdiag, mm): ?bdiag 说明：

‘bdiag()’ 的返回值继承自 ‘CsparseMatrix’ 类，而 ‘.bdiag()’ 返回一个 ‘TsparseMatrix’。

(前者是有序的压缩列向量形式，后者则是三元组形式：?"TsparseMatrix-class" 说 '一旦[创建了一个三元组矩阵]，但通常会将矩阵强制转换为 ‘CsparseMatrix’ 进行进一步操作。')

?bdiag 还有一个 注:

这个函数是针对数量相对较少的块矩阵的情况编写的，并且这些块矩阵本身通常是稀疏的。

所以，这个解决方案肯定比你现有的更好，但可能还有改进的余地。

感谢您提供这个富有洞见的答案！ - yrx1702
4
vec   <- 1:9
group <- rep(1:3, c(2,4,3))

我建议直接构建所需的行和列索引，然后将它们提供给稀疏矩阵的构造函数。
i <- unlist(split(vec, group), use.names = F)
j <- vapply(split(vec, group), length, numeric(1))
Matrix::sparseMatrix(i=i,
                     j=rep(1:length(j), j),
                     x=vec[i])

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

当组不是单调的时候，这个方法仍然有效：

vec   <- 1:9
group <- c(5:1, 2:5)

9 x 5 sparse Matrix of class "dgCMatrix"
               
 [1,] . . . . 1
 [2,] . . . 2 .
 [3,] . . 3 . .
 [4,] . 4 . . .
 [5,] 5 . . . .
 [6,] . 6 . . .
 [7,] . . 7 . .
 [8,] . . . 8 .
 [9,] . . . . 9

然而，当组是单调的时候，可以使用rle进行优化（如评论中所述）：

vec   <- 1:9
group <- rep(1:3, c(2,4,3))

j <- rle(group)$length
Matrix::sparseMatrix(i=1:length(group),
                     j=rep(1:length(j), j),
                     x=vec)

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
我认为 j 也是 rle(group)$lengths 吗？ - Ben Bolker
如果组是单调递增的话，使用split就可以了。这并不需要严格要求它们是单调递增的（但如果它们是严格单调递增的，rle可能是更好的解决方案）。 - CJR
谢谢你的回答，我接受了它，因为你的解决方案在速度方面表现得非常好，可以看到我上面添加的基准测试！ - yrx1702
另一个简短的评论：在高维情况下，使用rle（group）$lengths可以进一步加快速度超过50％，非常好！ - yrx1702
如果您可以确保您的组是单调的，那么这是一个很好的优化，这是有道理的。 - CJR
2
你可以尝试下面的代码。
> Matrix(`[<-`(M, cbind(seq_along(group), group), vec))
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

基准测试

microbenchmark(
  ben = {
    mm <- lapply(split(vec, group), matrix)
    Matrix::.bdiag(mm)
  },
  CJR = {
    i <- unlist(split(vec, group), use.names = F)
    j <- vapply(split(vec, group), length, numeric(1))
    Matrix::sparseMatrix(
      i = i,
      j = rep(1:length(j), j),
      x = vec[i]
    )
  },
  TIC = {
    Matrix(`[<-`(M, cbind(seq_along(group), group), vec))
  },
  check = "equivalent"
)

显示
Unit: microseconds
 expr   min     lq    mean median    uq    max neval
  ben 564.0 599.55 662.640 654.25 686.9 1213.5   100
  CJR 523.1 564.70 643.524 619.65 675.0 1222.6   100
  TIC 165.5 191.90 217.537 208.00 234.7  520.1   100
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接