fastPLS: Mathematical Methods, Implementations, and User Guide

Singular Value Decomposition Algorithms

fastPLS separates the PLS model from the singular-vector algorithm used inside the fit. fastPLS uses its native randomized SVD (rSVD) implementation following the randomized range-finder strategy of Halko, Martinsson, and Tropp (2011). In the compiled CPU implementation, a Gaussian test matrix is generated, the target matrix is sketched, optional power iterations are applied, the sketch is orthonormalized by an economic QR decomposition, and the final singular vectors are obtained from a much smaller projected matrix. When the randomized sketch dimension reaches the smaller matrix dimension, or when one matrix dimension is very small, the code falls back to an exact SVD because a truncated randomized approximation would no longer reduce the computation.

Native rSVD is the only public SVD solver. It is approximate: finite factors alone do not establish accurate singular directions or predictive agreement. Standalone fastsvd() reports singular-triplet checks and effective controls. If ordinary CPU sketches fail their checks, native rSVD retries with wider fresh sketches and additional power iterations. Failed or over-budget recovery raises an error. The recovery budget concerns major numerical buffers, not complete-process memory.

Sequential PLS fitting reports its seed, effective controls, execution profile, finite-factor checks and direction rule. These diagnostics do not establish equivalence to an independent PLS implementation. Compare several seeds and, when feasible, an independent reference before interpreting sensitive coefficients or latent subspaces. Increasing power or oversample can improve an approximation but is not a guarantee of agreement. Native rSVD is the only solver exposed by this package. PLS-SVD and standalone fastsvd() use oversample = 32, power = 5 by default. Accelerated SIMPLS, OPLS, and kernel-PLS use 32 oversampling directions and five power iterations for ordinary matrix shapes. Numeric responses with at least 64 columns and a response-to-sample ratio of at least 0.2 use 48 directions and six iterations. Classification with at least 32 classes and no more than 20 samples per class uses 64 directions and seven iterations. When the explicit predictor-response cross-covariance would exceed 512 MiB, the massive-matrix profile requests oversample = 12 and power = 1; executed controls and refresh width are recorded separately. CPU, Metal, and ordinary CUDA SIMPLS-family routes generate a seeded sketch from the current deflated cross-covariance. On a massive problem, CPU, Metal, and CUDA float64 use a rank-one randomized direction at each refresh. CUDA float32 can obtain up to eight independently initialized candidates per refresh through component 64 before returning to rank-one refreshes; CUDA keeps the sequential state on the device. For eligible large dummy-coded classification, CPU, CUDA, and Metal can instead obtain a candidate block of up to 64 directions. Candidates are accepted sequentially through the SIMPLS orthogonalization and deflation equations. Because the cross-covariance changes after every accepted component, the block route is an approximate sequential estimator and is not claimed to equal a fresh dominant-direction solve at every step. The explicit and matrix-free routes evaluate the same cross-covariance operator.

For large predictor-response cross-products, fastPLS can avoid forming S = X^T Y. This matrix-free route is called xprod in the code. The SVD backend evaluates products with S and S^T only when needed, for example Sv = X^T(Yv) and S^Tu = Y^T(Xu). This preserves the same mathematical operator while reducing memory pressure. For classification, labels can be used as PLS-DA responses without explicitly materializing a large dense one-hot matrix in the label-aware paths.

Partial Least Squares Algorithms

method = "plssvd" implements the direct PLS-SVD formulation, in which the dominant singular subspace of the predictor-response cross-covariance matrix is computed once and reused for all requested component counts. This makes PLS-SVD efficient when the response rank is the limiting dimension. The implementation adds compact latent prediction factors: instead of always storing and applying a full coefficient array, prediction can be performed through low-rank factors based on the latent projection and response-side coefficients. The cached score-Gram route forms the score cross-product once and reuses its leading blocks for the requested prefixes. If the centred cross-covariance is approximated by \(S \approx UDV^\top\), scores are \(T_a=XU_a\) and the response map solves \((T_a^\top T_a)H_a=D_aV_a^\top\). Predictions are \(X_{\mathrm{new}}U_aH_a\), with the training response mean restored. Singular directions alone are therefore not regression coefficients. PLS-SVD can also use the matrix-free xprod route for randomized SVD when the cross-covariance matrix would be large.

method = "simpls" is the optimized SIMPLS implementation and is the default general-purpose PLS core. The statistical structure follows de Jong’s SIMPLS algorithm (de Jong, 1993): each component is obtained from the current cross-covariance state, converted into a score and loading direction, and followed by SIMPLS orthogonalization and deflation. The rSVD route generates candidates from a newly randomized sketch at each refresh of the current deflated operator; a refresh can supply one or several candidates. Like pls::simpls.fit, one fit provides the sequential path through all components up to ncomp; fastPLS does not claim this standard path construction as a novelty. Cross-covariance deflation is updated after each accepted component. The component-wise CPU rSVD route derives a reproducible component-specific seed and generates a new Gaussian sketch at each refresh. Metal initializes one reproducible random sequence per fit and draws a new sketch for each component. CUDA generates newly initialized randomized candidates in resident device workspaces. Eligible large dummy-coded classification fits can generate up to 64 directions per CPU, CUDA, or Metal decomposition and consume them sequentially. The CPU route also evaluates candidate scores and predictor loadings in matrix blocks. Eligibility depends on predictor, response and sample dimensions, and the requested component count; it is not determined by sample count alone. This batched approximate route need not reproduce a fresh dominant-direction solve after each deflation. The ordinary CPU regression route remains component-wise; CUDA float32 additionally has a bounded massive-response block route. The C++ core also caches the rank-one deflation product, updates coefficient and fitted-value paths incrementally, and uses compact prediction factors. For low-dimensional high-sample problems, the code can cache cross-products such as X^T X and X^T Y; for large response or predictor dimensions, it can instead use matrix-free xprod products. These routes are selected by matrix shape: they reduce avoidable work or storage in their intended regimes but are not asserted to make every SIMPLS problem faster.

When the response dimension is no larger than the predictor dimension and is at most 512, the randomized direction calculation can use the smaller response-side Gram matrix \(S^\top S\). After each accepted component, fastPLS updates this matrix with the corresponding rank-one deflation product. This avoids repeatedly sketching the full tall cross-covariance while leaving the subsequent SIMPLS score, loading, orthogonalization, and deflation steps unchanged. The CPU implementation uses BLAS/LAPACK operations, CUDA keeps the Gram matrix and its updates in device memory for float32 and float64, and Metal uses its native float32 matrix operators.

method = "opls" implements orthogonal PLS following the orthogonal signal correction idea of Trygg and Wold (2002). It estimates response-orthogonal variation in X, subtracts this structured variation, stores the orthogonal filter, and then fits the optimized PLS core to the filtered matrix. The same filter is applied to new samples at prediction time.

method = "kernelpls" follows kernel PLS as described by Rosipal and Trejo (2001). Linear kernel PLS is dispatched to the ordinary linear PLS core to avoid unnecessary kernel materialization. For radial-basis and polynomial kernels, the training kernel is centered and the selected PLS core is fitted in kernel space. The same user-facing pls() function therefore covers linear PLS-SVD, SIMPLS, OPLS, and kernel PLS while sharing the accelerated low-rank machinery.

Classification Methods

Classification is handled as PLS-DA followed by one of two decoding heads. classifier = "argmax" uses the largest predicted dummy-response score and is the standard PLS-DA decoding rule. classifier = "lda" fits a regularized linear discriminant analysis model in the PLS score space following Fisher’s discriminant principle (Fisher, 1936).

Computational Backends

The backend argument selects where the main numerical work is performed. backend = "cpu" uses the package’s compiled, dependency-free C++ core with native rSVD. It is usually best for small and moderate datasets, low-dimensional problems, and machines without GPU support. The CPU kernels use Apple Accelerate by default on macOS. Linux and Windows prefer OpenBLAS when configuration finds it through OPENBLAS_ROOT, pkg-config, or Rtools; otherwise the package uses the BLAS/LAPACK supplied by R. Set FASTPLS_USE_OPENBLAS=1 before installation to require OpenBLAS. fastPLS_blas() reports the library selected when the installed package was compiled. Linux and Windows performance benchmarks should proceed only when this function returns "OpenBLAS". Linking against multithreaded OpenBLAS enables eligible matrix operations to use multiple CPU cores. For example, options(backend = "cpu", n.cores = 4L) selects CPU execution and requests four threads from the linked BLAS/OpenMP runtime. Explicit backend and n.cores arguments override their respective session options. The same resolved n.cores value is retained by nested folds, refits, permutations, and predictions. The sequential PLS deflation steps remain serial, so additional cores do not guarantee a speed improvement for every matrix shape. Additional BLAS threads help only when the matrix products are sufficiently large; sequential updates and thread-management overhead can limit or reverse the benefit. Compare timings on the intended workload rather than assuming that more threads are always faster.

For CPU self-products, fastPLS selects GEMM or the symmetric rank-k operation SYRK from the linked numerical library according to platform, precision, and matrix dimensions. SYRK computes one triangle; the other triangle is generated only when required by a later operation. Owned workspaces are aligned, and wide response matrices gathered for repeated fold computations use 64-byte column alignment. Compiled classification fitting and cross-validation products operate on compact labels and class sums, so they do not allocate a dense one-hot training-response matrix.

For eligible wide-response SIMPLS cross-validation, the compiled engine forms the uncentred sample-response Gram matrix once before the fold loop. Each training-fold Gram matrix is recovered from its principal submatrix and exactly double-centred, avoiding a materialized centred response matrix. CPU execution uses direct SSYRK or DSYRK; float32 CUDA uses cublasSsyrk; and the hybrid Metal route uses Accelerate SYRK for this host-visible operation. The raw Gram matrix retains only its lower triangle, and fold, output, and device workspaces are reused across the cross-validation call.

backend = "cuda" uses CUDA-native routines when the package is compiled with CUDA and a compatible NVIDIA runtime is available. After input transfer, CUDA paths perform preprocessing, cross-products, randomized direction extraction, component updates, low-rank prediction, and LDA scoring on the GPU for PLS-SVD, SIMPLS, OPLS, and kernel PLS. OPLS filtering and application of the stored orthogonal filter remain device-resident. Nonlinear kernel PLS constructs and centers both training and test Gram matrices on the GPU. Returned arrays and R model-object assembly are host operations. CUDA is most useful when dense matrix work offsets transfer and kernel-launch overhead. Nonlinear kernels are also subject to an explicit device-memory guard because their Gram matrix uses quadratic storage in the training sample count.

backend = "metal" provides a float32 route for PLS-SVD, SIMPLS, OPLS, and kernel PLS on Apple silicon unified-memory systems. It always enters the same operation-split algorithm rather than selecting a complete CPU or Metal estimator from the dataset dimensions. The CPU performs preprocessing, QR and reduced decompositions, sequential orthogonalization, deflation bookkeeping, model assembly, LDA, and prediction. Persistent Metal workspaces perform fitting products involving the training sample matrix. These include explicit cross-covariance formation, fused score/loading products, and randomized range-finder products when the cross-covariance is represented implicitly. In the transpose part of that implicit operator, consecutive products X V and Y^T (X V) share one Metal command buffer and one final synchronization before the CPU applies the response-centering correction. Reusable buffers and Metal matrix objects are retained for this sequence during the fit. The smaller cross-covariance and component-state operations remain on the CPU. This mathematical assignment is fixed and does not switch according to dataset dimensions. The route requires Metal support and float32 input and never silently converts a float64 request to CPU.

Bioconductor Data Containers

pls() and predict() accept a Biobase::ExpressionSet directly. Expression measurements are stored as variables by samples in an ExpressionSet, so fastPLS transposes the assay internally to its samples-by-variables layout. The following example predicts the case-control phenotype in the Biobase example data while keeping the test samples separate from model fitting.

data("sample.ExpressionSet", package = "Biobase")
y_bioc <- factor(Biobase::pData(sample.ExpressionSet)$type)

set.seed(42)
test_bioc <- sample(seq_len(ncol(sample.ExpressionSet)), 8)
fit_bioc <- pls(
    sample.ExpressionSet[, -test_bioc],
    y_bioc[-test_bioc],
    sample.ExpressionSet[, test_bioc],
    y_bioc[test_bioc],
    ncomp = 1:2,
    method = "simpls",
    backend = "cpu",
    seed = 42
)
fit_bioc$accuracy
#> ncomp=1 ncomp=2 
#>   0.875   0.750

Choosing a Starting Model

The public classification choices are "argmax" and "lda"; regression ignores classifier. LDA uses the fixed scale-normalized Cholesky fallback sequence described above and does not expose a ridge-tuning argument.

Fit And Predict A Classifier

The method argument selects the PLS algorithm, backend selects the implementation, and classifier selects the classification head. The default classification head is classifier = "argmax", which predicts the class with the largest PLS-DA response score.

fit_cls <- pls(
    Xtrain,
    Ytrain_cls,
    Xtest,
    Ytest_cls,
    ncomp = 1:2,
    fit = TRUE,
    return_variance = FALSE,
    seed = 101
)

fit_cls$accuracy
#>   ncomp=1   ncomp=2 
#> 0.6666667 0.7666667

Classification models can also be fitted once and predicted later. predict() returns only the predicted class by default and can optionally return ranked classes by setting top to a positive integer. The ordinary prediction in Ypred is always the rank-1 class. When top > 1, Ypred_top contains the ordered candidate labels for each sample: rank1 is the predicted class, rank2 is the next most likely class, and so on. If available, Ypred_top_score contains the corresponding class scores used to create that ranking. For both float64 and float32 models, the default raw_scores = FALSE evaluates ranked predictions in bounded row blocks and retains only the requested ranks. Set raw_scores = TRUE only when the complete class-score arrays are needed, because those arrays can be substantially larger.

fit_cls_train_only <- pls(
    Xtrain,
    Ytrain_cls,
    ncomp = 1:2,
    classifier = "lda",
    fit = TRUE,
    return_variance = FALSE,
    seed = 101
)

pred_cls_later <- predict(
    fit_cls_train_only,
    Xtest,
    Ytest = Ytest_cls,
    top = 2,
    raw_scores = TRUE
)

pred_cls_later$accuracy
#> ncomp=1 ncomp=2 
#>       1       1
pred_cls_later$metrics$metrics
#>          n accuracy no_information_rate lift_accuracy
#> ncomp=1 30        1           0.3333333             3
#> ncomp=2 30        1           0.3333333             3
#>         balanced_accuracy macro_precision macro_recall macro_f1
#> ncomp=1                 1               1            1        1
#> ncomp=2                 1               1            1        1
#>         kappa
#> ncomp=1     1
#> ncomp=2     1
head(pred_cls_later$Ypred_top[["ncomp=2"]])
#>      rank1        rank2       
#> [1,] "virginica"  "versicolor"
#> [2,] "virginica"  "versicolor"
#> [3,] "setosa"     "versicolor"
#> [4,] "versicolor" "virginica" 
#> [5,] "versicolor" "virginica" 
#> [6,] "versicolor" "virginica"

Evaluate Classification Predictions

Use evaluate() to summarize predicted class labels, ranked labels, class-score matrices, or a complete object returned by predict(). The task and requested rank are inferred from these inputs. For classification, lift_accuracy is accuracy divided by the no_information_rate, the accuracy obtained by always predicting the most frequent observed class. Values above one therefore improve on this simple majority-class baseline. For classification, the complete output includes global metrics, per-class metrics, and the confusion matrix.

eval_cls_path <- evaluate(
    observed = Ytest_cls,
    predicted = pred_cls_later
)

eval_cls <- eval_cls_path$by_component[["ncomp=2"]]
eval_cls
#> $task
#> [1] "classification"
#> 
#> $metrics
#>    n accuracy no_information_rate lift_accuracy balanced_accuracy
#> 1 30        1           0.3333333             3                 1
#>   macro_precision macro_recall macro_f1 kappa
#> 1               1            1        1     1
#> 
#> $metric_definitions
#> $metric_definitions$accuracy
#> [1] "Proportion of observed labels predicted correctly."
#> 
#> $metric_definitions$balanced_accuracy
#> [1] "Unweighted mean of class-specific recalls."
#> 
#> 
#> $per_class
#>        class support precision recall f1
#> 1     setosa      10         1      1  1
#> 2 versicolor      10         1      1  1
#> 3  virginica      10         1      1  1
#> 
#> $confusion
#>             observed
#> predicted    setosa versicolor virginica
#>   setosa         10          0         0
#>   versicolor      0         10         0
#>   virginica       0          0        10
#> 
#> $topk
#>   k accuracy
#> 1 1        1
#> 2 2        1

The same complete evaluation is also stored automatically in a fitted PLS object when observed test responses are supplied. Results are grouped first by data role and then by component count.

fit_cls$metrics$test[["ncomp=2"]]$metrics
#>    n  accuracy no_information_rate lift_accuracy balanced_accuracy
#> 1 30 0.7666667           0.3333333           2.3         0.7666667
#>   macro_precision macro_recall  macro_f1 kappa
#> 1       0.7765568    0.7666667 0.7612958  0.65

Rows of the confusion matrix are predicted labels and columns are observed labels. The confusion matrix is returned as an ordinary R table.

eval_cls$confusion
#>             observed
#> predicted    setosa versicolor virginica
#>   setosa         10          0         0
#>   versicolor      0         10         0
#>   virginica       0          0        10

When class-score matrices are available, evaluate() can also report top-k accuracy. Top-k accuracy asks whether the true class appears anywhere among the first k ranked labels. Thus top-1 accuracy is ordinary classification accuracy, while top-5 accuracy gives credit when the correct class appears among the five highest-scoring alternatives.

score_last <- pred_cls_later$LDA_scores[
    , , dim(pred_cls_later$LDA_scores)[3L]
]

evaluate(
    observed = Ytest_cls,
    predicted = score_last
)
#> $task
#> [1] "classification"
#> 
#> $metrics
#>    n accuracy no_information_rate lift_accuracy balanced_accuracy
#> 1 30        1           0.3333333             3                 1
#>   macro_precision macro_recall macro_f1 kappa
#> 1               1            1        1     1
#> 
#> $metric_definitions
#> $metric_definitions$accuracy
#> [1] "Proportion of observed labels predicted correctly."
#> 
#> $metric_definitions$balanced_accuracy
#> [1] "Unweighted mean of class-specific recalls."
#> 
#> 
#> $per_class
#>        class support precision recall f1
#> 1     setosa      10         1      1  1
#> 2 versicolor      10         1      1  1
#> 3  virginica      10         1      1  1
#> 
#> $confusion
#>             observed
#> predicted    setosa versicolor virginica
#>   setosa         10          0         0
#>   versicolor      0         10         0
#>   virginica       0          0        10
#> 
#> $topk
#>   k accuracy
#> 1 1        1
#> 2 2        1
#> 3 3        1

Classification Heads

The same pls() interface exposes two classification-specific heads for factor responses: argmax PLS-DA and latent-space LDA. They are decoders applied after the PLS model has produced class-response scores or latent scores. For regression, leave classifier at its default; numeric responses are predicted directly as continuous values.

fit_cls_plssvd <- pls(
    Xtrain,
    Ytrain_cls,
    Xtest,
    Ytest_cls,
    ncomp = 1:2,
    method = "plssvd",
    seed = 100
)

head(fit_cls_plssvd$Ypred)
#>     ncomp=1    ncomp=2
#> 1 virginica versicolor
#> 2 virginica  virginica
#> 3    setosa     setosa
#> 4 virginica  virginica
#> 5 virginica versicolor
#> 6 virginica  virginica

evaluate(
    observed = Ytest_cls,
    predicted = fit_cls_plssvd$Ypred[["ncomp=2"]]
)$confusion
#>             observed
#> predicted    setosa versicolor virginica
#>   setosa         10          0         0
#>   versicolor      0          5         2
#>   virginica       0          5         8

For PLS-DA with an LDA prediction head, use classifier = "lda". On systems with GPU support, the backend is selected through backend = "cuda" or backend = "metal" where available:

fit_cls_lda_gpu <- pls(
    Xtrain,
    Ytrain_cls,
    Xtest,
    Ytest_cls,
    ncomp = 1:2,
    method = "plssvd",
    backend = "cuda",
    classifier = "lda"
)

An unavailable CUDA or Metal selection raises an error; fastPLS never changes the request silently to CPU. To run on the CPU, select it explicitly:

fit_cls_lda_cpu <- pls(
    Xtrain,
    Ytrain_cls,
    Xtest,
    Ytest_cls,
    ncomp = 1:2,
    method = "plssvd",
    backend = "cpu",
    seed = 100,
    classifier = "lda"
)

head(fit_cls_lda_cpu$Ypred)
#>      ncomp=1    ncomp=2
#> 1  virginica  virginica
#> 2  virginica  virginica
#> 3     setosa     setosa
#> 4 versicolor versicolor
#> 5 versicolor versicolor
#> 6 versicolor versicolor

Kernel PLS For Classification

Kernel PLS changes the representation of the samples before the inner PLS fit. The linear kernel is equivalent to an ordinary inner-product representation and is useful as a fast baseline. The rbf kernel uses a radial-basis similarity; the poly kernel uses polynomial interactions among features.

kernel_fits <- lapply(c("linear", "rbf", "poly"), function(k) {
    pls(
    Xtrain, Ytrain_cls, Xtest, Ytest_cls,
    ncomp = 1:2,
    method = "kernelpls",
    kernel = k,
    degree = 2,
    seed = 102
    )
})
names(kernel_fits) <- c("linear", "rbf", "poly")

kernel_accuracy <- vapply(kernel_fits, function(fit) {
    mean(fit$Ypred[["ncomp=2"]] == Ytest_cls)
}, numeric(1))
kernel_accuracy
#>    linear       rbf      poly 
#> 0.7666667 0.9000000 0.9333333

Classification Score Plots

plot() can visualize stored PLS score maps. Refit with fit = TRUE to store training scores, or predict with proj = TRUE to store test scores. Ellipses can be ordinary confidence ellipses or Hotelling T2 ellipses.

plot(
    fit_cls,
    groups = Ytrain_cls,
    ellipse = TRUE,
    ellipse.type = "confidence"
)

The next example uses the two-class gene-expression data from the Biobase sample.ExpressionSet object introduced above. The same plot() call works for method = "plssvd", "simpls", "opls", and "kernelpls".

X_two <- t(Biobase::exprs(sample.ExpressionSet))
y_two <- factor(Biobase::pData(sample.ExpressionSet)$type)
train_two <- setdiff(seq_along(y_two), test_bioc)
fit_two <- pls(
    X_two[train_two, , drop = FALSE], y_two[train_two],
    X_two[test_bioc, , drop = FALSE], y_two[test_bioc],
    ncomp = 1:2, fit = TRUE, proj = TRUE, seed = 200
)

plot(fit_two, score.set = "train", groups = y_two[train_two],
        ellipse = TRUE, main = "Biobase expression data")

The second score-plot example fits OPLS to the three-class iris split defined at the beginning of this section, with two orthogonal components (north = 2).

opls_three <- pls(
    Xtrain,
    Ytrain_cls,
    Xtest,
    Ytest_cls,
    ncomp = 1:2,
    method = "opls",
    fit = TRUE,
    proj = TRUE,
    north = 2L,
    seed = 201
)

old_par <- par(mfrow = c(1, 2), mar = c(4.2, 4.2, 2.2, 0.8))
plot(
    opls_three,
    score.set = "train",
    groups = Ytrain_cls,
    ellipse = TRUE,
    ellipse.type = "hotelling",
    main = "OPLS training, north = 2"
)
plot(
    opls_three,
    score.set = "test",
    groups = opls_three$Ypred[["ncomp=2"]],
    xlim = c(-0.16, 0.42),
    main = "OPLS test prediction"
)

par(old_par)

Fit And Predict A Regression Model

The simplest workflow is to provide both the training data and an independent test set directly to pls(). In this case the fitted object already contains test-set predictions and, if Ytest is supplied, predictive metrics.

fit_reg <- pls(
    Xreg_train,
    Ytrain_reg,
    Xreg_test,
    Ytest_reg,
    ncomp = 1:3,
    fit = TRUE,
    return_variance = FALSE
)

fit_reg$Q2Y
#>   ncomp=1   ncomp=2   ncomp=3 
#> 0.7564006 0.7593961 0.8634626

Float32 Input

pls() also accepts predictor and response matrices from the float package. Supplying a float::float32 object requests single-precision execution. CUDA keeps PLS-SVD, SIMPLS, OPLS, nonlinear kernel construction, prediction, and LDA on the selected GPU. Metal accepts float32 only and uses the fixed CPU/Metal operation split described above; CUDA supports float32 and float64. Nonlinear kernel PLS requires an explicit n by n Gram matrix, so the package checks the estimated live device or unified-memory storage before fitting. Extreme-response and nonlinear-kernel routes should be checked against float64 when the scientific conclusion is sensitive to small numerical differences.

Float64 remains the baseline precision. Float32 approximately halves raw matrix storage, but paired tests show that runtime, incremental host RSS, and device memory may increase or decrease depending on the route. The package therefore warns once for measured-risk regimes. This includes precision-sensitive SIMPLS and linear kernel-PLS classification, nonlinear kernels and multivariate regression with at least 10,000 response columns and 50 components. In that last regime, runtime, memory, and numerical behavior must be checked against float64 before the float32 result is used for scientific interpretation. Unsupported combinations stop before allocation and are never silently promoted to float64.

On Windows, the standard R toolchain does not expose the single-precision BLAS/LAPACK symbols used by the Unix-like compiled kernels. The Windows CPU route therefore combines float-package rSVD and Cholesky operations with portable C++ float kernels. It supports PLS-SVD, SIMPLS, OPLS, linear and nonlinear kernel PLS, argmax, and latent-space LDA when backend = "cpu" and rSVD is selected automatically. Model factors, LDA buffers, kernel matrices, and predictions remain float32, but this portable route can be slower than the native Unix-like implementation. Windows float32 accelerator requests stop with an unsupported-combination error rather than silently converting data to float64.

Xreg32 <- float::fl(as.matrix(Xreg_train))
Yreg32 <- float::fl(matrix(Ytrain_reg, ncol = 1))
fit_reg32 <- pls(
    Xreg32,
    Yreg32,
    float::fl(as.matrix(Xreg_test)),
    float::fl(matrix(Ytest_reg, ncol = 1)),
    ncomp = 1:2
)
fit_reg32$Q2Y
#>   ncomp=1   ncomp=2 
#> 0.7564006 0.7593962

For standard double-precision regression, Ypred is a numeric prediction array with one slice for each requested number of components. For float32 input, Ypred is kept as a named list of float::float32 prediction matrices. In the standard example below, the last array slice corresponds to the largest requested component count (ncomp = 3). The metric vectors are named by component count, for example fit_reg$Q2Y["ncomp=3"].

reg_component <- dim(fit_reg$Ypred)[3L]
pred_mpg <- fit_reg$Ypred[, , reg_component]
plot(
    Ytest_reg,
    pred_mpg,
    pch = 21,
    bg = "#4E79A7",
    col = "black",
    xlab = "Observed mpg",
    ylab = "Predicted mpg",
    main = "Regression: observed vs predicted"
)
abline(0, 1, col = "#D55E00", lwd = 2)

The alternative workflow is to fit the model once without a test set, then call predict() later. This is useful when the same model must be applied to several independent datasets. predict() automatically applies the centering/scaling stored in the fitted object. When Ytest is supplied, it also passes the predictions to evaluate() and stores the complete evaluation under metrics.

fit_reg_train_only <- pls(
    Xreg_train,
    Ytrain_reg,
    ncomp = 1:3,
    fit = TRUE,
    return_variance = FALSE
)

pred_reg_later <- predict(
    fit_reg_train_only,
    Xreg_test,
    Ytest = Ytest_reg,
    proj = TRUE
)

pred_reg_later$Q2Y
#>   ncomp=1   ncomp=2   ncomp=3 
#> 0.7564006 0.7593961 0.8634626
pred_reg_later$metrics$metrics
#>             n        R2        Q2     RMSD     RMSE      MAE
#> component=1 8 0.7505150 0.7564006 3.066116 3.066116 2.739049
#> component=2 8 0.7535829 0.7593961 3.047205 3.047205 2.765666
#> component=3 8 0.8601637 0.8634626 2.295494 2.295494 1.929671
#>                  bias MRE_percent MAPE_percent      RPD Pearson_r
#> component=1 0.4938232    16.70827     14.92070 2.140295 0.8700688
#> component=2 0.6250298    17.87218     15.26995 2.153578 0.8743339
#> component=3 0.6784321    12.48956     10.76821 2.858815 0.9340810
#>             Spearman_r
#> component=1  0.8072875
#> component=2  0.8072875
#> component=3  0.8072875
head(pred_reg_later$Ttest)
#>             [,1]        [,2]       [,3]
#> [1,]  0.09841041 -0.01046739 -0.1900955
#> [2,] -0.10644825  0.19820127  0.1563451
#> [3,]  0.03394295  0.21190609 -0.2208999
#> [4,] -0.36597585  0.27972780 -0.2105406
#> [5,]  0.27550586  0.09182286  0.3494802
#> [6,] -0.25506247  0.31745909  0.2627156

Evaluate Regression Predictions

For numeric regression, evaluate() reports R2, Q2, RMSD/RMSE, MAE, bias, median relative error percentage, RPD, and correlations. If ytrain is supplied, Q2 is calculated relative to the training-set response mean, which is the preferred setting for independent test-set evaluation. For a response matrix, each response column is centered on its corresponding training mean before the denominator sums are aggregated. If ytrain is omitted, Q2 is returned as NA; it is not silently replaced by an R2 calculation.

For wide multivariate responses, pls(), pls.single.cv(), and pls.double.cv() return aggregate evaluation metrics by default. Set bycol = TRUE to also calculate the full response-wise metric table; this is useful for inspecting individual spectral bins, but can be expensive for NMR data with many response columns.

eval_reg <- evaluate(
    observed = Ytest_reg,
    predicted = pred_mpg,
    ytrain = Ytrain_reg
)

eval_reg$task
#> [1] "regression"
eval_reg$metrics
#>   n        R2        Q2     RMSD     RMSE      MAE      bias
#> 1 8 0.8601637 0.8634626 2.295494 2.295494 1.929671 0.6784321
#>   MRE_percent MAPE_percent      RPD Pearson_r Spearman_r
#> 1    12.48956     10.76821 2.858815  0.934081  0.8072875
lapply(eval_reg$metric_definitions, strwrap, width = 56)
#> $R2
#> [1] "Observed-set R2; each response is centered on its"   
#> [2] "observed mean before sums of squares are aggregated."
#> 
#> $Q2
#> [1] "Independent-test Q2; each response is centered on its"
#> [2] "training-response mean before sums of squares are"    
#> [3] "aggregated."
eval_reg$per_response
#>   response n        R2        Q2     RMSD     RMSE      MAE
#> 1       Y1 8 0.8601637 0.8634626 2.295494 2.295494 1.929671
#>        bias MRE_percent MAPE_percent      RPD Pearson_r Spearman_r
#> 1 0.6784321    12.48956     10.76821 2.858815  0.934081  0.8072875

OPLS For Regression

OPLS is accessed through the same pls() function with method = "opls". It is often used to separate predictive variation from response-orthogonal variation. ncomp counts predictive components, separately from the orthogonal components removed by north. Each removed direction reduces the available predictor rank. For example, with eight independent predictors and one removed direction, at most seven predictive components remain. A request exceeding the remaining rank raises an error; the same restriction applies within each CV training fold.

fit_opls <- pls(
    Xreg_train,
    Ytrain_reg,
    Xreg_test,
    Ytest_reg,
    ncomp = 1:2,
    method = "opls",
    seed = 101
)

class(fit_opls)
#> [1] "fastPLSOpls" "fastPLS"
fit_opls$Q2Y
#>   ncomp=1   ncomp=2 
#> 0.7593961 0.8634626

Ordinary k-fold CV

In ordinary k-fold CV, samples are split directly into kfold folds. The example below validates a fixed two-component SIMPLS classifier with five folds.

cv_kfold <- pls.single.cv(
    Xdata = Xtrain,
    Ydata = Ytrain_cls,
    ncomp = 2,
    kfold = 5,
    seed = 103
)

cv_kfold$metrics$cross_validated[["ncomp=2"]]$metrics
#>     n  accuracy no_information_rate lift_accuracy balanced_accuracy
#> 1 120 0.8333333           0.3333333           2.5         0.8333333
#>   macro_precision macro_recall  macro_f1 kappa
#> 1       0.8422619    0.8333333 0.8320269  0.75

metrics$cross_validated contains the complete evaluate() output for each tested component count, while selection_metrics is the compact internal table used to select the best setting. When fit = TRUE, metrics$fitted evaluates the corresponding model fitted to the full dataset.

Grouped k-fold CV with constrain

The constrain argument controls grouped splitting. It is a vector with one entry per sample; samples with the same value are assigned to the same fold. In practice, this prevents leakage when multiple rows come from the same patient, subject, batch, or technical replicate. For example, if two spectra come from the same patient, giving them the same patient identifier in constrain ensures that both spectra are placed either in the training set or in the test set, never one in each.

patient_id <- rep(
    seq_len(ceiling(nrow(Xtrain) / 2)), each = 2
)[seq_len(nrow(Xtrain))]

cv_grouped <- pls.single.cv(
    Xdata = Xtrain,
    Ydata = Ytrain_cls,
    constrain = patient_id,
    ncomp = 2,
    kfold = 4,
    seed = 104
)

c(n_folds = length(unique(cv_grouped$fold)),
    n_patient_groups = length(unique(patient_id)))
#>          n_folds n_patient_groups 
#>                4               60

Leave-one-out CV and leave-one-group-out CV

Leave-one-out CV is requested with kfold = "loocv". When constrain is not supplied, each sample is held out once. When constrain is supplied, LOOCV becomes leave-one-constraint-group-out CV, so a whole patient, subject, batch, or replicate group is held out together at each iteration. Numeric kfold values greater than or equal to the number of constraint groups are also treated as leave-one-group-out CV.

cv_loocv <- pls.single.cv(
    Xdata = Xtrain,
    Ydata = Ytrain_cls,
    constrain = patient_id,
    ncomp = 1,
    kfold = "loocv",
    seed = 105
)

c(n_loocv_folds = length(unique(cv_loocv$fold)),
    n_patient_groups = length(unique(patient_id)))
#>    n_loocv_folds n_patient_groups 
#>               60               60

Component and hyperparameter optimization CV

pls.single.cv() repeats the same CV splitting strategy over several candidate component counts and returns the best value according to the predictive metric: accuracy for classification and RMSD for regression. Fold-aware Q2Y, full-data fitted R2Y, balanced accuracy, and the other task-specific metrics documented below can instead be requested explicitly through selection. Predictive arguments can also be supplied as vectors to tune a compact grid. Fold-control arguments such as kfold remain single settings for the whole CV run.

For LDA classification, a rare class may be absent from an inner training fold. fastPLS fits LDA to the classes represented in that fold and maps the predictions back to the original factor levels. Held-out observations from an absent class are retained in the accuracy and balanced-accuracy calculations, although that class cannot be predicted by a model that did not observe it. When otherwise identical argmax and LDA configurations are requested together, fastPLS calculates the PLS component path and fold projection once. The two classification heads retain separate predictions, metric paths, and selected component counts. For sufficiently tall SIMPLS classification problems, full-data predictor and class moments are also computed once. Each training fold obtains its moments by subtracting the held-out contribution before applying fold-specific centering and scaling; the PLS model, LDA model, and predictions are still fitted independently within every fold. This avoids repeatedly materializing large training-score matrices without leaking held-out information into model fitting. The same sufficient statistics can assemble eligible PLS-SVD folds through the fold-specific predictor Gram and class cross-product, without materializing the training-score matrix. On the Metal route, small reduced products are evaluated with Accelerate while the large sample-matrix products assigned to Metal remain on the GPU.

cv_opt <- pls.single.cv(
    Xdata = Xreg_train,
    Ydata = Ytrain_reg,
    ncomp = 1:3,
    kfold = 5
)

cv_opt$best_ncomp
#> [1] 1

Interpreting R2Y, Q2Y, and RMSD

For regression models, R2Y, Q2Y, and RMSD answer different questions and should not usually be identical. In pls(), R2Y is training-set R2 and independent-test Q2Y uses the training-response mean in its denominator. For multivariate responses, each response is centered separately before the sums of squares are aggregated. In pls.single.cv(), each held-out prediction is evaluated relative to the corresponding fold-training response means. pls.double.cv() applies the same rule to the outer folds. RMSD is also calculated from held-out predictions and is reported on the response scale, so lower values are better. R2Y is a training-set explained-variance estimate from one additional model fitted on the full dataset; set fit = FALSE to skip this extra fit when only cross-validated performance is needed. For classification, Q2Y is calculated from held-out dummy-coded PLS-DA response scores using fold-training class proportions, accuracy reports decoded-label accuracy, and R2Y is calculated from the full-data PLS-DA fit on the dummy-coded response scores. Dummy-response Q2Y and R2Y are not classification accuracy. The exact convention used by each function is also reported in its metrics$definitions element. RMSD is not used for classification.

Use selection = "R2Y" to select from the full-data fitted-response path or selection = "Q2Y" to select from out-of-fold predictions standardized by each fold’s training-response mean. Selecting R2Y automatically enables and returns the fitted path even if fit = FALSE; because it is a training criterion, Q2Y or another held-out metric is usually preferable for choosing model complexity. The former names "r2" and "q2" are rejected because they did not identify these two different quantities unambiguously.

Classification settings can be selected by accuracy, balanced_accuracy, lift_accuracy, macro_precision, macro_recall, macro_f1, kappa, R2Y, or Q2Y. Use balanced accuracy when unequal class frequencies make majority-class accuracy misleading. Regression settings can be selected by R2Y, Q2Y, RMSD, MAE, MAPE_percent, RPD, Pearson_r, or Spearman_r. These additional aggregate metrics are useful for high-dimensional multivariate responses such as spectra. RMSD, MAE, and MAPE_percent are minimized; the other criteria are maximized. Signed bias and signed MRE_percent remain available from evaluate() but are not tuning criteria because they have no unambiguous one-sided optimization direction. The definitions match evaluate(), and an incompatible choice such as classification accuracy for a regression response raises an error before model fitting.

data.frame(
    ncomp = cv_opt$ncomp,
    training_R2Y = round(cv_opt$R2Y, 3),
    heldout_Q2Y = round(cv_opt$Q2Y, 3),
    heldout_RMSD = round(cv_opt$RMSD, 3)
)
#>         ncomp training_R2Y heldout_Q2Y heldout_RMSD
#> ncomp=1     1        0.742       0.687        3.471
#> ncomp=2     2        0.743       0.673        3.547
#> ncomp=3     3        0.774       0.461        4.555

Permutation-test p-values

fastPLS provides two permutation-test procedures. In pls(), the permutation test is a single train/test procedure: the rows of Xtrain are randomly permuted, the model is refitted, and the permuted test-set Q2Y values are compared with the observed Q2Y values component by component. The returned pval is the corrected Monte Carlo upper-tail value (b + 1) / (B + 1), where b counts successful null fits at least as extreme as observed and B counts successful null fits. It can therefore never be zero. Because pls() has no grouping argument, one training row is the permutation unit. The full permutation table is stored in permutation and can be visualized with plot.permutation(), where the x-axis is the correlation between the original and permuted response structure and the y-axis shows R2 and Q2.

perm_fit <- pls(
    Xtrain = Xreg_train,
    Ytrain = Ytrain_reg,
    Xtest = Xreg_test,
    Ytest = Ytest_reg,
    ncomp = 2,
    fit = TRUE,
    perm.test = TRUE,
    return_variance = FALSE,
    seed = 108
)

perm_fit$pval
#>    ncomp=2 
#> 0.00990099
plot.permutation(perm_fit, ncomp = 2)

In pls.double.cv(), the permutation test repeats the complete nested cross-validation workflow. Independent rows are exchanged individually. When constrain identifies repeated observations, complete constraint blocks are exchanged only between groups with the same number of rows. Thus group sizes, within-group response structure, and class frequencies are preserved exactly. The observed and permuted analyses use identical outer and inner folds and the same randomized-SVD seeds, so null variation reflects the exchange operation rather than new folds or sketches. Model selection and permutation inference use the same selection. For imbalanced classification, selection = "balanced_accuracy" therefore selects the model by mean class-specific recall and tests that statistic against its permutation distribution; it does not substitute dummy-response Q2. The output records the statistic in permutation_metric, permutation_observed, and permutation_sampled. Larger predictive metrics use the upper permutation tail, whereas losses such as RMSD use the lower tail. Both tails use (b + 1) / (B + 1). Failed null fits are stored in permutation_errors, omitted from B, and summarized by permutation_completed and permutation_failed.

dcv_perm <- pls.double.cv(
    Xdata = Xreg_train[1:20, ],
    Ydata = Ytrain_reg[1:20],
    ncomp = 1:2,
    kfold_inner = 2,
    kfold_outer = 2,
    perm.test = TRUE,
    seed = 109
)

data.frame(
    permutation_metric = dcv_perm$permutation_metric,
    observed = dcv_perm$permutation_observed,
    p_value = dcv_perm$p.value,
    completed = dcv_perm$permutation_completed,
    failed = dcv_perm$permutation_failed
)
#>   permutation_metric observed    p_value completed failed
#> 1               RMSD 3.293571 0.00990099       100      0

For an imbalanced classification analysis, the corresponding call is:

dcv_balanced <- pls.double.cv(
    Xdata = Xtrain,
    Ydata = Ytrain_cls,
    ncomp = 1:5,
    kfold_inner = 5,
    kfold_outer = 5,
    constrain = patient_id,
    selection = "balanced_accuracy",
    perm.test = TRUE,
    times = 100,
    seed = 109
)
dcv_balanced$balanced_accuracy
dcv_balanced$metrics$permutation

The selected CV object can then be passed directly to pls(). The final model is refitted on the full training set using the selected component count and selected tuning settings, then applied to the independent test set.

cv_select <- pls.single.cv(
    Xdata = Xtrain,
    Ydata = Ytrain_cls,
    ncomp = 1:3,
    kfold = 5,
    seed = 106
)

fit_selected <- pls(
    cv_select,
    Xtest = Xtest,
    Ytest = Ytest_cls,
    return_variance = FALSE
)

data.frame(
    best_ncomp = cv_select$best_ncomp,
    test_accuracy = mean(fit_selected$Ypred[[1]] == Ytest_cls)
)
#>   best_ncomp test_accuracy
#> 1          3     0.8333333

For example, kernelpls can select the best combination of component count and kernel setting. The selected values are returned in best_parameters, and the same CV object can be passed to pls() for the final refit.

cv_kernel <- pls.single.cv(
    Xdata = Xtrain,
    Ydata = Ytrain_cls,
    ncomp = 1:3,
    kfold = 5,
    method = "kernelpls",
    kernel = c("linear", "rbf"),
    gamma = c(0.1, 1),
    seed = 107
)

fit_kernel <- pls(
    cv_kernel,
    Xtest = Xtest,
    Ytest = Ytest_cls,
    return_variance = FALSE
)

data.frame(
    best_ncomp = cv_kernel$best_parameters$ncomp,
    best_kernel = cv_kernel$best_parameters$kernel,
    test_accuracy = mean(fit_kernel$Ypred[[1]] == Ytest_cls)
)
#>   best_ncomp best_kernel test_accuracy
#> 1          3         rbf             1

Nested or double CV

Double cross-validation is a nested validation design for separating model optimization from the final estimate of predictive performance. This separation is especially important for PLS-DA because the number of latent variables and other modelling choices can otherwise be tuned on the same samples used to report performance, producing optimistic accuracy or Q2 estimates. In the terminology of Szymanska et al. (2012), the inner CV loop (CV1) is used to optimize model complexity, such as the number of latent variables, whereas the outer CV loop (CV2) holds out samples that are not used during optimization and therefore provides the performance estimate of the complete modelling strategy. In fastPLS, pls.double.cv() follows this idea using a reproducible grouped- fold plan constructed in R. CPU and Metal use a compiled coordinator for inner selection, outer refitting, prediction, and metric accumulation. CUDA uses an R coordinator around native CUDA single-CV and outer-fit kernels. Both kfold_inner and kfold_outer can be ordinary fold counts or "loocv", and both respect constrain.

cv_double <- pls.double.cv(
    Xdata = Xtrain,
    Ydata = Ytrain_cls,
    constrain = patient_id,
    ncomp = 1:2,
    kfold_inner = 3,
    kfold_outer = 3,
    seed = 104
)

data.frame(
    selected_ncomp_mode = cv_double$bcomp,
    outer_metric = cv_double$metric_name,
    outer_accuracy = cv_double$accuracy,
    outer_Q2Y = cv_double$Q2Y,
    outer_R2Y = cv_double$R2Y
)
#>   selected_ncomp_mode outer_metric outer_accuracy outer_Q2Y
#> 1                   2     accuracy      0.8083333 0.5477306
#>   outer_R2Y
#> 1 0.5684508

The same interface can tune more than the component count. In this example the inner loop chooses the best combination of ncomp, kernel, and gamma for kernelpls. The outer loop then uses the selected combination for each held-out outer fold.

cv_double_kernel <- pls.double.cv(
    Xdata = Xtrain,
    Ydata = Ytrain_cls,
    constrain = patient_id,
    ncomp = 1:2,
    kfold_inner = 3,
    kfold_outer = 3,
    method = "kernelpls",
    kernel = c("linear", "rbf"),
    gamma = c(0.1, 1),
    seed = 108
)

cv_double_kernel$results[[1]]$best_parameters
#> [[1]]
#> [[1]]$ncomp
#> [1] 2
#> 
#> [[1]]$kernel
#> [1] "rbf"
#> 
#> [[1]]$gamma
#> [1] 1
#> 
#> 
#> [[2]]
#> [[2]]$ncomp
#> [1] 2
#> 
#> [[2]]$kernel
#> [1] "rbf"
#> 
#> [[2]]$gamma
#> [1] 1
#> 
#> 
#> [[3]]
#> [[3]]$ncomp
#> [1] 2
#> 
#> [[3]]$kernel
#> [1] "rbf"
#> 
#> [[3]]$gamma
#> [1] 1

The most commonly used pls.double.cv() output fields are:

Important elements returned by pls.double.cv() are:

• results: one entry per repeated outer CV run. Each entry stores Ypred and pred for that run, the outer fold assignment, the best_ncomp selected inside each outer fold, the full best_parameters selected inside each outer fold, the complete inner-CV objects in inner, the run-level metric_name and metric_value, and the fitted backend and method.
• Ypred: the final cross-validated prediction for each sample. For classification, repeated runs are combined by voting; for regression, predictions are averaged across repeated runs.
• acc_tot: classification-only text summary of the total number and percentage of correctly classified samples.
• conf: classification-only confusion matrix. Entries are printed as counts and column percentages so that class-wise errors can be inspected.
• vote_counts: classification-only matrix with one row per sample and one column per class, showing how many repeated outer-CV runs voted for each class.
• accuracy, Q2Y, RMSD, and R2Y: one value per repeated outer CV run. For classification, accuracy reports decoded-label accuracy, Q2Y reports held-out Q2 on dummy-coded PLS-DA response scores, R2Y reports the mean training-fit R2 of the selected outer-fold PLS-DA models, and RMSD is not used. For regression, Q2Y reports held-out Q2, RMSD reports held-out RMSD, and R2Y reports the mean training-fit R2 of the selected outer-fold models.
• metric_name: the metric used for inner/outer model selection. It is based on held-out predictions except when selection = "R2Y" explicitly requests the fitted-response criterion.
• medianQ2Y, CI95Q2Y, medianR2Y, CI95R2Y, medianRMSD, and CI95RMSD: summaries across repeated outer CV runs, returned only when runn > 1. They are omitted for the default runn = 1 output.
• bcomp: the most frequently selected number of components across all outer folds and repeated runs.
• backend and method: the default backend and PLS method used by the call. If vector-valued methods or backends are tuned in the inner loop, the selected fold-level values are stored in results[[run]]$best_parameters.
• selection_metric: the criterion used by the inner CV loop. The default "auto" means accuracy for classification and RMSD-style prediction error for multivariate regression. All task-compatible selection criteria listed above can also drive the nested permutation test.

Singular-vector backends

The PLS algorithms in fastPLS are separated from the inner singular-vector backend. The compiled implementation uses native randomized SVD, estimating the leading singular triplets needed by the requested PLS components. The randomized SVD backend follows the randomized range-finder framework described by Halko, Martinsson, and Tropp (2011): it draws a Gaussian test matrix, forms a low-dimensional sketch of the target matrix, optionally applies power iterations, orthonormalizes the sketch, and solves the final SVD in the small projected space. In pls() and the cross-validation helpers, these controls are supplied directly through ...; in stand-alone SVD calls they are supplied directly to fastsvd().

Standalone accelerator fastsvd() routes are unavailable because their reduced QR and SVD stages are not fully device-native for every matrix shape. The function raises an error rather than returning a partial result. CUDA acceleration remains device resident within PLS, while Metal acceleration uses the fixed CPU/Metal operation split described below.

The PLS Metal backend, pls(..., backend = "metal"), supports float32 PLS-SVD, SIMPLS, OPLS, and kernel PLS with Metal Performance Shaders. It returns host-accessible model factors and always uses a fixed mathematical division: Metal performs fitting products involving the training sample matrix, including explicit cross-covariance formation, fused score/loading products, and implicit randomized range-finder products. Where the implicit transpose operator requires consecutive X V and Y^T (X V) products, both are encoded in one command buffer and synchronized once before the CPU response-centering correction. CPU code performs the reduced decomposition, sequential PLS algebra, LDA, prediction, and model assembly. Dataset dimensions never select a complete CPU or Metal route. The executed path, storage decision, and randomized controls are reported in diagnostics$metal_operation_split. Ordinary R double matrices produce an explicit error rather than silently falling back to the CPU; use float::fl() to request Metal float32 execution. Repeated fits in one R session reuse initialized device and pipeline state, whereas an isolated first call includes that setup cost.

For CUDA builds, GPU execution is selected with backend = "cuda", so the numerical fitting and prediction stages are performed by resident GPU-native PLS routines after input transfer. When the cross-covariance \(S = X^T Y\) would be large, selected backends can use the matrix-free xprod route. In that route, the code avoids explicitly storing S and evaluates the products required by the SVD through identities such as \(Sv = X^T(Yv)\) and \(S^Tu = Y^T(Xu)\). For classification problems with many classes, label-aware products can avoid materializing a dense one-hot response matrix. Compact prediction instead retains the PLS latent factors and evaluates only the requested test-score blocks, avoiding dense coefficient and prediction paths when they are unnecessary. fastPLS chooses the matrix-free route automatically only when the cross-product is large enough to justify it.

For pls.single.cv(), seeded grouped fold construction and the numerical CV path are compiled. For pls.double.cv(), CPU and Metal nested loops are compiled after R constructs the fixed fold plan; CUDA retains R-level nested coordination around native CUDA kernels. Final R-object assembly remains host-side. CUDA SIMPLS CV uploads the task once and keeps fold gathering, training-only standardization, fitting, prediction, and metric reduction on the device. CUDA PLS-SVD uses the same resident strategy for large multivariate regression responses; classification and smaller responses use a fold-local CUDA route because it is faster for those workloads. Metal applies the fixed CPU/Metal operation split used by pls(). Unsupported accelerator combinations stop rather than falling back to a CPU fit.

PLS-SVD

method = "plssvd" implements the direct cross-covariance SVD formulation of PLS, within the same broad PLS regression and PLS-DA framework described by Wold, Sjostrom, and Eriksson (2001), Barker and Rayens (2003), and Boulesteix and Strimmer (2007). After preprocessing, the algorithm computes the dominant singular subspace of \(X^T Y\) and reuses that subspace to generate the requested component path. The maximum meaningful number of PLS-SVD components is capped by min(nrow(X), ncol(X), ncol(Y)); if a larger value is requested, the package uses the valid internal rank and records the capped component count. For an implicit cross-covariance, numerical rank is determined on the same singular-value scale as the direct decomposition. When a reduced Gram shortcut does not contain the complete requested retained rank, PLS-SVD uses the direct reduced decomposition rather than silently returning fewer directions.

The main innovation relative to a simple textbook PLS-SVD implementation is the execution strategy. fastPLS combines randomized SVD, compact low-rank prediction factors, CUDA-native fitting, and the matrix-free xprod route. For large classification tasks, the label-aware PLS-SVD path computes class-wise cross-products from the factor labels and avoids dense one-hot Y unless the legacy dense route is explicitly needed.

SIMPLS

method = "simpls" is the optimized fastPLS SIMPLS core. Statistically, it keeps the de Jong (1993) SIMPLS idea: one supervised latent direction is appended at a time, the score and loading vectors are computed, an orthogonal loading direction is constructed, and the cross-covariance is deflated before the next component.

At each refresh, fastPLS extracts fresh candidates from the current deflated cross-covariance. Depending on the active compiled route, a refresh produces one direction or a small candidate block. Accepted candidates are always applied one at a time through the standard SIMPLS orthogonalization and deflation equations, and the next refresh uses the resulting deflated operator. CUDA keeps the SIMPLS state on the device. PLS-SVD uses 32 oversampling directions and five power iterations. Accelerated SIMPLS, OPLS, and kernel-PLS use 32 directions and five iterations for ordinary matrix shapes; when the explicit cross-covariance would exceed 512 MiB, their automatic massive-matrix profile requests oversample = 12 and power = 1, while the executed controls and refresh width are recorded separately. Explicit controls supplied through ... take precedence.

The accepted direction is normalized, converted into a SIMPLS component, orthogonalized, and followed by SIMPLS deflation. A candidate block is consumed before the next sketch is generated. The innovation is therefore not a claim of exact de Jong estimator preservation: it is a backend-specialized randomized sequential estimator that trades strict component equivalence for lower computational cost while preserving the supervised deflation structure.

The retained optimizations reorganize work without altering this sequence. The cached-deflation update computes \(v^T S\) once and applies \(S <- S - v(v^T S)\). Coefficients are updated as \(B_a = B_{a-1} + r_a q_a^T\), and fitted values are updated similarly when requested. Cross-products can be cached for suitable matrix shapes, prediction can use compact latent factors, and matrix-free products avoid materializing a large cross-covariance. PLS-SVD instead uses a one-shot cross-covariance SVD; OPLS and kernel PLS reuse the same sequential SIMPLS core where applicable.

The complete component path is also a standard feature of reference SIMPLS software such as pls::simpls.fit, which returns coefficient and fitted-value arrays for components 1 through ncomp. The fastPLS contribution is therefore not the existence of that path, but its compiled incremental execution, compact prediction, matrix-free products, and backend-specific storage strategy. Incremental updates avoid rebuilding \(R_a Q_a^T\) and \(T_a Q_a^T\) from scratch at requested prefixes; compact prediction stores latent factors instead of all dense coefficient prefixes; conditional \(X^T X\) caching can benefit tall matrices; and matrix-free products trade repeated operator applications for avoiding an explicit \(p\)-by-\(q\) cross-covariance. Because these trade-offs depend on \(n\), \(p\), \(q\), and the requested component set, they are memory and execution optimizations rather than a claim of universal acceleration.

OPLS

method = "opls" follows the orthogonal PLS idea introduced by Trygg and Wold (2002). It first removes response-orthogonal variation from X and then fits the optimized PLS core to the filtered matrix. The number of orthogonal components is controlled by north. The filtering step estimates orthogonal scores and loadings, subtracts the corresponding variation from X, and stores the orthogonal filter so the same transformation is applied at prediction time. With float32 input, stored OPLS matrices use single precision. CUDA performs OPLS in float32 or float64, while Metal performs it in float32. CUDA keeps orthogonal filtering, the predictive core, and prediction on device. Metal uses GPU sample-matrix products and CPU orthogonal filtering, reduced factorizations, sequential PLS updates, and prediction. Users select OPLS with method = "opls"; the lower-level implementation details are handled internally.

Kernel PLS

method = "kernelpls" follows the kernel PLS formulation of Rosipal and Trejo (2001). It constructs a kernel representation before applying the selected inner PLS model. The supported kernels are "linear", "rbf", and "poly". For nonlinear kernels, the training kernel matrix is centered and the PLS model is fitted in kernel space. For the linear kernel, fastPLS avoids unnecessary kernel materialization and uses the ordinary linear PLS path. For RBF and polynomial kernels, CUDA constructs, centers, and retains the Gram representation on device through fitting and prediction. Metal evaluates Gram and sample-matrix products through Metal and completes reduced factorizations, sequential updates, and prediction on the CPU. Because this representation requires quadratic storage in the number of training observations, an explicit memory guard rejects unsafe workloads before device allocation.

Classification heads

For classification, factor responses are treated as PLS-DA targets. The following heads are classification-specific decoders and are not used for numeric regression responses. The default classifier, classifier = "argmax", predicts the class with the largest predicted dummy-response score. This is the classical PLS-DA decision rule described in chemometric discrimination work such as Barker and Rayens (2003) and reviewed by Boulesteix and Strimmer (2007).

The LDA head, classifier = "lda", fits a regularized linear discriminant model on the PLS latent scores, following the linear discriminant principle of Fisher (1936). For each class, fastPLS estimates a score-space mean, an empirical prior, and a common pooled within-class covariance matrix. If T contains the training scores, the pooled covariance is computed without materializing class-centered score blocks as (T^T T - sum_c n_c mu_c mu_c^T) / max(1, n - C). The class systems are solved by Cholesky factorization and triangular solves; fastPLS neither explicitly inverts the covariance nor uses Gauss-Jordan elimination. Numerical recovery is deterministic. With s = trace(Sigma) / q (or one when this scale is invalid), the implementation tries lambda = rho * s for rho = 1e-8, 1e-6, 1e-5, 1e-4, 1e-3, 1e-2, increasing the value only when Cholesky fails. Prediction uses the standard discriminant t^T w_c - 0.5 mu_c^T w_c + log(n_c / n). Float32 inputs keep LDA moments, factorization, coefficients, and prediction buffers in float32. CUDA retains the score-space covariance, Cholesky factors, triangular-solve workspaces, coefficients, and prediction buffers on device across component prefixes. The Metal operation split performs LDA on the CPU using the host-accessible latent scores returned by the fitting route.