Whittaker Smoothing

The ideas behind this smoother go back to Edmund Taylor Whittaker (1873 to 1956), a British mathematician and astronomer who served as Royal Astronomer of Ireland before spending most of his career at the University of Edinburgh. In 1923 he published On a New Method of Graduation in the Proceedings of the Edinburgh Mathematical Society [1]. "Graduation" was the actuaries' word for smoothing mortality tables, which are full of random fluctuations when the number of deaths in each age group is small. Whittaker's proposal was to treat smoothing as a compromise: find the curve that stays close to the observed values while keeping the sum of its squared differences, a measure of roughness, small.

The formulation was elegant, but in 1923 solving the resulting system of equations by hand was laborious. The method lived on mainly among actuaries, as Whittaker-Henderson graduation, after the American actuary Robert Henderson showed in 1924 how to carry out the computation more practically [6]. The same idea was later rediscovered by economists: the Hodrick-Prescott filter, widely used to separate trends from business cycles, is exactly a Whittaker smoother with second differences. It circulated as a working paper from 1980 and was published only in 1997, which is why the reference carries that date [5].

Eight decades after the original paper, Paul H. C. Eilers brought the method to chemometrics. His 2003 paper A Perfect Smoother in Analytical Chemistry [2] showed that, with sparse matrix techniques, the smoother takes only a few lines of code, runs in a time proportional to the number of points, handles missing data and uneven spacing, and is controlled by one intuitive parameter $\lambda$. Eilers had already worked with difference penalties when he and Brian Marx introduced P-splines (penalised B-splines) in 1996 [3]. Today Whittaker smoothing and its weighted variants underpin popular baseline correction methods such as asymmetric least squares, as well as general denoising.

The moving average, Gaussian and Savitzky-Golay filters all work locally: each smoothed value is computed from a small window of neighbours. Whittaker smoothing takes a different view. It asks a single question about the whole spectrum at once:

Which curve stays close to the data while being as smooth as possible?

Those two wishes pull in opposite directions, and the method is nothing more than a precise way of balancing them.

1. Fidelity: the smoothed curve $\mathbf{z}$ should stay close to the measured data $\mathbf{y}$.
2. Smoothness: $\mathbf{z}$ should not wiggle.

Perfect fidelity means $\mathbf{z}=\mathbf{y}$, keeping all the noise. Perfect smoothness means a featureless line, losing every band. Whittaker smoothing minimises a weighted sum of both:

where:

• $\mathbf{y}$ is the measured, noisy data ($n$ points),
• $\mathbf{z}$ is the smoothed curve we are looking for,
• $\lambda$ (lambda) is the smoothing parameter,
• ${\Delta }^d$ is the difference operator of order $d$, usually 2.

If this reminds you of least squares, it should. The first term is exactly the least squares misfit. The second is a penalty that punishes roughness. This is the same structure as ridge regression and many other regularised methods: fit the data, but pay a price for complexity.

is the sum of squared differences between data and curve. It is small when the curve follows the data closely.

The difference operator measures how much consecutive values change:

• First difference ($d=1$): $\Delta z_i=z_{i+1}-z_i$, a discrete slope.
• Second difference ($d=2$): ${\Delta }^2{z}_i=z_{i+2}-2z_{i+1}+z_i$, a discrete curvature.
• Third difference ($d=3$): ${\Delta }^3{z}_i=z_{i+3}-3z_{i+2}+3z_{i+1}-z_i$.

Squaring and summing these differences gives one number that is large for a jagged curve and small for a smooth one.

• $\lambda \to 0$: roughness is free, so $\mathbf{z}\to \mathbf{y}$ and nothing is smoothed.
• $\lambda \to \infty$: roughness is forbidden, so $\mathbf{z}$ becomes the least squares polynomial of degree $d-1$ (a straight line for $d=2$).
• In between: noise is removed while bands are kept.

Useful values cover many orders of magnitude, typically from about 1 to $10^8$ depending on the noise, the number of points and $d$. That is why $\lambda$ is almost always chosen on a logarithmic scale.

Despite looking like a complicated optimisation, the problem has a closed-form solution. Write the differences as a matrix product $\mathbf{Dz}$, where $\mathbf{D}$ is the difference matrix (shown below). Then

Just as in the derivation of the normal equations, we set the gradient with respect to $\mathbf{z}$ to zero:

which rearranges to

where $\mathbf{I}$ is the identity matrix. Smoothing a spectrum is therefore solving one linear system:

For second differences and $n=6$ points:

D = [ 1  -2   1   0   0   0 ]
    [ 0   1  -2   1   0   0 ]
    [ 0   0   1  -2   1   0 ]
    [ 0   0   0   1  -2   1 ]

Each row computes one second difference, $z_i-2z_{i+1}+z_{i+2}$. The matrix has $n-d$ rows. For first differences:

D = [ -1   1   0   0   0   0 ]
    [  0  -1   1   0   0   0 ]
    [  0   0  -1   1   0   0 ]
    [  0   0   0  -1   1   0 ]
    [  0   0   0   0  -1   1 ]

These matrices are almost entirely zeros, which is exactly why the computation is cheap. In code they are usually built by differencing the identity matrix $d$ times.

The panel below smooths a noisy spectrum with adjustable $\lambda$ and $d$. The metrics show the two terms of the objective. The lower plot is the trade-off curve: every dot is one value of $\lambda$, from almost no smoothing (top left: tiny misfit, large roughness) to heavy smoothing (bottom right: large misfit, almost no roughness). The highlighted dot is your current choice.

Try this:

1. Start at $\lambda \approx 5$ and observe moderate smoothing: the noise is much reduced and the narrow band keeps most of its height. At $\lambda \approx 30$ the noise is gone, but the narrow band has already lost about a sixth of its height.
2. Increase ${\log }_{10}\lambda$ to 6 and watch the curve stiffen until even the broad band is flattened.
3. Go down to $\lambda =0.1$ and see the curve follow the noise.
4. Switch between $d=1,2,3$ at a high $\lambda$ and see the curve approach a constant, a line and a parabola. The parabola needs the largest values: with $d=3$, go up to $\lambda =10^{11}$.
5. Watch the error against the true signal: it first falls, reaches a minimum, and then rises again as real features are smoothed away.

The trade-off curve also gives a useful landmark. If the curve $\mathbf{z}$ were exactly the true signal, the residuals would be pure noise, and the misfit would be about $n{\sigma }^2$, the number of points times the noise variance. In this panel $n=200$ and $\sigma =0.06$, so $n{\sigma }^2\approx 0.72$ (about $-0.14$ on the log scale). A misfit much smaller than that means the curve is still following noise; a misfit much larger means real signal is being removed.

The best curve, however, sits a little to the left of that landmark, because the smoother also absorbs part of the noise. The expected misfit is then roughly ${\sigma }^2(n-\text{ED})$, where ED, the effective number of parameters (also called effective degrees of freedom), is the trace (the sum of the diagonal elements) of the smoother matrix $\mathbf{H}=(\mathbf{I}+\lambda {\mathbf{D}}^T\mathbf{D}{)}^{-1}$, the matrix that turns $\mathbf{y}$ into $\hat{\mathbf{z}}$. In this panel, the lowest error against the true signal is reached near $\lambda \approx 5$, with ED about 49 and a misfit of about 0.5, while the λ that brings the misfit up to 0.72 gives an error about a third larger. So when the noise level of your instrument is known, use $n{\sigma }^2$ as an upper bound for the misfit, not as a target. The trade-off curve in the panel marks both the $n{\sigma }^2$ landmark and the λ with the lowest error, so you can see the gap.

For a long, evenly spaced signal, the Whittaker smoother acts like a filter that keeps a wave of frequency $f$ (in cycles per point) with gain

Slow waves (small $f$) have a gain close to 1 and pass; fast waves are strongly damped. The gain never becomes negative, so, unlike a moving average, it never flips noise components. Larger $\lambda$ moves the cut-off to lower frequencies, and larger $d$ makes the transition between kept and removed frequencies steeper.

This formula turns the abstract λ into something you can picture. The gain drops to one half when $\lambda (2\sin \pi f{)}^{2d}=1$. For $d=2$ and $\lambda =100$ this happens at $f_{50}=\frac{1}{\pi }\arcsin ({\lambda }^{-1/4}/2)\approx 0.05$ cycles per point: waves longer than about 20 points are mostly kept, shorter ones are mostly removed. A handy approximation: waves with a period of $P$ points are halved when $\lambda \approx (P/2\pi {)}^{2d}$, which for $d=2$ and $P=20$ gives again about 100.

A powerful extension gives each point a weight $w_i$:

whose solution is

with $\mathbf{W}$ a diagonal matrix holding the weights. Points with a large weight pull the curve strongly; points with a small weight barely matter. A weight of zero means the point is ignored entirely, which lets the smoother fill gaps: inside a gap only the penalty acts, so the curve continues as smoothly as possible [2].

Inside the gap the result is a smooth polynomial curve joining the two sides (a straight line for $d=1$, a cubic curve for $d=2$, a quintic for $d=3$). The kind of curve depends only on $d$; $\lambda$ changes how it meets the data at the edges of the gap. The grey dotted line in the panel is a plain straight-line interpolation, for comparison. The smoother cannot invent a band that was never measured, but it bridges the gap more naturally than a straight line when $d\ge 2$.

Weights are also the route to baseline correction. In asymmetric least squares, points above the current curve (likely to be bands) get a small weight and points below get a large weight, so after a few iterations the curve follows the baseline under the bands [4]. See Baseline correction with AsLS.

There is no universal value. The best $\lambda$ depends on:

• the noise level,
• the width of the narrowest feature you want to keep, measured in data points,
• the difference order $d$ (higher orders need much larger $\lambda$ for similar smoothing),
• your application: exploration or quantitative modelling.

Scan on a log scale. Try $\lambda =1,10,100,\dots$ and look at the results and the residuals.

Visual inspection:

1. Plot the original and smoothed signals.
2. Increase $\lambda$ gradually.
3. Stop when the noise is gone but the narrowest bands still keep their height.

Residual analysis:

1. Compute the residuals $r_i=y_i-z_i$.
2. If they look like random noise, the smoothing is appropriate.
3. If they show structure (band shapes, trends), you are over-smoothing.

Cross-validation: leave out one point at a time, smooth the rest, and see how well the left-out point is predicted. Eilers showed that for this smoother the leave-one-out error can be computed without refitting: the leave-one-out residual of point $i$ is $r_i/(1-h_{ii})$, where $h_{ii}$ are the diagonal elements of the "hat" matrix $\mathbf{H}=(\mathbf{W}+\lambda {\mathbf{D}}^T\mathbf{D}{)}^{-1}\mathbf{W}$, so $\lambda$ can be chosen automatically [2]. Garcia later combined a closely related criterion, generalised cross-validation, with a fast solver based on the discrete cosine transform, a popular way to pick λ automatically [7], and Weinert described efficient algorithms for the computation itself [8]. Be aware that cross-validation can choose too little smoothing when the noise is correlated from point to point.

For calibration models: smooth with several $\lambda$ values, build the model on each, and choose the value that gives the best cross-validated prediction error.

import numpy as np
import matplotlib.pyplot as plt
from scipy import sparse
from scipy.sparse.linalg import spsolve


def whittaker_smooth(y, lam=1e2, d=2, w=None):
    """
    Whittaker smoother (Eilers, 2003).

    y   : 1D array, the signal to smooth
    lam : smoothing parameter lambda (larger = smoother)
    d   : order of the difference penalty (usually 2)
    w   : optional weights (0 = ignore the point); default all ones
    """
    y = np.asarray(y, dtype=float)
    n = len(y)
    if w is None:
        w = np.ones(n)

    # d-th order difference matrix, shape (n - d, n)
    D = sparse.eye(n, format='csr')
    for _ in range(d):
        D = D[1:] - D[:-1]

    W = sparse.diags(w)
    A = (W + lam * D.T @ D).tocsc()
    return spsolve(A, w * y)


# Noisy synthetic spectrum
rng = np.random.default_rng(42)
wavelength = np.linspace(400, 800, 400)
true_signal = (2 * np.exp(-((wavelength - 500)**2) / 500) +
               1 * np.exp(-((wavelength - 600)**2) / 300))
noisy = true_signal + rng.normal(0, 0.1, len(wavelength))

plt.figure(figsize=(10, 6))
plt.plot(wavelength, noisy, color='gray', alpha=0.5, label='Noisy')
for lam in [1e1, 1e3, 1e5]:
    plt.plot(wavelength, whittaker_smooth(noisy, lam=lam), linewidth=2,
             label=f'λ = {lam:g}')
plt.plot(wavelength, true_signal, 'k--', label='True signal')
plt.xlabel('Wavelength (nm)')
plt.ylabel('Absorbance')
plt.legend()
plt.show()

# Gap filling: give a stretch of points weight 0
w = np.ones(len(noisy))
w[180:230] = 0
filled = whittaker_smooth(noisy, lam=1e3, d=2, w=w)

function z = whittaker_smooth(y, lam, d, w)
    % Whittaker smoother (Eilers, 2003)
    %   y   : signal (vector)
    %   lam : smoothing parameter lambda (larger = smoother)
    %   d   : order of the difference penalty (usually 2)
    %   w   : optional weights (0 = ignore the point)
    if nargin < 3, d = 2; end
    y = y(:);
    n = length(y);
    if nargin < 4, w = ones(n, 1); end

    D = diff(speye(n), d);            % sparse (n-d) x n difference matrix
    W = spdiags(w(:), 0, n, n);
    z = (W + lam * (D' * D)) \ (w(:) .* y);   % sparse Cholesky behind the backslash
end

% Example (in a separate script)
rng(42);
wavelength = linspace(400, 800, 400)';
true_signal = 2 * exp(-((wavelength - 500).^2) / 500) + ...
              1 * exp(-((wavelength - 600).^2) / 300);
noisy = true_signal + 0.1 * randn(size(wavelength));

figure; plot(wavelength, noisy, 'Color', [0.6 0.6 0.6]); hold on;
for lam = [1e1, 1e3, 1e5]
    plot(wavelength, whittaker_smooth(noisy, lam, 2), 'LineWidth', 2);
end
plot(wavelength, true_signal, '--k', 'LineWidth', 1.5);
legend('Noisy', '\lambda = 10', '\lambda = 10^3', '\lambda = 10^5', 'True');
xlabel('Wavelength (nm)'); ylabel('Absorbance'); grid on;

library(Matrix)  # sparse matrices

whittaker_smooth <- function(y, lambda = 1e2, d = 2, w = rep(1, length(y))) {
  # Whittaker smoother (Eilers, 2003)
  n <- length(y)
  E <- Diagonal(n)
  D <- diff(E, differences = d)          # sparse (n - d) x n difference matrix
  A <- Diagonal(x = w) + lambda * crossprod(D)
  as.vector(solve(A, w * y))
}

set.seed(42)
wavelength <- seq(400, 800, length.out = 400)
true_signal <- 2 * exp(-((wavelength - 500)^2) / 500) +
               1 * exp(-((wavelength - 600)^2) / 300)
noisy <- true_signal + rnorm(length(wavelength), sd = 0.1)

plot(wavelength, noisy, type = "l", col = "gray",
     xlab = "Wavelength (nm)", ylab = "Absorbance", main = "Whittaker smoothing")
cols <- c("blue", "red", "purple")
lams <- c(1e1, 1e3, 1e5)
for (k in seq_along(lams)) {
  lines(wavelength, whittaker_smooth(noisy, lams[k]), col = cols[k], lwd = 2)
}
lines(wavelength, true_signal, lty = 2, lwd = 1.5)
legend("topright", legend = c("Noisy", paste("lambda =", lams), "True"),
       col = c("gray", cols, "black"), lty = c(1, 1, 1, 1, 2))

It is an excellent choice when:

• You want very good band preservation controlled by one parameter.
• Your data have gaps, missing points or unreliable regions.
• You plan to use baseline correction methods built on the same machinery.
• You need automatic smoothing, with $\lambda$ chosen by cross-validation.

Prefer another method when:

• You need derivatives: use Savitzky-Golay.
• You need the simplest possible explanation: use a moving average.

• Spectroscopy: denoising NIR, Raman and IR spectra, and as the core of baseline correction.
• Chromatography: smoothing detector traces and estimating drifting baselines.
• Remote sensing: smoothing and gap-filling vegetation index time series from satellites, where clouds create missing values.
• Economics: the Hodrick-Prescott trend filter is the same smoother.
• Actuarial science: graduation of mortality tables, where it all started.

• λ too small: noise remains.
• λ too large: bands are flattened and broad features bend towards a line.
• Different λ for calibration and test data: inconsistent preprocessing.
• Using d = 1 without a reason: band tops and slopes get flattened.
• Using dense matrices: memory and time grow with $n^2$ or worse; always use sparse matrices.
• Not checking the residuals: they are the quickest test of over-smoothing.

• Whittaker smoothing finds the curve that minimises misfit plus $\lambda$ times roughness, with roughness measured by squared differences of order $d$.
• The solution is one sparse linear system, $(\mathbf{W}+\lambda {\mathbf{D}}^T\mathbf{D})\mathbf{z}=\mathbf{Wy}$, solved in time proportional to the number of points.
• $\lambda$ is explored on a logarithmic scale; very large values give a polynomial of degree $d-1$.
• Weights let the smoother ignore missing points and form the basis of baseline correction.

• Moving average: the simplest smoother.
• Savitzky-Golay smoothing: local polynomial fits and derivatives.
• Gaussian smoothing: weighted averaging with a bell-shaped kernel.
• Baseline correction with AsLS: Whittaker smoothing with asymmetric weights.

For most spectroscopic work, Whittaker and Savitzky-Golay are the two workhorses. Whittaker gives precise control through $\lambda$ and handles gaps; Savitzky-Golay excels at smoothed derivatives. Try both on your own data.