Thoughts about the integrating optimizations into the analysis stack

[LTLA]

Thinking about integration with BiocSingular in particular.

I think the best option would be to write a separate package with a S4 generic that takes an instance of a matrix-like-object and returns a BiocSingularParam object that provides the "best" choice of algorithm (in terms of the speed/accuracy trade-off). This approach has several advantages:

• Isolate implementation (in BiocSingular) from decisions about the choice of algorithm, which makes it easier to maintain as I only have to implement things rather than explicitly choose between them.
• Provide users with greater control - if they don't like the generated BiocSingularParam object, or if they know better (e.g., about their file system access speeds), they can just use their own.
• Simplify extensions for community-defined matrix representations, as anyone can write methods for the generic if they know that their representation is fast/slow at being touched.

The downside of the above strategy is that the choice of algorithm is not entirely transparent in the analysis stack, as the user has to actively call the new function that generates the new BiocSingularParam object. But I would argue that the choice of algorithm would not have been transparent in the first place, e.g., switching between IRLBA and RSVD can give results that differ beyond numerical precision, while switching between approximate and exact SVD will result in changes to the random seed.

Plus, if you write a separate package, you can also put in functions to perform empirical benchmarking for people who are really interested in optimizing their SVD speeds. Those won't go into BiocSingular.

[kasperdanielhansen]

I have not spent much time thinking about the practical implementation of a algorithmic selector. To some extent you have started doing this yourself with the auto selection of doing a cross product for wide/tall matrices. I agree that any implementation should have the following features

• Any choice should be user-overridable
• It should be transparent what choices are being made
  But for most users they also want convenience. One could have a separate package which analyses the problem and sets up the proposed approach. But I think BiocSingular should be able to act on it, so that for most users they don't need to write tons of code.

But for now my mind is on benchmarking. I ran an analysis script from Davide (330k cells) on a network mounted drive (ie. much slower than the SSDs we take advantages of in our Macs). Result: 4-7x slowdown. Specifically the PCA he did completed in 2.5h vs. 16h. Same number of cores (probably not same speed), but much slower IO. It may very well be that IO is by far the biggest bottleneck here.

[LTLA]

It wouldn't necessarily be tons of code. The best idea would be that a separate package could define its own BiocSingularParam class, and a runSVD method for it:

x <- runSVD(some_data, k=10, BSPARAM=YourOptimalParam())

... and the method would be able to choose the appropriate algorithm given the input matrix:

setMethod("runSVD", "YourOptimalParam",
    function (x, k, nu = k, nv = k, center = NULL, scale = NULL, 
        BPPARAM = SerialParam(), ..., BSPARAM) {
    # Can use class(x), dim(x), k, to make this choice.
    if (is(x, "HDF5Matrix")) {
        best <- RandomParam()
    } else {
        best <- IrlbaParam() 
    }

    # Re-dispatch.
    runSVD(x, k=k, nu=nu, nv=nv, center=center, scale=scale, 
        BPPARAM=BPPARAM, BSPARAM=best)
})

This new package can depend on whatever it wants (e.g., HDF5Array, other packages implementing matrix classes like bigmemory) without increasing the loading footprint of the basic BiocSingular package. It is also as easy to use as it can be - users just swap the BSPARAM= parameter, and they're optimizing. If the users don't want to optimize, then they shouldn't use YourOptimalParam(). Nice and simple and obvious.

It's an elegant separation of concerns on the user and the developer level. The latter is especially important as it means that I don't have to worry about your choices and you don't have to worry about mine.

[LTLA]

Incidentally, I have set the default for the cross-product determination to fold=Inf. For typical sparse matrices, computing the cross-product prior to irlba is actually pretty damaging because it needs to go via a DelayedMatrix. This realizes it into a dense intermediate unless deferred=TRUE.

So, I've just turned off the cross-product calculation by default, and leave it to your package to determine whether or not a cross-product is worth computing for a given input matrix. Depending on how it turns out, I may simplify the interface by replacing fold=<real> with a crossprod=<bool> argument.