bnprks
accurate and pass-effecient SVD method
Hi,
Nice work.
This is the author of PCAone and the R package. I am wondering if you're considering the accurate and much more pass-efficient SVD method in addition to the IRAM. The difference would be 10 vs 200 passes. I'm happy to contribute.
Best, Zilong.
Hi @Zilong-Li, thanks for reaching out! Your paper and R package seem like an interesting SVD method. I've previously had success implementing randomized SVD from Halko Algs 4.4 and 5.1 in pure R, building on top of the existing BPCells R interface functions. Your algorithm seems like a promising variant.
I think the most natural starting direction might be if you add BPCells compatibility in your PCAoneR package, i.e. implementing a pcaone.IterableMatrix method, and maybe adding BPCells to the package "suggests" list without introducing any hard dependencies. Then, once you've added that support we can link to your package from the BPCells documentation site so interested users will have an easier time finding your work.
I'm also happy to provide some advice along the way about what the performance considerations are for efficient matrix subsets with BPCells objects or getting multi-threaded sparse-dense matrix multiplies working. We have some plans for making C++-level BPCells extensions stable, but I think in this situation routing through the R interface will be a better initial approach as it is somewhat unlikely to cause performance bottlenecks for this kind of algorithm.
Reducing the number of matrix passes definitely has a lot of potential for improving disk-backed PCA performance, and it looks like CSV file formats probably hindered your performance in the paper benchmarks. I'd be very interested to see how performance and accuracy turn out on some of our own favorite benchmarking datasets when you can use a BPCells matrix back end hooked up to your SVD solver.
Happy to continue the thread here, or if you want to set up a call or have some more detailed conversations I'm also happy to communicate via email (bparks [at] stanford.edu)
-Ben
Hi @bnprks , thanks for being interested! I added winSVD function in pcaone package v1.5.0, which is the pure R implementation of the sliding window based randomized SVD. I followed the basic tutorial to generate the matrix mat_norm and tested the validity of winSVD for it. To clarify first, the winSVD is motivated by that learning the PCs of a large matrix through sliding window or mini-batchs in the first few passes can accelerate the process of minimizing the theoretical error. Similar idea recently proposed by Feng et. al 2024 employs a dynamic shifts of the singular values in each power iteration to accelerate it, which is also implemented in pcaone package called dashSVD. With the small matrix (1000 x 2600) from the tutorial, you may find dashSVD performs quite well and winSVD is less performant, which is expected since winSVD is designed for very large matrices or at least one dimension of it is way larger than the other. With large matrices, the default setting of winSVD , i.e. B=64, p=10, can have very accurate results in our experiments. Due to the similarity between winSVD and dashSVD, it`s possible to make a hybrid out of them. But I didn't see a clear scenario for that yet.
svd <- BPCells::svds(mat_norm, k=50)
s1 <- pcaone::winSVD(t(mat_norm), k=50, p = 20, s = 20, B = 16, perm = 0)
s2 <- pcaone::winSVD(t(mat_norm), k=50, p = 20, s = 30, B = 16, perm = 0)
s3 <- pcaone::winSVD(t(mat_norm), k=50, p = 20, s = 40, B = 16, perm = 0)
s4 <- pcaone::dashSVD(mat_norm, k=50, p = 20, s = 40)
par(mar = c(5, 5, 2, 1), pch = 16)
plot(svd$d - s1$d, xlab = "PC index", ylab = "Error of singular values", cex = 1.5, cex.lab = 2, cex.axis = 2)
points(svd$d - s2$d, col = "blue2", cex = 1.5)
points(svd$d - s3$d, col = "red2", cex = 1.5)
points(svd$d - s4$d, col = "green2", cex = 1.5)
legend("topleft", legend = paste0(c("winSVD, s=20, B=16", "winSVD, s=30, B=16", "winSVD, s=40, B=16", "dashSVD, s=40"), ", p=20" ), pch = 16,col = c("black", "blue2", "red2", "green2"), horiz = F, cex = 1.9, bty = "n" )
Currently, the pcaone function is the interface to Eigen matrix and operations with supports for dgCMatrix and dgRMatrix. To support BPCells at C++ level, definitely knowing the API for matrix slicing and multiplication is necessary. will email you about the details, but let me know your experiences with the pure R version of winSVD.
-Zilong
Hi @Zilong-Li, I haven't seen an email arrive from you -- if you sent one maybe there was an email address mixup?
I've been busy finalizing the BPCells preprint, but will have some time to check this out over the coming weeks. I'll try checking out accuracy and performance of the dashSVD implementation on some of the benchmark matrices from the BPCells manuscript I recently posted (link). That covers matrices up to 20M cells in the main benchmarking set, so it should be a pretty good test.
I think your R implementation of winSVD would benefit from a little BPCells-related customization to make sure that the matrix can be efficiently subset along the windowing axis. The storage_order() function in BPCells says this, and in some cases it is required to call transpose_storage_order() to create a matrix copy that supports efficient subsetting along the required axis.
It should be possible to test if I just manually make sure to pass the favorable storage ordering, but it's something to be aware of if you want to have more robust support for BPCells matrix objects.
-Ben