An Initial Exploration of Bayesian Model Calibration for Estimating the Composition of Rocks and Soils on Mars
Claire-Alice Hébert, Earl Lawrence, Kary Myers, James P. Colgan, Elizabeth J. Judge
RReceived: 00 Month 0000 Revised: 00 Month 0000 Accepted: 00 Month 0000DOI: xxx/xxxx
RESEARCH ARTICLE
An Initial Exploration of Bayesian Model Calibration forEstimating the Composition of Rocks and Soils on Mars
Claire-Alice Hébert* | Earl Lawrence | Kary Myers | James P. Colgan | Elizabeth J. Judge Department of Applied Physics, StanfordUniversity, 348 Via Pueblo, Stanford, CA94305, USA Statistical Sciences, Los Alamos NationalLaboratory, New Mexico, USA Physics and Chemistry of Materials, LosAlamos National Laboratory, New Mexico,USA Chemical Diagnostics and Engineering,Los Alamos National Laboratory, NewMexico, USA
Correspondence *Claire-Alice Hébert. Email:[email protected]
Summary
The Mars Curiosity rover carries an instrument, ChemCam, designed to measure thecomposition of surface rocks and soil using laser-induced breakdown spectroscopy(LIBS). The measured spectra from this instrument must be analyzed to identify thecomponent elements in the target sample, as well as their relative proportions. Thisprocess, which we call disaggregation, is complicated by so-called matrix effects,which describe nonlinear changes in the relative heights of emission lines as anunknown function of composition due to atomic interactions within the LIBS plasma.In this work we explore the use of the plasma physics code ATOMIC, developed atLos Alamos National Laboratory, for the disaggregation task. ATOMIC has recentlybeen used to model LIBS spectra and can robustly reproduce matrix effects fromfirst principles. The ability of ATOMIC to predict LIBS spectra presents an excitingopportunity to perform disaggregation in a manner not yet tried in the LIBS commu-nity, namely via Bayesian model calibration. However, using it directly to solve ourinverse problem is computationally intractable due to the large parameter space andthe computation time required to produce a single output. Therefore we also explorethe use of emulators as a fast solution for this analysis. We discuss a proof of conceptGaussian process emulator for disaggregating two-element compounds of sodiumand copper. The training and test datasets were simulated with ATOMIC using a Latinhypercube design. After testing the performance of the emulator, we successfullyrecover the composition of 25 test spectra with Bayesian model calibration.
KEYWORDS:
Bayesian calibration, laser-induced breakdown spectroscopy, LIBS, ATOMIC, Gaussian process, emula-tion, disaggregation, modular calibration
One of the main scientific drivers of the Mars rover Curiosityis to determine whether Mars has ever been host to forms oflife [13]. ChemCam, one of the instruments on board, devel-oped by Los Alamos National Laboratory and L’Institut deRecherche en Astrophysique et Planétologie, is designed torecord detailed data about the surface soil and rocks of the planet. The instrument uses laser-induced breakdown spec-troscopy (LIBS) to measure the abundance of all chemicalelements by firing a laser onto a small patch of rock or soilsurface, producing a plasma. As the high-temperature plasmacools, it emits light that ChemCam records via a spectrometerand CCD camera, producing a detailed spectrum over a rangeof wavelengths. Intensity peaks in such a spectrum can be usedto identify the presence and relative abundance of chemicalspecies in the sample of rock or soil. The absence or presence a r X i v : . [ s t a t . A P ] A ug Hébert
ET AL of certain elements could hold important evidence about thepossibility of life on the rocky planet.Curiosity can obtain hundreds of these spectra every day, butanalyzing them remains a slow manual process. To success-fully analyze, or disaggregate , each spectrum, we must answertwo questions: what elements were present in the sample, andwhat are their relative proportions or abundances?An expert can often answer the first question, identifyingcomponent elements by the presence of signature peaks inthe spectrum at specific wavelengths. Estimating the relativeabundance of these constituent species, however, can be quitedifficult due to atomic interactions within the plasma. Due tothese interactions, the relative heights of emissions lines oftwo elements can have a nonlinear dependence on the rela-tive abundance of those species. These matrix effects — where matrix refers to the components of the target rather than tothe mathematical concept — complicate our disaggregationproblem: the spectrum of a multi-species target is not sim-ply the linear combination of the spectra for each individualelement. These nonlinear effects pose significant challengesfor disaggregation of LIBS data, particularly because we don’thave a closed form expression of the effects.Clegg et al. [2] introduced the use of multivariateapproaches, such as partial least squares regression, to improvethe ability to determine the elements present in a sample andestimate their relative proportions in the presence of matrixeffects. More recently, scientists at Los Alamos National Lab-oratory adapted a first-principles plasma physics code calledATOMIC [12] to provide a forward model for the emissionfrom LIBS plasmas, including their matrix effects [3, 4, 8].The new existence of this forward model presents an excitingopportunity to explore the use of Bayesian model calibration [9] to compute estimates, with associated uncertainties, of thecomponents of a target and their relative abundances. That is,we want to solve the inverse problem to determine the inputparameters for ATOMIC, which include the elements presentand their proportions, that produce the simulated spectrum thatis most like an unknown measured spectrum.Bayesian model calibration in the context of high-dimensional outputs [7], of which LIBS spectra are anexample, has a successful history in a variety of scientificapplications. For example, in materials science, these methodshave been applied to estimate strength parameters of aluminumalloys in hydrodynamic shock experiments [16].We draw inspiration from Judge et al. [8], who demonstratedthe impact of matrix effects on sodium peak heights in a sim-ple two-element mix of sodium and copper measured with Note that we use the term disaggregation for our work to allow for eventualconsideration of other types of problems, such as estimating the devices drawingpower from a household given a single measurement of the household’s total powerusage. The chemistry community uses the term calibration for this task, which alsointroduces some vocabulary overloading with our use of Bayesian model calibration.
LIBS. The authors also demonstrated the ability of ATOMICto replicate those effects from first principles. In particular,their experimental observations, supported by ATOMIC’s the-oretical calculations and modeling, showed that the sodiumlines increased significantly in emission intensity as more cop-per was added to the target. This effect is explained by anincrease in electron density, due to the copper, leading toincreased recombination within the plasma.In this paper we will also use simplified two-element targetsof sodium and copper to begin to explore our ability to useBayesian model calibration with ATOMIC to perform disag-gregation. As we look forward to more complex targets drawnfrom the large parameter space of all possible elements, andparticularly in the context of the high data-collection rate ofChemCam on Curiosity, we recognize that using ATOMICdirectly in this framework will be computationally intractable.We therefore also present our results building and evaluating emulators to provide fast approximations to the computation-ally expensive ATOMIC runs.Emulators are well-established tools in the context of slowcomputer models and their use for modeling spectra has beendemonstrated in the field of cosmology. Large N-body simula-tions of the universe are prohibitively expensive, and emulat-ing the matter power spectrum of the universe on cosmologicalscales was shown to be an effective solution [10]. Althoughthese spectra have very different physical origins from those ofLIBS, these prior results provide some context for our effortspresented here.Our team’s ongoing and preliminary work [1] includessuccess at estimating the plasma temperature and density ofsimple, fixed compounds, taking into account a structured dis-crepancy between ATOMIC simulations and measured LIBSdata. While Judge et al. [8] showed that ATOMIC can repli-cate the matrix effects in experimental data, the discrepancystudy in Bhat et al. [1] and our own explorations showed someas yet unresolved challenges when comparing ATOMIC andmeasured LIBS spectra in the quantitative way required to sup-port Bayesian model calibration. Therefore we focus here onthe methodology of performing disaggregation in the presenceof matrix effects in a simplified scenario, using a test set ofsimulated ATOMIC spectra rather than of measured spectraas would be our ultimate intent. This allows us to demonstratethe feasibility of using modular Bayesian model calibration toperform disaggregation with LIBS data.We present an overview of the ATOMIC simulations, inputs,and outputs in the following section. In Section 3 we outlinethe statistical framework of emulation and Bayesian calibrationused in this work, and we discuss their application to simula-tion data in Section 4. We conclude with further discussion ofthe relevance and context of our results and future directions. ébert
ET AL Computer model calibration relies on a small number of runsof a high-fidelity simulation, tiled across parameter space. Thisset of simulations can then be used to build emulators whichcan quickly approximate the computer model at new param-eter settings. In the context of the analysis of LIBS data, themodel of interest is the ATOMIC forward model, a generalpurpose plasma modeling and kinetics code. It was devel-oped to simulate the emission spectra of chemical compoundsusing first principles theoretical atomic physics [12], in par-ticular, emission or absorption spectra from plasmas either inlocal-thermodynamic equilibrium (LTE) or in non-LTE.ATOMIC simulations require a few primary inputs: the tem-perature and density of the plasma, and a model describing theatomic structure and scattering data of the material(s) consti-tuting the plasma. These last, which include quantities such asenergy levels and transition probabilities, are generated fromthe Los Alamos suite of atomic physics codes [6]. The resultsin the simulations discussed here were generated from theCATS code [5] with modifications made for plasmas generatedfrom LIBS [4]. For a given temperature and density, ATOMICthen models the emissivity of the plasma by computing itsaverage ionization.As motivated above, the simulations discussed in this paperare of sodium copper (NaCu) plasmas with varying ratios ofsodium (Na) to copper (Cu), inspired by Judge et al. [8]. Inaddition, we vary the input plasma temperature 𝑇 and the massdensity 𝜌 of the plasma. In particular, we run simulations overthe following space of these three parameters:1. plasma temperature 𝑇 , in units electron volts (eV) andrange [0 . , .
2. mass density 𝜌 , in units of g/cm and range [−7 , −4] onthe log scale3. composition of the plasma, defined as the proportion ofone of the two elements in our compound. We have arbi-trarily chosen to use the proportion of sodium in theplasma, % Na, in the range [0 , . % Cu can be retrievedby using
Na.A few comments are in order about the third input parame-ter, % Na. It is by solving the inverse problem for this parameterthat we are doing disaggregation — i.e., identifying the ele-ments present and estimating their proportions. The validityof this particular parameterization, where % Cu = 1 − % Na,holds only because our simulations are run in an artificiallysimplified scenario without atmosphere, components of which,such as carbon and oxygen, would usually be present in theplasma in unknown quantities. Thus we know a priori that theonly elements that could be present are sodium and copper. This enormously simplifies the first part of the disaggregationproblem: identifying which elements are present. Our thinkingis that success in this very constrained regime will establish afoundation for addressing the challenges faced by ChemCamon Mars, such as the presence of atmosphere and much largersets of candidate elements.The output of the ATOMIC simulation is a spectrum withintensity as a function of wavelength for a particular set ofthe above inputs, i.e. a spectrum for a particular compound atsome plasma density and temperature. The simulation providesintensity as power per volume per photon energy per unit solidangle. ATOMIC predictions do not account for any effects aris-ing from the spectrometer, but do include matrix effects. Thecomputation time for a single ATOMIC run depends on thechemical complexity of the compound of interest, and typ-ically varies from minutes to hours on a high-performancecomputing system.The ATOMIC simulations used in this analysis wereselected using two Latin hypercube designs over the three inputparameters: a training set of 500 simulations, and an indepen-dent test design of 25 points to which we added noise beforeanalysis, as described in Section 4.2. The training set parame-ter design, as well as a few of the resulting spectra, are shown inFigure 1 . We have labeled the peak locations for two sodiumpeaks and one copper peak along the top wavelength axes, andwill show these markers throughout this work when relevant.These are not the complete set of peaks for sodium or copperand are simply meant to provide some reference to the eye.We note, in the left panel of Figure 1 , the orders of mag-nitude over which the simulation output varies. This motivatesa rescaling to the log scale, as shown in the middle panel,before statistical analysis, due to anticipated difficulty captur-ing small, yet potentially important, details in the presenceof these large variations in amplitude. While this log trans-form is not standard practice in the spectroscopy community,exploratory analysis confirmed that emulators built for thespectra on their original scale performed significantly worsethan the results presented here. Additional motivations formodeling on the log scale are discussed in Bhat et al. [1].
The aim of this work is to demonstrate the application ofBayesian calibration methods to the problem of disaggregatingLIBS spectra — identifying the elements present and esti-mating their proportions. We start with a brief overview ofcomputer model calibration from Kennedy and O’Hagan [9].In this context, computer model calibration entails estimatingthe input parameters of a computer model (here, ATOMIC)that most likely generated some given observed experimental
Hébert
ET AL
FIGURE 1
Three examples of ATOMIC simulations used for training along with the design for the training set. Left panel:three example training spectra, with intensity in units of power per volume per photon energy per unit solid angle. Middle: Thesame three examples plotted as log intensity, thereby reducing the orders of magnitude variation across the data. Right panel:simulation input design for the 500 training spectra, scaled to the [0 , range for all parameters. The parameter settings for thethree examples shown in the previous panels are highlighted in color.or simulated data. A key assumption is that our observed data, 𝑦 , are a noisy version of the simulator output at some unknownparameter setting 𝜃 : 𝑦 = 𝜂 ( 𝜃 ) + 𝜖 (1)We will assume that the data 𝑦 have been centered and scaledaccording to the mean vector 𝜇 and scalar standard deviation 𝜎 of the training data.Following the well-established literature on calibration, wedenote the ATOMIC computer model as 𝜂 ( 𝑡 ) , which takes a 𝑝 -dimensional parameter vector 𝑡 as input to produce a LIBSspectrum. Here, 𝑝 = 3 for the three input parameters describedin Section 2. The vector 𝜃 represents the parameter valuesthat yield the model output 𝜂 that most closely resembles theobserved data. Note that we are working here with 𝑦 and 𝜂 ( 𝑡 ) as log scaled versions of the measured and modeled spectra.We estimate 𝜃 for a given observation 𝑦 by exploring theposterior 𝑝 ( 𝜃 | 𝑦 ) with Markov chain Monte Carlo (MCMC).This posterior is calculated, via Bayes rule, as the product ofthe data likelihood 𝑝 ( 𝑦 | 𝜃 ) and the parameter prior 𝑝 ( 𝜃 ) . Wewill first define the data likelihood, and the choice of prior isdiscussed below. Equation 1 describes how the data are gen-erated given the parameters 𝜃 . The noise 𝜖 determines thesampling function for the data. Here we take 𝜖 to be nor-mally distributed with mean and variance Σ 𝑦 , which givesthe following likelihood: 𝑦 | 𝜃 ∼ ( 𝜂 ( 𝜃 ) , Σ 𝑦 ) (2)This is in principle all we need (along with a prior) to per-form MCMC. There is a significant computational challenge, though, as the ATOMIC model 𝜂 is slow to compute: givensome vector 𝑡 , evaluating 𝜂 ( 𝑡 ) takes order of minutes or hours,rendering exploration via MCMC extremely slow. We will fol-low the standard approach to overcome this through use of anemulator: a statistical model that provides a fast approximationof the simulator output.In the following subsections we discuss first the approachof emulating the ATOMIC outputs using Gaussian processes,then details of Bayesian model calibration. As mentioned above, we define a statistical model to pro-vide fast approximations of the slow ATOMIC outputs. Thisemulator will be some unknown function conditioned on atraining set of 𝑚 simulator runs { 𝜂 ( 𝑡 ) , ..., 𝜂 ( 𝑡 𝑚 )} at fixedinputs 𝑡 , ..., 𝑡 𝑚 . For simplicity, we scale these inputs such that 𝑡 ∈ [0 , 𝑝 . Before we describe the details of the emula-tor, we address a second computational bottleneck due to thehigh-dimensionality of the data. Each simulator output 𝜂 has 𝑛 𝜂 = 32 , wavelength bins. Naive MCMC implementationsmight require the inversion of a , 𝑚 × 32 , 𝑚 matrixor , 𝑚 × 𝑚 matrices at each step to calculate the pos-terior. The solution to this problem, as developed in Higdonet al. [7], relies on using dimensionality reduction to find areduced set of basis vectors. This set of 𝑛 𝜂 -dimensional vec-tors { 𝑘 𝑖 , 𝑖 = 1 , ..., 𝑞 } describes the model for any input 𝑡 in thefollowing way: 𝜂 ( 𝑡 ) = 𝑞 ∑ 𝑖 =1 𝑘 𝑖 𝑤 𝑖 ( 𝑡 ) (3) ébert ET AL where the weights 𝑤 𝑖 ( 𝑡 ) hold the dependence on the input 𝑡 . Wedenote the number of components included in the reconstruc-tion by 𝑞 , with a maximum value of 𝑚 , the size of training set.Typically smaller values ∼ 10 suffice for good performance.This formalism reduces the computational complexity of theproblem enormously: we can now emulate and sample just 𝑞 ∼ 10 weights instead of 𝑛 𝜂 ∼ 10 wavelength bins.We find this new basis via a singular value decomposition(SVD) of the training simulation matrix 𝑋 . Each of the 𝑚 columns of 𝑋 holds a simulation of length , . The SVDfactorization of 𝑋 , in terms of orthogonal matrices 𝑈 , 𝑉 anddiagonal matrix 𝑆 of singular values of 𝑋 , can be written: 𝑋 = 𝑈 𝑆𝑉 𝑇 = 𝐾𝑊 (4)The second equality relates the SVD to Equation 3 in matrixform, with 𝐾 a column matrix of all the 𝑘 𝑖 s and the weights 𝑤 in row vector 𝑊 . We have defined 𝐾 = 𝑈 𝑆 ∕ √ 𝑚 and 𝑊 = 𝑉 𝑇 √ 𝑚 .In terms of the data described in Section 2, this decomposi-tion results in a set of principal components 𝐾 , each of whichhas length 𝑛 𝜂 = 32 , , that are common between all thespectra, and a set of weights which now hold the dependenceon the input parameters and which vary between each spec-trum. The first four of these principal components (PCs) areshown in Figure 2 along with the associated weights for the500 training and 25 test simulations.Our ATOMIC emulator will take the form of Equation 3,with 𝐾 calculated from the set of training simulations. Theweights 𝑤 𝑖 define a surface in parameter space for each com-ponent 𝑖 (see Figure 2 ): as is common in the literature, we willrepresent each of these by a Gaussian process (GP): 𝑤 𝑖 ( 𝑡 ) ∼ ( , 𝜎 𝑤𝑖 𝑅 𝑖 ( 𝑡 ) ) (5)where 𝜎 𝑤𝑖 is the marginal variance for weight 𝑖 and 𝑅 𝑖 ( 𝑡 ) isa correlation matrix with each entry given by the correlationfunction: corr ( 𝑡, 𝑡 ′ ; 𝑙 𝑖 ) = ∏ 𝑗 =1 exp ( − | 𝑡 𝑗 − 𝑡 ′ 𝑗 | 𝑙 𝑖𝑗 ) (6)where 𝑙 𝑖𝑗 is the length scale hyperparameter for weight 𝑖 andparameter 𝑗 .Whereas it is common to perform the estimation of the GPhyperparameters concurrently with the calibration, we treatthese in a modular way as in [11], fixing the hyperparam-eters by maximum likelihood estimation and keeping themfixed during the Bayesian calibration. This approach greatlyspeeds up the estimation process because we do not need torebuild and reinvert the covariance matrices at each step ofthe MCMC. During calibration, we use the GPs to predictthe value of each weight for parameter settings not present inthe training simulations. The prediction of weight 𝑖 for someparameters 𝜃 can be found using properties of conditional normal distributions: ̂𝑤 𝑖 | 𝜃 ∼ ( 𝑟 𝑖 ( 𝜃 ) 𝑅 −1 𝑖 𝑤 𝑖 , 𝜎 𝑤𝑖 [1 − 𝑟 𝑖 ( 𝜃 ) 𝑇 𝑅 −1 𝑖 𝑟 𝑖 ( 𝜃 )] ) (7)Here 𝑟 𝑖 ( 𝜃 ) is the 𝑚 × 1 vector found by applying the GPcovariance function in Equation 6 (with hyperparameters forweight 𝑖 ) to 𝜃 and the set of training simulation parameters { 𝑡 , ..., 𝑡 𝑚 } . The 𝑚 × 𝑚 matrix 𝑅 𝑖 is the correlation matrix ofthe training parameters, and 𝑤 𝑖 is the 𝑚 × 1 vector of trainingweights for principal component 𝑖 . We now describe the process of estimating the input parame-ters of the simulations via Bayesian model calibration. As seenin Equation 2, we assume our data 𝑦 is given by a simulatorrun 𝜂 ( 𝜃 ) , with some fixed, diagonal covariance which can beparameterized by a precision 𝜆 𝑦 : Σ 𝑦 = 𝜆 −1 𝑦 𝐼 . 𝑦 | 𝜂 ( 𝜃 ) ∼ ( 𝜂 ( 𝜃 ) , 𝜆 −1 𝑦 𝐼 ) (8)We have seen in the previous section how to express ourdata in a new basis of principal components 𝐾 and parameter-dependent weights 𝑊 . The emulator output 𝜂 ( 𝜃 ) , given param-eters 𝜃 , is expressed in the new basis using the weightspredicted by the GPs (see Equation 7): 𝜂 ( 𝜃 ) = 𝐾 ̂𝑤 ( 𝜃 ) (9)To cast a new observed spectrum 𝑦 into this basis, we use: 𝑤 𝑜𝑏𝑠 = ̃𝐾𝑦 (10)where ̃𝐾 = ( 𝐾 𝑇 𝐾 ) −1 𝐾 𝑇 , and we define 𝑤 𝑜𝑏𝑠 as the vector of 𝑞 weights corresponding to the observation 𝑦 .We use Equations 8, 10, and 9 and properties of Gaus-sian distributions to find the sampling of the observed weightsgiven ̂𝑤 ( 𝜃 ) , the GP predictions: 𝑤 𝑜𝑏𝑠 | ̂𝑤 ( 𝜃 ) ∼ ( ̃𝐾𝐾 ̂𝑤 ( 𝜃 ) , ̃𝐾𝜆 −1 𝑦 𝐼 ̃𝐾 𝑇 ) (11) = ( ̂𝑤 ( 𝜃 ) , ( 𝜆 𝑦 𝐾 𝑇 𝐾 ) −1 ) (12)Finally, since each predicted value ̂𝑤 𝑖 is drawn from a Gaussianprocess, these weights are normally distributed given parame-ters 𝜃 , as in Equation 7. This means that the observed weights 𝑤 𝑜𝑏𝑠 , given some parameter vector 𝜃 , are given by: 𝑤 𝑜𝑏𝑠 | 𝜃 ∼ ( 𝜇 𝑤 , ( 𝜆 𝑦 𝐾 𝑇 𝐾 ) −1 + Σ 𝑤 ) , (13)where 𝜇 𝑤 is a vector with entry 𝑖 given by 𝑟 𝑖 ( 𝜃 ) 𝑖 𝑅 −1 𝑤 𝑖 , 𝑖 =1 , ..., 𝑞 and Σ 𝑤 is a 𝑞 × 𝑞 diagonal matrix with element 𝑖𝑖 givenby 𝜎 𝑤𝑖 [1 − 𝑟 𝑖 ( 𝜃 ) 𝑇 𝑅 −1 𝑖 𝑟 𝑖 ( 𝜃 )] , where the 𝑖 denotes that these arecalculated from the covariance matrix with hyperparametersfor weight 𝑖 . Equation 13 is the likelihood we will use forMCMC exploration. Hébert
ET AL
FIGURE 2
Exploring the singular value decomposition of the LIBS simulations. The left panel shows the first four basiscomponents and the right panel shows the associated weights (in grey for each of the 500 training simulations; in black for 25noise-added test examples) as functions of the three input simulation parameters, which are each scaled to lie in [0 , . We applied the methods above using the Python Scikit-learnimplementation of Gaussian processes [14] and the MonteCarlo sampling in PyMC3 [15]. We chose a value of 𝑞 = 15 principal components in the emulator used to generate all theresults discussed in this section. The SVD reconstruction withthese 15 components explains > . of the variance ofthe training set. In addition, the emulator and calibration errorsdid not significantly improve with added components beyondthis point. Indeed, when a larger value is chosen for 𝑞 , theadditional weights added are less and less constraining forcalibration since their associated sampling variances (see thelikelihood in Equation 13) increase with weight index. Addingmore weights, then, is not expected to result in more accuratecalibration results.Both emulation performance and calibration results are cal-culated based on a test set of 25 ATOMIC simulations run ona Latin hypercube design independent from the training set, towhich we added noise. Figure 3 summarizes the emulator performance for the set-tings listed at the start of this section. The performance iscalculated on 25 test simulations.We quantify emulator performance by two metrics: 𝑅 andpercent error, both calculated point-wise for each wavelengthmodeled. We calculate these by: • 𝑅 = 1− 𝜎 𝑟𝑒𝑠 ∕ 𝜎 𝑟𝑎𝑤 , where 𝜎 𝑟𝑎𝑤 is the variance of the testsimulations around their mean, and 𝜎 𝑟𝑒𝑠 is the varianceof the residuals (emulator output - truth). • % error = 100× 𝜎 ( 𝜂 ( 𝜃 )− 𝑦 ) 𝜎𝑦 + 𝜇 , calculated for each test exampleand where 𝜇, 𝜎 are the mean vector and scalar standarddeviation of the training set.Overall, these errors lie between ±1% , with the exceptionof the test run plotted in yellow which we discuss below. Weconclude that the emulator is performing well.We expect that the information most important for our even-tual goal of disaggregation should lie near the peaks associatedwith constituent elements. For example, the locations of someimportant sodium and copper peaks are indicated in Figure 4.1by tick marks and vertical dashed lines. We see that the per-centage error is greater at peak locations than at most other ébert ET AL FIGURE 3
Summary of the emulator performance. Left: The 𝑅 and percent error as functions of wavelength (see Section 4.1for details). The percent error for each of the 25 test runs is overlaid in the bottom panel. The yellow curve shows the percenterror for the test run with the largest percent median absolute error, as shown in the right panel. Right: The test set design coloredby percent median (over all wavelengths) absolute error of the emulator prediction for each point. The training design is shownin grey.points along the wavelength axis. The dataset contains a widerange of sodium concentrations, implying high data varianceat these peak locations. In absolute terms, the percent errorshows that this variation is harder for the emulator to capturefully. However, the spikes in 𝑅 value at those locations are notsignificantly worse than at many other wavelengths, indicatingthat the emulator is performing well relative to the variance ofthe data. In summary, the emulator is able to capture impor-tant variations in the data which we expect to be indicative ofsodium concentration.The right panel in Figure 4.1 shows the test design, eachpoint colored according to the median percentage absoluteerror of the emulator predictions. The test run in yellow standsout as having the highest median error by far. Given its loca-tion on the very edge of parameter space for all three inputs, itis not surprising that the emulator has trouble with it. The per-centage error of this test run is highlighted in yellow in the lefthand panel, and we will discuss it further below.In summary, both the percentage errors and the 𝑅 values ofthe emulator predictions are indicative of strong performance,so we proceed with calibration. The Bayesian calibration was performed for noise-added ver-sions of the 25 test spectra, following Equation 8 with precision 𝜆 𝑦 = 4 . In addition to the 25 two-element NaCu compoundswith varying proportions of Na, we also ran the calibration for50 spectra of single-element simulations, with the composition for 25 fixed at and 25 fixed at Na. These single-element simulations were run with the same input design forplasma temperature and density as the 25 two-element testsimulations. During analysis these single-element spectra weretreated identically to the two-element test examples.The summary of these calibration results is shown in Figure4 . The left three panels of this figure show very encourag-ing calibration results for the two-target simulations, with thepoints colored according to the emulator errors as in the rightpanel of Figure 3 . The leftmost panel shows successful dis-aggregation for all 25 test spectra. Of particular interest is thetest run with highest emulator error that we discussed ear-lier, indicated in yellow. Recall that its true parameters are thehighest %Na in the design, relatively high temperature 𝑇 , andvery low density 𝜌 . We find that the calibration results for thischallenging test run are on par with the rest of the test set.Our single-element simulation results are summarized bythe histograms in the last panel of 4 . Our success in recog-nizing pure Na and pure Cu demonstrates the ability of theemulator trained on the two-element NaCu simulations to cor-rectly identify examples of single elements. We had had someconcern that the difference between a target with a vanishinglysmall amount of one element and one with none at all wouldlook like a discontinuity of some sort that would be difficult tocapture with our emulators, and these results suggest that ouremulators are robust to that difference.Figure 5 shows bivariate marginal distributions of theMCMC samples alongside the log scaled test spectrum for twotwo-element test examples. Results for the example with high Hébert
ET AL
FIGURE 4
Summary of calibration results. The left three panels show the posterior mean of 15,000 MCMC samples for25 noise-added test simulations, using the same color scheme as in the right panel of Figure 3 to indicate emulator errors.Panels show results for estimated %Na, temperature 𝑇 , and density 𝜌 , compared to the true parameter value. Right: histogramof estimated %Na from single target simulations (see Section 4.2), for sodium-only targets in the upper panel and copper-onlyin the lower panel.emulator error discussed previously are displayed in the toppanel, while a randomly selected example was chosen for thelower panel. Both spectra and true parameter values are shownin colors corresponding to their median absolute percent errorsdiscussed in Section 4.1. All three parameter combinationsshow some evidence of correlations between the samples, mostnotably between density and temperature in the upper panel. In this work we show promising first results toward the goalof disaggregating LIBS data from ChemCam with Bayesianmodel calibration. The emulation and Bayesian calibrationmethods discussed here successfully perform disaggregationin the presence of matrix effects in noise-added simulations oftwo-element compounds.An important next step is to test this approach on measuredLIBS data. There are a number of challenges to overcomebefore this can succeed. While ATOMIC simulations comparevery well to experimental data on many metrics and for otherapplications [3, 4, 8], the wavelength-by-wavelength matchingthat drives the emulation in this work is not one of them. Forexample, we have found that the peak shapes, widths, and tosome extent even locations, differ between ATOMIC outputsand measured LIBS spectra. These small effects have not beenstraightforward to accommodate with the methods presentedhere.In addition, LIBS spectra measured “in the wild” on Mars oreven in a laboratory setting typically contain extraneous peaksthat don’t provide information about the target of interest. For instance, when measured in air, the plasmas include contribu-tions from the atmosphere in unknown proportions that appearas peaks in the spectra. And laboratory targets of even simpletwo-element compounds like NaCu typically are formed usingbinders such as stearic acid that introduce still more peakswhen measured via LIBS. These extraneous components of theplasmas exacerbate matrix effects and complicate our analy-sis. A possible path forward that we are exploring is to limitthe set of wavelengths in our analyses to those around expert-identified peaks of interest, rather than considering the entirespectrum shown in the work here.In addition, instrument response and other effects shouldbe accounted for when comparing simulations to measure-ments. Standard practice for ChemCam data is to correct forsome of these effects, and we have begun exploring how thesecorrections impact our ability to perform disaggregation.While we continue those explorations toward comparingsimulations with measurements as a separate line of research,we will build on our successes with simulated data that wepresented here in order to provide a proof of principle fordisaggregating increasingly complex multi-element targets.
ACKNOWLEDGMENTS
Research presented in this paper was supported by theLaboratory Directed Research and Development programof Los Alamos National Laboratory under project number20180097ER. CAH was supported by the Department ofEnergy Computational Science Graduate Fellowship (DE-FG02-97ER25308). ébert
ET AL FIGURE 5
Example bivariate marginal distribution for twotwo-element test examples. 15,000 MCMC samples are shownfor each. The log-scale test spectrum and true parameter valuesare shown, colored to match the associated emulator error fromFigure 3 ; for example, the upper triangle plot shows samplesfor the test example which demonstrated the highest emulatorerrors. The vertical range is the same for both spectra shownhere, with ticks at log intensity of 20 and 22.5.
Conflict of interest
The authors declare no potential conflict of interests.
References [1] Bhat, K. S., K. Myers, E. Lawrence, J. Colgan, andE. Judge, 2020: Structured Discrepancy in BayesianModel Calibration for ChemCam on the Mars CuriosityRover. in press .[2] Clegg, S. M., E. Sklute, M. D. Dyar, J. E. Barefield, andR. C. Wiens, 2009: Multivariate analysis of remote laser-induced breakdown spectroscopy spectra using partialleast squares, principal component analysis, and relatedtechniques.
Spectrochimica Acta Part B: Atomic Spec-troscopy , , no. 1, 79–88.[3] Colgan, J., E. Judge, H. Johns, D. Kilcrease, J. Bare-field II, R. McInroy, P. Hakel, R. Wiens, and S. Clegg,2015: Theoretical modeling and analysis of the emis-sion spectra of a ChemCam standard: Basalt BIR-1A. Spectrochimica Acta Part B: Atomic Spectroscopy , ,20–30.[4] Colgan, J., E. Judge, D. Kilcrease, and J. Barefield II,2014: Ab-initio modeling of an iron laser-inducedplasma: Comparison between theoretical and experimen-tal atomic emission spectra. Spectrochimica Acta Part B:Atomic Spectroscopy , , 65–73.[5] Cowan, R. D., 1981: The theory of atomic structure andspectra . Univ of California Press.[6] Fontes, C., H. Zhang, J. Abdallah Jr, R. Clark, D. Kil-crease, J. Colgan, R. Cunningham, P. Hakel, N. Magee,and M. Sherrill, 2015: The Los Alamos suite of relativis-tic atomic physics codes.
Journal of Physics B: Atomic,Molecular and Optical Physics , , no. 14, 144014.[7] Higdon, D., J. Gattiker, B. Williams, and M. Right-ley, 2008: Computer model calibration usinghigh-dimensional output. Journal of the Ameri-can Statistical Association , , no. 482, 570–583,doi:10.1198/016214507000000888.[8] Judge, E. J., J. Colgan, K. Campbell, J. E. Bare-field, H. M. Johns, D. P. Kilcrease, and S. Clegg,2016: Theoretical and experimental investigationof matrix effects observed in emission spectra ofbinary mixtures of sodium and copper and magne-sium and copper pressed powders. SpectrochimicaActa Part B: Atomic Spectroscopy , , 142 – 148,doi:https://doi.org/10.1016/j.sab.2016.06.004.[9] Kennedy, M. and A. O’Hagan, 2001: Bayesian Calibra-tion of Computer Models. Journal of the Royal StatisticalSociety. Series B (Statistical Methodology) , , no. 3,425–464. Hébert
ET AL [10] Lawrence, E., K. Heitmann, J. Kwan, A. Upadhye,D. Bingham, S. Habib, D. Higdon, A. Pope, H. Finkel,and N. Frontiere, 2017: The mira-titan universe. II. matterpower spectrum emulation.
The Astrophysical Journal , , no. 1, 50, doi:10.3847/1538-4357/aa86a9.URL https://doi.org/10.3847%2F1538-4357%2Faa86a9[11] Liu, F., M. Bayarri, J. Berger, et al., 2009: Modulariza-tion in Bayesian analysis, with emphasis on analysis ofcomputer models. Bayesian Analysis , , no. 1, 119–150.[12] Magee, N. H., J. Abdallah, J. Colgan, P. Hakel, D. Kil-crease, S. Mazevet, M. Sherrill, C. J. Fontes, andH. Zhang, 2004: Los Alamos opacities: Transition fromLEDCOP to ATOMIC. AIP Conference Proceedings ,London, UK, volume 730, 168âĂŞ179.[13] NASA, 2019: Curiosity Rover Mission Overview. https://mars.nasa.gov/msl/mission/overview/ .[14] Pedregosa, F., G. Varoquaux, A. Gramfort, V. Michel,B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer,R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cour-napeau, M. Brucher, M. Perrot, and E. Duchesnay, 2011:Scikit-learn: Machine learning in Python.
Journal ofMachine Learning Research , , 2825–2830.[15] Salvatier, J., T. V. Wiecki, and C. Fonnesbeck, 2016:Probabilistic programming in Python using PyMC3. PeerJ Computer Science , 2:e55, doi:10.7717/peerj-cs.55.[16] Walters, D. J., A. Biswas, E. C. Lawrence, D. C. Fran-com, D. J. Luscher, D. A. Fredenburg, K. R. Moran, C. M.Sweeney, R. L. Sandberg, J. P. Ahrens, and C. A. Bolme,2018: Bayesian calibration of strength parameters usinghydrocode simulations of symmetric impact shock exper-iments of al-5083.
Journal of Applied Physics ,124