# Communication-avoiding micro-architecture to compute Xcorr scores for
peptide identification

Sumesh Kumar, Fahad Saeed  Knight Foundation School of Computing and
Information Sciences
Florida International University (FIU)
Miami, FL USA 33199
{sumesh.kumar<EMAIL_ADDRESS>
###### Abstract

Database algorithms play a crucial part in systems biology studies by
identifying proteins from mass spectrometry data. Many of these database
search algorithms incur huge computational costs by computing similarity
scores for each pair of sparse experimental spectrum and candidate theoretical
spectrum vectors. Modern MS instrumentation techniques which are capable of
generating high-resolution spectrometry data require comparison against an
enormous search space, further emphasizing the need of efficient accelerators.
Recent research has shown that the overall cost of scoring, and deducing
peptides is dominated by the communication costs between different hierarchies
of memory and processing units. However, these communication costs are seldom
considered in accelerator-based architectures leading to inefficient DRAM
accesses, and poor data-utilization due to irregular memory access patterns.
In this paper, we propose a novel communication-avoiding micro-architecture to
compute cross-correlation based similarity score by utilizing efficient local
cache, and peptide pre-fetching to minimize DRAM accesses, and a custom-
designed peptide broadcast bus to allow input reuse. An efficient bus
arbitration scheme was designed, and implemented to minimize synchronization
cost and exploit parallelism of processing elements. Our simulation results
show that the proposed micro-architecture performs on average 24x better than
a CPU implementation running on a 3.6 GHz Intel i7-4970 processor with 16GB
memory.

###### Index Terms:

cross-correlation, protein identification, SEQUEST, accelerator, micro-
architecture

††publicationid: pubid: 978-1-5386-5541-2/18/$31.00 ©2018 IEEE

##  I Introduction

Mass-spectrometry based analysis has been the preferred method for
identification of proteins from complex biological samples [1]. The last two
decades have seen tremendous developments in data acquisition and analysis
techniques which have enabled many powerful proteomic applications. Database
search algorithms such as SEQUEST[7], X!Tandem[6], and MSFragger [9] can now
search high resolution mass-spectrometry data against an ever increasing
protein database to produce high quality matches. This has drastically
increased the compute load for existing implementations of the database search
algorithms. In this regard, several studies have used parallelization
strategies using high-performing compute clusters[10, 18, 16], GPUs[13, 11,
3], and FPGAs[5, 4, 19, 15] to speed up the computation process. However, a
recent study suggests that the major bottleneck in mass-spectrometry based
analysis is the cost of communication i.e cost of moving input and output data
between different hierarchies of a system[17]. Thus, even though CPUs are
operating at a much higher frequency, their performance gain for proteomics
studies relies on efficiently utilizing system cache or some other input reuse
technique [2] to minimize the number of DRAM accesses. Consequently, the
implementation of Crux [12], state-of-the-art software for computing cross-
correlation (Xcorr) scores, utilizes processor registers to store peptide
fragment ions to allow peptide reuse. While this allows one-side data reuse,
the cost of accessing experimental spectra from main memory is not minimized
as generally CPU registers are not large enough to hold the entire
experimental spectrum. On the other hand custom architectures using FPGAs can
achieve better performance for memory bound applications by utilizing the
abundant on-chip RAM resources and custom-designed communication minimizing
pipelines to allow experimental spectrum reuse[14].

In this paper, we propose a communication-avoiding micro-architecture to
accelerate the Xcorr score computation which achieves two-side data reuse by
utilizing the on-chip RAM to cache an entire experimental spectrum and a
peptide broadcast bus to decrease the number of DRAM accesses. Our experiments
show that these optimizations result in 600x reduction in the average number
of DRAM accesses compared with a no-caching approach and 24x times speed-up
over Crux. The main contributions of this paper are as follows:

  1. 1.

We implemented a block RAM based cache of size 2kB to store experimental
spectra and minimize redundant DRAM accesses.

  2. 2.

We pre-sorted the peptide database which allows the use of binary search to
search candidate peptides. The search operation needs to be performed only
once per spectrum as next peptide can be pre-fetched, hence achieving input
locality.

  3. 3.

To allow input reuse, we designed a peptide broadcast bus to make it
accessible to all the processing elements.

  4. 4.

We implemented a first-come first-serve (FCFS) based bus arbitration scheme to
minimize the synchronization time of processing elements sharing the system
bus.

The rest of the paper is organized as follows. Section II explains the
background of Xcorr computation problem and related work. Section III
describes the proposed architecture. Section IV presents the experimental
results. Section V concludes the paper.

##  II Background

###  II-A Xcorr theoretical formulation

The Xcorr score between a theoretical spectrum vector $\bm{X}$ and an
experimental spectrum vector $\bm{Y}$ of length $n$ is defined in [8] as,

$X_{corr}=\sum_{i=0}^{n-1}X[i]Y[i]-\frac{1}{151}\sum_{i=0}^{n-1}\sum_{\tau=-75}^{\tau=75}X[i]Y[i-\tau]$
(1)

where $\tau$ is the amount by which vector is being serially shifted. However,
SEQUEST implementation performs an optimization by pre-processing the
experimental spectrum to perform dot product only once as described in[7] and
summarized below,

$\bm{Y_{P}}=\sum_{i=0}^{n-1}\left(Y[i]-\frac{1}{151}\sum_{\tau=-75}^{\tau=75}Y[i-\tau]\right)$
(2)

using (2) reduces the Xcorr computation to

$X_{corr}=\sum_{i=0}^{n-1}X[i]Y_{P}[i]$

###  II-B Related work

A significant amount of work has been done on acceleration of peptide
deduction algorithms using FPGA based architectures. Bogdan, Coca and
Beynon[5] designed a FPGA based accelerator for peptide deduction using
Profound[20] algorithm by instantiating 48 parallel search processors to
achieve 950$\times$ speed-up over a single core processor. Another study
accelerated the scoring process of X!Tandem[15] by instantiating 6 score
generation modules and one fragment ion generation module to achieve
17$\times$ speed-up over a CPU only implementation. In [19], a complete CPU-
FPGA system which was based on[15] but included support for multiple FPGAs and
achieved 10-fold speed-up for the entire search operation. To the best of our
knowledge, there hasn’t been any studies performed on accelerating the
SEQUEST[3] algorithm, which is widely used in proteomics.

##  III Proposed Architecture

The architectural setting of the heterogeneous computational system for Xcorr
is shown in Fig. 1. Host CPU communicates with PCIe DMA via the PCIe link to
transfer the experimental spectra from host memory to FPGA memory. A set of
directly accessible core control registers, hold computation parameters, and
control the operation. Each step of the algorithm takes place inside the
processing element (PE) i.e. reading experimental spectrum vectors one by one,
searching for candidate peptides, generating theoretical spectrum, computing
dot product scores, and writing the results back to main memory. The system
allows deployment of multiple PEs which execute the computations in a parallel
and asynchronous manner. Since all the PEs share the same memory bus, we
implemented a first come first serve (FCFS) based bus arbitration scheme to
achieve maximum bandwidth utilization. The detailed view of the bus
arbitration scheme is described in Fig. 3.

Figure 1: Complete system architecture shows host CPU communicates with FPGA
RAM via PCIe DMA bridge which is connected to Intel’s Avalon memory mapped
bus. Core registers module contains the computation parameters and is also
used for FPGA-CPU communication. To allow efficient use of the Avalon memory
mapped bus, all PEs are connected to FCFS based bus arbiter which is in turn
connected to Avalon memory mapped bus.

###  III-A Processing Element construction

Each PE takes over the computation of a single spectrum with all the candidate
peptides. At the heart of a PE, sits a controller which determines the flow of
computation as shown in Fig. 2. The controller copies an experimental spectrum
in the form of m/z and ion intensity values from the external memory into on-
chip RAM and starts the computation. Once the scores have been computed, they
are collected in the on-chip RAM and a request for bus access is generated
again to copy the scores into the DRAM.

Figure 2: Detailed internal construction of a single processing element. At
the heart is the control logic which controls the function of all the sub-
modules in the figure. Binary search module fetches a candidate peptide and
stores it in a peptide FIFO. Ion generator reads the peptide and generates
fragment ions. A 512 bit packet containing 16 32-bit ion mz and intensity pair
values along with a 32-bit theoretical ion and intensity pair values are fed
to the ion-matching kernel which finds the matching peak and stores the
partial score in on-chip RAM.

###  III-B Bus Arbiter Design

To ensure load-balancing among the PEs, we designed a bus-arbiter shown in
Fig. 3 which aimed to minimize total wait time for all the PEs. All the PEs
requiring access to the bus connect with the ”bus request” signal which is
connected to a wait counter register. The wait counter register keeps track of
the wait time of every bus master so that the decision of contest for bus
access is based on fairness i.e. access is granted to a master which has been
wait for the longest.

Figure 3: Brief logic description of bus arbitration module. Bus request lines
from all the PEs are coming into the arbiter. When a PE is denied service, its
wait count register is incremented, dynamically increasing its priority for
the next turn. Find max module is a comparator tree which finds which
registers has the maximum value and grants access to the corresponding PE.

###  III-C Ion-Matching Kernel

To compute the dot product scores between experimental spectra and a candidate
peptide, the processing element moves a 64-byte word from the on-chip RAM and
a theoretically generated ion-pair from the candidate peptide to the peak-
matching circuit. Each 64-byte word has 16 ion-pairs (using 16-bit floating
point representation for intensity and 16 bit binary representation for m/z)
from the experimental spectrum which is stored in 16 32-bit registers inside
the peak-matching circuit. The m/z value of the theoretical ion is compared
with all the experimental m/z values using a set of 16 parallel comparators as
shown in Fig. 4 and the corresponding matching intensity value of the peak is
multiplied and accumulated in the output register. Once all the ions are
traversed, the final Xcorr score is sent to the local on-chip RAM and the
process repeats for the next candidate peptide.

Figure 4: The ion-matching circuit receives a 512 bit packet containing 16
experimental ions which are all compared with a theoretical ion in one cycle.
The matched ions are multiplied and accumulated in the score register. If the
theoretical ion is outside the range of current experimental ions, next packet
is requested from the onchip-RAM by incrementing the counter.

##  IV Experiments

###  IV-A Methodology

We designed the entire hardware using Intel Quartus Pro and Qsys system
builder for Intel Stratix 10 FPGA. VHDL description was compiled using Quartus
Pro to verify the maximum operable frequency of 200MHz. To evaluate the timing
performance of our design, we implemented a cycle accurate simulator in python
which mimicked the exact timing response of the hardware. In our simulator, we
modeled each sub-module as a class whose data objects represented the internal
and external signals of the module and a clock-event() method which updated
the signals whenever a clock edge occurred.

For our experiments, we used the PXD000612 dataset from PRIDE database which
contained 90494 experimental spectra to score against human proteome dataset
containing 669964 peptides. The experimental spectra were stored in the
compressed sparse row (CSR) format with ion m/z value as the data index and
ion intensity value as data element. Ion m/z values were stored in a 16 bit
binary format and ion intensity values were represented using 16 bit half-
precision floating point format.

(a)

(b)

(c)

(d)

Figure 5: (a), Total computation time vs number of instantiated PEs for cache
size of 512B, 1kB, 2kB, and 4kB is shown. The search time decreases
consistently with increasing number PEs for cache size 2kB and 4kB but its
saturates for 512B and 1kB after 6 PEs due to increased memory requests. (b)
Average synchronization time per PE in the system. This is the time spent by
PE waiting for memory bus access. The wait time increases exponentially for
cache size below 2kB. (c) Average computation time per PE in the system. This
is the time spent on actual computation which decreases with the increasing
number of PEs but is not affected by cache size. (d) Average I/O time per PE
is also higher for cache size below 2kB. However there is little difference in
performance between 2kB and 4kB.

###  IV-B Results

The performance gains in our design come from a combination of optimizations
which minimize DRAM accesses and allow input reuse by using an on-chip RAM as
a local cache. To find the optimal cache size, we performed a design space
exploration for four different cache sizes along with the number of
instantiated PEs in the design. The results of these experiments are presented
in Fig. 5.

We analyzed the performance of our design by elaborating the total processing
time spent on computation and communication. To understand the effect of cache
size, we further divide the communication time in terms of I/O and waiting
time. We define these terms as Average computation time: average time each PE
spends on computing dot product. Average I/O time: average time spent by each
PE on DRAM read/write operations. Average waiting time: average time each PE
spends on waiting to get access to system bus. The total processing time for
dot product computations is shown in Fig. 5a, it is evident that increasing
the number of PEs from 1 to 31 displayed significant speed up until 15
processing elements for a cache size of 2kB and 4kB, while for cache size
below 2kB the speed up is plateaued after 6 processing elements.

Fig. 5b, 5c, and 5d show the breakdown of the total processing time in terms
of computation, I/O and waiting time for a single PE. Fig. 5b shows that
waiting time is zero for 1 PE as memory bus is not being shared. As the number
of PEs are increased, there is an almost exponential increase in the average
waiting time of a PE for cache-sizes below 2kB. The waiting times for 2kB and
4kB cache stay constant even when 31 PEs are instantiated. The average
computation time per PE is not impacted by the size of cache, but it decreases
sharply when processing elements are increased to 11. Fig. 5d shows that the
total I/O time i.e. total number of DRAM accesses are orders of magnitude
greater for cache size below 2kB. Table I further illustrates that increasing
cache-size from 1kB to 2kB results in 600$\times$ reduction in the average I/O
time.

TABLE I: Effect of cache size on average I/O and synchronization time while using 16 PEs Cache | I/O | Waiting | Total communication
---|---|---|---
size | time | time | time
512B | 1.01s | 11.57s | 12.58s
1kB | 0.52s | 5.19s | 5.71s
2kB | 0.86ms | 2.2ms | 3.06ms
4kB | 0.84ms | 2.1ms | 2.95ms

Based on our experiments, in our final design we instantiated 16 processing
elements to achieve maximum performance from the system. We compared the total
search time of our design with Crux[12] for 6 different values of precursor
mass window. Table II presents the total run-time of Crux running on a 3.6GHz
Intel i7-4970 processor with 16GB of system memory and the run-time of our
proposed hardware accelerator running at 200MHz clock frequency.

TABLE II: Run-time Comparison of FPGA accelerator with Crux running on 3.6GHz Intel i7-4970 using 8 threads and 16GB memory. Precursor mass | Dot product |  | Hardware | Relative
---|---|---|---|---
Tolerance (Da) | operations | Crux | Accelerator | Speed-up
1.5 | 162.79M | 53.2s | 1.25s | 42
3 | 325.41M | 75s | 2.45s | 30
5 | 541.42M | 86s | 4.10s | 21
10 | 1.07B | 139s | 7.782s | 20
25 | 1.99B | 304s | 20.75s | 15
50 | 3.076B | 648s | 39.6s | 16

##  V Conclusion

In this paper we designed, and developed an efficient communication-avoiding
micro-architecture. By using extensive experimentation, we demonstrated the
applicability of custom hardware design approach to accelerate crucial memory
bound problems in MS based omics. We presented optimizations for input reuse
at all stages of the computation including cache implementation, pre-fetching,
and input broadcasting. Although the system was designed for SEQUEST, it can
easily be applied for other scoring techniques which involve dot product
computation with little modification. Our simulation results suggest that our
design is scalable for up-to 32 PEs with linear speed-ups. In future, we plan
to extend this work to include the complete peptide identification process
involving protein digestion, and build towards a general purpose proteomics
processor.

## Acknowledgment

Research reported in this paper was supported by NIGMS of the National
Institutes of Health under award number: R01GM134384. Fahad Saeed was further
supported by the National Science Foundations (NSF) under the Award Numbers
NSF CAREER OAC-1925960. The content is solely the responsibility of the
authors and does not necessarily represent the official views of the National
Institutes of Health or the National Science Foundation.

## References

  * [1] Ruedi Aebersold and Matthias Mann.  Mass spectrometry-based proteomics.  Nature, 422(6928):198–207, 2003.
  * [2] Shuichi Asano, Tsutomu Maruyama, and Yoshiki Yamaguchi.  Performance comparison of fpga, gpu and cpu in image processing.  In 2009 international conference on field programmable logic and applications, pages 126–131. IEEE, 2009.
  * [3] Lydia Ashleigh Baumgardner, Avinash Kumar Shanmugam, Henry Lam, Jimmy K Eng, and Daniel B Martin.  Fast parallel tandem mass spectral library searching using gpu hardware acceleration.  Journal of proteome research, 10(6):2882–2888, 2011.
  * [4] Fabiola Casasopra, Gea Bianchi, Gianluca C Durelli, and Marco D Santambrogio.  Parallel protein identification using an fpga-based solution.  In 2016 IEEE International Parallel and Distributed Processing Symposium Workshops (IPDPSW), pages 295–299. IEEE, 2016.
  * [5] Daniel Coca, Istvan Bogdan, and Robert J Beynon.  Parallel fpga search engine for protein identification.  Bioinformatics: High Performance Parallel Computer Architectures, page 313, 2010.
  * [6] Robertson Craig and Ronald C Beavis.  Tandem: matching proteins with tandem mass spectra.  Bioinformatics, 20(9):1466–1467, 2004.
  * [7] Jimmy K Eng, Bernd Fischer, Jonas Grossmann, and Michael J MacCoss.  A fast sequest cross correlation algorithm.  Journal of proteome research, 7(10):4598–4602, 2008.
  * [8] Jimmy K Eng, Ashley L McCormack, and John R Yates.  An approach to correlate tandem mass spectral data of peptides with amino acid sequences in a protein database.  Journal of the american society for mass spectrometry, 5(11):976–989, 1994.
  * [9] Andy T Kong, Felipe V Leprevost, Dmitry M Avtonomov, Dattatreya Mellacheruvu, and Alexey I Nesvizhskii.  Msfragger: ultrafast and comprehensive peptide identification in mass spectrometry–based proteomics.  Nature methods, 14(5):513–520, 2017.
  * [10] Chuang Li, Kenli Li, Tao Chen, Yunping Zhu, and Qiang He.  Sw-tandem: a highly efficient tool for large-scale peptide identification with parallel spectrum dot product on sunway taihulight.  Bioinformatics, 35(19):3861–3863, 2019.
  * [11] You Li, Hao Chi, Leihao Xia, and Xiaowen Chu.  Accelerating the scoring module of mass spectrometry-based peptide identification using gpus.  BMC bioinformatics, 15(1):1–11, 2014.
  * [12] Sean McIlwain, Kaipo Tamura, Attila Kertesz-Farkas, Charles E Grant, Benjamin Diament, Barbara Frewen, J Jeffry Howbert, Michael R Hoopmann, Lukas Käll, Jimmy K Eng, et al.  Crux: rapid open source protein tandem mass spectrometry analysis.  Journal of proteome research, 13(10):4488–4491, 2014.
  * [13] Jeffrey A Milloy, Brendan K Faherty, and Scott A Gerber.  Tempest: Gpu-cpu computing for high-throughput database spectral matching.  Journal of proteome research, 11(7):3581–3591, 2012.
  * [14] Eriko Nurvitadhi, Jaewoong Sim, David Sheffield, Asit Mishra, Srivatsan Krishnan, and Debbie Marr.  Accelerating recurrent neural networks in analytics servers: Comparison of fpga, cpu, gpu, and asic.  In 2016 26th International Conference on Field Programmable Logic and Applications (FPL), pages 1–4. IEEE, 2016.
  * [15] Jin Qiu, Ping Kang, Li Ding, Yipeng Yuan, Wenbo Yin, and Lingli Wang.  Fpga acceleration of the scoring process of x! tandem for protein identification.  In 2017 27th International Conference on Field Programmable Logic and Applications (FPL), pages 1–4. IEEE, 2017.
  * [16] Rovshan G Sadygov, Jimmy Eng, Eberhard Durr, Anita Saraf, Hayes McDonald, Michael J MacCoss, and John R Yates.  Code developments to improve the efficiency of automated ms/ms spectra interpretation.  Journal of proteome research, 1(3):211–215, 2002.
  * [17] Fahad Saeed.  Communication lower-bounds for distributed-memory computations for mass spectrometry based omics data, 2020.  unpublished.
  * [18] Leheng Wang, Wenping Wang, Hao Chi, Yanjie Wu, You Li, Yan Fu, Chen Zhou, Ruixiang Sun, Haipeng Wang, Chao Liu, et al.  An efficient parallelization of phosphorylated peptide and protein identification.  Rapid Communications in Mass Spectrometry, 24(12):1791–1798, 2010\.
  * [19] Moucheng Yang, Tao Chen, Xuegong Zhou, Liang Zhao, Yunping Zhu, and Lingli Wang.  A complete cpu-fpga architecture for protein identification with tandem mass spectrometry.  In 2019 International Conference on Field-Programmable Technology (ICFPT), pages 295–298. IEEE, 2019.
  * [20] Wenzhu Zhang and Brian T Chait.  Profound: an expert system for protein identification using mass spectrometric peptide mapping information.  Analytical chemistry, 72(11):2482–2489, 2000.


# A Roadmap to Decipher Ultrafast Photophysics in Two-Dimensional
Nanomaterials

Pieter Schiettecatte  Physics and Chemistry of Nanostructures, Department of
Chemistry, Ghent University, 9000 Ghent, Belgium
<EMAIL_ADDRESS>Zeger Hens  Physics and Chemistry of
Nanostructures, Department of Chemistry, Ghent University, 9000 Ghent, Belgium
Pieter Geiregat  Physics and Chemistry of Nanostructures, Department of
Chemistry, Ghent University, 9000 Ghent, Belgium<EMAIL_ADDRESS>
###### Abstract

Atomically thin two-dimensional (2D) semiconductors are extensively
investigated for opto-electronic applications that require strong light-matter
interactions. In view of such applications, it is essential to understand how
(photo)excitation alters the non-linear optical response of these materials
under high carrier density conditions. Broadband transient absorption (TA)
spectroscopy is by now a widely used tool to study semiconductor physics in
such highly excited systems. However, the complex interplay between different
many-body interactions in 2D materials produces highly congested spectral
information and an ensuing non-trivial nonlinear photo-response, thereby
masking the desired intrinsic photophysics. Herein, we outline a concise
roadmap for analyzing such congested datasets based on examples of TA analysis
of various 2D materials. In particular, we emphasize the synergy between an
initial qualitative understanding of the transient photo-response based on
line shapes and their derivatives, and a consequent quantitative spectral
deconvolution backed by such insights.

###### keywords:

two-dimensional materials, transition metal dichalcogenides, transient
absorption, exciton, pump-probe

Center for Nano and Biophotonics, Ghent University, 9000 Ghent, Belgium Center
for Nano and Biophotonics, Ghent University, 9000 Ghent, Belgium Center for
Nano and Biophotonics, Ghent University, 9000 Ghent, Belgium

##  1 Introduction

Advancements in synthetic chemistry and material characterization over the
last decades have prompted the discovery of strongly confined two-dimensional
(2D) semiconductors with atomically thin dimensions. The most notable examples
include the large family of transition metal di-chalcogenides (TMDs), cadmium
chalcogenide nanoplatelets, and the recently emerging library of layered 2D
perovskites (Figure 1a).1, 2, 3, 4 In all these examples, a strong out-of-
plane carrier confinement alongside a reduced dielectric screening of the
Coulomb interaction enhances many-body interactions between photo-excited
charge carriers. As a consequence, otherwise unbound electrons and holes
condense into tightly bound excitons with unparalleled exciton binding
energies at room temperature.5, 6, 7, 8 At higher densities, these excitons
can further condense into bi-excitonic molecules, as observed in several 2D
materials, even at room temperature. Compared to continuum transitions
involving unbound carriers, light-matter interaction at the exciton transition
energy is greatly enhanced, causing a pronounced absorption of light with
large optical absorption coefficients in the order of
$\mathrm{10^{6}\>cm^{-1}}$.1, 9 In summary, the envisaged prospect of strong
absorption and, in parallel, efficient emission through excitons at room
temperature, alongside the tantalizing possibility to explore multi-exciton
physics, is what stimulates research into these atomically thin
semiconductors.10, 11, 12, 13, 14

Figure 1: (a) (top) Exemplary optical absorbance spectra of two-dimensional
semiconductors including (from top to bottom) $\mathrm{MoS}_{2}$ and
$\mathrm{WS}_{2}$,15 2D $\mathrm{CdSe}$ nanoplatelets,16 and 2D hybrid organic
inorganic perovskites.17 The optical absorbance is characterized by a sharp
exciton line superimposed on a 2D DOS-like continuum. (bottom) Schematic
depicting the absorbance spectrum of a two-dimensional excitonic semiconductor
with an optical band gap $E_{\mathrm{g,o}}$ and an exciton transition at
binding energy $E_{\mathrm{b,X}}$ superimposed on a step-like continuum at
energy $E_{g}$. (b) Schematic depicting the principles behind ultrafast pump-
probe spectroscopy. A specimen is pumped with an ultrafast laser pulse,
inducing changes in the absorbance. A delayed probe pulse monitors changes in
absorbance. Generally, experiments are mapped by tracking changes in the
differential absorbance ($\Delta A$) as a function of (ordinate) pump probe
delay and (abscissa) wavelength/energy, of which an example color map is shown
for colloidal MoS2. Adapted with permission from Reference 18. Copyright 2019
American Chemical Society. (c) Pictorial representation of the most common
transient effects after photo-excitation, including (from left to right)
changes in the integrated absorbance brought by state filling and Coulomb
screening, shifts in the transition energy due to band gap renormalization and
binding energy reduction, and linewidth changes induced by collisional
broadening through carrier-carrier and carrier-phonon scattering. (d) Spectral
slices of $\Delta A$ evaluated at the $A$ (660 nm) and $B$ (610 nm) transition
energies, extracted from the color map in panel c, highlighting the difficulty
to extract reliable information based on slices at a fixed probe wavelength.

The study of the underlying ultrafast phenomena dictating exciton formation
and their dissociation or association in strongly excited systems heavily
relies on various forms of pump-probe spectroscopy, e.g., the increasingly
prevalent transient absorption (TA) spectroscopy where visible/near-infrared
light pulses are used as pump and/or probe. Through these methods exciton
properties can be disentangled at otherwise inaccessible time and energy
scales (Figure 1c).19 Focusing on 2D nanomaterials, initial pump-probe studies
were often limited to single color/wavelength probing due to the technological
limitations of measuring single flakes with broadband probes using
microscopy.20, 21, 22 Even though such an approach captures the charge carrier
dynamics, looking to the spectrally strongly varying linear and transient
absorption in Figures 1a and 1b, one can understand that this often impedes a
consistent interpretation. Lately, however, the development of scalable
production methods – liquid-phase exfoliation (LPE), colloidal synthesis, and
large area chemical vapor deposition (CVD) – together with increasing
availability of TA microscopy equipment has resulted in a shift toward
broadband TA spectroscopy.23, 24, 18, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34,
35, 36, 37, 38

Such broadband TA studies highlighted that the photophysics of 2D materials is
far from trivial to interpret based on spectral information. Figure 1b depicts
an example of a 2D time-wavelength map obtained after photo-excitation of
MoS2, a typical 2D material, featuring a complex interplay between spectral
shifts, linewidth changes and changes in the integrated absorbance of the
exciton resonances. These effects, summarized in Figure 1c, on the line shape
and probability of transitions are generally rationalized based on the many-
body theory of quasi-2D quantum wells (QWs) which predicts drastic changes in
the optical absorption of a QW when additional (free) carriers are present in
the system.39, 40 Most notably, photo-excited carriers are predicted to modify
the optical band gap through band gap renormalization (BGR) and exciton
binding energy reduction (BER), broaden the exciton linewidth due to
collisional broadening, and reduce its oscillator strength through state
filling (Pauli blocking) and Coulomb screening of the e-h attraction (Figure
1c). As these effects typically occur simultaneously, the spectral information
is often highly congested and non-trivial to interpret, and kinetic traces can
widely differ depending on which probe wavelength they are evaluated at
(Figure 1d). In particular, these concerns raise questions on how exciton
parameters can be reliably extracted from such spectrally complex $\Delta A$
maps.

Herein, we outline a concise roadmap for analyzing spectrally complex
transient absorbance ($\Delta A$) maps based on recent literature examples.
The eventual goal is to obtain a rational methodology to go from congested 2D
time-wavelength maps to clear kinetic information on exciton parameters as the
exciton energy change $\Delta E$, the linewidth change $\Delta w$ and the
change of oscillator strength or the integrated change in absorbance $\Delta
S$. To this end, we first develop a series of qualitative insights for each
observable into the shape of transient spectra based on a derivative analysis
of the well-understood linear absorption spectrum. In particular, we develop
this approach for spectral shifts ($\Delta E$) and linewidth changes ($\Delta
w$). Based on these insights, we can more rationally optimize global fitting
procedures to attain the net $\Delta S$, often directly interpreted as the
’carrier dynamics’.

###  1.1 Spectral Shifts

As depicted in Figure 1a, the optical band gap $E_{\mathrm{g,o}}$ of an
excitonic semiconductor depends on the free particle band gap $E_{g}$ and the
exciton binding energy $E_{\mathrm{b,X}}$, such that $E_{\mathrm{g,o}}$ equals
$E_{\mathrm{g}}-E_{\mathrm{b,X}}$.5 Compared to an unexcited semiconductor,
carrier creation and/or subsequent heat generation through photo-excitation
alters the free particle band gap and the exciton binding energy.41 A photo-
induced shrinkage of $E_{g}$, referred to as band gap renormalization, results
from screening of repulsive Coulomb and exchange-correlation interactions.42
Similarly, photo-induced changes in lattice temperature and non-equilibrium
phonon populations after cooling of hot carriers affect the free particle band
gap, where typically a higher lattice temperature decreases $E_{g}$.30 As
shown in Figure 2a, both effects manifest themselves as a shrinkage of the
optical band gap, which concomitantly red-shifts transition energies.
Additionally, in excitonic systems, photo-excited e-h pairs dynamically screen
(or reduce) the electron-hole attraction that binds the exciton, thereby
reducing the exciton binding energy and reducing its oscillator strength. As
illustrated in Figure 2b, binding energy reduction (BER) shifts the exciton
absorption toward higher energy (blue-shift). Both mechanisms work in
spectrally opposite directions, yet typically do not fully compensate in
quasi-two dimensional systems. 43 As such, the direction of spectral shifts
observed experimentally depends on an intricate balance between BER and BGR.44

Figure 2: Schematics depicting (a) a spectral red-shift due to band gap
renormalization and (b) a spectral blue-shift due to binding energy
renormalization. $E_{\mathrm{g,o}}$, $E_{\mathrm{B,X}}$ and $\Delta$E denote
the optical band gap, the exciton binding energy and the spectral shift,
respectively. (c,d) Schematics depicting the effect of a spectral (c) red-
shift and (d) blue-shift on the transient absorbance $\Delta A$. A spectral
shift results in a $\Delta A$ spectrum resembling a first-derivative-like
spectrum $dA_{0}/dE$. (e) The transient response of colloidal molybdenum
disulfide nanosheets as an exemplary color map of a red-shift.18 Particularly
at longer time delays, $\Delta A$ shows multiple bands of bleach and PIA
mirrored around the $C$-band, $B$ and $A$ exciton. Adapted with permission
from Reference 18. Copyright 2019 American Chemical Society. (f) Traces of
$\Delta A$ probed at different excitation energies, extracted from color maps
similar to those shown in panel e, together with the first derivative of the
linear absorbance to the energy (grey) $dA_{0}/dE$. The spectral shape of
$\Delta A$ is independent on the excitation energy and resembles $dA_{0}/dE$,
as is indicative of a spectral red-shift of all transition energies due to BGR
dominating over BER. (g) The transient response of bilayer lead iodide
perovskites as a didactic color map and spectral slice showing a clear blue-
shift of the exciton transition.45 $\Delta A$ shows an alternate sequence of
PIA and bleach mirrored around the $1S$ exciton. Reprinted with permission
from Reference 45. Copyright 2015 American Chemical Society.

#### Effect on $\Delta A$ and Derivative Analysis

Figures 2c and 2d illustrate the effect of a change in transition energy on
the transient absorbance $\Delta A$ through the shift of an isolated Gaussian
absorption band to lower energy (red-shift) and higher energy (blue-shift),
respectively. Clearly, such a spectral shift reduces the absorbance at one
side of the absorption band, while increasing it at the other side, without
changing the integrated absorbance. A spectral shift therefore yields a
$\Delta A$ spectrum that resembles the first derivative of the absorbance
spectrum, consisting of an alternated sequence between bleach (negative
$\Delta A$, blue) and photo-induced absorption (PIA, positive $\Delta A$, red)
for the example shown. Note that the specific bleach/PIA sequence depends on
the direction of the transition’s energy shift.

A particularly useful approach to further analyse spectral shifts involves
combining such a conceptual interpretation with an analysis based on the
derivatives of the absorbance spectrum given that, for small enough shifts,
the ratio of $\Delta A$ and $dA_{0}/dE$ is directly proportional to the
spectral shift $\Delta E$ according to:18, 46

$\Delta E\approx\frac{\Delta A}{(dA_{0}/dE)}$ (1)

We illustrate this point in Figure 2f by plotting traces of $\Delta A$,
extracted from similar color maps to the one shown for $\mathrm{MoS}_{2}$ in
Figure 2e, together with $dA_{0}/dE$. Clearly, features in the TA spectrum
coincide with local minima and maxima in $dA_{0}/dE$. From Equation 1, we
could further estimate the magnitude of the spectral red-shift, which in this
case amounted to $\approx\mathrm{20\>meV}$.18 Importantly, this number closely
matched values extracted from a spectral deconvolution through a global fit
(see Figure 4), making this approach valuable to validate numerical fitting
approaches.

#### Literature Examples

Throughout literature, spectral red-shifts have been reported for various
transition metal di-chalcogenides produced by exfoliation, colloidal synthesis
and chemical vapor deposition.23, 38, 28, 24, 18, 25, 26, 27, 37 Often, those
spectral red-shifts were accompanied by an exciton linewidth broadening,24,
23, 25, 18, 27, 30 see further. The color map recorded on few-layer
colloidally dispersed molybdenum disulfide (MoS_2) in Figure 2e is a typical
example of such a red-shifted $\Delta A$ spectrum.18 Particularly at time
delays beyond $10\>\mathrm{ps}$, the transient absorbance features a
consecutive sequence of bleach and PIA bands centered around the excitonic
resonances, which reflects band gap renormalization.26 Similarly, blue-shifted
exciton transitions have been observed in quantum wells43, 47, 48, 49,
atomically-thin TMDs,28, 29, 30, 31, 32 and 2D perovskites.45 In
$\mathrm{MoS}_{2}$29 and $\mathrm{WS_{2}}$28, 30, 31, 32, blue-shifts at early
delays are typically accompanied by a broadening of the exciton band due to
carrier scattering, see further section on broadening. Cross-overs from blue-
to-red-shifts have been observed in $\mathrm{WS_{2}}$,$\mathrm{WSe_{2}}$, and
$\mathrm{ReS_{2}}$, evidencing the intricate balance between BGR and BER on
the excitation density.33, 28, 31, 34, 35 A revealing case of a relatively
uncongested spectral blue-shift is provided in Figure 2g.45 Here, the $\Delta
A$ map of a bilayer lead iodide perovskite shows a first-derivative-like
signature, mirrored around the fundamental exciton at 570 nm, with the PIA
band being at the high energy side and bleach band at the lower energy side of
the exciton line. Take note that the sequence of PIA/bleach is opposite from
the one shown in Figures 2e-f and thus characterizes a blue-shift.

###  1.2 Linewidth Changes

The homogeneous linewidth of an exciton transition is determined by the de-
phasing rate or coherence time of the exciton state(s).50 A loss of exciton
coherence can occur due to (in)elastic scattering of excitons with carriers or
lattice vibrations, generally referred to as carrier-carrier and carrier-
phonon scattering. Collisional broadening due to an increase in carrier
population or lattice temperature is therefore the dominant broadening
mechanism of absorption (and emission) line(s) in semiconductors, see Figure
3a.43, 29 In reality, transitions can also inhomogeneously broaden due to
local potentials arising from defects.51, 52 While being less prominent in the
literature, photo-induced line narrowing has been reported as a result of
carrier-induced screening of parasitic local defect-induced potentials or
externally applied and built-in electric fields, see Figure 3b.53, 39, 34, 54,
55, 56, 57, 58

Figure 3: Schematics depicting (a) linewidth broadening $\Delta w>0$ due to
collisional broadening and (b) linewidth narrowing $\Delta w>0$ due to
screening of an internal electric field. $w$ and $\Delta\tau$ respectively
denote the width of the exciton band, and the exciton de-phasing time. (c,d)
Schematics depicting the effect of a linewidth (c) broadening and (d)
narrowing on the transient absorbance $\Delta A$. Linewidth changes lead to
$\Delta A$ spectra that resemble the second derivative of the linear
absorbance $A_{0}$ to the energy $E$ ($d^{2}A_{0}/dE^{2}$). (e) Overview of
the transient absorbance $\Delta A$ recorded on a film of few-layer
$\mathrm{ReS}_{2}$, showing a complex interplay between several spectral
effects. The dashed line represents the linear absorbance spectrum $A_{0}$.
(f) Traces of the $\Delta A$ map in panel e sliced at specific pump-probe time
delays, highlighting the evolution of (i) a line broadening compounded by a
blue-shift into (ii-iii) a line narrowing (iv-v) congested by a red-shift at
an even longer delay. (g,h) Overlay plots of $\Delta A$ with
$d^{2}A_{0}/dE^{2}$, highlighting a one-to-one correspondence of a
broadened/narrowed $\Delta A$ spectrum with a second derivative spectrum. The
grey and black lines represent the derivatives of $A_{0}$ and the isolated
exciton profile, respectively. Reprinted and adapted with permission from
Reference 34. Copyright 2021 American Chemical Society.

#### Effect on $\Delta A$ and Derivative Analysis

Figures 3c and 3d illustrate the effect of a linewidth change on the transient
absorbance $\Delta A$ by, respectively, broadening and narrowing an isolated
Gaussian absorption band. Clearly, broadening (narrowing) a resonance
decreases (increases) the absorbance maximum, while increasing (decreasing)
the intensity at the tails of the absorbance band. A central bleach sided by
two positive side bands of photoinduced absorption therefore reflects line
broadening, whereas the pattern is reversed for narrowed lines. Generally, a
linewidth change can thus be best described as a second-derivative-like
$\Delta A$ spectrum. For Gaussian line shapes and small enough width changes,
the ratio of $\Delta A$ and $d^{2}A_{0}/dE^{2}$ is directly proportional to
the change in width $\Delta\sigma=\sigma-\sigma_{0}$ according to:34

$\Delta\sigma\approx\frac{\Delta A}{(d^{2}A_{0}/dE^{2})\sigma_{0}}$ (2)

#### Literature Examples

The $\Delta A$ map of liquid-phase exfoliated (LPE) $\mathrm{ReS}_{2}$ in
Figure 3e serves to illustrate broadened and narrowed resonances, and their
relation to the second derivative of the linear absorbance to the energy.59,
34 A more detailed description of the peculiar photophysics observed in LPE
$\mathrm{ReS}_{2}$ can be found in Ref. 34. To analyse this map conceptually,
we present traces of $\Delta A$ at specific pump-probe time delays in Figure
3f and overlay plots of $\Delta A$ with $d^{2}A_{0}/dE^{2}$ in Figures 3g-h.
Shortly after photoexcitation ($t$=200 fs, spectrum i), the transient
resembles a second derivative-like spectrum (see Figure 3g), and therefore
predominantly reflects a broadened absorption profile. Notwithstanding the
overlap of $\Delta A$ with $d^{2}A_{0}/dE^{2}$, the bleach center of $\Delta
A$ is slightly off-set toward lower energy (see vertical line), likely due to
a superimposed blue-shift. From 1 ps onwards (spectra ii and iii), the
transient has inverted and mimics $-d^{2}A_{0}/dE^{2}$ (see Figure 3h), and
thus reflects a narrowed resonance. At longer time delays (spectra iv-v), a
red-shift grows in, resulting in a transient with a first derivative-like
shape. This example points out how a conceptual interpretation of the color
map based on simple line shapes, in addition to comparing slices of $\Delta A$
to first and second derivatives, can greatly contribute to decrypt TA maps
qualitatively, prior to attempting more elaborate curve fits.

###  1.3 Extracting Exciton Parameters through Curve Fits

Figure 4: Overview of the general principle behind a spectral deconvolution of
a $\Delta A$ map. (a) The transient absorbance $\Delta A$ is summed together
with the linear absorbance $A_{0}$, yielding the nonlinear absorbance $A$. (b)
$A_{0}$ and $A$ are fit to (filled curve) a Gaussian function and (blue dashed
line) a polynomial background, yielding the central position $E_{X,(0)}$, the
area $S_{(0)}$ and the width $w_{1/2,(0)}$ of the exciton band. (c) Repeating
step b as a function of pump-probe delay maps the time-dependent evolution of
these fit parameters. This is illustrated by fitting the $\Delta A$ map of
colloidal $\mathrm{MoS}_{2}$ in Figure 2e to a sum of two Gaussians,
accounting for the (red markers) $A$ and (blue markers) $B$ exciton, and a
polynomial background.18 Adapted with permission from Reference 18. Copyright
2019 American Chemical Society.

In the examples listed above, we emphasized how $\Delta A$ spectra often
result from a complex interplay between spectral shifts ($\Delta E$) and
linewidth changes ($\Delta w$) of the exciton band(s). Additionally, these
effects are supplemented by changes of the integrated absorbance of the
exciton transition ($\Delta S$). Such changes are linked to a reduction in
oscillator strength, which occurs both due to state filling/Pauli-blocking or
due to screening of the Coulomb potential.44 Since $\Delta S$ mostly scales
linearly with carrier density in low density regimes, it’s dynamics are a
direct metric for the carrier lifetime. As such, $\Delta S$ is of most
interest, e.g. in charge transfer studies or to evaluate trapping dynamics.60,
61 The sheer complexity of the resulting $\Delta A$ spectra made researchers
devise data analysis techniques to deconvolute and quantitatively analyse the
time-dependent evolution of these metrics. Throughout the literature, spectral
deconvolution approaches have been used to unravel charge carrier decay in a
variety of TMDs, including $\mathrm{MoS}_{2}$, $\mathrm{ReS}_{2}$, and
$\mathrm{WS_{2}}$.25, 34, 62, 63, 28, 23, 30, 27, 32 Generally, these methods
rely on a curve fit to the excited-state absorbance $A$ (or reflectance $R$).
We illustrate the basic principle behind this approach in Figure 4a starting
from a $\Delta A$ spectrum that results from a broadened, red-shifted and
bleached Gaussian absorption profile superimposed on a background – a shape
commonly encountered in the literature.37, 18, 23, 25, 30 First, the excited-
state or nonlinear absorbance ($A$) is computed by summing the transient
absorbance ($\Delta A$) to the ground-state or linear absorbance ($A_{0}$),
i.e., $A(E,t)=A_{0}(E)+\Delta A(E,t)$. Both ground-state and excited-state
absorbance are consequently fit to Gaussian (or where adequate Lorentzian)
line-shapes $G(E)$, describing the exciton band(s), and an eventual background
absorbance $C(E)$, accounting for the continuum energy transitions and/or
scattering:

$\displaystyle A(E)$ $\displaystyle=$ $\displaystyle G(E)+C(E)$ (3)
$\displaystyle=$ $\displaystyle
h\exp((-E-E_{X})^{2}/w^{2})+\sum_{j}^{m}K_{j}E^{j}$ (4)

Here, $h$, $E_{X}$ and $w$ are, respectively, the amplitude, central position
and linewidth of the Gaussian (see Figure 4b for a graphical representation),
whereas the area $S$ scales with the product of the $w$ and $h$. The
background absorbance $C(E)$, on the other hand, is typically implemented as a
polynomial with order $m$ and pre-factors $K_{j}$.

Figure 5: Summary of the roadmap to interpret congested $\Delta A$ spectra and
obtain kinetic information on intrinsic variables as exciton energy shifts
$\Delta E$, linewidth changes $\Delta w$ and changes in the integrated
absorbance $\Delta S$ of the exciton resonance. This toolbox includes a
continuous feedback loop between (left) a qualitative interpretation based on
line shapes and derivatives, and (right) a quantitative analysis based on a
spectral deconvolution through rational curve fits of the full absorbance
$A_{0}$+$\Delta A$. The transient absorbance $\Delta A$ is represented by the
black curve filled in red (positive $\Delta A$) and blue (negative $\Delta
A$). The filled grey curve is the linear absorbance $A_{0}$ and the red curve
is the nonlinear absorbance $A$ computed as $A_{0}$+$\Delta A$.

In this way, by comparing these parameters before and after photo-excitation,
the transient absorbance is deconvolved into a clear set of parameters: photo-
induced areal changes, spectral shifts and linewidth changes. Repeating this
procedure for each pump-probe delay maps the time-dependent evolution of these
fit parameters. We illustrate this method in Figure 4c by fitting the $\Delta
A$ map of $\mathrm{MoS}_{2}$ found in Figure 2e. From such a deconvolution, it
follows that the transient absorbance resulted from a reduced oscillator
strength and broadening of the $A$ and $B$ excitons and a red-shift of the
entire absorbance spectrum; a result conform with a qualitative assessment
based on line shapes and derivatives, see Figure 4. As in particular the $A$
exciton bleach exhibited a rapid first-order decay, one could propose that
non-radiative recombination proceeded through trapping of photoexcited holes
in mid-gap states.18 This example shows how exciton parameters, and their
dynamics, can be extracted through curve fits, backed up by insights obtained
from qualitative interpretation based on line shapes and derivatives (see
overview of the roadmap in Figure 5). Combining such ’blind’ output of a
fitting algorithm with qualitative insights allows us to guide the former in
the right direction. For example, based on the derivative insights we expected
redshifts in the order of 20 meV. If the fitting routine would produce much
smaller shifts, or even blueshifts, we had to re-adjust our numerical approach
by, for example, restricting certain parameters.

##  2 Conclusion

In summary, we outlined a concise roadmap for analyzing spectrally complex
transient absorbance maps based on recent literature examples. The roadmap,
summarized in Figure 5, offers a rational methodology to go from congested 2D
time-wavelength maps to clear kinetic information on exciton parameters as the
exciton energy change $\Delta E$, the linewidth change $\Delta w$ and the
change of oscillator strength or the integrated change in absorbance $\Delta
S$. In particular, we emphasized the synergy between an initial qualitative
understanding of the transient photo-response based on line shapes and their
derivatives, and a consequent quantitative spectral deconvolution backed by
such insights.

##  3 Acknowledgements

P.S. acknowledges the FWO-Vlaanderen (1S40117N) and VLAIO
(AIO.ONO.2020.0007.01) for funding. P.G. acknowledges support from FWO-
Vlaanderen (12K8216N). Z.H. acknowledges the Research Foundation Flanders
(research projects 17006602 and G0F0920N) and Ghent University (GOA no.
01G01513) for funding.

## References

  * Mak et al. 2010 Mak, K. F.; Lee, C.; Hone, J.; Shan, J.; Heinz, T. F. Atomically Thin $MoS_{2}$: a New Direct-gap Semiconductor. _Phys. Rev. Lett._ 2010, _105_ , 136805
  * Splendiani et al. 2010 Splendiani, A.; Sun, L.; Zhang, Y.; Li, T.; Kim, J.; Chim, C.-Y.; Galli, G.; Wang, F. Emerging Photoluminescence in Monolayer $MoS_{2}$. _Nano Lett._ 2010, _10_ , 1271–1275
  * Ithurria et al. 2011 Ithurria, S.; Tessier, M.; Mahler, B.; Lobo, R.; Dubertret, B.; Efros, A. L. Colloidal Nanoplatelets with Two-dimensional Electronic Structure. _Nat. Mater._ 2011, _10_ , 936–941
  * Shi et al. 2018 Shi, E.; Gao, Y.; Finkenauer, B. P.; Coffey, A. H.; Dou, L., et al. Two-Dimensional Halide Perovskite Nanomaterials and Heterostructures. _Chem. Soc. Rev._ 2018, _47_ , 6046–6072
  * Wang et al. 2018 Wang, G.; Chernikov, A.; Glazov, M. M.; Heinz, T. F.; Marie, X.; Amand, T.; Urbaszek, B. Colloquium: Excitons in Atomically Thin Transition Metal Dichalcogenides. _Rev. Mod. Phys._ 2018, _90_ , 021001
  * Naeem et al. 2015 Naeem, A.; Masia, F.; Christodoulou, S.; Moreels, I.; Borri, P.; Langbein, W. Giant Exciton Oscillator Strength and Radiatively Limited Dephasing in Two-dimensional Platelets. _Phys. Rev. B_ 2015, _91_ , 121302
  * Ugeda et al. 2014 Ugeda, M. M.; Bradley, A. J.; Shi, S.-F.; Da Jornada, F. H.; Zhang, Y.; Qiu, D. Y.; Ruan, W.; Mo, S.-K.; Hussain, Z.; Shen, Z.-X. et al. Giant Bandgap Renormalization and Excitonic Effects in a Monolayer Transition Metal Dichalcogenide Semiconductor. _Nat. Mat._ 2014, _13_ , 1091–1095
  * Wang et al. 2018 Wang, Q.; Liu, X.-D.; Qiu, Y.-H.; Chen, K.; Zhou, L.; Wang, Q.-Q. Quantum Confinement Effect and Exciton Binding Energy of Layered Perovskite Nanoplatelets. _AIP Adv._ 2018, _8_ , 025108
  * Bernardi et al. 2013 Bernardi, M.; Palummo, M.; Grossman, J. C. Extraordinary Sunlight Absorption and One Nanometer Thick Photovoltaics Using Two-Dimensional Monolayer Materials. _Nano Lett._ 2013, _13_ , 3664–3670
  * Chen et al. 2014 Chen, Z.; Nadal, B.; Mahler, B.; Aubin, H.; Dubertret, B. Quasi-2D Colloidal Semiconductor Nanoplatelets for Narrow Electroluminescence. _Adv. Funct. Mater._ 2014, _24_ , 295–302
  * Grim et al. 2014 Grim, J. Q.; Christodoulou, S.; Di Stasio, F.; Krahne, R.; Cingolani, R.; Manna, L.; Moreels, I. Continuous-Wave Biexciton Lasing at Room Temperature Using Solution-Processed Quantum Wells. _Nat. Nanotechnol._ 2014, _9_ , 891–895
  * Ye et al. 2015 Ye, Y.; Wong, Z. J.; Lu, X.; Ni, X.; Zhu, H.; Chen, X.; Wang, Y.; Zhang, X. Monolayer Excitonic Laser. _Nat. Photonics_ 2015, _9_ , 733
  * Ross et al. 2014 Ross, J. S.; Klement, P.; Jones, A. M.; Ghimire, N. J.; Yan, J.; Mandrus, D.; Taniguchi, T.; Watanabe, K.; Kitamura, K.; Yao, W. et al. Electrically Tunable Excitonic Light-Emitting Diodes Based on Monolayer $WSe_{2}$ p–n Junctions. _Nat. Nanotechnol._ 2014, _9_ , 268–272
  * Cheng et al. 2014 Cheng, R.; Li, D.; Zhou, H.; Wang, C.; Yin, A.; Jiang, S.; Liu, Y.; Chen, Y.; Huang, Y.; Duan, X. Electroluminescence and Photocurrent Generation from Atomically Sharp $WSe_{2}$/$MoS_{2}$ Heterojunction p–n Diodes. _Nano Lett._ 2014, _14_ , 5590–5597
  * Li et al. 2014 Li, Y.; Chernikov, A.; Zhang, X.; Rigosi, A.; Hill, H. M.; Van Der Zande, A. M.; Chenet, D. A.; Shih, E.-M.; Hone, J.; Heinz, T. F. Measurement of The Optical Dielectric Function of Monolayer Transition-metal Dichalcogenides: $MoS_{2}$, $MoSe_{2}$, $WS_{2}$, and $WSe_{2}$. _Phys. Rev.B_ 2014, _90_ , 205422
  * Di Giacomo et al. 2020 Di Giacomo, A.; Roda, C.; Khan, A. H.; Moreels, I. Colloidal Synthesis of Laterally Confined Blue-Emitting 3.5 Monolayer CdSe Nanoplatelets. _Chem. Mater_ 2020, _32_ , 9260–9267
  * Song et al. 2020 Song, B.; Hou, J.; Wang, H.; Sidhik, S.; Miao, J.; Gu, H.; Zhang, H.; Liu, S.; Fakhraai, Z.; Even, J. et al. Determination of Dielectric Functions and Exciton Oscillator Strength of Two-Dimensional Hybrid Perovskites. _ACS Mater. Lett._ 2020, _3_ , 148–159
  * Schiettecatte et al. 2019 Schiettecatte, P.; Geiregat, P.; Hens, Z. Ultrafast Carrier Dynamics in Few-Layer Colloidal Molybdenum Disulfide Probed by Broadband Transient Absorption Spectroscopy. _J. Phys. Chem. C_ 2019, _123_ , 10571–10577
  * Maiuri et al. 2019 Maiuri, M.; Garavelli, M.; Cerullo, G. Ultrafast Spectroscopy: State of the Art and Open Challenges. _J. Am. Chem. Soc._ 2019, _142_ , 3–15
  * Wang et al. 2014 Wang, H.; Zhang, C.; Rana, F. Ultrafast Dynamics of Defect-Assisted Electron–Hole Recombination in Monolayer MoS2. _Nano Lett._ 2014, _15_ , 339–345
  * Wang et al. 2012 Wang, R.; Ruzicka, B. A.; Kumar, N.; Bellus, M. Z.; Chiu, H.-Y.; Zhao, H. Ultrafast and Spatially Resolved Studies of Charge Carriers in Atomically Thin Molybdenum Disulfide. _Phys. Rev. B_ 2012, _86_ , 045406
  * Sun et al. 2014 Sun, D.; Rao, Y.; Reider, G. A.; Chen, G.; You, Y.; Brezin, L.; Harutyunyan, A. R.; Heinz, T. F. Observation of Rapid Exciton–Exciton Annihilation in Monolayer Molybdenum Disulfide. _Nano Lett._ 2014, _14_ , 5625–5629
  * Zhou et al. 2019 Zhou, P.; Tanghe, I.; Schiettecatte, P.; Van Thourhout, D.; Hens, Z.; Geiregat, P. Ultrafast Carrier Dynamics in Colloidal $WS_{2}$ Nanosheets Obtained Through a Hot Injection Synthesis. _J. Chem. Phys._ 2019, _151_ , 164701
  * Chernikov et al. 2015 Chernikov, A.; Ruppert, C.; Hill, H. M.; Rigosi, A. F.; Heinz, T. F. Population inversion and giant bandgap renormalization in atomically thin WS 2 layers. _Nature Photonics_ 2015, _9_ , 466–470
  * Cunningham et al. 2016 Cunningham, P. D.; McCreary, K. M.; Hanbicki, A. T.; Currie, M.; Jonker, B. T.; Hayden, L. M. Charge Trapping and Exciton Dynamics in Large-Area CVD Grown MoS2. _J. Phys. Chem. C._ 2016, _120_ , 5819–5826
  * Pogna et al. 2016 Pogna, E. A. A.; Marsili, M.; De Fazio, D.; Dal Conte, S.; Manzoni, C.; Sangalli, D.; Yoon, D.; Lombardo, A.; Ferrari, A. C.; Marini, A. et al. Photo-Induced Bandgap Renormalization Governs the Ultrafast Response of Single-Layer MoS2. _ACS Nano_ 2016, _10_ , 1182–1188
  * Sie et al. 2017 Sie, E. J.; Steinhoff, A.; Gies, C.; Lui, C. H.; Ma, Q.; Rosner, M.; Schonhoff, G.; Jahnke, F.; Wehling, T. O.; Lee, Y.-H. et al. Observation of Exciton Redshift–Blueshift Crossover in Monolayer WS2. _Nano Lett._ 2017, _17_ , 4210–4216
  * Cunningham et al. 2017 Cunningham, P. D.; Hanbicki, A. T.; McCreary, K. M.; Jonker, B. T. Photoinduced Bandgap Renormalization and Exciton Binding Energy Reduction in $WS_{2}$. _ACS Nano_ 2017, _11_ , 12601–12608
  * Sim et al. 2013 Sim, S.; Park, J.; Song, J.-G.; In, C.; Lee, Y.-S.; Kim, H.; Choi, H. Exciton Dynamics in Atomically Thin MoS2: Interexcitonic Interaction and Broadening Kinetics. _Phys. Rev. B_ 2013, _88_ , 075434
  * Ruppert et al. 2017 Ruppert, C.; Chernikov, A.; Hill, H. M.; Rigosi, A. F.; Heinz, T. F. The Role of Electronic and Phononic Excitation in the Optical Response of Monolayer $WS_{2}$ after Ultrafast Excitation. _Nano Lett._ 2017, _17_ , 644–651
  * Calati et al. 2021 Calati, S.; Li, Q.; Zhu, X.; Stähler, J. Ultrafast Evolution of the Complex Dielectric Function of Monolayer $\mathrm{WS_{2}}$ after Photoexcitation. _Phys. Chem. Chem. Phys._ 2021, _23_ , 22640–22646
  * Trovatello et al. 2022 Trovatello, C.; Katsch, F.; Li, Q.; Zhu, X.; Knorr, A.; Cerullo, G.; Dal Conte, S. Disentangling Many-Body Effects in the Coherent Optical Response of 2D Semiconductors. _Nano Letters_ 2022, _22_ , 5322–5329
  * Aivazian et al. 2017 Aivazian, G.; Yu, H.; Wu, S.; Yan, J.; Mandrus, D. G.; Cobden, D.; Yao, W.; Xu, X. Many-body Effects in Nonlinear Optical Responses of 2D Layered Semiconductors. _2D Materials_ 2017, _4_ , 025024
  * Schiettecatte et al. 2021 Schiettecatte, P.; Poonia, D.; Tanghe, I.; Maiti, S.; Failla, M.; Kinge, S.; Hens, Z.; Siebbeles, L. D.; Geiregat, P. Unraveling the Photophysics of Liquid-Phase Exfoliated Two-Dimensional ReS2 Nanoflakes. _J. Phys. Chem. C_ 2021, _125_ , 20993–21002
  * Calati et al. 2022 Calati, S.; Li, Q.; Zhu, X.; Stähler, J. Dynamic Screening of Quasiparticles in $\mathrm{WS_{2}}$ Monolayers. _arXiv preprint arXiv:2204.02125_ 2022,
  * Taank et al. 2022 Taank, P.; Karmakar, R.; Sharma, R.; Yadav, R. K.; Shrivastava, M.; Maurya, N. C.; Maji, T. K.; Karmakar, D.; Adarsh, K. An Insightful Picture of Multi-Particle Recombination in Few-Layer MoS2 Nanosheets. _J. Phys. Chem. C_ 2022, _126_ , 416–422
  * Karmakar et al. 2021 Karmakar, M.; Mukherjee, S.; Ray, S. K.; Datta, P. K. Electron-Hole Plasma Formation Dynamics Observed Through Exciton-Plasma Interactions in Transition Metal Dichalcogenides. _Phys. Rev. B_ 2021, _104_ , 075446
  * Bera et al. 2021 Bera, S. K.; Shrivastava, M.; Bramhachari, K.; Zhang, H.; Poonia, A. K.; Mandal, D.; Miller, E. M.; Beard, M. C.; Agarwal, A.; Adarsh, K. Atomlike Interaction and Optically Tunable Giant Band-Gap Renormalization in Large-Area Atomically Thin $\mathrm{MoS_{2}}$. _Phys. Rev. B_ 2021, _104_ , L201404
  * Schmitt-Rink et al. 1989 Schmitt-Rink, S.; Chemla, D.; Miller, D. Linear and Nonlinear Optical Properties of Semiconductor Quantum Wells. _Advances in Physics_ 1989, _38_ , 89–188
  * Haug and Koch 2009 Haug, H.; Koch, S. W. _Quantum Theory of the Optical and Electronic Properties of Semiconductors_ ; World Scientific Publishing Company, 2009
  * Schmitt-Rink et al. 1986 Schmitt-Rink, S.; Ell, C.; Haug, H. Many-body Effects in the Absorption, Gain, and Luminescence Spectra of Semiconductor Quantum-well Structures. _Phys. Rev. B_ 1986, _33_ , 1183
  * Kalt and Rinker 1992 Kalt, H.; Rinker, M. Band-Gap Renormalization in Semiconductors With Multiple Inequivalent Valleys. _Phys. Rev. B_ 1992, _45_ , 1139
  * Wake et al. 1992 Wake, D.; Yoon, H.; Wolfe, J.; Morkoc, H. Response of Excitonic Absorption Spectra to Photoexcited Carriers in GaAs Quantum Wells. _Phys. Rev. B_ 1992, _46_ , 13452
  * Steinhoff et al. 2014 Steinhoff, A.; Rosner, M.; Jahnke, F.; Wehling, T. O.; Gies, C. Influence of Excited Carriers on the Optical and Electronic Properties of $MoS_{2}$. _Nano Lett._ 2014, _14_ , 3743–3748
  * Wu et al. 2015 Wu, X.; Trinh, M. T.; Zhu, X.-Y. Excitonic Many-body Interactions in Two-dimensional Lead Iodide Perovskite Quantum Wells. _J. Phys. Chem. C_ 2015, _119_ , 14714–14721
  * Geiregat et al. 2014 Geiregat, P.; Houtepen, A.; Justo, Y.; Grozema, F. C.; Van Thourhout, D.; Hens, Z. Coulomb Shifts upon Exciton Addition to Photoexcited PbS Colloidal Quantum Dots. _J. Phys. Chem. C_ 2014, _118_ , 22284–22290
  * Schultheis et al. 1986 Schultheis, L.; Kuhl, J.; Honold, A.; Tu, C. Ultrafast Phase Relaxation of Excitons via Exciton-Exciton and Exciton-Electron Collisions. _Phys. Rev. Lett._ 1986, _57_ , 1635
  * Leite et al. 1969 Leite, R.; Shah, J.; Gordon, J. Effect of Electron-Exciton Collisions on the Free-Exciton Linewidth in Epitaxial GaAs. _Phys. Rev. Lett._ 1969, _23_ , 1332
  * Hulin et al. 1986 Hulin, D.; Mysyrowicz, A.; Antonetti, A.; Migus, A.; Masselink, W.; Morkoc, H.; Gibbs, H.; Peyghambarian, N. Well-size Dependence of Exciton Blue Shift in GaAs Multiple-quantum-well Structures. _Phys. Rev. B_ 1986, _33_ , 4389
  * Moody et al. 2015 Moody, G.; Dass, C. K.; Hao, K.; Chen, C.-H.; Li, L.-J.; Singh, A.; Tran, K.; Clark, G.; Xu, X.; Berghäuser, G. et al. Intrinsic Homogeneous Linewidth and Broadening Mechanisms of Excitons in Monolayer Transition Metal Dichalcogenides. _Nat. Commun._ 2015, _6_ , 1–6
  * Stoneham 1969 Stoneham, A. Shapes of Inhomogeneously Broadened Resonance Lines in Solids. _Rev. Mod. Phys._ 1969, _41_ , 82
  * Raja et al. 2019 Raja, A.; Waldecker, L.; Zipfel, J.; Cho, Y.; Brem, S.; Ziegler, J. D.; Kulig, M.; Taniguchi, T.; Watanabe, K.; Malic, E. et al. Dielectric Disorder in Two-Dimensional Materials. _Nat. Nanotechnol._ 2019, _14_ , 832–837
  * Miller et al. 1985 Miller, D. A.; Chemla, D.; Damen, T.; Gossard, A.; Wiegmann, W.; Wood, T.; Burrus, C. Electric Field Dependence of Optical Absorption Near the Band Gap of Quantum-well Structures. _Phys. Rev. B_ 1985, _32_ , 1043
  * Shikanai et al. 2000 Shikanai, A.; Deguchi, T.; Sota, T.; Kuroda, T.; Tackeuchi, A.; Chichibu, S.; Nakamura, S. A Pump and Probe Study of Photoinduced Internal Field Screening Dynamics in an AlGaN/GaN Single-Quantum-Well Structure. _Appl. Phys. Lett._ 2000, _76_ , 454–456
  * Chen et al. 2003 Chen, F.; Cheung, M.; Sweeney, P. M.; Kirkey, W.; Furis, M.; Cartwright, A. Ultrafast differential transmission spectroscopy of excitonic transitions in InGaN/GaN multiple quantum wells. _Journal of applied physics_ 2003, _93_ , 4933–4935
  * Chen et al. 2003 Chen, F.; Kirkey, W.; Furis, M.; Cheung, M.; Cartwright, A. Excitonic Field Screening and Bleaching in InGaN/GaN Multiple Quantum Wells. _Solid State Commun._ 2003, _125_ , 617–622
  * Omae et al. 2003 Omae, K.; Kawakami, Y.; Fujita, S.; Narukawa, Y.; Mukai, T. Effects of Internal Electrical Field on Transient Absorption in $In_{x}Ga_{1-x}N$ Thin Layers and Quantum Wells with Different Thickness by Pump and Probe Spectroscopy. _Phys. Rev. B._ 2003, _68_ , 085303
  * Omae et al. 2002 Omae, K.; Kawakami, Y.; Fujita, S.; Yamada, M.; Narukawa, Y.; Mukai, T. Effects of Internal Electric Field and Carrier Density on Transient Absorption Spectra in a Thin GaN Epilayer. _Phys. Rev. B._ 2002, _65_ , 073308
  * Schiettecatte et al. 2020 Schiettecatte, P.; Rousaki, A.; Vandenabeele, P.; Geiregat, P.; Hens, Z. Liquid-Phase Exfoliation of Rhenium Disulfide by Solubility Parameter Matching. _Langmuir_ 2020, _36_ , 15493–15500
  * Maiti et al. 2021 Maiti, S.; Poonia, D.; Schiettecatte, P.; Hens, Z.; Geiregat, P.; Kinge, S.; Siebbeles, L. D. Generating Triplets in Organic Semiconductor Tetracene Upon Photoexcitation of Transition Metal Dichalcogenide $\mathrm{ReS_{2}}$. _J. Phys. Chem. Lett._ 2021, _12_ , 5256–5260
  * Markeev et al. 2022 Markeev, P. A.; Najafidehaghani, E.; Samu, G. F.; Sarosi, K.; Kalkan, S. B.; Gan, Z.; George, A.; Reisner, V.; Mogyorosi, K.; Chikan, V. et al. Exciton Dynamics in $\mathrm{MoS_{2}}$-Pentacene and $\mathrm{WSe_{2}}$-Pentacene Heterojunctions. _ACS Nano_ 2022,
  * Cunningham et al. 2019 Cunningham, P. D.; Hanbicki, A. T.; Reinecke, T. L.; McCreary, K. M.; Jonker, B. T. Resonant Optical Stark Effect in Monolayer $\mathrm{WS_{2}}$. _Nat. Commun._ 2019, _10_ , 1–8
  * Zhao et al. 2020 Zhao, J.; Zhao, W.; Du, W.; Su, R.; Xiong, Q. Dynamics of Exciton Energy Renormalization in Monolayer Transition Metal Disulfides. _Nano Res._ 2020, _13_ , 1399–1405


# Real-time Anomaly Detection for Multivariate Data Streams

Kenneth Odoh  Microsoft Corporation<EMAIL_ADDRESS>
(2018)

###### Abstract.

We present a real-time multivariate anomaly detection algorithm for data
streams based on the Probabilistic Exponentially Weighted Moving Average
(PEWMA). Our formulation is resilient to (abrupt transient, abrupt
distributional, and gradual distributional) shifts in the data. The novel
anomaly detection routines utilize an incremental online algorithm to handle
streams. Furthermore, our proposed anomaly detection algorithm works in an
unsupervised manner eliminating the need for labeled examples. Our algorithm
performs well and is resilient in the face of concept drifts.

signal processing, online learning

††copyright: acmcopyright††journalyear: 2018††doi:
10.1145/1122445.1122456††conference: Atlanta ’22: 31st ACM International
Conference on Information and Knowledge Management; October 17–22, 2022;
Atlanta, USA††price: 15.00††isbn: 978-1-4503-XXXX-X/18/06††ccs: Computing
methodologies Machine learning algorithms

##  1\. Introduction

Anomaly detection is the task of classifying patterns that depict abnormal
behavior. Outliers can arise due to (human/equipment) errors, faulty systems,
and others. Anomaly detection is well-suited for unbalanced data, where the
ideal scenario is to predict the behavior of the minority class. There are
numerous applications in detecting default on loans, fraud detection, and
network intrusion detections. An anomaly detection algorithm can work in
Unsupervised, Supervised, or hybrid mode.

There are different types of anomalies described as follows.

  * •

Point Anomaly: the algorithm decide a single instance as an anomaly concerning
the entire data set.

  * •

Contextual Anomaly: a data instance can be anomalous based on the context
(attributes and position in the stream) and proximity of the chosen anomaly.
This anomaly type is ideal for multivariate data, e.g. in the snapshot reading
of a machine, an attribute of a single data point may seem abnormal but can be
normal behavior based on consideration of the entire data.

  * •

Collective Anomaly: the algorithm decide a set of data points that are
anomalies as a group, but individually these data points exhibit normal
behaviors.

Anomaly detection algorithms can operate in the following setting:

  * •

Static: These algorithms work in static datasets. Every item is loaded into
memory at one time to perform computation.

  * •

Online: These algorithms work in real-time data streams. Items are
incrementally loaded into memory and processed in chunks.

  * •

Static + Online: The model may operate in two stages as initial parameters get
estimated in the static setting. The parameters are incrementally updated as
more data arrives. Our work is of this type.

PEWMA was introduced in the work (Carter and Streilein, 2012) for online
anomaly detection on univariate time series. Drawing inspiration from their
work, we have provided extensions to support real-time anomaly detection of a
multivariate data stream.

111Source code: https://github.com/kenluck2001/anomalyMulti, Blog:
https://kenluck2001.github.io/blog_post/real-
time_anomaly_detection_for_multivariate_data_stream.html

##  2\. Background

An anomaly detection algorithm can identify outliers in several signal changes
in time-dependent data. Signal changes are common in time-dependent data. They
include abrupt transient shift, abrupt distributional shift, and gradual
distributional shift (Carter and Streilein, 2012) labeled as ”A”, ”B”, and ”C”
in Figure 1 respectively.

Figure 1. Different types of signal changes: abrupt transient, abrupt
distributional shift, and gradual distributional shift (Carter and Streilein,
2012) (from left to right).

Online algorithms are useful for real-time applications, as they operate
incrementally on data streams. These algorithms incrementally receive input
and decide based on an updated parameter that conveys the current state of the
data stream. This philosophy contrasts with offline algorithms that assume the
entire data is available in memory. The issue with an offline algorithm is
that the data must fit in memory. The online algorithm should be both time and
space-efficient.

Anomaly detection algorithm can work in modes such as diagnosis and
accommodation (Hodge and Austin, 2004). Firstly, the diagnosis method finds
the outlier in the data and removes it from the data sample to avoid skewing
the distribution. This method is applicable when the distribution of expected
behaviors is known. The outliers get excluded when the estimation of the
parameters of the distribution (Hodge and Austin, 2004). Secondly, the
accommodation method finds the outliers in the data and includes estimating
the parameters of the statistical model. The accommodation method is suitable
for data streams that account for the effect of concept drift (Ditzler and
Polikar, 2013).

Muth (Muth, 1960) laid the foundation for exponential smoothing by showing
that it provides optimal estimates using random walk with some noise. Further
works were done to provide a statistical framework for exponential smoothing
leading to the creation of linear models (Box and Jenkins, 1976; Roberts,
1982). Unfortunately, this does apply to nonlinear models. Yule postulated a
formulation of time series as a realization of a stochastic process (Yule,
1927). Exponential Weighted Moving Average (EWMA) is ideal for keeping a set
of running moments in the data stream. EWMA is a smoothing technique that adds
a forgetting parameter to modify the influence of a recent item in the data
stream as shown in Equation 1. This smoothing causes volatility in abrupt
transient changes and is unfit for distribution shifts.

(1) $\mu_{t+1}=\alpha\mu_{t-1}+(1-\alpha)X_{t}$

EWMA limitations motivated the discovery of Probabilistic Exponentially
Weighted Moving Average (PEWMA). PEWMA (Carter and Streilein, 2012) improves
on EMWA by adding parameter which includes the probability of the data in the
model as shown in Equation 2. PEWMA works for every shift including abrupt
transient shift, abrupt distributional shift, and gradual distributional shift
respectively. PEWMA and EWMA make use of a damped window (Salehi and Rashidi,
2018).

(2) $\mu_{t+1}=\alpha(1-\beta P_{t})\mu_{t-1}+(1-\alpha(1-\beta P_{t}))X_{t}$

##  3\. Method

Our formulation provided an implementation of the online covariance matrix in
Subsection 3.2, alongside an online inverse covariance matrix based on
Sherman$-$Morrison formula in Subsection 3.3, and PEWMA in Subsection 3.1.

Our implementation of the online covariance matrix builds on work (Igel et
al., 2006). We simplify the algorithm by ignoring the details of evolutionary
computation in the paper. Our work adapted evolution as described in the paper
(Igel et al., 2006) as a transition from one generation to the next; is
equivalent to moving from one state to another state. This is analogous to how
online algorithms work with dynamic changes as new data enters the stream.
Cholesky decomposition is used extensively in the algorithms.

###  3.1. Probabilistic Exponentially Weighted Moving Average (PEWMA)

PEWMA (Carter and Streilein, 2012) algorithm works in the accommodation mode.
The routine shown in Algorithm [1] allows for concept drift (Ditzler and
Polikar, 2013), which occurs in data streams by updating the set of parameters
that convey the state of the stream.

Algorithm 1 Probabilistic Exponential Weighted Moving Average (Carter and
Streilein, 2012)

$X_{t},\hat{X_{t}},\hat{\alpha_{t}},T,t$

$P_{t},\hat{X_{t+1}},\hat{\alpha_{t+1}}$

    //incremental Z score

$Z_{t}\leftarrow\frac{X_{t}-\hat{X_{t}}}{\hat{\alpha_{t}}}$

    //probability density function

$P_{t}\leftarrow\frac{Z_{t}}{\sqrt{2\pi}}e^{\frac{Z_{t}}{2}}$

if $t<T$ then

    //increment standard deviation (training phase)

$\alpha_{t}\leftarrow 1-1/t$

else

    //increment standard deviation

$\alpha_{t}\leftarrow(1-\beta P_{t})\alpha$

end if

    //moving average

$s_{1}\leftarrow\alpha_{t}s_{1}+(1-\alpha_{t})X_{t}$

$s_{2}\leftarrow\alpha_{t}s_{2}+(1-\alpha_{t})X_{t}^{2}$

    //incremental mean

$\hat{X_{t+1}}\leftarrow s_{1}$

    //incremental standard deviation

$\hat{\alpha_{t+1}}\leftarrow\sqrt{s_{2}-s_{1}^{2}}$

The parameters of the anomaly detection algorithm consist of $X_{t}$ the
current data, $\mu_{t}$ the mean of the data, $\hat{X_{t}}$ is the mean of the
data, $\hat{\alpha_{t}}$ the current standard deviation, $P_{t}$ the
probability density function, $\hat{X_{t+1}}$ the mean of the next data
(incremental aggregate), $\hat{\alpha_{t+1}}$ the next standard deviation
(incremental aggregates), $T$ the data size, and $t$ a point in $T$.
Initialize the process by setting the initial data for training the model
$s_{1}=X_{1}$ and $s_{2}=X_{1}^{2}$.

Our work made use of the following parameters $\alpha=0.98,\beta=0.98$, and
$\tau=0.0044$. These thresholds are chosen based on the criteria that outliers
are $\geq 3$ times the standard deviation in normally distributed data.

###  3.2. Online Covariance matrix

  1. (1)

Estimate covariance matrix for initial data, $X\in R^{n\times m}$. Initial
covariance matrix, $C$ where $C\in R^{n\times m}$, $n$ is the number of
samples, $m$ is the number of dimensions as shown in Equation 3.

(3) $C=X*{X}^{T}$

  2. (2)

Perform Cholesky factorization on the initial covariance matrix (positive-
definite), $C$ as shown in Equation 4.

(4) $C_{t}=A_{t}*{A_{t}}^{T}$

  3. (3)

The updated covariance in the presence of new data is equivalent to the
weighted average of the past covariance without the updated data and
covariance of the transformed input as shown in Equation 5.

(5) $C_{t+1}=\alpha*C_{t}+\beta*v_{t}*{v_{t}}^{T}$

Where $v_{t}=A_{t}*z_{t}$ and $z_{t}\in R^{m}$ is understood in our
implementation is the current data. $\alpha$ and $\beta$ are positive scalar
values.

  4. (4)

Increment the Cholesky factor of the covariance matrix as shown in Equation 6.

(6)
$A_{t+1}=\sqrt{\alpha}*A_{t}+\frac{\sqrt{\alpha}}{\Big{\|}z_{t}\Big{\|}^{2}}*\left(\sqrt{1+\frac{\beta*\Big{\|}z_{t}\Big{\|}^{2}}{\alpha}}-1\right)*v_{t}*z_{t}$

  5. (5)

There are difficulties with setting the values of $\alpha$ and $\beta$
respectively. $\alpha+\beta=1$ as an explicit form of exponential moving
average coefficients. The author chose to set the values of $\alpha$, $\beta$
using the statistics of the data stream as shown in Equation 7.

The parameters are set as $\alpha={C_{a}}^{2}$, $\beta=1-{C_{a}}^{2}$ and $n$
is the size of the original data in the static settings. Where
${C_{a}}=\sqrt{1-C_{cov}}$ and $C_{cov}=\frac{2}{{n^{2}}+6}$.

(7)
$A_{t+1}={C_{a}}*A_{t}+\frac{{C_{a}}}{\Big{\|}z_{t}\Big{\|}^{2}}*\left(\sqrt{1+\frac{(1-{C_{a}}^{2})*\Big{\|}z_{t}\Big{\|}^{2}}{{C_{a}}^{2}}}-1\right)*v_{t}*z_{t}$

###  3.3. Online Inverse Covariance matrix

  1. (1)

Estimate covariance matrix for initial data, $X\in R^{n\times m}$. Initial
covariance matrix, $C$ where $C\in R^{n\times m}$, $n$ is the number of
samples, $m$ is the number of dimensions as shown in Equation 8.

(8) $C=X*{X}^{T}$

Inverse the covariance matrix, $C^{-1}$ as shown in Equation 9.

(9) $C^{-1}=\left(X*{X}^{T}\right)^{-1}$

  2. (2)

Perform Cholesky factorization on initial covariance matrix, $C$ as shown in
Equation 10.

(10) $C_{t}=A_{t}*{A_{t}}^{T}$

  3. (3)

Increment the Cholesky factor of the covariance matrix

(11) $C_{t+1}^{-1}=(\alpha*C_{t}+\beta*v_{t}*{v_{t}}^{T})^{-1}$

(12)
$C_{t+1}^{-1}=\alpha^{-1}*(C_{t}+\frac{\beta*v_{t}*{v_{t}}^{T}}{\alpha})^{-1}$

Let us fix, $\hat{v_{t}}=\frac{\beta*v_{t}}{\alpha}$. The resulting
simplification using Sherman$-$Morrison Formula reduces the expression to

(13)
$C_{t+1}^{-1}=\frac{1}{\alpha}*\left({{C_{t}}^{-1}}-\frac{{{C_{t}}^{-1}}*\hat{v_{t}}*{v_{t}}^{T}*{{C_{t}}^{-1}}}{1+(\hat{v_{t}}*{{C_{t}}^{-1}}*{v_{t}}^{T})}\right)$

###  3.4. Online Multivariate Anomaly Detection

The probability density function utilizes ideas from hypothesis testing for
Deciding on a threshold to set the confidence level for determining the
acceptance and rejection regions of the Gaussian distribution curve.

  1. (1)

Use the covariance matrix, $C_{t+1}$ and inverse covariance matrix,
${C_{t+1}}^{-1}$.

  2. (2)

We increment the mean vector, $\mu$ as new data arrives. It is possible to
simplify the Covariance matrix, $C$, which will capture a number of the
dynamics of the system. Let $n$ represent the current count of data before new
data has arrived. Also, $\hat{x}$: is the new data, $\mu_{t+1}$: moving
average as shown in Equation 14.

(14) $\mu_{t+1}=\frac{(n*\mu_{t})+\hat{x}}{n+1}$

  3. (3)

Set a threshold to determine the acceptance and rejection regions. Items in
the acceptance region are considered to be normal behavior as shown in
Equation 15.

(15)
$p(x)=\frac{1}{\sqrt{(2\pi)^{m}|C|}}\exp\left(-\frac{1}{2}(x-\mu)^{T}{C}^{-1}(x-\mu)\right)$

Where $\mu$ is mean vector, $C$ is covariance matrix, $|C|$ is the determinant
of $C$ matrix, $x\in R^{m}$ is data vector, and $m$ is the dimension of $x$
respectively.

##  4\. Experiment

Furthermore, we have provided a set of detailed experiments on the proposed
algorithms in different realistic scenarios. However, we maintain the
statistical framework provided by the work (Carter and Streilein, 2012) with
theoretical guarantees.

We have experimented to determine the usefulness of our algorithm by creating
a simulation with 10000000 random vectors with dimensions of 15. The repeated
trial shows that our algorithm is not sensitive to initialization seeds and
dimensions of the matrix. This requirement was a deciding factor in the choice
of the evaluation metric as shown in Equation 16. We have provided more
information on the metric in Section 5.

Our experiment will check the effect of varying the size of the initial static
window versus the update window as shown in Subsection 4.1 and Subsection 4.2.

Figure 2. Compare the threshold of static vs incremental impact performance of
anomaly detection (Version 1)

###  4.1. Experiment 1

We evaluate the trade-off between the static window and the update window. The
experiment setup is as follows:

  * •

Split the data into 5 segments train on 1st segment(static), update covariance
on 2nd (online), compare with static covariance, and calculate the error.

  * •

Train on 1, 2 segments (static), update covariance on 3rd (online), compare
with static covariance and calculate the error.

  * •

Train on 1, 2, 3 segments (static), update covariance on 4th (online), compare
with static covariance and calculate the error.

  * •

Train on 1, 2, 3, 4 segments (static), update covariance on 5th (online),
compare with static covariance and calculate the error.

Figure 2 contains the experimental result.

Figure 3. Compare the threshold of static vs incremental impact performance of
anomaly detection (Version 2)

###  4.2. Experiment 2

The experiment setup is as follows:

  * •

Split the data into 5 segments.

  * •

Train on 1st segment(static), update covariance on remaining segments
(2,3,4,5) (online), compare with static covariance and calculate error on
segments (2,3,4,5)

  * •

Train on 1, 2 segments (static), update covariance on remaining segments
(3,4,5) (online), compare with static covariance and calculate error on
segments (3,4,5)

  * •

Train on 1, 2, 3 segments (static), update covariance on remaining segments
(4,5) (online), compare with static covariance and calculate error on segments
(4,5)

  * •

Train on 1, 2, 3, 4 segment(static), update covariance on remaining segments
(5) (online), compare with static covariance and calculate error on segments
(5)

Figure 3 contains the experimental result.

##  5\. Result Analysis

Our random matrices get flattened to a vector and utilized as input. The
length of the flattened vector is used as a normalization factor to make the
loss metric that is agnostic to the dimension of the matrix. The loss function
used in the evaluation is Absolute Average Deviation (AAD) because it gives a
tighter bound on the error compared to mean squared error (MSE) or mean
absolute deviation (MAD) as shown in Equation 16. We take the average of the
residuals divided by the ground truth for every sample in our evaluation set.
If the residual is close to zero, we contribute almost nothing to the measure.
On the contrary, if the residual is high, we want to know the difference from
the ground truth.

(16) $AAD=\sum_{i=1}^{n}\left|\frac{\hat{Y_{i}}-Y_{i}}{Y_{i}}\right|$

Where $\hat{Y_{i}}$ is the predicted value, $Y_{i}$ is the ground truth, and
$n$ is the length of the flattened matrix.

We can observe that building your model with more data in the init (static)
phase leads to lower errors compared to having fewer data in the init phase
and using more of the data for an update. The observation matches our
intuition because when you operate in an online mode, you tend to use smaller
storage space. However, there is still a performance trade-off when compared
to batch mode.

The error at the beginning of our training is significant in both charts. This
insight shows that rather than performing the expensive operation of
converting a covariance matrix to have the property of positive definiteness,
it is better to use random matrices that are positive definite. More data
would help us get to convergence as more data arrives.

The success of the experiments has given us the confidence that our
multivariate anomaly detection algorithm would have similar characteristics to
the univariate case described in the work (Carter and Streilein, 2012).

##  6\. Conclusion

There is no generic anomaly detection that works for every possible task. The
underlying assumption in this work is that the features in use capture
relevant information about the underlying dynamic of the system. Our proposed
implementation is an anomaly detection algorithm for handling multivariate
streams even with challenging shifts. In future work, we will make extensions
to support non-stationary distributions in multivariate data streams.

## References

  * (1)
  * Box and Jenkins (1976) George E. P. Box and Gwilym M. Jenkins. 1976.  _Time Series Analysis: Forecasting and Control_.  Holden-Day.
  * Carter and Streilein (2012) Kevin M. Carter and William W. Streilein. 2012.  Probabilistic reasoning for streaming anomaly detection. In _Proceedings of the Statistical Signal Processing Workshop_. 377–380.
  * Ditzler and Polikar (2013) Gregory Ditzler and Robi Polikar. 2013.  Incremental Learning of Concept Drift from Streaming Imbalanced Data.  _IEEE Transactions on Knowledge and Data Engineering_ 25, 10 (2013), 2283–2301.
  * Hodge and Austin (2004) Victoria Hodge and Jim Austin. 2004.  A Survey of Outlier Detection Methodologies.  _Artificial Intelligence Review_ 22, 2 (2004), 85–126.
  * Igel et al. (2006) Christian Igel, Thorsten Suttorp, and Nikolaus Hansen. 2006\.  A computational efficient covariance matrix update and a (1+1)-CMA for evolution strategies. In _Proceedings of the Genetic and Evolutionary Computation Conference, GECCO_. Washington, USA.
  * Muth (1960) John F. Muth. 1960\.  Optimal Properties of Exponentially Weighted Forecasts.  _J. Amer. Statist. Assoc._ 55, 290 (1960), 299–306.
  * Roberts (1982) S. A. Roberts. 1982\.  A General Class of Holt-Winters Type Forecasting Models.  _Journal of Management Science_ 28, 7 (1982), 808–820.
  * Salehi and Rashidi (2018) Mahsa Salehi and Lida Rashidi. 2018.  A Survey on Anomaly Detection in Evolving Data: With Application to Forest Fire Risk Prediction.  _SIGKDD Explorations Newsletter_ 20, 1 (2018), 13–23.
  * Yule (1927) G. Udny Yule. 1927\.  On a Method of Investigating Periodicities in Disturbed Series, with Special Reference to Wolfer’s Sunspot Numbers.  _Philosophical Transactions of the Royal Society of London_ 226 (1927), 267–298.