---
title: "fastPLS User Guide"
output: rmarkdown::html_vignette
header-includes:
    - \usepackage{array}
vignette: >
    %\VignetteIndexEntry{fastPLS User Guide}
    %\VignetteEngine{knitr::rmarkdown}
    %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
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:

```{r install-cran, eval = FALSE}
install.packages("fastPLS")
```

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

```{r install-github, eval = FALSE}
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 require-metal, eval = FALSE}
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 require-openblas-linux, eval = FALSE}
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:

```{r require-openblas-windows, eval = FALSE}
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:

```{r verify-compiled-libraries, eval = FALSE}
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:

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

```{r chunk-002}
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.

```{r chunk-003}
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.

```{r chunk-004}
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.

```{r chunk-005}
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.

```{r chunk-005a}
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.

```{r chunk-006}
eval_cls$confusion
```

<div style="page-break-before: always;"></div>

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.

```{r chunk-007}
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.

```{r chunk-008}
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:

```{r chunk-009, eval = FALSE}
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:

```{r chunk-010}
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.

```{r chunk-012}
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.

```{r chunk-013, fig.width = 5, fig.height = 4}
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.

```{r iris-three-class-opls, fig.width = 7.5, fig.height = 4}
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`.

```{r chunk-016}
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.

```{r chunk-017}
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.

```{r chunk-018}
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"]`.

```{r chunk-019, fig.width=4.8, fig.height=4}
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`.

```{r chunk-020}
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.

```{r chunk-021}
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.

```{r multivariate-regression}
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.

```{r multivariate-regression-cv}
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.

```{r chunk-022}
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.

```{r chunk-023}
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.

```{r chunk-024}
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.

```{r chunk-025}
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`.

```{r chunk-026}
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.

```{r chunk-027}
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.

```{r chunk-028}
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`.

```{r chunk-029}
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:

```{r balanced-permutation-example, eval=FALSE}
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.

```{r chunk-030}
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.

```{r chunk-031}
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`.

```{r chunk-032}
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.

```{r nested-regression-examples}
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.

```{r chunk-033}
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:

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

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

```{r chunk-034}
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.

```{r chunk-038}
has_cuda()
has_metal()
```

```{r chunk-039, eval = FALSE}
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"
    )
}
```

```{r chunk-040, eval = FALSE}
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.

```{r chunk-041}
C <- fastcor(Xtrain, byrow = FALSE, diag = FALSE)
dim(C)
```

```{r chunk-042}
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.

- Barker, M. and Rayens, W. (2003). Partial least squares for discrimination.
    *Journal of Chemometrics*, 17, 166-173.
- Boulesteix, A.-L. and Strimmer, K. (2007). Partial least squares: a versatile
    tool for the analysis of high-dimensional genomic data. *Briefings in
    Bioinformatics*, 8, 32-44.
- Chong, I.-G. and Jun, C.-H. (2005). Performance of some variable selection
    methods when multicollinearity is present. *Chemometrics and Intelligent
    Laboratory Systems*, 78, 103-112.
- de Jong, S. (1993). SIMPLS: an alternative approach to partial least squares
    regression. *Chemometrics and Intelligent Laboratory Systems*, 18, 251-263.
- Fisher, R. A. (1936). The use of multiple measurements in taxonomic problems.
    *Annals of Eugenics*, 7, 179-188.
- Geladi, P. and Kowalski, B. R. (1986). Partial least-squares regression: a
    tutorial. *Analytica Chimica Acta*, 185, 1-17.
- Halko, N., Martinsson, P.-G., and Tropp, J. A. (2011). Finding structure with
    randomness: probabilistic algorithms for constructing approximate matrix
    decompositions. *SIAM Review*, 53, 217-288.
- Mevik, B.-H. and Wehrens, R. (2007). The pls package: principal component and
    partial least squares regression in R. *Journal of Statistical Software*,
        18,
    1-23.
- Rosipal, R. and Trejo, L. J. (2001). Kernel partial least squares regression
    in reproducing kernel Hilbert space. *Journal of Machine Learning Research*,
    2, 97-123.
- Szymanska, E., Saccenti, E., Smilde, A. K., and Westerhuis, J. A. (2012).
    Double-check: validation of diagnostic statistics for PLS-DA models in
    metabolomics studies. *Metabolomics*, 8, S3-S16.
        doi:10.1007/s11306-011-0330-3.
- Trygg, J. and Wold, S. (2002). Orthogonal projections to latent structures
    (O-PLS). *Journal of Chemometrics*, 16, 119-128.
- Wold, S., Sjostrom, M., and Eriksson, L. (2001). PLS-regression: a basic tool
    of chemometrics. *Chemometrics and Intelligent Laboratory Systems*, 58,
    109-130.

## Session Information

```{r session-info}
sessionInfo()
```