Classification of Copper Minerals by Handheld Laser-Induced Breakdown Spectroscopy and Nonnegative Tensor Factorisation

Laser-induced breakdown spectroscopy (LIBS) analysers are becoming increasingly common for material classification purposes. However, to achieve good classification accuracy, mostly noncompact units are used based on their stability and reproducibility. In addition, computational algorithms that require significant hardware resources are commonly applied. For performing measurement campaigns in hard-to-access environments, such as mining sites, there is a need for compact, portable, or even handheld devices capable of reaching high measurement accuracy. The optics and hardware of small (i.e., handheld) devices are limited by space and power consumption and require a compromise of the achievable spectral quality. As long as the size of such a device is a major constraint, the software is the primary field for improvement. In this study, we propose a novel combination of handheld LIBS with non-negative tensor factorisation to investigate its classification capabilities of copper minerals. The proposed approach is based on the extraction of source spectra for each mineral (with the use of tensor methods) and their labelling based on the percentage contribution within the dataset. These latent spectra are then used in a regression model for validation purposes. The application of such an approach leads to an increase in the classification score by approximately 5% compared to that obtained using commonly used classifiers such as support vector machines, linear discriminant analysis, and the k-nearest neighbours algorithm.


Introduction
Laser-induced breakdown spectroscopy (LIBS) [1] is a remote sensing technique used for both qualitative and quantitative analysis of various materials. The operational principle is to use a high pulse energy laser to instantaneously heat the matter to evaporate a small amount of the substrate and eject it as a plasma plume. Then, the light emitted by the plasma is dispersed and registered by a camera. After the specified time of continuous wavelength radiation, a time window exists with the quickly cooling plasma, whereas the individual spectral lines representing the elemental material composition can be recorded.
To shed light on the given classification problem and understand the general similarities of the minerals, they were all grouped in accordance with the Strunz Classes [23] given in Table 1. The first three characters of each mineral name (indicated with bold font in the above list) together with a Strunz Class number constitute short names used in the remainder of this article (with one exception for cornetite being co0-8). The exact Strunz Class assignment for each mineral is given in Table S1 in the supplementary data. The material suppliers were chosen from among different countries and/or geographical regions to differentiate the population of rocks with specific minerals and make the database more versatile. The origins of the samples include six continents, 29  The mineral shapes and sizes attached to the rocks were extremely diverse, from perfect crystals spread on the surface of the base rock to thin layers, often combined with the base rock, with unknown percentage distributions of both. Figure 1 presents examples of such mineral shapes. Examples are for chalcopyrite, azurite, and malachite, where each pair of photos (Figure 1a-f, respectively) shows two cases: (1) a well-built and clear crystal structure without impurities attached to the base rock, and (2) evenly distributed mineral over the rock surface with an unknown percentage mix with the base rock. Moreover, the attached crystal sizes differed from large, as in the case of chalcopyrite (Figure 1a), to very small as in the case of azurite (Figure 1c).
Sensors 2020, 20  A total of 127 rock samples were taken into consideration. The statistics for the number of rock samples per mineral are presented in Figure S2 in the supplementary data.

Element Composition
As the chemical bonds are destroyed during the plasma creation, the analysed spectral signal is primarily dependent on the elements that existed in the material before the laser action. All 62 minerals are described in Table S1 in the supplementary data by their empirical element composition taken from [24], where a total of 27 elements were found. Then, based on the elemental composition, an investigation of the similarities between them was conducted with the use of PCA to find theoretically undistinguishable minerals by LIBS. In this way, a new N-dimensional (N = 27) space was created with new orthogonal variables called principal components (PCs) that describe the dataset sorted from the highest to the lowest variance [25]. Even with such PCA transformation, it was impossible to find good representative PCs to determine how the minerals are distinguishable from each other in a lower dimensional space (i.e., 3D representation).
In that case, N-dimensional measures are used to find mineral pairs that are most like each other. For this purpose, the Matlab ® pdist function was used with four selected distance measures: [26]. The final computation was an average of those four measures as they gave similar yet different orders of hard to distinguish mineral pairs. Figure 2a presents the least distinguishable 21 pairs of minerals. The worst cases were the mal-5 and azu-5 groups as well as a group consisting of lan-7, pos-7, and bro-7. These two barely reach 2% of the maximum distance between minerals in the PCA space. Figure 2b presents the worst-case groups of minerals extracted from Figure 2a with the other minerals distributed within the space of the first two PCs. Those groups are the already mentioned lan-7/pos-7/bro-7 and mal-5/azu-5 as well as co0-8/lib-8/pse-8, cli-8/cor-8/oli-8, dio-9/chr-9/ajo-9, hal-9/all-9, gla-2/kru-2, and the last one, tsu-8/vau-7, which is actually the first pair that shows close similarity above the Strunz Class division; hence, there is no need to search for another similar groups. Further analysis of the elemental composition values ( Figure 3) clearly indicated that malachite and azurite as well as the group of langite, posnjakite, and brochantite cannot be distinguished from each other. As a result, new mineral labels were created: a+m for the mal-5/azu-5 and blp for the other three. These minerals were combined and labelled together for further classification. The full heatmap with all 27 elements and 62 minerals is shown in Figure S3 in the supplementary data.  Further analysis of the elemental composition values ( Figure 3) clearly indicated that malachite and azurite as well as the group of langite, posnjakite, and brochantite cannot be distinguished from each other. As a result, new mineral labels were created: a+m for the mal-5/azu-5 and blp for the other three. These minerals were combined and labelled together for further classification. The full heatmap with all 27 elements and 62 minerals is shown in Figure S3 in the supplementary data. Further analysis of the elemental composition values ( Figure 3) clearly indicated that malachite and azurite as well as the group of langite, posnjakite, and brochantite cannot be distinguished from each other. As a result, new mineral labels were created: a+m for the mal-5/azu-5 and blp for the other three. These minerals were combined and labelled together for further classification. The full heatmap with all 27 elements and 62 minerals is shown in Figure S3 in the supplementary data.

Experimental
The measurements were performed using a handheld LIBS device (Z-300, SciAps). The DPSS Nd:YAG laser emits radiation at a wavelength of 1064 nm, with a repetition rate up to 50 Hz, a pulse energy of 5-6 mJ, and a pulse duration of 1-2 ns. The spectrometer covers a range of 190-950 nm. Immediately before and during the measurement, the measuring region was purged with Ar gas to remove ambient air and enhance the LIBS signals. Each specific measurement point (MP) on the rock sample consisted of 64 single shots in an 8 × 8 grid covering an area of roughly 2-4 mm 2 . It is worth noting that the device did not save spectra below a specific intensity threshold. Therefore, in some cases, fewer than 64 spectra were recorded per MP.
The mineral crystals on the surface were uniform, and the area covered by the 8 × 8 grid often extended beyond the crystal of interest. For that reason, all MPs had to be manually labelled.
An example of such labelling is shown in Figure 4 and Table 2.

Experimental
The measurements were performed using a handheld LIBS device (Z-300, SciAps). The DPSS Nd:YAG laser emits radiation at a wavelength of 1064 nm, with a repetition rate up to 50 Hz, a pulse energy of 5-6 mJ, and a pulse duration of 1-2 ns. The spectrometer covers a range of 190-950 nm. Immediately before and during the measurement, the measuring region was purged with Ar gas to remove ambient air and enhance the LIBS signals. Each specific measurement point (MP) on the rock sample consisted of 64 single shots in an 8 × 8 grid covering an area of roughly 2-4 mm 2 . It is worth noting that the device did not save spectra below a specific intensity threshold. Therefore, in some cases, fewer than 64 spectra were recorded per MP.
The mineral crystals on the surface were uniform, and the area covered by the 8 × 8 grid often extended beyond the crystal of interest. For that reason, all MPs had to be manually labelled. An example of such labelling is shown in Figure 4 and Table 2.   When an MP seems to cover more than 90% of the interesting mineral, it is flagged as a possible reference, and can then be used for both training and validation. If not flagged with logic 1 in the Reference column, then the MP remains as validation only. In the Coverage column of the description table, a percentage value of the mineral of interest within the 8 × 8 grid is provided. These values are rough assessments (limited to a difference of +/−5%) made just after the LIBS measurement, based on the pattern of laser spots on the sample and descriptions of the minerals of interest made by geologists. This will assist with further automated validation of MPs (because if they reach that percentage value, it means that classification succeeded). Even if the coverage description may in some cases be slightly inaccurate, this does not favour our proposed NTF method as the bias remains the same for all classification algorithms used.
During further classification, the MP is a single entity and the 8 × 8 grid spots will not be separated or analysed exclusively. In total, 458 MPs were recorded, of which 311 were flagged as possible reference and 147 as validation only. The full statistics of MPs per mineral are given in Figure  S1 in the supplementary data.

Latent Spectrum Extraction
The latent spectrum extraction in its principal form is used to select some underlying (latent) spectra from observed LIBS spectra that are considered as mixtures of multiple latent components generated by multiple spectral sources. The latent spectra are more frequent and have a common   Figure 4.

MP Mineral
Reference Coverage   1  azurite  0  80  2  azurite  0  20  3  azurite  1  90  4  azurite  0  30  5  azurite  0  60 When an MP seems to cover more than 90% of the interesting mineral, it is flagged as a possible reference, and can then be used for both training and validation. If not flagged with logic 1 in the Reference column, then the MP remains as validation only. In the Coverage column of the description table, a percentage value of the mineral of interest within the 8 × 8 grid is provided. These values are rough assessments (limited to a difference of +/−5%) made just after the LIBS measurement, based on the pattern of laser spots on the sample and descriptions of the minerals of interest made by geologists. This will assist with further automated validation of MPs (because if they reach that percentage value, it means that classification succeeded). Even if the coverage description may in some cases be slightly inaccurate, this does not favour our proposed NTF method as the bias remains the same for all classification algorithms used.
During further classification, the MP is a single entity and the 8 × 8 grid spots will not be separated or analysed exclusively. In total, 458 MPs were recorded, of which 311 were flagged as possible reference and 147 as validation only. The full statistics of MPs per mineral are given in Figure S1 in the supplementary data.

Latent Spectrum Extraction
The latent spectrum extraction in its principal form is used to select some underlying (latent) spectra from observed LIBS spectra that are considered as mixtures of multiple latent components generated by multiple spectral sources. The latent spectra are more frequent and have a common pattern across all observed LIBS spectra. In an ideal case, after the extraction, one latent spectrum should resemble an artificial spectrum of the desired source. In practise, we obtain few latent components that approximately represent the true spectra of the analysed minerals.
In the investigated case, each object is probed with several to a dozen measuring points.
Thus, each MP contains a maximum of 64 spectra. Let y + be the i 1 -th spectrum of the m-th object, measured in the i 2 -th MP. Each MP is assumed to provide I is usually lower than 64 because some shots that correspond to the spectra of a low variance (below a threshold) must be neglected. The spectral resolution is determined by the number of samples in each spectrum, that is, the number I (m) 3 . Because the spectral resolution is the same for each observed spectrum, then ∀m : can be regarded as a superposition of latent spectra that could be pure spectra of analysed minerals (endmembers) or other unwanted or perturbing spectra. The latent spectra for the m-th object can be collected into the matrix U (m, 3) where J m is the number of latent spectra in the m-th object. Considering the above, the spectrum y (m) i 1 ,i 2 can be expressed by the following superposition rule: where the coefficient ξ i 1 ,i 2 ,j m ≥ 0 determines the contribution of the j m -th latent spectrum to the i 1 -th observed spectrum of the i 2 -th MP in the m-th object. The coefficient ξ i 1 ,i 2 ,j m can then be factorised as ≥ 0 represents the contribution from the i 1 -th shot and u 2 . Sweeping over the indices i 1 and i 2 , be a 3-way array (3-modal tensor) created from a set of spectra y (m) i 1 ,i 2 for the m-th object. It is thus easy to notice that model (1) takes the form where the symbol • denotes the outer product. Model (2) can be equivalently expressed in the form Sensors 2020, 20, 5152 is a superdiagonal identity tensor, and the symbol × n stands for the tensor-matrix product across the n-th mode.
Factor U (m, 3) contains J m latent spectra, and U (m,1) and U (m,2) represent the contribution coefficients (concentrations) of the latent spectra to observations across the first and second modes of the tensor Y (m) , respectively. The latent spectra can thus be obtained by performing the NTF of Y (m) , given the assumed number J m .
There are numerous computational strategies for NTF, and nearly all of them are based on an alternating optimisation scheme with unfolding imposed on each mode. Model (3) expressed in the unfolded version takes the form where is a matrix obtained by the unfolding tensor Y (m) along its n-th mode; where n = 1, 2, 3, and the symbol stands for the Khatri-Rao product. Note that the system in (4) is considerably overdetermined because p n I (m) p J n . To alleviate the problem of scaling ambiguity in the NTF, the columns in matrices U (m,2) and U (m,3) are normalised to the unit l 1 -norm. The system of linear equations in (4) can be solved with numerous linear solvers subject to nonnegativity constraints. In our study, we used the hierarchical alternating least-squares (HALS) algorithm proposed in [29], and then computationally improved in [30]. It belongs to a family of block coordinate descent The graphical representation of the NTF model is presented in Figure S4 of the supplementary data.

Regression Model Using Latent Spectra
To estimate the percentage rate of minerals in a newly measured MP, the latent spectra are extracted from the known MPs labelled with the so-called strong reference. We assume that J m latent spectra are extracted from the m-th labelled object using the NTF. The latent spectra after being postprocessed can be regarded as regressors for predicting the percentage of minerals in such an unknown MP. Any unknown spectrum y ∈ R I 3 + is assumed to be approximated by a linear regression model, Coefficient α where λ ≥ 0 is a regularisation parameter that controls the overfitting. In this study, problem (6) was solved using the interior-point least-squares algorithms for regularised box-constrained problems, implemented in the function lsqlin in Matlab ® 2016a. Note that coefficients in vector α * can also be regarded as unknown percentages of expected minerals in the analysed MP.

Determining the Number of Latent Spectra per Mineral
The last problem to be solved is how to determine the correct (sufficient) J m number of latent spectra for each m-th object. The problem would be trivial if only pure mineral samples were analysed, since only a single latent spectrum may represent the desired mineral. However, in our measurements, impurities in the spectra resulting from base rocks and other minerals will occur. Moreover, the NTF is sensitive to the difference in light intensity distribution between spectrometer channels, which is not constant because the light is propagated inside the device and thus can cause additional perturbations of the desired m-th object spectra. In that case, the assumption was made that we will search for J m latent spectra in which most of them represent the scattered spectra of the desired mineral. The less contributing examples after the so-called self-regression will be labelled as unknown spectra (U).
In such a case, we developed an iterative process in which we start from J m = 2 (J m = 1 will be similar to a weighted average of the data) and increment the number until reaching a break loop condition. The loop breaking condition is based on two parameters: level L and ratio R, which will be set up prior to this procedure.
The L-condition is superior and relates to the cumulative contribution given by the sum of coefficients α (m) j m of the first K latent spectra (j m = 1, . . . , K), sorted in descending order of their contributions. Note that number K cannot be equal to J m , which means that at least one latent spectrum should be classified as undefined U. If the cumulative contribution of K latent spectra in a given iterative step is equal to or greater than the L value, then the first K spectra are assigned as the desired mineral, and the last J m -K spectra are labelled as U. Otherwise, J m is incremented and again the L-condition is checked.
If the L-condition is satisfied, then the second R-condition should be checked. This condition requires the minimum ratio between the last (K-th) spectrum assigned as a mineral to the first (J m -K+1) spectrum assigned as U. This relation is very important, as for a higher value of J m , the distribution of latent spectra contributions might end up equal. Therefore, it may be difficult to assess the U spectra, as the difference between the last labelled as mineral and first labelled as U may be only a few percentage points. If the R-condition is not met, J m is incremented, and the L-condition is checked.
After meeting these two (L and R) conditions, we achieve the solution for which we have the desired level of contributions and the ratio of the contributions of the last mineral spectrum to the first U spectrum large enough to assume the U spectra are really the unwanted signals. An example of such a loop operation with given L/R conditions is shown in Figure S5 of the supplementary data.
Obviously, for some combination of R and L in a given m-th object, such a pair of conditions can never be met. To avoid an infinite loop, parameter J , then a given pair of L/R conditions is excluded from consideration for all m objects. Figure 5 presents an example of such latent spectra extraction for malachite with R set to 1.8 and L to 0.96. In this case, three components were extracted as scattered in the m-th object and one as the U spectrum.
For the final validation process of a new unknown MP, the scattered K latent spectra of each m-th object become different regressors. Finally, the predicted percentage contribution for the m-th object in the new MP becomes the sum of the single K contributions of the m-th object scattered latent spectra. The classifying label is then decided on the mineral that has the top percentage contribution. never be met. To avoid an infinite loop, parameter ( ) is introduced, which is the maximum Jm number to be incremented in the loop. If the conditions are not met after = ( ) , then a given pair of L/R conditions is excluded from consideration for all m objects. Figure 5 presents an example of such latent spectra extraction for malachite with R set to 1.8 and L to 0.96. In this case, three components were extracted as scattered in the m-th object and one as the U spectrum.

Setup
The number of MPs flagged as validation only differs among minerals, some do not even have such or have only one. To make the proportion of training and validation data equal among the minerals, some of the MPs flagged as possible reference were also moved to be validation data. This new division can be observed in Figure S6 of the supplementary data.
The parameters that were set for our algorithm were a maximum number of latent spectra (J (max) m = 20), ratios R (from 1.1 to 2.0), and levels L (from 0.80 to 0.98). The following sections present 10 × 10 heatmaps with classification measures from which we can select the best parameters for the analysis of the confusion matrix.
Because we deal with an imbalanced dataset, there is a risk of the results being overwhelmed by the outcome of the larger mineral classes. Thus, the metrics of precision and recall are introduced together with their bounding metric called F-measure [31]: where tp is the true positive rate, fp is the false positive rate, and fn is the false negative rate. The proposed NTF-based classification also introduces the U class (Section 4.3), which should be included in some way in the metrics in (7). For that purpose, we propose to differentiate two cases: Uin and Uex.
The Uin case assumes that we include MPs labelled as U in the classification score calculation. However, the metrics in (7) require, among others, true positive rates, which in the case of U class, will never exist (U is not a true class). This will result in precision and recall being zero all the time and independently from the number of MPs labelled as U, the F-measure (if calculated) for U class will always be zero and, therefore, lower the total F-measure for all the classes together. The only way to introduce the U class into the F-measure, therefore, is assuming that if any MP is going to be assigned as U class, it will appear as a false negative value for the original true class (i.e., azurite labelled as U).
The Uex case totally excludes MPs labelled as U from the classification score calculation. This decision was made on the assumption that U labels give us the information that the measurement outcome is uncertain and the user should repeat the measurement on such a sample for more confidence (this is not yet an error at that point but restrains us from introducing false positive rates into another class).
Again, as U is not a true class, the false positive rates will never exist for it, so, actually, the only difference between Uin and Uex scores will be held by the recall part of the equation in (7) (the precision will stay the same for both).
The NaN values within the heat maps are related to the L/R pairs with which the algorithm could not find a latent extraction solution for at least one mineral. For improved clarity, each heatmap has the top three L/R solutions listed in its title.
As the confusion matrices in the case of our 59 classes were very large, they were added as supplementary data (Figures S8-S11).

Analysis of Training Data
The first portion of the results is focussed on finding the ideal combination of parameters for our proposed method. To do this, the model was trained and validated using the same training data ( Figure S6-supplementary data). This is reasonable, as while the proposed algorithm extracts the latent spectra, it does it within a single mineral class and is unaware of the existence of other classes. Because of this, it is unlikely for the model with so many classes to reach a 100% score even when trained and validated with the same data (which is different in the case of the SVM that reached 100% under the same conditions). However, this gives us an opportunity to revise the model on the basis of training data and determine the best combination of R and L values for the given dataset.
From Figure 6, we see that the F-measure score increased with increasing L value. Moreover, it is clear that we reached the parameterisation boundary from each side of the 10 × 10 matrices, as it was impossible to deliver results for L equal or greater than 0.98 and R equal to 1.8 (or greater). If the L values were taken from the range [0.8, 0.82], we had few possible solutions and much lower scores. Finally, it was futile to use R values smaller than 1.1 because 1.0 would designate equality.
The best results for both F Uin (Figure 6a) and F Uex (Figure 6b) are a combination of R = 1.5 and L = 0.96, and those parameters will be selected as a final solution with the validation dataset.
of training data and determine the best combination of R and L values for the given dataset.
From Figure 6, we see that the F-measure score increased with increasing L value. Moreover, it is clear that we reached the parameterisation boundary from each side of the 10 × 10 matrices, as it was impossible to deliver results for L equal or greater than 0.98 and R equal to 1.8 (or greater). If the L values were taken from the range [0.8, 0.82], we had few possible solutions and much lower scores. Finally, it was futile to use R values smaller than 1.1 because 1.0 would designate equality. The best results for both FUin (Figure 6a) and FUex (Figure 6b) are a combination of R = 1.5 and L = 0.96, and those parameters will be selected as a final solution with the validation dataset.

Analysis of Validation Data
For this section, analyses were taken with separate training and validation data, in accordance with Figure S6 in the supplementary data. Although we already selected parameters in the previous section, we decided to perform the same parameterisation of R and L, resulting in 10 × 10 matrices to verify that the selection was accurate.
Similar to Figure 6, the results in Figure 7 indicate (with few exceptions) that the F-measure score increased with an increase in L. The boundary parameters also remained the same as the learned model, as in the previous section. Figure 7 presents the F-measure results of the parameterisation, but apart from the scores for the Uin and Uex cases, it also presents the results in comparison to the best selected classifier ( Figure  7c,d), which, in our case, was SVM scoring 67.22%. It is important to note here that this was the highest value that was possible to reach for an SVM trying different kernel functions and their parameterisation. In fact, the basic linear SVM scored the best among all SVM variants that was investigated. The additional results for NTF covering the F-measure's partial scores-precision and recall-are presented in Figure S7 in the supplementary data.

Analysis of Validation Data
For this section, analyses were taken with separate training and validation data, in accordance with Figure S6 in the supplementary data. Although we already selected parameters in the previous section, we decided to perform the same parameterisation of R and L, resulting in 10 × 10 matrices to verify that the selection was accurate.
Similar to Figure 6, the results in Figure 7 indicate (with few exceptions) that the F-measure score increased with an increase in L. The boundary parameters also remained the same as the learned model, as in the previous section. Figure 7 presents the F-measure results of the parameterisation, but apart from the scores for the Uin and Uex cases, it also presents the results in comparison to the best selected classifier (Figure 7c,d), which, in our case, was SVM scoring 67.22%. It is important to note here that this was the highest value that was possible to reach for an SVM trying different kernel functions and their parameterisation. In fact, the basic linear SVM scored the best among all SVM variants that was investigated. The additional results for NTF covering the F-measure's partial scores-precision and recall-are presented in Figure S7 in the supplementary data.
Although the best results in the case of Figure 7 were not for R = 1.5 and L = 0.96, these results were selected as final because they could be foreseen (Section 5.2) and did not differ much from the other high scores, especially in the case of the highest results such as F Uex . Table 3 presents the discussed measures obtained for the NTF-based method and three standard classifiers for the analysed validation data. The best built-in solution was SVM, which reached an F-measure of approximately 17.5% better than that associated with KNN and LDA. The proposed algorithm reached 1.87% (Uin) and 5.02% (Uex) of the F-measure score with respect to the second-best SVM, followed by increases in every other case where R Uex had the highest gain of 6.92%.   Table 3; (d) FUex-scores with U class excluded from measure in accordance with the best SVM result from Table 3.
Although the best results in the case of Figure 7 were not for R = 1.5 and L = 0.96, these results were selected as final because they could be foreseen (Section 5.2.) and did not differ much from the other high scores, especially in the case of the highest results such as FUex. Table 3 presents the discussed measures obtained for the NTF-based method and three standard classifiers for the analysed validation data. The best built-in solution was SVM, which reached an Fmeasure of approximately 17.5% better than that associated with KNN and LDA. The proposed algorithm reached 1.87% (Uin) and 5.02% (Uex) of the F-measure score with respect to the secondbest SVM, followed by increases in every other case where RUex had the highest gain of 6.92%. The confusion matrices for the built-in classifiers present a trend to seek a host for the more likely hard-to-assess samples. In the case of the best out of three SVM ( Figure S8 Table 3; (d) F Uex -scores with U class excluded from measure in accordance with the best SVM result from Table 3.
The confusion matrices for the built-in classifiers present a trend to seek a host for the more likely hard-to-assess samples. In the case of the best out of three SVM ( Figure S8-supplementary data), 23 mineral classes had a 100% score, and the classes that caused the most errors were blp (19 false positives) and a+m (12 false positives). For LDA ( Figure S9-supplementary data), only 11 mineral classes had a 100% score, and the most confusing was all-5 (16 false positives). The KNN (Figure S10-supplementary data) performed similarly to LDA for the final F-measure score, while 13 mineral classes were perfectly assessed, and the classes that caused the most errors were a+m (19 false positives) and ajo-5 (13 false positives). Using the NTF method ( Figure S11-supplementary data), we managed to classify the top F-measure score and the top 25 mineral classes without error. The U class perfectly took the top host position for the hardest to assess samples (17 false positives), followed by cha-2 (14 false positives) and cup-4 (8 false positives), which is very reasonable.
The CPU time and disk space usage of the above methods are compared in Table S12 in the supplementary data. It is clear from this table that the proposed NTF-based method requires less disk space and needs half the time for data validation compared to the competing SVM.

Example of Mineral Contribution for Selected MPs
The regression method proposed, apart from classifying the MPs, also gives the mineral percentage contribution followed by their geometrical distribution within the 8 × 8 shooting grid. Figure 8 presents the results for MP1 from the selected chalcopyrite sample. The bar plot presents the single contributions of the regressors (J m latent spectra) summed for each mineral. In this MP case, the top count goes to chalcopyrite (almost 80%), leaving all the other minerals far behind, so the given classification label is correct. The geometrical distributions of the four main minerals (m-th objects) and the unknown class were equal within the 8 × 8 grid. Both can be confirmed with the MP photo, as the surface looks like an equal mineral distribution without any sign of the base rock.
pattern. This fact is clearly visible in the mineral distribution heatmaps of a+m, cha-2, and chr-9, where in the case of the proper azurite (a+m) class in the bottom-left quarter, we observe a higher contribution of that mineral, while in the same area on the cha-2 and chr-9 heatmaps, there is almost zero contribution. The total sum of a+m contribution on the bar plot is almost 30%, which perfectly covers one-quarter of the crystal that fits within the 8 × 8 shooting grid plus some average error of the latent spectra contribution.  However, as mentioned in Section 2, there were also samples that did not have a consistent mineral distribution on the surface, and the crystal size was even smaller than the 8 × 8 shooting grid. An example of such a situation is the selected azurite sample presented in Figure 9. From the photo, we observe that the azurite crystal actually covers less than one-quarter of the visible, burned in laser pattern. This fact is clearly visible in the mineral distribution heatmaps of a+m, cha-2, and chr-9, where in the case of the proper azurite (a+m) class in the bottom-left quarter, we observe a higher contribution of that mineral, while in the same area on the cha-2 and chr-9 heatmaps, there is almost zero contribution. The total sum of a+m contribution on the bar plot is almost 30%, which perfectly covers one-quarter of the crystal that fits within the 8 × 8 shooting grid plus some average error of the latent spectra contribution.  In both cases, as the classification labels were correctly set, the U class did not play an important role in the final percentage rate within those MPs. Both MPs are missing some spectra.

Conclusions
In this research, we hypothesised that it is possible to extract artificial (latent) spectra for each of In both cases, as the classification labels were correctly set, the U class did not play an important role in the final percentage rate within those MPs. Both MPs are missing some spectra.

Conclusions
In this research, we hypothesised that it is possible to extract artificial (latent) spectra for each of the investigated minerals and use them as predictors in a linear regression model. For this purpose, NTF was used. Because of the heterogeneity of the mineral samples and weak reproducibility of the spectra acquired with the use of a handheld LIBS device, the procedure for proper latent spectra selection was proposed. In such a procedure, a parameterisation of ratio R and level L variables is required. The results show that these variables can be limited, with high accuracy, to one selection by performing validation with the use of the same data as for the model training.
The NTF-based classification performed well, reaching higher F-measure scores of around 1.9% (when the Unknown class was included in the measures) and around 5.0% (when the Unknown class was excluded from the measures), both in accordance with the best SVM classifier. The standard methods seem to find a host class for the hardest-to-assess samples, so using the method that already contains an Unknown class inside its model was even more reasonable here.
In addition to the output labelling required for classification purposes, a percentage contribution of the minerals within measuring points is given. Such additional information on the MPs allows the creation of 2D mapping of the shooting area 8 × 8 grids with a smooth distribution of the mineral classes among them.
The final regression model can be stored using low disk space, and the regression function is not as memory-and CPU-intensive (compared to the commonly used classifiers SVM, LDA, and KNN). Such a model may even be implemented in mobiles and other handheld LIBS devices and thus increase their functionality as fast, on-site mineral analysers.
Future research will be devoted to applying the method to minerals other than those containing copper as their primary element and verifying its universal usability. The NTF was confirmed as a method for extracting artificial spectra for the mineral classes and was successfully used as a regressor in a linear model.

Supplementary Materials:
The following are available online at http://www.mdpi.com/1424-8220/20/18/5152/s1, Table S1: detailed list of 62 minerals used in research; Figure S1: the number of measurement points (MPs) taken per mineral together with their partitions into possible reference and validation only ones; Figure S2: the number of different (exclusive) rock samples on which MPs were taken-listed per mineral; Figure S3: heatmap representing elements composition among all 62 minerals in accordance with the mineral's empirical formula [24]; Figure S4: NTF model; Figure S5: example of J m number of latent spectra decision for one m-th object; Figure S6: dataset partitioning of the training and validation data using some of the possible reference MPs as validation data to evenly match the dataset among mineral classes; Figure S7: classification precision and recall output for the parametrisation of the NTF method with R and L variables (validation dataset), Figure S8: confusion matrix of SVM classifier for validation dataset; Figure S9: confusion matrix of KNN classifier for validation dataset; Figure S10: confusion matrix of LDA classifier for validation dataset; Figure S11: confusion matrix of NTF based classifier for validation dataset; Table S12: CPU time and disk space usage of the compared methods.