fastPLS Public API and Implementation

SVD Algorithms

fastPLS separates the PLS model from the singular-vector algorithm used inside the fit. The default svd.method = "rsvd" is a randomized SVD backend 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 numerical rank of the target matrix, 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.

The alternative svd.method = "irlba" is CPU-only and uses the bundled implicitly restarted Lanczos bidiagonalization implementation, following the IRLBA approach of Baglama and Reichel (2005). The package calls the internal C IRLBA operator with a working subspace, maximum iteration count, convergence tolerance, epsilon threshold, and singular-value tolerance. These controls can be supplied through ... in pls() or directly in fastsvd(). The IRLBA backend is also used in a matrix-free form for selected xprod paths, where the algorithm receives functions that compute products with the cross-product operator rather than a precomputed dense matrix.

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.

PLS 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 PLSSVD 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. PLSSVD can also use the matrix-free xprod route for randomized SVD or bundled IRLBA 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 fastPLS implementation improves the execution strategy. Candidate directions can be refreshed from the current deflated cross-covariance using randomized/block-style subspace ideas related to modern randomized low-rank methods (Halko, Martinsson, and Tropp, 2011), then consumed sequentially so that the accepted components still follow the SIMPLS deflation geometry. The C++ core also uses cached rank-one deflation, warm starts from the previous latent direction, optional reorthogonalization of deflation vectors, and incremental prediction updates. 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.

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 PLSSVD, SIMPLS, OPLS, and kernel PLS while sharing the accelerated low-rank machinery.

Classifiers

Classification is handled as PLS-DA followed by one of three 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). classifier = "cknn" reranks centroid-selected candidate classes by nearest-neighbour evidence in the supervised latent space, following the nearest-neighbour principle of Cover and Hart (1967).

Backends

The backend argument selects where the main numerical work is performed. backend = "cpu" is the portable compiled C++/RcppArmadillo backend and supports both rsvd and irlba. It is usually best for small and moderate datasets, low-dimensional problems, and machines without GPU support.

backend = "cuda" uses CUDA-native routines when the package is compiled with CUDA and a compatible NVIDIA runtime is available. CUDA paths perform the main PLSSVD and SIMPLS randomized-SVD matrix products, selected xprod products, low-rank prediction, LDA scoring, and candidate-kNN scoring on the GPU. CUDA is most useful for large dense matrices, high component counts, many classes, or many response variables, where GPU matrix multiplication offsets transfer and kernel-launch overhead.

backend = "metal" provides an Apple Silicon backend. The Metal SVD and PCA paths offload the large randomized-SVD matrix multiplications to Apple Metal Performance Shaders while keeping QR and the final small SVD on CPU. The PLS Metal route extends the same idea to PLSSVD, SIMPLS, OPLS, and kernel PLS by dispatching large cross-products, score projections, deflation matrix products, and prediction products through Metal where implemented. Metal is intended for portability and exploratory acceleration on Apple hardware; CPU remains the most robust fallback for small examples.

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.8333333

Classification models can also be fitted once and predicted later. predict() can optionally return ranked classes by setting top or top5 = TRUE. 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.

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
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 or class-score matrices. For classification, the complete output includes global metrics, per-class metrics, and the confusion matrix.

eval_cls <- evaluate(
  observed = Ytest_cls,
  predicted = fit_cls$Ypred[["ncomp=2"]]
)

eval_cls
#> $task
#> [1] "classification"
#> 
#> $metrics
#>    n  accuracy balanced_accuracy macro_precision macro_recall macro_f1 kappa
#> 1 30 0.8333333         0.8333333       0.8498168    0.8333333 0.829497  0.75
#> 
#> $per_class
#>        class support precision recall        f1
#> 1     setosa      10 1.0000000    1.0 1.0000000
#> 2 versicolor      10 0.8571429    0.6 0.7058824
#> 3  virginica      10 0.6923077    0.9 0.7826087
#> 
#> $confusion
#>             observed
#> predicted    setosa versicolor virginica
#>   setosa         10          0         0
#>   versicolor      0          6         1
#>   virginica       0          4         9
#> 
#> $topk
#> NULL

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          6         1
#>   virginica       0          4         9

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,
  top_k = c(1, 2)
)
#> $task
#> [1] "classification"
#> 
#> $metrics
#>    n accuracy balanced_accuracy macro_precision macro_recall macro_f1 kappa
#> 1 30        1                 1               1            1        1     1
#> 
#> $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

Classification Heads

The same pls() interface exposes three classification-specific heads for factor responses: argmax PLS-DA, latent-space LDA, and candidate-kNN. 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  virginica
#> 2 virginica  virginica
#> 3    setosa     setosa
#> 4 virginica versicolor
#> 5 virginica versicolor
#> 6 virginica versicolor

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

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"
)

When a GPU backend is unavailable, the compiled CPU fallback is:

fit_cls_lda_cpu <- pls(
  Xtrain,
  Ytrain_cls,
  Xtest,
  Ytest_cls,
  ncomp = 1:2,
  method = "plssvd",
  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

The candidate-kNN head is selected with classifier = "cknn". It is most useful for larger multiclass data, but the same syntax works on smaller examples: The main tuning parameters control how strongly the local neighbourhood can modify the global PLS class ranking. k is the number of same-class neighbours used for each candidate class. tau is a positive temperature: small values make the local score behave almost like the single best neighbour, whereas larger values smooth the evidence over the top neighbours. alpha weights the centroid/prototype score that is added to the local kNN score. Larger alpha keeps predictions closer to the global class-centroid ranking; smaller alpha lets local neighbours dominate the reranking. top_m controls how many centroid-ranked classes are passed to this reranker.

fit_cls_cknn <- pls(
  Xtrain,
  Ytrain_cls,
  Xtest,
  Ytest_cls,
  ncomp = 1:2,
  classifier = "cknn",
  k = 5,
  top_m = 3,
  seed = 100
)

evaluate(
  observed = Ytest_cls,
  predicted = fit_cls_cknn$Ypred[["ncomp=2"]]
)$confusion
#>             observed
#> predicted    setosa versicolor virginica
#>   setosa         10          0         0
#>   versicolor      0          6         4
#>   virginica       0          4         6

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.8333333 1.0000000 0.9666667

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 bundled two-class colon gene-expression data. The same plot() call works for method = "plssvd", "simpls", "opls", and "kernelpls"; here we keep the code short and show one fitted model.

data("colon", package = "fastPLS")
set.seed(200)
X_two <- as.matrix(colon$X)
y_two <- colon$y
train_two <- unlist(tapply(seq_along(y_two), y_two, function(i) {
  sample(i, floor(0.7 * length(i)))
}))

fit_two <- pls(
  X_two[train_two, , drop = FALSE], y_two[train_two],
  X_two[-train_two, , drop = FALSE], y_two[-train_two],
  ncomp = 1:2, fit = TRUE, proj = TRUE, seed = 200
)

plot(fit_two, score.set = "train", groups = y_two[train_two],
     ellipse = TRUE, main = "Colon, two classes")

The second score-plot example uses a bundled three-subtype breast-cancer expression example and fits OPLS with two orthogonal components (north = 2). The example uses predefined training and test sets.

data("breast", package = "fastPLS")
X_three <- as.matrix(breast$X_train)
y_three <- breast$y_train
X_three_test <- as.matrix(breast$X_test)
y_three_test <- breast$y_test

opls_three <- pls(
  X_three,
  y_three,
  X_three_test,
  y_three_test,
  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 = y_three,
  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,
  svd.method = "irlba",
  fit = TRUE,
  return_variance = FALSE
)

fit_reg$Q2Y
#>   ncomp=1   ncomp=2   ncomp=3 
#> 0.7505150 0.7535829 0.8601637

Float32 Input

pls() also accepts predictor and response matrices from the float package. When a float::float32 object is supplied, the float32 route keeps centering, scaling, cross-products, low-rank SVD steps, latent scores, and predictions in 32-bit arithmetic. This is currently available for method = "plssvd" or method = "simpls". CPU and Metal support svd.method = "rsvd" or svd.method = "irlba", while CUDA supports float32 randomized SVD. For classification, float32 models support classifier = "argmax", classifier = "lda", and classifier = "cknn". Unsupported combinations stop with an error rather than silently converting the data to double.

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.7505149 0.7535830

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.

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.7505150 0.7535829 0.8601637
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.

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

eval_reg
#> $task
#> [1] "regression"
#> 
#> $metrics
#>   n        R2        Q2     RMSD     RMSE      MAE      bias MRE_percent
#> 1 8 0.8601637 0.8634626 2.295494 2.295494 1.929671 0.6784321    12.48956
#>   MAPE_percent      RPD Pearson_r Spearman_r
#> 1     10.76821 2.858815  0.934081  0.8072875
#> 
#> $per_response
#>   response n        R2        Q2     RMSD     RMSE      MAE      bias
#> 1       Y1 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

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.

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.7535829 0.8601637

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,
  svd.method = "rsvd",
  seed = 103
)

cv_kfold$metrics
#>   ncomp_index metric_name metric_value
#> 1           1    accuracy        0.825

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,
  svd.method = "rsvd",
  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",
  svd.method = "rsvd",
  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, Q2 for univariate regression, and RMSD-style prediction error for multivariate regression outputs. 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.

cv_opt <- pls.single.cv(
  Xdata = Xreg_train,
  Ydata = Ytrain_reg,
  ncomp = 1:3,
  kfold = 5,
  svd.method = "irlba",
  maxit = 1000L
)

cv_opt$best_ncomp
#> [1] 3

Interpreting R2Y, Q2Y, and RMSD

For regression models, R2Y, Q2Y, and RMSD answer different questions and should not usually be identical. Q2Y is calculated by default from held-out cross-validation predictions and measures predictive performance on samples that were not used to fit the corresponding fold model. 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, accuracy reports decoded-label accuracy, and R2Y is calculated from the full-data PLS-DA fit on the dummy-coded response scores. RMSD is not used for classification.

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.683        3.290
#> ncomp=2     2        0.743       0.682        3.292
#> ncomp=3     3        0.774       0.684        3.283

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 empirical fraction mean(Q2Y_permuted > Q2Y_observed), without a +1 correction. 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
#> [1] 0
plot.permutation(perm_fit, ncomp = 2)

In pls.double.cv(), the permutation test repeats the complete nested cross-validation workflow after shuffling the rows of Xdata. The observed statistic is the median outer-CV Q2Y, and Q2Ysampled stores the permuted median values. For metrics where larger is better, the p-value is mean(Q2Y_permuted >= Q2Y_observed); for loss metrics where smaller is better, it is mean(loss_permuted <= loss_observed).

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

data.frame(
  observed_median_Q2Y = median(dcv_perm$Q2Y, na.rm = TRUE),
  p_value = dcv_perm$p.value
)
#>   observed_median_Q2Y p_value
#> 1           0.7161204       0

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,
  svd.method = "rsvd",
  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",
  svd.method = "rsvd",
  kernel = c("linear", "rbf"),
  gamma = c(0.1, 1),
  seed = 107
)

fit_kernel <- pls(
  cv_kernel,
  Xtest = Xtest,
  Ytest = Ytest_cls,
  return_variance = FALSE
)
#> Warning in (function (Xtrain, Ytrain, Xtest = NULL, Ytest = NULL, ncomp = 2, :
#> NAs introduced by coercion

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          2      linear     0.8333333

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: the inner loop calls pls.single.cv() to select the best number of components and any vector-valued predictive hyperparameters, while the outer loop refits the selected model on the corresponding outer training set and predicts the outer test fold. 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,
  svd.method = "rsvd",
  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 outer_R2Y
#> 1                   2     accuracy      0.8333333 0.5709992 0.5808943

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",
  svd.method = "rsvd",
  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] "linear"
#> 
#> [[1]]$gamma
#> [1] NA
#> 
#> 
#> [[2]]
#> [[2]]$ncomp
#> [1] 2
#> 
#> [[2]]$kernel
#> [1] "linear"
#> 
#> [[2]]$gamma
#> [1] NA
#> 
#> 
#> [[3]]
#> [[3]]$ncomp
#> [1] 2
#> 
#> [[3]]$kernel
#> [1] "linear"
#> 
#> [[3]]$gamma
#> [1] NA

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 held-out metric used for inner/outer model selection.
• 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.

Singular-vector backends

The PLS algorithms in fastPLS are separated from the inner singular-vector backend. The compiled CPU implementation supports two truncated SVD routes: irlba and randomized SVD. The irlba backend follows the implicitly restarted Lanczos bidiagonalization strategy of Baglama and Reichel (2005). It estimates only the leading singular triplets needed by the requested PLS components, rather than computing a full dense SVD. 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().

On Apple Silicon builds where has_metal() returns TRUE, fastsvd(backend = "metal", method = "rsvd") and pca(backend = "metal", method = "rsvd") can use the Metal backend. This backend keeps the same randomized SVD mathematics as the CPU backend, but sends the large sketching matrix multiplications to Apple Metal Performance Shaders. The QR decomposition and final small SVD remain on the CPU. Because matrices are copied between R and Metal buffers, this backend is most useful for larger matrices where the GPU multiplication work is large enough to dominate transfer overhead.

The experimental PLS Metal backend, pls(..., backend = "metal"), extends the same idea to PLSSVD, SIMPLS, OPLS, and kernel PLS. The large cross-products, score projections, deflation matrix products, and prediction matrix products are dispatched through Metal. Small decompositions, scalar reductions, and some kernel-centering operations remain CPU-side until dedicated Metal kernels are added.

For CUDA builds, GPU execution is selected with backend = "cuda", so the main matrix operations are performed by GPU-native PLS routines. 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. The low-rank streamed prediction path is inspired by the FlashSVD streaming low-rank inference design of Shao et al. (2025), but is adapted here to the PLS prediction factors exposed by fastPLS. In ordinary use, leave xprod = NULL; fastPLS then chooses the matrix-free route automatically only when the cross-product is large enough to justify it.

PLSSVD

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 PLSSVD 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.

The main innovation relative to a simple textbook PLSSVD implementation is the execution strategy. fastPLS can use truncated IRLBA or randomized SVD, compact low-rank prediction factors, CUDA-native fitting, and the matrix-free xprod route. For large classification tasks, the label-aware PLSSVD 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.

The implementation differs from a naive classical SIMPLS loop in how it obtains new singular directions. Instead of always calling the SVD backend for one fresh vector after every deflation, fastPLS can refresh a small block of candidate directions from the current deflated cross-covariance. The block is consumed sequentially: each candidate direction is normalized, converted into a SIMPLS component, orthogonalized, and followed by the usual SIMPLS deflation before the next component is accepted. The cached-deflation update explicitly computes \(v^T S\) and applies the rank-one update \(S <- S - v(v^T S)\), reducing temporary allocations and repeated projection work. The incremental refresh uses the previous latent direction as a warm start for the next randomized refresh, which can reduce the number of iterations needed to recover the next dominant direction.

These block-refresh, warm-start, and cached-deflation optimizations are SIMPLS-specific because SIMPLS is a sequential deflation method. PLSSVD uses a different one-shot cross-covariance SVD strategy. OPLS and kernel PLS reuse the same optimized compiled building blocks where their internal fitting step is linear PLS-like.

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. 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.

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. A small ridge term scaled by the average covariance diagonal is added for numerical stability. Prediction uses the standard LDA discriminant score and returns the largest-scoring class. The CUDA and Metal routes are selected automatically from backend when available.

The cKNN head, classifier = "cknn", is the short public name for the candidate-kNN classifier used in large multiclass PLS-DA. It first scores classes by centroids in supervised PLS score space, keeps the top top_m candidate classes, and then reranks only those candidates by same-class nearest neighbours in the PLS score space, using the nearest-neighbour classification principle of Cover and Hart (1967) after a supervised PLS compression step. For each candidate class, the neighbour similarities are combined with a temperature-smoothed top-k score. The tau parameter controls this smoothing: lower values emphasize the strongest neighbour, while higher values make the local evidence closer to an average of the top neighbours. The final candidate score is the smoothed local kNN evidence plus alpha times the centroid score, so alpha controls the balance between local correction and the global PLS class-centroid ranking. The top_m parameter is internal to the cKNN decision rule: it controls how many centroid-ranked classes are passed to the kNN reranker. This is different from the prediction argument top, which only controls how many ranked output labels are returned to the user. The default tuning parameters are k = 10, tau = 0.2, alpha = 0.75, and top_m = 20. predict() can also return ranked predictions by using top or top5 = TRUE.

For large multiclass data, cKNN can be memory-intensive because it uses PLS scores for both training and test samples. The cknn_memory argument controls how these score matrices are materialized. cknn_memory = "standard" uses the historical one-pass path. cknn_memory = "blocked" keeps the same classifier but predicts test samples in blocks, avoiding a full test-score matrix in memory. cknn_memory = "streaming" additionally builds the training candidate-score cache in row blocks for scalar ncomp, which is useful for ImageNet-scale feature matrices. The default, cknn_memory = "auto", selects a memory-aware strategy from the size of the latent-score cache. These modes do not change the cKNN scoring rule; they only change when intermediate matrices are created.