Neural Density Estimation and Uncertainty
Quantification for Laser Induced Breakdown
Spectroscopy Spectra

Katiana Kontolati
Johns Hopkins University
Baltimore, MD 21218
kontolati@jhu.edu

Natalie Klein
Los Alamos National Laboratory
Los Alamos, NM 87545
neklein@lanl.gov

Nishant Panda
Los Alamos National Laboratory
Los Alamos, NM 87545
nishpan@lanl.gov

Diane Oyen
Los Alamos National Laboratory
Los Alamos, NM 87545
doyen@lanl.gov

Abstract
Constructing probability densities for inference in high-dimensional spectral data
is often intractable. In this work, we use normalizing flows on structured spectral
latent spaces to estimate such densities, enabling downstream inference tasks. In
addition, we evaluate a method for uncertainty quantification when predicting
unobserved state vectors associated with each spectrum. We demonstrate thecapability of this approach on laser-induced breakdown spectroscopy data collectedby the ChemCam instrument on the Mars rover Curiosity. Using our approach, weare able to generate realistic spectral samples and to accurately predict state vectors
with associated well-calibrated uncertainties. We anticipate that this methodology
will enable efficient probabilistic modeling of spectral data, leading to potential
advances in several areas, including out-of-distribution detection and sensitivity
analysis.

1

Introduction

The ChemCam instrument on the Mars rover Curiosity uses laser-induced breakdown spectroscopy(LIBS), a type of an atomic emission spectroscopy, to remotely analyze Martian rocks (Wienset al., 2012). Spectral information is used to extract the qualitative and quantitative chemical content(QQCC) of a material sample, where the QQCC can be seen as an unobserved state vector representing
the sample. Both linear and nonlinear supervised learning techniques have been applied for mapping
LIBS spectra to QQCC with good accuracy (Forni et al., 2013; Boucher et al., 2015; Castorena et al.,2021). In this work, we build upon these efforts by proposing a framework for constructing theprobability density function (PDF) of LIBS spectra. In addition, unlike previous methods which
produce point estimates of QQCC, we propose and evaluate a method for uncertainty quantification(UQ) on point predictions of QQCC.
Many real-world spectra, including LIBS spectra, are characterized by a large number of features,complicating the construction of the data PDF due to the curse of dimensionality. We develop a novel
framework for constructing low-dimensional PDFs suitable for downstream inference (e.g., sampling,
density estimation, outlier detection, or unsupervised representation learning) using state-of-the-art
Fourth Workshop on Machine Learning and the Physical Sciences (NeurIPS 2021).

neural density estimators with normalizing flows (NF) on spectral latent spaces. This framework
allows us to generate realistic spectral samples on a reduced space and using an inverse-transformation
project them back to the physically interpretable space. Furthermore, we propose a bootstrapping
approach for quantifying uncertainty in predictions of unobserved state vectors corresponding toeach spectrum. We demonstrate the capabilities of the proposed approach to construct the PDF of a
LIBS spectral data set, and to learn a mapping to the known QQCC with uncertainty. The validated
framework can be then employed for QQCC prediction and direct UQ of novel samples, such as
artificial samples generated by the NF model as well as spectra collected directly on Mars.
To the best of our knowledge, this work is the first time normalizing flows are constructed on spectral
latent spaces and can be readily employed for any kind of spectroscopy data. We show that the
proposed framework provides a straightforward way to perform downstream inference tasks and direct
UQ for high-dimensional spectral data. Codes can be found at https://github.com/katiana22/
uq4chemcam.

2

Methods

2.1

Problem statement

Let’s assume y ∈ RM is an M −dimensional random vector with non-negative elements, and a true
data distribution pY , which represents the spectral signals. Our goal is to learn an invertible, stable
mapping between the approximate data distribution p̂Y and a latent distribution pZ (e.g., Gaussian)that will allow for fast evaluation of various inference tasks. However, estimating the full-joint densityof very high-dimensional spectra is a challenging and often intractable task. Therefore, we introduce
a second mapping, that transforms y ∈ RM to x ∈ RL where L ≪ M to discover the spectral latent
representation of the signals. Next, we learn an invertible mapping between x ∼ p̂X (spectral latentvariable) and z ∼ pZ (latent variable). This framework allows us to generate novel samples on the
reduced spectral latent space x ∼ p̂X and use the inverse transformation to map back to the original
space RL → RM and therefore approximate the true data distribution.
We also want to estimate an unobserved state vector v ∈ RC where C is the vector dimensionality.
For the ChemCam application we consider, this represents the QQCC, an 8-dimensional vector withthe relative weight percentages of 8 major oxides commonly found on Mars. Given a training datasetof LIBS spectra y ∈ RM and associated compositions v ∈ RC (samples generated on Earth), weare interested in constructing a surrogate of the mapping f : RM → RC , that will allow us to makepredictions of the chemical concentration of novel samples. To calculate uncertainties related to thesepredictions, we propose an approach based on bootstrapping that allows us to quantify both modeland data uncertainties and thus assign measures of accuracy to sample estimates. This approach canbe then employed for UQ of data generated by the normalizing flow model.
2.2

Spectral NMF latent space

Consider N observations of the random vector y ∈ RM and let the data matrix be Y =N ×M[y1 , y2 , .., yN ]T ∈ R≥0. We use non-negative matrix factorization (NMF) to decompose Y
×Linto a product of a non-negative basis matrix X ∈ RN
and a non-negative coefficient matrix
≥0
PL
L×M
V ∈ R≥0 , such that Y ≈ XV (or equivalently yj ≈ i=1 xi Vij ) (Paatero and Tapper, 1994;
Wang and Zhang, 2012). NMF decomposes each data point into the linear combination of the basis
vectors. The NMF optimization problem consists of minimizing the Frobenius norm between Y andXV .
NMF is appropriate for non-negative data and a powerful method for feature selection; and thus allows
us to ignore the non-informative LIBS dimensions and enable interpretability of results (Rammelkamp
et al., 2020).
2

2.3

Inference via a spectral normalizing flow

We propose the construction of a normalizing flow model on the latent space of the LIBS spectra,
obtained by the NMF decomposition, to learn the underlying probability distribution of the spectral
latent variable x ∈ RL . The inverse NMF mapping introduced in the previous section can be used toproject generated samples back to the physically interpretable space (i.e., RL → RM ).
Normalizing flows are a powerful class of likelihood generative models which transform a basedensity into a target density by a series of deterministic and invertible transformations (Kobyzev et al.,2020). Consider the base density pZ (z) (spectral latent variable), the more complex density pX (x)(latent variable) and an invertible mapping x = f (z). Under the change of variables formula we cancompute the log-likelihood of x as
∂f (x|θ)(1)
log(pX (x|θ)) = log(pZ (z|θ)) + log det∂xT
where θ represents the trainable parameters of the flow. To train the NF model the negative loglikelihood (NLL) of Eq.(1) is minimized.
Here, we parameterize the normalizing flow with a sequence of real-valued non-volume preserving(RealNVP) transformations (Dinh et al., 2016). The RealNVP model, composes two types of
invertible transformations: additive coupling layers and rescaling. The model uses the so-calledaffine coupling layers for the coupling flows, which are simple and computationally efficient. The
transformation can be written as
x1:d = z1:d
xd+1:D = zd+1:D ⊙ exp(fα (z1:d )) + fµ (z1:d ),(2)where ⊙ is the Hadamard product or element-wise product and the exp(·) is applied to each elementof α. The above transformation performs a ‘1-1’ mapping to the first d elements and scales and shifts
the remaining D − d. By incorporating coupling layers into the flow, the elements are permuted
across layers so that a different set of elements is copied each time. Here we model fα , fµ as neuralnetworks. Once the PDF is learned, downstream inference tasks can be performed straightforwardly.
In the next Section, we are interested in predicting the elemental composition of novel samples, with
their associated uncertainty.
2.4

Uncertainty quantification via bootsrapping

We now aim to construct a mapping between the LIBS signal signatures and QQCC. We train shallow
neural networks, one for each oxide element. The models are formed asv (i) = σ(wT x(i) + b),

i = 1, .., t(3)

where t is the total number of oxides to be determined, w, b denote the trainable weights and bias andσ is the activation function. We learn the parameters of the models with a training set {(yi , vi )}N
i=1 ,
where N is the total number of samples. The main advantage of such models is that they are both fast
to train and result in very good accuracy scores (see Results).
To quantify uncertainties related to predictions of elemental compositions we use bootstrapping(Kumar and Srivastava, 2012), a statistics resampling method which allow us to assign measures ofaccuracy to a sample estimate. In general, bootstrapping performs as well as parametric prediction
intervals and its implementation is straightforward. For our application, it does not result in highcomputational cost given the choice of simplistic bootstrapped shallow neural networks. In case of
more complex models, methods that leverage the last-layer of the network can be employed as they
have shown good performance (Brosse et al., 2020). Given a new observation y0 ∈ RM , we canwrite
v0 = v̂n (y0 ) + r(y0 )(4)where v̂n (y0 ) represents the model estimate at the n-th bootstrap iteration (model uncertainty) andr(y0 ) the predicted residual between true and predicted values which can been modeled for a training
3

Figure 2: Regression results of the composition of each of the 8 major oxides for a holdout set of
140 samples. Accuracy is measured with the coefficient of determination (R2 score) and x,y axes
represent the true and predicted oxide wt.% respectively.
Table 1: Coverage results for 95% confidence intervals and 140 validation samples.
oxidecoverage (%)

SiO2
84.89

TiO2
98.56

Al2 O3
86.33

FeOT
86.33

MgO
86.33

CaO
96.40

Na2 O
93.53

K2 O
89.93

dataset with a regression model (data uncertainty). To measure the quality of prediction intervals, we
compute the coverage of validation samples (the rate at which the actual values fall within the rangeof the prediction interval). A summary of the proposed approach is presented in Algorithm 1.
Algorithm 1 Proposed methodology
Require: Input spectra y ∈ RM , state vectors v ∈ RC , embedding dimension L
Ensure: y ≥ 0
PLPerform NMF yj ≈ i=1 xi Vij , and obtain data on the latent space, x ∈ RL
Construct density of latent variable pX (x) with RealNVP
for i=1,t do
Train a shallow NN as v (i) = σ(wT x(i) + b)
▷ Map reduced spectra to state vectors
end forGiven a new observation y0 ∈ RM , prediction is v0 = v̂n (y0 ) + r(y0 ), n : bootstrap iteration
Compute coverage to evaluate the confidence intervals

3

Results

Consider Y ∈ RN ×M the LIBS spectra matrix with N = 426 andM = 5606. We perform NMF using L = 15 (selected by 5-foldcross-validation) with transformed data matrix X ∈ RN ×L . Next weconstruct a NF model based on the RealNVP architecture with 5 coupling layers and a Gaussian distribution as a base density (Kostrikov(2019)). We should highlight here the computational advantagesof this approach. Constructing a NF model on the 15-dimensional
latent space is extremely fast as the training process required lessthan 1 minute of CPU time. In Figure 1, we show a novel random
sample generated by the NF model which is transformed back to the Figure 1: Random sample generated by the normalizing flow
original space with inverse NMF.
model and transformed to the
We consider 8 oxides and therefore we construct t = 8, single hidden original space.
layer neural network models, with ReLU activation function, trainedwith stochastic gradient descent (SGD). To measure the accuracy of results we compute the coefficient4

of determination, calculated as R2 = 1 − RSS
TSS , where RSS represents the sum of squares of residuals
and TSS the total sum of squares. We show the accuracy of models in Figure 2 for a holdout set of
140 samples. The results reported show that point estimates are close to the optimal regression (1:1
line) and overall the response is comparable to state-of-the-art deep CNN approaches (Castorena
et al., 2021). For TiO2 and K2 O elements, estimates show larger deviation, which are not consideredsignificant due to the small oxide wt.% values at these regions.
Finally, we perform bootstrapping for a dataset of LIBS
spectra samples collected on Earth (with associated groundtruth) and we compute the prediction intervals for the same
holdout dataset of 140 samples. In Figure 3, we plot the
distribution of bootstrap predictions with the ground truth
for a random Earth sample where box plots represent thefirst quartile to the third quartile and green stars are theground truth. Table 1 shows the coverage calculated for
95% prediction intervals for all holdout samples, and we
see that the intervals appear to nearly achieve the nominal
coverage. The validated approach can be therefore used Figure 3: Bootstrap results for a random
for making predictions with uncertainty for novel LIBS sample.
spectra, generated either by the normalizing flow modelor for Martian samples directly collected from ChemCam.

4

Conclusions

We showed that the proposed framework provides a straightforward way to perform downstream
inference tasks for high-dimensional spectral data by identifying a spectral latent space, estimated asa parsimonious representation of the data, and constructing a spectral normalizing flow model on
the reduced space. The proposed approach is ideal for modeling high-dimensional data and enables
the learning of complex distributions in a fast and efficient way. Beyond general high-dimensional
inference, the proposed UQ approach allows for predictions of state vectors associated with novel
out-of-distribution data or data generated by the trained normalizing flow. We believe that this work
does not have future societal or ethical consequences.

Acknowledgments and Disclosure of Funding
This project was supported by the Laboratory Directed Research and Development program of Los
Alamos National Laboratory under project number LDRD-20210043DR. Research was performed
while K.K. was an Applied Machine Learning Summer Research Fellow at LANL.

5

Checklist• Do the main claims made in the abstract and introduction accurately reflect the paper’s
contributions and scope?: Yes• Have you read the ethics review guidelines and ensured that your paper conforms to them?:
Yes• Did you discuss any potential negative societal impacts of your work?: As we mention in thelast sentence of Section 4 we believe that this work does not have negative societal impacts.
• Did you describe the limitations of your work?: The main limitation of this work is related tothe high computational cost of the bootstrapping method in cases of complex neural networkmodels. We mention this limitation and also propose an alternative solution in Section 2.4.
• Did you include the code, data, and instructions needed to reproduce the main experimental
results (either in the supplemental material or as a URL)?: Yes
5• Did you specify all the training details (e.g., data splits, hyperparameters, how they werechosen)?: Yes• Did you report error bars (e.g., with respect to the random seed after running experiments
multiple times)?: Yes• Did you include the amount of compute and the type of resources used (e.g., type of GPUs,
internal cluster, or cloud provider)?: We have used CPUs and run the codes on a personal
computer. Therefore we do not discuss about specific computational resources.
• If your work uses existing assets, did you cite the creators?: Yes.
• Did you mention the license of the assets?: Yes.
• Did you include any new assets either in the supplemental material or as a URL?: Yes.
• Did you discuss whether and how consent was obtained from people whose data you’reusing/curating?: bluen/a.
• Did you discuss whether the data you are using/curating contains personally identifiable
information or offensive content?: The data we are using do not contain such content.

References
Wiens, R. C.; Maurice, S.; Barraclough, B.; Saccoccio, M.; Barkley, W. C.; Bell, J. F.; Bender, S.; Bernardin, J.;
Blaney, D.; Blank, J., et al. The ChemCam instrument suite on the Mars Science Laboratory (MSL) rover:
Body unit and combined system tests. Space science reviews 2012, 170, 167–227.
Forni, O.; Maurice, S.; Gasnault, O.; Wiens, R. C.; Cousin, A.; Clegg, S. M.; Sirven, J.-B.; Lasue, J. Independent
component analysis classification of laser induced breakdown spectroscopy spectra. Spectrochimica Acta Part
B: Atomic Spectroscopy 2013, 86, 31–41.
Boucher, T. F.; Ozanne, M. V.; Carmosino, M. L.; Dyar, M. D.; Mahadevan, S.; Breves, E. A.; Lepore, K. H.;
Clegg, S. M. A study of machine learning regression methods for major elemental analysis of rocks using
laser-induced breakdown spectroscopy. Spectrochimica Acta Part B: Atomic Spectroscopy 2015, 107, 1–10.
Castorena, J.; Oyen, D.; Ollila, A.; Legget, C.; Lanza, N. Deep spectral CNN for laser induced breakdown
spectroscopy. Spectrochimica Acta Part B: Atomic Spectroscopy 2021, 178, 106125.
Paatero, P.; Tapper, U. Positive matrix factorization: A non-negative factor model with optimal utilization oferror estimates of data values. Environmetrics 1994, 5, 111–126.
Wang, Y.-X.; Zhang, Y.-J. Nonnegative matrix factorization: A comprehensive review. IEEE Transactions on
knowledge and data engineering 2012, 25, 1336–1353.
Rammelkamp, K.; Gasnault, O.; Forni, O.; Lasue, J.; Maurice, S. Optimization of clustering analyses forclassification of ChemCam data from Gale Crater, Mars. European Planetary Science Congress. 2020; pp
EPSC2020–867.
Kobyzev, I.; Prince, S.; Brubaker, M. Normalizing flows: An introduction and review of current methods. IEEE
Transactions on Pattern Analysis and Machine Intelligence 2020,
Dinh, L.; Sohl-Dickstein, J.; Bengio, S. Density estimation using real NVP. arXiv preprint arXiv:1605.08803
2016,
Kumar, S.; Srivastava, A. Bootstrap prediction intervals in non-parametric regression with applications toanomaly detection. Proc. 18th ACM SIGKDD Conf. Knowl. Discovery Data Mining. 2012.
Brosse, N.; Riquelme, C.; Martin, A.; Gelly, S.; Moulines, É. On last-layer algorithms for classification:
Decoupling representation from uncertainty estimation. arXiv preprint arXiv:2001.08049 2020,
Kostrikov, I. Pytorch Flows (MIT License). https://github.com/ikostrikov/pytorch-flows, 2019.

6