Performance impact with sparse matrix

[yoda-vid]

First off, thanks for this package, which has saved me countless hours of debugging in R.

I've encountered a perplexing issue where accessing sparse matrices appears to be slower or even to hang within a tryCatchLog block. I can reproduce it by generating a somewhat large sparse matrix with column names, which prints fine on its own but appears to hang when printing it within a tryCatchLog block.

# generate sparse matrix with columns
mat <- matrix(rnorm(1e6), ncol=10)
mat.sparse <- Matrix::Matrix(mat, sparse=TRUE)
colnames(mat.sparse) <- rep("a", ncol(mat.sparse))

# print matrix, which completes quickly
print(mat.sparse)

# print matrix in tryCatchLog, which hangs for me
tryCatchLog::tryCatchLog(print(mat.sparse))

For some reason, I only see this issue when including column names (the situation where I first encountered this). Even stranger is if I load a library such as Seurat (where I discovered this issue), the problem is only apparent when printing the matrix within a list. Here is a comparison with and without tryCatchLog for matrices of increasing size, where the time differences become more dramatic as size increases.

testSize <- function(size) {
  mat <- matrix(rnorm(size), ncol=10)
  mat.sparse <- Matrix::Matrix(mat, sparse=TRUE)
  colnames(mat.sparse) <- rep("a", ncol(mat.sparse))
  
  message("Running outside of tryCatchLog")
  start.time <- Sys.time()
  print(mat.sparse)
  times <- list(base=difftime(Sys.time(), start.time, units="mins"))
  
  tryCatchLog::tryCatchLog({
    message("Running inside tryCatchLog")
    start.time <- Sys.time()
    print(mat.sparse)
    times[["try"]] <- difftime(Sys.time(), start.time, units="mins")
  })
  
  return(times)
}

elapsed <- list()
for (n in c(1e3, 1e4, 1e5, 2e5, 3e5, 4e5, 5e5)) {
  elapsed[[paste0("Size_", n)]] <- testSize(n)
}
times <- matrix(unlist(elapsed), nrow=2)
dimnames(times) <- list(c("Baseline", "tryCatchLog"), names(elapsed))
print(times)

Outputs for increasing matrix sizes (in minutes):

               Size_1000   Size_10000 Size_1e+05 Size_2e+05 Size_3e+05
Baseline    0.0001921495 0.0009166161 0.02194637 0.02247133 0.02464442
tryCatchLog 0.0004020850 0.0029591680 0.13794237 0.39542582 0.80269897
            Size_4e+05 Size_5e+05
Baseline    0.02131515  0.0225001
tryCatchLog 1.36883435  2.2514157

Would this have something to do with serializing the matrix in the tryCatchLog block?

Thanks for your help and for this package!

sessionInfo() output:

R version 4.0.2 Patched (2020-08-10 r79000)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Catalina 10.15.7

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

loaded via a namespace (and not attached):
[1] compiler_4.0.2

[aryoda]

Thanks for reporting this issue together with a reproducible example and the benchmark code!

My first thought is that the established condition handlers could cost a lot of time but fired conditions may not always be visible (eg. custom condtions).

Which version of tryCatchLog did you use (the sessionInfo does not show it)?

[aryoda]

FYI:

I can reproduce the behaviour on the current release candidate of tryCatchLog!

If I call the underlying tryCatch of base R directly (tryCatch(print(mat.sparse))) there is no performance impact so something must happen inside of tryCatchLog. I can see a multiple messages that are thrown for some reason when using tryCatchLog:

INFO [2021-09-19 17:46:23]    [[ suppressing 10 column names ‘a’, ‘a’, ‘a’ ... ]]

Compact call stack:
  1 tryCatchLog::tryCatchLog(print(mat.sparse))
  2 tryCatchLog.R#353: tryCatch(withCallingHandlers(expr, condition = cond.handler), ..., finally = finally)
  3 tryCatchLog.R#353: withCallingHandlers(expr, condition = cond.handler)
  4 tryCatchLog.R#353: print(mat.sparse)

Full call stack:
  1 tryCatchLog::tryCatchLog(print(mat.sparse))
  2 tryCatchLog.R#353: tryCatch(withCallingHandlers(expr, condition = cond.handler), ..., finally = finally)
  3 tryCatchList(expr, classes, parentenv, handlers)
  4 tryCatchLog.R#353: withCallingHandlers(expr, condition = cond.handler)
  5 tryCatchLog.R#353: print(mat.sparse)
  6 print.default(mat.sparse)
  7 (new("standardGeneric", .Data = function (object) 
    standardGeneric("show"), generic = "show", package = "methods", group
  8 (new("standardGeneric", .Data = function (object) 
    standardGeneric("show"), generic = "show", package = "methods", group
  9 printSpMatrix2(object)
  10 formatSpMatrix(x[seq_len(nr1), , drop = FALSE], digits = digits, maxp = maxp, zero.print = zero.print, col.names = col.names
  11 .formatSparseSimple(m, asLogical = logi, digits = digits, col.names = col.names, note.dropping.colnames = note.dropping.coln
  12 emptyColnames(cx, msg.if.not.empty = note.dropping.colnames)
  13 message(if (lc > 3) gettextf("   [[ suppressing %d column names %s ... ]]", nc, paste(sQuote(cn[1:3]), collapse = ", ")) els
  14 withRestarts({
        signalCondition(cond)
        defaultHandler(cond)
    }, muffleMessage = function() NULL)
  15 withOneRestart(expr, restarts[[1]])
  16 doWithOneRestart(return(expr), restart)
  17 signalCondition(cond)



INFO [2021-09-19 17:46:25]    [[ suppressing 10 column names ‘a’, ‘a’, ‘a’ ... ]]


Compact call stack:
  1 tryCatchLog::tryCatchLog(print(mat.sparse))
  2 tryCatchLog.R#353: tryCatch(withCallingHandlers(expr, condition = cond.handler), ..., finally = finally)
  3 tryCatchLog.R#353: withCallingHandlers(expr, condition = cond.handler)
  4 tryCatchLog.R#353: print(mat.sparse)

Full call stack:
  1 tryCatchLog::tryCatchLog(print(mat.sparse))
  2 tryCatchLog.R#353: tryCatch(withCallingHandlers(expr, condition = cond.handler), ..., finally = finally)
  3 tryCatchList(expr, classes, parentenv, handlers)
  4 tryCatchLog.R#353: withCallingHandlers(expr, condition = cond.handler)
  5 tryCatchLog.R#353: print(mat.sparse)
  6 print.default(mat.sparse)
  7 (new("standardGeneric", .Data = function (object) 
    standardGeneric("show"), generic = "show", package = "methods", group = list(
  8 (new("standardGeneric", .Data = function (object) 
    standardGeneric("show"), generic = "show", package = "methods", group = list(
  9 printSpMatrix2(object)
  10 printSpMatrix(x[seq_len(nr1), , drop = FALSE], digits = digits, maxp = maxp, zero.print = zero.print, col.names = col.names, note.dr
  11 formatSpMatrix(x, digits = digits, maxp = maxp, cld = cld, zero.print = zero.print, col.names = col.names, note.dropping.colnames = 
  12 .formatSparseSimple(m, asLogical = logi, digits = digits, col.names = col.names, note.dropping.colnames = note.dropping.colnames, dn
  13 emptyColnames(cx, msg.if.not.empty = note.dropping.colnames)
  14 message(if (lc > 3) gettextf("   [[ suppressing %d column names %s ... ]]", nc, paste(sQuote(cn[1:3]), collapse = ", ")) else gettex
  15 withRestarts({
        signalCondition(cond)
        defaultHandler(cond)
    }, muffleMessage = function() NULL)
  16 withOneRestart(expr, restarts[[1]])
  17 doWithOneRestart(return(expr), restart)
  18 signalCondition(cond)



INFO [2021-09-19 17:46:27]    [[ suppressing 10 column names ‘a’, ‘a’, ‘a’ ... ]]


Compact call stack:
  1 tryCatchLog::tryCatchLog(print(mat.sparse))
  2 tryCatchLog.R#353: tryCatch(withCallingHandlers(expr, condition = cond.handler), ..., finally = finally)
  3 tryCatchLog.R#353: withCallingHandlers(expr, condition = cond.handler)
  4 tryCatchLog.R#353: print(mat.sparse)

Full call stack:
  1 tryCatchLog::tryCatchLog(print(mat.sparse))
  2 tryCatchLog.R#353: tryCatch(withCallingHandlers(expr, condition = cond.handler), ..., finally = finally)
  3 tryCatchList(expr, classes, parentenv, handlers)
  4 tryCatchLog.R#353: withCallingHandlers(expr, condition = cond.handler)
  5 tryCatchLog.R#353: print(mat.sparse)
  6 print.default(mat.sparse)
  7 (new("standardGeneric", .Data = function (object) 
    standardGeneric("show"), generic = "show", package = "methods", group = list(
  8 (new("standardGeneric", .Data = function (object) 
    standardGeneric("show"), generic = "show", package = "methods", group = list(
  9 printSpMatrix2(object)
  10 printSpMatrix(tail(x, nr2), digits = digits, maxp = maxp, zero.print = zero.print, col.names = col.names, note.dropping.colnames = n
  11 formatSpMatrix(x, digits = digits, maxp = maxp, cld = cld, zero.print = zero.print, col.names = col.names, note.dropping.colnames = 
  12 .formatSparseSimple(m, asLogical = logi, digits = digits, col.names = col.names, note.dropping.colnames = note.dropping.colnames, dn
  13 emptyColnames(cx, msg.if.not.empty = note.dropping.colnames)
  14 message(if (lc > 3) gettextf("   [[ suppressing %d column names %s ... ]]", nc, paste(sQuote(cn[1:3]), collapse = ", ")) else gettex
  15 withRestarts({
        signalCondition(cond)
        defaultHandler(cond)
    }, muffleMessage = function() NULL)
  16 withOneRestart(expr, restarts[[1]])
  17 doWithOneRestart(return(expr), restart)
  18 signalCondition(cond)

[aryoda]

First analysis:

The performance loss happens in the function limitedLabelsCompact when the call stack in value (a sys.calls object) is converted into a string via as.character (eg. here https://github.com/aryoda/tryCatchLog/blob/master/R/limited_Labels_Compact.R#L72):

Somehow the print function in call 7 and 8 of the above full call stack has a huge string length (nchar > 1 Mio)
which is then cut back to maxwidth via strtrim.

The goal is now to avoid creating the huge strings (just to cut them later) but this is done in C code inside of (Primitive("as.character"))...

Edit 1: The R-internal function limitedLabels seems to have the same problem (maybe worth an upstream bug report to base R)...

[yoda-vid]

Thank you for checking this and finding the source!

Sorry about the missing version number. Here is sessionInfo after loading everything, with the latest release of tryCatchLog from CRAN.

R version 4.0.2 Patched (2020-08-10 r79000)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Catalina 10.15.7

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

loaded via a namespace (and not attached):
[1] compiler_4.0.2       tryCatchLog_1.2.4    Matrix_1.3-4        
[4] formatR_1.11         futile.logger_1.4.3  lambda.r_1.2.4      
[7] futile.options_1.0.1 grid_4.0.2           lattice_0.20-44     

When I first ran into this issue while using Seurat, I initially got around it by quoteing the sparse matrix here, which avoided the hang during the do.call step on the list holding the matrix. It was later when I realized that I could only reproduce the original issue in tryCatchBlock that I looked here. I tried to quote the value in the function you mentioned but did not see improvement, though I do not fully understand the usage of quote and whether it would be applicable here.

Thanks for your help with this!

[aryoda]

Bug fix implementation details

Implementation see the new branch:

https://github.com/aryoda/tryCatchLog/tree/bug/68_performance_with_matrix

How to install the bug fix branch for testing

# install.packages("devtools")
library(devtools)
install_github("aryoda/tryCatchLog", ref = "bug/68_performance_with_matrix")

Minimal reproducible example

You can use this code snippet to reproduce the problem and set a breakpoint into the limitedLabelsCompact function:

# generate sparse matrix with columns
mat <- matrix(rnorm(1e3), ncol=10)
mat.sparse <- Matrix::Matrix(mat, sparse=TRUE)
colnames(mat.sparse) <- rep("a", ncol(mat.sparse))
tryCatchLog::tryCatchLog(print(mat.sparse))

Background

limitedLabelsCompact is passed the calls stack from sys.calls in the "value" argument.
sys.calls returns a (pair)list of all active calls and frames.
Each element of the list contains an (unevaluated) call of a function which consists of
of the function and the already evaluated (!) arguments ("data" like the matrix).

Each call element of the list in value is a language object with sub elements (like an AST):

Browse[2]> str(value)
List of 18
 $ : language tryCatchLog::tryCatchLog(print(mat.sparse))
 $ : language tryCatch(withCallingHandlers(expr, condition = cond.handler), ..., finally = finally)
  ..- attr(*, "srcref")= 'srcref' int [1:8] 353 3 365 3 3 3 1230 1242
  .. ..- attr(*, "srcfile")=Classes 'srcfilealias', 'srcfile' <environment: 0x560d846aaa20> 
 $ : language tryCatchList(expr, classes, parentenv, handlers)
 $ : language withCallingHandlers(expr, condition = cond.handler)
  ..- attr(*, "srcref")= 'srcref' int [1:8] 353 3 365 3 3 3 1230 1242
  .. ..- attr(*, "srcfile")=Classes 'srcfilealias', 'srcfile' <environment: 0x560d846aaa20> 
 $ : language print(mat.sparse)
  ..- attr(*, "srcref")= 'srcref' int [1:8] 353 3 365 3 3 3 1230 1242
  .. ..- attr(*, "srcfile")=Classes 'srcfilealias', 'srcfile' <environment: 0x560d846aaa20> 
 $ : language print.default(mat.sparse)
 $ : language (new("standardGeneric", .Data = function (object)  standardGeneric("show"), generic = "show", package = "methods"| __truncated__ ...
 $ : language (new("standardGeneric", .Data = function (object)  standardGeneric("show"), generic = "show", package = "methods"| __truncated__ ...
 $ : language printSpMatrix2(object)
 $ : language printSpMatrix(x, cld = cld, digits = digits, maxp = maxp, zero.print = zero.print, col.names = col.names, note.dr| __truncated__
 $ : language formatSpMatrix(x, digits = digits, maxp = maxp, cld = cld, zero.print = zero.print, col.names = col.names, note.d| __truncated__
 $ : language .formatSparseSimple(m, asLogical = logi, digits = digits, col.names = col.names, note.dropping.colnames = note.dr| __truncated__
 $ : language emptyColnames(cx, msg.if.not.empty = note.dropping.colnames)
 $ : language message(if (lc > 3) gettextf("   [[ suppressing %d column names %s ... ]]", nc, paste(sQuote(cn[1:3]), collapse =| __truncated__
 $ : language withRestarts({     signalCondition(cond) ...
 $ : language withOneRestart(expr, restarts[[1L]])
 $ : language doWithOneRestart(return(expr), restart)
 $ : language signalCondition(cond)

Now you can find the complete data within value at the 7th call stack element as 2nd sub element (= 1st argument) of the function call:

Browse[2]> value[[7]][[2]]
10000 x 10 sparse Matrix of class "dgCMatrix"
   [[ suppressing 10 column names ‘a’, ‘a’, ‘a’ ... ]]
   [[ suppressing 10 column names ‘a’, ‘a’, ‘a’ ... ]]
                                                                                                                                
 [1,] -0.27101375  1.27095210  0.27920829  0.113388502  0.683061498 -0.26678504  0.15789238 -1.37377665 -1.676342271  0.39381950
 [2,]  0.89699289 -1.21024923  2.31602159 -1.072096302 -0.186674139 -0.13325089  2.20758412  0.39885638  2.307467052  0.22826009
 [3,]  0.05281388  0.97670370 -1.09529982  0.004114273  0.477866129  1.17097650 -1.07072814 -0.34552467  0.806597933  0.63766978
 [4,]  0.12681865  1.38273777  2.35218714  1.043297314  0.582625477  2.53524418  1.13970371 -0.55565279 -2.037443232  0.62122915
 [5,] -0.08898476 -0.04354125 -0.21621703 -0.711896155  0.566070531 -0.07279835 -0.30158637 -0.34781138 -2.202115348 -0.15463500
 [6,] ...

Since limitedLabelsCompact tries to convert the matrix data into a full character string it takes a lot of time (and RAM memory):
The larger the matrix the more time and memory is consumed:

Browse[2]> object.size(value[[7]][[2]])
1201776 bytes
as.character(value[[7]][[2]])  # may take very long...

Implemented fix

To avoid this the conversion of the full value call stack into character must be avoid.

Luckily the implementation of as.character for language objects is implemented via the deparse function
which allows to limit the the output to a maximum number of lines which can then be converted into character.

Now it is easy to change the implementation to extract only the first line (compact = TRUE)
or the first n characters of each call (compact = FALSE with maxwidth = getOption("width") - 5L):

  if (compact == TRUE)
    max.rows <- 1 else
    max.rows <- 30    # nobody wants to see more than 30 rows for each call in the call stack ;-)

  value <- lapply(value, function(x) {
                            if (is.language(x))
                               paste(deparse(x, width.cutoff = 500, nlines = max.rows), collapse = "\n") else
                               as.character(x)
                            })

Open problems

The stack trace output is shorter now than before because the internal implementation in C uses a different
width.cutoff than the default of 60 characters.

See R's internal implementation:

# Max cutoff is bufsize - 12 => 500
https://github.com/wch/r-source/blob/trunk/src/main/deparse.c#L120
# SEXP attribute_hidden do_deparse(SEXP call, SEXP op, SEXP args, SEXP rho)
https://github.com/wch/r-source/blob/trunk/src/main/deparse.c#L165

Expected best solution:

With width.cutoff = 60 (R's default) nine unit tests fail!

I will hard code the internal width.cutoff (500) in tryCatchLog and hope it will never change in R internally
(otherwise the original as.character output of a call has a different character count than tryCatchLog's implementation
or may even fail if R throws an error instead of a warning internally in the future!).

There are "only" three failing unit tests in test_build_log_entry.R and test_build_log_output.R then
which may be caused by different indentions now.

Next steps

• I will further examine the output differences in a few days to estimate the severity...
• @yoda-vid It would be great if you could install the bug fix as described above to test if the perfomance issues are gone now. THX :-)

[yoda-vid]

This is really amazing! I confirm that the update fixes the issue on my end. I see no lag or hang with large sparse matrices now.

Here is the output from the original benchmark test on this new branch:

               Size_1000  Size_10000 Size_1e+05 Size_2e+05 Size_3e+05
Baseline    0.0002154509 0.001738803 0.04303877 0.04058193 0.04426195
tryCatchLog 0.0004967690 0.001989186 0.04322071 0.04099272 0.04334397
            Size_4e+05 Size_5e+05
Baseline    0.03875228  0.0414661
tryCatchLog 0.03971055  0.0429032

Times are now on par with and without tryCatchLog.

Thank you for the fix and the thorough explanation! I've learned a ton about R internals from your posts. Very clever to use deparse directly!

[aryoda]

I have analyzed the three failing unit tests now (all related to the call stack output).

Analysis results:

There are two reasons for the failing unit tests:

1. The number of lines per call is now limited to 30 (instead of unlimited before). So some code lines are missing. This is OK.
2. The deparsed code differs in some cases between direct deparse and as.character calls:

   $ diff test.txt expected_full_stack_trace.txt
   46c46
   <   5 tryCatchOne(expr, names, parentenv, handlers[[1L]])
   ---
   >   5 tryCatchOne(expr, names, parentenv, handlers[[1]])
   

Problem 1 is OK, reducing the number or code lines in the stack trace is intentional.

Problem 2 is due inconsistent default arguments in the implementation of as.character in R + C vs. deparse in R.

A direct deparse call in R uses control = c("keepNA", "keepInteger", "niceNames", "showAttributes") as default (see ?.deparseOpts for details).

The as.character implementation in the C code of base R calls the internal deparse C function with another default for .deparseOpts: The SIMPLEDEPARSE C constant which corresponds to control = NULL.

Solution

I will apply the following changes:

1. Add control = NULL to the deparse call in the function limitedLabelsCompact:

   paste(deparse(x, width.cutoff = 500, nlines = max.rows, control = NULL), collapse = "\n")
   

This is reasonable since ?.deparseOpts says:

For the most readable (but perhaps incomplete) display, use control = NULL.

2. Adjust the expected stack trace output in the three failing unit tests to the new line output limit (which I will further reduce to minimize the output size - R does print only one line per call in an error call stack!).

Details

The function call path of R's as.character into R's C code is:

> as.character
function (x, ...)  .Primitive("as.character")

1. #define AS_CHARACTER(x) Rf_coerceVector(x,STRSXP)
2. #define coerceVector Rf_coerceVector
3. SEXP coerceVector(SEXP v, SEXPTYPE type)
4. ans = coercePairList(v, type);
5. SET_STRING_ELT(rval, i, STRING_ELT(deparse1line(CAR(vp), 0), 0));
6. SEXP deparse1line(SEXP call, Rboolean abbrev) { return deparse1line_(call, abbrev, SIMPLEDEPARSE); }
7. SEXP deparse1line_(SEXP call, Rboolean abbrev, int opts)

The SIMPLEDEPARSE option is defined as follows in C:

#define KEEPINTEGER 		1
...
#define KEEPNA			64
...
#define NICE_NAMES             	1024
/* common combinations of the above */
#define SIMPLEDEPARSE		0
#define DEFAULTDEPARSE		1089 /* KEEPINTEGER | KEEPNA | NICE_NAMES, used for calls */
...

Compare: The stack trace of an error is printed with different deparse options and only one line per call!

Call path:

./errors.c: PrintWarnings() | warningcall_dflt() ... -> deparse1s() -> deparse1WithCutoff()

deparse1WithCutoff(call, FALSE, DEFAULT_Cutoff, backtick, DEFAULTDEPARSE, /* nlines = */ 1);

[aryoda]

BTW:

1. I have limited the max number of output lines per call in the stack trace output to 10 lines.
2. I have also added an option tryCatchLog.max.lines.per.call to change the maximum number of source code rows printed in the full call stack per call.

[yoda-vid]

This is really awesome, continuing to learn a ton of R here. Everything looks good from my initial usage of the latest updates on this branch, including the initial case with the large sparse matrix.

[aryoda]

I am preparing the CRAN release new (planned end of October 2021)

[aryoda]

BTW: Core R had a similar performance problem in try() which was fixed in R4.1.1patched (Sept. 17, 2021):

https://stat.ethz.ch/pipermail/r-devel/2021-October/081159.html

[yoda-vid]

This is great, I can't thank you enough for it! Been using the fix daily without any issues. Glad Core R has benefited from a similar fix as well!

[aryoda]

FYI: Released at CRAN today as version 1.3.1:

https://cran.r-project.org/package=tryCatchLog

@yoda-vid Many thanks for your productive and valuable contribution and support!

[yoda-vid]

Installed the latest release, and it's working great.

Thank you again for this fix (and your detailed analysis) and putting out such an excellent and useful package!