fastPLS
Fast Partial Least Squares for High-Dimensional Data

Package index
Search the fastPLS package
Vignettes
Functions
559
Source code
16
Man pages
14
Home
/
CRAN
/
fastPLS
/
fastPLS User Guide

fastPLS User Guide
In fastPLS: Fast Partial Least Squares for High-Dimensional Data 

knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
options(width = 68)
library(fastPLS)

Overview

fastPLS combines partial least squares (PLS) modelling with compiled numerical implementations and optional accelerator execution. This vignette is a practical guide to choosing a model, compiling the package, fitting and evaluating models, and running single or nested cross-validation. Detailed derivations and pseudocode are reserved for the accompanying software-methods manuscript.

Choosing a Model

fastPLS provides four related PLS families through one fitting interface. The choice should follow the scientific question and the size of the data rather than the execution backend. The direct PLS-SVD method uses singular value decomposition (SVD) to obtain its latent directions.

| Method | Use it when | Main consideration | |:--|:--|:--| | plssvd | A direct low-rank model is suitable, especially for multivariate responses. | One decomposition supplies every requested component prefix. | | simpls | A general linear PLS model is required for regression or classification. | This is the optimized SIMPLS-family estimator used by the package. | | opls | Response-orthogonal variation in the predictors should be separated before prediction. | The orthogonal filter adds sequential work and must be applied to new observations. | | kernelpls | A nonlinear relation is plausible and the number of observations is moderate. | Nonlinear kernels require a training Gram matrix with quadratic storage in the sample count. |

Classification can decode PLS scores by the largest dummy-response score (classifier = "argmax") or by linear discriminant analysis (LDA) in the retained score space (classifier = "lda"). LDA is the default. It can improve discrimination when class boundaries are not well represented by the largest raw PLS score.

All public methods use the package's native randomized SVD. It is an approximation, so the fitted object records the seed and effective numerical controls.

Installation and Compilation

Install the released package with:

install.packages("fastPLS")

To compile the development version from source, install remotes and use a fresh R session:

install.packages("remotes")
remotes::install_github(
    "tkcaccia/fastPLS",
    upgrade = "never",
    force = TRUE,
    build_vignettes = TRUE
)

The package always builds a CPU backend. CUDA and Metal are optional and are included only when their toolchains are available at compilation time. An unavailable accelerator produces an error when requested; it does not silently fall back to the CPU.

macOS

Install Apple's command-line developer tools before a source build:

```{sh macos-command-line-tools, eval = FALSE} xcode-select --install

macOS uses Apple Accelerate for CPU BLAS/LAPACK. Metal is detected automatically
on supported Apple systems. To require Metal and stop installation if it cannot
be compiled, set the following before the source installation:

```r
Sys.setenv(FASTPLS_USE_METAL = "1")

In the benchmarks conducted on the Apple M3 system used for package testing, the Metal backend did not provide a meaningful computational-speed improvement over the macOS CPU backend. Performance may differ with the matrix dimensions, model configuration, and Apple hardware generation.

Ubuntu and Debian

Install the compiler toolchain and OpenBLAS development files before building the package:

```{sh ubuntu-openblas, eval = FALSE} sudo apt update sudo apt install build-essential gfortran pkg-config libopenblas-dev

Then require OpenBLAS during the source build. This prevents an unnoticed
fallback to the BLAS/LAPACK supplied by R:

```r
Sys.setenv(FASTPLS_USE_OPENBLAS = "1")
install.packages("fastPLS", type = "source")

Ubuntu and Debian releases may provide different OpenBLAS versions, compilation options, threading implementations, and CPU kernels. Consequently, identical fastPLS calls can have substantially different runtimes across distributions or OpenBLAS installations without indicating a change in the statistical model. For reproducible timing, use a current build compiled for the target processor and record the version and selected kernel reported by fastPLS_blas(). Detecting the OpenBLAS family alone does not establish that an architecture-appropriate kernel is active.

Fedora

Install the corresponding development packages, then use the same R command shown above:

```{sh fedora-openblas, eval = FALSE} sudo dnf install gcc gcc-c++ gcc-gfortran make pkgconf-pkg-config \ openblas-devel

### Windows x86-64

Install the Rtools release that matches the installed R version. A source build
can use R's BLAS/LAPACK, but OpenBLAS is recommended for the Linux and Windows
performance routes evaluated with fastPLS. One option is to install an x86-64
OpenBLAS archive from an MSYS2 UCRT64 terminal:

```{sh windows-openblas, eval = FALSE}
pacman -S --needed mingw-w64-ucrt-x86_64-openblas

Point the package configuration to the matching prefix and require OpenBLAS:

Sys.setenv(
    FASTPLS_USE_OPENBLAS = "1",
    OPENBLAS_ROOT = "C:/msys64/ucrt64"
)
remotes::install_github(
    "tkcaccia/fastPLS",
    upgrade = "never",
    force = TRUE,
    build_vignettes = TRUE
)

The OpenBLAS library must match the target architecture. In particular, a Windows ARM64 build must not use an x86-64 Rtools or MSYS2 archive. If a compatible ARM64 OpenBLAS development package is unavailable, leave automatic detection enabled and allow the build to use R's BLAS/LAPACK instead.

Windows OpenBLAS distributions can also differ in version, compiler options, threading implementation, and processor-specific kernel. These differences can materially affect runtime even when the fastPLS version and model settings are unchanged. For reproducible benchmarks, record the OpenBLAS version, selected kernel, thread count, and library path reported by fastPLS_blas().

If OpenBLAS is installed in a nonstandard location on Linux or Windows, set OPENBLAS_ROOT to its installation prefix. Set FASTPLS_USE_OPENBLAS = "0" only when an R-supplied BLAS build is explicitly desired.

Verify the Compiled Libraries

After installation, restart R and verify the selected CPU library and optional accelerators:

library(fastPLS)

fastPLS_blas()
has_cuda()
has_metal()

fastPLS_blas() reports the backend together with the OpenBLAS version, configuration, selected CPU core, parallel runtime, active thread count, and resolved library path when available. Linux and Windows benchmarks should proceed only when fastPLS_blas()$backend is "OpenBLAS" and the expected version and core are active. Use fastPLS_blas(details = FALSE) when only the former scalar backend name is required. The campaign tools in fastPLS-extra validate the detailed report before running fastPLS timing stages. The separate installation vignette contains CUDA toolkit requirements, environment variables, and troubleshooting for architecture or linker errors.

Backend and CPU Configuration

CPU is the default backend. After choosing a mathematical model, users can set one execution backend for the R session when several calls should use the same hardware:

options(backend = "cuda")

fit <- pls(X, y, method = "simpls")    # CUDA session default
fit.cpu <- pls(X, y, backend = "cpu")  # explicit override

Batch jobs may use FASTPLS_BACKEND=cpu|cuda|metal. Precedence is an explicit function argument, options(backend = ...), FASTPLS_BACKEND, and finally CPU. For CPU execution, options(n.cores = 4L) requests four threads from the linked BLAS/OpenMP runtime. Eligible matrix operations may use these threads, but sequential PLS deflation remains serial and multicore acceleration depends on matrix shape and the installed numerical library. An explicit n.cores = argument takes precedence over the option. The first operation that resolves an unavailable option or environment value raises an error; CPU is never substituted. Prediction must use the backend that fitted the model; therefore, a model fitted with CUDA or Metal requires the same explicit or session-level backend selection during prediction.

Supported Computational Routes

The terms in this table distinguish implementation coverage and data residency. Tested denotes a route covered by fixed-seed numerical tests. Approximate identifies rSVD execution. Native means that preprocessing, cross-products, decomposition, component updates, prediction, and an optional LDA head execute on the selected device after input transfer. Unavailable combinations stop; they do not silently fall back to CPU.

The table uses T for tested, A for approximate, N for device-native, H for an explicit host/device hybrid, and U for unavailable.

| Model or operation | CPU | CUDA | Metal | |---|---:|---:|---:| | PLS-SVD, float64 | T/A | T/A/N | U | | SIMPLS family, float64 | T/A | T/A/N | U | | OPLS, float64 | T/A | T/A/N | U | | Linear kernel PLS, float64 | T/A | T/A/N | U | | Nonlinear kernel PLS, float64 | T/A | T/A/N | U | | PLS-SVD, float32 | T/A | T/A/N | T/A/H | | SIMPLS family, float32 | T/A | T/A/N | T/A/H | | OPLS, float32 | T/A | T/A/N | T/A/H | | Linear kernel PLS, float32 | T/A | T/A/N | T/A/H | | Nonlinear kernel PLS, float32 | T/A | T/A/N | T/A/H | | Argmax PLS-DA | T | T/N | T/H | | Latent-space LDA | T | T/N | T/H |

Float32 support is route and platform dependent, as described in the Float32 Input section. CUDA uses a device-native route, whereas Metal uses the fixed CPU/Metal operation split. Metal float64 requests stop because the hardware does not provide native double-precision arithmetic. Nonlinear kernel workloads that exceed the guarded device-memory budget stop before allocation and never fall back to another method or backend.

Public API and Data Inputs

Users select the mathematical method and implementation through the public functions; the internal C++, CUDA, Metal, and benchmarking helpers are not called directly.

| Function | Purpose | |---|---| | pls() | Fit PLS models for regression or classification. | | predict() | Predict from fitted fastPLS, OPLS, or kernel PLS models. | | plot() | Plot PLS scores and optional ellipses. | | plot.permutation() | Plot R2/Q2 diagnostics from a PLS permutation test. | | pls.single.cv() | Select components by grouped cross-validation. | | pls.double.cv() | Run nested cross-validation. | | evaluate() | Evaluate classification or regression predictions. | | fastsvd() | Run the stand-alone native CPU randomized SVD. | | fastcor() | Compute fast Pearson-style correlations. | | ViP() | Compute variable-importance-in-projection trajectories. | | fastPLS_blas() | Report the CPU library, version, configuration, and runtime core. | | has_cuda() | Check whether CUDA-native fastPLS support is available. | | has_metal() | Check for Apple Metal support. |

Choosing a Starting Model

| Goal | Suggested call | |---|---| | Fast exploratory modelling | SIMPLS-family estimator with CPU rSVD | | Repeatable CPU analysis | backend = "cpu", seed = 1 | | Direct cross-covariance PLS-SVD | method = "plssvd" | | Remove response-unrelated structured variation | method = "opls" | | Nonlinear relationships | Kernel PLS with an RBF or polynomial kernel | | Classification: latent-space discriminant analysis (default) | classifier = "lda" | | Classification: argmax PLS-DA decoding | classifier = "argmax" | | Small or moderate data | backend = "cpu" | | Large dense matrices | CUDA when available | | Apple Silicon acceleration | Metal when available | | Float32 CPU input | float::fl(X), backend = "cpu" |

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.

Classification Tasks

For classification, responses are supplied as factors. fastPLS handles the PLS-DA response encoding internally and returns predicted class labels. The classifier argument is used only for this type of task and selects the classification head: argmax or lda. These options are not regression models and are not used for numeric responses. The examples in this section use iris for compact multiclass classification.

set.seed(100)
X <- as.matrix(iris[, 1:4])
Y_cls <- iris$Species
cls_test_id <- sample(seq_len(nrow(X)), 30)
Xtrain <- X[-cls_test_id, , drop = FALSE]
Xtest <- X[cls_test_id, , drop = FALSE]
Ytrain_cls <- Y_cls[-cls_test_id]
Ytest_cls <- Y_cls[cls_test_id]

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 = "lda", which fits linear discriminant analysis in the retained PLS score space. Use classifier = "argmax" explicitly to predict 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

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
pred_cls_later$metrics$metrics
head(pred_cls_later$Ypred_top[["ncomp=2"]])

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

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

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

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
)

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)

evaluate(
    observed = Ytest_cls,
    predicted = fit_cls_plssvd$Ypred[["ncomp=2"]]
)$confusion

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)

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

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 following example fits OPLS to the three-class iris split defined at the beginning of this section, with two orthogonal components (north = 2). The same plot() interface is available for PLS-SVD, the SIMPLS-family estimator, OPLS, and kernel PLS.

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)

Regression Tasks

For regression, responses are supplied as numeric vectors or matrices. The examples in this section use mtcars to demonstrate univariate regression, multivariate regression, and OPLS. Regression does not use classifier = "argmax" or classifier = "lda": those are classification heads for factor responses. For numeric responses, pls() returns continuous predictions and regression metrics such as R2Y, Q2Y, and RMSD.

set.seed(100)
Xreg <- as.matrix(mtcars[, c("disp", "hp", "wt", "qsec", "drat")])
Y_reg <- mtcars$mpg
reg_test_id <- sample(seq_len(nrow(Xreg)), 8)
Xreg_train <- Xreg[-reg_test_id, , drop = FALSE]
Xreg_test <- Xreg[reg_test_id, , drop = FALSE]
Ytrain_reg <- Y_reg[-reg_test_id]
Ytest_reg <- Y_reg[reg_test_id]

Univariate Regression

For a univariate task, supply the response as a numeric vector. The simplest workflow is to provide both the training data and an independent test set directly to pls(). The fitted object then contains test-set predictions and, when 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

Float64 and Float32

Standard R numeric matrices use float64, whereas a float::float32 object requests float32 execution. Float64 stores each value in eight bytes and retains approximately 15--16 decimal digits of precision. Float32 stores each value in four bytes and retains approximately seven decimal digits. Float32 can therefore halve the representation size of a matrix and reduce data-transfer and memory-bandwidth costs, but it also introduces more rounding error.

| Property | float64 | float32 | |:--|:--|:--| | R input | Standard numeric matrix | float::fl() matrix | | Storage per value | 8 bytes | 4 bytes | | Approximate decimal precision | 15--16 digits | 7 digits | | fastPLS backends | CPU and CUDA | CPU, CUDA, and Metal | | Suggested role | Baseline and confirmatory analyses | Reduced-storage or accelerated analyses after numerical comparison |

Reduced representation size does not guarantee a faster fit or lower peak process memory. Temporary workspaces, decomposition costs, host-device transfers, and backend-specific kernels can dominate the calculation. Runtime, incremental host memory, and device memory may therefore increase or decrease when float32 replaces float64.

The lower precision can affect nearly tied latent directions, selected component counts, predictions, or classifications close to a decision boundary. For a new scientific workflow, compare float32 and float64 using the same split, seed, component grid, and model settings, and examine prediction differences as well as the final performance metric. Float64 remains the baseline precision when small numerical differences could alter the scientific conclusion.

CUDA keeps PLS-SVD, the SIMPLS-family estimator, 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.

The package warns once for measured-risk regimes. These include precision-sensitive SIMPLS-family and linear kernel-PLS classification, nonlinear kernels, and multivariate regression with at least 10,000 response columns and 50 components. In the latter 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, the SIMPLS-family estimator, 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

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
pred_reg_later$metrics$metrics
head(pred_reg_later$Ttest)

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
eval_reg$metrics
lapply(eval_reg$metric_definitions, strwrap, width = 56)
eval_reg$per_response

Multivariate Regression

For multivariate regression, supply a numeric response matrix with one column per outcome. The following example predicts mpg, qsec, and drat simultaneously from a common predictor matrix. bycol = TRUE requests both the aggregate multivariate metrics and a response-wise summary.

Xmulti <- as.matrix(mtcars[, c("disp", "hp", "wt", "gear", "carb")])
Ymulti <- as.matrix(mtcars[, c("mpg", "qsec", "drat")])
Xmulti_train <- Xmulti[-reg_test_id, , drop = FALSE]
Xmulti_test <- Xmulti[reg_test_id, , drop = FALSE]
Ymulti_train <- Ymulti[-reg_test_id, , drop = FALSE]
Ymulti_test <- Ymulti[reg_test_id, , drop = FALSE]

fit_multi <- pls(
    Xtrain = Xmulti_train,
    Ytrain = Ymulti_train,
    Xtest = Xmulti_test,
    Ytest = Ymulti_test,
    ncomp = 1:3,
    method = "plssvd",
    bycol = TRUE,
    return_variance = FALSE,
    seed = 102
)

fit_multi$metrics$test[["ncomp=3"]]$metrics
fit_multi$metrics$test[["ncomp=3"]]$per_response

Component selection uses the same interface. With a numeric response matrix, the default criterion is aggregate held-out RMSD across all response columns.

cv_multi <- pls.single.cv(
    Xdata = Xmulti_train,
    Ydata = Ymulti_train,
    ncomp = 1:3,
    kfold = 3,
    method = "plssvd",
    fit = FALSE,
    bycol = TRUE,
    return_splits = TRUE,
    seed = 102
)

cv_multi$best_ncomp
cv_multi$best_metric_value
head(cv_multi$split_index)

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)
fit_opls$Q2Y

Cross-Validation

fastPLS provides two cross-validation helpers. Use pls.single.cv() for grouped k-fold or leave-one-group-out validation; pass a scalar ncomp for a fixed-component CV or a vector of candidates when the number of components should be selected from a grid. Use pls.double.cv() for nested validation, where an inner CV chooses the number of components and an outer CV estimates predictive performance. Both helpers support regression and classification and dispatch to compiled CPU or accelerator routes when available. CUDA SIMPLS-family CV keeps task data and fold workspaces on the device. CUDA library handles are created once per CV call and reused across folds rather than reconstructed for every fit. CUDA PLS-SVD CV also uses a resident route for large multivariate regression responses, while classification and smaller-response tasks retain the faster fold-local CUDA route. This choice is automatic and does not change fold assignment, fold-specific preprocessing, requested components, or prediction semantics. For large multivariate responses, predictor and response marginal moments are computed once. Each training fold obtains its own moments by subtracting the held-out contribution before fold-specific centering and scaling. The large cross-covariance remains matrix-free, so this reuse does not introduce a dense predictor-by-response cache. When its cost and storage model is favorable, the compiled CPU engine can also form a bounded response Gram matrix once for wide-response SIMPLS-family problems and extract an exactly double-centred training submatrix for each fold. The Gram kernel is selected between GEMM and SYRK by platform, precision, and matrix shape; a SYRK route retains one triangle until a complete fold matrix is required. Fold-specific preprocessing, seeds, and model fitting remain unchanged. Aggregate regression metrics are computed directly from the native out-of-fold predictions before those predictions are formatted for R.

On CPU and Metal, the complete nested coordinator is compiled when inner selection uses accuracy, balanced accuracy (or its identical macro-recall definition), dummy-response Q2Y, or regression RMSD/Q2Y. Other documented criteria, such as macro F1, MAE, RPD, or correlation, retain R-level outer coordination around compiled single-CV and fitting kernels because their full metric paths are assembled by evaluate(). CUDA likewise uses R-level nested coordination around its native CUDA kernels. These routes return the same documented fold structure, but the fully compiled coordinator avoids repeated R calls and model-object assembly and is therefore preferable when its metric matches the scientific objective.

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-family classifier with five folds.

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

cv_kfold$metrics$cross_validated[["ncomp=2"]]$metrics
head(cv_kfold$split_index)

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. With return_splits = TRUE, split_index has one row per input sample and one column per fold; entries identify whether that sample was used for training or testing. Nested CV uses the same optional field for all outer and inner splits; outer_test marks samples that are unavailable to a given inner CV cycle.

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,
    return_splits = TRUE,
    seed = 104
)

c(n_folds = length(unique(cv_grouped$fold)),
    n_patient_groups = length(unique(patient_id)))

head(data.frame(
    sample_index = seq_len(nrow(Xtrain)),
    patient_id = patient_id,
    cv_grouped$split_index,
    check.names = FALSE
), 8)

patient_rows <- split(seq_along(patient_id), patient_id)
patients_kept_together <- vapply(
    seq_len(ncol(cv_grouped$split_index)),
    function(fold) {
        all(vapply(patient_rows, function(rows) {
            length(unique(cv_grouped$split_index[rows, fold])) == 1L
        }, logical(1L)))
    },
    logical(1L)
)

stopifnot(all(patients_kept_together))
data.frame(
    fold = colnames(cv_grouped$split_index),
    patients_kept_together = patients_kept_together
)

The displayed rows show that the two samples from each patient receive the same training or test label within a fold. The executable stopifnot() check applies this rule to every patient and every fold, so the vignette build fails if grouped splitting ever places samples from one patient on opposite sides of a split.

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

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.

Some grouped or highly imbalanced folds have lower effective rank than the requested component path. In that case, fastPLS fits only the estimable score prefix and repeats its prediction and metric for larger requested counts. The requested path is preserved, while effective_ncomp records the prefix used in each fold. For LDA, a fold with no estimable PLS direction uses its empirical training-class priors and finite log-prior scores; a training fold containing one class predicts that class. These fallbacks are reported through the fold status field rather than surfacing as a matrix-dimension error.

The same rule applies to a direct regression fit. If fewer directions are estimable than requested, effective_ncomp records the usable prefix and later paths repeat the last estimable prediction and coefficient matrix. If the training response is constant, no direction is estimable: predictions equal the training-response mean, coefficient paths are zero, and response-variance metrics with a zero denominator are NA.

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

cv_opt$best_ncomp

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, AUROC, lift_accuracy, macro_precision, macro_recall, macro_f1, kappa, R2Y, or Q2Y. Binary AUROC pools continuous held-out class scores across folds and treats the second factor level as positive. 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)
)

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

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

Cross-validation results do not retain the training matrices. To fit the final model, combine the selected configuration with best_parameters and pass the original training data explicitly to pls(). This keeps the CV object compact and makes the data used for refitting unambiguous.

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

selected <- utils::modifyList(
    cv_select$tuning_config,
    cv_select$best_parameters
)
svd_controls <- selected$svd_dots
selected$svd_dots <- NULL
fit_selected <- do.call(
    pls,
    c(
        list(
            Xtrain = Xtrain,
            Ytrain = Ytrain_cls,
            Xtest = Xtest,
            Ytest = Ytest_cls,
            return_variance = FALSE
        ),
        selected,
        svd_controls
    )
)

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

For example, kernelpls can select the best combination of component count and kernel setting. The selected values are returned in best_parameters and can be combined with tuning_config for the final explicit 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
)

selected_kernel <- utils::modifyList(
    cv_kernel$tuning_config,
    cv_kernel$best_parameters
)
svd_controls <- selected_kernel$svd_dots
selected_kernel$svd_dots <- NULL
fit_kernel <- do.call(
    pls,
    c(
        list(
            Xtrain = Xtrain,
            Ytrain = Ytrain_cls,
            Xtest = Xtest,
            Ytest = Ytest_cls,
            return_variance = FALSE
        ),
        selected_kernel,
        svd_controls
    )
)

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

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. For the eligible selection criteria described above, 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:3,
    kfold_inner = 3,
    kfold_outer = 3,
    method = "simpls",
    classifier = "lda",
    selection = "balanced_accuracy",
    backend = "cpu",
    return_splits = TRUE,
    seed = 104
)

data.frame(
    selected_ncomp_mode = cv_double$bcomp,
    outer_metric = cv_double$metric_name,
    outer_accuracy = cv_double$accuracy,
    outer_balanced_accuracy = cv_double$balanced_accuracy,
    outer_Q2Y = cv_double$Q2Y,
    outer_R2Y = cv_double$R2Y
)

data.frame(
    outer_fold = seq_along(cv_double$results[[1]]$best_ncomp),
    inner_selected_ncomp = cv_double$results[[1]]$best_ncomp
)

head(cv_double$split_index)

Here, all rows from one illustrative patient remain together in both loops. If an inner training fold contains one class, fastPLS predicts that class and uses a finite constant score for every requested component. These predictions remain in the pooled inner metric. If an outer training partition contains one class, the same rule supplies its held-out prediction and the fold is marked as not estimable for discrimination. The returned degenerate_inner_folds, constant_classifier_fallback, minimum class counts, component_selection_informative, and outer_discrimination_estimable fields make these cases explicit.

The same nested workflow accepts a numeric vector for univariate regression or a numeric matrix for multivariate regression. The default inner-loop criterion for both tasks is held-out RMSD; for a response matrix it is aggregated across response columns.

cv_double_uni <- pls.double.cv(
    Xdata = Xreg_train,
    Ydata = Ytrain_reg,
    ncomp = 1:2,
    kfold_inner = 2,
    kfold_outer = 2,
    method = "simpls",
    backend = "cpu",
    seed = 105
)

cv_double_multi <- pls.double.cv(
    Xdata = Xmulti_train,
    Ydata = Ymulti_train,
    ncomp = 1:2,
    kfold_inner = 2,
    kfold_outer = 2,
    method = "plssvd",
    backend = "cpu",
    bycol = TRUE,
    seed = 105
)

data.frame(
    task = c("univariate regression", "multivariate regression"),
    selected_ncomp = c(cv_double_uni$bcomp, cv_double_multi$bcomp),
    outer_RMSD = c(cv_double_uni$RMSD, cv_double_multi$RMSD)
)

For each outer fold, the inner loop selects the component count by balanced accuracy and the selected LDA model predicts samples that were not involved in that choice. Component counts may therefore differ among outer folds; bcomp reports their most frequent value, while balanced_accuracy and Ypred summarize the outer held-out predictions.

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

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

| Field | Meaning | |---|---| | bcomp | Most frequently selected number of components. | | accuracy, Q2Y, RMSD | Outer-CV predictive performance. | | Ypred | Final cross-validated prediction for each sample. | | conf | Classification confusion matrix. | | metrics | Detailed outer-CV evaluation results. | | results | Detailed per-run and per-fold information. |

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

SVD Utility

fastsvd() provides direct access to the native CPU randomized SVD used by PLS. It returns the left singular vectors (u), singular values (d), and right singular vectors (v) for users who want a truncated decomposition outside a PLS model. A standard R numeric matrix selects float64 computation. A matrix created with float::fl() selects the native float32 path automatically; its left and right singular vectors remain in float32 format. The function exposes no solver selector because rSVD is its only algorithm. Stand-alone CUDA and Metal SVD routes are not exposed; their reduced decompositions are not device-native for every matrix shape.

s64 <- fastsvd(Xtrain, ncomp = 3, seed = 104)
s32 <- fastsvd(float::fl(Xtrain), ncomp = 3, seed = 104)
names(s64)
s64$d
inherits(s32$u, "float32")

CUDA and Metal Availability

CUDA and Apple Metal are optional. Use has_cuda() before selecting backend = "cuda" and has_metal() before selecting backend = "metal". An unavailable accelerator request stops with an informative error and is not redirected to CPU.

has_cuda()
has_metal()

if (has_cuda()) {
    fit_gpu <- pls(
    Xtrain,
    Ytrain_cls,
    Xtest,
    Ytest_cls,
    ncomp = 1:2,
    backend = "cuda"
    )

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

if (has_metal()) {
    Xtrain_metal <- float::fl(as.matrix(Xtrain))
    Xtest_metal <- float::fl(as.matrix(Xtest))
    fit_metal <- pls(
        Xtrain_metal,
        Ytrain_cls,
        Xtest_metal,
        Ytest_cls,
        ncomp = 1:2,
        backend = "metal"
    )
}

Helper Functions

fastcor() computes fast Pearson-style correlations. ViP() returns variable importance in projection trajectories for direct SIMPLS-family fits, following the standard VIP interpretation used for PLS variable ranking (Wold, Sjostrom, and Eriksson, 2001; Chong and Jun, 2005). VIP is most useful when the columns of X are interpretable predictors, such as genes, metabolites, spectral bins, or clinical variables; larger values indicate stronger contribution to the fitted latent predictive model. Linear kernel PLS uses the same direct SIMPLS-family path and is also supported. PLS-SVD scores are not generally orthogonal and its stored response singular vectors are not the final latent regression map; OPLS applies an additional predictor filter; and nonlinear kernel weights index training samples rather than original predictors. ViP() therefore stops for these three cases instead of returning a quantity with a misleading variable- importance interpretation. The supported matrix calculation, like the metrics returned by evaluate(), is implemented in the dependency-free C++ core and is independent of CUDA or Metal availability; the R functions validate inputs and format the returned objects.

C <- fastcor(Xtrain, byrow = FALSE, diag = FALSE)
dim(C)

vip <- ViP(fit_reg)
dim(vip)

Methodological Scope

This vignette is a practical guide to fitting, validating, and predicting with fastPLS. The mathematical derivations and executable pseudocode for native rSVD, PLS-SVD, the SIMPLS-family estimator, OPLS, kernel PLS, compact prediction, and sufficient-statistics cross-validation are maintained in the future Journal of Statistical Software manuscript. Keeping those details in the methods article avoids repeating a long technical description in the package vignette while preserving a single place in which the equations can be checked against the shared C++ core.

The fitted object still records the information needed to audit an analysis, including the executed method, backend, precision, randomized controls, effective component count, and route diagnostics. These fields should be saved with the analysis whenever numerical reproducibility is important.

Method References

The implementation in fastPLS is not a line-by-line copy of the papers below; rather, these papers define the statistical algorithms or numerical building blocks that the package implements and accelerates.

Session Information

sessionInfo()


Try the fastPLS package in your browser

Any scripts or data that you put into this service are public.

fastPLS documentation built on Sept. 29, 2026, 1:06 a.m.

R Package Documentation

rdrr.io home R language documentation Run R code online

Browse R Packages

CRAN packages Bioconductor packages R-Forge packages GitHub packages

We want your feedback!

Note that we can't provide technical support on individual packages. You should contact the package authors for that.
GitHub issue tracker
ian@mutexlabs.com
Personal blog

 
Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.

Close