# Scalable multicomponent spectral analysis
for high-throughput data annotation

R. Patrick Xian1,2∗ Ralph Ernstorfer1 Philipp M. Pelz3,4

1Fritz Haber Institute of the Max Planck Society, 14195 Berlin, Germany.
2Department of Neurobiology, Northwestern University, Evanston 60208, IL, USA.
3National Center for Electron Microscopy Facility, Molecular Foundry,
Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA.
4Department of Materials Science and Engineering,
University of California, Berkeley, CA 94720, USA.
∗Correspondence authors<EMAIL_ADDRESS>

> Orchestrating parametric fitting of multicomponent spectra at scale is an
> essential yet underappreciated task in high-throughput quantification of
> materials and chemical composition. To automate the annotation process for
> spectroscopic and diffraction data collected in counts of hundreds to
> thousands, we present a systematic approach compatible with high-performance
> computing infrastructures using the MapReduce model and task-based
> parallelization. We implement the approach in software and demonstrate
> linear computational scaling with respect to spectral components using
> multidimensional experimental materials characterization datasets from
> photoemission spectroscopy and powder electron diffraction as benchmarks.
> Our approach enables efficient generation of high-quality data annotation
> and online spectral analysis and is applicable to a variety of analytical
> techniques in materials science and chemistry as a building block for
> closed-loop experimental systems.

## I. Introduction

Real-time understanding of experimental data in materials and chemical
characterization is a pursuit for the coming age of automation [1]. Building
machine learning models for experimental techniques that produce spectral or
spectrum-like data requires tools for spectrum annotation, which pairs data
with quantitative information extracted by fitting approximate lineshape or
structural models [2]. Achieving computational efficiency in these tasks
allows processing large amounts of data with negligible human intervention as
well as execution of online or on-the-fly analysis [3, 4]. In materials
characterization measurements, changes in sample and environmental conditions
could appear in minutes to hours, while experimental campaigns often run
continuously for days or longer. Moreover, the data acquisition process often
involves tuning of external variables [5] such as the spatial locations, time,
temperature, concentration, etc., resulting in multidimensional datasets of
multicomponent spectra (each component being a spectral lineshape or
background model). As a realistic example, fitting one thousand 10-component
spectra, with 3-4 parameters for each component would result in a large-scale
optimization problem involving on the order of $10^{4}$ parameters. To assess
the experimental outcome on this scale is not uncommon in modern high-
throughput experimentation [6]. The dataset shifts during the experiment also
make it challenging to apply modern machine learning frameworks without
sufficient human interaction [3, 4]. Therefore, data annotation methods with
efficient computational scaling with respect to spectral complexity,
quantified by the number of components, becomes crucial for automating
materials characterization.

In this work, we leverage high-performance computing (HPC) to achieve linear
scaling for a class of spectrum annotation task – multicomponent spectral
analysis (MSA) [7] or fitting. Broadly speaking, the MSA problem encompasses
two complementary scenarios: (1) Fitting an unknown number of peaks at fixed
spectral locations with varying amplitudes; (2) Fitting a fixed number of
peaks at varying spectral locations. Ultimately, MSA extracts information such
as the shape and height of spectral components and disentangles signal from
background. Mathematically, MSA may be mapped to separable nonlinear least
squares problems [8] present in various branches of science and engineering.
In analytical chemistry, where scenario (1) is frequently encountered for
understanding spectroscopic data, matrix factorization techniques without [9]
or with limited prior knowledge [10] of pure spectral components (i.e. spectra
from single sources or species) are often adopted for data analysis. In
comparison, scenario (2) remains a recurring challenge in upscaling the
analysis pipeline of quantitative electron- and X ray-based characterization
techniques, when large, annotated, empirical databases are yet to exist for
training machine learning models. Moreover, the equivalents in scenario (2) of
pure spectral components are generally unattainable. Currently, relevant
softwares for materials characterization usually only provide the option for
sequentially fitting individual spectra [11] or pattern in the case of
diffraction [12]. Therefore, we have developed pesfit [13], a software to
tackle the MSA problem in (2) at scale to facilitate online and offline
analysis and experimental feedback. In the following, we discuss the software
architecture, validate the computational scaling on batch (or distributed)
spectrum fitting and present use cases in materials characterization and
discuss the insights gained from data.

## II. Results

### I. Scalable software architecture

Fitting a multicomponent spectrum requires iterative calculation of its
components, which are typically modelled as unimodal probability distributions
or analytical functions. The mathematical context is provided in Appendix B.
In practice, a complex spectrum may comprise tens of overlapping components or
may be divided into segments of a similar size. Maintaining a reasonable
computational performance with a large number of spectra with increasing
complexities requires both resource and architectural upgrade on almost all
steps of single-spectrum fitting that remains the mainstay in existing data
analysis workflows [11, 12].

Figure 1: Schematic of the scalable software architecture. Each fitting task
for a multicomponent spectrum (left) requires iterative evaluation of an
approximate spectral lineshape model. A MapReduce approach implemented in
MultipeakModel (see Appendix B) is used to execute the calculation of all
components. During optimization, the fitting parameters are updated to refine
the outcome, guided by a least squares error function involving the data and
the fit. The entire pool of fitting tasks (middle) formed by aggregating
individual ones is distributed over the multicore CPUs in a high-performance
computing infrastructure (right) and different task schedulers may be used.
The cartoon shows 16 worker processes each assigned with 6 fitting tasks.
Task-based parallelization with the asynchronous scheme shown here (right)
allows to balance the computational loads between fitting tasks.

We have realized linear scaling in pesfit to satisfy the computational needs
and alleviate the manual task load in (1) multicomponent spectrum computation,
(2) batch execution, as well as the associated steps of (3) initial condition
generation and (4) fitting result collection. To minimize redundancy in (1),
we constructed a generating class (MultipeakModel, see Appendix B) to
represent multicompoent spectra, making use of a MapReduce structure [14] for
evaluating the spectra given the constituents (see Fig. 1 left). Our approach
extends the popular spectral lineshape fitting framework lmfit [15] to allow
batch fitting a multitude of spectra with an arbitrary number of components.
The MapReduce approach splits up (“map” step) the evaluation of multicomponent
spectra into parallelizable computation of the components followed by merging
(“reduce” step) of the results (see Fig. 1 left) and allows to generate an
arbitrary complex spectral shape from existing lineshape models or user-
defined functions without hard-coding each type of multicomponent spectrum
separately. Each task thus constructed executes fitting of a single spectrum
(see Eq. (2) in Appendix B). To enable (2), the compartmentalized fitting
tasks (see Fig. 1 middle) are executed concurrently through dask [16], parmap
[17], and torcpy [18] on an HPC infrastructure with multicore CPUs (see Fig. 1
right). Resource sharing (including the numerical data and the multicomponent
spectrum model) between tasks is implemented in distributed fitting
(DistributedFitter, see Appendix B) to reduce the memory footprint. To
simplify (3), the pesfit package includes convenient methods to
programmatically generate initial conditions for a large number of fitting
parameters (see Appendix B). The initial conditions are assembled into a
nested key-value pairs (Python dictionary) to deliver to the optimizer. In
step (4), the fitting results in the form of parameter lists are collected by
a DataFrame data structure from pandas, whose expansive functionality allows
further statistical analysis of the outcome [19]. The fitting results and the
initial conditions may be exported for reproducible examination of the
specific fitting tasks. A schematic of the computational workflow is shown in
SM Fig. 1.

### B. Performance evaluation with use cases

We demonstrate the use of batch fitting on experimental data from two popular
materials characterization techniques: photoemission spectroscopy and powder
electron diffraction. Both of these techniques generate multidimensional
datasets involving multicomponent spectra or spectrum-like signals. Spectrum
fitting is carried out with predefined, domain-specific multicomponent models
(see Appendix C).

Figure 2: Computational evaluation on experimental datasets. Performance
comparison between the sequential fitting procedure and the distributed
fitting procedure on benchmark datasets: a, photoemission spectroscopy (PES)
of the valence band (around K point in momentum space) of WSe2 with up to 14
energy bands considered; b, Time-resolved powder electron diffraction (PED) of
polycrystalline Pt [20], with up to 10 diffraction peaks and a background
profile considered. The computing times ($T_{\mathrm{sequential}}$ and
$T_{\mathrm{distributed}}$), here used as the evaluation metric, are
normalized by both the number of components and spectra (or diffraction
peaks).

In momentum-resolved photoemission spectroscopy (PES), extracting electronic
band dispersion from the 3D dataset requires batch fitting of momentum-
dependent intensity-valued energy spectra. We use photoemission band mapping
data of bulk tungsten diselenide (WSe2), which was measured for about an hour
and contains a maximum of 14 energy bands within the probed energy range of
$\sim$ 8 eV, to benchmark the fitting algorithms. The electron diffraction
dataset was recorded for polycrystalline platinum (Pt) with a pump-probe
scheme – pump light induces electronic and lattice dynamics that are probed by
electron diffraction – to obtain time- and fluence-dependent 1D diffraction
patterns [20], resulting in 3D dataset. The data was acquired continuously
over 2 days with up to 10 diffraction peaks of Pt within the measured
reciprocal space. Further details of the experiments are provided in Appendix
A.

In both use cases, the batch fitting is initiated with realistic guesses in
their respective domains (see Appendix C), the optimization processes using
both sequential and distributed approaches are timed for comparison. We
compare the computing time (see Fig. 2) that estimates the duration between
the start of the fitting until the results are collected, which accounts for
the accumulated wall time of steps (1), (2), and (4) mentioned earlier. The
computing times are normalized to the per-spectrum and per-component level to
allow comparison between datasets from different experimental modalities. As
shown in Fig. 2, within the tested cases, the distributed approach generally
exhibits a linear dependence on the number of components in the data, while
the sequential approach, as is common practice in the field, has at least a
quadratic scaling. The disparity in normalized computing time between the two
use cases reflects characteristics of the fitting tasks (spectrum model,
initialization, hyperparameters, etc). Nevertheless, in both cases,
distributed fitting of a significant portion of the data is sufficiently
faster than the respective experimental duration, further discussed in
Appendix C. The polynomial speed-up therefore allows multicomponent spectral
analysis on the fly during the experimental measurements, with only the added
cost on the computing hardware. Moreover, the tuning of the initialization to
improve the fitting outcome may be conducted interactively using the same
software, which allows to design custom approaches to deal with complex
spectra. Some intuitive approaches to achieve high performance and consistency
within the fitting results are discussed in the SM and schematized in SM Fig.
1.

Specifically for the PES data, SM Fig. 2 shows that the batch fitting method
provides rapid access to the local dispersion in an extended region ($\sim$
0.6 Å $\times$ 0.6 Å covered by 1600 energy spectra) around K point of the
surface Brillouin zone of WSe2, therefore uncovers the different strength of
the trigonal warping (TW) feature in the energy bands visible by eye. TW is a
third-order angular modulation [21] of the ideal parabolic dispersion
(representing isotropic dispersion without angular distortion) often depicted
for semiconductors. It is a manifestation of the strong spin-orbit interaction
within the material [22] and has influence over the transport properties [23,
24, 25]. Although theoretical descriptions exist [22, 21], so far it has not
yet been densely reconstructed directly from experimental data. The much
improved throughput of our method now allow us to further study the details of
TW in other energy bands of WSe2 and compare with theoretical calculations.

## III. Discussion and conclusion

We have described a linear-scaling software architecture for mapping the
multicomponent spectral analysis problem to distributed hardware. Our approach
fully exploits the HPC infrastructure to provide a convenient toolkit for
speeding up data annotation in materials characterization without requiring
repeated data loading and human intervention. In the use case demonstrations,
although we have assumed prior knowledge for the number of components, the
requirement may be achieved using numerical estimation [26, 27]. The modular
architecture of the associated software, pesfit, allows convenient upgrade of
its components individually in the future. For example, it is possible to plug
in alternative optimizers for batch fitting, such as brute-force search or
variable projection [28], which have respective speed-accuracy trade-offs in
solving large-scale nonlinear least squares problems. The computational
performance demonstrated here provides a baseline for future algorithm
development using such annotations for training supervised learning
algorithms. Besides the use cases discussed here, we envision the method and
software to be applicable to techniques such as transient vibrational
spectroscopy [29], X-ray fluorescence imaging [30], to name a few, where
spectral shift along with the changes in spectral width and amplitude have
unique physical and chemical implications, and where multicomponent spectral
analysis are adopted for analyzing the measured multidimensional datasets. The
software may be integrated with existing computational workflows [12, 31],
thereby functioning as a building block for closed-loop materials or chemical
characterization and analytics [3, 4] to facilitate experimental optimization
and discovery.

### Appendix A: Materials characterization datasets

Both experimental datasets are measured with home-built instruments at the
Fritz Haber Institute. The photoemission spectroscopy experiment is conducted
using pulsed light source based on high harmonic generation as described
previously [32]. The sample under study is single crystalline tungsten
diselenide (WSe2). Extreme UV pulses with an energy of 21.7 eV liberate
electrons from the WSe2 sample surface, which are subsequently captured as
single events in 3D (energy and two momenta) by a time-of-flight delay-line
detector (METIS 1000, SPECS GmbH). A custom computational pipeline is used to
preprocess the single events to produce the 3D photoemission band mapping data
[31]. The time-resolved powder electron diffraction experiment [20] contains
fluence-dependent measurements of polycrystalline Pt diffraction pattern with
the laser excitation at 0.70 eV and 70 kV electron bunches generated from a
gold photocathode [33]. The 2D diffraction patterns are recorded with a CMOS
camera (TemCam-F416, TVIPS GmbH) at a frame integration time of 5 s. The
measured diffraction rings are collapsed into 1D diffraction peaks via radial
averaging.

### Appendix B: Software details and mathematical context

The main software modules in pesfit are lineshape, fitter and metrics. The
lineshape module contains the MultipeakModel class. The fitter module contains
funtions to automate the construction of fitting tasks and initializations as
well as simple visualization methods for displaying fitting results. The
metrics module includes evaluation metrics for assessing the quality of fits.
The sequential and distributed fitting are achieved through the PatchFitter
and DistributedFitter classes, respectively, both lie within the fitter
module. To fit higher-dimensional ($>$ 2D) spectral data, the data are first
partially transformed to 2D before the fitting tasks are farmed out to
processors. The initial conditions can be generated programmatically using
init_generator and supplied to a multicomponent model via varsetter, both are
integrated into the fitter module classes discussed here, but are also
invocable separately and individually. The ordering of fitting tasks is
tracked by a spectrum index (spec_id) to accommodate asynchronous concurrent
execution, which permits index scrambling.

Both sequential and distributed fitting of the benchmark datasets uses the
default Levenberg-Marquardt algorithm for minimization with a spectrum-wise
least-squares loss function and bound constraints. We formulate the single
multicomponent spectral analysis task as

$\tilde{\textbf{S}}=\mathrm{argmin}\,\left\|I(\omega)-\sum_{n}f_{n}(\omega,\textbf{s}_{n})\right\|_{\omega}^{2}.$
(1)

Here, $I$ is the intensity-valued spectrum data with the spectral dimension
$\omega$. $\textbf{s}_{n}=(s_{n,1},s_{n,2},...)$ is the set of parameters of
the single-component spectral lineshape or background model $f_{n}$, with $n$
being the index of spectral components. The symbol
$\left\|\cdot\right\|_{\omega}$ denotes the Euclidean norm with respect to the
spectral dimension.
$\tilde{\textbf{S}}=\left[\tilde{\textbf{s}}_{1}\,\tilde{\textbf{s}}_{2}\,...\right]$
is the concatenated list of component-wise fitting parameters after
optimization. For example, a 10-component spectrum with 4 parameters for each
component results in an $\tilde{\textbf{S}}$ of size 40 for a single-spectrum
fitting task. Based on this notation, batch fitting of multidimensional data
is written as

$\\{\tilde{\textbf{S}}^{\\{a_{i}\\}}\\}=\mathrm{argmin}\,\left\|I(\omega,\\{a_{i}\\})-\sum_{n}f_{n}^{\\{a_{i}\\}}(\omega,\textbf{s}_{n}^{\\{a_{i}\\}})\right\|_{\omega}^{2},$
(2)

where the spectrum indices $i$ are positive integers. Now, $I$ becomes
multidimensional with the coordinates of other dimensions represented by
$\\{a_{i}\\}=(a_{1},a_{2},...)$. A scalar-valued spectrum index is assigned to
each multidimensional coordinate $\\{a_{i}\\}$ according to the row-major
order. On the left-hand side of Eq. (2),
$\\{\tilde{\textbf{S}}^{\\{a_{i}\\}}\\}=\\{\left[\tilde{\textbf{s}}_{1}\,\tilde{\textbf{s}}_{2}\,...\right]^{\\{a_{i}\\}}\\}$
involves a double concatenation (firstly row-wise over the inner brackets
$\left[\,\cdots\right]$ and ordered by the spectral components, secondly
column-wise over the outer brackets $\\{\,\cdots\\}$ and ordered by the
spectrum index), which represents the entirety of spectrum-wise optimized
parameter sets gathered over all $\\{a_{i}\\}$ coordinates. This arrangement
turns the fitting results into a 2D array regardless of the dimensionality of
the data. The nonlinear least-squares problems in Eqs. (1)-(2) are considered
separable in the sense that the spectrum model may be written as a sum of
components [8, 28].

### Appendix C: Computational benchmarks

All batch fitting benchmark use cases are run on an on-premises computing
server (Dell PowerEdge R840), equipped with four Intel Xeon Gold 6150
multicore CPUs. The sequential fitting tasks are run with a single logical
core, while the distributed fitting tasks are run using the designated number
of processes mapped onto logical cores. Throughout the runs, the RAM size
doesn’t pose a limit on the computing time of the benchmark fitting tasks.
Fitting of the photoemission data of WSe2 uses the density functional theory
electronic structure calculation available within the NOMAD repository (DOI:
10.17172/NOMAD/2020.03.28-1) along with a rigid global shift to all band
energies as the initial condition for the band positions (in both sequential
and distributed approaches), while fitting the diffraction peak positions of
polycrystalline Pt uses the estimate positions from a reference spectrum [20].
The trends of the computational scaling of the two approaches persist
regardless of the number of spectra used for benchmarking. For photoemission
(see Fig. 2a), the scaling was estimated by fitting 900 spectra. Each spectrum
was fit with up to 14 Voigt lineshapes in combination. For experimental
monitoring, band dispersion information from a region around a high-symmetry
point of the material suffices. As a time reference, batch fitting of 1600
4-component photoemission spectra (see SM Fig. 2) costs around 6 mins.

The electron diffraction data was fit with up to 10 Voigt lineshapes and an
exponential-decay background model. A total of 225 1D diffraction patterns
were obtained by averaging from the over 4600 patterns measured in the 2-day
experimental campaign. For computational scaling estimation (see Fig. 2b),
fitting the entire set of 225 11-component 1D diffraction patterns costs $<$ 3
mins. Since in practice, only a fraction of the data (both in range and in
quantity) is needed to monitor the sample condition, the speed boost is
sufficient to enable real-time interaction with experiments. Running the
distributed fitting consistently at the optimal speed requires tuning of the
task parallelization parameters – the number of worker processes and the
number of tasks per process. An extended example that illustrate the
interdependence between these two parameters is shown in SM Fig. 3.
Nevertheless, in practice, only a few trials guided by heuristics on subsets
of data are needed to find a reasonable trade-off that minimizes the computing
time.

### Acknowledgments

We thank S. Dong, S. Beaulieu and L. Rettig for providing the photoemission
spectroscopy data, D. Zahn for providing the electron diffraction data and
helpful discussions. We thank S. Schülke and G. Schnapka at Gemeinsames
Netzwerkzentrum (GNZ) and L. Rettig at the Fritz Haber Insitute in Berlin for
support on the computing infrastructure. R.P.X. thanks S. Oller at
Universitätsklinikum Hamburg-Eppendorf (UKE) for helpful discussion on code
parallelization. R.P.X. and R.E. acknowledge the support by BiGmax, the Max
Planck Society’s Research Network on Big-Data-Driven Materials-Science, the
European Research Council (ERC) under the European Union’s Horizon 2020
research and innovation program (Grant No. ERC-2015-CoG-682843). P.M.P.
acknowledges financial support from STROBE: A National Science Foundation
Science & Technology Center under Grant No. DMR 1548924.

### Data Availability

The data presented in the paper is available upon reasonable request from the
corresponding author.

### Author contributions

R.P.X. conceived the project and developed the software with help from P.M.P..
R.P.X. wrote the first draft of the paper. All authors contributed to
discussion and revision of the manuscript to its final form.

### Competing interests

The authors declare no competing interests in the content of the article.

> Supplementary Material
>  Scalable multicomponent spectral analysis for high-throughput data
> annotation

R. Patrick Xian1,2∗ Ralph Ernstorfer1 Philipp M. Pelz3,4

1Fritz Haber Institute of the Max Planck Society, 14195 Berlin, Germany.
2Department of Neurobiology, Northwestern University, Evanston 60208, IL, USA.
3National Center for Electron Microscopy Facility, Molecular Foundry,
Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA.
4Department of Materials Science and Engineering,
University of California, Berkeley, CA 94720, USA.
∗Correspondence authors<EMAIL_ADDRESS>

##  S1 Modes of operation

Although the problem mapping and software realization presented in this work
offer a scalable and automatable solution for batch fitting of multicomponent
spectra, situations will always exist when a subset of the fitting tasks yield
unfeasible answers due to experimental resolution, noise in the data, the
sloppiness of the model used for fitting [34], etc. This also applies to
offline data analysis where a more accurate model is often required for
fitting to extract specific parameters, while in online analysis, approximate
models with less computational overhead are favored. In the following, we
offer here practical solutions outlined in the computational workflow (see SM
Fig. 1) to resolve potential outstanding cases while making use of the
existing functionalities offered by pesfit. In the demonstrations below, we
discuss primarily in the context of photoemission data due to its complexity
compared with diffraction data, but a similar line of reasoning also applies
to other materials characterization methods that produce spectra or spectrum-
like data.

###  S1.1 Tuning of hyperparameters and initialization

To arrive at a meaningful outcome, batch fitting of multicomponent spectra
requires suitable initialization with knowledge from theoretical calculations
or educated guesses. We treat rigid shifts applied to theory or reference in
initialization as hyperparameters in the fitting procedure, which may be tuned
according to the data characteristics. We discuss here three possible ways
that may be realized using pesfit.

SM Figure 1: Proposed multicomponent spectral analysis workflow. The
computational workflow starts from the dataset organized into a
multidimensional array. Then, the multicomponent spectrum model is constructed
using domain knowledge. If theoretical calculation or a reference spectrum is
available, the workflow progresses to the stage of tuning task parallelization
parameters. Otherwise, an initial guess fit of single spectra at distinct
locations is required, along with interpolation of the spectral peak locations
to other instances without guess fitting results. Using theory or reference
spectrum positions as input, hyperparameters tuning is needed to compensate
for the offsets from experimental data, but this is often unnecessary if guess
fitting results are available. Finally, the batch fitting routine may be
deployed for offline analysis or online monitoring of experimental results,
using user-defined metrics constructed with the spectrum fitting outcome.

Theory-based approach. When theoretical calculations or reference spectra are
present, (i) The random_varshift method in the fitter module of pesfit enables
automated tuning of initial condition using random search within user-defined
values to optimize a goodness-of-fit metric. (ii) Alternatively, if segments
of the spectra may be isolated with no overlap with the rest (e.g. when
spectral intensity goes to baseline level at the two ends of the segment), a
(1) divide-and-conquer approach may be adopted to batch fit and tune the
hyperparameters for a segment at a time. This applies to diffraction data and
the photoemission data of materials with relatively flat and isolated bands,
such as organic molecular crystals or monolayers [35]. If the isolation
between segments is less perfect,

SM Figure 2: Semi-automated fine-tuning of hyperparameters. Tuning of the
hyperparameters in batch fitting multicomponent spectra with a, the divide-
and-conquer approach and b, the expanding window approach, illustrated with
synthetic spectra. The two segments are numbered in dashed boxes with
different colors. Demonstration of the tuning uses c. the photoemission
dataset (containing 1600 multicomponent spectra) measured near the
$\overline{\mathrm{K}}$ point of WSe2. Other important landmarks within the
projected first Brillouin zone (1BZ) include the zone center
($\overline{\Gamma}$) and zone edge ($\overline{\mathrm{K}}$). The overline
signifies the surface projection (e.g. $\overline{\mathrm{K}}$ in
photoemission corresponds to K). d. Comparison of the band dispersions
reconstructed using batch fitting without (upper row) and with band-wise fine-
tuning (middle row) of the shift hyperparameters of band structure
calculation. The lower row shows the batch fitting results using data-derived
initialization. The fourth band is not resolved here due to strong overlap
with the third band at around K. The band mapping data used for fitting has
been symmetrized to balance the matrix elements. The colormap for each energy
band is scaled separately for visualization.

one may adopt an (2) expanding window approach to batch fit and tune an
increasing number of spectral components (thus the “expanding window”),
retaining the hyperparameters used for the previous segments while fitting in
an increasingly broader spectral window. These two approaches are illustrated
in SM Fig. 3a-b. During the tuning process, the fitting outcome may be
appraised using either goodness-of-fit metrics or via visual inspection. An
example showing the results before and after hyperparameter tuning is shown in
SM Fig. 3c-d using photoemission data from WSe2. The results compared include
batch fitting initialized with a 0.2 eV rigid shift of all bands (3c) and that
initialized with a band-wise tuned rigid shift of 0.2 eV, 0.18 eV, 0.22 eV and
-0.24 eV for each energy band (3d), respectively.

Data-driven approach. When theoretical calculations or reference spectra are
not available, (iii) one can resort to a data-driven approach. A selected
number of distinctive spectra may be fit manually (guess fit) to determine the
anchor points as reference for interpolation, which generates the initial
conditions for the entire dataset. The data-derived initial conditions may
then be used in batch fitting. An example using 16 equally-spaced anchor
points to fit 1600 (40 $\times$ 40 grid in $k_{x}$ – $k_{y}$ coordinates)
multicomponent spectra is shown in SM Fig. 3d (lower row). Here, we used the
Clough-Tocher interpolator [36] implemented in scipy [37] to obtain the data-
derived initial conditions for all bands, although other multivariate
nonlinear interpolators should suffice in general. The approach is realized
with the InteractiveFitter class within pesfit. In comparison, the data-driven
approach shows higher consistency among neighboring spectra due to the closer
resemblance of data-derived initial condition than the theoretical calculation
to the final outcome.

###  S1.2 Tuning of model complexity

In cases where strictly isolating a group of spectral features using data
slicing is challenging, there may exist the situation where a subset of the
data requires a more complex model (secondary model) to fit than the rest
(primary model). Alternatively, the tuning may also be realized by changing
the data range used in fitting. This appears commonly in cases where spectral
shift is present, such as in photoemission data due to the existence of band
dispersion and matrix element effects, which results in the dependence of
spectral features on momentum coordinate (in momentum-resolved measurements),
sample position (in spatially-resolved measurements), etc. Comparatively, in
diffraction datasets (such as in [38]) involving phase transition due to
changes in macroscopic parameters or dynamics, the diffraction peak position
changes, if any, are relatively small with respect to the peak width. After
batch fitting with the primary model, the less well-fit spectra may be
selected using goodness-of-fit metrics. These outlier spectra may be batch fit
again using a secondary spectrum model with more (or less) components than the
primary model.

SM Figure 3: Effects of task parallelization parameters on computing time.
Linear-scale heatmap of the computing time for the distributed multicomponent
fitting procedure presented in this work. Each entry represents the computing
time (in seconds) elapsed in a numerical experiment, where a distinct pair of
parameters (processes and tasks per process) is selected to execute the same
number of multicomponent spectral analysis tasks (900 spectra in total, each
with 4 components). The benchmark uses the photoemission data from around the
K point of WSe2. The purple boxes represent the possible cases to try in the
performance tuning guided by heuristics.

###  S1.3 Tuning of task parallelization parameters

For any batch fitting task, there are generally two key parameters to tune to
minimize the overall computing time by balancing the number of worker
processes ($n_{\mathrm{W}}$) and the number of tasks per worker process
($n_{\mathrm{T}}$), or chunk size in this context [39]. Their trade-off
depends on the complexity of the multicomponent spectrum model and the quality
of the initial condition provided to the optimizer. We illustrate the trade-
off in the example in SM Fig. 3. To facilitate the use of our approach, we
provide here a list of potential choices for these parameters for reference.
From our experience, batch fitting hundreds to thousands of spectra generally
requires $n_{\mathrm{W}}>10$ for performance optimization, although for low-
complexity problems (e.g. with very few spectral components), a smaller number
of processes will suffice. it’s generally. when there are only 2-3 components
in the multicomponent spectrum model, $n_{\mathrm{T}}$ should be a few
hundreds to maximize the computational efficiency. When there are 4-6
components per spectrum, $n_{\mathrm{T}}$ should be below 100. When there are
8-10 components per spectrum, $n_{\mathrm{T}}$ should be less than 30-40.
Beyond 10 components per spectrum, it is advised not to use $n_{\mathrm{T}}$
larger than 10-20. As a rule of thumb, the product of $n_{\mathrm{W}}$ and
$n_{\mathrm{T}}$ should be equal or slightly larger than the total number of
fitting tasks. We illustrate the trial cases (as purple boxes) in the example
in SM Fig. 3 for tuning the parallelization parameters using the heuristics
discussed here.

## References

  * [1] Steven R. Spurgeon et al.  “Towards data-driven next-generation transmission electron microscopy”  In _Nature Materials_ , 2020  DOI: 10.1038/s41563-020-00833-z
  * [2] Phil Woodruff  “Modern Techniques of Surface Science”  Cambridge: Cambridge University Press, 2016  DOI: 10.1017/CBO9781139149716
  * [3] A. Kusne et al.  “On-the-fly closed-loop materials discovery via Bayesian active learning”  In _Nature Communications_ 11.1, 2020, pp. 5966  DOI: 10.1038/s41467-020-19597-w
  * [4] Peter M. Attia et al.  “Closed-loop optimization of fast-charging protocols for batteries with machine learning”  In _Nature_ 578.7795, 2020, pp. 397–402  DOI: 10.1038/s41586-020-1994-5
  * [5] Philip Hofmann  “In operando devices studied by angle-resolved photoemission spectroscopy”  In _arXiv_ , 2020, pp. 2011.11490  arXiv: http://arxiv.org/abs/2011.11490
  * [6] Mitsutaro Umehara et al.  “Analyzing machine learning models to accelerate generation of fundamental materials insights”  In _npj Computational Materials_ 5.1, 2019, pp. 34  DOI: 10.1038/s41524-019-0172-5
  * [7] J.. Blackburn  “Computer Program for Multicomponent Spectrum Analysis Using Least Squares Method.”  In _Analytical Chemistry_ 37.8, 1965, pp. 1000–1003  DOI: 10.1021/ac60227a013
  * [8] Axel Ruhe and Per Åke Wedin  “Algorithms for Separable Nonlinear Least Squares Problems”  In _SIAM Review_ 22.3, 1980, pp. 318–337  DOI: 10.1137/1022057
  * [9] Anna Juan, Joaquim Jaumot and Romà Tauler  “Multivariate Curve Resolution (MCR). Solving the mixture analysis problem”  In _Anal. Methods_ 6.14, 2014, pp. 4964–4976  DOI: 10.1039/C4AY00571F
  * [10] E. Kriesten et al.  “Identification of unknown pure component spectra by indirect hard modeling”  In _Chemometrics and Intelligent Laboratory Systems_ 93.2, 2008, pp. 108–119  DOI: 10.1016/j.chemolab.2008.05.002
  * [11] Giorgio Speranza and Roberto Canteri  “RxpsG a new open project for Photoelectron and Electron Spectroscopy data processing”  In _SoftwareX_ 10, 2019, pp. 100282  DOI: 10.1016/j.softx.2019.100282
  * [12] Laurent P. René de Cotret, Martin R. Otto, Mark J. Stern and Bradley J. Siwick  “An open-source software ecosystem for the interactive exploration of ultrafast electron scattering data”  In _Advanced Structural and Chemical Imaging_ 4.1, 2018, pp. 11  DOI: 10.1186/s40679-018-0060-y
  * [13] Rui Patrick Xian  “pesfit”, https://github.com/mpes-kit/pesfit URL: https://github.com/mpes-kit/pesfit
  * [14] Ralf Lämmel  “Google’s MapReduce programming model — Revisited”  In _Science of Computer Programming_ 70.1, 2008, pp. 1–30  DOI: 10.1016/j.scico.2007.07.001
  * [15] Matt Newville et al.  “lmfit/lmfit-py 1.0.0”  In _Zenodo_ , https://doi.org/10.5281/zenodo.3588521, 2019  DOI: 10.5281/ZENODO.3588521
  * [16] Dask Development Team  “Dask: Library for dynamic task scheduling”, 2016  URL: https://dask.org
  * [17] Sergio Oller  “parmap”, https://github.com/zeehio/parmap URL: https://github.com/zeehio/parmap
  * [18] P.E. Hadjidoukas et al.  “torcpy: Supporting task parallelism in Python”  In _SoftwareX_ 12, 2020, pp. 100517  DOI: 10.1016/j.softx.2020.100517
  * [19] Wes McKinney  “Data Structures for Statistical Computing in Python”  In _Proceedings of the 9th Python in Science Conference_ , 2010, pp. 56–61  DOI: 10.25080/Majora-92bf1922-00a
  * [20] Daniela Zahn, Hélène Seiler, Yoav William Windsor and Ralph Ernstorfer  “Ultrafast lattice dynamics and electron-phonon coupling in laser-excited platinum”  In _arXiv_ , 2020, pp. 2012.10428  arXiv: http://arxiv.org/abs/2012.10428
  * [21] Andor Kormányos et al.  “k · p theory for two-dimensional transition metal dichalcogenide semiconductors”  In _2D Materials_ 2.2, 2015, pp. 022001  DOI: 10.1088/2053-1583/2/2/022001
  * [22] Andor Kormányos et al.  “Monolayer MoS ¡math display=”inline”¿ ¡msub¿ ¡mrow/¿ ¡mn¿2¡/mn¿ ¡/msub¿ ¡/math¿ : Trigonal warping, the ¡math display=”inline”¿ ¡mi¿$\Gamma$¡/mi¿ ¡/math¿ valley, and spin-orbit coupling effects”  In _Physical Review B_ 88.4, 2013, pp. 045416  DOI: 10.1103/PhysRevB.88.045416
  * [23] Hongyi Yu et al.  “Nonlinear Valley and Spin Currents from Fermi Pocket Anisotropy in 2D Crystals”  In _Physical Review Letters_ 113.15, 2014, pp. 156603  DOI: 10.1103/PhysRevLett.113.156603
  * [24] Wen-Yu Shan, Jianhui Zhou and Di Xiao  “Optical generation and detection of pure valley current in monolayer transition-metal dichalcogenides”  In _Physical Review B_ 91.3, 2015, pp. 035402  DOI: 10.1103/PhysRevB.91.035402
  * [25] A.. Kalameitsev, V.. Kovalev and I.. Savenko  “Valley Acoustoelectric Effect”  In _Physical Review Letters_ 122.25, 2019, pp. 256801  DOI: 10.1103/PhysRevLett.122.256801
  * [26] Robert R. Lindner et al.  “Autonomous Gaussian decomposition”  In _The Astronomical Journal_ 149.4, 2015, pp. 138  DOI: 10.1088/0004-6256/149/4/138
  * [27] Sebastian E. Ament et al.  “Multi-component background learning automates signal detection for spectroscopic data”  In _npj Computational Materials_ 5.1, 2019, pp. 77  DOI: 10.1038/s41524-019-0213-0
  * [28] Gene Golub and Victor Pereyra  “Separable nonlinear least squares: the variable projection method and its applications”  In _Inverse Problems_ 19.2, 2003, pp. R1–R26  DOI: 10.1088/0266-5611/19/2/201
  * [29] Michael P. Grubb et al.  “Translational, rotational and vibrational relaxation dynamics of a solute molecule in a non-interacting solvent”  In _Nature Chemistry_ 8.11, 2016, pp. 1042–1046  DOI: 10.1038/nchem.2570
  * [30] M. Pushie et al.  “Elemental and Chemically Specific X-ray Fluorescence Imaging of Biological Systems”  In _Chemical Reviews_ 114.17, 2014, pp. 8499–8541  DOI: 10.1021/cr4007297
  * [31] R. Xian et al.  “An open-source, end-to-end workflow for multidimensional photoemission spectroscopy”  In _Scientific Data_ 7.1, 2020, pp. 442  DOI: 10.1038/s41597-020-00769-8
  * [32] M. Puppin et al.  “Time- and angle-resolved photoemission spectroscopy of solids in the extreme ultraviolet at 500 kHz repetition rate”  In _Review of Scientific Instruments_ 90.2, 2019, pp. 023104  DOI: 10.1063/1.5081938
  * [33] Lutz Waldecker, Roman Bertoni and Ralph Ernstorfer  “Compact femtosecond electron diffractometer with 100 keV electron bunches approaching the single-electron pulse duration limit”  In _Journal of Applied Physics_ 117.4, 2015, pp. 044903  DOI: 10.1063/1.4906786
  * [34] Mark K. Transtrum et al.  “Perspective: Sloppiness and emergent theories in physics, biology, and beyond”  In _The Journal of Chemical Physics_ 143.1, 2015, pp. 010901  DOI: 10.1063/1.4923066
  * [35] Yasuo Nakayama, Satoshi Kera and Nobuo Ueno  “Photoelectron spectroscopy on single crystals of organic semiconductors: experimental electronic band structure for optoelectronic properties”  In _Journal of Materials Chemistry C_ 8.27, 2020, pp. 9090–9132  DOI: 10.1039/D0TC00891E
  * [36] Peter Alfeld  “A trivariate clough—tocher scheme for tetrahedral data”  In _Computer Aided Geometric Design_ 1.2, 1984, pp. 169–181  DOI: 10.1016/0167-8396(84)90029-3
  * [37] Pauli Virtanen et al.  “SciPy 1.0: fundamental algorithms for scientific computing in Python”  In _Nature Methods_ 17.3, 2020, pp. 261–272  DOI: 10.1038/s41592-019-0686-2
  * [38] Sanjay Singh et al.  “Robust Bain distortion in the premartensite phase of a platinum-substituted Ni2MnGa magnetic shape memory alloy”  In _Nature Communications_ 8.1, 2017, pp. 1006  DOI: 10.1038/s41467-017-00883-z
  * [39] Navtej Singh, Lisa-Marie Browne and Ray Butler  “Parallel astronomical data processing with Python: Recipes for multicore machines”  In _Astronomy and Computing_ 2, 2013, pp. 1–10  DOI: 10.1016/j.ascom.2013.04.002