fast (and lean) RDA, CCA, dbRDA

[jarioksa]

New CRAN package fastrda promises faster RDA for large data sets. It uses C++ backend with OMP for parallel computation. The package is mainly intended for large genomic data with dimensions in thousands or tens of thousands. I had a quick look at the function and ran some comparisons against vegan – not because I doubt fastrda claims but to see if we can learn something from it and why would that be faster. I used ARM-based Mac M1 which may not be the best platform for fastrda due to poor or lacking OMP support in macOS. Mac ARM platforms also have optimized and very fast BLAS in a matrix algebra suprocessor ("neural engine"). While constrained ordination is written in R, it uses BLAS for all calculations and the R overhead is pretty small so that I do not expect much speed-up from switching from R to compiled code – most of the time is already spent in compiled code of BLAS/Lapack. All below is based on my particular platform, and results can be different with different hardware and OS. Parallel analyses are welcom.

vegan::rda was much faster in the first quick test with mite data. This was expected because fastrda::fastrda is designed for large data sets. The proper test used bryceveg & brycesite data of the labdsv package These are more realistic larger vegetation data with dimensions 160 × 169 though not yet large and pretty small when compared to genomic data. I used model matrix ~ elev + annrad + pos (rank 6). Here microbenchmark::microbenchmark got fastrda as a clear winner when not using M1 optimized BLAS: vegan::rda 23.3 msec, fastrda::fastrda 6.6 msec median times. With optimized BLAS both rda and fastrda used 6.0 msec.

The C++ backend is not the only thing that fastrda has for speed-up. It calculates (default) only three first constrained eigenvalues, and it does not calculate residual unconstrained solutions. So I had a look where vegan::rda spends its time with Rprof for a loop of 100 rda with the same data. Here a selection of relevant entries:

function	total used
rda	100.0%
ordConstrained	98.1%
svd	70.0%
ordResid	63.8%
ordConstrain	33.8%
qr.resid	22.4%

These numbers are with the slower "reference BLAS". Residual unconstrained solution uses nearly two times as much time as the constrained solutions. Moreover, qr.resid is only used in ordConstrain to prepare matrix passed on to ordResid so that the time cost of unconstrained component is 63.8 + 22.4 = 86.2% of running time.

Do we always need the unconstrained component? fastrda never calculates that component and that is one explanation why it is faster. We do know the sum of residual eigenvalues (residual inertia) without explicit calculation of eigenvalues and eigenvectors (and do this in vegan::simpleRDA2 that is used in vegan::varpart). Moreover, dropping residual component would make result object smaller: the constrained scores have dimensions 160 × 6, but in this case the unconstrained scores have dimensions 160 × 148. Especially with large genomic data sets with a huge number of columns, the cost of calculating unconstrained scores is huge and object sizes get really big. If we want to suppress the calculation of unconstrained residual component, this should be optional and the default should be to calculate.

You perhaps noticed that unconstrained scores had only 148 columns although the input data had dims 160 × 169. This happened because the residual matrix had no full rank but several last eigenvalues were practically zero and dropped from the vegan::rda result. Therefore vegan::rda reports the residual rank as 148, but fastrda as 153 (160 – 6 constraints - 1). This you only know if you evaluate the rank. vegan:::ordResid does this by looking at the eigenvalues, but this could also be estimated faster from the rank of the QR decomposition of the residual matrix (in ordResid(Y) with qr(Y)$rank). However, calculation of residual axis scores takes much time, and dropping their calculation can make a big speed-up: according to microbenchmark running svd(Y) took 10.4 msec, but svd(Y, nv=0, nu=0) 3.7 msec or nearly three times speed-up. So if we want to keep the eigenvalues we could optionally suppress only calculating scores. This alone could mean great savings both in time and object sizes.

See also Discussion #776 and Issue #527.

[jarioksa]

Small update to my initial message: I profiled rda with the optimized M1 BLAS of my Mac. With optimized BLAS the function was ca. 7 times faster (from 4.2 sec profiling time per 100 runs to 6.3 sec per 1000 runs). The same computer, the same R code, but using different BLAS. These comparisons are not general but they depend on your platform and I would really appreciate parallel analyses. The profiling table also changed:

function	total used
rda	100.0%
ordConstrained	92.1%
svd	66.8%
ordConstrain	45.6%
ordResid	39.2%
qr.resid	5.4%

Now the time cost of residual component is much smaller, with qr.resid 44.6%. What remains is that svd takes much time. According to microbenchmark calculating score vectors takes about half of the time so that calculating only singular values / eigenvalues would give a considerable speed-up also with optimized BLAS. Optional suppression of calculating unconstrained scores (but calculating eigenvalues) could make sense.