Sparse matrix from vector

user_5055647 0
Viewed 123

I have a vector with values (val) and a vector indicating group membership (group):

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

Let's say we have K groups and a total of N values, hence both vectors have length N. The goal is to efficiently construct a sparse 'block-diagonal' matrix where the first column holds the values of group 1, the second column holds the values of group 2, and so on. However, the values should not 'overlap' in the sense that there should be only one value per row, see the solution below. I need to do this thousands of times with K and N very large. Hence, the following, loop-based solution is not efficient enough:

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

For memory reasons, it would be ideal to directly construct a sparse matrix without the intermediate step based on the dense N times K matrix.

EDIT

For the small example given above, it turns out that the loop based solution is quite efficient:

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

However, when turning to high-dimensional examples (N = 10,000 and K = 1,000), the solution by CJR is the winner speed-wise:

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
sparse-matrix r matrix
3 Answers

Matrix::.bdiag() will let you construct a block-diagonal (sparse) matrix directly from a list of matrices:

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

.bdiag(mm) is approximately equivalent to do.call(Matrix::bdiag, mm): ?bdiag says

The value of ‘bdiag()’ inherits from class ‘CsparseMatrix’, whereas ‘.bdiag()’ returns a ‘TsparseMatrix’.

(the former is a sorted compressed column-oriented form, the latter is in triplet form: ?"TsparseMatrix-class" says 'once [a triplet-oriented matrix] is created, however, the matrix is generally coerced to a ‘CsparseMatrix’ for further operations.')

?bdiag also has a Note:

This function has been written and is efficient for the case of relatively few block matrices which are typically sparse themselves.

So, this solution will definitely be better than what you've got, but further improvements might be possible.

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

I would suggest building the row & column indices you need directly and then giving them to the sparse constructor.

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

This works when groups are not monotonic:

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

But when groups are monotonic, it can be optimized using rle (as noted in comments):

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
user_73409480

You can try the code below

> 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

Benchmark

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"
)

shows

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
user_121587570
Related