Whitened Repack Invariance: A Generalized Least-Squares Extension of the Lagrangian Repack-Correction Algorithm for Diffuse Reflection Clutter Suppression in NIR Calibration

Column #251, Part II of a series continuing "Development of an Algorithm That Minimizes/Eliminates Sampling Repack Variation" (Column #250)

Abstract

The preceding column (#250) in this series derived, from Lagrange's method of undetermined multipliers, a constrained multiple linear regression (MLR) algorithm whose coefficients satisfy Σ b_iΔA_i = 0, rendering the predicted constituent value insensitive to a single, assumed-proportional direction of spectral change caused by sample repacking. That single-constraint assumption is a special case of a richer phenomenon: diffuse-reflection variability from particle size, packing density, surface texture, and specular contamination is generally multidimensional and only approximately proportional across wavelengths. Independently, Gallagher and Watson formalized clutter suppression for classical least squares (CLS) through generalized least squares (GLS, covariance weighting) and extended least squares (ELS, explicit clutter basis vectors), including a clutter-suppression preprocessing step for inverse least squares (ILS) models such as MLR and partial least squares (PLS). This paper unifies the two approaches. We define a clutter covariance matrix W_c from replicate repack spectra, use it to whiten both the calibration absorbance matrix and the repack-difference vector or vectors, and show that applying the earlier Lagrangian bordered-normal-equations algorithm to the whitened data yields a diffuse-reflection-invariant calibration that reduces exactly to the original algorithm when the clutter covariance is isotropic and a single repack direction is used, extends naturally to multiple, non-proportional clutter directions through a multi-multiplier bordered system, has a soft, diagnostic-rich counterpart based on an ELS estimator, and, via an EPO-filter route that instead reproduces DR-GLS exactly, is compatible with PCR and PLS. We provide the full derivation, a compact worked illustration, a proposed experimental validation protocol using the wheat repack data referenced in the earlier column, and a discussion connecting the clutter subspace to Kubelka-Munk scattering theory and to established scatter-correction methods such as MSC, SNV, and EPO.

Keywords: diffuse-reflection correction, repack effect, clutter suppression, generalized least squares, extended least squares, external parameter orthogonalization, Lagrange multipliers, NIR calibration, multivariate statistics

Nomenclature and Notation

Our preceding column¹ and Gallagher and Watson's white paper² each grew their own notation, and the two do not always agree: both use C, for example, for two different quantities. Rather than force the reader to keep two cheat sheets, we reconcile them here: we follow the preceding column wherever the two sources agree, and follow Gallagher and Watson where they do not. Matrices and vectors are set in italic capital and lowercase roman letters, following both source documents; boldface is reserved for emphasis in prose. Matrix equations are written in bracket notation with rows separated by a semicolon (for example, [X; Y] denotes a matrix with row X above row Y), since the underlying operations are identical to, but more compact than, the fully expanded scalar form we used for illustration last time.

Introduction

Two columns ago we asked a question that has floated around the NIR community since the earliest days of the technology: could a calibration algorithm be built that is mathematically blind to "repack" (the systematic but uncontrollable change in a diffuse reflectance spectrum that shows up when a powdered or granular sample is removed from its cup and reloaded) while staying fully alert to genuine compositional differences? Working from Lagrange's method of undetermined multipliers,³ that column¹ derived a constrained MLR algorithm whose regression coefficients b satisfy a single linear constraint tied to one measured repack-difference vector ΔA. The resulting bordered normal-equations system, solved exactly like ordinary MLR, gives a solution that is, by construction, unaffected by spectral changes proportional to ΔA.

That derivation leaned on one simplifying assumption, and we said so plainly at the time: that repack-induced changes at different wavelengths stay in constant proportion to one another for a given sample, that is, that a single direction ΔA adequately describes the clutter. We even showed a figure of repack spectra from 31 ground hard red wheat samples in which that proportionality holds only roughly, and we admitted that stretching the method to PCR or PLS, or to more complicated, non-proportional clutter, was an open problem. Our exact words were that it made "a good project for the Chemometrics community to tackle."

As it happens, somebody had already been tackling a closely related problem, just not our problem. Working independently, and for a different application altogether (target detection and quantification in hyperspectral and process spectroscopy), Neal Gallagher and Nathanial Watson² formalized the general question of separating a signal of interest from structured interference, which they call clutter, inside the classical least squares (CLS) family of models. Their generalized least squares (GLS) and extended least squares (ELS) estimators use a clutter covariance matrix, or an explicit clutter basis, to down-weight or orthogonalize a fit against directions of interference that carry no target signal. Buried in their concluding section is a short but important remark: the same clutter-suppression weighting extends to inverse least squares (ILS) models such as MLR and PLS, simply by preprocessing the predictor spectra with the inverse square root of the clutter covariance.

That remark is what caught our attention, because it is essentially the missing piece of our own puzzle. This column shows that the two developments are not just cousins in spirit but are formally the same idea seen from two angles, and that the connection resolves the open question we left dangling: a single proportional repack direction turns out to be the degenerate, rank-one case of a general clutter covariance; our original Lagrange-multiplier hard constraint is one of two natural ways of removing a clutter subspace from an ILS calibration, the other being Gallagher and Watson's softer, estimated-coefficient ELS filter; and the combined algorithm works with any linear calibration method, PLS and PCR included, because it lives entirely in a preprocessing (whitening) step rather than in the regression itself.

We will not re-review the broader literature on the repack effect, on multivariate calibration generally, or on the statistical tools used below; that background is already gathered in the compendia and primary sources this series has leaned on before,^4-9 and we cite those again only where a specific point is drawn from them.

Background

A Quick Recap: The Lagrangian Repack Constraint

We will not repeat the full scalar derivation here; it took the better part of the last column,¹ worked out in detail for a three-wavelength model. The short version goes like this. Write the ordinary MLR model for the ith wavelength coefficient as

C = b₀ + b₁A₁ + b₂A₂ + b₃A₃ + … + error, [1]

and require that two readings of the same sample, one before repacking and one after, give the same predicted constituent value (C = C′). Because b0 cancels, that requirement reduces to the linear constraint b1ΔA1 + b2ΔA2 + b3ΔA3 = 0, where ΔAi = Ai - A′i. Appending this constraint to the least-squares objective by Lagrange's method and differentiating gives a bordered (KKT) linear system that reads equation 2:

which is identical in form to the ordinary MLR normal equations, except that the coefficient matrix is bordered with the constraint direction ΔA and an extra unknown, λ. Its solution follows the same block-elimination procedure worked out in full in the Appendix for the general, whitened case; setting the clutter covariance to the identity there recovers the unweighted solution used here, no new machinery required. We also defined a diagnostic statistic, D = Σ(bΔA), which should land at zero for a perfectly repack-invariant model, and recommended that the ΔA used in the constraint be the largest, most repack-dominated, least noise-dominated difference vector available in the data set.

Classical, Weighted, Generalized, and Extended Least Squares

Gallagher and Watson² start from the classical least-squares (CLS) linear-mixture model X = CSᵀ + E and its weighted single-measurement form x = Sc + e, minimizing the weighted objective eᵀW⁻¹e. Four choices of the weighting matrix W give four named models: ordinary CLS (W = σ²I, uncorrelated and equal-variance noise), weighted least squares or WLS (W diagonal, unequal variance by channel), generalized least squares or GLS (W = W_c, a full, generally non-diagonal clutter covariance matrix estimated from measurements Xc that contain clutter but no target signal), and extended least squares or ELS, which instead of weighting the objective adds explicit clutter basis vectors P as extra regressors, x = [S P][c; t] + e, and estimates the clutter loadings t directly. They also show that the ELS estimator can be rewritten as a GLS-type estimator using the projector

W_ELS⁻¹ = I − P(P^TP)⁻¹P^T [3]

(equation 3), which will be familiar to many readers as the External Parameter Orthogonalization (EPO) filter of Roger, Chauchard, and Bellon-Maurel.¹⁰ Tucked into their concluding section, Gallagher and Watson note that the same weighting extends to inverse least squares (ILS) models such as MLR and PLS, y = Xb + e, simply by preprocessing X with the inverse square root of the clutter covariance, y = XW_c^−1/2b + e; and they flag that the regression coefficients obtained in that preprocessed space are generally not the same coefficients one would get on raw data. That distinction turns out to matter a great deal for what follows.

A word on credit, though, since none of this sprang from nowhere. GLS and WLS themselves are old news: Aitken had weighted and generalized least squares in essentially modern form back in 1936,⁴ and Gallagher and Watson's own bibliography pushes the lineage further back still, to work they cite as their refs. 5 and 6 (their numbering, not ours). So what exactly do we think Gallagher and Watson added? Not covariance weighting itself, but a single framework that puts CLS, WLS, GLS, and ELS on the same footing as special cases of one weighted objective, finally letting the family be compared, and extended to ILS, side by side. That framework is what we build on here.

The ILS extension has a track record, too. Gallagher and Watson trace that reasoning back further still, to their ref. 19, and a covariance-weighted preprocessing step much like this one has quietly sat inside commercial chemometrics software, Eigenvector Research's PLS_Toolbox among others, for on the order of two decades, handling general spectral clutter the way we're asking it to handle repack here. What's new in this column isn't the weighting mechanism. It's marrying that mechanism explicitly to the Lagrangian repack constraint from our earlier column, and building out the multidirectional, PLS/PCR-compatible version we derive below.

Theory: A Unified Diffuse-Reflection Clutter-Suppression Algorithm

Diffuse-Reflection Variability as Clutter

Following Gallagher and Watson,² we define clutter as any systematic, wavelength-correlated spectral signal that carries no information about the constituent being calibrated for. Diffuse-reflection variability (repack, particle-size distribution, packing density, surface texture, and the specular, or first-surface, reflectance component) is clutter in exactly this sense. It is real, physically caused, systematic, and correlated across wavelengths (our own repack conditions¹ state that it should be linear with respect to itself across wavelengths, and that spectra from repacks of one sample should not cross each other), yet by construction it carries no compositional signal.

Last time¹ we characterized this clutter with a single vector, ΔA, taken from one pair of repack measurements. Here we generalize that idea by characterizing the clutter statistically, from a whole population of repack, particle-size, or packing-density replicate measurements, exactly the way Gallagher and Watson characterize clutter for target detection. Let a representative subset of the calibration samples be measured in replicate, two or more repacks each. For the mth such sample, define the repack-difference vector

ΔA_m,i [4]

with entries ΔA_m,i = A_m,i − A′_m,i at each wavelength i, and stack the difference vectors as the rows of a clutter matrix Xc (Mc × N, one row per replicate pair or per additional repack beyond the first). This is simply the multi-sample version of our earlier single ΔA, and it is exactly what Gallagher and Watson mean by measurements that include clutter signal but no target signal: repack replicates of the same sample have, by definition, ΔC = 0.

The Clutter Covariance Matrix

From Xc we estimate the diffuse-reflection clutter covariance, using the non-centered form that Gallagher and Watson identify as appropriate for reflectance spectra:

When Wc is ill-conditioned or rank-deficient, which it typically will be when Mc is modest relative to N (a common situation in a repack study), we regularize it exactly as Gallagher and Watson recommend, by ridge-type inflation toward the isotropic case:

W_c → W_c + α I, α ≈ σ_x² [6]

with α on the order of the mean single-channel noise variance. This one regularization step is what lets the method fall back gracefully to ordinary, unweighted MLR as the amount of characterized clutter shrinks toward zero: precisely the isotropic W = σ²I case in which Gallagher and Watson's GLS collapses back to plain CLS.

Where more than one physically distinct source of diffuse-reflection variability is at play (particle size and packing density, for instance, behave quite differently), principal component analysis of Xc yields an orthonormal basis P (N × L) for the dominant L directions of the clutter subspace, Xc ≈ TPᵀ. We come back to this multidirectional case in a later section; the physical interpretation of the loadings in P is taken up in the Discussion.

Weighted Clutter and a Link to Time-Series GLS

Equation 5 is not the only way to build a clutter covariance. Time-series applications often construct GLS clutter from differences between consecutive measurements, on the logic that whatever drifts slowly from one reading to the next is probably instrumental or environmental noise rather than real signal. Repack has a close analogue. Instead of treating every replicate pair in X_c as equally informative, weight the difference vectors by how close the two measurements' reference values are to each other. Repacks of samples with nearly identical reference values get weighted up, since a difference between them is more likely to be pure clutter; repacks whose reference values sit farther apart get weighted down, because some of that gap could be real composition, not repack. This reference-value-weighted W_c has been used in practice for some time, and works well in at least some instances, though we have not seen it published or tied explicitly to the CLS/GLS framework the way the simpler, unweighted W_c is above. Consider it a candidate refinement, worth testing alongside the baseline W_c of equation 5 in the validation protocol below.

Whitening and the Generalized Lagrangian Constraint

Here is the key step, and the original contribution of this paper. Gallagher and Watson² already showed how to apply ILS clutter-suppression preprocessing to a calibration absorbance matrix. What we add is applying that same transformation to the repack-difference vector or vectors that define our earlier Lagrangian constraint,¹ so that the constraint and the data it acts on are expressed in the same coordinate system.

Define the whitening transform using the symmetric inverse square root of the regularized clutter covariance:

à = A W_c^−1/2, ΔÃ_m = W_c^−1/2 ΔA_m [7]

and reparameterize the regression coefficients as b = W_c^−1/2g, so that Ãg = Ab for any g. Because W_c^−1/2 is symmetric, our original repack-invariance requirement, bᵀΔA = 0, becomes, in the whitened coordinates,

g^TΔÃ = 0 [8]

which is formally identical to our original constraint, just written for g and ΔÃ instead of b and ΔA. That means the entire bordered-normal-equations machinery from the earlier column carries over unchanged; only the data it operates on has been pre-whitened. Writing the bordered system in real matrix form,

whose closed-form, Schur-complement solution (derived in full in the Appendix, following the same block-inversion logic Gallagher and Watson use to pass from their equation 21 to equation 22) is just the classical restricted-least-squares estimator¹¹ applied to the whitened data:

ĝ = (Ã^TÃ)⁻¹Ã^TC − (Ã^TÃ)⁻¹ΔÃ[ΔÃ^T(Ã^TÃ)⁻¹ΔÃ]⁻¹ΔÃ^T(Ã^TÃ)⁻¹Ã^TC [10]

with the associated Lagrange multiplier

λ = [ΔÃ^T(Ã^TÃ)⁻¹ΔÃ]⁻¹ΔÃ^T(Ã^TÃ)⁻¹Ã^TC [11]

Finally, we return the coefficient vector to the original absorbance units, so it can be applied to new, un-whitened spectra at prediction time:

b̂ = W_c^−1/2 ĝ [12]

We call this the Diffuse Reflection Generalized Least Squares (DR-GLS) estimator, and two things about it are worth spelling out. First, when the clutter is isotropic (W_c = σ²I, meaning no directional clutter has been characterized beyond ordinary noise), W_c^−1/2 is just a scalar multiple of the identity, Ã ∝ A and ΔÃ ∝ ΔA, and equation 10 collapses exactly, term for term, into our original bordered solution. DR-GLS, in other words, is not a replacement for the earlier algorithm but a strict generalization of it. Second, because the clutter weighting enters only through the whitening step in equation 7, the same reasoning that justifies equation 10 for MLR justifies using it as a preprocessing step ahead of PCR or PLS. Ã and ΔÃ can be handed to any linear calibration method, which is exactly the extension we could not offer last time.

Multiple, Non-Proportional Clutter Directions

Our earlier repack conditions¹ assumed a single dominant direction of proportional change. Where diffuse-reflection clutter has more than one physically distinct mode (a particle-size-driven baseline-and-tilt component together with an independent packing-density offset, say), the PCA loadings P (N × L) from the previous section replace the single vector ΔA. Whitening P by the same transform, P̃ = Wc^−1/2P, lets the bordered system of equation 9 generalize directly to L simultaneous constraints and L Lagrange multipliers,

solved by the identical block-matrix procedure used for equations 10-11, with ΔÃ replaced by the matrix P̃. This is exactly the "obvious" expansion we predicted would be needed for more wavelengths; it turns out the same expansion also handles more constraints.

A Softer Alternative: Diffuse Reflection Extended Least Squares (DR-ELS)

The Lagrange approach above imposes an exact algebraic constraint: the fitted model is forced to have zero sensitivity to the specified clutter direction or directions, whatever that costs the fit elsewhere. Gallagher and Watson's ELS estimator offers a gentler alternative, in which the clutter loadings are estimated as ordinary regression coefficients alongside the target coefficients, rather than nulled outright by a multiplier. Because the design matrix must have one row per calibration sample, the clutter block enters through the whitened calibration spectra's own projections onto the clutter directions, ÃP̃ (M × L), rather than through the loadings P̃ themselves (N × L). Applied to whitened diffuse-reflection data, the model becomes

with ordinary least-squares solution

A related way to remove the clutter contribution is by filtering à with an EPO-type projector before an ordinary fit:

ĝ = (W_ELS⁻¹Ã^TÃW_ELS⁻¹)⁻¹W_ELS⁻¹Ã^TC, W_ELS⁻¹ = I − P̃(P̃^TP̃)⁻¹P̃^T [16]

The augmented-design estimator of equations 14-15 gives up the hard repack-invariance guarantee of DR-GLS in exchange for two things our original algorithm could never offer. The estimated clutter coefficients t̂ are themselves a diagnostic (a large t̂ flags a sample whose measured spectrum is dominated by diffuse-reflection artifact rather than composition), and the correction degrades gracefully, rather than catastrophically, if the assumed clutter subspace P turns out to be only approximately right. The equation-16 filter, in contrast, gives up nothing: because it solves the identical equality-constrained least-squares problem as the bordered system of equations 9-13, it reproduces the DR-GLS solution exactly, with no t̂ diagnostic and no relaxation of the invariance guarantee; its advantage is purely computational, since it is what extends the correction to PLS and PCR below. Our own preference is DR-GLS, computed either via the bordered system of equations 9-13 or, equivalently, via the equation-16 filter, where an exact invariance guarantee is required, for example in a regulatory or referee-facing calibration, and the augmented-design estimator of equations 14-15 where diagnostic value and tolerance for an imperfectly characterized clutter subspace matter more.

Combining a Known Interferent with Adaptive Clutter Weighting

Gallagher and Watson point out, in their equation 26, that a known, fixed interferent filter and a separately estimated, updatable GLS weighting can be nested together. The diffuse-reflection version of that idea is handy when, say, cup geometry or a specular-reflectance component is known and stable across an instrument fleet (so P̃ and WELS⁻¹ stay fixed) while the finer repack and particle-size clutter covariance Wc is re-estimated periodically from fresh repack replicates:

b̂ = W_c^−1/2W_ELS⁻¹(W_ELS⁻¹W_c^−1/2A^TAW_c^−1/2W_ELS⁻¹)⁻¹W_ELS⁻¹W_c^−1/2A^TC [17]

As Gallagher and Watson note for the analogous hyperspectral case, this is exactly the sense in which clutter-suppression correction is "adaptive": the fixed geometric filter need not be recomputed every time new repack data come in, even though Wc can be.

One more credit check, on a smaller point. Gallagher and Watson's equation 26 nests a fixed interferent filter inside an adaptively re-estimated GLS weighting; our equation 17 is that same nesting, specialized to a fixed geometric or specular filter combined with a repack-and-particle-size covariance that gets updated periodically. We could not find published work applying this fixed-plus-adaptive combination to diffuse-reflection or repack correction specifically. Treat it, accordingly, as part of the proposed algorithm rather than as settled practice.

Gallagher and Watson's equation 26 formalizes the concept. It was shown in a poster on multivariate curve resolution applied to surface-enhanced Raman spectroscopy (SERS) measurements of 2,4-dinitrotoluene (2,4-DNT)¹² and used algorithmically in an automated atmospheric-compensation and target-detection algorithm for standoff hyperspectral detection.¹³

Diagnostics

Two diagnostics fall directly out of the estimation-error and residual expressions Gallagher and Watson give for the general weighted CLS, GLS, and ELS family, specialized here to the ILS case. The weighted residual for the mth calibration or validation sample,

q_m = ẽ_m^Tẽ_m, ẽ = C − Ã ĝ [18]

flags samples whose whitened residual is unexpectedly large, which tells us either that the assumed clutter subspace does not fully describe that sample's diffuse-reflection behavior, or that we are looking at a genuine outlier. This is the same caution Gallagher and Watson raise: estimates with high residuals should not be trusted, because they may be biased. Second, our earlier repack-invariance statistic generalizes directly to a per-sample check across the whole calibration and validation sets,

which should stay small for every sample with characterized repack replicates, not merely for the single pair we used originally. A systematic, nonzero Dm across a subset of samples is itself diagnostic: it tells us that those samples' diffuse-reflection behavior lies outside the clutter subspace spanned by Wc or P, and gives us a direct, quantitative reason to enlarge L or collect additional repack replicates for that sample type.

A Visual Illustration: Applying the Correction to Simulated Repack Spectra

Before working through the scalar case, it helps to see what the whitening step of equation 7 actually does to a set of repack spectra. Figure 1 shows six simulated repack measurements of a single sample, built from one underlying absorbance spectrum with two overlapping NIR bands, each repack perturbed by an independent multiplicative scattering factor and a small linear tilt against wavelength, the same two clutter modes discussed above (a particle-size-driven scattering multiplier and a packing-density-driven tilt). These are, deliberately, simulated data meant only to illustrate the mechanism; they are not a substitute for the wheat repack data validation proposed below.

Figure 1. Conceptual illustration of diffuse-reflection clutter suppression. Left: six simulated repack spectra of the same sample before correction, showing the multiplicative and tilt-like scatter variation typical of repack, particle-size, and packing-density effects. Right: the same six spectra after projecting out the two-dimensional clutter subspace spanned by the multiplicative and tilt modes, the same conceptual operation performed by the Wc^−1/2 whitening of equation 7 ahead of the bordered-system solve of equation 9. The mean across-repack standard deviation drops roughly 16-fold after correction, illustrating why a repack-invariant calibration built on the corrected spectra should be far less sensitive to which particular repack was measured.

One figure only goes so far. Readers with their own repack or clutter data might want to build the same before-and-after picture for their own instrument and sample type; that, along with a worked case on the wheat repack data itself, is on the list of extensions below. And it's worth saying plainly: a clutter-corrected spectrum doesn't always look, to the eye, like an obviously "better" or more natural spectrum. The correction pulls a spectrum toward its net analyte signal, and net analyte signal is not the same thing as visually interpretable, even though it is the representation most tightly tied to the constituent being predicted, a point Gallagher and Watson draw from still earlier work of their own. If a correction looks a little strange, that alone doesn't mean something went wrong.

Illustrative Case: A Diagonally Weighted Three-Wavelength Model

To connect equations 7-12 to the fully expanded scalar derivation we walked through last time,¹ consider the same three-wavelength illustration, but now suppose repack replicates show that channel 2 is roughly twice as clutter-prone as channels 1 and 3 (perhaps because it sits nearer a particle-size-sensitive combination band). A diagonal clutter covariance capturing this is Wc = diag(σ², 4σ², σ²), so that Wc^−1/2 = diag(σ⁻¹, ½σ⁻¹, σ⁻¹). Whitening simply rescales each absorbance channel, and the corresponding entry of ΔA, by these factors:

with ΔÃ1 = ΔA1/σ, ΔÃ2 = ΔA2/2σ, and ΔÃ3 = ΔA3/σ transformed the same way. Substituting into the constraint of equation 8 shows the mechanism plainly: the whitened constraint

g₁ΔÃ₁ + g₂ΔÃ₂ + g₃ΔÃ₃ = 0 [21]

is identical in form to our original constraint,¹ but because channel 2 has been shrunk by a factor of two relative to channels 1 and 3, the least-squares split of the constraint burden across wavelengths that comes out of solving equation 10 automatically leans less on channel 2 to satisfy repack invariance. That is exactly the behavior we want when channel 2 is known, from the repack replicates, to be a noisier or more clutter-dominated carrier of the constraint. Un-whitening at the end, equation 12, returns coefficients in ordinary absorbance units, b1 = ĝ1/σ, b2 = ĝ2/2σ, b3 = ĝ3/σ, ready to apply to a new, raw, unwhitened sample spectrum. The scalar bookkeeping is otherwise identical, term for term, to the derivatives we worked through in full last time; only the constants multiplying each Ai,k differ.

Proposed Algorithm

Here is the procedure in practical terms: DR-GLS is the hard-constraint route, and, in brackets, we note where DR-ELS (the softer, diagnostic-rich alternative) departs from it, for a calibration set of M samples at N wavelengths.

Discussion

Physical Interpretation of the Clutter Subspace

Last time¹ we deliberately avoided committing to a physical theory of diffuse reflection, noting that neither physical theories of fundamental optics nor optical theories based on diffusion properties of light have fully explained the phenomenon, and pointing to the long lineage from Schuster^{14,\ 15} through Kubelka and Munk¹⁶ to the Representative Layer theory of Dahm and Dahm^17-27 and the discontinuum treatments reviewed there. We are not going to settle that debate here, but the clutter-covariance framework does give it a precise statistical role: whatever combination of scattering coefficient S(λ), absorption coefficient K(λ), packing fraction, and particle-size distribution produces the repack effect, its net contribution to the measured spectrum is exactly what Xc captures empirically, with no need to know the underlying physics in closed form. This mirrors the Kubelka-Munk relation itself,

equation 22, in which S, the scattering coefficient primarily responsible for particle-size and packing-density sensitivity, enters multiplicatively and in a wavelength-dependent way. That is consistent with what we saw in the Illustrative Case above: different channels can be differently clutter-prone, and therefore benefit from unequal weighting rather than the single uniform proportionality we assumed for simplicity last time. In practice we would expect the leading PCA loadings of Xc to correspond to recognizable physical modes: a broad, low-frequency loading dominated by S(λ), reflecting particle-size and packing offset and tilt; a narrower loading associated with specular or first-surface reflectance contamination; and, for granular or fibrous materials, higher-order loadings associated with preferred orientation effects on repacking. This lines up with how Gallagher and Watson describe PCA loadings of a clutter set as spanning the systematic variance in a hyperspectral background, and with the physical loadings sought in the Representative Layer theory's treatment of surface geometry.

A Rayleigh-Quotient View of the Objective Function

There's a second reason the Lagrangian route of our earlier column¹ and Gallagher and Watson's covariance-weighted route might be the same idea, and we have not chased it down to a full derivation. Flagging it anyway seems worthwhile. The CLS/GLS/ELS objective functions can be written as a ratio of quadratic forms in the fitted coefficients, i.e., as a Rayleigh quotient. Rayleigh quotients are stationary at the eigenvectors of a generalized eigenvalue problem, and that eigenvalue problem is itself usually derived by imposing a normalization constraint through Lagrange's method of undetermined multipliers, the same tool we reached for last time, just applied there to a simpler, single-constraint version of the problem. So the bordered normal-equations system of equation 9 might not just resemble a constrained eigenvalue problem; it might be one, in its degenerate, single-mode limit. We have not checked whether the connection actually reproduces equation 10 exactly, and that strikes us as good follow-up work. If it holds up, it would mean two sets of authors arrived, independently and without realizing it, at the same underlying eigenvalue problem from opposite directions.

The generalized eigenvector problem (GEP) mentioned above is associated with two square symmetric matrices that can be transformed to a symmetric eigenvector problem that can be solved with an SVD.

Relationship to Established Scatter-Correction Methods

DR-GLS and DR-ELS are best understood as a population-level, covariance-aware generalization of two families of scatter-correction preprocessing already in wide NIR use. Multiplicative scatter correction, or MSC,²⁸ regresses each individual sample spectrum against a reference, usually the mean, spectrum and removes the fitted multiplicative and additive terms. It is, in effect, a per-sample, rank-one ELS correction in which the clutter basis is estimated from a single reference spectrum rather than from replicate measurements, and it does not use any information about how the sample's own repacks vary. Standard normal variate transformation, or SNV,²⁹ row-standardizes each spectrum to zero mean and unit variance; it can be viewed as an extreme, purely diagonal, per-sample special case of W_c in which every wavelength is treated as equally and independently clutter-prone, with no cross-wavelength correlation exploited at all. External Parameter Orthogonalization,¹⁰ by contrast, is architecturally identical to DR-ELS (equation 16) and, as Gallagher and Watson note, already requires that the interfering parameter vary within the calibration set used to build P. DR-GLS and DR-ELS differ from all three chiefly in using a population, replicate-based clutter covariance or basis rather than a single reference spectrum or a diagonal assumption, and, in the DR-GLS case, in offering an exact rather than approximate invariance guarantee for the specific clutter direction or directions characterized.

One clarification, since the two are easy to conflate: clutter suppression, in every form discussed above, deals with spectral shape. It reweights or reprojects the wavelength-to-wavelength correlation structure of the interference. It does not correct magnitude effects, overall scaling, offset, or path-length differences, the kind of thing MSC, SNV, and plain normalization are built for; the two families of preprocessing target different aspects of diffuse-reflection variability and typically get used together, not as substitutes for each other. There's also a limiting case worth knowing about at the opposite end from MSC and SNV. As the estimated clutter covariance W_c becomes more nearly diagonal, meaning the clutter across wavelength channels grows more statistically independent, DR-GLS collapses toward an ordinary, per-channel weighted least squares, the WLS case defined earlier. A full, non-diagonal Wc earns its keep only when clutter is correlated across wavelengths, which for diffuse-reflection variability is the rule, not the exception.

Answering the Open Question We Left Last Time

We closed the previous column¹ with an honest admission: it was unknown whether a similar variation could be applied to calibrations based on PCR or PLS. Because DR-GLS and DR-ELS act entirely as a preprocessing transformation of the predictor matrix, Ã, or à with P̃ appended, rather than as a modification of the regression step itself, we think that question now has a direct answer: yes, by construction. Any linear calibration method that accepts a preprocessed spectral matrix (PCR, PLS, ridge regression, or plain MLR) can be substituted for the ordinary-least-squares step inside equation 10, equation 13, or equation 15 without altering the derivation. For DR-GLS specifically, the hard Lagrangian constraint of equation 8 remains available only when the downstream method is itself a linear regression solved from normal equations, such as MLR, ridge, or PCR. PLS's iterative deflation does not accept an equality-constrained inner solve in closed form, so for PLS we recommend the equation-16 filter instead, which requires only that Wc^−1/2 be applied to, and the clutter subspace projected out of, the spectra before the usual PLS algorithm runs, unmodified; because that filter reproduces the DR-GLS solution exactly, PLS calibrated this way retains the same exact repack-invariance guarantee as MLR or PCR under DR-GLS.

Recommendations for Experimental Validation

This paper is, deliberately, a theoretical extension in the same spirit as the column it follows,¹ whose own closing section listed confirmation on actual data as an open item and left it there. We have not fabricated or simulated numerical results here. Instead, we propose a concrete validation protocol using the same repack data we referenced previously (two repacks each of 31 ground hard red wheat samples) so that a direct, apples-to-apples comparison with the original Lagrangian algorithm is possible.

1. Baseline replication. Reproduce our original single-constraint result on the wheat data set (protein or moisture, whichever constituent value is available) to confirm D is small, establishing the performance floor against which DR-GLS and DR-ELS are compared.

2. Clutter characterization. If more than two repacks per sample are available, or can be collected, estimate Wc (equation 5) from the full replicate set rather than a single pair. Otherwise, construct Wc from the 31 available pairwise ΔA vectors, which already gives a population estimate that the original single-vector approach could never use.

3. Head-to-head comparison. Using an independent, held-out subset of samples and their own repacks, compare the original estimator against DR-GLS and DR-ELS on three fronts: the repack-invariance statistic D over all held-out repack pairs, not just the pair used to build the constraint; standard error of prediction, or SEP, on the held-out reference values; and the weighted residual diagnostic q (equation 18). As we cautioned previously, SEP and D should be reported and interpreted separately, since SEP reflects several error sources beyond repack.

4. Multi-direction check. Run principal component analysis on the 31 ΔA vectors. If more than one component carries appreciable variance, repeat the comparison using the L-direction DR-GLS or DR-ELS extension against the single-direction baseline, to test directly whether the single-proportionality assumption is what limits performance on this data set.

5. Cross-instrument or cross-particle-size generalization. Because Wc and P are estimated independently of the regression step, test whether a clutter covariance estimated on one particle-size grind, or one instrument, usefully transfers, through equation 17, to another. That is a question of real practical interest, and one the original single-pair constraint could not address at all.

Future Extensions of These Studies

1. Empirical validation of DR-GLS and DR-ELS on the wheat repack data set and on additional matrices, such as ground grain of varying particle size, pharmaceutical powders, and minerals, following the protocol above.

2. A systematic comparison of the estimated PCA loadings of Xc against Kubelka-Munk-predicted scattering-coefficient behavior and against the Representative Layer theory's treatment of surface geometry, to test whether the empirically derived clutter subspace has a recoverable physical interpretation.

3. Extension of the PLS-compatible DR-ELS preprocessing to a worked PLS case study, including comparison with EPO-PLS¹⁰ on the same data.

4. An investigation of whether the clutter covariance Wc itself drifts with instrument aging or environmental conditions, and whether the adaptive combination of equation 17 is an adequate model for that drift.

5. As we said last time, the field could use a standardized definition of repack and a community consensus protocol for generating repack replicates, which would materially improve the comparability of Wc estimated across laboratories.

Conclusions

Last time¹ we showed that a single, carefully chosen linear constraint, derived from Lagrange's two-hundred-year-old method of undetermined multipliers, can make an MLR calibration exactly insensitive to one measured direction of repack-induced spectral change. Gallagher and Watson,² working independently and for a different purpose, showed that structured spectral interferences, or clutter, can be suppressed in least-squares models either by covariance weighting (GLS) or by explicit orthogonalization against an estimated clutter basis (ELS), and sketched how the same weighting extends to inverse least squares models such as MLR and PLS. This paper has argued that these are not two unrelated tools but two views of one problem. Our earlier Lagrangian repack constraint is the exact, hard-constraint, rank-one special case of a clutter-whitened, generally multidirectional correction that we have now derived in closed form, DR-GLS, together with a softer, diagnostic-rich counterpart, DR-ELS, built directly from Gallagher and Watson's ELS estimator, together with an equation-16 EPO-filter route that instead reproduces DR-GLS exactly. Because the correction is implemented as a whitening preprocessing step rather than as a change to the regression itself, it works with MLR, PCR, and, by way of the equation-16 EPO-filter route, PLS, which finally answers the extension question we left open. We want to be explicit about what we have not done: we have not demonstrated empirical superiority over the original algorithm on real data. That demonstration, on the same wheat repack data set referenced throughout this series, is the necessary next step, and we hope to bring it to you soon.

Acknowledgement

The authors acknowledge the contributions of Drs. Neal Gallagher and Nathanial Watson to the concepts found in the article and also to their insightful review comments.

References

(1) Mark, H.; Workman, J., Jr. Development of an Algorithm That Minimizes/Eliminates Sampling Repack Variation. Spectroscopy 2026. https://www.spectroscopyonline.com/view/development-of-an-algorithm-that-minimizes-eliminates-sampling-repack-variation (accessed 2026-07-24).

(2) Gallagher, N. B.; Watson, N. Classical Least Squares and Clutter Suppression; Eigenvector Research, Inc.: Manson, WA. https://www.researchgate.net/publication/395446881_Classical_Least_Squares_and_Clutter_Suppression (accessed 2026-07-24). Revision Published as Classical Least Squares and Clutter Suppression. Spectroscopy Online, August 5, 2026. Available at: https://www.spectroscopyonline.com/view/classical-least-squares-and-clutter-suppression ((accessed 2026-08-05).

(3) Taylor, A. E. Advanced Calculus; Ginn and Company: New York, 1955.

(4) Aitken, C. On Least Squares and Linear Combination of Observations. Proc. R. Soc. Edinburgh 1936, 55, 42-48. DOI: 10.1017/S0370164600014346.

(5) Martens, H.; Næs, T. Multivariate Calibration, 2nd ed.; Wiley: Chichester, U.K., 1989.

(6) Mark, H.; Workman, J., Jr. Data Transforms in Chemometric Calibrations: Variation in MLR, Part 3: Reducing Sensitivity to Repack Σbi = 0. Spectroscopy 2023, 38 (10), 14-18. DOI: 10.56530/spectroscopy.jm1087r1.

(7) Mark, H.; Workman, J., Jr. Effect of Repack on Calibrations Produced for Near-Infrared Reflectance Analysis. Anal. Chem. 1986, 58 (7), 1454-1459. DOI: 10.1021/ac00298a041.

(8) ASTM International. Standard Practices for Infrared Multivariate Quantitative Analysis; ASTM E1655-05; West Conshohocken, PA, 2005. https://www.astm.org/DATABASE.CART/HISTORICAL/E1655-05.htm (accessed 2026-07-24).

(9) Ciurczak, E. W.; Igne, B.; Workman, J., Jr.; Burns, D. A., Eds. Handbook of Near-Infrared Analysis, 4th ed.; CRC Press: Boca Raton, FL, 2021. DOI: 10.1201/b22513.

(10) Roger, J. M.; Chauchard, F.; Bellon-Maurel, V. EPO-PLS External Parameter Orthogonalisation of PLS Application to Temperature-Independent Measurement of Sugar Content of Intact Fruits. Chemom. Intell. Lab. Syst. 2003, 66 (2), 191-204. DOI: 10.1016/S0169-7439(03)00051-000051-0).

(11) Rao, C. R. Linear Statistical Inference and Its Applications, 2nd ed.; Wiley: New York, 1973. DOI: 10.1002/9780470316436.

(12) Gallagher, N. B.; Lee, S. J.; Meinhart, C. D.; Moskovits, M.; Piorek, B. D.; Wise, B. M. Multivariate Curve Resolution Applied to SERS Measurements of 2,4-DNT; poster presented at APACT'13, Chester, Cheshire, U.K., April 24-26, 2013. DOI: 10.13140/RG.2.2.23274.38089.

(13) Forland, B. M.; Gallagher, N. B.; Johnson, T. J. Automated Atmospheric Compensation and Target Detection Algorithm for Standoff Detection. In Algorithms, Technologies, and Applications for Multispectral and Hyperspectral Imaging XXXI; Velez-Reyes, M., Messinger, D. W., Eds.; Proc. SPIE 13455; SPIE: Bellingham, WA, 2025; Art. No. 134550C; SPIE Defense + Commercial Sensing, Orlando, FL, April 13-17, 2025. DOI: 10.1117/12.3054470.

(14) Schuster, A. XXII. The Influence of Radiation on the Transmission of Heat. Philos. Mag. 1903, 5, 243-257. DOI: 10.1080/14786440309462919.

(15) Schuster, A. Radiation Through a Foggy Atmosphere. Astrophys. J. 1905, 21, 1-22. DOI: 10.1086/141186.

(16) Kubelka, P.; Munk, F. Ein Beitrag zur Optik der Farbanstriche. Z. Tech. Phys. 1931, 12, 593-601.

(17) Dahm, D. J.; Dahm, K. D. Interpreting Diffuse Reflectance and Transmittance: A Theoretical Introduction to Absorption Spectroscopy of Scattering Materials; IM Publications: Chichester, West Sussex, U.K., 2007.

(18) Dahm, K. D.; Dahm, D. J. Principles of Diffuse Reflectance Spectroscopy. In Handbook of Near-Infrared Analysis, 4th ed.; Ciurczak, E. W., Igne, B., Workman, J., Jr., Burns, D. A., Eds.; CRC Press: Boca Raton, FL, 2021; Chapter 2, pp 17-40.

(19) Dahm, D. J.; Dahm, K. D. Letter: Math Pretreatment of NIR Reflectance Data: Log(1/R) vs F(R). J. Near Infrared Spectrosc. 1995, 3, 53-59. https://opg.optica.org/jnirs/abstract.cfm?uri=jnirs-3-1-53 (accessed 2026-07-24).

(20) Dahm, D. J.; Dahm, K. D. Representative Layer Theory for Diffuse Reflectance. Appl. Spectrosc. 1999, 53 (6), 647-654. DOI: 10.1366/0003702991947298.

(21) Dahm, D. J.; Dahm, K. D. Bridging the Continuum-Discontinuum Gap in the Theory of Diffuse Reflectance. J. Near Infrared Spectrosc. 1999, 7, 47-53.

(22) Dahm, D. J.; Dahm, K. D.; Norris, K. H. Test of the Representative Layer Theory of Diffuse Reflectance Using Plane Parallel Samples. J. Near Infrared Spectrosc. 2000, 8 (3), 171-181. DOI: 10.1255/jnirs.276.

(23) Dahm, D. J.; Dahm, K. D. The Physics of Near-Infrared Scattering. In Near-Infrared Technology in the Agricultural and Food Industries; Williams, P., Norris, K. H., Eds.; AACC International: St. Paul, MN, 2001; pp 1-17.

(24) Dahm, D. J.; Dahm, K. D.; Norris, K. H. J. Near Infrared Spectrosc. 2002, 10, 53-59.

(25) Dahm, D. J.; Dahm, K. D. Discontinuum Theory of Diffuse Reflection. In Handbook of Vibrational Spectroscopy; Chalmers, J. M., Griffiths, P. R., Eds.; John Wiley and Sons: Chichester, U.K., 2002; Vol. 2, pp 1149-1153.

(26) Dahm, D. J. NIR News 2004, 15 (4), 6-10.

(27) Dahm, D. J. NIR News 2006, 17 (1), 4-7; and 17 (2), 12-15.

(28) Geladi, P.; MacDougall, D.; Martens, H. Linearization and Scatter-Correction for Near-Infrared Reflectance Spectra of Meat. Appl. Spectrosc. 1985, 39 (3), 491-500. DOI: 10.1366/0003702854248656.

(29) Barnes, R. J.; Dhanoa, M. S.; Lister, S. J. Standard Normal Variate Transformation and De-Trending of Near-Infrared Diffuse Reflectance Spectra. Appl. Spectrosc. 1989, 43 (5), 772-777. DOI: 10.1366/0003702894202201.

(30) Golub, G. H.; Van Loan, C. F. Matrix Computations, 4th ed.; Johns Hopkins Studies in the Mathematical Sciences; Johns Hopkins University Press: Baltimore, MD, 2013.

(31) Petersen, K. B.; Pedersen, M. S. The Matrix Cookbook, Nov 15, 2012 version. https://www.math.uwaterloo.ca/~hwolkowi/matrixcookbook.pdf (accessed 2026-07-24).

Appendix. Derivation of the Restricted (Bordered) GLS Estimator

Here is equation 10 through equation 11 derived in full, following the same block-elimination style Gallagher and Watson use to pass from their extended-mixture-model estimator, their equation 21, to its EPO-filtered equivalent, their equation 22, using the standard block-matrix, or Schur-complement, inversion identities^{30,\ 31} and matching the classical restricted-least-squares result.¹¹

Let M = ÃᵀÃ (N × N, assumed invertible), let v = ΔÃ (N × 1) denote the whitened constraint direction, and let r = ÃᵀC denote the ordinary cross-product vector, where C keeps its usual meaning as the reference constituent-value vector. The bordered system of equation 9 is, written out by block,

Mg + vλ = r (row 1) [A.1]

v^Tg = 0 (row 2) [A.2]

From row 1, since M is invertible,

g = M⁻¹r − λM⁻¹v [A.3]

Substituting equation A.3 into row 2, equation A.2:

v^T(M⁻¹r − λM⁻¹v) = 0 [A.4]

v^TM⁻¹r − λ(v^TM⁻¹v) = 0 [A.5]

and, since vᵀM⁻¹v is a scalar, solving for λ gives us equation 11:

λ = (v^TM⁻¹v)⁻¹ v^TM⁻¹r [A.6]

Substituting equation A.6 back into equation A.3 gives equation 10:

g = M⁻¹r − M⁻¹v(v^TM⁻¹v)⁻¹ v^TM⁻¹r = [I − M⁻¹v(v^TM⁻¹v)⁻¹ v^T] M⁻¹r [A.7]

The bracketed term in equation A.7 is itself an oblique projector that annihilates the component of M⁻¹r along v; left-multiplying equation A.7 by vᵀ and substituting M, v, and r confirms vᵀg = 0 identically, as required. Reinstating M = ÃᵀÃ, v = ΔÃ, and r = ÃᵀC recovers equation 10 exactly.

The multi-constraint system of equation 13 follows from the identical elimination, with v replaced by the N × L matrix P̃ and λ by an L × 1 vector. Row 2 becomes P̃ᵀg = 0, an L-equation system, and the same substitution gives

g = [I − M⁻¹P̃(P̃^TM⁻¹P̃)⁻¹ P̃^T] M⁻¹r [A.8]

which is the direct L-dimensional analogue of equation A.7. It is structurally similar to, though not identical with, Gallagher and Watson's EPO filter, their equation 22, shown here as equation 16: both are oblique projectors that annihilate P̃, and despite their different algebraic forms — equation A.8 explicitly weights by M⁻¹ = (ÃᵀÃ)⁻¹, while the EPO or ELS filter of equation 16 does not — the two give the identical fitted coefficient g: any g satisfying P̃ᵀg = 0 can be written as g = WELS⁻¹h for some h, reducing the constrained minimization of ‖C − Ãg‖² to the unconstrained fit equation 16 already performs; by uniqueness of the restricted-least-squares estimator (cf. ref. 11), that fit must agree with equation A.8. Equation A.8 and equation 16 are accordingly two computational routes to the same hard, exact constraint — imposing gᵀP̃ = 0 on the fitted coefficients — rather than a hard-versus-soft distinction; the genuinely soft alternative is the augmented-design estimator of equations 14-15, which does not satisfy gᵀP̃ = 0 exactly.

About the Authors

Howard Mark and Jerome Workman, Jr. are internationally recognized as authorities in spectroscopy, chemometrics, and analytical instrumentation. Their collaboration began at Technicon Instrument Corporation and has spanned four decades, producing influential books, technical papers, and the long-running Statistics in Spectroscopy and Chemometrics in Spectroscopy columns. Mark, a physicial chemist, spectroscopist, and chemometrician, has made major contributions to near-infrared spectroscopy, statistical analysis, and instrumentation. Workman, a biochemist and analytical spectroscopist, inventor, and educator, has contributed extensively to spectroscopy, chemometrics, instrument development, and multivariate analysis. Together, they have authored numerous books and hundreds of technical publications and have received many of the field's distinguished professional honors. Their work has helped establish and advance the practical foundations of modern chemometrics and quantitative spectroscopy. Together they received the 2020 NYSAS Gold Medal Award.