Data Reduction

A near-infrared spectrum can contain two thousand numbers. A hyperspectral image can contain hundreds of thousands per pixel. It is tempting to feed all of them to a model, on the grounds that more data can only help. In practice, it often does not. Many of those numbers repeat what their neighbours already say, some record only noise, and some are dominated by effects that have nothing to do with the chemistry we care about. Data reduction is the craft of keeping the information and dropping the rest.

The idea that a signal can contain far fewer "real" pieces of information than raw numbers has deep roots. In 1948, Claude Shannon published "A Mathematical Theory of Communication", which defined information mathematically and showed that redundant parts of a message can be removed without losing anything [1]. Shannon was thinking about telephone lines, but the principle applies just as well to a spectrum whose neighbouring wavelengths carry almost the same information.

In chemometrics, variable selection became an active research topic in the 1990s, as instruments produced more and more variables. Wold, Johansson and Cocchi described the variable importance in projection (VIP) for PLS models in 1993 [2]. Centner and colleagues proposed uninformative variable elimination (UVE) in 1996, which compares each variable's regression coefficients with those of deliberately added random variables and removes the ones that behave no better than noise [3]. In 2000, Nørgaard and colleagues introduced interval PLS (iPLS), which splits the spectrum into windows and asks which windows predict best [4]. Together with PCA compression and plain common sense about which spectral regions are useful, these ideas form the modern toolkit.

Before building a multivariate model, it is worth asking: do I need all of these variables? In spectroscopy the answer is usually no. A typical NIR spectrum has around a thousand or more wavelength channels, but the chemistry behind it may involve only a few absorbing species, each with a handful of broad bands. There are five main reasons to reduce:

Computational cost. PCA and PLS are fast for ordinary data sets, but hyperspectral imaging, high-frequency process monitoring and multi-way data can produce matrices with hundreds of thousands of variables. Fewer variables make computation and storage manageable.

Noise. Variables that record mostly noise add noise to the model without adding signal. They make the calibration statistics worse and make it harder for the model to find the real patterns. Removing them complements noise reduction by smoothing.

Interpretability. A model that uses a few selected regions can be explained in chemical terms: "the model relies on the C-H band near 1730 nm and the O-H band near 1940 nm". This matters for validation, for regulatory acceptance and for scientific understanding.

Overfitting. When there are far more variables $p$ than samples $n$ (the "$p\gg n$" situation), a model has so much freedom that it can fit the calibration samples perfectly and still fail on new ones. Fewer variables mean fewer opportunities to fit noise.

The curse of dimensionality. A related, more geometric problem is captured by this expression, coined by Richard Bellman in the late 1950s in his work on dynamic programming [5]. As the number of dimensions grows, the volume of the space grows so quickly that any realistic number of samples becomes sparse, and distances between points become more and more alike, so that the "nearest" neighbour is hardly nearer than the farthest [6],[7]. Methods based on distances or densities (nearest neighbours, clustering) suffer most.

There are four broad families of reduction, from the simplest to the most data-driven: excluding regions you know are useless, merging neighbouring variables (binning), choosing variables with a statistical criterion (variable selection), and replacing the variables with a few new ones (PCA compression).

The simplest and often the most effective reduction is to leave out spectral regions that are known to be uninformative. It requires domain knowledge, but when the excluded regions are genuinely useless it is both safe and easy to explain.

In NIR spectroscopy of moist samples, two regions often deserve attention [8]:

• Around 1450 nm, the first overtone of the O-H stretching vibration of water.
• Around 1940 nm, the combination band of O-H stretching and bending.

In samples with a lot of water, these bands can be so intense that the measured signal becomes noisy, nonlinear or even saturated, and they may overwhelm the weaker bands of the analyte. In other applications, of course, water is the analyte, and these are the most important regions of all.

Other common candidates are the edges of the detector range, where sensitivity drops and noise rises, regions affected by atmospheric water vapour or carbon dioxide in mid-infrared spectra, and regions where nothing varies between samples.

How to decide what to exclude:

• Physical knowledge. Known interferences (water, CO₂, a cuvette or fibre absorption) can be excluded before any modelling.
• Visual inspection. Overlay all spectra. Regions with very high noise, flat-topped (saturated) peaks, or no variation between samples are candidates.
• Model diagnostics. Large loadings or regression coefficients in regions that make no chemical sense suggest the model is fitting noise or artefacts there.

Binning is the simplest data-driven reduction: group $b$ neighbouring points and replace each group by its average. A spectrum of 2000 points binned with $b=4$ becomes a spectrum of 500 points. Formally, for bin $k$,

where $p$ is the original number of variables and $\lfloor \cdot \rfloor$ means rounding down.

Binning does two things at once. It reduces the number of variables by a factor $b$, and, because each value is an average of $b$ points, it reduces independent random noise by a factor of about $\sqrt{b}$. The price is resolution: features narrower than a few bins get flattened. Try it:

Binning is appropriate when:

• The instrument samples the spectrum much more finely than the features require (the data are oversampled).
• The bands of interest are broad compared with the sampling interval, as in most NIR spectra.
• You want a quick, transparent way to reduce data size.

It is not appropriate when:

• The spectrum contains sharp features (Raman lines, gas-phase or high-resolution IR bands) whose height, width or position matters.
• The sampling interval is already matched to the width of the features.

Variable selection keeps a subset of the original variables (wavelengths, mass channels, chromatographic time points) and discards the rest. Unlike PCA, the variables that remain are still real, physically meaningful measurements [9],[10].

The classical strategies come from statistics:

Forward selection starts with no variables and adds them one at a time. At each step it adds the variable that most improves the model, for example the one that gives the largest drop in cross-validated prediction error. It stops when no remaining variable helps.

Backward elimination starts with all variables and removes them one at a time, each time dropping the variable whose removal hurts least (or helps most).

Both are greedy: they make the best choice at each step without looking ahead. Forward selection can miss pairs of variables that are useless alone but powerful together, and both become slow and unstable with thousands of correlated wavelengths. For spectra, methods that exploit the continuous nature of the signal usually work better.

Genetic algorithms search the space of variable subsets in a less greedy way. They start from a population of random subsets, keep the ones that give the lowest cross-validated error, and combine and randomly modify them over many generations, in imitation of natural selection. Riccardo Leardi and colleagues adapted them to spectroscopic calibration and showed how to use them without overfitting [11].

Penalised regression offers yet another route: the lasso sets most coefficients exactly to zero and so selects variables while it fits.

Interval PLS [4] builds on a simple observation: in a spectrum, information comes in regions, not in isolated wavelengths.

The bar plot of RMSECV per interval is a map of where the information lives in the spectrum, which makes iPLS very easy to explain. An extension called synergy iPLS tests combinations of two, three or four intervals together, to find complementary regions that are weak alone but strong together [12]. Moving-window PLS slides a window of fixed width along the spectrum one point at a time instead of using fixed, non-overlapping intervals, so the result does not depend on where the interval borders happen to fall [13].

VIP scores summarise how much each variable contributes to a PLS model [2],[14]. They combine two pieces of information: how much each PLS component explains of the response, and how much each variable contributes to that component (see the PLS article for the weights and scores used below):

Here $p$ is the number of variables, $A$ the number of PLS components, $w_{ja}$ the weight of variable $j$ in component $a$, ${\mathbf{t}}_a$ the scores of component $a$, and $q_a$ the coefficient that links those scores to the response. $S{S}_a$ is the part of the response's sum of squares explained by component $a$.

The formula is built so that the average of the squared VIP scores is exactly 1. This gives the usual reading:

• ${\text{VIP}}_j>1$: above-average contribution, usually considered important.
• ${\text{VIP}}_j<0.8$: usually considered unimportant and a candidate for removal (some authors use 0.5).
• In between: a grey zone that depends on the application.

Plotted against wavelength, VIP scores look like a spectrum whose peaks show where the model gets its information, which can often be matched to known absorption bands.

The panel below runs both methods on simulated spectra with the same recipe as the code example further down: three constituents, of which the second is correlated with the one we predict.

CARS [15] repeats PLS many times, each time on a different random subset of the samples (a Monte Carlo scheme). After each run, variables with small regression coefficients (in absolute value) are less likely to be kept, and the number of variables is reduced step by step following an exponentially decreasing function (EDF): many variables are dropped in the first runs and fewer and fewer later. In the spirit of natural selection, variables "compete" to survive. At the end, the subset with the lowest cross-validation error is chosen. CARS usually ends with far fewer variables than VIP and is popular in NIR applications, but its results can vary from run to run because of the random sampling. Fix the random seed so that a result can be reproduced, and repeat the selection with different seeds to look for variables that are selected consistently.

Variable selection keeps some variables. PCA compression replaces all of them with a few new ones. The data matrix $\mathbf{X}$ ($n\times p$) is written as

where $\mathbf{T}$ is the $n\times A$ score matrix, $\mathbf{P}$ the $p\times A$ loading matrix, and $\mathbf{E}$ the residuals. Using the scores $\mathbf{T}$ as input for later steps (regression, classification, clustering) reduces the number of variables from $p$ to $A$. The PCA article explains how the decomposition is computed.

How good is the compression? One natural measure is the reconstruction error, how well the original data can be rebuilt from $A$ components:

When the error falls to the level of the instrumental noise, the remaining components contain mostly noise, and keeping them stores noise rather than information. The panel below shows this for a small set of simulated spectra:

Three common approaches:

• Explained variance. Plot the cumulative explained variance and stop where it levels off, or at a threshold such as 99%. The threshold is arbitrary: in the panel above, 99% keeps one component too many (five, at 99.02%) and 95% one too few (three, at 95.2%), while the data contain four sources of variation.
• Scree plot. Plot the eigenvalues and look for the "elbow" where they change from steep to flat.
• Cross-validation. If the scores will feed a regression (principal component regression), choose $A$ to minimise the cross-validated prediction error. This optimises what you actually care about.

The examples simulate NIR-like spectra of three constituents, run an interval PLS scan, and compute VIP scores. The three languages draw different random numbers and split the samples into different cross-validation folds, so their numbers differ slightly, but the picture is the same.

import numpy as np
from sklearn.cross_decomposition import PLSRegression
from sklearn.model_selection import cross_val_predict
import matplotlib.pyplot as plt


def best_rmsecv(X, y, max_components=8, cv=10):
    """Lowest RMSECV over 1..max_components PLS components."""
    best = np.inf
    for nc in range(1, min(max_components, X.shape[1], X.shape[0] - 1) + 1):
        y_pred = cross_val_predict(PLSRegression(n_components=nc, scale=False), X, y, cv=cv)
        best = min(best, np.sqrt(np.mean((y - y_pred.ravel()) ** 2)))
    return best


def interval_pls(X, y, n_intervals=20, max_components=8):
    """Interval PLS: RMSECV for every interval and for the full spectrum."""
    edges = np.linspace(0, X.shape[1], n_intervals + 1).astype(int)
    rmsecv = [best_rmsecv(X[:, a:b], y, max_components) for a, b in zip(edges[:-1], edges[1:])]
    return np.array(rmsecv), best_rmsecv(X, y, max_components), edges


def compute_vip(X, y, n_components):
    """VIP scores of a PLS model (sklearn's x_weights_ have unit length)."""
    pls = PLSRegression(n_components=n_components, scale=False).fit(X, y)
    W = pls.x_weights_                     # (p, A)
    T = pls.x_scores_                      # (n, A)
    q = pls.y_loadings_[0]                 # (A,)
    ss = q**2 * np.sum(T**2, axis=0)       # y variance explained per component
    p = W.shape[0]
    return np.sqrt(p * (W**2 @ ss) / ss.sum())


# Simulated NIR data: three constituents, we predict the first one
rng = np.random.default_rng(42)
n_samples, n_wl = 80, 500
wavelengths = np.linspace(1000, 2500, n_wl)
comp1 = np.exp(-((wavelengths - 1200) ** 2) / 5000)
comp2 = np.exp(-((wavelengths - 1720) ** 2) / 3000)
comp3 = np.exp(-((wavelengths - 2100) ** 2) / 4000)

conc = rng.random((n_samples, 3))
conc[:, 1] = 0.5 * conc[:, 0] + 0.5 * rng.random(n_samples)   # correlated constituents
X = conc @ np.vstack([comp1, comp2, comp3]) + rng.normal(0, 0.02, (n_samples, n_wl))
y = conc[:, 0]

rmsecv_int, rmsecv_full, edges = interval_pls(X, y, n_intervals=20)
vip = compute_vip(X, y, n_components=4)

best = np.argmin(rmsecv_int)
print(f"Full-spectrum RMSECV: {rmsecv_full:.4f}")
print(f"Best interval: {wavelengths[edges[best]]:.0f}-{wavelengths[edges[best + 1] - 1]:.0f} nm, "
      f"RMSECV {rmsecv_int[best]:.4f}")
print(f"Wavelengths with VIP > 1: {np.sum(vip > 1)} of {n_wl}")

fig, axes = plt.subplots(2, 1, figsize=(11, 7))
centers = wavelengths[((edges[:-1] + edges[1:]) // 2).clip(max=n_wl - 1)]
width = wavelengths[1] - wavelengths[0]
colors = ['#2563eb' if k == best else '#94a3b8' for k in range(len(rmsecv_int))]
axes[0].bar(centers, rmsecv_int, width=np.diff(edges) * width, color=colors, edgecolor='white')
axes[0].axhline(rmsecv_full, color='#dc2626', linestyle='--', label=f'Full spectrum ({rmsecv_full:.4f})')
axes[0].set_xlabel('Wavelength (nm)'); axes[0].set_ylabel('RMSECV')
axes[0].set_title('iPLS: RMSECV per interval'); axes[0].legend()

axes[1].plot(wavelengths, vip, 'k-', linewidth=1)
axes[1].fill_between(wavelengths, vip, alpha=0.3, color='#2563eb')
axes[1].axhline(1.0, color='#dc2626', linestyle='--', label='VIP = 1')
axes[1].set_xlabel('Wavelength (nm)'); axes[1].set_ylabel('VIP score')
axes[1].set_title('Variable importance in projection'); axes[1].legend()

plt.tight_layout()
plt.show()

% Simulated NIR data: three constituents, we predict the first one
rng(42);
n_samples = 80;
n_wl = 500;
wavelengths = linspace(1000, 2500, n_wl);
comp1 = exp(-((wavelengths - 1200).^2) / 5000);
comp2 = exp(-((wavelengths - 1720).^2) / 3000);
comp3 = exp(-((wavelengths - 2100).^2) / 4000);

conc = rand(n_samples, 3);
conc(:, 2) = 0.5 * conc(:, 1) + 0.5 * rand(n_samples, 1);
X = conc * [comp1; comp2; comp3] + randn(n_samples, n_wl) * 0.02;
y = conc(:, 1);

% --- Interval PLS: 10-fold RMSECV for each of 20 intervals ---
n_int = 20;
edges = round(linspace(0, n_wl, n_int + 1));
rmsecv = zeros(n_int, 1);
for k = 1:n_int
    idx = edges(k) + 1 : edges(k + 1);
    [~, ~, ~, ~, ~, ~, mse] = plsregress(X(:, idx), y, 8, 'CV', 10);
    rmsecv(k) = sqrt(min(mse(2, 2:end)));     % row 2 = response, column 1 = 0 components
end
[~, ~, ~, ~, ~, ~, mse] = plsregress(X, y, 8, 'CV', 10);
rmsecv_full = sqrt(min(mse(2, 2:end)));

% --- VIP scores (following the plsregress documentation) ---
nc = 4;
[XL, YL, XS, ~, ~, ~, ~, stats] = plsregress(X, y, nc);
W0 = stats.W ./ sqrt(sum(stats.W.^2, 1));     % unit-length weights
p = size(XL, 1);
sumSq = sum(XS.^2, 1) .* sum(YL.^2, 1);       % y variance explained per component
vip = sqrt(p * sum(sumSq .* (W0.^2), 2) ./ sum(sumSq, 2));

[~, best] = min(rmsecv);
fprintf('Full-spectrum RMSECV: %.4f\n', rmsecv_full);
fprintf('Best interval: %.0f-%.0f nm, RMSECV %.4f\n', ...
        wavelengths(edges(best) + 1), wavelengths(edges(best + 1)), rmsecv(best));
fprintf('Wavelengths with VIP > 1: %d of %d\n', sum(vip > 1), n_wl);

figure;
subplot(2, 1, 1);
centers = wavelengths(round((edges(1:end-1) + edges(2:end)) / 2));
hb = bar(centers, rmsecv, 'FaceColor', 'flat'); hold on;
hb.CData = repmat([0.58 0.64 0.72], n_int, 1);
hb.CData(best, :) = [0.15 0.39 0.92];            % best interval in blue
yline(rmsecv_full, '--r', 'Full spectrum');
xlabel('Wavelength (nm)'); ylabel('RMSECV'); title('iPLS: RMSECV per interval');

subplot(2, 1, 2);
plot(wavelengths, vip, 'k-'); hold on;
yline(1, '--r');
xlabel('Wavelength (nm)'); ylabel('VIP score'); title('Variable importance in projection');

library(pls)

# Simulated NIR data: three constituents, we predict the first one
set.seed(42)
n_samples <- 80
n_wl <- 500
wavelengths <- seq(1000, 2500, length.out = n_wl)
comp1 <- exp(-((wavelengths - 1200)^2) / 5000)
comp2 <- exp(-((wavelengths - 1720)^2) / 3000)
comp3 <- exp(-((wavelengths - 2100)^2) / 4000)

conc <- matrix(runif(n_samples * 3), ncol = 3)
conc[, 2] <- 0.5 * conc[, 1] + 0.5 * runif(n_samples)
X <- conc %*% rbind(comp1, comp2, comp3) +
     matrix(rnorm(n_samples * n_wl, sd = 0.02), nrow = n_samples)
y <- conc[, 1]

# Lowest 10-fold RMSECV over 1..max_comp components
best_rmsecv <- function(Xi, y, max_comp = 8) {
  df <- data.frame(y = y)
  df$Xi <- I(Xi)
  nc <- min(max_comp, ncol(Xi), nrow(Xi) - 2)
  fit <- plsr(y ~ Xi, data = df, ncomp = nc, validation = "CV", segments = 10)
  min(RMSEP(fit, estimate = "CV", intercept = FALSE)$val)
}

# --- Interval PLS ---
n_int <- 20
edges <- round(seq(0, n_wl, length.out = n_int + 1))
rmsecv <- sapply(seq_len(n_int), function(k) best_rmsecv(X[, (edges[k] + 1):edges[k + 1]], y))
rmsecv_full <- best_rmsecv(X, y)

# --- VIP scores ---
compute_vip <- function(X, y, n_components) {
  df <- data.frame(y = y)
  df$X <- I(X)
  fit <- plsr(y ~ X, data = df, ncomp = n_components)
  W <- fit$loading.weights                 # unit-length weights
  ss <- fit$Yloadings[1, ]^2 * colSums(fit$scores^2)
  sqrt(ncol(X) * colSums(t(W^2) * ss) / sum(ss))
}
vip <- compute_vip(X, y, n_components = 4)

best <- which.min(rmsecv)
cat(sprintf("Full-spectrum RMSECV: %.4f\n", rmsecv_full))
cat(sprintf("Best interval: %.0f-%.0f nm, RMSECV %.4f\n",
            wavelengths[edges[best] + 1], wavelengths[edges[best + 1]], rmsecv[best]))
cat(sprintf("Wavelengths with VIP > 1: %d of %d\n", sum(vip > 1), n_wl))

par(mfrow = c(2, 1), mar = c(4.5, 4.5, 2.5, 1))
centers <- wavelengths[round((edges[-1] + edges[-length(edges)]) / 2)]
barplot(rmsecv, names.arg = round(centers), col = ifelse(seq_len(n_int) == best, "#2563eb", "#94a3b8"),
        xlab = "Interval centre (nm)", ylab = "RMSECV", main = "iPLS: RMSECV per interval")
abline(h = rmsecv_full, col = "#dc2626", lty = 2)

plot(wavelengths, vip, type = "l", xlab = "Wavelength (nm)", ylab = "VIP score",
     main = "Variable importance in projection")
abline(h = 1, col = "#dc2626", lty = 2)

Run the example and compare the two panels. iPLS finds that the interval between 1150 and 1222 nm (in blue) alone predicts about as well as the full spectrum (in Python, RMSECV 0.0046 for both), while the intervals with no analyte signal are dozens of times worse. VIP also flags the analyte band, but it flags a second region near 1720 nm as well, the band of the second constituent. That constituent is only correlated with the analyte, and iPLS shows that its region predicts poorly on its own. This is the warning above in action: VIP measures how much a variable is used by the model, not whether it is chemically causal.

Reduction usually helps with:

• Many variables and few samples, such as 2000 wavelengths and 50 samples, where overfitting is a real risk.
• Noisy or saturated regions, whose removal makes the model more robust.
• A need for interpretation, when the model must be explained in terms of chemical bands.
• Tight computational budgets, as in hyperspectral imaging or real-time process monitoring.

It calls for caution when:

• Information is subtle. Aggressive selection can remove variables that carry small but real chemical information. Always validate on independent data.
• Selection is tuned on the same samples used to judge the model. The selected variables then overfit; use nested cross-validation or a separate test set.
• The model is PLS on a well-behaved data set. PLS already compresses; variable selection may add complexity without benefit.
• Robustness matters. A model built on a few wavelengths can be more sensitive to instrument drift, temperature or wavelength shifts than a full-spectrum model, because it relies on fewer points [10].

Data reduction is about keeping information and dropping redundancy, noise and irrelevant effects. Start with the cheapest and most transparent steps: exclude regions you know are useless and bin oversampled spectra. Move on to statistical variable selection (iPLS, VIP, CARS) when you need interpretable or leaner models, and to PCA compression when you need a few uncorrelated inputs for later steps. Whatever you do, judge the result on data that played no part in the reduction.