"""
The module :mod:`chemotools.regression._pls_regression` implements PLS regression
backed by the fast and numerically stable Improved Kernel PLS algorithms from the
`ikpls <https://github.com/sm00thix/ikpls>`_ package, with automatic explained
variance ratios for both X and Y spaces computed from fitted scores and loadings.
"""
# Authors: Pau Cabaneros, Ole-Christian Galbo Engstrøm
# License: MIT
from numbers import Integral
from typing import Tuple
import numpy as np
from ikpls.sklearn import PLS as _IkplsPLS
from sklearn.utils._param_validation import Interval, Options
from chemotools._doc_mixin import DocLinkMixin
[docs]
class PLSRegression(DocLinkMixin, _IkplsPLS):
"""PLS regression via Improved Kernel PLS with automatic explained variance.
This estimator wraps the fast, exact Improved Kernel PLS algorithms [1]_
from the `ikpls <https://github.com/sm00thix/ikpls>`_ package [2]_ in a
scikit-learn-conformant regressor + transformer, and automatically
calculates explained variance ratios for both X-space and Y-space after
fitting, making it easy to use with diagnostic plots and following the
same API as PCA.
Compared to ``sklearn.cross_decomposition.PLSRegression`` (NIPALS) [3]_ and [4]_,
the Improved Kernel PLS algorithms are faster [5]_ while still being numerically
stable [6]_. Predictions and regression coefficients agree with NIPALS;
individual weight/score/loading vectors are identical up to a per-component sign,
which is arbitrary in PLS.
Additional capabilities inherited from ikpls:
- :meth:`predict_all_components` returns predictions for every number of
components ``1..n_components`` in a single call (shape
``(n_components, n_samples, n_targets)``), making component selection
via (cross-)validation cheap.
- Fine-grained preprocessing on the backend: the ikpls model supports
independent ``center_X`` / ``center_Y`` / ``scale_X`` / ``scale_Y`` and a
configurable ``ddof``. This wrapper exposes only a single ``scale`` flag
(``X`` and ``Y`` are always mean-centered) and fixes ``ddof = 1``, to
mirror the previous scikit-learn-backed ``PLSRegression``. A maintainer can
enable the independent flags or a different ``ddof`` on the inner model,
but they change the fitted model and the explained-variance computation.
- Sample weights: the backend's ``fit`` accepts ``sample_weight`` (when
given, the inner ``X`` / ``Y`` means and standard deviations become their
weighted variants). This wrapper does not expose ``sample_weight``;
whether to enable it is left to the maintainer.
Parameters
----------
n_components : int, default=2
Number of components to keep. Should be in
[1, min(n_samples, n_features)].
scale : bool, default=True
Whether to scale ``X`` and ``Y`` to unit standard deviation before
fitting. ``X`` and ``Y`` are always mean-centered, matching the previous
scikit-learn-backed ``PLSRegression``. Internally this maps onto the
ikpls backend as ``center_X = center_Y = True`` and
``scale_X = scale_Y = scale``.
algorithm : int, default=1
Improved Kernel PLS algorithm to use, either 1 or 2. Algorithm 1 uses
``X`` directly, while algorithm 2 builds ``X.T @ X`` and is typically
faster for tall matrices (many more samples than features). Both
algorithms give the same results.
copy : bool, default=True
Whether to copy ``X`` and ``Y`` in fit before applying centering and
potentially scaling.
dtype : type, default=numpy.float64
Floating point dtype used for the computations.
Attributes
----------
explained_x_variance_ratio_ : ndarray of shape (n_components,)
Fraction of the total preprocessed training X sum of squares explained
by each component (not cumulative). See Notes for the calculation and
numerical caveats.
explained_y_variance_ratio_ : ndarray of shape (n_components,)
Fraction of the total preprocessed training Y sum of squares explained
by each component (not cumulative).
For multiple targets, sums of squares are aggregated across all targets.
These ratios describe training-data variance, not predictive performance
on unseen data. See Notes for the calculation and numerical caveats.
x_scores_ : ndarray of shape (n_samples, n_components)
The training X-scores ``T`` (the projection of the preprocessed
training ``X`` onto ``x_rotations_``), equal to ``transform(X)`` on
the training data. (There is no ``y_scores_`` attribute; the Y-scores
are available on demand via ``transform(X, y)`` -- see :meth:`transform`
for the caveat on how they relate to scikit-learn.)
All other fitted attributes (``x_weights_``, ``y_weights_``,
``x_loadings_``, ``y_loadings_``, ``x_rotations_``, ``y_rotations_``,
``coef_``, ``intercept_``, ``n_features_in_``, ``feature_names_in_``) are
inherited unchanged from the ikpls backend (``ikpls.sklearn.PLS``); see the
ikpls documentation (https://ikpls.readthedocs.io/en/latest/) for their
definitions.
References
----------
.. [1] Dayal, B. S., & MacGregor, J. F. (1997).
Improved PLS algorithms.
Journal of Chemometrics, 11(1), 73-85.
https://doi.org/10.1002/(SICI)1099-128X(199701)11:1%3C73::AID-CEM435%3E3.0.CO;2-%23
.. [2] Engstrøm, O.-C. G., Dreier, E. S., Jespersen, B. M., & Pedersen,
K. S. (2024).
IKPLS: Improved Kernel Partial Least Squares and Fast Cross-Validation
Algorithms for Python with CPU and GPU Implementations Using NumPy and
JAX.
Journal of Open Source Software, 9(99), 6533.
https://doi.org/10.21105/joss.06533
.. [3] sklearn.cross_decomposition.PLSRegression
https://scikit-learn.org/stable/modules/generated/sklearn.cross_decomposition.PLSRegression.html
.. [4] Wegelin, J. A. (2000).
A Survey of Partial Least Squares (PLS) Methods,
with Emphasis on the Two-Block Case.
Technical Report No. 371, Department of Statistics,
University of Washington, Seattle, WA
.. [5] Alin, A. (2009).
Comparison of PLS algorithms when number of objects is much larger
than number of variables.
Statistical Papers, 50, 711-720.
https://doi.org/10.1007/s00362-009-0251-7
.. [6] Andersson, M. (2009).
A comparison of nine PLS1 algorithms.
Journal of Chemometrics, 23(10), 518-529.
https://doi.org/10.1002/cem.1248
.. [7] Abdi, H. (2003).
Partial Least Squares (PLS) Regression.
In Lewis-Beck M., Bryman A., Futing T. (Eds.),
Encyclopedia of Social Sciences Research Methods.
Thousand Oaks (CA): Sage.
Examples
--------
**Basic usage with automatic variance calculation**
>>> from chemotools.regression import PLSRegression
>>> import numpy as np
>>>
>>> # Generate sample data
>>> X = np.random.randn(100, 50)
>>> y = X[:, 0] + 2*X[:, 1] + np.random.randn(100)*0.1
>>>
>>> # Fit model
>>> pls = PLSRegression(n_components=5)
>>> pls.fit(X, y)
>>>
>>> # Variance ratios are automatically available!
>>> print(
... f"LV1 explains {pls.explained_y_variance_ratio_[0]*100:.1f}%"
... )
>>> print(f"Total Y variance: {pls.explained_y_variance_ratio_.sum()*100:.1f}%")
>>>
>>> # Predictions for ALL component counts 1..5 in a single call
>>> all_predictions = pls.predict_all_components(X)
>>>
>>> # Use with plotting
>>> from chemotools.plotting import ExplainedVariancePlot
>>> plot = ExplainedVariancePlot(pls.explained_y_variance_ratio_)
>>> plot.show()
Notes
-----
**Variance Calculation:**
Component sums of squares are computed directly from squared X-score and
loading norms [7]_. Both X and Y contributions use X-scores, consistent
with regression-mode PLS.
Let ``t_a``, ``p_a``, and ``q_a`` be column ``a`` of ``x_scores_``,
``x_loadings_``, and ``y_loadings_``, respectively. Let ``X_c`` and ``Y_c``
denote the training blocks after applying the fitted means and, when
``scale=True``, the fitted standard deviations (``ddof=1``). Then::
explained_x_variance_ratio_[a] = sum(t_a**2) * sum(p_a**2) / sum(X_c**2)
explained_y_variance_ratio_[a] = sum(t_a**2) * sum(q_a**2) / sum(Y_c**2)
Each sum is over all elements of its vector or matrix. Thus the denominators
are the total sums of squares of the preprocessed blocks, not the residual
sums of squares after earlier components.
- For orthogonal fitted X-scores, these contributions equal the reduction
in residual sum of squares from sequential PLS deflation, up to numerical
precision [4]_ [7]_. No explicit deflation is performed here.
- X ratios sum to approximately 1 when the fitted components reconstruct
the preprocessed X block. Y ratios do not need to sum to 1: variation in
Y that the X-scores cannot explain remains in the residual.
- The same calculation is used for both algorithms. Ratios are non-negative
for finite fitted scores and loadings and positive total variance.
There is no explicit rank thresholding, clipping, or renormalization;
components near or beyond numerical rank depend on the backend's
numerical behavior.
- A block with zero total variance has undefined ratios (``NaN``).
**Differences from the previous scikit-learn (NIPALS) backend:**
- The NIPALS-specific ``max_iter`` and ``tol`` parameters do not exist:
Improved Kernel PLS is an exact, non-iterative solver.
- The single ``scale`` parameter is retained with the same meaning (``X``
and ``Y`` are always mean-centered, and ``scale`` toggles unit-variance
scaling of both). The ikpls backend's independent ``center_X`` /
``center_Y`` / ``scale_X`` / ``scale_Y`` are not exposed here.
- Score/weight/loading vectors may differ from NIPALS by a per-component
sign; predictions, coefficients, and explained variances are unaffected.
See Also
--------
sklearn.cross_decomposition.PLSRegression : scikit-learn's NIPALS-based PLS.
chemotools.plotting.ExplainedVariancePlot : Visualization for explained variance.
"""
_parameter_constraints: dict = {
"n_components": [Interval(Integral, 1, None, closed="left")],
"scale": ["boolean"],
"algorithm": [Options(Integral, {1, 2})],
"copy": ["boolean"],
"dtype": "no_validation",
}
def __init__(
self,
n_components: int = 2,
*,
scale: bool = True,
algorithm: int = 1,
copy: bool = True,
dtype: type = np.float64,
):
# This wrapper deliberately mirrors the previous scikit-learn-backed
# PLSRegression: it exposes a single ``scale`` flag (fixing
# ``center_X = center_Y = True`` and ``scale_X = scale_Y = scale``) and
# fixes ``ddof = 1`` (Bessel's correction, as scikit-learn uses). The
# ikpls backend supports independent center_X/center_Y/scale_X/scale_Y
# and a configurable ddof; a maintainer can enable them on the inner
# model, but they are not part of the public API here -- they affect the
# fitted model and the explained-variance computation (see the
# limitation below).
#
# ``ikpls.sklearn.PLS.__init__`` only stores the preprocessing/solver
# parameters as attributes; we set our five public ones directly (so
# ``get_params`` / ``clone`` / ``check_estimator`` see exactly the public
# parameters) and expose the fixed/derived backend flags as read-only
# properties for the inherited ``fit`` to consume.
self.n_components = n_components
self.scale = scale
self.algorithm = algorithm
self.copy = copy
self.dtype = dtype
# The inherited ikpls ``fit`` reads ``center_X`` / ``center_Y`` / ``scale_X``
# / ``scale_Y`` / ``ddof``; expose them as fixed (centering always on,
# ``ddof = 1``) / derived (scaling tied to ``scale``) read-only properties.
@property
def center_X(self) -> bool: # noqa: N802 - mirrors the ikpls attribute name
return True
@property
def center_Y(self) -> bool: # noqa: N802 - mirrors the ikpls attribute name
return True
@property
def scale_X(self) -> bool: # noqa: N802 - mirrors the ikpls attribute name
return self.scale
@property
def scale_Y(self) -> bool: # noqa: N802 - mirrors the ikpls attribute name
return self.scale
@property
def ddof(self) -> int:
return 1
[docs]
def fit(self, X: np.ndarray, y: np.ndarray) -> "PLSRegression": # type: ignore[ty:invalid-method-override] # sample_weight intentionally not exposed; ikpls PLS.fit's optional sample_weight is deliberately dropped
"""Fit model to data and compute explained variance ratios.
This method extends ``ikpls.sklearn.PLS.fit`` by storing the training
X-scores and calculating explained variance ratios. See the Notes
section of :class:`PLSRegression` for calculation details.
Parameters
----------
X : array-like of shape (n_samples, n_features)
Training vectors. Accepts numpy arrays, pandas DataFrames.
y : array-like of shape (n_samples,) or (n_samples, n_targets)
Target vectors. Accepts 1D (univariate) or 2D (multivariate) targets.
Returns
-------
self : PLSRegression
Fitted estimator with populated variance attributes:
``explained_x_variance_ratio_`` and ``explained_y_variance_ratio_``.
"""
self._validate_params()
# The inner PLS centers/scales X and y IN PLACE when copy=False, so
# snapshot the raw inputs first. Both the score computation (algorithm 2
# re-projects via transform) and the variance normalization re-apply
# the fitted preprocessing, so they must start from the original data, not
# data already preprocessed in place. With copy=True the inputs are left
# untouched, so no snapshot is needed.
if self.copy:
X_raw, y_raw = X, y
else:
X_raw = np.array(X, dtype=self.dtype)
y_raw = np.array(y, dtype=self.dtype)
super().fit(X, y)
# ikpls stores the training X-scores in inner_.T for algorithm 1 only
# (algorithm 2 leaves it None), so recompute them there. Use the inner
# model's own transform (not super().transform) so this internal recompute
# does not re-run scikit-learn's feature-name check: with copy=False X_raw
# is a label-stripped snapshot, which would otherwise emit a spurious "X
# does not have valid feature names" warning for a DataFrame input.
# np.asarray/np.copy so x_scores_ is a fresh array that never aliases the
# inner model's stored T.
self.x_scores_ = (
np.asarray(
self.inner_.transform(
X=np.asarray(X_raw, dtype=self.dtype),
n_components=self.n_components_,
)
)
if self.algorithm == 2
# inner_.T is the stored X-scores (non-None) in the algorithm-1 branch.
else np.copy(self.inner_.T)
)
# Calculate explained variance ratios automatically (from the raw inputs)
(
self.explained_x_variance_ratio_,
self.explained_y_variance_ratio_,
) = self._calculate_explained_variance(X_raw, y_raw)
return self
def _calculate_explained_variance(self, X, y) -> Tuple[np.ndarray, np.ndarray]:
"""Calculate explained variance ratios from fitted scores and loadings.
Requires fitted training scores, loadings, and preprocessing statistics.
See the Notes section of :class:`PLSRegression` for the formulas and
numerical caveats.
Parameters
----------
X : array-like of shape (n_samples, n_features)
Original training predictors, before centering and scaling.
y : array-like of shape (n_samples,) or (n_samples, n_targets)
Original training responses, before centering and scaling.
Returns
-------
x_ratios : ndarray of shape (n_components,)
Fraction of total preprocessed X variance explained by each
component.
y_ratios : ndarray of shape (n_components,)
Fraction of total preprocessed Y variance explained by each
component.
"""
# Convert to arrays and ensure y is 2D (handles pandas DataFrame/Series).
X = np.asarray(X, dtype=self.dtype)
y_array = np.asarray(y, dtype=self.dtype)
y = np.atleast_2d(y_array).T if y_array.ndim == 1 else y_array
# Apply the fitted preprocessing
inner = self.inner_
X_centered = X - inner.X_mean
y_centered = y - inner.Y_mean
if self.scale:
X_centered /= inner.X_std
y_centered /= inner.Y_std
# Component sums of squares from squared score and loading norms.
scores_ss = np.sum(self.x_scores_**2, axis=0)
explained_x_ss = scores_ss * np.sum(self.x_loadings_**2, axis=0)
explained_y_ss = scores_ss * np.sum(self.y_loadings_**2, axis=0)
# Zero total variance give undefined ratios (NaN)
with np.errstate(divide="ignore", invalid="ignore"):
x_ratios = explained_x_ss / np.sum(X_centered**2)
y_ratios = explained_y_ss / np.sum(y_centered**2)
return x_ratios, y_ratios
