Multivariate Curve Resolution - Alternating Least Squares

Multivariate Curve Resolution is a family of methods intended to recover scientifically meaningful component contributions from multivariate mixture data through an additive bilinear, or later extended multilinear, representation. de Juan, Jaumot, and Tauler define MCR as methods that provide a chemically meaningful additive bilinear model of pure contributions from a mixed data matrix. MCR is the methodology. It is not itself one algorithm.

MCR-ALS is the combination of that bilinear MCR model with an Alternating Least Squares optimizer. Jaumot, Gargallo, de Juan, and Tauler describe ALS as an iterative calculation of the contribution matrix $\mathbf{C}$ and the response matrix ${\mathbf{S}}^{\top }$ that fit the experimental matrix $\mathbf{D}$ for a proposed number of components, starting from an initial estimate of either factor matrix, with optional constraints applied during the optimization.

This page teaches classical two-way MCR-ALS. Multiway, weighted, nonlinear, kinetic hard-modeling, and correlation constraints are later extensions and are not part of the main derivation.

Lawton and Sylvestre formulated self-modeling curve resolution as the recovery of two overlapping functions from observed additive mixtures of those functions, a problem they associated with spectrophotometry, chromatography, and related fields. Under mild restrictions the feasible bands of functions can collapse to unique curves; without those restrictions a range of feasible profiles remains. That historical setting is the mixture-analysis problem. It is not the modern MCR-ALS algorithm.

In the spectroscopic bilinear setting used throughout this page, rows of $\mathbf{D}$ are mixture measurements and the goal is to recover contribution-like profiles and spectral-like responses whose product reconstructs the mixed data. The algorithm does not assign molecular names. Identity requires reference spectra, chemical knowledge, or other analytical evidence.

Jaumot et al. write the MCR bilinear model as

with $\mathbf{D}\in {\mathbb{R}}^{I\times J}$, $\mathbf{C}\in {\mathbb{R}}^{I\times N}$, ${\mathbf{S}}^{\top }\in {\mathbb{R}}^{N\times J}$, and $\mathbf{E}\in {\mathbb{R}}^{I\times J}$. Tauler used the equivalent storage $\mathbf{S}\in {\mathbb{R}}^{J\times N}$, so that ${\mathbf{S}}^{\top }$ has size $N\times J$. $N$ is the selected number of resolved components. $I$ is the number of observations (samples, times, or process conditions). $J$ is the number of variables (wavelengths or other channels).

For spectroscopic measurements, Jaumot et al. associate rows of $\mathbf{D}$ with measured spectra, columns of $\mathbf{C}$ with concentration-like profiles, and rows of ${\mathbf{S}}^{\top }$ with resolved spectral profiles. Those $\mathbf{C}$ values are not automatically absolute concentrations. Absolute scale requires additional information such as known spectra, calibration, or a mass-balance constraint. Equation (1) is the multiwavelength matrix form of a Lambert-Beer-type bilinear law when that chemical model is appropriate.

The residual matrix is $\mathbf{E}=\mathbf{D}-\mathbf{C}{\mathbf{S}}^{\top }$. Jaumot et al. state that $\mathbf{E}$ should ideally approach experimental error. It is the structure not reproduced by the selected bilinear model. It is not claimed to equal experimental noise exactly.

Conventional MCR-ALS requires an initial estimate of $\mathbf{C}$ or of ${\mathbf{S}}^{\top }$ together with a chosen $N$. Jaumot et al. mention Evolving Factor Analysis and SIMPLISMA-derived methods as possible sources of those estimates. They are not the only methods. de Juan, Jaumot, and Tauler add that known pure spectra, when available, can also be used.

$N$ is not discovered automatically by ALS. Rank analysis, SVD or PCA, local-rank methods such as EFA, and chemical knowledge can all inform the choice. There is no universal 95% variance rule. Different starting estimates can produce different resolved profiles when rotational ambiguity or local optimization effects are present. Marcel Maeder, reviewing the 2005 toolbox, noted that the selection of starting values is crucial in many instances.

SIMPLISMA and EFA are initialization or local-rank methods with their own literature. They are named here only as possible initial-estimate sources. Their equations are not derived on this page. The educational code accepts a deterministic initial ${\mathbf{S}}^{\top }$.

For unconstrained least-squares fitting of the bilinear model, ALS targets

where the Frobenius norm is $\Vert \mathbf{A}{\Vert }_F=\sqrt{{\sum }_{i,j}{a}_{ij}^2}$. Constrained MCR-ALS is not presented as a jointly convex global optimization problem. ALS alternates two linear least-squares subproblems. That alternation is not a one-step formula for a unique chemical resolution.

With ${\mathbf{S}}^{\top }$ held fixed, the contribution update is

subject to any selected constraints. In the unconstrained case this is the least-squares solution of $\mathbf{C}{\mathbf{S}}^{\top }\approx \mathbf{D}$. Tauler wrote the equivalent pseudoinverse form $\mathbf{C}=\mathbf{D}({\mathbf{S}}^{\top }{)}^{+}$. If ${\mathbf{S}}^{\top }\in {\mathbb{R}}^{N\times J}$, then $({\mathbf{S}}^{\top }{)}^{+}\in {\mathbb{R}}^{J\times N}$ and $\mathbf{D}({\mathbf{S}}^{\top }{)}^{+}$ has size $I\times N$, matching $\mathbf{C}$. Educational Python uses numpy.linalg.lstsq on the transposed system rather than an explicit inverse.

With $\mathbf{C}$ held fixed, the spectral update is

Unconstrained, Tauler's corresponding form is ${\mathbf{S}}^{\top }={\mathbf{C}}^{+}\mathbf{D}$. If $\mathbf{C}\in {\mathbb{R}}^{I\times N}$, then ${\mathbf{C}}^{+}\in {\mathbb{R}}^{N\times I}$ and ${\mathbf{C}}^{+}\mathbf{D}$ has size $N\times J$, matching ${\mathbf{S}}^{\top }$. Some implementations, including Tauler 1995, optionally replace $\mathbf{D}$ by a PCA-reproduced matrix of rank $N$. Jaumot et al. report lack of fit against both the raw matrix and a PCA reproduction. The educational code on this page operates on the input matrix $\mathbf{D}$.

Constraints introduce chemically or mathematically justified information into the resolution. They can reduce the set of feasible solutions. They are optional. Incorrect or overly strict constraints can distort the resolution. de Juan, Jaumot, and Tauler recommend comparing models with and without a doubtful constraint: a genuine constraint should not produce a large drop in explained variance.

Non-negativity forces selected profiles to contain nonnegative values. It is often reasonable for concentration-like profiles and for absorbance-type spectra. It is not required for every spectroscopy. Derivative spectra and difference signals can contain negative values. Jaumot et al. implement non-negativity by nonnegative least squares (Lawson and Hanson), by fast nonnegative least squares (Bro and de Jong), or by replacing negatives with zeros. Those three operations are not interchangeable. The educational code uses true nonnegative least squares in the regression step. It does not clip an unconstrained solution.

Unimodality allows a single maximum per constrained profile under the selected implementation. It can be useful for some elution or monotonic reaction profiles. Chromatographic peaks are not always unimodal in practice. This page does not invent a universal unimodality algorithm, and the educational code does not implement it.

Closure expresses a mass-balance relation when that relation is scientifically appropriate. Jaumot et al. allow a constant or a known varying total, and equality or inequality forms. Closure does not mean that every row of $\mathbf{C}$ always sums to one. That is one possible normalization convention, not the definition. Selectivity and local-rank constraints encode known absence or unique presence of components in selected windows. Equality constraints can fix known spectra or known contribution values. Hard kinetic models and correlation (calibration) constraints are later extensions and are not derived here.

Classical two-way MCR-ALS is a bilinear factorization, two alternating least-squares updates, optional constraints, and source-defined fit statistics. It is not a closed-form unique chemical resolution.

Equations (1) to (3) define ordinary two-way ALS. Equation (4) is the Jaumot et al. lack-of-fit convention against the input matrix. Equation (5) is rotational ambiguity. sklearn-style clipping of negatives is not part of this mathematics.

The procedure follows the Jaumot et al. ALS cycle: initialize one factor matrix, alternate least-squares updates of $\mathbf{C}$ and ${\mathbf{S}}^{\top }$, reconstruct, and stop on a residual-based relative change or an iteration limit. Unimodality, closure, selectivity, and profile normalization are omitted from the educational loop. When non-negativity is selected, it is incorporated in the least-squares solve, not applied afterwards by setting negatives to zero.

Jaumot et al. declare convergence when, in two consecutive iterative cycles, the relative difference in the standard deviation of residuals between experimental and ALS-calculated values is smaller than a user-selected threshold. They note 0.1% as a commonly chosen value in their implementation, and they note that users may change it. 0.1% is not a universal MCR-ALS constant. The educational code uses that residual-SD relative change, with 0.1% as an example threshold. de Juan, Jaumot, and Tauler also describe stopping on a preset iteration count or on a threshold for fit improvement between consecutive iterations. Those stopping rules are related but are not claimed to be identical.

Jaumot et al. define percentage lack of fit as

where $d_{ij}$ is an entry of the input matrix and $e_{ij}$ is the corresponding residual. Their variance-explained quantity is

That is a fit statistic under their MCR-ALS convention. It is not ordinary centered regression $R^2$. Their residual standard deviation is

with denominator equal to the number of entries of $\mathbf{D}$. No extra degrees-of-freedom correction is introduced.

A low lack of fit, or a high explained variance, indicates that $\mathbf{C}{\mathbf{S}}^{\top }$ reconstructs $\mathbf{D}$ well. It does not by itself prove that the resolved profiles are chemically pure, that component identities are correct, or that the solution is unique.

The bilinear factorization is generally non-unique. de Juan, Jaumot, and Tauler distinguish permutation ambiguity, intensity (scale) ambiguity, and rotational ambiguity. Component order is arbitrary: swapping corresponding columns of $\mathbf{C}$ and rows of ${\mathbf{S}}^{\top }$ leaves the product unchanged. Intensity ambiguity means that multiplying one profile by a constant and dividing the paired profile by the same constant leaves the product unchanged, which is why recovered $\mathbf{C}$ values are in arbitrary units unless reference intensity information is used. Rotational ambiguity is

for an invertible $N\times N$ matrix $\mathbf{T}$ whose transformed profiles still obey the selected constraints. Tauler, Smilde, and Kowalski focused on rotational and intensity ambiguities and on the role of selectivity and local rank in uniqueness. ALS alone does not guarantee unique pure profiles. Constraints, selectivity, and multiset information can reduce the feasible set. They do not automatically produce uniqueness.

Jaumot et al. describe simultaneous analysis of several experiments that share spectral responses as a column-wise augmented bilinear model

with a common ${\mathbf{S}}^{\top }$ and experiment-specific contribution blocks. They also describe row-wise augmentation when several techniques share concentration profiles. Simultaneous analysis can add information and may reduce some ambiguities or rank-deficiency issues. Jaumot et al. write that it may eventually remove rotational ambiguities and rank-deficiency problems in favorable structures. That is not a guarantee of uniqueness for every multiset. Trilinearity and other multiway constraints are later extensions and are not derived here. The educational code remains two-way.

Jaumot et al. list chemical reactions, industrial processes, chromatographic elutions, spectroscopic images, and environmental data among bilinear systems analyzed by MCR-ALS, with spectroscopic, electrochemical, and composition-profile responses. de Juan, Jaumot, and Tauler add process monitoring, HPLC-DAD, hyperspectral imaging, and related bilinear measurements. Those are sourced application classes, not a claim that MCR-ALS is the default method for every mixture data set.

PCA finds orthogonal latent directions ordered by variance. MCR-ALS seeks bilinear component profiles that can carry chemical constraints. MCR-ALS factors are not required to be orthogonal. MCR-ALS is not PCA with constraints. SVD supplies an orthogonal matrix factorization and can help estimate rank; SVD components are not pure spectra. PLS is a supervised latent-variable regression relating predictor and response blocks. MCR-ALS resolves bilinear mixture contributions from the mixture matrix and selected constraints. It is not a substitute for calibration by PLS or MLR.

Camp describes pyMCR as an MCR alternating-regression library in which ordinary or nonnegative least squares is one possible regressor choice. That software is related to MCR-ALS but is not claimed to be byte-for-byte equivalent to the Barcelona MCR-ALS implementation. The educational listings below implement the two-way ALS updates taught on this page.

The first Python listing implements two-way ALS from a supplied initial ${\mathbf{S}}^{\top }$: unconstrained numpy.linalg.lstsq updates, optional scipy.optimize.nnls when non-negativity is selected, reconstruction, Jaumot lack of fit, and stopping on relative change of residual standard deviation. The example threshold tol_percent=0.1 follows the commonly selected value in Jaumot et al. It is not a universal default. Unimodality, closure, selectivity, and SIMPLISMA are not implemented. The second listing only calls that educational function with non-negativity. Outputs are text: profiles, residuals, iteration count, and fit statistics.

MATLAB uses transpose-backslash least squares, which matches the same unconstrained updates, and official lsqnonneg for nonnegative solves. The Barcelona MCR-ALS GUI (historically als2004, later mcr_main) is a specialized MATLAB implementation from that research group. It is not an official MathWorks toolbox and is not reproduced here.

import numpy as npfrom scipy.optimize import nnls  def _validate_matrix(A, name):    A = np.asarray(A, dtype=float)    if A.ndim != 2:        raise ValueError(f"{name} must be a 2D array.")    if min(A.shape) == 0:        raise ValueError(f"{name} must have at least one row and one column.")    if not np.all(np.isfinite(A)):        raise ValueError(f"{name} must contain only finite values.")    return A  def frobenius_sse(D, C, ST):    D = _validate_matrix(D, "D")    C = _validate_matrix(C, "C")    ST = _validate_matrix(ST, "ST")    if C.shape[0] != D.shape[0] or ST.shape[1] != D.shape[1]:        raise ValueError("C and ST must match the dimensions of D.")    if C.shape[1] != ST.shape[0]:        raise ValueError("C and ST must share the component dimension N.")    return float(np.sum((D - C @ ST) ** 2))  def lack_of_fit_percent(D, E):    D = _validate_matrix(D, "D")    E = np.asarray(E, dtype=float)    if E.shape != D.shape:        raise ValueError("E must have the same shape as D.")    den = float(np.sum(D ** 2))    if den == 0:        raise ValueError("The data matrix D must have nonzero sum of squares.")    return 100.0 * float(np.sqrt(np.sum(E ** 2) / den))  def variance_explained(D, E):    D = _validate_matrix(D, "D")    E = np.asarray(E, dtype=float)    if E.shape != D.shape:        raise ValueError("E must have the same shape as D.")    den = float(np.sum(D ** 2))    if den == 0:        raise ValueError("The data matrix D must have nonzero sum of squares.")    return float((den - np.sum(E ** 2)) / den)  def residual_sd(E):    E = np.asarray(E, dtype=float)    if E.size == 0:        raise ValueError("E must be nonempty.")    return float(np.sqrt(np.sum(E ** 2) / E.size))  def update_C(D, ST, nonneg=False):    D = _validate_matrix(D, "D")    ST = _validate_matrix(ST, "ST")    if ST.shape[1] != D.shape[1]:        raise ValueError("ST must have shape (N, J) matching D.")    if ST.shape[0] < 1:        raise ValueError("The number of components N must be at least 1.")    if nonneg:        C = np.empty((D.shape[0], ST.shape[0]), dtype=float)        A = ST.T        for i in range(D.shape[0]):            C[i], _ = nnls(A, D[i])        return C    CT, *_ = np.linalg.lstsq(ST.T, D.T, rcond=None)    return CT.T  def update_ST(D, C, nonneg=False):    D = _validate_matrix(D, "D")    C = _validate_matrix(C, "C")    if C.shape[0] != D.shape[0]:        raise ValueError("C must have shape (I, N) matching D.")    if C.shape[1] < 1:        raise ValueError("The number of components N must be at least 1.")    if nonneg:        ST = np.empty((C.shape[1], D.shape[1]), dtype=float)        for j in range(D.shape[1]):            ST[:, j], _ = nnls(C, D[:, j])        return ST    ST, *_ = np.linalg.lstsq(C, D, rcond=None)    return ST  def mcr_als(D, ST0, max_iter=50, tol_percent=0.1, nonneg=False):    D = _validate_matrix(D, "D")    ST = _validate_matrix(ST0, "ST0")    if ST.shape[1] != D.shape[1]:        raise ValueError("Initial ST must have shape (N, J) matching D.")    history = []    sigma_prev = None    converged = False    C = None    for it in range(1, int(max_iter) + 1):        C = update_C(D, ST, nonneg=nonneg)        ST = update_ST(D, C, nonneg=nonneg)        Dhat = C @ ST        E = D - Dhat        sigma = residual_sd(E)        lof = lack_of_fit_percent(D, E)        sse = frobenius_sse(D, C, ST)        history.append(            {"iteration": it, "sigma": sigma, "lof": lof, "sse": sse}        )        if sigma_prev is not None and sigma_prev > 0:            rel = 100.0 * abs(sigma - sigma_prev) / sigma_prev            if rel <= float(tol_percent):                converged = True                break        sigma_prev = sigma    Dhat = C @ ST    E = D - Dhat    return {        "C": C,        "ST": ST,        "E": E,        "Dhat": Dhat,        "n_iter": history[-1]["iteration"],        "converged": converged,        "lof": lack_of_fit_percent(D, E),        "variance_explained": variance_explained(D, E),        "sigma": residual_sd(E),        "history": history,    } 

1. 1.

   Lawton, W. H., & Sylvestre, E. A. (1971). Self Modeling Curve Resolution. Technometrics, 13(3), 617-633.

   doi:10.1080/00401706.1971.10488823
2. 2.

   Tauler, R. (1995). Multivariate curve resolution applied to second order data. Chemometrics and Intelligent Laboratory Systems, 30(1), 133-146.

   doi:10.1016/0169-7439(95)00047-X
3. 3.

   Tauler, R., Smilde, A., & Kowalski, B. R. (1995). Selectivity, local rank, three-way data analysis and ambiguity in multivariate curve resolution. Journal of Chemometrics, 9, 31-58.

   doi:10.1002/cem.1180090105
4. 4.

   Jaumot, J., Gargallo, R., de Juan, A., & Tauler, R. (2005). A graphical user-friendly interface for MCR-ALS: a new tool for multivariate curve resolution in MATLAB. Chemometrics and Intelligent Laboratory Systems, 76, 101-110.

   doi:10.1016/j.chemolab.2004.12.007
5. 5.

   de Juan, A., Jaumot, J., & Tauler, R. (2014). Multivariate Curve Resolution (MCR). Solving the mixture analysis problem. Analytical Methods, 6, 4964-4976.

   doi:10.1039/C4AY00571F
6. 6.

   de Juan, A., & Tauler, R. (2021). Multivariate Curve Resolution: 50 years addressing the mixture analysis problem - A review. Analytica Chimica Acta, 1145, 59-78.

   doi:10.1016/j.aca.2020.10.051
7. 7.

   Jaumot, J., de Juan, A., & Tauler, R. (2015). MCR-ALS GUI 2.0: New features and applications. Chemometrics and Intelligent Laboratory Systems, 140, 1-12.

   doi:10.1016/j.chemolab.2014.10.003
8. 8.

   Camp, C. H., Jr. (2019). pyMCR: A Python Library for Multivariate Curve Resolution Analysis with Alternating Regression (MCR-AR). Journal of Research of the National Institute of Standards and Technology, 124, 124018.

   doi:10.6028/jres.124.018