Row-wise ECDF on SVT_SparseMatrix

[rcastelo]

Hi,

I'm interested in calculating the values of a so-called empirical cumulative distribution function (ECDF) for every row of a SVT_SparseMatrix object. This kind of ECDF calculation is used in the GSVA package in the first step of the GSVA algorithm. The current version of GSVA implements this calculation for dense base R matrices and sparse dgCMatrix objects. When the input is a sparse dgCMatrix, GSVA can calculate two different types of outputs: a dense output, as if one would coerce the input to a dense matrix and calculate the ECDF for every value of every row, or a sparse output, where GSVA calculates the ECDF for the nonzero values (only) of every row. It would be great if the SparseArray package would offer those two kinds of calculations for SVT_SparseMatrix objects.

ECDF values are calculated in two steps: (1) build the ECDF; and (2) evaluate the ECDF at the values of interest. I'm interested in evaluating the ECDF at the same values that I use to build it. The ECDF is calculated by first sorting the numerical values, then tabulating their occurrences, and finally calculating from the lowest to the highest value the running fraction of values smaller or equal than each value. The base R ecdf() function does this, e.g.:

x <- c(4, 3, 2, 2)
f <- ecdf(x)
f
Empirical CDF 
Call: ecdf(x)
 x[1:3] =      2,      3,      4
f(3)
[1] 0.75
f(x)
f(x)
[1] 1.00 0.75 0.50 0.50

In this toy example, the 75% of values are equal or below 3, 100% of values are equal or below 4, and 50% of values are equal or below 2. Instead of using the base R ecdf() function, we can also do this calculation as follows:

vals <- sort(unique(x))
mt <- match(x, vals)
tab <- tabulate(mt, nbins=length(vals))
ecdfvals <- cumsum(tab) / length(x)
ecdfvals[mt]
[1] 1.00 0.75 0.50 0.50

We can do this calculation for every row of a dense matrix with the following function:

## using base R 'ecdf()' function
## input and output are a dense matrix
ecdfvals_dense_to_dense <- function(x) t(apply(x, 1, function(rx) ecdf(rx)(rx)))

If the input matrix is a sparse dgCMatrix and the output is dense, this can be implemented as follows:

## input a sparse dgCMatrix, output a dense matrix
## ecdf values are calculated on both zero and nonzero values
ecdfvals_sparse_to_dense <- function(x) {
    stopifnot(is(x, "dgCMatrix")) ## QC
    res <- matrix(0, nrow(x), ncol(x))
    nc <- ncol(x)
    for (i in 1:nrow(x)) {
        rx <- x[i, , drop=FALSE]   ## drop=FALSE to keep w/ dgCMatrix
        vals <- unique(rx@x)       ## fetch unique nonzero values
        if (nnzero(rx) < nc)       ## if there is at least one zero
            vals <- c(0, vals)     ## add it to the unique values
        vals <- sort(vals)         ## sort unique values
        mt <- match(rx@x, vals)    ## tabulate unique nonzero values
        tab <- tabulate(mt, nbins=length(vals))
        whz <- which(tab == 0)     ## add counts of zero values
        if (nnzero(rx) < nc)       ## if present
            tab[whz] <- nc - nnzero(rx)
        ecdfvals <- cumsum(tab)/nc ## calculate ECDF values and put
        if (nnzero(rx) < nc)       ## them where they belong to
            res[i, ] <- rep(ecdfvals[whz], nc)
        res[i, which(diff(rx@p) > 0)] <- ecdfvals[mt]
        res
    }
    res
}

The following is a comparison of the performance of these two functions, which give the same result, and their current C implementations in the GSVA package:

library(Matrix)
library(GSVA)

nr <- 500
nc <- 1000
p <- nr*nc
fnz <- 0.05

x <- numeric(p)
x[sample(1:p, size=ceiling(fnz*p), replace=FALSE)] <- rnorm(ceiling(fnz*p))
m <- matrix(x, nrow=nr, ncol=nc)
M <- Matrix(m, sparse=TRUE)

bench::mark(d2d=ecdfvals_dense_to_dense(m),
            d2d_c=GSVA:::.ecdfvals_dense_to_dense(m, verbose=FALSE),
            s2d=ecdfvals_sparse_to_dense(M),
            s2d_c=GSVA:::.ecdfvals_sparse_to_dense(M, verbose=FALSE))
# A tibble: 4 × 13
  expression      min   median `itr/sec` mem_alloc `gc/sec` n_itr  n_gc
  <bch:expr> <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl> <int> <dbl>
1 d2d         79.47ms  87.37ms     11.2    40.13MB     5.60     6     3
2 d2d_c       19.07ms  24.56ms     38.0    25.04MB    13.3     20     7
3 s2d        783.16ms 783.16ms      1.28     1.9GB    35.8      1    28
4 s2d_c        4.55ms   5.83ms    133.      5.41MB     9.90    67     5
# ℹ 5 more variables: total_time <bch:tm>, result <list>, memory <list>,
#   time <list>, gc <list>
Warning message:
Some expressions had a GC in every iteration; so filtering is disabled. 

While in the dense-input and dense-output ("dense_to_dense") case, the C implementation is "only" about four times faster and consumes about 40% less memory, in the sparse-input and dense-output ("sparse_to_dense") case, the C implementation is two orders of magnitude more performant in terms of both speed and memory consumption.

An R implementation for the sparse-input and sparse-output could be the following:

## input and output are a sparse dgCMatrix
## ecdf values are calculated on nonzero values only
ecdfvals_sparse_to_sparse <- function(x) {
    stopifnot(is(x, "dgCMatrix")) ## QC
    nc <- ncol(x)
    for (i in 1:nrow(x)) {
        rx <- x[i, , drop=FALSE]   ## drop=FALSE to keep w/ dgCMatrix
        vals <- unique(rx@x)       ## fetch unique nonzero values
        vals <- sort(vals)         ## sort unique values
        mt <- match(rx@x, vals)    ## tabulate unique nonzero values
        tab <- tabulate(mt, nbins=length(vals))
                                   ## calculate ECDF values and put
        ecdfvals <- cumsum(tab)/nnzero(rx)
                                   ## them where they belong to
        x[i, which(diff(rx@p) > 0)] <- ecdfvals[mt]
    }
    x
}

The benchmarking of this R implementation and the current C implementation in GSVA is as follows:

bench::mark(s2s=ecdfvals_sparse_to_sparse(M),
            s2s_c=GSVA:::.ecdfvals_sparse_to_sparse(M, verbose=FALSE))
# A tibble: 2 × 13
  expression      min   median `itr/sec` mem_alloc `gc/sec` n_itr  n_gc
  <bch:expr> <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl> <int> <dbl>
1 s2s        282.35ms 500.34ms      2.00  174.68MB     5.00     2     5
2 s2s_c        3.45ms   4.04ms    216.      1.83MB     3.97   109     2
# ℹ 5 more variables: total_time <bch:tm>, result <list>, memory <list>,
#   time <list>, gc <list>
Warning message:
Some expressions had a GC in every iteration; so filtering is disabled. 

So, again the C implementation is two orders of magnitude more performant than the R implementation in both speed and memory consumption. In the C implementation, when the input is a sparse dgCMatrix, a RsparseMatrix copy of the input is created first, and then used to efficiently go through the rows.

If you consider that such calculations for SVT_SparseMatrix objects don't belong to this package, I'd appreciate if you could sketch me how could I implement this efficiently for SVT_SparseMatrix objects.

[vjcitn]

Hi @rcastelo is this still an ongoing concern?

[rcastelo]

Hi @vjcitn thanks for asking, I was about to write again here now that the release is over. In the last devel cycle I have provided in the GSVA package a C implementation of this calculation for SVT_SparseMatrix objects, and it would be great if @hpages or you would have some hint or advice on how to improve the performance of that implementation.

Let me explain my implementation using R code, since the implementation in the GSVA package is simply that R code re-implemented in C. First, I have implemented a function to fetch the nonzero values for a given (1-based) row of an SVT_SparseArray matrix object:

fetch_row_nzvals <- function(svt, i) {
    i_offsets <- lapply(svt@SVT, "[[", 2)
    i_offsets_mask <- i-1 == unlist(i_offsets) ## SVT offsets are 0-based
    i_offsets_mask <- relist(i_offsets_mask, i_offsets)
    i_nzvalues <- unlist(mapply(function(v, j) v[[1]][j], svt@SVT, i_offsets_mask))
    i_cols <- which(sapply(i_offsets_mask, any))
    list(x=i_nzvalues, j=i_cols)
}

Using this function, this would be the implementation of the ECDF calculation for an input SVT_SparseMatrix object and an output dense matrix:

## input a sparse SVT_SparseMatrix, output a dense matrix
## ecdf values are calculated on both zero and nonzero values
ecdfvals_svt_to_dense <- function(x) {
    stopifnot(is(x, "SVT_SparseMatrix")) ## QC
    res <- matrix(0, nrow(x), ncol(x))
    nc <- ncol(x)
    for (i in 1:nrow(x)) {
        rx <- fetch_row_nzvals(x, i)
        vals <- unique(rx$x)       ## fetch unique nonzero values
        if (length(rx$x) < nc)     ## if there is at least one zero
            vals <- c(0, vals)     ## add it to the unique values
        vals <- sort(vals)         ## sort unique values
        mt <- match(rx$x, vals)    ## tabulate unique nonzero values
        tab <- tabulate(mt, nbins=length(vals))
        whz <- which(tab == 0)     ## add counts of zero values
        if (length(rx$x) < nc)     ## if present
            tab[whz] <- nc - length(rx$x)
        ecdfvals <- cumsum(tab)/nc ## calculate ECDF values and put
        if (length(rx$x) < nc)     ## them where they belong to
            res[i, ] <- rep(ecdfvals[whz], nc)
        res[i, rx$j] <- ecdfvals[mt]
        res
    }
    res
}

This is the benchmarking of the this R implementation compared with the currently available C implementation of this same function in the GSVA package, as well as the C implementation of the function with a dgCMatrix as input:

S <- as(M, "SVT_SparseArray")

bench::mark(dgc2d_c=GSVA:::.ecdfvals_sparse_to_dense(M, verbose=FALSE),
            svt2d_c=GSVA:::.ecdfvals_svt_to_dense(S, verbose=FALSE),
            svt2d=ecdfvals_svt_to_dense(S))
# A tibble: 3 × 13
  <bch:expr> <bch:tm> <bch:t>     <dbl> <bch:byt>    <dbl> <int> <dbl>   <bch:tm> <list>   <list>     <list>     <list>
1 dgc2d_c      4.86ms  4.93ms    21.0      5.42MB     2.74    23     3      1.09s <dbl[…]> <Rprofmem> <bench_tm> <tibble>
2 svt2d_c     21.57ms 22.16ms    42.5      5.07MB     0       22     0   518.22ms <dbl[…]> <Rprofmem> <bench_tm> <tibble>
3 svt2d         3.16s   3.16s     0.316  168.44MB     2.21     1     7      3.16s <dbl[…]> <Rprofmem> <bench_tm> <tibble>
Warning message:
Some expressions had a GC in every iteration; so filtering is disabled.

So, not surprisingly, the C implementation with an SVT_SparseMatrix input is about two orders of magnitude faster than the R implementation, but it is about 5 times slower than the C implementation for dgCMatrix objects, which makes first a copy of the input dgCMatrix into a (row-wise) dgRMatrix that of course it is what speeds up the calculation.

The R implementation when both input and output are sparse as SVT_SparseMatrix objects is as follows:

## ecdf values are calculated on nonzero values only
ecdfvals_svt_to_svt <- function(x) {
    stopifnot(is(x, "SVT_SparseMatrix")) ## QC
    nc <- ncol(x)
    for (i in 1:nrow(x)) {
        rx <- fetch_row_nzvals(x, i)
        vals <- unique(rx$x)       ## fetch unique nonzero values
        vals <- sort(vals)         ## sort unique values
        mt <- match(rx$x, vals)    ## tabulate unique nonzero values
        tab <- tabulate(mt, nbins=length(vals))
                                   ## calculate ECDF values and put
        ecdfvals <- cumsum(tab)/length(rx$x)
                                   ## them where they belong to
        x[cbind(i, rx$j)] <- ecdfvals[mt]
    }
    x
}

This is the benchmarking between this R implementation and the corresponding C implementation in the GSVA package for both inputs, dgCMatrix and SVT_SparseArray matrix objects:

bench::mark(dgc2dgc_c=GSVA:::.ecdfvals_sparse_to_sparse(M, verbose=FALSE),
            svt2svt_c=GSVA:::.ecdfvals_svt_to_svt(S, verbose=FALSE),
            svt2svt=ecdfvals_svt_to_svt(S),
            check=FALSE) ## necessary since dgCMatrix and SVT_SparseArray results are different
# A tibble: 3 × 13                                                                                                       
  expression      min   median `itr/sec` mem_alloc `gc/sec` n_itr  n_gc total_time result memory     time       gc       
  <bch:expr> <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl> <int> <dbl>   <bch:tm> <list> <list>     <list>     <list>   
1 dgc2dgc_c    3.84ms   4.23ms   207.       1.83MB     3.98   104     2   502.12ms <NULL> <Rprofmem> <bench_tm> <tibble> 
2 svt2svt_c   19.72ms  20.76ms    44.3       1.5MB     0       23     0   518.65ms <NULL> <Rprofmem> <bench_tm> <tibble> 
3 svt2svt       3.71s    3.71s     0.269  169.37MB     2.15     1     8      3.71s <NULL> <Rprofmem> <bench_tm> <tibble> 
Warning message:                                                                                                         
Some expressions had a GC in every iteration; so filtering is disabled.                                                  

Again, the C implementation is two orders of magnitude faster than the R implementation, but 5 times slower than the C implementation when the input is a dgCMatrix object.

I would be interested in your advice on how to make the implementation for SVT_SparseMatrix objects equally performant, or even faster, than the implementation for dgCMatrix objects.

Another avenue of improvement might arise from reducing the memory consumption by "delaying" the division step, i.e., the dividing by the number of elements which is a common denominator for the whole row, by keeping only the tallies of the different values in a integer matrix with an additional integer vector with the denominator value for each row, as a row metadata column, so that the division could be done later on in the algorithm when needed.

[rcastelo]

Hi @hpages , following the discussion we had in issue #22 and here, I have been able to improve the performance of two different tasks that require scanning the rows of a SVT_SparseMatrix object:

1. Calculating the ranges of values and ranges of nonzero values for each row.
2. Calculating the ECDF values (see above for definition) for each row.

A fundamental difference between these two tasks is that the first one amounts to calculate the minimum and maximum values for each row, and following the current implementations of rowMins() and rowMaxs() as you explained in issue #22, you can implement it quite easily by walking on the columns, since you only need to keep track of a vector of minimum (maximum) values, and update it as you progress to the next column.

On the other hand, the second task requires identifying the unique values of a row, tally their occurrences and put them where the corresponding value was located. I guess there could be a way to carry on for each row the occurrences of unique values, and where do they occur, and be able to walk on the columns, but I thought this would require more time than I wanted to invest on it at the moment, and found the following way to improve row-wise operations like this.

My C implementation of the function that fetches the nonzero values and their column indices from a given row, , i.e. the C version of the fetch_row_nzvals() function above, walks on the columns, but because for each column it scans the vector of offsets until it finds the sought row, it leads to the observation you made by which its performance degrades with the rank of the row. Now, here is the trick, I'm not using this function to fetch arbitrary rows, but so far I'm only using it to fetch all the rows from the first one to the last one of the matrix. This means that using a vector with as many values as columns, initialized to zeroes, I can keep track of the last position of the sought offset for each column, updating it at each call of the function, and use it to speed up fetching the following offset. This leads to about a 30% improvement in performance for the task number 2.

For testing purposes, I have made a private R function that calls this C function with an optional third parameter corresponding to this vector. The following code compares using this function with and without this vector that tells which positions had the sought offset in the previous row. By the way, the following code uses the last devel version of GSVA 2.5.13, and if anyone reads this in the future, some of these private functions I made for testing purposes, may disappear.

library(SparseArray)
library(GSVA)

## coerce sparse to dense by fetching all rows independently
indrows <- function(X) {
    res <- matrix(0, nrow=nrow(X), ncol=ncol(X))
    for (i in seq_len(nrow(X))) {
        r <- GSVA:::.fetch_row_nzvals(X, i, NULL)
        res[i, r[[2]]+1] <- r[[1]] ## r[[2]] are 0-based indices of nonzero values 
    }
    res
} 

## coerce sparse to dense by fetching all rows with information
## dependent on the previous row
deprows <- function(X) {
    res <- matrix(0, nrow=nrow(X), ncol=ncol(X))
    whimin1 <- integer(ncol(X))
    for (i in seq_len(nrow(X))) {
        r <- GSVA:::.fetch_row_nzvals(X, i, whimin1)
        res[i, r[[2]]+1] <- r[[1]] ## r[[2]] are 0-based indices of nonzero values
        whimin1 <- r[[3]]
    }
    res
}

nr <- 1000
nc <- 10000
p <- nr*nc
fnz <- 0.05 ## fraction of nonzero values
x <- numeric(p)
x[sample(1:p, size=ceiling(fnz*p), replace=FALSE)] <- rnorm(ceiling(fnz*p))
m <- matrix(x, nrow=nr, ncol=nc)
S <- SparseArray(m)

bench::mark(indrows=indrows(S),
            deprows=deprows(S),
            native=as.matrix(S))
# A tibble: 3 × 13
  expression      min   median `itr/sec` mem_alloc `gc/sec` n_itr  n_gc total_time result   memory     time
  <bch:expr> <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl> <int> <dbl>   <bch:tm> <list>   <list>     <list>
1 indrows     708.9ms  708.9ms      1.41   202.5MB     1.41     1     1      709ms <dbl[…]> <Rprofmem> <bench_tm [1]>
2 deprows     432.1ms  437.6ms      2.29   202.5MB     1.14     2     1      875ms <dbl[…]> <Rprofmem> <bench_tm [2]>
3 native       26.8ms   27.3ms     31.0     76.3MB    11.6     16     6      516ms <dbl[…]> <Rprofmem> <bench_tm [16]>
# ℹ 1 more variable: gc <list>

As you can see, using the information to fetch the next row leads to about (708.9-432.1) / 708.9 ~ 39% performance improvement. I included the native way to coerce a sparse SVT_SparseMatrix object to a dense one just for comparison.

Using this improved implementation to fetch all the (rank-increasing) rows of a matrix, my initial naive implementation of calculating the (nonzero) ranges for each row also improved, but is still an order of magnitude slower than walking on the columns as we discussed.

library(Matrix)

M <- Matrix(m, sparse=TRUE)

bench::mark(dgCMtr_C=GSVA:::.rowNzRanges_dgCMatrix(M),
            SVTrow_C=GSVA:::.rowNzRanges_SVT_SparseArray_byrow(S),
            SVTtrp_C=GSVA:::.rowNzRanges_SVT_SparseArray_transpose_C(S),
            SVTrxc_R=GSVA:::.rowNzRanges_SVT_SparseArray_rowbycols_R(S),
            SVTrxc_C=GSVA:::.rowNzRanges_SVT_SparseArray(S))
# A tibble: 5 × 13
  expression      min   median `itr/sec` mem_alloc `gc/sec` n_itr  n_gc total_time result   memory     time
  <bch:expr> <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl> <int> <dbl>   <bch:tm> <list>   <list>     <list>
1 dgCMtr_C     6.19ms    6.3ms    149.      5.76MB     2.05    73     1      489ms <dbl[…]> <Rprofmem> <bench_tm [74]>
2 SVTrow_C   317.84ms 319.99ms      3.13    31.3KB     0        2     0      640ms <dbl[…]> <Rprofmem> <bench_tm [2]>
3 SVTtrp_C    10.44ms  10.53ms     94.4     5.88MB     2.05    46     1      487ms <dbl[…]> <Rprofmem> <bench_tm [47]>
4 SVTrxc_R     8.65ms   8.76ms    114.     177.3KB     0       57     0      500ms <dbl[…]> <Rprofmem> <bench_tm [57]>
5 SVTrxc_C     4.94ms   4.96ms    201.      31.3KB     0      101     0      502ms <dbl[…]> <Rprofmem> <bench_tm [101]>
# ℹ 1 more variable: gc <list>

Here I have included in the comparison the calculation when the input is a dgCMatrix, labeled dgCMtr_C, in which case the function first coerces it in to a dgRMatrix, which is also one of the strategies I first implemented for SVT_SparseMatrix object, here labeled SVTtrp_C, and as expected both have a fairly good performance but with the largest memory consumption. Then the row-wise strategy I described above, labeled SVTrow_C, is the slowest one with a minimal memory consumption, and finally the walking-on-the-columns strategy, implemented in both R (label SVTrxc_R) and C (label SVTrxc_C), are the fastest, but the C implementation has obviously a smaller memory footprint.

So, finally, this is the performance for the ECDF values calculation (generating again the simulated data to compare it the performances reported previously in the conversation):

nr <- 500
nc <- 1000   
p <- nr*nc    
fnz <- 0.05    

x <- numeric(p)
x[sample(1:p, size=ceiling(fnz*p), replace=FALSE)] <- rnorm(ceiling(fnz*p))
m <- matrix(x, nrow=nr, ncol=nc)
M <- Matrix(m, sparse=TRUE)
S <- SparseArray(m)

bench::mark(dgc2d_c=GSVA:::.ecdfvals_sparse_to_dense(M, verbose=FALSE),
            svt2d_c=GSVA:::.ecdfvals_svt_to_dense(S, verbose=FALSE))
# A tibble: 2 × 13
  expression      min   median `itr/sec` mem_alloc `gc/sec` n_itr  n_gc total_time result   memory     time
  <bch:expr> <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl> <int> <dbl>   <bch:tm> <list>   <list>     <list>
1 dgc2d_c      4.44ms   4.61ms     210.     5.36MB     6.63    95     3      452ms <dbl[…]> <Rprofmem> <bench_tm [98]>
2 svt2d_c     13.71ms  13.83ms      72.1    5.07MB     0       37     0      513ms <dbl[…]> <Rprofmem> <bench_tm [37]>
# ℹ 1 more variable: gc <list>

bench::mark(dgc2dgc_c=GSVA:::.ecdfvals_sparse_to_sparse(M, verbose=FALSE),
            svt2svt_c=GSVA:::.ecdfvals_svt_to_svt(S, verbose=FALSE),
            check=FALSE) ## necessary since dgCMatrix and SVT_SparseArray results are different
# A tibble: 2 × 13
  expression      min  median `itr/sec` mem_alloc `gc/sec` n_itr  n_gc total_time result memory     time       gc
  <bch:expr> <bch:tm> <bch:t>     <dbl> <bch:byt>    <dbl> <int> <dbl>   <bch:tm> <list> <list>     <list>     <list>
1 dgc2dgc_c    3.69ms  3.83ms     261.     1.83MB     2.16   121     1      464ms <NULL> <Rprofmem> <bench_tm> <tibble>
2 svt2svt_c   12.51ms 12.61ms      79.1     1.5MB     0       40     0      506ms <NULL> <Rprofmem> <bench_tm> <tibble>

So, here I'm comparing, as I did previously, the implementation for dgCMatrix objects, which coerces it into a dgRMatrix to efficiently scan rows, with the one scanning the rows of a SVT_SparseMatrix object, but this last performance test is using the improved version of row scanning, and in fact one can see that in the previous implementation, the dgCMatrix based approach was 4.43 (sparse-to-dense) and 5.13 (sparse-to-sparse) times faster, but using the new row-scanning implementation this has reduced to 3.09 and 3.39, respectively, which is about a 30% improvement.

In summary, for the first task of computing the (nonzero) ranges per row, we have improved the implementation in more than one order of magnitude (317.84/4.94=64.34 times faster), and for the second task I have now an implementation that is about 30% faster, and while this cannot be seen in this toy benchmarking, with a substantial reduction in memory consumption with respect to the approach based on "transposing" a sparse matrix.

If you have no further thoughts about this, you may close the issue, thank you very much for your help!!