Thoughts on algorithmic improvements? (comparison to latest Rspectra)

[LTLA]

I was toying around with RSpectra, and damn, its truncated SVD got pretty good since you did your comparison a decade-ish ago.

y <- Matrix::rsparsematrix(10000, 5000, density=0.1)

system.time(iout <- irlba::irlba(y, nv=20, nu=20))
##    user  system elapsed 
##   7.862   0.004   7.866 

system.time(sout <- RSpectra::svds(y, k=20))
##    user  system elapsed 
##   4.656   0.000   4.657 

str(iout)
## List of 5
##  $ d    : num [1:20] 53.9 53.8 53.8 53.8 53.7 ...
##  $ u    : num [1:10000, 1:20] -0.02029 0.00435 0.01116 -0.01035 -0.0011 ...
##  $ v    : num [1:5000, 1:20] -0.02178 -0.0053 -0.0031 -0.00887 0.00627 ...
##  $ iter : int 172
##  $ mprod: int 1200

str(sout)
## List of 5
##  $ d    : num [1:20] 53.9 53.8 53.8 53.8 53.7 ...
##  $ u    : num [1:10000, 1:20] -0.02029 0.00435 0.01116 -0.01035 -0.0011 ...
##  $ v    : num [1:5000, 1:20] -0.02178 -0.0053 -0.0031 -0.00887 0.00627 ...
##  $ niter: num 26
##  $ nops : num 834

Almost twice as fast, and pretty much the same results once you take out the indeterminate sign. AFAICT this isn't a quality-of-implementation issue as my C++ port is about the same speed as your original R package and it uses different code/libraries for every step (matrix multiplication, the internal SVD, etc.). So there is some algorithmic improvement in Spectra's svds() that causes it to converge much faster, as evidenced by the lower niter and nops.

Is there any scope for similar improvements to IRLBA? It's not a big deal, but I'd like to try to squeeze a bit more performance out of my C++ library if I can. I don't know enough math to have any good ideas here (my best guess would be to put a cap on the number of iterations in the Lancoz process to avoid a quadratic increase in the number of matrix-vector multiplications with increasing nv or nu) but if you know of any relevant advancements, I'd be happy to implement them and try them out.

Session information

Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 22.04.5 LTS

Matrix products: default
BLAS:   /home/luna/Software/R/trunk/lib/libRblas.so 
LAPACK: /home/luna/Software/R/trunk/lib/libRlapack.so;  LAPACK version 3.12.1

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

time zone: Australia/Sydney
tzcode source: system (glibc)

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

loaded via a namespace (and not attached):
[1] compiler_4.6.0  Matrix_1.7-5    tools_4.6.0     Rcpp_1.1.1     
[5] RSpectra_0.16-2 grid_4.6.0      irlba_2.3.7     lattice_0.22-9 

[LTLA]

Some extra thoughts after some investigation:

As ?irlba suggests, bumping up work can speed things up, but I didn't realize exactly how much effect it would have. A few timings in libscran/irlba#21 indicates that much of the performance gap to svds() can be cut by just setting work to be proportional to the requested number of singular triplets. A natural choice is to just set the number of extra dimensions to be equal to the number of requested dimensions (work = 2 * nv in the irlba() context), and anecdotally this seems to be near the optimum work for each tested dataset.

So, I think there's a decent argument in favor of bumping up the default work=. Memory efficiency probably isn't a major concern; work is pretty small in the grand scheme of things, and if someone had enough memory to store work dimensions, they probably have enough memory to store 2 * work dimensions. In fact, poking around in svds()'s documentation indicates that their workspace defaults to 2 * number + 1 (omitting some of their constraints for simplicity), so if the effect is comparable, that would explain the differences in default behavior.

Another observation is that the implementation has a moderate effect on the performance. irlba.tests is a small R package that wraps my C++ port for easier testing (disclaimer: not intended for production use!), and it is faster than the original irlba package despite ostensibly performing the same calculations. I'd guess that this is because the matrix multiplication operations in the original package's C code are not as fast as those in Eigen, maybe because there's more overhead to call the relevant CHOLMOD routines from Matrix compared to Eigen's potentially-inlined multiplication.

set.seed(99)
y <- Matrix::rsparsematrix(10000, 5000, density=0.1)

# Setting tol= to be the same as the same default in svds.
# Not really sure that these arguments are comparable, but whatever.
library(irlba.tests)
system.time(iout <- irlba.tests::run_irlba(y, number=20, work=20, tol=1e-10))
##    user  system elapsed 
##   4.225   0.000   4.224 

library(irlba)
system.time(iout2 <- irlba::irlba(y, nv=20, nu=20, work=40, tol=1e-10))
##    user  system elapsed 
##   6.703   0.000   6.703 

library(RSpectra)
system.time(sout <- RSpectra::svds(y, k=20))
##    user  system elapsed 
##   4.735   0.000   4.734 

I also noticed that, generally speaking, IRLBA (regardless of the implementation) performs fewer matrix multiplications but requires more iterations compared to svds(). If we consider comparable Eigen-based implementations, this means that sometimes run_irlba() is faster and sometimes svds() is faster, depending on the relative speed of each iteration versus the multiplication. I'd speculate that IRLBA would win out at lower nu where the per-iteration SVD and multiplication of the working W/V matrices are cheaper, or for dense matrices where the multiplication is relatively expensive.

str(iout)
## List of 6
##  $ U              : num [1:10000, 1:20] -0.00335 -0.00861 -0.01101 -0.00307 0.00213 ...
##  $ V              : num [1:5000, 1:20] -0.00643 0.00151 -0.00704 -0.00951 0.00649 ...
##  $ D              : num [1:20] 54 53.9 53.8 53.8 53.7 ...
##  $ converged      : logi TRUE
##  $ iterations     : int 49
##  $ multiplications: int 830

str(sout)
## List of 5
##  $ d    : num [1:20] 54 53.9 53.8 53.8 53.7 ...
##  $ u    : num [1:10000, 1:20] -0.00335 -0.00861 -0.01101 -0.00307 0.00213 ...
##  $ v    : num [1:5000, 1:20] -0.00643 0.00151 -0.00704 -0.00951 0.00649 ...
##  $ niter: num 27
##  $ nops : num 848

# For comparison, let's ramp up the number of requested singular triplets.
y <- Matrix::rsparsematrix(10000, 5000, density=0.1)

system.time(iout.big <- irlba.tests::run_irlba(y, number=200, work=200, tol=1e-10))
##    user  system elapsed 
##  23.616   0.000  23.618 

system.time(sout.big <- RSpectra::svds(y, k=200))
##    user  system elapsed 
##  18.883   0.000  18.884 

str(iout.big)
## List of 6
##  $ U              : num [1:10000, 1:200] -0.0021 -0.0049 -0.00335 0.00499 -0.00295 ...
##  $ V              : num [1:5000, 1:200] 0.00108 -0.00791 0.00908 0.02931 -0.00244 ...
##  $ D              : num [1:200] 54.1 54.1 53.9 53.8 53.8 ...
##  $ converged      : logi TRUE
##  $ iterations     : int 19
##  $ multiplications: int 2432

str(sout.big)
## List of 5
##  $ d    : num [1:200] 54.1 54.1 53.9 53.8 53.8 ...
##  $ u    : num [1:10000, 1:200] 0.0021 0.0049 0.00335 -0.00499 0.00295 ...
##  $ v    : num [1:5000, 1:200] -0.00108 0.00791 -0.00908 -0.02931 0.00244 ...
##  $ niter: num 8
##  $ nops : num 2680

But in any case, the timings are so similar it doesn't seem to matter all that much.