Automatic Algorithm Selection in bigPLSR

When does each win?

• XtX (SIMPLS): moderate \(p\) (fits \(p^2\) in RAM). Fast BLAS-3; excellent for \(n \gg p\).
• XXᵗ (wide-kernel): moderate \(n\) (fits \(n^2\)). Great when \(p\gg n\) (wide problems).
• NIPALS / streaming: both \(p^2\) and \(n^2\) exceed budget; supports file-backed scores and large-scale chunked BLAS.

Sanity check

set.seed(1)
n <- 1e5; p <- 200
X <- matrix(rnorm(n*p), n, p)
y <- X[,1]*0.5 + rnorm(n)
bmX <- bigmemory::as.big.matrix(X)
bmy <- bigmemory::as.big.matrix(matrix(y, n, 1))

options(bigPLSR.memory_budget_bytes = 2L*1024^3)
fit <- pls_fit(bmX, bmy, ncomp=3, backend="bigmem", scores="r")
fit$algorithm

Overview

bigPLSR::pls_fit() can automatically choose an algorithm based on problem shape and a user-configurable memory budget:

• SIMPLS (XtX route) when forming the p × p cross-product fits in memory.
• SIMPLS (XXt / kernel route) when XtX does not fit but XXt (n × n) does.
• NIPALS (streaming) when neither XtX nor XXt comfortably fit.

This selection only applies when algorithm = "auto" (the default). Any explicit algorithm = overrides the decision.

The decision rule

Let the memory budget be B bytes (defaults to 8 GB, configurable via options(bigPLSR.mem_budget_gb = ...)). With doubles (8 bytes), we estimate the size of each symmetric matrix as:

• need_XtX = 8 * p^2
• need_XXt = 8 * n^2

Then:

  if (can_XtX && shape_XtX) { algo_in <- "simpls"}.    # XtX
  if (can_XXt && shape_XXt) { algo_in <- "widekernelpls"}. XXt (a.k.a. "kernel" route)
  else                      { algorithm <- "nipals"}         # streaming

Configuring the memory budget

# Use 16 GB as selection budget
options(bigPLSR.mem_budget_gb = 16)

This does not change R’s actual memory limit; it only controls the selection.

Reproducibility knobs

For tight numerical parity in tests:

set.seed(1)
if (requireNamespace("RhpcBLASctl", quietly = TRUE)) {
  RhpcBLASctl::blas_set_num_threads(1L)
  RhpcBLASctl::omp_set_num_threads(1L)
}
# otherwise, you can try environment variables:
# Sys.setenv(OPENBLAS_NUM_THREADS = "1", OMP_NUM_THREADS = "1")

Examples

library(bigPLSR)

n <- 2e3; p <- 5e2
X <- matrix(rnorm(n*p), n, p)
y <- X[,1] - 0.5*X[,2] + rnorm(n)

# Auto will likely pick SIMPLS (XtX) here
fit <- pls_fit(X, y, ncomp = 10, algorithm = "auto")
fit$algorithm  # "simpls"

Wide case:

n <- 200; p <- 4000
X <- matrix(rnorm(n*p), n, p)
y <- rnorm(n)

# If budget is small, auto picks kernel (XXt) or NIPALS
options(bigPLSR.mem_budget_gb = 2)  # small budget
fit <- pls_fit(X, y, ncomp = 5, algorithm = "auto")
fit$algorithm  # "kernelpls" or "nipals" depending on n^2 vs budget

Big-matrix streaming:

library(bigmemory)
n <- 1e6; p <- 50
# (example only; allocate according to your RAM)
# bmX <- big.matrix(n, p, type = "double")
# bmy <- big.matrix(n, 1, type = "double")
# fit <- pls_fit(bmX, bmy, ncomp = 10, backend = "bigmem", algorithm = "auto")
# fit$algorithm  # "simpls" or "nipals"

References

• de Jong, S. (1993). SIMPLS: An alternative approach to partial least squares regression. Chemometrics and Intelligent Laboratory Systems, 18(3), 251–263.
• Dayal, B., & MacGregor, J. F. (1997). Improved PLS algorithms. Journal of Chemometrics, 11(1), 73–85.
• Rosipal, R., & Trejo, L. J. (2001). Kernel Partial Least Squares Regression in Reproducing Kernel Hilbert Space. Journal of Machine Learning Research, 2, 97–123.
• Wold, H. (1966, 1985). NIPALS algorithm (original PLS formulation).

Appendix: streaming Gram math

For column blocks \(J\), \[ K \approx \sum_{J} X_{[:,J]} X_{[:,J]}^\top,\quad (Kv) \leftarrow (Kv) + X_{[:,J]} \big(X_{[:,J]}^\top v\big). \]

For row blocks \(B\), \[ K \approx \sum_{B} X_B X^\top,\quad (Kv) \leftarrow (Kv) + X_B \big(X^\top v\big)_B. \]

Center on the fly: \(H K H v = K v - \tfrac{1}{n}\mathbf{1}\mathbf{1}^\top K v - \tfrac{1}{n}K\mathbf{1}\mathbf{1}^\top v + \tfrac{1}{n^2}\mathbf{1}\mathbf{1}^\top K \mathbf{1}\,\mathbf{1}^\top v\). Maintain the needed aggregated vectors once per pass.