Classical Least Squares and Clutter Suppression

Abstract

Classical least squares (CLS) is based on a linear mixture model. It is the underlying model for multivariate curve resolution and many target detection algorithms. CLS is also useful in process modeling and classification. However, there is a misconception that CLS cannot account for interferences. This white paper shows the derivation of the CLS and weighted least squares models and discusses two forms of the model that account for interferences: generalized least squares and extended least squares. Although the interferences need not be known explicitly, the key is that the interference signal must be characterized. This restriction is similar to the requirements for inverse least squares models such as partial least squares. This white paper provides derivations and model concepts for four different forms of CLS: traditional CLS, weighted CLS, generalized least squares (GLS) and extended least squares (ELS). GLS and ELS are used to account for interferences and perform clutter suppression in CLS and inverse least squares models.

Introduction: The Classical Least Squares Model

For a given set of basis vectors S and measurements X, the classical least squares (CLS) model can be used to estimate contributions C. The CLS model is

X = CS^T + E      [1]

where X is a M × N matrix of measurements. For example, each row of X can be represented as a N × 1 column vector x_m (and x_m^T is a 1 × N vector), C is a M × K matrix of contributions for K analytes and S is a N × K matrix where each column corresponds to a pure component spectrum. The units of C depend on the units of S, and the coefficients in C correspond to the contribution of each basis vector to the signal in X. The CLS model is also known as a linear-mixture-model, forward least squares and multicomponent Beer’s law.

In spectroscopy, the columns of S correspond to the spectrum of each analyte contributing to the signal in X, and the least-squares estimator for C is given by

Ĉ = XS( S^TS )^-1      [2]

where Ĉ is the estimate of C, and S( S^TS )^-1 is a Moore-Penrose pseudoinverse. This white paper shows the details of the derivation of [2] from [1]. The first derivation provides an estimator for a single row of X, and the second, which provides an estimator for all rows of X using matrix identities, is shown in Appendix A. However, to provide a more general approach, the derivations are shown for weighted versions of CLS. For a more complete list of identities see¹ and for a detailed discussion of derivations of least squares estimators see.² See Appendix B for a list of definitions.

Derivation of Weighted CLS for a Single Measurement

For the m^th measurement (the m^th row of X as a column vector), the CLS model is

x_m = Sc_m+e_m      [3]

and the weighted objective function to minimize is

O( c_m ) = e_m^TW^-1e_m = ( x_m - Sc_m )^TW^-1( x_m - Sc_m )      [4]

where W^-1 is a N × N weighting matrix discussed further below. Dropping the subscript x_m for clarity and taking the derivative with respect to c^T gives

(dO( c )/dc^T) = 2( x - Sc )^TW^-1(d( x - Sc )/dc^T) = - 2( x - Sc )^TW^-1S      [5]

where it is noted that (d( c )/dc^T) = I. Setting the result to zero, taking the transpose and rearranging gives

S^TW^-1Sc = S^TW^-1x.      [6]

Multiplying both sides on the left by ( S^TW^-1S )^-1 gives the CLS estimator

ĉ = ( S^TW^-1S )^-1S^TW^-1x      [7]

where ^ is used to indicate an estimation.

An estimate of the estimation error is obtained by taking the differential of the estimator:

dĉ = ( S^TW^-1S )^-1S^TW^-1dx      [8]

The outer product of equation 8 is used to approximate the error covariance as

dĉdĉ^T = ( S^TW^-1S )^-1S^TW^-1dxdx^TW^-1S( S^TW^-1S )^-1.      [9]

In the asymptotic limit dĉdĉ^T → V( ĉ ) K × K and dxdx^T → V( x ) and N × N where V( ) is the error covariance. An assumption of the weighted CLS model is that V( x )=W and that W^(-1/2)e is normally distributed with mean zero. Substitution of W =dxdx^T into equation 9 gives

V( ĉ ) = ( S^TW^-1S )^-1S^TW^-1S( S^TW^-1S )^-1 = ( S^TW^-1S )^-1.      [10]

The weighted residuals, q, are given by

q = e^TW^-1e=( e^TW^(-1/2) )( W^(-1/2)e )      [11]

where e = x - Sĉ. The residuals q, and the residual contributions ( W^(-1/2)e ) can be used as diagnostics for goodness of model fit.

For completeness, a derivation for multiple measurements, X, is given in Appendix A.

There are four common CLS model forms that differ based on the definition of W. In each of the four forms it is assumed that q is mean zero and normally distributed, and the next two sections discuss the four models.

Classical Least Squares (CLS) and Weighted Least Squares (WLS)

Suppose that X_d is measured on dark measurements (i.e., noise background) with mean x̅_d then

W_d = [1/(M_d - 1)]( X_d- 1x̅_d^T )^T( X_d- 1x̅_d^T ).      [12]

In this case, the mean background is subtracted from the measurements such that x → x - x̅_d. If the measured noise in X_d is uncorrelated and of similar magnitude for each of the N variables, then W_d can be represented asymptotically as W_CLS=σ_x²I and substitution of W_CLS for W in equation 7 results in equation 2 giving the CLS model. For CLS, σ_x² is estimated from

σ_x² = [1/(N - K)]e^Te.      [13]

The CLS model is often used for laboratory and process measurements where all the basis vectors in S are known.³

Next, assume the measured noise in X_d is uncorrelated and not of similar magnitude for each of the N variables. In this case, W_d can be represented asymptotically by a diagonal matrix W_WLS where the diagonal elements are given by the variance for each measurement σ_x,n² > 0. Substitution of W_WLS for W in equation 7 gives the WLS model. The variance σ_x,n² is estimated from

σ_x,n² = [1/(M_d - 1)]ε_n^Tε_n      [14]

where ε_n is a N × 1 column of X_d- 1x̅_d^T corresponding to the mean-centered dark current on the n^th channel. However, it should be noted that the diagonal elements of W_WLS can be defined differently depending on model requirements. In which case, the residuals e^TW_WLS^-1e are interpreted accordingly. The WLS model is often used in spectroscopic applications where all the basis vectors in S are known.⁴

CLS and WLS are useful models, but what happens when, after subtracting the mean background signal, W_d is not diagonal? I.e., what happens if interference signals are present in the background that are not accounted for? In this case, the weighted residuals are not normally distributed and the estimates can be biased. To account for interferences requires the definition of “clutter” and two models designed to provide “clutter suppression” are discussed in the next section.

Generalized Least Squares (GLS) and Extended Least Squares (ELS)

The basis vectors in S correspond to signal in X that is of interest – it is target signal. In contrast, interference and noise signal in X corresponds to signal not of interest. This interference signal is defined as “clutter” and generalized least squares (GLS) is a weighted version of CLS that accounts for the interferences i.e., GLS provides clutter suppression. The weighting in GLS, W_GLS, is a non-diagonal symmetric matrix that is used for clutter suppression. As will be shown, the interferences need not be known explicitly but the interference signal must be characterized. This restriction parallels the requirements for inverse least squares models (ILS) such as partial least squares. For ILS models it is often stated that interferences must vary in the calibration for the model to be robust to the interference. Clutter is characterized by measurements that include clutter signal but do not contain target signal, X_c. For example, clutter can be characterized from a pooled estimated of intra-class variance where target signal does not vary within each class, but the interferences do vary. An example of a class might be a set of replicate measurements where the same sample has been reintroduced to the measurement system multiple times. GLS is also called the Aitken estimator⁵ or matched filter;⁶ it is discussed in detail.⁷

In the additive model,⁸ the measured signal is modeled as the sum of target signal and the mean clutter signal, x̅_c, where

x = Sc +x̅_c+e      [15]

and the clutter covariance is given by

W_GLS = W_c = [1/(M_c - 1)]( X_c- 1x̅_c^T )^T( X_c- 1x̅_c^T ).      [16]

The estimator is given as

ĉ = ( S^TW_GLS^-1S )^-1S^TW_GLS^-1( x -x̅_c )      [17]

If the clutter is not additive, a non-centered version calculates the clutter covariance as

W_GLS = W_c = (1/M_c)X_c^TX_c      [18]

with corresponding estimator given as

ĉ = ( S^TW_GLS^-1S )^-1S^TW_GLS^-1x.      [19]

The non-centered model, and the replacement model discussed below, are appropriate for reflectance spectra.⁷

A difficulty with the GLS model is that W_c is often ill-conditioned or even not full rank and regularization is required when calculating the inverse. One consideration on how to regularize assumes that as the interference signal → 0, W_GLS→σ_x²I (or W_WLS) and a simple approach is based on ridge regression. In this case, W_GLS is replaced by W_GLS→W_GLS+αI where α ≈ σ_x². Other approaches are also used.⁷

Another useful model used to account for clutter uses extended mixture model⁹ also known as a replacement model⁸ given by

x = [S P][c t]+e      [20]

where P are basis vectors of X_c and t are coefficients to be estimated. The columns of P can consist of loadings from principal components analysis^10,11 or any set of basis vectors that span the systematic variance in X_c. The extended least squares (ELS) estimator is given by

[ĉ t̂] = ( [S P]^T[S P] )^-1[S P]^Tx.      [21]

ELS is directly related to orthogonal subspace projection,¹² orthogonal background subtraction,¹³ and external parameter orthogonalization (EPO).¹⁴ It can be shown with block inverse identities that the estimator for ĉ can be given by

ĉ = ( S^TW_ELS^-1S )^-1S^TW_ELS^-1x      [22]

where W_ELS^-1 = I - P( P^TP )^-1P^T is an EPO filter. An advantage of equation 21 is that it provides estimates of t̂ that may be useful in exploring the clutter signal.

Discussion and Conclusions

The estimator for weighted CLS models is given generally by equation 7 and can also be written as

ĉ = ( S̃^TS̃ )^-1S̃^Tx̃.      [23]

where S̃ = W^(-1/2)S and x̃ = W^(-1/2)x. For GLS and ELS these can be considered filtered target and measurements, and provide a motivation for using clutter suppression in inverse least squares modeling such as partial least squares.¹⁴ Equation 23 can be interpreted as performing CLS on filtered data. When GLS is used, it is often referred to as “whitening” or “weighting”. When ELS is used, it is referred to as orthogonalization.

The estimation error is given by equation 10 rewritten as

V( ĉ ) = ( S^TW^-1S )^-1 = ( S̃^TS̃ )^-1,      [24]

with residuals q from equation 11:

q=( e^TW^(-1/2) )( W^(-1/2)e ) = ẽ^Tẽ.      [25]

The residuals, q, and residual contributions ẽ are used as diagnostics for goodness of model fit. High residuals are an indication that S does not span the space of the measured x, and/or the clutter was not appropriately identified. This has implications for target detection⁷ and stepwise regression modeling^15,16 because estimates with high residuals should not be trusted as they may be biased.

The text above showed four different types of CLS models depending on the definition of weighting W. The weighting to use depends on the sensing scenario of the measurements and can be combined to provide the desired estimator. For example, if a measurement system is expected to always include a known interferent P, it can be used to develop an EPO filter W_ELS^-1. Additionally, if that same process also includes unknown interference that can be characterized a separate GLS weighting W_GLS can also be obtained and the estimator is given by

ĉ = ( S^TW_GLS^(-1/2)W_ELS^-1W_GLS^(-1/2)S )^-1S^TW_GLS^(-1/2)W_ELS^-1W_GLS^(-1/2)x.      [26]

An advantage of this approach is that while W_ELS^-1 may remain constant W_GLS can be easily updated as new measurements are obtained.¹⁷ This is the reason that using GLS filtering is often referred to as an “adaptive” approach.

The CLS-based modeling approach has two distinct advantages. The first is that S is interpretable and can be easily modified. The second is that W_GLS and W_ELS are easily updated with new data. For example, in hyperspectral imaging they can be determined on an image specific basis image-to-image.¹⁸

For ILS models, the definition x̃ = W^(-1/2)x used in equation 23 can also be used for clutter suppression. For example, an ILS model is defined as

y=Xb + e      [27]

where y are the reference values and b is a regression vector. Preprocessing the predictor variables X with the inverse square root of the clutter covariance provides clutter suppression in ILS:

y=XW^(-1/2)b + e =X̃b + e      [28]

where b and e are generally different in equation 27 than compared to those in equation 28.¹⁹

References

1. Petersen, K. B.; Pedersen, M. S. The Matrix Cookbook; November 15, 2012. Available at: http://www.math.uwaterloo.ca/~hwolkowi/matrixcookbook.pdf.
2. Magnus, J. R.; Neudecker, H. Matrix Differential Calculus with Applications in Statistics and Econometrics; Wiley: New York, 1999.
3. Ayesa, U.; Li, S.; Dance, Z. E. X. Harnessing the Power of Chemometric Modeling and Plate Readers for High-Throughput Analysis of Biocatalytic Reactions. Anal. Chem. 2025, 97 (22), 11760–11766. DOI: https://doi.org/10.1021/acs.analchem.5c01207.
4. Chu, P. M.; Guenther, F. R.; Rhoderick, G. C.; Lafferty, W. J. The NIST Quantitative Infrared Database. J. Res. Natl. Inst. Stand. Technol. 1999, 104 (1), 59–81. DOI: https://doi.org/10.6028/jres.104.004.
5. Aitken, C. On Least-Squares and Linear Combinations of Observations. Proc. R. Soc. Edinburgh 1936, 55, 42–48. DOI: https://doi.org/10.1017/S0370164600014346.
6. Turin, G. L. An Introduction to Matched Filters. IEEE Trans. Inf. Theory 1960, 6 (3), 311–329. DOI: https://doi.org/10.1109/TIT.1960.1057571.
7. Gallagher, N. B.; Goyetche, R.; Rubio, J. M. A.; Kucheryavskiy, S. Extended Least Squares (ELS) and Generalized Least Squares (GLS) for Clutter Suppression in Hyperspectral Imaging: A Theoretical Discussion. Chemom. Intell. Lab. Syst. 2024, 244, 105032. DOI: https://doi.org/10.1016/j.chemolab.2023.105032.
8. Theiler, J.; Ziemann, A.; Matteoli, S.; Diani, M. Spectral Variability of Remotely Sensed Target Materials: Causes, Models, and Strategies for Mitigation and Robust Exploitation. IEEE Geosci. Remote Sens. Mag. 2019, 7 (2), 8–30. DOI: https://doi.org/10.1109/MGRS.2019.2890997.
9. Martens, H.; Næs, T. Multivariate Calibration, 2nd ed.; Wiley: Chichester, U.K., 1989.
10. Bro, R.; Smilde, A. K. Principal Components Analysis. Anal. Methods 2014, 6, 2812–2831. DOI: https://doi.org/10.1039/C3AY41907J.
11. Jackson, J. E. A User's Guide to Principal Components; Wiley: New York, 1991.
12. Harsanyi, J. C.; Chang, C.-I. Hyperspectral Image Classification and Dimensionality Reduction: An Orthogonal Subspace Projection Approach. IEEE Trans. Geosci. Remote Sens. 1994, 32 (4), 779–785. DOI: https://doi.org/10.1109/36.298007.
13. Hayden, A.; Niple, E.; Boyce, B. Determination of Trace-Gas Amounts in Plumes by the Use of Orthogonal Digital Filtering of Thermal-Emission Spectra. Appl. Opt. 1996, 35 (16), 2802–2809. DOI: https://doi.org/10.1364/AO.35.002802.
14. 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: https://doi.org/10.1016/S0169-7439(03)00051-0.
15. Thompson, C.; Gallagher, N. B.; Hughey, K. D.; Dunlap, M. K.; Meyers, T. L.; Johnson, T. J. An Interactive Spectral Analysis Tool for Chemical Identification and Quantification of Gas-Phase Species. Appl. Spectrosc. 2023, 77 (6), 557–568. DOI: https://doi.org/10.1177/00037028231169304.
16. Gallagher, N. B. Classical Least Squares for Detection and Classification. In Hyperspectral Imaging; Amigo, J. M., Ed.; Data Handling in Science and Technology; Elsevier: Amsterdam, 2019; Vol. 32, pp 231–246. DOI: https://doi.org/10.1016/B978-0-444-63977-6.00011-0.
17. 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: https://doi.org/10.13140/RG.2.2.23274.38089.
18. Forland, B. M.; Gallagher, N. B.; Johnson, T. J. Application of Automated Iterative Target Detection for Standoff Hyperspectral Imaging. J. Appl. Remote Sens. 2025, 19 (1), 016509. DOI: https://doi.org/10.1117/1.JRS.19.016509.
19. Martens, H.; Høy, M.; Wise, B. M.; Bro, R.; Brockhoff, P. B. Pre-Whitening of Data by Covariance-Weighted Pre-Processing. J. Chemom. 2003, 17 (3), 153–165. DOI: https://doi.org/10.1002/cem.780.

Appendix A. Derivation of the CLS Estimator for Multiple Measurements

For completeness, the derivation for multiple measurements is provided and for clarity the non-weighted version is shown. For a set of measurements, X, the CLS model is

X = CS^T + E      [A.1]

and the objective is to minimize

O( C ) = tr( E^TE ) = vec( E )^Tvec( E ) = vec( X - CS^T )^Tvec( X - CS^T )

= [ vec( X ) - vec( CS^T ) ]^T[ vec( X ) - vec( CS^T ) ]

= [ vec( X ) - ( S⊗I )vec( C ) ]^T[ vec( X ) - ( S⊗I )vec( C ) ].      [A.2]

Taking the derivative of A.2 with respect to vec( C )^T gives

(dO( C )/dvec( C )^T) = 2[ vec( X ) - ( S⊗I )vec( C ) ]^T(d[ vec( X ) - ( S⊗I )vec( C ) ]/dvec( C )^T)

= - 2[ vec( X ) - ( S⊗I )vec( C ) ]^T( S⊗I )      [A.3]

where (d[ vec( C ) ]/dvec( C )^T) = I and using the identity vec( ABC ) = ( C^T⊗ A )vec( B ). Setting the result in A.3 to zero, noting the identity ( A⊗ B )^T( C⊗ D ) = A^TC⊗B^TD and rearranging gives

0 = - 2[ vec( X ) - ( S⊗I )vec( C ) ]^T( S⊗I ),

vec( X )^T( S⊗I ) = vec( C )^T( S⊗I )^T( S⊗I ),

vec( X )^T( S⊗I ) = vec( C )^T( S^TS⊗I ).

Taking the transpose and swapping the equality yields

( S^TS⊗I )vec( C ) = ( S^T⊗I )vec( X ).      [A.4]

Using the identity vec( ABC ) = ( C^T⊗ A )vec( B ) A.4 can be written as

C( S^TS ) = XS.      [A.5]

Multiplying both sides on the right by ( S^TS )^-1 gives the CLS estimator

Ĉ = XS( S^TS )^-1.      [A.6]

where ^ is used to indicate an estimation.

The residuals, q, are given by

q = diag( EE^T )      [A.7]

where E = X - ĈS^T and diag( ) takes the diagonal elements.

Appendix B. Variable and Operator Definitions

Variable Definitions Table

Matrices and Vectors

Scalars and Parameters

Functions and Operations

About the Authors

Neal B. Gallagher, Ph.D.
Vice President and co-founder of Eigenvector Research, Inc. established January 1995.
www.eigenvector.com
Neal Gallagher earned a Ph.D. in Chemical Engineering with a mathematics minor from the University of Arizona in 1992, an M.S. in Chemical Engineering from the University of Washington in 1987 and B.S. degrees in Chemical Engineering and Engineering Physics from the University of Colorado in 1985.
Recent research includes novel algorithms for iterative target detection in hyperspectral imaging and shift-invariant tri-linearity (SIT) for modeling hyphenated chromatography data. Since founding Eigenvector in January 1995, Neal has been intimately involved in chemometrics consulting, teaching short courses and software development including algorithms for detection, classification and quantification. Specific interests include hyperspectral imaging, process modeling, multi-variate curve resolution and classical least squares modeling.

Nathanial E. Watson, PhD.

Principal Data Scientist and Software Developer, Eigenvector Research, Inc.

www.eigenvector.com

Nate Watson holds BS, MS and PhD degrees in analytical chemistry and other applied disciplines earned between 1999 and 2017. He has nearly three decades of experience using chemometric methods to answer analytical questions.

His research interests are varied. He particularly enjoys bridging complicated scientific and mathematical concepts into language that is easily understood by end users and technicians.

He is married and has four children. He currently resides in Western Montana.