mlcollyer
commented
7 months ago Dear Ryan,
Let me give this some thought. I agree the function is not optimized and to be honest, it is now kind of old and could use some updating. I don’t think the use of the sapply function is the issue, and switching to the mclapply approach in your case might not help things too much. I think the problem is the enormous matrices involved in the calculations, both in terms of the memory needed to store them in the function environment plus make the calculations. When I originally contrived this function and its support functions, looking at options in ape and phytools as guidance, the goal was to have a matrix calculation that would be quicker than cycling through a recursive algorithm that finds the node-by-node ancestral character state estimates using PICs, variable-by-variable in a highly multivariate data set. This is what phytools::fastAnc does, using ape::ace, which relies on C to do the computational work. It is a nice approach with a single variable, but if done variable by variable, it suffers a similar problem that you ponder, why keep going to the same node multiple times to do similar calculations? (The way the algorithm works is to reroot the tree, node by node, calculate all PICs but retain just the first PIC as the ancestral state for the target node.) There is certainly a way to make an algorithm that alternatively reroots the tree at each node and then calculates the first PIC for each of many variables — that might be the best solution.
I think the anc.BM function we use is proficient to a point, but when the matrices become too large, the algebra needs to be abandoned for algorithms. Your example is elucidating this as an issue. The algorithms might also be slow but would not be as prohibitive. Ideally, writing C++ code that does what the ape::ace function does with PICs but for multivariate data would be the way to go, but that is quite an undertaking and one that would not be a simple package update. However, I think I can try to find a better solution than what we have now. It maybe will take a couple of weeks though.
Do you need to obtain the ancestral states for your analysis? Or is this just a roadblock for what you are trying to do? There might be faster ways to perform ordinations by avoiding gm.prcomp and using RRPP::ordinate, if you did not explicitly need to put a tree in the plot.
Thanks for the alert!
Mike
On Apr 15, 2024, at 5:13 PM, Ryan McMinds @.***> wrote:
Hello -
I'm experimenting with the phylogenetic PCA / PACA to analyze microbial communities. My datasets often have thousands of species. I have left gm.prcomp running for hours on what is actually a relatively small dataset for me, and have no idea how much longer it may take, so I started looking into ways to speed things up.
Within the gm.prcomp function, the first hurdle appears to be 'anc.BM'. I first thought the 'sapply' function within could be simply change to something like mclapply to at least parallelize things, but I tried it and it's still taking quite a while, and now takes a lot of memory. It looks to me like there is a lot of redundant computation here - for simple ancestral state prediction, if you start at the tips and move toward the root, then you should only have to visit each node once while accumulating the ancestral states, right? Whereas a function like sapply is going to process each node individually - and the nodes closer to the root have to re-calculate all the more-tipwise results.
I can probably sit and code up something myself, but I'm not sure if I'm missing something important! I'm a bit concerned that even if I put the time into speeding up the ANC portion, later parts are still going to be intractable (e.g. the naive covariance matrix inversion in phylo.mat... which might be possible to speed up with tricks like those used in MCMCglmm?).
Let me know if you have thought about these issues, have any ideas to optimize things for larger datasets, or if you'd be happy for me to try to contribute! Thank you!
— Reply to this email directly, view it on GitHub https://github.com/geomorphR/geomorph/issues/48, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABUNU4WEJFWXAUD7OEWUMITY5Q7IRAVCNFSM6AAAAABGIC5M4OVHI2DSMVQWIX3LMV43ASLTON2WKOZSGI2DINRTGA4TANY. You are receiving this because you are subscribed to this thread.
rmcminds
Hello -
I'm experimenting with the phylogenetic PCA / PACA to analyze microbial communities. My datasets often have thousands of species. I have left gm.prcomp running for hours on what is actually a relatively small dataset for me, and have no idea how much longer it may take, so I started looking into ways to speed things up.
Within the gm.prcomp function, the first hurdle appears to be 'anc.BM'. I first thought the 'sapply' function within could be simply change to something like mclapply to at least parallelize things, but I tried it and it's still taking quite a while, and now takes a lot of memory. It looks to me like there is a lot of redundant computation here - for simple ancestral state prediction, if you start at the tips and move toward the root, then you should only have to visit each node once while accumulating the ancestral states, right? Whereas a function like sapply is going to process each node individually - and the nodes closer to the root have to re-calculate all the more-tipwise results.
I can probably sit and code up something myself, but I'm not sure if I'm missing something important! I'm a bit concerned that even if I put the time into speeding up the ANC portion, later parts are still going to be intractable (e.g. the naive covariance matrix inversion in phylo.mat... which might be possible to speed up with tricks like those used in MCMCglmm?).
Let me know if you have thought about these issues, have any ideas to optimize things for larger datasets, or if you'd be happy for me to try to contribute! Thank you!