Partial Least Squares Regression (PLS)

Few methods can claim to have shaped a whole discipline. Partial least squares regression (PLS) is one of them. If you have ever seen a near-infrared analyser predict protein in grain, moisture in a tablet or octane in fuel, there is a very good chance a PLS model was doing the work.

PLS keeps the strategy of principal component regression: compress the many correlated wavelengths into a few latent variables, then regress. It changes one thing. Instead of compressing the spectra blindly, it lets the property you want to predict guide the compression. That single change is the subject of this entry.

The story begins in econometrics, with the Swedish statistician Herman Wold. In the 1960s he developed iterative algorithms that estimate models by alternating simple least squares regressions, starting with principal components [1]. In the 1970s he extended them into a framework for models with latent variables, which he called soft modelling, deliberately avoiding strong distributional assumptions [2]. Wold first called these procedures NILES (nonlinear iterative least squares); by the mid-1970s they had become NIPALS (nonlinear iterative partial least squares), and that is where the name comes from.

The jump to chemistry came through his son, Svante Wold, a chemist at Umeå University and one of the founders of chemometrics, together with Harald Martens, who was working on near-infrared analysis of food and agricultural products in Norway. NIR spectroscopy produced exactly the data that classical regression could not handle: hundreds of strongly correlated wavelengths measured on a few dozen samples. In 1983 Svante Wold, Martens and Herman Wold described how PLS solves this multivariate calibration problem, in a chapter of the proceedings of a conference on matrix pencils [3]. A year later, Wold and co-workers analysed the method from the numerical analysis side and connected it with the conjugate gradient method [7].

What brought PLS to the bench chemist was the 1986 tutorial by Paul Geladi and Bruce Kowalski in Analytica Chimica Acta [4], one of the most cited papers in chemometrics. It walked through the algorithm step by step, with examples and geometric pictures, at the moment when personal computers were putting multivariate methods within reach of any laboratory with a spectrometer.

The name has always caused some confusion, because "partial least squares" does not say much about what the method does. Svante Wold later proposed reading PLS as projection to latent structures [5],[16], which describes it better: the data are projected onto a few latent directions. Both readings share the acronym.

You have a set of samples (tablets, grain, fuel) and for each one you measure a spectrum, the rows of $\mathbf{X}$, and a reference property, the elements of $\mathbf{y}$. Typically there are far more wavelengths than samples, and neighbouring wavelengths are almost perfectly correlated. As we saw in the previous entries:

• Multiple linear regression has no unique solution when $p\ge n$ (it can fit the calibration samples perfectly and still predict new ones badly), and is wildly unstable when $p<n$ but the predictors are collinear.
• Choosing a few wavelengths by hand throws away most of the spectrum and has to be redone for every new application.
• PCR handles both problems, but builds its components without looking at $\mathbf{y}$, so it may spend its first components on large variations that have nothing to do with the analyte.

PLS keeps the compression of PCR and fixes its blind spot.

The core of PLS fits in one sentence:

Find directions in X-space along which the samples vary a lot and which are strongly related to y.

Formally, for mean-centred $\mathbf{X}$ and $\mathbf{y}$, the first PLS direction is the unit vector $\mathbf{w}$ that maximizes the covariance between the scores $\mathbf{t}=\mathbf{Xw}$ and the response:

The decomposition on the right is the key to the whole method. The standard deviation of $\mathbf{y}$ does not depend on $\mathbf{w}$, so maximizing the covariance means finding a direction that is both well correlated with the response and carries a decent amount of spectral variation. Compare the three related criteria:

Maximizing the variance alone ignores the response. Maximizing the correlation alone is dangerous with spectra: with many variables you can almost always find some tiny, noisy direction that correlates perfectly with $\mathbf{y}$ by chance. Covariance is the compromise. And it has a remarkably simple solution: for centred data, the best direction is just proportional to ${\mathbf{X}}^T\mathbf{y}$, the vector of covariances between each wavelength and the response.

For a single response, the maximum-correlation direction is exactly the MLR regression vector (whenever MLR has a unique solution), so the table reads as a continuum: PCR looks only at variance, MLR only at correlation, and PLS sits in between with covariance. This is why MLR overfits wide data so easily, and why PLS, which also rewards spectral spread, is much less prone to it. Stone and Brooks turned this continuum into a single method with a tuning knob, continuum regression [18].

The demo below compares the first PCA direction with the first PLS direction on two wavelengths. Move the slider to change where the analyte information lies.

When the analyte information lies along the big spread of the data, the two directions agree. When it lies along the small spread, PC1 keeps following the big spread, while the PLS direction tilts towards the colour gradient and its scores predict the analyte much better, although a single PLS component still leaves most of $\mathbf{y}$ unexplained here (compare its $R^2$ with the MLR tile); the second component will pick up the rest. Notice that the PLS direction rarely points exactly along the relevant direction either: it is pulled partly towards the large spectral variation. That is the covariance compromise at work.

PLS describes both blocks of data with the same scores:

Here $\mathbf{T}$ ($n\times A$) contains the scores of the $A$ components, $\mathbf{P}$ ($p\times A$) the X-loadings, $\mathbf{q}$ ($A\times 1$) the y-loadings, and $\mathbf{E}$ and $\mathbf{f}$ the residuals. The structure looks just like PCR. The difference is entirely in how $\mathbf{T}$ is computed: with the help of $\mathbf{y}$.

The most widely taught algorithm is NIPALS. Here is the version for a single response (PLS1). Start with mean-centred data, ${\mathbf{X}}_0=\mathbf{X}$ and ${\mathbf{y}}_0=\mathbf{y}$, and for each component $a=1,2,\dots ,A$:

Every step is a small, well-conditioned least squares problem, and the solution is confined to an $A$-dimensional subspace, which is why the method works even when ${\mathbf{X}}^T\mathbf{X}$ (the matrix of the normal equations) is hopelessly singular. The deflation also makes the scores of different components mutually orthogonal, just like in PCR. The weight vectors ${\mathbf{w}}_a$ are orthonormal too; the loadings ${\mathbf{p}}_a$ in general are not.

After $A$ components, the model can be written as a single regression vector on the original wavelengths:

where $\mathbf{W}=[{\mathbf{w}}_1,\dots ,{\mathbf{w}}_A]$, $\mathbf{P}=[{\mathbf{p}}_1,\dots ,{\mathbf{p}}_A]$ and $\mathbf{q}=[q_1,\dots ,q_A{]}^T$. The matrix ${\mathbf{W}}^{\ast }=\mathbf{W}({\mathbf{P}}^T\mathbf{W}{)}^{-1}$ is needed because the weights of later components act on deflated data, not on the original spectra; this product translates them back, so that $\mathbf{T}=\mathbf{X}{\mathbf{W}}^{\ast }$. Software calls it ${\mathbf{W}}^{\ast }$ or $\mathbf{R}$ (SIMCA, PLS_Toolbox), or x_rotations_ (scikit-learn).

PLS1 predicts one response at a time: protein, moisture, fat, octane number. It is by far the most common case in spectroscopic calibration, and it is the algorithm above.

PLS2 predicts several responses at once, using a response matrix $\mathbf{Y}$ instead of a vector. The weights are chosen to maximize the covariance with the whole of $\mathbf{Y}$, so the components are shared by all the responses [9].

PLS2 is attractive when the responses are related and share an underlying structure, or when one compact model is preferred. Separate PLS1 models are usually better when the responses are unrelated, need different preprocessing, or need different numbers of components; for a single property, a dedicated PLS1 model typically predicts at least as well. In practice many chemometricians build one PLS1 model per property.

The number of components $A$ is the critical tuning parameter of PLS, exactly as in PCR. Too few and the model misses real structure; too many and it starts to fit noise in the calibration set and predicts new samples worse.

The standard tool is cross-validation, and it works exactly as described for PCR: leave samples out, predict them with a model that never saw them, and compute the root mean square error of cross-validation (RMSECV) for each $A$. The only difference is that the PLS model, weights included, is rebuilt inside every fold.

The demo below uses the same twenty simulated NIR spectra as in the PCR entry. Compare PLS with PCR (the grey curve), and then make the analyte signal smaller.

Look at the green bars and switch between the two methods. With the default settings PLS explains more than half of the variance of $\mathbf{y}$ with two components, while the first three PCR components explain about 2 %: they spend themselves on the interferents and the baseline and leave $\mathbf{y}$ for a later component. Now make the analyte signal tiny. For the first four components both error curves stay high, because the interferents dominate everything. Then look at components 5 to 10: the red PLS curve reaches its plateau at about six components, while the grey PCR curve needs nine or ten to get down to the same error. Both methods end at a similar minimum error; PLS simply gets there with fewer components. The green dot marks the suggested number of components (the fewest within 5 % of the minimum), and a grey ring marks the absolute minimum when it is somewhere else. Look also at the regression vector: as components are added, it develops the shape needed to read the analyte band while cancelling the overlapping interferents. Add too many components and it becomes noisy, a visual sign of overfitting.

A typical curve has three parts:

1. A rapid decrease over the first components, as the main predictive structure is captured.
2. A minimum or plateau, where additional components bring little.
3. A slow rise, when components start describing noise that happens to correlate with $\mathbf{y}$ in the calibration set.

A large gap between RMSEC and RMSECV or RMSEP is a classic sign of overfitting: the model has learned the calibration samples rather than the chemistry.

A PLS model is not a black box. Several plots show what it has learned and whether it is reliable.

Score plots. Plotting ${\mathbf{t}}_1$ against ${\mathbf{t}}_2$ shows the samples in the latent space, like a PCA score plot but with axes oriented towards prediction. Use it to spot outliers, clusters, drifts over time, and to check that calibration and validation samples cover the same region.

Weights and loadings. The weights ${\mathbf{w}}_a$ show which wavelengths are used to build each component; the loadings ${\mathbf{p}}_a$ show the spectral pattern that each component describes. Plotted against wavelength, they should highlight chemically meaningful regions. If they do not, the model may be relying on a spurious correlation.

Regression vector. The final ${\mathbf{b}}_{\text{PLS}}$, plotted against wavelength, is the overall recipe the model applies to a spectrum. Positive and negative lobes typically reflect the analyte bands and the corrections for overlapping interferents.

VIP scores. The variable importance in projection summarizes the importance of each wavelength over all components [5],[10],[16]:

where $S{S}_a$ is the part of the sum of squares of $\mathbf{y}$ explained by component $a$. The squared VIP values average to 1, so "VIP greater than 1" is a common rule of thumb for important variables. Treat it as a guide rather than a test: its behaviour depends on the correlation structure of the data [11].

Predicted versus reference. Finally, plot the predictions against the reference values for calibration and validation samples. A good model gives points scattered tightly around the 1:1 line; curvature or offsets point to nonlinearity or bias.

The examples simulate NIR spectra of mixtures of a weak analyte and two interferents, one of them overlapping the analyte band, with a baseline that shifts and tilts from sample to sample and some error in the reference values. Thirty samples are used for calibration and fifty for validation. The code chooses the number of components by cross-validation and reports RMSEC, RMSECV and RMSEP. You should see RMSEC < RMSECV < RMSEP, and an RMSECV curve that falls steeply, bottoms out at about five components (analyte, two interferents, offset and tilt) and then rises slowly.

import numpy as np
from sklearn.cross_decomposition import PLSRegression
from sklearn.model_selection import cross_val_predict, KFold
from sklearn.metrics import mean_squared_error
import matplotlib.pyplot as plt

# --- Simulate NIR spectra: a weak analyte band overlapping an interferent ---
rng = np.random.default_rng(42)
n_samples, n_wavelengths = 80, 200
wavelengths = np.linspace(1000, 2500, n_wavelengths)
pure_spectra = np.array([
    np.exp(-((wavelengths - 1750) ** 2) / 10000),   # analyte
    np.exp(-((wavelengths - 1900) ** 2) / 20000),   # interferent 1
    np.exp(-((wavelengths - 1700) ** 2) / 15000),   # interferent 2, overlaps the analyte
])
concentrations = rng.uniform(0, 10, (n_samples, 3))
y = concentrations[:, 0] + rng.normal(0, 0.2, n_samples)              # reference method error
X = (concentrations * [0.3, 1, 1]) @ pure_spectra                       # weak analyte signal
offset = rng.normal(0, 0.2, (n_samples, 1))                             # baseline shift
slope = rng.normal(0, 0.2, (n_samples, 1)) * (wavelengths - 1750) / 750  # baseline tilt
X = X + offset + slope + rng.normal(0, 0.02, (n_samples, n_wavelengths))

X_cal, X_val, y_cal, y_val = X[:30], X[30:], y[:30], y[30:]

# --- Choose the number of components by cross-validation ---
# scale=False: spectra are mean-centred but not scaled to unit variance
max_components = 12
kf = KFold(n_splits=10, shuffle=True, random_state=42)
rmsecv = []
for a in range(1, max_components + 1):
    y_cv = cross_val_predict(PLSRegression(n_components=a, scale=False), X_cal, y_cal, cv=kf)
    rmsecv.append(np.sqrt(mean_squared_error(y_cal, y_cv)))
rmsecv = np.array(rmsecv)
print("RMSECV:", np.round(rmsecv, 3))

# Parsimony: the fewest components within 2 % of the minimum RMSECV
n_opt = int(np.argmax(rmsecv <= 1.02 * rmsecv.min())) + 1
print(f"Chosen number of components: {n_opt}")

# --- Final model ---
pls = PLSRegression(n_components=n_opt, scale=False).fit(X_cal, y_cal)
rmsec = np.sqrt(mean_squared_error(y_cal, pls.predict(X_cal).ravel()))
rmsep = np.sqrt(mean_squared_error(y_val, pls.predict(X_val).ravel()))
print(f"RMSEC = {rmsec:.3f}, RMSECV = {rmsecv[n_opt - 1]:.3f}, RMSEP = {rmsep:.3f}")
# about 0.19, 0.24 and 0.26: the error of the reference values (0.2) sets the floor

# --- VIP scores ---
def vip_scores(model):
    t = model.x_scores_                 # n x A
    w = model.x_weights_                # p x A (unit-norm columns)
    q = model.y_loadings_               # 1 x A
    p = w.shape[0]
    ss = (q[0] ** 2) * np.sum(t ** 2, axis=0)              # sum of squares of y explained per component
    w_norm = w / np.linalg.norm(w, axis=0)
    return np.sqrt(p * (w_norm ** 2 @ ss) / ss.sum())

vip = vip_scores(pls)

fig, axes = plt.subplots(1, 3, figsize=(16, 4))
axes[0].plot(range(1, max_components + 1), rmsecv, "o-")
axes[0].axvline(n_opt, color="red", linestyle="--")
axes[0].set_xlabel("Number of components"); axes[0].set_ylabel("RMSECV")
axes[1].plot(wavelengths, pls.coef_.ravel())
axes[1].set_xlabel("Wavelength (nm)"); axes[1].set_ylabel("Regression coefficient")
axes[2].plot(wavelengths, vip)
axes[2].axhline(1, color="red", linestyle="--", alpha=0.6)
axes[2].set_xlabel("Wavelength (nm)"); axes[2].set_ylabel("VIP")
for ax in axes:
    ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

% PLS regression with plsregress (Statistics and Machine Learning Toolbox)
rng(42);
n_samples = 80;
n_wavelengths = 200;
wavelengths = linspace(1000, 2500, n_wavelengths);
pure_spectra = [
    exp(-((wavelengths - 1750).^2) / 10000);   % analyte
    exp(-((wavelengths - 1900).^2) / 20000);   % interferent 1
    exp(-((wavelengths - 1700).^2) / 15000)];  % interferent 2, overlaps the analyte

concentrations = 10 * rand(n_samples, 3);
y = concentrations(:, 1) + 0.2 * randn(n_samples, 1);            % reference method error
X = (concentrations .* [0.3 1 1]) * pure_spectra;                  % weak analyte signal
offset = 0.2 * randn(n_samples, 1);                                % baseline shift
slope = 0.2 * randn(n_samples, 1) .* (wavelengths - 1750) / 750;   % baseline tilt
X = X + offset + slope + 0.02 * randn(n_samples, n_wavelengths);

X_cal = X(1:30, :);    y_cal = y(1:30);
X_val = X(31:end, :);  y_val = y(31:end);

% --- Cross-validation for 1 to 12 components in one call ---
% plsregress mean-centres X and y. MSE(2, a+1) is the CV error with a components.
max_comp = 12;
[~, ~, ~, ~, ~, ~, MSE] = plsregress(X_cal, y_cal, max_comp, 'cv', 10);
rmsecv = sqrt(MSE(2, 2:end));
disp(round(rmsecv, 3));

% Parsimony: the fewest components within 2 % of the minimum RMSECV
n_opt = find(rmsecv <= 1.02 * min(rmsecv), 1);
fprintf('Chosen number of components: %d\n', n_opt);

% --- Final model (BETA includes the intercept as its first element) ---
[~, ~, ~, ~, BETA] = plsregress(X_cal, y_cal, n_opt);
y_cal_pred = [ones(size(X_cal, 1), 1), X_cal] * BETA;
y_val_pred = [ones(size(X_val, 1), 1), X_val] * BETA;
fprintf('RMSEC = %.3f, RMSECV = %.3f, RMSEP = %.3f\n', ...
    sqrt(mean((y_cal - y_cal_pred).^2)), rmsecv(n_opt), sqrt(mean((y_val - y_val_pred).^2)));

figure;
subplot(1, 3, 1); plot(1:max_comp, rmsecv, 'o-'); xline(n_opt, 'r--');
xlabel('Number of components'); ylabel('RMSECV'); grid on;
subplot(1, 3, 2); plot(wavelengths, BETA(2:end));
xlabel('Wavelength (nm)'); ylabel('Regression coefficient'); grid on;
subplot(1, 3, 3); scatter(y_val, y_val_pred, 'filled'); hold on;
plot([0 10], [0 10], 'k--'); xlabel('Reference'); ylabel('Predicted'); grid on;

library(pls)
set.seed(42)

n_samples <- 80
n_wavelengths <- 200
wavelengths <- seq(1000, 2500, length.out = n_wavelengths)
pure_spectra <- rbind(
  exp(-((wavelengths - 1750)^2) / 10000),   # analyte
  exp(-((wavelengths - 1900)^2) / 20000),   # interferent 1
  exp(-((wavelengths - 1700)^2) / 15000)    # interferent 2, overlaps the analyte
)
concentrations <- matrix(runif(n_samples * 3, 0, 10), ncol = 3)
y <- concentrations[, 1] + rnorm(n_samples, 0, 0.2)                 # reference method error
X <- concentrations %*% diag(c(0.3, 1, 1)) %*% pure_spectra           # weak analyte signal
offset <- rnorm(n_samples, 0, 0.2)                                    # baseline shift
slope <- outer(rnorm(n_samples, 0, 0.2), (wavelengths - 1750) / 750)  # baseline tilt
X <- X + offset + slope +
     matrix(rnorm(n_samples * n_wavelengths, 0, 0.02), nrow = n_samples)

dat <- data.frame(y = y)
dat$X <- I(X)
cal_idx <- 1:30
val_idx <- 31:80

# --- PLS with 10-fold cross-validation (plsr centres X; no scaling by default) ---
pls_model <- plsr(y ~ X, data = dat, subset = cal_idx,
                  ncomp = 12, validation = "CV", segments = 10)
rmsecv <- RMSEP(pls_model, estimate = "CV")$val[1, 1, -1]   # drop intercept-only model
print(round(rmsecv, 3))

# Parsimony: the fewest components within 2 % of the minimum RMSECV
n_opt <- unname(which(rmsecv <= 1.02 * min(rmsecv))[1])
cat("Chosen number of components:", n_opt, "\n")

y_cal_pred <- drop(predict(pls_model, ncomp = n_opt, newdata = dat[cal_idx, ]))
y_val_pred <- drop(predict(pls_model, ncomp = n_opt, newdata = dat[val_idx, ]))
cat(sprintf("RMSEC = %.3f, RMSECV = %.3f, RMSEP = %.3f\n",
            sqrt(mean((y[cal_idx] - y_cal_pred)^2)), rmsecv[n_opt],
            sqrt(mean((y[val_idx] - y_val_pred)^2))))

par(mfrow = c(1, 3))
plot(rmsecv, type = "b", xlab = "Number of components", ylab = "RMSECV")
abline(v = n_opt, lty = 2, col = "red")
plot(wavelengths, coef(pls_model, ncomp = n_opt), type = "l",
     xlab = "Wavelength (nm)", ylab = "Regression coefficient")
plot(y[val_idx], y_val_pred, pch = 16, xlab = "Reference", ylab = "Predicted")
abline(0, 1, lty = 2)

PLS is the default choice for:

• NIR, MIR and Raman calibration: predicting protein, moisture, fat, octane number and countless other properties from spectra.
• Process analytical technology: real-time monitoring of manufacturing processes with in-line spectroscopy.
• Any regression with many correlated variables, or with more variables than samples.

Consider alternatives when:

• you have few, nearly uncorrelated predictors: ordinary least squares is simpler and just as good;
• you need strict variable selection: PLS uses every variable, so consider the lasso or sparse PLS;
• the relationship is strongly nonlinear: standard PLS is a linear method;
• the goal is classification: use PLS-DA instead.

Strengths:

• Handles wide data: works even when $p\gg n$.
• Handles collinearity: correlated wavelengths are an asset rather than a problem, since the latent variables average them.
• Uses the response: usually reaches its best performance with fewer components than PCR.
• Fast: every step is a small least squares problem.
• Interpretable: scores, loadings, weights, VIP and the regression vector all have a meaning you can check against chemistry.

Limitations:

• Linear: nonlinear effects (strong scattering, detector saturation) are not captured.
• The number of components matters: too many overfit, too few underfit, so careful validation is essential.
• Sensitive to outliers: extreme samples can distort the latent space; check for them before modelling.
• Keeps every variable: irrelevant spectral regions add noise unless you remove them (see data reduction for iPLS and other variable selection methods).
• Needs chemical judgement: a model that predicts well for the wrong reasons will eventually fail. Always check that it makes chemical sense.

In head-to-head comparisons, PLS, PCR and ridge regression often reach similar prediction errors when each is tuned properly [14]. PLS owes its popularity as much to its interpretability and its efficient use of few components as to its accuracy.

• PLS-DA (discriminant analysis) codes class membership as dummy responses and adapts PLS to classification, for example to authenticate foods or detect counterfeit drugs [12].
• OPLS (orthogonal PLS) separates the variation in $\mathbf{X}$ that is related to $\mathbf{y}$ from the variation that is not, which makes the model easier to interpret [13].
• Multi-block PLS handles data from several sources (for example NIR, Raman and process variables) while keeping their block structure.
• Kernel and nonlinear PLS extend the method to curved relationships.
• Sparse PLS adds L1 penalties to the weights so that each component uses only some of the variables [17].

PLS compresses the spectra into a few latent variables, like PCR, but it chooses each latent direction to maximize the covariance with the response. Because covariance is correlation times spread, that means directions that both carry real spectral variation and relate to what you want to predict. The NIPALS algorithm builds them one at a time with small least squares regressions and deflation, and the resulting model can always be written as a single regression vector over the wavelengths. Seen from least squares, PLS1 is conjugate gradient on the normal equations stopped early: the number of components is how early you stop.

The number of components is the one parameter that really matters, and it must be chosen by validation, never by the calibration fit. Choose it well, check that the model makes chemical sense, and PLS will handle exactly the kind of data that broke ordinary least squares. For a classic introduction, see Geladi and Kowalski [4]; for the chemometric perspective, Wold, Sjöström and Eriksson [5] and Martens and Næs [8].