---
title: "fastPLS Public API and Implementation"
output: BiocStyle::html_document
vignette: >
  %\VignetteIndexEntry{fastPLS Public API and Implementation}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
```

## Overview

`fastPLS` exposes a compact public API. Users do not need to call the internal
C++, CUDA, or benchmarking helpers directly; all model families and backends are
selected through the user-facing functions below.

| Function | Purpose |
|---|---|
| `pls()` | Fit PLSSVD, SIMPLS, OPLS, or kernel PLS for regression or classification. |
| `predict()` | Predict from fitted `fastPLS`, OPLS, or kernel PLS models. |
| `plot()` | Plot PLS or PCA score maps, with optional confidence or Hotelling T2 ellipses. |
| `plot.permutation()` | Plot R2/Q2 permutation-test diagnostics from `pls(..., perm.test = TRUE)`. |
| `pls.single.cv()` | Run grouped cross-validation and choose the best component count. |
| `pls.double.cv()` | Run nested/double cross-validation with inner component optimization. |
| `evaluate()` | Compute classification and regression performance metrics from predictions. |
| `fastsvd()` | Run a stand-alone truncated SVD through the bundled backends. |
| `pca()` | Compute PCA using the same SVD backend selector. |
| `fastcor()` | Compute fast Pearson-style correlations. |
| `ViP()` | Compute variable-importance-in-projection trajectories. |
| `has_cuda()` | Check whether CUDA-native fastPLS support is available. |
| `has_metal()` | Check whether Apple Metal-native fastPLS support is available. |

A practical starting point is:

| Goal | Suggested call |
|---|---|
| General regression or PLS-DA | `method = "simpls"`, `backend = "cpu"`, `svd.method = "rsvd"` |
| Direct cross-covariance PLS-SVD | `method = "plssvd"` |
| Remove response-unrelated structured variation | `method = "opls"` |
| Try nonlinear sample relationships | `method = "kernelpls"`, `kernel = "rbf"` or `"poly"` |
| Classification: argmax PLS-DA decoding | `classifier = "argmax"` |
| Classification: latent-space discriminant analysis | `classifier = "lda"` |
| Classification: candidate-neighbour classification | `classifier = "cknn"` |
| Small or moderate data | `backend = "cpu"` |
| Large dense matrices or many classes/responses | `backend = "cuda"` when available |
| Apple Silicon exploratory acceleration | `backend = "metal"` when `has_metal()` is `TRUE` |
| Float32 CPU input | `float::fl(X)`, `backend = "cpu"`, `svd.method = "rsvd"` |

## Implementation Summary

### SVD Algorithms

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

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

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

### PLS Algorithms

`method = "plssvd"` implements the direct PLS-SVD formulation, in which the
dominant singular subspace of the predictor-response cross-covariance matrix is
computed once and reused for all requested component counts. This makes PLSSVD
efficient when the response rank is the limiting dimension. The implementation
adds compact latent prediction factors: instead of always storing and applying
a full coefficient array, prediction can be performed through low-rank factors
based on the latent projection and response-side coefficients. PLSSVD can also
use the matrix-free `xprod` route for randomized SVD or bundled IRLBA when the
cross-covariance matrix would be large.

`method = "simpls"` is the optimized SIMPLS implementation and is the default
general-purpose PLS core. The statistical structure follows de Jong's SIMPLS
algorithm (de Jong, 1993): each component is obtained from the current
cross-covariance state, converted into a score and loading direction, and
followed by SIMPLS orthogonalization and deflation. The fastPLS implementation
improves the execution strategy. Candidate directions can be refreshed from the
current deflated cross-covariance using randomized/block-style subspace ideas
related to modern randomized low-rank methods (Halko, Martinsson, and Tropp,
2011), then consumed sequentially so that the accepted components still follow
the SIMPLS deflation geometry. The C++ core also uses cached rank-one deflation,
warm starts from the previous latent direction, optional reorthogonalization of
deflation vectors, and incremental prediction updates. For low-dimensional
high-sample problems, the code can cache cross-products such as `X^T X` and
`X^T Y`; for large response or predictor dimensions, it can instead use
matrix-free `xprod` products.

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

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

### Classifiers

Classification is handled as PLS-DA followed by one of three decoding heads.
`classifier = "argmax"` uses the largest predicted dummy-response score and is
the standard PLS-DA decoding rule. `classifier = "lda"` fits a regularized
linear discriminant analysis model in the PLS score space following Fisher's
discriminant principle (Fisher, 1936). `classifier = "cknn"` reranks
centroid-selected candidate classes by nearest-neighbour evidence in the
supervised latent space, following the nearest-neighbour principle of Cover and
Hart (1967).

### Backends

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

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

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

## 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`, `lda`, or `cknn`. These options are not
regression models and are not used for numeric responses. The examples in this
section use the `iris` data for compact classification examples and the bundled
`colon` and `breast` omics datasets for score plots.

```{r chunk-002}
library(fastPLS)

set.seed(100)
X <- as.matrix(iris[, c("Sepal.Width", "Petal.Length", "Petal.Width")])
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 = "argmax"`, which predicts 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()`
can optionally return ranked classes by setting `top` or `top5 = TRUE`. The
ordinary prediction in `Ypred` is always the rank-1 class. When `top > 1`,
`Ypred_top` contains the ordered candidate labels for each sample: `rank1` is
the predicted class, `rank2` is the next most likely class, and so on. If
available, `Ypred_top_score` contains the corresponding class scores used to
create that ranking.

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

### Evaluate Classification Predictions

Use `evaluate()` to summarize predicted class labels or class-score matrices.
For classification, the complete output includes global metrics, per-class
metrics, and the confusion matrix.

```{r chunk-005}
eval_cls <- evaluate(
  observed = Ytest_cls,
  predicted = fit_cls$Ypred[["ncomp=2"]]
)

eval_cls
```

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

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,
  top_k = c(1, 2)
)
```

### Classification Heads

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

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

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

```{r chunk-010}
fit_cls_lda_cpu <- pls(
  Xtrain,
  Ytrain_cls,
  Xtest,
  Ytest_cls,
  ncomp = 1:2,
  method = "plssvd",
  seed = 100,
  classifier = "lda"
)

head(fit_cls_lda_cpu$Ypred)
```

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

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

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

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

```{r colon-binary-pls-plots, fig.width = 5.5, fig.height = 4.5}
data("colon", package = "fastPLS")
set.seed(200)
X_two <- as.matrix(colon$X)
y_two <- colon$y
train_two <- unlist(tapply(seq_along(y_two), y_two, function(i) {
  sample(i, floor(0.7 * length(i)))
}))

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

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

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

```{r breast-tcga-three-class-opls, fig.width = 7.5, fig.height = 4}
data("breast", package = "fastPLS")
X_three <- as.matrix(breast$X_train)
y_three <- breast$y_train
X_three_test <- as.matrix(breast$X_test)
y_three_test <- breast$y_test

opls_three <- pls(
  X_three,
  y_three,
  X_three_test,
  y_three_test,
  ncomp = 1:2,
  method = "opls",
  fit = TRUE,
  proj = TRUE,
  north = 2L,
  seed = 201
)

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

## Regression Tasks

For regression, responses are supplied as numeric vectors or matrices. The
examples in this section use `mtcars` to demonstrate univariate regression and
OPLS. Regression does not use `classifier = "argmax"`, `classifier = "lda"`, or
`classifier = "cknn"`: 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 <- matrix(mtcars$mpg, ncol = 1)
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, , drop = FALSE]
Ytest_reg <- Y_reg[reg_test_id, , drop = FALSE]
```

### Fit And Predict A Regression Model

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

```{r chunk-017}
fit_reg <- pls(
  Xreg_train,
  Ytrain_reg,
  Xreg_test,
  Ytest_reg,
  ncomp = 1:3,
  svd.method = "irlba",
  fit = TRUE,
  return_variance = FALSE
)

fit_reg$Q2Y
```

### Float32 Input

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

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

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

```{r chunk-021}
eval_reg <- evaluate(
  observed = Ytest_reg,
  predicted = pred_mpg,
  ytrain = Ytrain_reg
)

eval_reg
```

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

```{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 the same compiled CPU/CUDA core when available.

The cKNN memory strategy can also be used during cross-validation. For example,
`pls.single.cv(..., classifier = "cknn", cknn_memory = "blocked")` evaluates
the same candidate-kNN classifier while reducing fold-level prediction memory.
For very large scalar-component cKNN validations, `cknn_memory = "streaming"`
can further reduce peak RAM by constructing the training score cache in blocks
inside each fold.

### Ordinary k-fold CV

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

```{r chunk-023}
cv_kfold <- pls.single.cv(
  Xdata = Xtrain,
  Ydata = Ytrain_cls,
  ncomp = 2,
  kfold = 5,
  svd.method = "rsvd",
  seed = 103
)

cv_kfold$metrics
```

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

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

### 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",
  svd.method = "rsvd",
  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, Q2 for univariate regression, and RMSD-style
prediction error for multivariate regression outputs. Predictive arguments can
also be supplied as vectors to tune a compact grid. Fold-control arguments such
as `kfold` remain single settings for the whole CV run.

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

cv_opt$best_ncomp
```

### Interpreting `R2Y`, `Q2Y`, and `RMSD`

For regression models, `R2Y`, `Q2Y`, and `RMSD` answer different questions and
should not usually be identical. `Q2Y` is calculated by default from held-out
cross-validation predictions and measures predictive performance on samples that
were not used to fit the corresponding fold model. `RMSD` is also calculated
from held-out predictions and is reported on the response scale, so lower values
are better. `R2Y` is a training-set explained-variance estimate from one
additional model fitted on the full dataset; set `fit = FALSE` to skip
this extra fit when only cross-validated performance is needed. For
classification, `Q2Y` is calculated from held-out dummy-coded PLS-DA response
scores, `accuracy` reports decoded-label accuracy, and `R2Y` is calculated from
the full-data PLS-DA fit on the dummy-coded response scores. `RMSD` is not used
for classification.

```{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 empirical fraction `mean(Q2Y_permuted > Q2Y_observed)`, without a
`+1` correction. The full permutation table is stored in `permutation` and can
be visualized with `plot.permutation()`, where the x-axis is the correlation
between the original and permuted response structure and the y-axis shows R2
and Q2.

```{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 after shuffling the rows of `Xdata`. The observed
statistic is the median outer-CV `Q2Y`, and `Q2Ysampled` stores the permuted
median values. For metrics where larger is better, the p-value is
`mean(Q2Y_permuted >= Q2Y_observed)`; for loss metrics where smaller is better,
it is `mean(loss_permuted <= loss_observed)`.

```{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,
  svd.method = "rsvd",
  perm.test = TRUE,
  seed = 109
)

data.frame(
  observed_median_Q2Y = median(dcv_perm$Q2Y, na.rm = TRUE),
  p_value = dcv_perm$p.value
)
```

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

```{r chunk-030}
cv_select <- pls.single.cv(
  Xdata = Xtrain,
  Ydata = Ytrain_cls,
  ncomp = 1:3,
  kfold = 5,
  svd.method = "rsvd",
  seed = 106
)

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

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

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

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

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

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: the inner loop calls
`pls.single.cv()` to select the best number of components and any vector-valued
predictive hyperparameters, while the outer loop refits the selected model on
the corresponding outer training set and predicts the outer test fold. Both
`kfold_inner` and `kfold_outer` can be ordinary fold counts or `"loocv"`, and
both respect `constrain`.

```{r chunk-032}
cv_double <- pls.double.cv(
  Xdata = Xtrain,
  Ydata = Ytrain_cls,
  constrain = patient_id,
  ncomp = 1:2,
  kfold_inner = 3,
  kfold_outer = 3,
  svd.method = "rsvd",
  seed = 104
)

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

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",
  svd.method = "rsvd",
  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, depending on task type. |
| `Ypred` | Final cross-validated prediction for each sample. |
| `conf` | Classification confusion matrix. |
| `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 held-out metric used for inner/outer model selection.
- `medianQ2Y`, `CI95Q2Y`, `medianR2Y`, `CI95R2Y`, `medianRMSD`, and
  `CI95RMSD`: summaries across repeated outer CV runs, returned only when
  `runn > 1`. They are omitted for the default `runn = 1` output.
- `bcomp`: the most frequently selected number of components across all outer
  folds and repeated runs.
- `backend` and `method`: the default backend and PLS method used by the call.
  If vector-valued methods or backends are tuned in the inner loop, the selected
  fold-level values are stored in `results[[run]]$best_parameters`.
- `selection_metric`: the criterion used by the inner CV loop. The default
  `"auto"` means accuracy for classification and RMSD-style prediction error for
  multivariate regression.

## PCA and SVD Utilities

`fastsvd()` provides direct access to the same SVD backends 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. `pca()` builds PCA scores and loadings from those backends, following
the classical principal-component construction of Pearson (1901). The fitted
PCA object stores the training centering, scaling, and loading matrix, so the
same projection can be applied to independent data either during fitting with
`xtest` or later with `predict()`.

```{r chunk-034}
s <- fastsvd(Xtrain, ncomp = 3, seed = 104)
names(s)
s$d
```

```{r chunk-035, fig.width = 5, fig.height = 4}
pc <- pca(Xtrain, ncomp = 3, xtest = Xtest, seed = 105)
plot(pc, groups = Ytrain_cls, ellipse = TRUE)
```

```{r chunk-036}
head(pc$scores_test)
head(predict(pc, Xtest, ncomp = 2))
```

The same workflow can be used with the bundled breast cancer example. Here PCA
is fitted only on the training matrix, the training scores are plotted first,
and the independent test samples are projected with the stored training
centering, scaling, and loadings before being overlaid on the same score map.

```{r chunk-037, fig.width = 5, fig.height = 4}
data(breast)

pp <- pca(breast$X_train)
plot(pp, bg = breast$y_train)

qq <- predict(pp, breast$X_test)
points(qq, bg = breast$y_test)
```

## CUDA and Metal Availability

CUDA and Apple Metal are optional. Use `has_cuda()` before selecting
`backend = "cuda"` and `has_metal()` before selecting `backend = "metal"`.

```{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()) {
  fit_metal <- pls(
    Xtrain,
    Ytrain_cls,
    Xtest,
    Ytest_cls,
    ncomp = 1:2,
    backend = "metal"
  )

  pca_metal <- pca(
    Xtrain,
    backend = "metal"
  )
}
```

## Helper Functions

`fastcor()` computes fast Pearson-style correlations. `ViP()` returns variable
importance in projection trajectories for fitted PLS models, 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. These functions are ordinary R-level helpers and are
independent of CUDA or Metal availability.

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

```{r chunk-042}
vip <- ViP(fit_reg)
dim(vip)
```

## Implementation of the Main Methods

### Singular-vector backends

The PLS algorithms in `fastPLS` are separated from the inner singular-vector
backend. The compiled CPU implementation supports two truncated SVD routes:
`irlba` and randomized SVD. The `irlba` backend follows the implicitly
restarted Lanczos bidiagonalization strategy of Baglama and Reichel (2005). It
estimates only the leading singular triplets needed by the requested PLS
components, rather than computing a full dense SVD. The randomized SVD backend
follows the randomized range-finder framework described by Halko, Martinsson,
and Tropp (2011): it draws a Gaussian test matrix, forms a low-dimensional
sketch of the target matrix, optionally applies power iterations,
orthonormalizes the sketch, and solves the final SVD in the small projected
space. In `pls()` and the cross-validation helpers, these controls are supplied
directly through `...`; in stand-alone SVD calls they are supplied directly to
`fastsvd()`.

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

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

For CUDA builds, GPU execution is selected with `backend = "cuda"`, so the main
matrix operations are performed by GPU-native PLS routines. When the
cross-covariance \(S = X^T Y\) would be large, selected backends can use the
matrix-free `xprod` route. In that route, the code avoids explicitly storing
`S` and evaluates the products required by the SVD through identities such as
\(Sv = X^T(Yv)\) and \(S^Tu = Y^T(Xu)\). For classification problems with many
classes, label-aware products can avoid materializing a dense one-hot response
matrix. The low-rank streamed prediction path is inspired by the FlashSVD
streaming low-rank inference design of Shao et al. (2025), but is adapted here
to the PLS prediction factors exposed by `fastPLS`.
In ordinary use, leave `xprod = NULL`; fastPLS then chooses the matrix-free
route automatically only when the cross-product is large enough to justify it.

### PLSSVD

`method = "plssvd"` implements the direct cross-covariance SVD formulation of
PLS, within the same broad PLS regression and PLS-DA framework described by
Wold, Sjostrom, and Eriksson (2001), Barker and Rayens (2003), and Boulesteix
and Strimmer (2007). After preprocessing, the algorithm computes the dominant
singular subspace of \(X^T Y\) and reuses that subspace to generate the
requested component path.
The maximum meaningful number of PLSSVD components is capped by
`min(nrow(X), ncol(X), ncol(Y))`; if a larger value is requested, the package
uses the valid internal rank and records the capped component count.

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

### SIMPLS

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

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

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

### OPLS

`method = "opls"` follows the orthogonal PLS idea introduced by Trygg and Wold
(2002). It first removes response-orthogonal variation from `X` and then fits
the optimized PLS core to the filtered matrix. The number of orthogonal
components is controlled by `north`. The filtering step estimates orthogonal
scores and loadings, subtracts the corresponding variation from `X`, and stores
the orthogonal filter so the same transformation is applied at prediction time.
Users select OPLS with `method = "opls"`; the lower-level implementation details
are handled internally.

### Kernel PLS

`method = "kernelpls"` follows the kernel PLS formulation of Rosipal and Trejo
(2001). It constructs a kernel representation before applying the selected
inner PLS model. The supported kernels are `"linear"`, `"rbf"`, and `"poly"`.
For nonlinear kernels, the training kernel matrix is centered and the PLS model
is fitted in kernel space. For the linear kernel, fastPLS avoids unnecessary
kernel materialization and uses the ordinary linear PLS path.

### Classification heads

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

The LDA head, `classifier = "lda"`, fits a
regularized linear discriminant model on the PLS latent scores, following the
linear discriminant principle of Fisher (1936). For each class, fastPLS
estimates a score-space mean, an empirical prior, and a common pooled
within-class covariance matrix. A small ridge term scaled by the average
covariance diagonal is added for numerical stability. Prediction uses the
standard LDA discriminant score and returns the largest-scoring class. The CUDA
and Metal routes are selected automatically from `backend` when available.

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

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


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

- Baglama, J. and Reichel, L. (2005). Augmented implicitly restarted Lanczos
  bidiagonalization methods. *SIAM Journal on Scientific Computing*, 27,
  19-42.
- 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.
- Cover, T. and Hart, P. (1967). Nearest neighbor pattern classification.
  *IEEE Transactions on Information Theory*, 13, 21-27.
- 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.
- Hotelling, H. (1931). The generalization of Student's ratio. *Annals of
  Mathematical Statistics*, 2, 360-378.
- Johnson, W. B. and Lindenstrauss, J. (1984). Extensions of Lipschitz mappings
  into a Hilbert space. *Contemporary Mathematics*, 26, 189-206.
- 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.
- Pearson, K. (1901). On lines and planes of closest fit to systems of points
  in space. *Philosophical Magazine*, 2, 559-572.
- 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.
- Shao, Z., Wang, Y., Wang, Q., Jiang, T., Du, Z., Ye, H., Zhuo, D., Chen, Y.,
  and Li, H. (2025). FlashSVD: memory-efficient inference with streaming for
  low-rank models. arXiv:2508.01506.
- 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()
```