# Using separable non-negative matrix factorization techniques for the analysis of time-resolved Raman spectra

###### Abstract

The key challenge of time-resolved Raman spectroscopy is the identification of the constituent species and the analysis of the kinetics of the underlying reaction network. In this work we present an integral approach that allows for determining both the component spectra and the rate constants simultaneously from a series of vibrational spectra. It is based on an algorithm for non-negative matrix factorization which is applied to the experimental data set following a few pre-processing steps. As a prerequisite for physically unambiguous solutions, each component spectrum must include one vibrational band that does not significantly interfere with vibrational bands of other species. The approach is applied to synthetic “experimental” spectra derived from model systems comprising a set of species with component spectra differing with respect to their degree of spectral interferences and signal-to-noise ratios. In each case, the species involved are connected via monomolecular reaction pathways. The potential and limitations of the approach for recovering the respective rate constants and component spectra are discussed.

## 1 Introduction

Raman spectroscopy is a versatile tool to probe molecular structure changes that are associated with the temporal evolution of chemical or physical processes [21, 3, 20]. Since time-resolved Raman spectroscopy is applicable in a wide dynamic range down to the femtosecond time scale, it is capable to monitor quite different events, including intramolecular rearrangements in the excited state as well as chemical reactions in the ground state. The individual Raman spectra measured as a function of time represent a superposition of the intrinsic spectra of the individual species or molecular states that are involved in the reaction sequence. The relative contributions of the various spectra to each measured spectrum then reflect the actual composition of the sample at the respective time, and hence the entire series of experimental spectra represents the kinetics of the underlying processes. While just this information can also be provided by other transient optical techniques, the unique advantage of time-resolved Raman spectroscopy resides in the fact that the molecular structure of the individual species is encoded in the respective component spectra. This allows for identifying intermediate states and characterizing their structural and electronic properties as a prerequisite for determining the molecular reaction mechanism. The central task is, therefore, to disentangle the series of time-resolved Raman spectra in terms of the individual component spectra and their temporal evolution.

In the past decades, different concepts were designed for analyzing sets of complex Raman spectra [8, 14, 15, 22]. In most cases, these efforts were spurred by the concomitant vivid development in Raman and IR spectroscopic investigation of cellular systems for biological and medical applications.[23, 4] These studies are dedicated to classify and to distinguish microorganisms or to identify pathological tissues. To achieve these objectives, it is necessary to determine spectral signatures of a complex ensemble of biomolecules in a certain environment that are characteristic of a specific type or state of the cellular system. In these cases, analytical methods for pattern recognition are required which, in view of the large number of experimental spectra, are based on statistical procedures, such as principal component analysis, factor analysis, or singular value decomposition[8, 14, 15, 22, 23, 4, 27, 28]. Such approaches provide mathematical solutions but typically not the intrinsic spectra of the large number of pure components. Thus, extending multivariate analyses to a series of spectra reflecting physical or chemical changes of well-defined molecular species may just afford the number of the species involved but not necessarily their component spectra.

Such systems, on the other hand, are frequently treated by least-square methods in which either single Lorentzian/Gaussian bands and complete component spectra (component analysis) are employed to achieve a global fit to all experimental spectra[8]. The component analysis is quite robust as the number of degrees of freedom in the fitting process just corresponds to the number of components, given that all component spectra are known a priori. If, however, this is not the case, the “intuition” factor gains weight and thus the overall error increases substantially with the number of unknown component spectra.

In this work, we tried to overcome these drawbacks by developing an unsupervised analytical method that is based on non-negative matrix factorization (NMF). Such NMF techniques have rarely been used for factor analysis, with two notable exceptions[1, 19]. Unlike these previous approaches, our method takes advantage of the so-called separability inherent to the measurement data.

This paper is organized as follows. In Section 2 we present the mathematical model and describe the so-called separability condition that is derived from the specific properties of complex Raman spectra composed by a finite number of component spectra. Section 3 then describes the numerical method for computing the NMF and extracting the reaction coefficients. In Section 4 we illustrate on a sequence of artificial first-order reactions the reliability of our method under interference among the component spectra and under measurement noise. Concluding remarks are given in Section 5.

## 2 Mathematical background

### 2.1 Model for time-resolved Raman spectra of chemical reactions

Mathematically, the acquired measurements correspond to certain convex combinations of sums of Lorentz functions or Lorentzians that constitute the Raman spectra of the individual reactants. We will formalize our notion of the model in this section.

A Lorentzian is a non-negative “peak function” with its maximum at the base point (corresponding to the frequency of the normal mode), the width at half-height and the intensity , which is defined by

(1) |

For simplicity we will usually skip the parameters and simply write .

Consider a chemical reaction with reactant species. Then the Raman spectrum of each reactant can be modeled as a non-negative sum of Lorentzians, so that

(2) |

We assume that all base points (or “peaks”) are located in the finite interval . Note that one could also implement a model where the individual spectra are given by sums of Lorentzians and Gaussians, or Gaussians only.

Now we will consider the relative concentrations of the reactant species. We denote the concentration function of species by

so that corresponds to the relative concentration of species at time of the reaction. Consequently, at each time the concentrations sum to . The functions represent the reaction kinetics.

In this functional setting, the time-resolved vibrational Raman spectrum of the reaction can be modeled as

(3) |

Discretizing (3) over time through in time steps, and over frequencies through in frequencies, models the measured data from the reaction. We denote the resulting measurement matrix by .

Now we discretize the functions over the frequencies and obtain the vectors

as well as the matrix . Similarly, we discretize the functions over time and obtain the vectors

as well as the matrix

Then the measurement matrix , which includes all experimental spectra, can be written as

(4) |

i.e., the entry-wise non-negative matrix is the product of the two entry-wise non-negative matrices and .

### 2.2 NMF for Raman experimental data

The task of analyzing time-resolved Raman spectroscopic data consists of (1) determining the spectra of the individual species (component spectra), and (2) identifying the underlying reaction kinetics. Using the notation from Section 2.1, the corresponding mathematical problem is as follows: Given the non-negative data matrix , and assuming that the reaction involves species, find non-negative matrices and such that

(5) |

Note that a factorization of the form (5) where and/or have negative entries has no physical meaning, as neither a measured intensity nor a relative concentration can be negative.

The mathematical task of finding a non-negative matrix factorization (NMF) of is one of formidable difficulty: Without further assumptions on the given data, the problem is ill-posed and its solutions are non-unique in general. From a complexity point of view, NMF is NP-hard[26]. Moreover, an exact factorization as in (5) is more a theoretical desire than achievable in practice. The presence of noise or other forms of data uncertainty may simply rule out the existence of such a factorization.

The theory and computation of non-negative matrix factorizations is a very active research topic, whose popularity gained much from an article by Lee and Seung on the use of NMF for feature extraction[17]; an earlier work on NMF dates back to 1974[24]. A useful overview is given in the recent survey by Gillis[11].

A practically much better justified view is adopted by considering NMF as an approximation problem rather than an exact factorization. The usual way to give a formal definition of this approximation problem is as follows: Let be a matrix norm. Given the non-negative data matrix , we seek non-negative matrices and such that

is small. The norm plays the role of a “distance function” and measures how close the eventually found factors , reproduce the given data. In this work we mostly use the Frobenius norm, which for any (rectangular) matrix is defined by .

When treating time-resolved Raman data on chemical reactions, another problem of interest arises. Usually one does not know the number of intermediate species in the reaction. Thus the required information to set up an NMF problem, , is missing. To overcome this shortcoming, one may assume that is not very large (in practice it is often between two and ten). It then is a simple matter of trying different values for and selecting the best solution (see Section 4.5). Another heuristic estimate is available by the number of “large” singular values of the matrix [18].

### 2.3 The separability condition

While the general NMF problem introduced in Section 2.2 is very difficult in general, there is a very important special case, the so-called separable NMF problem [7, 2]. As before, the measurement matrix is denoted by , and the entries of each row of sum to . (This can always be achieved by applying a diagonal scaling matrix from the left.) Algebraically, the data matrix is called -separable, if it can be written in the form

(6) |

where is the -by- identity matrix, and is a permutation matrix. Separability implies that all the rows of can be reconstructed by using only rows of (these constitute the factor ) by convex combinations with weights given through .

The separability condition in equation (6) can be interpreted in the model for time-resolved Raman spectroscopy data from the previous section as follows. The measurement matrix is separable when each species , represented by the -th column of , has a characteristic frequency at which (see (2)), but for all other species . If such a frequency is present for each species, the measurement matrix contains rows that are equal to the rows of the sought-for kinetic matrix . Thus, the primary task is to determine the characteristic frequencies. For example, if all the species involved in the reaction contain a Lorentz band that does not interfere with a Lorentz band of any other species, then the separability condition is satisfied.

However, recall from the definition (1) that the intensity of a Lorentzian vanishes reciprocally to a quadratic function in the distance from the base frequency . In particular, the separability condition as explained above cannot be satisfied in any case, since for every species and frequency . Consequently, the factorization (4) will not exactly be separable (in the theoretical framework of Section 2.1).

While an exact interference-free set of characteristic frequencies is in general impossible, the interference can be numerically small or even zero if the corresponding base frequencies of the characteristic bands are sufficiently far apart. Algebraically, it means that instead of the exact separability condition (6), our problem is properly described by a near-separable problem, meaning that

(7) |

Here we interpret the matrix as noise originating from interference of the species at frequencies at which some species is strongly dominant in comparison with the other species. If is small, the original factors may be well approximated by applying an algorithm for separable NMF to . Quantitative investigations on the allowable noise for some algorithms have been persued in a purely mathematical context.[13, 9, 12]

In Section 4.3, we provide a numerical study on a model problem that shows the effect of growing on the overall approximation quality.

So far we considered only ideal, noise free data measurements. Real experimental data involve measurement noise, and hence

(8) |

where we assume the noise to be purely additive, i.e. is componentwise non-negative. Algebraically, the measurement noise is no different from the noise arising from interference. The effect of measurement noise is studied numerically in Section 4.4.

Note that the separability assumption is widely used in other fields, such as hyperspectral imaging, text mining, or other blind source separation applications [16, 5, 6]. In some of these contexts, the separability assumption is assumed to be satisfied with respect to the time axis (in our notation, for ). In our application this would mean that for each species there exists some point in time at which the relative concentration of the species is , which is highly unlikely for a typical reaction.

## 3 Computational method

The method we use for the identification of the reaction kinetics is based on the successive non-negative projection algorithm (SNPA) [10]. SNPA is an algorithm for computing the factor in (6) (the kinetics), provided that the problem at hand is near-separable. In the language of Section 2.3, we use SNPA to compute approximate characteristic frequencies of the species. We have chosen SNPA for its computational speed and robustness with respect to noise, but any other algorithm for separable NMF could be used as well.

While SNPA is at the center of our method, we also need to deal with a number of other computational tasks. The overall method is shown in Algorithm 1.

#### Step 1: Removing insignificant frequencies

In our approach it has turned out to be useful to remove all the experimental data (intensity-frequency pairs) at frequencies which did not display any significant intensities with respect to the measurement noise level. Here, we first estimate the noise level inherent to the measurement data by the standard deviation on the frequency having the least mean intensity. If we make the reasonable assumption that there are frequencies at which only noise is measured, such a frequency will have minimum mean intensity, and its standard deviation is a measure for the noise level. Algebraically, this filtering just leads to removal of some rows of . In order to simplify the notation, we will still assume that .

#### Step 2: Smoothing the data

An useful preprocessing step is to “smooth” the measurement data in direction of the time. In all numerical experiments in the following section that involve measurement noise, we smoothed the data using a running mean with a window size of . Algebraically, this smoothing just effects that each data entry of is replaced by the mean of the values for .

#### Step 3: Find characteristic frequencies

Subsequently, we use SNPA to find a set of approximate characteristic frequencies. If is the set of indices computed by SNPA, we obtain the pseudo-kinetic matrix by stacking the rows of indexed by (i.e., in Matlab-like notation). Note that the columns of may not sum to (approximately) , and hence cannot be interpreted as a reaction kinetic matrix.

#### Step 4: Scaling

In order to obtain a reaction kinetic matrix whose columns sum approximately to , we next compute a diagonal scaling matrix such that has this property. Here use that the approximation error of an NMF is invariant under such scalings, viz. . To find a suitable scaling matrix , we solve the non-negative least squares problem

where . After finding the optimal scaling values , we rescale the kinetic matrix . In our experiments described in the following section, we used Matlab’s lsqnonneg function.

#### Step 5: Compute corresponding spectra

#### Step 6: Extracting reaction coefficients

To extract the reaction coefficients (rate constants) from a given kinetic matrix , we restrict the analysis to the case that all reaction steps are of first order. Recall from Section 2.1 that if is the reaction kinetics matrix of the true kinetics function , and is the matrix of reaction coefficients, then , where is the initial concentration vector. If is not known, it can be recovered from by solving the nonlinear least squares problem

where denotes the -th column of . If has full row rank, then is uniquely determined by . In our experiments we used Matlab’s fminunc function.

## 4 Numerical study with synthetic data

In this section we illustrate the effectiveness of our method using a sequence of artificial first-order reactions involving five species. We first describe the model reactions and the component spectra of the involved species (“fingerprints”) in Section 4.1. The results in Section 4.2 show that in the low-interference, noiseless regime, both the kinetics and species fingerprints are perfectly recovered. Then, in Section 4.3, we study numerically the effect of increasing interference among the species, and in Section 4.4 we also add measurement noise to the data and study the recovery quality. Finally in Section 4.5, we address the question of determining the correct number of species.

### 4.1 Description of the model reaction

The reaction scheme is set-up by five species A, B, C, D and E which are inter-related by first-order reactions. These reactions are characterized by reaction coefficients (in arbitrary units of reciprocal time) as given as follows:

A \arrow(–bb)<=>[0.53][0.02] B \arrow<=>[0.43][0.25] C \arrow<=>[0.11][0.1] E \arrow(@bb–)->[*00.21][-90] D \schemestop

We let species A be the only educt in the reaction, resulting in the initial concentration vector . If we denote the reaction coefficient matrix corresponding to the above reaction scheme by

(9) |

the reaction kinetics are given as a function over time by

as displayed in Figure 1 (left). The corresponding kinetics matrix is obtained by discretization of at equidistant steps , so that .

The five component spectra are constructed by arbitrarily chosen combinations of Lorentzians, inspired by the Raman spectra of various organic compounds. These component spectra were constructed such that each species has at least one characteristic frequency, so that the imposed separability condition (see Section 2.3) is satisfied. The fingerprints constitute the columns of the matrix (Figure 1 (right)). Visual inspection already reveals that each of the five species has at least one characteristic frequency.

The resulting data matrix is obtained as the product of the kinetic- and spectra matrix, viz. . A visualization of the data matrix is given in Figure 2 (top).

### 4.2 Recovery in the noiseless case

Given the data corresponding to the artificial reaction scheme described in Section 4.1, our goal is now to recover both the component spectra of the five species and the reaction kinetics only from these data, i.e. to recover the matrices and using only the data in without any further information. By construction of the data, we know that is separable, and we will now apply the methods described in Section 2.3.

The results are displayed in Figure 3, demonstrating that both the reaction kinetics and all the component spectra are recovered to such a high level of accuracy, that they can hardly be distinguished visually from the original data.

Nevertheless, the computed factors are not identical to the true factors and . For the relative errors of and we find

Hence the relative error for both the kinetics and the species is less than 1%.

Finally, we will recover the reaction coefficient matrix from the computed kinetics matrix . Applying the methodology described in Section 3, we obtain

Compared with the true reaction coefficients in (9), we find that the largest error made in estimating any of the reaction coefficients from the computed kinetics is , related to the coefficient .

In summary, we find that all the data that constitute the reaction network described in Section 4.1 have been recovered quite accurately.

### 4.3 Effect of increased interference

In the model used in the previous section, the bands of the species involved were (visually) well separated from each other. In the other extreme, i.e., if all bands of a species interfere with those of other species, our approach will not be applicable for the recovery of the reaction kinetics and the species fingerprints, as the factors are no longer close to a separable factorization (the term in (7) becomes too large). Thus we now consider the case of “modest” interference.

We enforce an increased level of interference among the species by moving all the base points in all species towards three focal points (see Figure 4).

The result of Algorithm 1 being applied to these interference-rich data is shown in Figure 5. Because of the increased interference, the computed kinetics deviate slightly from the true ones, but all the species bands have been identified quite satisfactorily. For the relative errors of and we find

so the relative error for both the kinetics and the species spectra is not more than and , respectively.

If, however, the bands in the original species spectra are moved even closer to each other, our method will eventually fail to produce a qualitatively good solution.

### 4.4 Recovery under the influence of measurement noise

The data used in the previous section are highly idealized in the sense that they are free of measurement noise. In any practical setting, the experimentally acquired Raman measurements will be contaminated with noise from different sources, such as signal shot noise or background noise (e.g. fluorescence).

We will now simulate the effect of noise added to the measurement matrix in the sense of (8). We assume that all the noise from different sources taken together resemble additive Gaussian white noise. Overall we disturb the data matrix according to

where the entries of are drawn from the normal distribution and is the relative noise level. Further we will assume for our noise model, that any constant background has already been removed from the measurement data. In Figure 2 (bottom), we visualize the resulting noisy data .

The component spectra and kinetics were recovered as described in Section 3. Figure 6 shows the result of our algorithm applied to the noisy, interference-rich data. Because of the large distance of to being truly separable, the computed kinetics displays some deviations from the true data, but still provides a satisfactory description. Also the computed component spectra show a reasonable agreement with the true spectra.

### 4.5 Determination of the number of species

In the numerical experiments described above we have always assumed that the number of species is known (here ). Of course, in practical applications the correct determination of the value of is one of the greatest challenges. One approach for solving this problem is to use Algorithm 1 for different values of , which yields approximate species and kinetics , and then compute the (relative) data error .

In the following table we show the data errors for and the three experimental setups from Sections 4.2–4.4: noiseless case (first row of the table), increased interference (second), increased interference and measurement noise (third). In each case a significant drop in the data error occurs at the correct number of species, while no significant further reduction of the error is achieved when increasing the number of (suspected) species even further.

1 | 2 | 3 | 4 | 5 | 6 | 7 | |
---|---|---|---|---|---|---|---|

Sec. 4.2 | 0.97294 | 0.77203 | 0.18479 | 0.06055 | 0.00003 | 0.00003 | 0.00003 |

Sec. 4.3 | 0.94908 | 0.68692 | 0.20209 | 0.04698 | 0.00029 | 0.00029 | 0.00029 |

Sec. 4.4 | 0.83106 | 0.35383 | 0.25894 | 0.12683 | 0.08169 | 0.08107 | 0.08066 |

The results indicate that using the NMF for determining the number of species is an alternative to existing techniques in this context such as singular values of the data matrix. An extensive survey of the latter technique is given in [18].

## 5 Concluding remarks

We have presented an algorithm for the recovery of component spectra and reaction kinetics from data obtained through time resolved Raman spectroscopy. The key tool we used in our approach is the so-called separability condition in the non-negative matrix factorization (NMF). In terms of the component spectra this condition is approximately satisfied, if each of the species has at least one band in its spectrum that does not interfere too much with the bands of other species. Thus, it is conceptually similar to the classical Rayleigh criterion limit. Our approach combines a standard algorithm for separable NMF (“SNPA”) with a few pre-processing steps for the measurement data. A number of other separable NMF algorithms exist, but we did not yet pursue a detailed comparison among them for the present application. In our numerical study we have demonstrated that the component spectra and reaction kinetics can be recovered with a reasonable quality under modest measurement noise and interference among the component spectra. Whereas in terms of the kinetics the approach currently restricted to a network of first-order or pseudo-first-order reactions, it may equally well applied to other spectroscopic techniques such as IR or NMR spectroscopies.

#### Acknowledgements

This research was supported by the DFG cluster of excellence “UniCat”.

## References

- 1. M. Ando and H.-o. Hamaguchi, “Molecular Component Distribution Imaging of Living Cells by Multivariate Curve Resolution Analysis of Space-Resolved Raman Spectra”. J. Biomed. Opt. 2013. 19(1):011016.
- 2. S. Arora, R. Ge, R. Kannan, and A. Moitra, “Computing a Nonnegative Matrix Factorization—Provably”. In STOC’12—Proceedings of the 2012 ACM Symposium on Theory of Computing, pages 145–161, ACM, New York, 2012.
- 3. G. Balakrishnan, C. L. Weeks, M. Ibrahim, A. V. Soldatova, and T. G. Spiro, “Protein Dynamics from Time Resolved UV Raman Spectroscopy.” Curr. Opin. Struct. Biol. 2008. 18(5):623–629.
- 4. A. Beljebbar, O. Bouché, M. D. Diébold, P. J. Guillou, J. P. Palot, D. Eudes, and M. Manfait, “Identification of Raman Spectroscopic Markers for the Characterization of Normal and Adenocarcinomatous Colonic Tissues.” Crit. Rev. Oncol. Hematol. 2009. 72(3):255–264.
- 5. J. M. Bioucas-Dias, A. Plaza, N. Dobigeon, M. Parente, Q. Du, P. Gader, and J. Chanussot, “Hyperspectral Unmixing Overview: Geometrical, Statistical, and Sparse Regression-Based Approaches”. IEEE J. Sel. Top. Appl. Earth Observations Remote Sensing 2012. 5(2):354–379.
- 6. T.-H. Chan, W.-K. Ma, C.-Y. Chi, and Y. Wang, “A Convex Analysis Framework for Blind Separation of Non-Negative Sources”. IEEE Trans. Signal Process. 2008. 56(10, part 2):5120–5134.
- 7. D. Donoho and V. Stodden, “When Does Non-Negative Matrix Factorization give a Correct Decomposition into Parts?” In S. Thrun, L. Saul, and B. Schölkopf, editors, Advances in Neural Information Processing Systems 16, MIT Press, Cambridge, MA, 2004.
- 8. S. Döpner, P. Hildebrandt, A. G. Mauk, H. Lenk, and W. Stempfle, “Analysis of Vibrational Spectra of Multicomponent Systems. Application to pH-dependent Resonance Raman Spectra of ferricytochrome c”. Spectrochim. Acta, Part A 1996. 52(5):573 – 584.
- 9. N. Gillis, “Robustness Analysis of Hottopixx, a Linear Programming Model for Factoring Nonnegative Matrices”. SIAM J. Matrix Anal. Appl. 2013. 34(3):1189–1212.
- 10. N. Gillis, “Successive Nonnegative Projection Algorithm for Robust Nonnegative Blind Source Separation”. SIAM J. Imaging Sci. 2014. 7(2):1420–1450.
- 11. N. Gillis, “The Why and How of Nonnegative Matrix Factorization”. In Chapman & Hall/CRC Machine Learning & Pattern Recognition, pages 257–291, Chapman and Hall/CRC, 2014.
- 12. N. Gillis and R. Luce, “Robust Near-Separable Nonnegative Matrix Factorization using Linear Optimization”. J. Mach. Learn. Res. 2014. 15:1249–1280.
- 13. N. Gillis and S. A. Vavasis, “Fast and Robust Recursive Algorithms for Separable Nonnegative Matrix Factorization”. IEEE Trans. Pattern Anal. Mach. Intell. 2014. 36(4):698–714.
- 14. R. W. Hendler and R. I. Shrager, “Deconvolutions Based on Singular Value Decomposition and the Pseudoinverse: a Guide for Beginners.” J. Biochem. Biophys. Methods 1994. 28(1):1–33.
- 15. E. Henry and J. Hofrichter, “Singular Value Decomposition: Application to Analysis of Experimental Data”. In Numerical Computer Methods, volume 210 of Methods in Enzymology, pages 129 – 192, Academic Press, 1992.
- 16. A. Kumar, V. Sindhwani, and P. Kambadur, “Fast Conical Hull Algorithms for Near-Separable Non-Negative Matrix Factorization”. ICML 2013. (PART 1):231–239.
- 17. D. D. Lee and H. S. Seung, “Learning the Parts of Objects by Non-Negative Matrix Factorization”. Nature 1999. 401(6755):788–791.
- 18. E. R. Malinkowski, Factor Analysis in Chemistry. John Wiley & Sons, Inc., New York, 3rd edition, 2002.
- 19. K. Neymeyr, M. Sawall, and D. Hess, “Pure Component Spectral Recovery and Constrained Matrix Factorizations: Concepts and Applications”. J. Chemom. 2010. 24(2):67–74.
- 20. M. Ondrias, M. Simpson, and R. Larsen, “Time-Resolved Resonance Raman Spectroscopy”. In J. Laserna, editor, Modern Techniques in Raman Spectroscopy, Wiley, 1996.
- 21. S. K. Sahoo, S. Umapathy, and A. W. Parker, “Time-Resolved Resonance Raman Spectroscopy: Exploring Reactive Intermediates”. Appl. Spectrosc. 2011. 65(10):1087–1115.
- 22. H. Shinzawa, K. Awa, W. Kanematsu, and Y. Ozaki, “Multivariate Data Analysis for Raman Spectroscopic Imaging”. J. Raman Spectrosc. 2009. 40(12):1720–1725.
- 23. A. C. S. Talari, C. A. Evans, I. Holen, R. E. Coleman, and I. U. Rehman, “Raman Spectroscopic Analysis Differentiates Between Breast Cancer Cell Lines”. J. Raman Spectrosc. 2015. 46(5):421–427.
- 24. L. Thomas, “Rank Factorization of Nonnegative Matrices (A. Berman)”. SIAM Rev. 1974. 16(3):393–394.
- 25. M. H. Van Benthem and M. R. Keenan, “Fast Algorithm for the Solution of Large-Scale Non-Negativity-Constrained Least Squares Problems”. J. Chemom. 2004. 18(10):441–450.
- 26. S. A. Vavasis, “On the Complexity of Nonnegative Matrix Factorization”. SIAM J. Optim. 2009. 20(3):1364–1377.
- 27. A. T. Weakley, Warwick, P. C. Temple, T. E. Bitterwolf, and D. E. Aston, “Multivariate Analysis of Micro-Raman Spectra of Thermoplastic Polyurethane Blends Using Principal Component Analysis and Principal Component Regression”. Appl. Spectrosc. 2012. 66(11):1269–1278.
- 28. A. Zhang, W. Zeng, T. M. Niemczyk, M. R. Keenan, and D. M. Haaland, “Multivariate analysis of Infrared Spectra for Monitoring and Understanding the Kinetics and Mechanisms of Adsorption Processes”. Appl. Spectrosc. 2005. 59(1):47–55.