Potential performance improvement in predict() via aggregate() optimization

[RestlessTail]

Describe the issue
When running predictions on multi-block datasets (especially with cross-validation), I noticed significant slowdowns (~10s) in predict.mixo_pls. Profiling revealed that WeightedPredict calculation was the bottleneck (as mentioned in #236). I revisited the weighted vote part of predict.mixo_pls and ran additional correctness and benchmark tests.

Currently, the code uses:

out$WeightedPredict = array(unlist(lapply(temp.all, function(x){
	apply(x, c(1,2), function(z){
		temp = aggregate(rowMeans(object$weights),list(z),sum)
		ind = which(temp[,2]== max (temp[,2]))# if two max, then NA
		if(length(ind) == 1)
		{
			res = temp[ind, 1]
		} else {
			res = NA
		}
		res
	})})), dim(Y.hat[[1]]), dimnames = list(rownames(newdata[[1]]), colnames(Y), paste0("dim", c(1:min(ncomp[-object$indY])))))

As discussed earlier, this construct:

• repeatedly computes rowMeans(object$weights) even though it is constant,
• calls aggregate() on floating point values (almost no aggregation occurs),
• results in very slow predictions on larger datasets (especially regression with CV/permutation).

Expected behavior
This block can be replaced by:

out$WeightedPredict = array(unlist(lapply(temp.all, function(x){
	x[, , which.max(rowMeans(object$weights))]
})), dim(Y.hat[[1]]), dimnames = list(rownames(newdata[[1]]), colnames(Y), paste0("dim", c(1:min(ncomp[-object$indY])))))
out$temp.all.double <- is.double(unlist(temp.all)) # For testing purposes, always TRUE

Screenshots

I tested 4 affected models (Multi-Block PLS, Multi-Block sPLS, Multi-Block PLS-DA, Multi-Block sPLS-DA), each using random matrices generated by rnorm and diagonal matrices to simulate sparse omics data.

In each model, the value of temp.all is a floating point number, and the resulting WeightedPredict and WeightedPredict.class (for DA) are no different from the original function.

I also benchmarked each model, and for regression tasks, the performance improvement was more significant.

	Modified (ms)	Unmodified (ms)
Multi-Block PLS	67.11455	13010.12190
Multi-Block sPLS	65.91465	6566.60865
Multi-Block PLS-DA	97.25735	247.86015
Multi-Block sPLS-DA	97.06685	248.61410

In addition, the modified function also passed the mixOmics test suite.

Test script

Output of sessionInfo() (if it is an issue):

# session info
Platform: x86_64-w64-mingw32/x64
Running under: Windows 11 x64 (build 22631)

Matrix products: default


locale:
[1] LC_COLLATE=Chinese (Simplified)_China.utf8 
[2] LC_CTYPE=Chinese (Simplified)_China.utf8   
[3] LC_MONETARY=Chinese (Simplified)_China.utf8
[4] LC_NUMERIC=C                               
[5] LC_TIME=Chinese (Simplified)_China.utf8    

time zone: Etc/GMT-8
tzcode source: internal

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

other attached packages:
[1] microbenchmark_1.4.10 mixOmics_6.32.0       ggplot2_3.5.2        
[4] lattice_0.22-6        MASS_7.3-60.2        

loaded via a namespace (and not attached):
 [1] Matrix_1.7-0        gtable_0.3.6        dplyr_1.1.4        
 [4] compiler_4.4.1      tidyselect_1.2.1    Rcpp_1.0.12        
 [7] ellipse_0.5.0       stringr_1.5.1       parallel_4.4.1     
[10] gridExtra_2.3       tidyr_1.3.1         rARPACK_0.11-0     
[13] scales_1.4.0        BiocParallel_1.38.0 R6_2.6.1           
[16] plyr_1.8.9          generics_0.1.3      igraph_2.0.3       
[19] ggrepel_0.9.5       tibble_3.3.0        pillar_1.11.0      
[22] RColorBrewer_1.1-3  rlang_1.1.4         stringi_1.8.4      
[25] cli_3.6.3           withr_3.0.2         magrittr_2.0.3     
[28] grid_4.4.1          lifecycle_1.0.4     vctrs_0.6.5        
[31] RSpectra_0.16-2     glue_1.7.0          farver_2.1.2       
[34] corpcor_1.6.10      codetools_0.2-20    purrr_1.0.2        
[37] reshape2_1.4.4      matrixStats_1.3.0   tools_4.4.1        
[40] pkgconfig_2.0.3    

# package version
[1] "C:/Users/18269/AppData/Local/R/win-library/4.4/mixOmics"
 package  * version date (UTC) lib source
 mixOmics * 6.32.0  2025-09-23 [1] Github (mixOmicsTeam/mixOmics@ad47493)

[brendanlu]

Hey, thanks for looking into this - and apologies for the delay. mixOmics does not have a dedicated maintainer, and might not have one for some while.

However, given that your proposed changes represent a substantial performance improvement (with benchmarking evidence), and are passing the test suite locally, would you like to open a PR implementing your proposed changes + adding any non-regression unit tests? IIRC, the repository should then run testing across all different OS.

I see there are two optimizations you have made, being on the order of approx 6x and 1000x from your testing with the same numerical output. Perhaps open a PR with the more efficient method? I did not work on this part of the library, so I can't be certain the aggregate() method was a concious choice, but I'm leaning towards no; regardless, a 1000x speedup is enticing and in this case I can't immediately imagine a design / compatibility decision compelling enough for us to avoid such an optimization.

In late November, I should have some time to review your work and merge your changes as there's some housekeeping work I have to do on the repo anyway; I will also have the time to look a bit closer to make sure there's no convincing reason the old logic must be preserved. Obviously, only if you have the time, I would prefer you opened a PR seeing as you have basically done all the investigative work, and it can count as your contribution. If you don't have time, let me know and I can also easily implement them.

[RestlessTail]

Thanks for your detailed reply and clarification!
I've prepared the optimized implementation and a non-regression unit test as discussed. To be honest, I'm not very familiar with unit test conventions, so I followed the style of the existing tests in the repository, but I’m not entirely sure if it’s the best approach. I’d greatly appreciate any feedback or suggestions for improving the test. I’ll open a PR shortly to make it easier to review and test on your end.
Thank you for your time and openness to contributions. It’s great to see the care put into maintaining this project.

[frohart]

Hey guys,
Jumping in here quickly: there seems to be a difference between the 2 codes and I think it should be checked a bit more before being approved.
The original code has a special case when there's multiple max and that's probably why the 'which.max' function wasn't used as it returns only 1 position. There might be a smarter way to do that than with the aggregate function though..

Agreed that proposed solution is much faster, but I'd suggest to look at the behaviour of the 2 pieces of code with an example where there's multiple max. Otherwise we might have something faster but a lot more errors..

[RestlessTail]

Hi frohart
Thank you for the thoughtful review and for pointing out the edge case with multiple maxima. You're absolutely right about the difference in behavior. I'd like to discuss the necessity of the current implementation given the performance impact. The code contains two floating-point equality checks that seem to guard against extremely rare scenarios:

1. The first occurs in aggregate() where Y.hat values (continuous probabilities) are used as grouping factors - exact duplicates are extremely rare
2. The second in which(temp[,2] == max(temp[,2])) is even less likely, as we're comparing sums of weights that accumulate rounding errors, making exact equality virtually impossible

The current implementation spends more time to handle these edge cases. I'm wondering if this level of protection is necessary given their practical improbability, especially when it significantly impacts user experience during cross-validation and permutations.

That said, I completely understand the value of robust code. If we prefer to maintain the exact same behavior, here's a theoretically faster fallback:

out$WeightedPredict = array(unlist(lapply(temp.all, function(x){
if (all(apply(x, c(1,2), anyDuplicated.default) == 0L) && anyDuplicated(rowMeans(object$weights)) == 0L)
x[, , which.max(rowMeans(object$weights))]
else
# fallback to current implementation
})), dim(Y.hat[[1]]), dimnames = list(rownames(newdata[[1]]), colnames(Y), paste0("dim", c(1:min(ncomp[-object$indY])))))

I haven't benchmarked this new fallback approach yet, but since anyDuplicated.default is a built-in function with simpler logic, it should theoretically be faster than the complex aggregate workflow.