【问题标题】：Sparse matrix from vector来自向量的稀疏矩阵
【发布时间】：2022-01-06 07:46:52
【问题描述】：

我有一个带有值的向量 (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() 如此糟糕真的很有趣，但我想这就是文档中的注释警告的内容。大概在 very 大尺寸下，它避免分配密集矩阵的事实至少会使其比原始矩阵更好......

标签： r matrix sparse-matrix

【解决方案1】：

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还有一个备注：

这个函数已经写好了，并且对于相对较少的块矩阵本身通常是稀疏的情况是有效的。

因此，此解决方案肯定会比您现有的解决方案更好，但可能会进一步改进。

【讨论】：

感谢您的深刻回答！

【解决方案2】：

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 进行优化（如 cmets 中所述）：

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 ?
@BenBolker 如果组是单调递增的。使用 split 不需要这样做（但如果它们严格单调增加rle 可能是更好的解决方案）。
感谢您的回答，我接受了它，因为您的解决方案在速度方面的扩展性非常好，请参阅我在上面添加的基准！
另一个简短的评论：使用rle(group)$lengths在高维情况下可以进一步加速超过50%，非常好！
这是有道理的，如果您可以确保您的组是单调的，那么这是一个很好的优化。

【解决方案3】：

你可以试试下面的代码

> 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

【讨论】：