Are MCRALS NonNegativity Constraints fully implemented yet?

[cedleau]

In the first place, thank you very much for your work on developing Spectrochempy, which I discovered very recently.

I installed developer version SpectroChemPy's API - v.0.6.6.dev13 and started playing aroud with the package.
Even though it was not mentionned in the tutorial, I noticed parameters in the sources to apply Non-Negativity Constraint, which seems to be implemented.

I ran the following code to test it hoping to see the effect of non-negativity, yet the two fittings return exactly the same results, with several negative spectrum values (circled in blue on the figure).

mcr = scp.MCRALS(log_level="INFO")
mcr.fit(X, St0)
_ = mcr.C.T.plot(cmap=None, color="red")

mcr = scp.MCRALS(log_level="INFO", nonnegSpec='all', nonnegConc='all' )
mcr.fit(X, St0)
_ = mcr.C.T.plot(cmap=None, color="red")

Hence I'd be curious to know if I'm missing something in the use of this functionnality, or if your developements are simply not completed yet (which would explain the lack of documentation on it).

Thank you for your help,,
Cedric

[fernandezc]

Hi,

Thanks for your interest and try on SpectroChemPy!

I'm not the principal developer of the MCRALS implementation (@atravert can probably answer better than me).

First, If you need feature information, you can look here https://www.spectrochempy.fr/latest/reference/generated/spectrochempy.MCRALS.html#spectrochempy.MCRALS, because the tutorial is indeed not the first thing we update when we develop.

Also, perhaps you could use solverSpec and solverConc to impose positive values.

[fernandezc]

I would add that it is not impossible to have bugs since the new implementation of MCRALS is quite recent. Perhaps you could attach the dataset you used for testing so that I can take a closer look to see if there are any problems.

[cedleau]

Hi,

Thank you for your quick answer,

Indeed, I missed the documentation that you shared, which is quite helpful. I managed to apply non-negativity to all spectra using 'nnls' for solverSpec and solverConc .

Yet, when using 'pnls' with applying negative constraint to all the columns, I face results that I did not expected. It looks like the negative condition was not properly applied. Here is some minimal example to reproduce, with dataset corresponding to a HPLC-DAD run, from Jaumot et al. Chemolab, 76 (2005), pp. 101-110 and Jaumot et al. Chemolab, 140 (2015) pp. 1-12, loaded from the “Multivariate Curve Resolution Homepage” at https://mcrals.wordpress.com/download/example-data-sets as ‘als2004dataset.MAT’. (This is the dataset from the tutorial.)

A = scp.read_matlab("matlabdata/als2004dataset.MAT")
X = A[-1]
St0 = A[3]

for param in [['lstsq', 'lstsq', 'all', 'all'],
              ['nnls', 'nnls', 'all', 'all'],
              ['pnnls', 'pnnls', [0,1,2,3], [0,1,2,3]],
              ['pnnls', 'lstsq', [0,1,2,3], 'all']]:
    mcr = scp.MCRALS(solverSpec=param[0], solverConc=param[1], nonnegSpec=param[2], nonnegConc=param[3])
    mcr.fit(X, St0)
    #_ = mcr.C.T.plot(cmap=None, color="red")
    for attr in ['solverSpec', 'solverConc', 'nonnegSpec', 'nonnegConc']:
        print((f'{attr}: {getattr(mcr,attr)}'))
    print(f'Number of negative concentration values: {(mcr.C.T.values<0).reshape(1,-1).squeeze().tolist().count(True)}')
    print(f'Number of negative spectral values: {(mcr.St.values<0).reshape(1,-1).squeeze().tolist().count(True)}')
    print('--')

returns

solverSpec: lstsq
solverConc: lstsq
nonnegSpec: [0, 1, 2, 3]
nonnegConc: [0, 1, 2, 3]
Number of negative concentration values: 27
Number of negative spectral values: 9
--
solverSpec: nnls
solverConc: nnls
nonnegSpec: [0, 1, 2, 3]
nonnegConc: [0, 1, 2, 3]
Number of negative concentration values: 0
Number of negative spectral values: 0
--
solverSpec: pnnls
solverConc: pnnls
nonnegSpec: [0, 1, 2, 3]
nonnegConc: [0, 1, 2, 3]
Number of negative concentration values: 23
Number of negative spectral values: 9
--
solverSpec: pnnls
solverConc: lstsq
nonnegSpec: [0, 1, 2, 3]
nonnegConc: [0, 1, 2, 3]
Number of negative concentration values: 27
Number of negative spectral values: 9
--

[atravert]

Thanks, I will look at this too.

[atravert]

Hello Cedric,

• first note that nonnegConc='all' is the default, so in principle there is no need to specify it.
• This behaviour is expected in this case because the ALS (Alternate Least Squares) loops are as follows (the + exponent indicates the matrix pseudo inverse, and the prime ' the matrix after application of the constraints:

1 - compute $C$ from $X$ and $S^t$ using least-squares: $C = X S^{t+}$
2 - apply constraints (including non-negativity) to $C$, yielding $C'$
3 - compute $S^t$ from $X$ and $C'$ using least-squares $S^t = C'^+ X$
4 - apply constraints to $S^t$ yielding ${S^t}'$
5 - compute $C$ from $X$ and ${S^t} '$ using least-squares $C = X {{S^t} '}^+$
6 - stop if converged else goto 1 using $S^t ={S^t}'$

The non-negative concentration profile ($C'$ above) is stored in the mcr.C_constrained attribute and you can plot it and see tthat it is indeed non-negative:

_ = mcr.C_constrained.T.plot()

The (slightly) negative points appear in the final concentration profile $C$. Lines 3-5 above are equivalent to:

$$ C = X {{S^t} '}^+ = X {({C'^+ X}) '}^+ $$

This formula shows that the final concentration profile $C$ is obtained from the constrained (non-negative) concentration profile $C'$ "filtered" by the data matrix $X$ and the optional constraints possibly applied on spectra ($'$). It is important to carry out these steps in order to check that the constraints you apply are indeed compatible with the data $X$. In this case, the negative points are - I think - not significant wrt the noise of the data matrix.

As indicated by @fernandezc , you can also directly enforce non-negativity using the nnls solver for concentrations and in this case, $C$ is indeed non-negative

mcr = scp.MCRALS(log_level="INFO")
mcr.solverConc = `nnls`
_ = mcr.fit(X, St0)
_ = mcr.C.T.plot()

Cuirrently nnls is not the default for compatibility with previous versions and also for efficiency ('nnls' increases the computation time by 40 to 50% wrt the default solver). In most cases - as here I think - the default parameters (lstsq solver + non-negativity constraints) allow reasonable enforcement of non negativity.

Finallly note that using the algorithm above yieds a concentration matrix $C$ (attribute mcr.C) directly obtained by least-squares and spectra matrix $St'$ ((attribute mcr.St) obtained by least-squares and possible application of constraints. We could have chosen the other way around ($St$ directly from least-squares and $C'$ aftert least-squares + constraints. If one prefers this solution, they can be obtained using the attributes mcr.C_constrained and mcr.St_unconstrained which store the results of lines 1 and 3 above.

I hope this helps. Let us know if you need further details: I will add this in the documentation, so your input and remarks are valuable to us !

Best regards
Arnaud

[cedleau]

Thank you very much for this extensive answer, very useful.
What I missed in the first place was that simple non-negativity constraint and nnls were two different things, it's much clearer now.
Best,