\section{Public Repository}

The dataset, code, and documentation for Exathlon\ are publicly available at \url{https://github.com/exathlonbenchmark/exathlon}.

\section{Details on Data Collection}
\label{appendix:data}

In this section, we provide additional details on our anomaly design and data collection.

\begin{figure*}[t]
	\centering
	\begin{tabular}{lcc}
		\hspace{-0.1in}
		\subfigure[\small{Metrics in an undisturbed trace}]
		{\label{fig:normal}\includegraphics[width=0.33\textwidth,height=3.6cm]{figures/undisturbed.png}}
		&
		\hspace{-0.1in}
		\subfigure[\small{Bursty input}]
		{\label{fig:burstyinput2}\includegraphics[width=0.33\textwidth,height=3.5cm]{figures/bursty_input.png}}
		&
		\hspace{-0.1in}
		\subfigure[\small{Bursty input until crash}]
		{\label{fig:burstycrash}\includegraphics[width=0.33\textwidth,height=3.5cm]{figures/bursty_input_crash.png}}
		\\
		\hspace{-0.2in}
		\subfigure[\small{Stalled input}]
		{\label{fig:stalledinput}\includegraphics[width=0.33\textwidth,height=3.5cm]{figures/stalled_input.png}}
		&
		\hspace{-0.1in}
		\subfigure[\small{CPU contention}]
		{\label{fig:cpucontention}\includegraphics[width=0.33\textwidth,height=3.5cm]{figures/cpu_contention.png}}
	\end{tabular}
\vspace{-0.2in}
\caption{\small Metrics observed in normal data and anomaly instances (a pair of red vertical bars marks a root cause event) }
\vspace{-0.1in}
\end{figure*}

\subsection{Spark Workload}
\label{appendix:spark-applications}

Our Spark workload consists of 10 stream processing applications that analyzed user click streams from the WorldCup 1998 website~\cite{LiMD+12} (replicated with a scale factor for our long-running applications). The applications processed data formatted as \texttt{(user\_id - - timestamp - - url)} records, using a batch interval of 5 seconds. Join data stored in HDFS was also available for the applications to use, in the form of \texttt{(page\_rank - - url)} records.

We summarize the operation performed by each application below, and provide their source code under \url{https://github.com/exathlonbenchmark/exathlon/tree/master/apps} for more details:

\noindent \textbf{App 1.} GROUP BY + COUNT last batch values.

\noindent \textbf{App 2.} GROUP BY + COUNT + FILTER values from the start.

\noindent \textbf{App 3.} GROUP BY + COUNT + FILTER values in a 30s sliding window (with 20s steps).

\noindent \textbf{App 4.} GROUP BY values in a 30s sliding window.

\noindent \textbf{App 5.} JOIN + GROUP-BY + SUM values in a 30s sliding window.

\noindent \textbf{App 6.} JOIN + GROUP-BY + SUM values in a 30s sliding window.

\noindent \textbf{App 7.} GROUP-BY + SIMULATED UDF for values in a 30s sliding window.

\noindent \textbf{App 8.} GROUP-BY + COUNT + FILTER values in a 10s jumping window.

\noindent \textbf{App 9.} GROUP-BY values in a 10s jumping window.

\noindent \textbf{App 10.} GROUP-BY + COUNT + FILTER values in a 10s jumping window.

The results for these operations were entirely saved on HDFS for applications 6, 8, 9 and 10. For other applications, the results were saved for the first batch only.

Overall, these applications commonly involved group-by aggregations (counts or sums), and varied in terms of %
window type (sliding vs. jumping), parameters (size and slide) and filtering conditions.
Two of the applications additionally involved join operations, and one of them used a user-defined function. As such, our workload aimed to capture a diverse mix of popular streaming primitives.

\subsection{Undisturbed and Disturbed Traces}
\label{appendix:traces}

\noindent
{\bf Undisturbed Traces.}
To collect traces characterizing the normal execution behavior of our Spark cluster, we continually ran experiments (5/10 randomly selected applications at a time) for a month, using input rate values and Spark parameters that suited the capacity of our cluster. Some traces affected by cluster downtime were manually pruned, leaving 59 undisturbed traces, around 15.3GB of data, to constitute the main bulk of our dataset. We consider these traces to be representative of the normal state because the cluster was not interrupted or stress tested by disruptive external events during their recording.

It is important to note that even undisturbed traces exhibit a great deal of variability in the recorded metrics.
For example, Figure~\ref{fig:normal} shows that both the {\em scheduling delay} (delay between the scheduling of a task and the start of its processing) and {\em processing time} (time taken to process an input batch of streaming data) of the recorded application could rise high beyond their usual ranges of values, while its total {\em number of processed records} steadily increased over time.
The spikes observed for the first two metrics were likely caused by other system operations, e.g., checkpointing in Spark or CPU usage by a DataNode in HDFS (Hadoop File System). Since such variations appear in almost every trace, they should be considered to be part of the normal state.

\noindent
{\bf Disturbed Traces.}
Disturbed traces were obtained by introducing anomalous events during an execution. Based on discussions with industry contacts from the Spark ecosystem, we came up with 6 types of anomalous events. When designing these, we considered that:
(i) they lead to a visible effect in the trace,
(ii) they do not lead to an instant crash of the application (since AD would be of little help in this case),
(iii) they can be tracked back to their root causes.

\minipb{Bursty Input (Type 1):}
For every application, the user expects a certain range of data input rates and configures Spark parameters to allocate sufficient resources for such input rates.  However,  on some special occasions (e.g., a special sales day for a retailer) the input rate can  increase significantly, and existing resources be insufficient for handling the higher data rate. Such phenomena can be reflected in the data, e.g., in increased values of processing time, memory usage, and scheduling delay, as shown in Figure~\ref{fig:burstyinput2}. It is because each batch of data (e.g., received in the past 10 seconds) increases dramatically in size, and the processing time of the batch increases accordingly. When the processing time exceeds the batch interval (e.g., 10 seconds), the application cannot process data fast enough; so data is put in memory by the receivers and such data experiences a higher delay before it is scheduled for processing. Early detection of bursty input is helpful because it can lead to corrective actions such as allocating more resources to the application.

To mimic input rate spikes, we ran a disruptive event generator (DEG) on the Data Senders to temporarily increase the input rate by a given factor for a duration of 15-30 minutes and then set the input rate back to its normal range. We repeated this pattern multiple times during a given trace, with sufficient gap between two instances. As such, we created a total of 29 instances of this anomaly type over 6 different traces.

\minipb{Bursty Input Until Crash (Type 2):}
Spark developers rank the out of memory (OOM) condition to be the number one reason for Spark application failures. To capture such phenomena, we extended the bursty input design:  instead of reducing the input rate back to a normal range, we kept the input at the same high rate until the application eventually crashed. Figure~\ref{fig:burstycrash} shows how the system behaviors during such an event are reflected in the data: the application starts with a normal input rate. At some point, the input rate rises high, and the processing time and scheduling delay build up until one of the executors fails due to an out-of-memory error.
A new executor is then launched to replace the failed one, but the sustained high rate causes more executors to fail. When the number of failed executors reaches a threshold, the application is killed by Spark, constituting an application crash.

To instrument this anomaly type, we ran the DEG  forever, first crashing the executors due to lack of memory, and eventually crashing the application. We injected this anomaly into 7 different traces.

\minipb{Stalled Input (Type 3):}
Another type of input-related anomaly is stalled input, i.e., no data is injected into the Spark application. This may  indicate a failure of the data source (e.g., Kafka or HDFS) and results in a waste of resources as the driver and executors are still running with no data to process. As shown in Figure~\ref{fig:stalledinput}, this type of anomaly can  be reflected in our collected metrics: since no data is coming from the receiver, the number of processed records, as well as other related metrics, do not increase, and the processing time is much lower than in the normal state.

To generate these anomalies, we ran a DEG that set the input rates to 0 for about 15 minutes, and then periodically repeated this pattern every few hours. This gave us a total of 16 anomaly instances across 4 different traces..

\minipb{CPU Contention (Type 4):}
A Spark cluster is usually run with a resource manager (e.g., YARN) that allocates resources to each application so that each of them has exclusive access to its own resources, including CPU and memory. However, the resource manager cannot prevent external programs from using the resources that it allocated previously, which is a common reason for performance issues in distributed environments with high concurrency. For instance, it is not uncommon to have a Hadoop DataNode using a large amount of CPU on a node also used by a Spark application. This generates a contention for resources that can often slow down the progress of the Spark application.

We reproduced this type of anomaly using a DEG that ran Python programs to consume all CPU cores available on a given Spark node. We created 26 such anomaly instances over 6 different traces. An example trace containing multiple instances of this type of anomaly is shown in Figure~\ref{fig:cpucontention}.
From this, we see there are several intervals during which the processing time is higher than normal due to  CPU contention on the affected node, meaning that one or more executors are using fewer resources than allocated. If the processing time is still lower than the batch interval, the CPU contention does not affect much the progress of the application;  otherwise, scheduling delay starts to build up. Besides processing time, CPU contention can have a more severe impact on the Spark application: if the driver is performing a lot of computation and CPU contention occurs on the driver node, it can cause it to crash.
Then the application master has to restart the driver, which will itself restart all the executors.

\minipb{Driver Failure (Type 5) and Executor Failure (Type 6):}
Other types of anomaly include abrupt failures in distributed systems. A common example is a node failure caused  by a  hardware fault or a maintenance operation. In this case, all the processes (drivers and/or executors) located on that node become unreachable. Such processes must be restarted on another node, which causes processing delays.
We created such anomalies by failing driver processes, causing cumulative metrics to drop back to their initial values and the application to stop until its driver is restarted in about 20 seconds.
We also created anomalies by failing executor processes, which get restarted 10 seconds after the failure, but whose effects on metrics such as processing delay may continue longer.
We created 9 driver failures and 10 executor failures over 11 different traces.

\subsection{Extended Effect Intervals}
\label{appendix:extended-effect}

In this section, we describe the way that we set the extended effect interval (see Table~\ref{tab:traces}) for each anomalous event.
Each anomaly instance was initially labeled with its known type and {\em root cause interval} (RCI).
However, we observed that some injected events could have a long-lasting effect on relevant metrics such as processing time and scheduling delay even after their end time.
To characterize such events, we therefore sought to extend their root cause interval to include this long-lasting effect.
More specifically, we used domain knowledge to set an {\em extended effect interval} (EEI) for each anomalous event, starting immediately after the end of its RCI, and ending
at a point after which we deem anomaly detection not helpful.
In practice, we either set the end of an EEI to when the application had fully restarted, or to when its main metrics had come back to normal. Here are the rules we used for each type of anomalous event:

\begin{enumerate}
\item Bursty Input: The end of the EEI is set to the point when highly related metrics such as the processing time and scheduling delay come back to normal.

\item Bursty Input Until Crash: The EEI is set to null, because the root cause event already ends at the time of the application crash.

\item Stalled Input: The end of the EEI is set to the point when the processing time comes back to normal (the application restarts processing data at its usual rate).

\item CPU Contention:
a) If the effect is increased processing time, we set the end of the EEI to the time when processing time and scheduling delay come back to normal.
b) If the effect is a driver crash, we set the end of the EEI to the time when the application restarts (typically 1 minute after the crash in practice).

\item Driver Failure:
The end of the EEI is set to when the application restarts.

\item Executor Failure:
If the executor is replaced, we observe an increase in scheduling delay during the replacement time. If it is not replaced, then the buildup in scheduling delay either lasts temporarily or until the application crashes. In all these cases, we set the end of the EEI to the point when the scheduling delay comes back to normal.
\end{enumerate}

\noindent The combined anomaly interval, RCI plus EEI, gives our benchmark more freedom to evaluate anomaly detection algorithms.

\section{Details on Benchmark Design}
\label{appendix:benchmark}

\begin{table}[t]
\centering
{\small
\begin{tabular}{|l|c|c|c|c|c|c|}
\hline
& \multicolumn{3}{c|}{\bf Precision} & \multicolumn{3}{c|}{\bf Recall} \\
& \multicolumn{3}{c|}{\bf Parameters} & \multicolumn{3}{c|}{\bf Parameters} \\
\hline
{\bf AD Functionality Level} & {\bf $\alpha$} & {\bf $\delta$} & {\bf $\gamma$} & {\bf $\alpha$} & {\bf $\delta$} & {\bf $\gamma$} \\
\hline
\hline
AD1: Anomaly Existence & 0 & Flat & 1 & 1 & N/A & N/A \\
\hline
AD2: Range Detection & 0 & Flat & 1 & 0 & Flat & 1 \\
\hline
AD3: Early Detection & 0 & Flat & 1 & 0 & Front & 1 \\
\hline
AD4: Exactly-Once Detection & 0 & Flat & 0 & 0 & Front & 0 \\
\hline
\end{tabular}
}
\caption{\small Range-based precision/recall parameter settings}
\label{tab:params}
\end{table}

\subsection{Anomaly Detection (AD) Functionality}
\label{appendix:ad-functionality}

\noindent
\textbf{Range-based Precision and Recall.}
Table \ref{tab:params} further shows how to set the tunable parameters of the range-based precision/recall framework~\cite{TatbulLZAG18} to capture the requirements of AD1-AD4 in general.
We provide a brief insight for these parameters here, and refer the reader to the original paper for further details \cite{TatbulLZAG18}.
The $\alpha$ parameter is to reward detecting the existence of an anomaly range and is only relevant here in the context of recall when we only care about anomaly existence. Therefore, it always gets a 0 value except when computing AD1's recall.
The $\delta$ parameter captures the positional bias, i.e., the reward to be earned from which portion of the anomaly range was successfully detected. It is the most relevant for AD3 and beyond, where detection latency is considered. For recall of AD3-4, $\delta$ should be set to ``Front'' to favor early detection, and ``Flat'' in all other cases, where detection position does not matter.
Finally, the $\gamma$ parameter is to penalize fragmented ranges (i.e., an anomaly range detected as multiple subranges) and is, therefore, directly relevant to exactly-once detection. For AD4, we set $\gamma$ to 0, indicating that fragmentation is strictly undesirable; for all others, it has no penalizing effect to the score.
We would like to note that there is a fourth parameter, $\omega$, which is to compute how much reward is earned from the size of the anomaly range that is correctly detected (relevant for AD2 and beyond). We use its default, additive definition in the original model with a minor normalization adjustment to ensure monotonicity. More specifically, we make the precision/recall scores decrease from AD2 to AD3 by defining the union of detected ranges under AD2 as the best positioned (i.e., the earliest) one under AD3 to achieve a monotonic behavior.

\noindent
\textbf{Learning Settings.}
By default, our benchmark offers only undisturbed traces  as training data for building an AD model, because it is the most practical for real-world use.
This can be seen as a {\em ``noisy'' semi-supervised anomaly detection} problem~\cite{ChandolaBK2009}, meaning that the training data is assumed mostly normal but subject to small amounts of noise.
In this context, Exathlon\ also tests how well a learned AD model generalizes to different workload characteristics, including the Spark application ($A$), input rate ($R$), and concurrency ($C$) characteristics.
We offer four alternative learning settings, considering the following two aspects:

\noindent
\underline{\emph{Modeling Subject (1-App vs. N-App Learning)}:} Spark applications differ in workload characteristics (e.g., CPU intensive versus I/O intensive). Such characteristics may lead to different runtime observations (e.g., a CPU-intensive application would be more sensitive to CPU contention anomalies). Hence detecting anomalies in traces of different applications requires more generalization power.
(i)~\emph{1-App learning} focuses on training and evaluating an AD model on a single application basis. As such, there is no need for the model to generalize to other applications.
(ii)~\emph{N-App learning} trains an AD model using multiple applications together, and hence requires it to jointly learn how to characterize different applications.

\noindent
\underline{\emph{Training Constraints (Many vs. Few Examples)}:}
Each of our traces contains data characterizing an application running with a given input rate in a given concurrent environment.
An AD model trained on undisturbed traces with certain ($A$, $R$, $C$) settings may later be subject to a new trace with a previously unseen ($A$, $R$, $C$) setting.
A well-learned model is supposed to generalize across these differing workload settings.
However, such generalization power is easier to achieve if the training data includes many examples of ($A$, $R$, $C$),
or many ($R$, $C$) examples in the {\em 1-App} learning setting.
By default, Exathlon\ reserves the disturbed traces entirely for testing.
As such, the training data includes only the 59 undisturbed traces, which is unlikely to cover most ($A$, $R$, $C$) values and hence constitutes the {\em Few Example}-training scenario. This  bears similarity with the few-shot learning problem~\cite{WangYKN20} recently studied in the ML community, and poses a great challenge for learning.
Exathlon\ also offers a simpler, yet less realistic, training setting, {\em Many Example}-training, in case some algorithms cannot learn from limited examples. Here, we allow the AD algorithm to also leverage an earlier segment (with normal data only) of each disturbed trace in training, while testing on the anomalies in the later segments of the same traces. This way, the model is given a chance to ``peek'' at the normal state (including all workload characteristics) of a particular test trace.
This simulates the scenario of a model trained on many ($A$, $R$, $C$) examples, some of which similar to the current test trace.

By combining the above options, we obtain four learning settings to experiment with: {\em LS1: 1-App, Many-Examples}; {\em LS2: N-App, Many-Examples}; {\em LS3: 1-App, Few-Examples}; and {\em LS4: N-App, Few-Examples}.

\subsection{Broader Applicability of AD Evaluation}
\label{appendix:ad-applicability}

\textbf{Applicability to Other AD Techniques.}
In its current form, the Exathlon benchmark considers all techniques that work by 1) training a model for data normality via learning a specific task on undisturbed traces, 2) assigning real-valued outlier scores to new test records as proportional to the error made by the model, and 3) deriving binary predictions for these records based on whether their scores exceeded a given threshold. In our experimental study, this third step was performed in the same way for the three compared methods (LSTM, AE, and BiGAN).

In practice, any technique that can assign an outlier score to each record of a test sequence based on a model learned from (mostly) normal data can directly be evaluated using this benchmark. Depending on the feature extraction step, the range of directly applicable techniques is hence quite broad, ranging from deep learning-based techniques like the ones presented in this paper to simpler statistical techniques~\cite{bianco01}, to other recent techniques for subsequence AD~\cite{seq2graph}, \cite{norm} and data drift detection~\cite{conf-constraints}.

Although not considered in this paper, our benchmark could also be used for problems focusing on sliding windows in disturbed traces only, with or without pre-training a model on undisturbed traces. This would extend the range of applicable techniques to also include the ones of the distance-based outlier detection for data streams (DODDS) literature~\cite{Tran2015}. In such streaming scenarios, where both CPU time and memory usage can be of primary concern, the maximum memory used to update the outlier information could be reported along with the inference time as part of our AD efficiency metrics, since this metric is frequently encountered in the DODDS literature.

\textbf{Applicability to Other Datasets.}
The AD metrics of Exathlon could in principle be used with any labeled time series anomaly test datasets similar to ours. For datasets where ground truth anomalies are provided as ranges, like those most frequently found in the discord discovery literature (e.g., the simulated engine disks (SED) dataset and the MIT-BIH Supraventricular Arrhythmia Database (MBA) recently used in \cite{seq2graph}), the four AD levels of the benchmark would be directly usable. For datasets where labels are provided as points, our benchmark would still be usable by resorting to the classical definitions of precision and recall metrics, which are special cases of the range-based ones \cite{TatbulLZAG18}. The availability of labels for multiple anomaly types (e.g., as in MBA \cite{physiobank}, \cite{mba-impact}) would also enable to directly report type-wise recall scores in addition to the global ones (where all anomaly types are considered the same).

Data granularity is use case dependent. That is, trace/application-level data separation may not make sense for all use cases. Depending on the use case, there may be other granularity levels or only a single one. Our AD evaluation metrics and methodology are equally applicable to such use case scenarios.

While we have not applied our AD metrics to other anomaly datasets beyond our Spark-based traces, the range-based precision/recall framework constituting their foundation was effectively tested on multiple real and artificial time series anomaly datasets from the AD literature \cite{TatbulLZAG18}. In fact, this framework is a generalization of the classical precision/recall metrics from point-based data to range-based data. As such, it is a superset that encompasses the classical precision/recall metrics and other similar metrics based on them (e.g., F-score), which are widely used in the literature for AD evaluation. As was shown in \cite{TatbulLZAG18} based on real datasets from the Numenta Anomaly Benchmark \cite{numenta:icmla15}, the range-based metrics can in fact be tuned to simulate the behavior of both the point-based metrics and Numenta's AD scoring metrics, as well as capturing additional behaviors that different AD applications may require. The range-based metrics also provide a modular and extensible formulaic structure. As such, we believe that our choice of AD metrics both provide good coverage of the existing metrics in the literature as well as being amenable to future extensions if necessary.

\subsection{Broader Applicability of ED Evaluation}
\label{appendix:ed-applicability}

\textbf{Applicability to Other ED Works.}
While a user study may be the best way to evaluate the usefulness of explanations, it is not always available and may come at a high cost. Therefore, our benchmark aims to provide {\em automated} evaluation of ED methods based on intuitive metrics. For selecting those, we surveyed multiple other ED works and summarized the main metrics they used (see Table \ref{tab:survey_measurement}). Among these, we found conciseness and accuracy to be the most general and intuitive: human users typically prefer small sets of important factors as explanations, and if these were to be used as an AD tool, they should ideally achieve good accuracy. We also added the concept of consistency, given that accuracy is not always applicable or meaningful for a given ED method. As an example, running LIME to explain the outlier score assigned to a given test sample might not be directly usable for binary prediction. Another issue might be the impact of hard-coded constants in explanation predicates. For example, a valid explanation to a human user could be that ``the CPU usage is too high''. However, this notion of ``too high'' might be a relative concept, meaning above 50\% in a given context but above 80\% in another. A lot of ED methods typically include such constants in their explanations, sometimes leading a ``good'' instance explanation to have poor accuracy when applied to other instances. Like conciseness, the concept of consistency has the benefit to rule out the effect of such constants in the evaluation of an ED method. Although the concept of these three metrics is common to all evaluated methods, the way we calculate them typically depends on how the method works. For instance, if the method takes a dataset as input to generate the explanation, like in~\cite{ZDM2017, macrobase:sigmod17, LakkarajuBL16, WuHPZ2018}, we consider the consistency implementation described in Section~\ref{sec:benchmark}. For other methods, e.g., ~\cite{Ribeiro0G16, Ribeiro0G18, LundbergL17}, the calculation of consistency is window-based, implemented for LIME~\cite{Ribeiro0G16} in our pipeline. If needed, other implementations can also be added by the user.

\textbf{Applicability to Other Datasets.}
The ED metrics of Exathlon\ are directly applicable to the datasets used in~\cite{ZDM2017, macrobase:sigmod17}. Our metrics can also be applied for evaluating explanations of classification results, which is relevant for the datasets of~\cite{LakkarajuBL16, WuHPZ2018, Ribeiro0G18}, balanced and used for classification problems. For the text and image datasets used in ~\cite{Ribeiro0G16, LundbergL17}, our metrics are applicable after an additional data preprocessing step. If an application scenario does not provide any detailed anomaly type information, one can consider all anomalies as of the same type. Finally, if there is no clear separation of a dataset into instances, the user may need to decide to consider either ED1 or ED2 based on the dataset semantics.

\begingroup
\setlength{\tabcolsep}{1pt}
\begin{table}[t]
	\centering
	\footnotesize
	\begin{tabular}{|c|c|c|c|}
		\hline
		 ED Work & Metrics & Applicability of Exathlon\ & Comment \\ \hline \hline
		\multirow{3}{*}{{EXstream~\cite{ZDM2017}}} &
		{Conciseness} & {\textbf{Yes}} & {ED2 Conciseness} \\
		&  {Consistency} & {No} & {Needs ground truth}  \\
		& {Accuracy} & {\textbf{Yes}} & {ED2 Accuracy} \\ \hline
		{{MacroBase~\cite{macrobase:sigmod17}}} &
		{Efficiency} & {\textbf{Yes}} & {Running time}  \\ \hline
		\multirow{2}{*}{{LIME~\cite{Ribeiro0G16}}} &
		{Coverage of gold set} & {No} & {Model-dependent}  \\
		& {Overlap with noise} & {No} & {Model-dependent}  \\ \hline
		\multirow{2}{*}{{Anchors~\cite{Ribeiro0G18}}} &
		{Coverage} & {No}& {Algorithm-specific}  \\
		& {Precision} & {\textbf{Yes}} & {ED1 Accuracy}  \\  \hline
		\multirow{3}{*}{{Decision S.~\cite{LakkarajuBL16}}} &
		{F1-score, AUROC} & {\textbf{Yes}} & {ED2 Accuracy} \\
		& {Fraction of 3 types} & {No} & {Algorithm-specific}  \\
		& {Rule width/\# of rules} & {\textbf{Yes}} & {ED2 Conciseness} \\ \hline
		{{SHAP~\cite{LundbergL17}}} &
		{Visually on MNIST} & {No} & {Visual inspection}  \\ \hline
		{{Tree Reg.~\cite{WuHPZ2018}}} &
		{AUROC} & {\textbf{Yes}} & {ED2 Accuracy}  \\ \hline
	\end{tabular}
	\caption{\small Metrics used in other ED works and applicability of the Exathlon\ metrics}
	\label{tab:survey_measurement}
\end{table}
\endgroup

\section{Details on Pipeline Design}
\label{appendix:pipeline}

\minipa{1. Data Partitioning.}
This initial phase takes as input the 93 raw traces, as described in \S \ref{sec:anomalies}. It first performs simple data cleaning, e.g., replacing missing data with a default value. It then performs data selection and partitioning of the 93 traces according to the selected learning setting among the ones described in \S \ref{sec:benchmark}:
(i) In the {\em 1-App} learning settings (LS1 and LS3), we consider each Spark application, $i$, separately. The data selection step collects all the traces relating to $i$. The data partitioning step then takes all undisturbed traces of $i$ as training data, $D_{\text{train}}$, and all disturbed traces of $i$ as test data, $D_{\text{test}}$. These two datasets will be used to run the full pipeline, learning the model for $i$ and conducting the final evaluation. We repeat this process for each Spark application.
(ii) In the {\em N-App} learning settings (LS2 and LS4), the data selection step collects all 93 traces. The data partitioning step then takes all  undisturbed traces as training data, $D_{\text{train}}$, and all disturbed traces as test data, $D_{\text{test}}$.
(iii) The above implementations correspond to the {\em Few-Examples} learning settings (LS3 and LS4). If the pipeline is configured to run under the  {\em Many-Examples} settings (LS2 and LS4), we further augment $D_{\text{train}}$, for {\em 1-App} and {\em N-App} learning, respectively, with an earlier segment of each relevant test trace (mostly normal data).

\minipa{7. ED Evaluation.}
After processing each test trace, we obtain a set of anomalies with their corresponding explanations. We then collect the explanations from all the test traces to run the final ED evaluation.
To compute our metrics, we implemented for each ED method:
(i) The extraction function $G_\mathbf{A} (F_{t,w})$, used to extract the feature set appearing in a given explanation, required for conciseness and consistency. We support extraction functions for common forms of explanations, including logical formulas, linear models, and decisions trees.
(ii) The subsampling procedure, applied to each anomaly $X_{t,w}$  when computing its ED1 stability and accuracy. If an ED method can explain an anomaly of arbitrary duration $w$, our sampling procedure performs a random split so that 80\% of the anomaly is used to construct its explanation (and its remaining 20\% is used along with neighboring normal data to compute its ED1 accuracy). For ED methods that can only explain anomalies of fixed size $s$, e.g., LIME~\cite{Ribeiro0G16}, we create samples of size $s$ such that their union entirely covers $X_{t,w}$.
The implementation of concordance and ED2 accuracy directly follows their definitions.

\section{Details of Our Experimental Study}
\label{appendix:experiments}

All of our experiments were performed in a cluster of 20 compute nodes, each with 2 Intel{\scriptsize{\textsuperscript{\textregistered}}} Xeon{\scriptsize{\textsuperscript{\textregistered}}} Gold 6130 16-core processors, 768GB of memory, and 64TB disk.

\subsection{Customized Feature Set}
\label{appendix:custom-set}

Our custom feature set was produced via the manual selection of relevant features resulting from domain knowledge and various visualizations. This feature set characterizes different aspects of a running application, with 9 features relating to the Spark application driver, 6 to the Spark executors, and 4 to the OS metrics. We give the complete list of features below:

\begin{itemize}
	\item \textbf{\underline{Driver features}}
	\begin{itemize}
	\item Processing delay, scheduling delay and total delay of the last completed batch.
	\item First-order difference ($f_t := f_{t+1} - f_t$, simply called \emph{difference} in the following) in number of completed batches, received records and processed records since the start of the application.
	\item Difference in number of records in the last received batch.
	\item Difference in total memory used.
	\item Difference in heap memory used by the JVM.
	\end{itemize}
	\item \textbf{\underline{Executor features}} -- All averaged across active Spark executors (before differencing).
	\begin{itemize}
	\item Difference in number of HDFS writing operations.
	\item Difference in CPU time and runtime.
	\item Difference in number of records read and written during shuffling operations.
	\item Difference in heap memory used by the JVM.
	\end{itemize}
	\item \textbf{\underline{OS features}} -- All for each of the four nodes of the mini-cluster the applications were run on.
	\begin{itemize}
	\item Difference in global CPU percentage idle time.
	\end{itemize}
\end{itemize}

In the following, we list the feature names with their indices, as used in the explanations shown in Figure~\ref{fig:ed_example_all}.

\begin{enumerate}[label={\arabic*}.]
	\item[0.] driver{\_}Streaming{\_}lastCompletedBatch{\_}processingDelay{\_}value
	\item driver{\_}Streaming{\_}lastCompletedBatch{\_}schedulingDelay{\_}value
	\item driver{\_}Streaming{\_}lastCompletedBatch{\_}totalDelay{\_}value
	\item 1{\_}diff{\_}driver{\_}Streaming{\_}totalCompletedBatches{\_}value
	\item 1{\_}diff{\_}driver{\_}Streaming{\_}totalProcessedRecords{\_}value
	\item 1{\_}diff{\_}driver{\_}Streaming{\_}totalReceivedRecords{\_}value
	\item 1{\_}diff{\_}driver{\_}Streaming{\_}lastReceivedBatch{\_}records{\_}value
	\item 1{\_}diff{\_}driver{\_}BlockManager{\_}memory{\_}memUsed{\_}MB{\_}value
	\item 1{\_}diff{\_}driver{\_}jvm{\_}heap{\_}used{\_}value
	\item 1{\_}diff{\_}node5{\_}CPU{\_}ALL{\_}Idle\%
	\item 1{\_}diff{\_}node6{\_}CPU{\_}ALL{\_}Idle\%
	\item 1{\_}diff{\_}node7{\_}CPU{\_}ALL{\_}Idle\%
	\item 1{\_}diff{\_}node8{\_}CPU{\_}ALL{\_}Idle\%
	\item 1{\_}diff{\_}avg{\_}executor{\_}filesystem{\_}hdfs{\_}write{\_}ops{\_}value
	\item 1{\_}diff{\_}avg{\_}executor{\_}cpuTime{\_}count
	\item 1{\_}diff{\_}avg{\_}executor{\_}runTime{\_}count
	\item 1{\_}diff{\_}avg{\_}executor{\_}shuffleRecordsRead{\_}count
	\item 1{\_}diff{\_}avg{\_}executor{\_}shuffleRecordsWritten{\_}count
	\item 1{\_}diff{\_}avg{\_}jvm{\_}heap{\_}used{\_}value
\end{enumerate}

\subsection{Details of AD Methods}
\label{appendix:ad-methods}

We next describe in more detail the methods we used for the normality modeling, outlier score derivation and threshold selection steps of AD modeling when conducting our experiments.

\textbf{LSTM}. Long Short-Term Memory (LSTM) \cite{DHochreiterS97} is a deep neural network structure designed for handling sequential data. In the context of AD, it is typically used as a forecasting-based method, where the outlier score of a given data point is derived from its forecasting error by the model. Like in a recent work~\cite{BontempsCML16}, our experiments used LSTM for AD by setting the outlier score of a test record as its relative forecasting error. Unlike in~\cite{BontempsCML16}, the produced outlier scores were however not further averaged but kept as is.

\textbf{AE}. An Autoencoder (AE) \cite{HS} is a deep neural network structure designed to map its input data to a corresponding latent representation in lower dimensional space (or \emph{coding}), by learning to accurately reconstruct it. In the context of AD, autoencoders are typically used as reconstruction-based methods, where the outlier score of a given data sample is derived from its reconstruction error by the model. In our experiments, we used a dense autoencoder architecture for AD by setting the outlier score of a test window as its reconstruction mean squared error (MSE), and the outlier score of a test record as the average score of its enclosing windows.

\textbf{BiGAN}. Generative Adversarial Networks (GAN)~\cite{GoodfellowPMXWOCB14} learn to generate samples with similar characteristics as their training data, by training a generator network to map random noise in latent space to samples that appear realistic to a discriminator network. Bidirectional Generative Adversarial Networks (BiGAN)~\cite{DonahueKD17} are used to automatically learn the inverse mapping of the generator, through the additional training of an encoder network, that can be later used in pair with the generator to reconstruct some (real) input samples. As such, BiGAN can also be used in the context of AD as a reconstruction-based method, with an additional flexibility in deriving the outlier score of a test window, by leveraging the architectures of all three networks~\cite{abs180206222}. In our pipeline, we experimented with network architectures allowing the usage of LSTM networks to better handle time sequences, like done in~\cite{abs180904758}. The outlier score of a test window was defined as the average of its MSE by the (encoder, generator) pair and its feature loss by the discriminator, as defined in \cite{abs180206222}. Like for AE, the outlier score of a test record was defined as the average score of its enclosing windows.

\textbf{Threshold Selection}. Throughout our experiments, threshold selection was performed for an AD method based on the distribution of the outlier scores it assigned to the samples of $D_{\text{train}}^2$. More precisely, we experimented with threshold definitions of the form:

\[
\texttt{threshold} = S_1 + c \cdot S_2
\]

Where $S_1$ and $S_2$ are two descriptive statistics of the $D_{\text{train}}^2$ sample outlier scores, and $c$ is a constant we refer to as the \emph{thresholding factor}. We also allowed this computation to be performed twice, removing for the second iteration the outlier scores that were above the threshold found for the first iteration, in order to drop any obvious outliers that could prevent us from finding a suitable threshold.

For every AD evaluation in our experiments, the following threshold selection methods were evaluated for $c \in \left\{1.5, 2, 2.5, 3\right\}$:

\begin{itemize}
    \item \textbf{STD} : $S_1 = \bar m$, $S_2 = s$, \emph{i.e.} sample mean and standard deviation.
    \item \textbf{MAD} : $S_1 = \text{median}, S_2 = \text{MAD}$, where $\text{MAD} = 1.4826 \cdot \text{median}(\lvert X - \text{median}(X) \rvert)$.
    \item \textbf{IQR} : $S_1 = Q_3$, $S_2 = Q_3 - Q_1 = \text{IQR}$.
\end{itemize}

This led to a total of 24 thresholding parameter combinations, corresponding to 24 final detection performances for each evaluated AD method, for which we reported the median values in Table~\ref{tab:custom-ls4-ad-detection}.

\subsection{Details of ED Methods}
\label{appendix:ed-methods}

We also integrated three recent ED methods into our pipeline: EXstream~\cite{ZDM2017} and MacroBase~\cite{macrobase:sigmod17} from the database community for outlier explanation in data streams, and LIME~\cite{Ribeiro0G16}, an influential ED method from the ML community.

{\bf EXstream}~\cite{ZDM2017} takes as input an anomaly dataset and a normal, reference, dataset. It first calculates the entropy-based single-feature reward, which indicates the ability of a feature to separate normal from anomalous data, and then uses several heuristics to find a small set of high-reward features to build the explanation. One of its original steps was false positive filtering, pruning the features only incidentally correlated with abnormality. Since this step required generating extra labeled data based on user annotations, unavailable in our problem setting, we did not include it in our implementation.

{\bf MacroBase}~\cite{macrobase:sigmod17} considers the same input format as EXstream. It first calculates the risk ratio for each feature, and then for each combination of features based on a frequent itemset mining procedure. It also employs several optimization techniques to improve efficiency. %
Since it was originally designed for categorical features, we added the extra step of assigning categorical values to our numerical features via equal width binning in our implementation.

{\bf LIME}~\cite{Ribeiro0G16} was originally designed to explain a given prediction of a classifier, but also provides an interface for regression models and temporal inputs. Internally, LIME first leverages Lasso~\cite{efron_2004} to identify $k$ important features, and then uses a customized loss function to learn a new linear model approximating the original model's output locally. In our implementation, we used LIME's RecurrentTabularExplainer interface, and set $k=5$.

\begin{figure*}[!b]
	\centering
	\begin{tabular}{lc}
		\hspace{-0.2in}
		\subfigure[\small{LSTM record-wise outlier scores}]
		{\label{fig:lm4-lstm-app1-cpu}\includegraphics[width=0.49\textwidth,height=2.0cm]{figures/lm4_lstm_app1_cpu.png}}
		&
		\subfigure[\small{AE record-wise outlier scores}]
		{\label{fig:lm4-ae-app1-cpu}\includegraphics[width=0.49\textwidth,height=2.0cm]{figures/lm4_ae_app1_cpu.png}} \\
	\end{tabular}
\vspace{-0.2in}
\caption{\small Record-wise outlier scores on a T4 trace of Application 1 (FS$_\text{custom}$)}
\label{fig:lm4-app1-scores}
\end{figure*}

\subsection{Additional AD Evaluation Results}
\label{appendix:ad-results}

As shown in Figure~\ref{fig:lm4-app1-scores} for a T4 trace of Application 1, the AE method typically produces {\em smooth} record-wise outlier scores, by taking averages over overlapping windows, while the ones of LSTM often exhibit discontinuous {\em spikes}.
The outlier scores of the last two records were removed for LSTM, as their higher values visually packed down the others.

In the rest of this section, we describe the two additional AD experiments we conducted, respectively relating to the impact of feature extraction and learning setting on the AD performance.

\minipa{Experiment 5 (AD2).}
In this experiment, we study the effect of feature extraction on the performance of the three methods, by comparing our custom feature selection (FS$_{\text{custom}}$) to applying PCA to the raw input features directly (FS$_{\text{pca}}$). From Table~\ref{tab:pca-ls4-ad2-global-separation}, we can see that using FS$_{\text{pca}}$ caused the global separation abilities of all methods to drop significantly (especially for LSTM), largely due to the fact that PCA was not able to preserve the right signals for detecting a major portion of our anomalies.

When using FS$_{\text{pca}}$, we are not averaging across active executor features anymore, which likely plays a part in the fact that separation abilities for T6 were either maintained or increased for AE and BiGAN. Another aspect of FS$_{\text{pca}}$ is that it selects features based on their variance, which means that features tending not to vary much in normal scenarios, like scheduling delay or input rate, will be far less represented in the output space than features having naturally more activity, like CPU- or memory-related metrics. This likely explains why T4 separation ability was also maintained or increased for AE and BiGAN. Like previously hinted by remark R3, we can see this benchmark can be used to analyze the effect of constructing different feature spaces on the overall AD down to the detection of specific anomaly types.

\begin{table}[t]
    \centering
    \small
    \begin{tabular}{|c|c|cccccc|}
    	\hline
    	 Method & Ave & \multicolumn{6}{c|}{AUPRC for Anomaly Types T1$\rightarrow$T6}\\ \hline \hline
   		\textbf{LSTM} & 0.14 & 0.15 & 0.14 & 0.09 & 0.14 & 0.21 & 0.12 \\
   		\textbf{AE} & 0.39 & 0.43 & 0.39 & 0.18 & 0.49 & 0.49 & 0.38 \\
   		\textbf{BiGAN} & 0.34 & 0.40 & 0.29 & 0.21 & 0.47 & 0.37 & 0.29 \\ \hline
    \end{tabular}
    \caption{\small Global separation results (FS$_{\text{pca}}$, AD2)}
    \label{tab:pca-ls4-ad2-global-separation}
\end{table}

\minipa{Experiment 6 (FS$_{\text{custom}}$, AD2).}
In this last experiment, we study the effect of using different learning settings for the three AD methods. Table~\ref{tab:custom-ls1:4-ad2-app-separation} reports their application-level separation abilities for each learning setting. From these, we see that AE and BiGAN followed the expected trend, with their average separation ability decreasing as we increase data variety and disable peeking at the beginning of test traces. On the contrary, LSTM seemed to perform better for LS2 and LS4 ({\em N-App} settings), which seems consistent with Experiment 3 showing that increasing the amount of training data for LSTM increases its performance, even when data variety is also increased. Although the changes in results appear as minor in these examples, we can see that training models with different learning settings can help draw some conclusions regarding their data needs and generalization ability ({\bf R8}).

\begin{table}[t]
    \centering
    \small
    \begin{tabular}{|c|c|c|cccccc|}
    	\hline
    	 LS & Method & Ave & \multicolumn{6}{c|}{AUPRC for Anomaly Types T1$\rightarrow$T6}\\ \hline \hline
   		\multirow{3}{*}{\textbf{LS1}} &
   		\textbf{LSTM} & 0.43 & 0.55 & 0.35 & 0.46 & 0.33 & 0.58 & 0.32 \\
   		& \textbf{AE} & 0.58 & 0.72 & 0.41 & 0.67 & 0.54 & 0.71 & 0.42 \\
   		& \textbf{BiGAN} & 0.55 & 0.82 & 0.35 & 0.37 & 0.53 & 0.77 & 0.47 \\ \hline \hline
   		\multirow{3}{*}{\textbf{LS2}} &
   		\textbf{LSTM} & 0.46 & 0.55 & 0.37 & 0.55 & 0.36 & 0.61 & 0.32 \\
   		& \textbf{AE} & 0.58 & 0.66 & 0.41 & 0.61 & 0.57 & 0.75 & 0.45 \\
   		& \textbf{BiGAN} & 0.54 & 0.72 & 0.33 & 0.65 & 0.56 & 0.68 & 0.31 \\ \hline \hline
   		\multirow{3}{*}{\textbf{LS3}} &
   		\textbf{LSTM} & 0.43 & 0.58 & 0.34 & 0.44 & 0.33 & 0.58 & 0.30 \\
   		& \textbf{AE} & 0.57 & 0.72 & 0.40 & 0.64 & 0.53 & 0.69 & 0.44 \\
   		& \textbf{BiGAN} & 0.52 & 0.83 & 0.31 & 0.33 & 0.51 & 0.72 & 0.43 \\ \hline \hline
   		\multirow{3}{*}{\textbf{LS4}} &
   		\textbf{LSTM} & 0.47 & 0.57 & 0.37 & 0.56 & 0.38 & 0.60 & 0.35 \\
   		& \textbf{AE} & 0.57 & 0.65 & 0.40 & 0.63 & 0.55 & 0.79 & 0.43 \\
   		& \textbf{BiGAN} & 0.52 & 0.81 & 0.36 & 0.25 & 0.54 & 0.69 & 0.48 \\ \hline
    \end{tabular}
    \caption{\small App-level separation results (LS1:4, FS$_{\text{custom}}$, AD2)}
    \label{tab:custom-ls1:4-ad2-app-separation}
\end{table}

\section{Conclusions and Future Directions} \label{sec:conclusions}

In this paper, we presented Exathlon\ -- a novel public benchmark for explainable AD, and demonstrated its utility through an experimental analysis of selected AD and ED algorithms from recent literature. Our AD results show that Exathlon's dataset is valuable for evaluating AD algorithms due to rich signals and diverse anomaly types included in the data. Yet more importantly, our results reveal the limitations of these AD methods for semi-supervised learning under noisy training data and
mixed anomaly types.
On the ED front, the literature lacked comparative analysis tools and studies. Our benchmark fills this gap by providing a common framework for analyzing the strengths and limitations of diverse ED methods in their conciseness, consistency, accuracy, and efficiency.
These results call for new research to advance the current state of the art of AD and ED, as well as integrated solutions to anomaly and explanation discovery. For a true integration, ED methods should first become capable of discovering range-based explanations, which is also a key step towards automated root cause analysis (a.k.a., ``why explanations''). Exathlon's dataset and extensible design are well-positioned to support research progress towards these goals in the long term. Going forward, we envision Exathlon\ to develop into a collaborative community platform for fostering reproducible research and experimentation in the area. We intend to actively maintain and extend this platform, as well as welcoming feedback and contributions from the AD and ED communities.

\section{#1}\vspace{-0.5mm}}
\newcommand{\ssubsection}[1]{\vspace{-2mm}\subsection{#1}\vspace{-0.5mm}}
\newcommand{Exathlon}{Exathlon}
\newcommand{\msbold}[1]{\boldsymbol{~#1}}

\begin{document}
\title{Exathlon: A Benchmark for Explainable Anomaly Detection \\ over Time Series}
\author{
Vincent Jacob$^\dagger$
\hbox{~~~~}Fei Song$^\dagger$
\hbox{~~~~}Arnaud Stiegler$^\dagger$
\hbox{~~~~}Bijan Rad$^\dagger$
\hbox{~~~~}Yanlei Diao$^\dagger$
\hbox{~~~~}Nesime Tatbul$^\ast$
}
\affiliation{\institution{
$^\dagger$Ecole Polytechnique, France
\hbox{~~~~}$^\ast$Intel Labs and MIT, USA}
}
\email{
{vincent.jacob, fei.song, arnaud.stiegler, bijan.rad, yanlei.diao}@polytechnique.edu, tatbul@csail.mit.edu
}

\begin{abstract}
Access to high-quality data repositories and benchmarks have been instrumental in advancing the state of the art in many experimental research domains.
While advanced analytics tasks over time series data have been gaining lots of attention, lack of such community resources severely limits scientific progress.
In this paper, we present Exathlon, the first comprehensive public benchmark for explainable anomaly detection over high-dimensional time series data. Exathlon\ has been systematically constructed based on real data traces from repeated executions of large-scale stream processing jobs on an Apache Spark cluster. Some of these executions were intentionally disturbed by introducing instances of six different types of anomalous events (e.g., misbehaving inputs, resource contention, process failures). For each of the anomaly instances, ground truth labels for the root cause interval as well as those for the extended effect interval are provided, supporting the development and evaluation of a wide range of anomaly detection (AD) and explanation discovery (ED) tasks.
We demonstrate the practical utility of Exathlon's dataset, evaluation methodology, and end-to-end data science pipeline design through an experimental study with three state-of-the-art AD and ED techniques.
\end{abstract}

\maketitle

\pagestyle{\vldbpagestyle}

\input{introduction}
\input{related-work}
\input{data-collection}
\input{benchmark}
\input{pipeline}
\input{experimental-setup}
\input{experiments}
\input{ed-experiments}
\input{conclusions}

\balance
\bibliographystyle{ACM-Reference-Format}


\section{Introduction}
An intriguing feature of colloids is their ability to self-assemble into ordered structures with interparticle distances commensurate with wavelengths of light \cite{Manoharan2015}.  Advances in  synthesis have produced a wide variety of anisootropic particles, such as polyhedra \cite{Henzie2012,Young2013, Gong2017}, dumbbells \cite{Forster2011}, spherocylinders \cite{Hosein2010}, superballs \cite{Meijer2017} and octapods \cite{Miszta2011}. Experiments and simulations have demonstrated a diverse range of close-packed superlattices whose structure depends on particle shape. Simulations of hard colloids (e.g., Refs.~\cite{Gang2011, Haji-Akbari2009, Kraft2012, Agarwal2011, Damasceno2012_science, Ni2012, Marechal2010, Chen2014, Wan2019}) predict complex crystals from an even larger variety of anisotropic shapes, which are versatile in terms of modification and functionalization.  Colloidal self-assembly is one route that has been explored to fabricate photonic crystals with photonic band gaps (PBGs) (e.g., Refs.~\cite{Ye2001,Vlasov2001, Hosein2010, Forster2011,Sowade2016}). Different from top-down design, this bottom-up method has advantages such as low cost and low energy consumption, and crystals can be produced over large areas \cite{Moon2010, Kim2011, Zhao2014}. Moreover, the particle size can be varied from tens of nanometers to micrometers in experiments, tuning the periodicity of crystals and consequently the PBG frequency \cite{Moon2010}. Besides the photonic crystals already realized in experiments, many other promising lattice structures from self-assembly approaches have been theoretically proposed (e.g., Refs.~\cite{Ding2014, Hynninen2007, Woldering2011, Pattabhiraman2017, Wang2017, Changizrezaei2017, Cersonsky2018, Lei2018, Cersonsky2019}). Previous studies of self-assembled photonic crystals have focused primarily on structure design and self-assembly pathways. How ``randomness" arising from the inevitable thermal noise in a self-assembled colloidal crystal affects PBGs has received little attention. At first glance one might expect noise to weaken or destroy photonic band gaps, but is this generally true?  If it is, it would suggest that thermodynamic self-assembly may not be a viable synthesis route for PBG materials.

Here we study the PBGs of self-assembled two-dimensional photonic colloidal crystals (Fig.~\ref{sa} shows some examples) using computer simulations. We consider self-assembled lattices of parallel dielectric rods of infinitely long length that interact through a hard core potential. We study rods of circular and square cross section and investigate a wide range of packing densities $\phi$ and dielectric constants
$\epsilon$. Surprisingly, we find that the widest PBG does not always appear at the packing fraction where the rods are perfectly ordered. Moreover, for rods with square cross section at intermediate packing densities, we find that the transverse magnetic (TM) band gap of the self-assembled system can be larger than that of its corresponding perfect system. Further, we show that by considering hollow rods and optimizing the internal radius, the PBG of the transverse electric (TE) mode can be substantially increased while that of the TM mode does not show obvious improvement. We discuss the possibility of engineering PBGs in self-assembled colloidal systems by controlling relevant factors.

\section{Methods}
To be consistent with previous work, we adapt the method used in Ref.~\cite{Michael2016} to generate self-assembled structures of cross sections of $N=200$ circular hard rods. For particles of other shapes, we performed Monte Carlo (MC) simulations with periodic boundary conditions using the hard particle Monte Carlo (HPMC) module in HOOMD-blue \cite{hoomd,Glaser2015}. Simulations were initialized at very low packing density ($\phi=0.01$) in a random configuration and slowly compressed to a target packing density using MC simulations. Here one MC step consists of $N+1$ trial moves including translation (plus rotation for rods of square cross section) of $N$ particles or rescaling of the box, where during each compression step the length of the sides of the box are rescaled to 0.9995 of their current value. At high packing densities, box rescaling can create  unphysical overlaps that allow us to reach higher densities; these overlaps are subsequently eliminated with isochoric Monte Carlo \cite{Haji-Akbari2009}. After the system reaches the targeted packing density, it is further equilibrated for $10^6$  MC steps. We use the supercell method \cite{Joannopoulos2008} implemented in the open source code MIT Photonic-Bands \cite{Johnson2001} to obtain the photonic band structure of equilibrated snapshots. All band gap sizes of self-assembled structures presented in this work are averaged over five independent runs.

\section{Results and discussion}

\begin{figure}
\centering
\includegraphics[width=3.5in]{fig1}
\caption{(Color online) Snapshots of cross section of a system of $N=200$ self-assembled rods described in the main text. (a) Rods of square cross section at packing density $\phi=1$ (homogeneous material, 196 rods shown), 0.85 and 0.4 (left to right).  (b) At $\phi=0.85$, equilibrium configurations of rods of circular, hexagonal and triangular cross sections (left to right).}
\label{sa}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=3.5in]{fig2}
\caption{(Color online) Photonic band structure and band gap for arrays of parallel dielectric rods of circular cross section with radius $r=0.5a$ ($a$ is the diameter and the length unit) embedded in air ($\epsilon=1$). (a-c) Rods in a perfect triangular lattice. (a) An example of the band structure with TE polarization, at $\phi=0.5$ and $\epsilon=20$. The inset shows the first Brillouin zone for the periodic structure studied, with the symmetry points indicated. The yellow area represents the band gap observed. (b) Gap size (the percentage of $\triangle \omega / \omega_{0}$) as a function of $\phi$ and $\epsilon$ for the TE mode. (c) Gap size for the TM mode. (d-f) Parallel results for snapshots of $N=200$ self-assembled rods. Parameters in (d) are the same as those in (a). The band gap (e-f) is obtained by averaging over five independent simulation snapshots. }
\label{disks}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=3.5in]{fig3}
\caption{(Color online) Gap size as a function of $\phi$ and $\epsilon$ for rods with square cross section in a perfect square lattice (a-b) and the self-assembled crystal (c-d). The white dashed line in (d) indicates the region where both the perfect square lattice and the self-assembled crystal have a non-vanishing TM gap, and the gap size of the self-assembled structure is wider than in the perfect structure. }
\label{squares}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=3.3in]{fig4}
\caption{(Color online)  Electric field distribution in the system of $N=500$ self-assembled rods with square cross section for TM polarization: $\phi=0.5$ and $\epsilon=20$. (a) Extended mode before the PBG. (b) Localized mode before the PBG. (c) Localized mode after the PBG. (d) Extended mode after the PBG.}
\label{modes}
\end{figure}

\begin{figure*}
\centering
\includegraphics[width=7in]{fig5}
\caption{(Color online) Band gap for an array of parallel dielectric rods of annular cross section (a-d) or hollow square cross section (e-h). (a) Annular cross section. (b) Band structure at optimal $r_{in}$ value (described in the text) for TE polarization. Color indicates the value of $r_{in}/r_{out}$. (c-d) Band gap at optimal $r_{in}$, as a function of $\phi$ and $\epsilon$, for the self-assembled structure. (e) Hollow square cross section. (f) The optimal $a_{in}$ value for TE polarization. Color indicates the value of $a_{in}/a_{out}$. (g-h) Band gap at optimal $a_{in}$, as a function of $\phi$ and $\epsilon$, for the self-assembled structure.}
\label{hollow}
\end{figure*}

We start from the photonic band structure for electromagnetic waves in a periodic array of parallel dielectric rods of circular cross section, whose intersection with a perpendicular plane form a perfect triangular lattice. This system was first studied in Ref.~\cite{Plihal1991}. The band structure can be calculated from a unit cell and Fig.~\ref{disks}(a) shows an example for TE polarization at $\phi=0.5$ and dielectric constant $\epsilon=20$. The yellow area indicates there is a band gap between the first and second bands. The relative gap size is defined as $\triangle \omega / \omega_{0}$, where $\triangle \omega$ is the width of the PBG and $\omega_{0}$ is the central frequency. Fig.~\ref{disks}(b and c) show the relative gap size as a function of $\phi$ and $\epsilon$. We investigated a wide range of $\phi$  and with $\epsilon$ in the range from 2 to 20.  Some low $\epsilon$  materials such as polystyrene and silica, and some high $\epsilon$ materials such as titania, selenium and amorphous silicon, fall in this range \cite{Moon2010}. For the TE mode, the largest PBG is at $\phi=0.6$ and $\epsilon=20$; for the TM mode, the largest PBG is at much lower packing density, i.e., $\phi=0.1$ and $\epsilon=20$. Furthermore, the largest band gap for the TM mode is much wider than that of the TE mode.
Fig.~\ref{disks}(d-f) show the equivalent results calculated for simulation snapshots of a self-assembled system of 200 rods (Fig.~\ref{sa}(b) (left)).
Fig.~\ref{disks}(d) shows the band structure with the same parameters as those in Fig.~\ref{disks}(a). The yellow area indicates the band gap between the $N$th and the $(N+1)$th band. The band structure shows that the central frequency remains nearly constant, but the gap size decreases. Fig.~\ref{disks}(e,f) show the relative PBG over the entire range of parameters.  Because the bands are dense near the band edge and the deviation is relatively small, a sample size of five independent configurations is sufficient to obtain good statistics (see the Supplimentary Material for more information). We see that for the TE mode, the maximal band gap of the self-assembled system (Fig.~\ref{disks}(e)) appears at about $\phi=0.6$, similar to that in the perfect triangular system (Fig.~\ref{disks}(b)). At $\phi=0.9$, where the system is most ordered (close to the densest packing fraction 0.907), there is no band gap. We observe that for the TM mode, the packing fraction at which the maximal band gap appears shifts from $\phi=0.1$ (Fig.~\ref{disks}(c)) to $\phi=0.3$ (Fig.~\ref{disks}(f)), and there is no obvious band gap at $\phi=0.1$ in the self-assembled system. This is understandable because at very low density the system is highly disordered, where we do not expect a PBG. Thus, surprisingly, the largest gap does not occur when  the system is the most ordered. However, at the same $\phi$ and $\epsilon$ values, the self-assembled system always has a smaller PBG compared to that of the perfect ordered system. This is consistent with Ref.~\cite{Michael2016} that used a hard disk system to generate a seed pattern and placed cylindrical rods with an arbitrary fixed radius at the points of the seed pattern. Here we consider the seed pattern as the self-assembled structure of infinitly long rods. In practice, when rods are long enough and at high packing densities, they will tend to align parallel to each other \cite{Anders2014PNAS, Anders2014, Manoharan2015, Wan2018}.  Experiments have been able to obtain monodomain films of highly aligned carbon nanotubes from suspension [43]. We expect that similarly, by taking advantage of entropic effects and using auxiliary experimental skills, it is possible to align long colloidal rods as well.

Particle shape affects the PBG of a self-assembled crystal in two ways: (1) it defines the region of dielectric materials and (2) it determines the assembled structure and densest packing structure \cite{Bernard2011, Engel2013, Anderson2017}. Some examples are shown in Fig.~\ref{sa}. As a simple example, although infinite rods of circular or hexagonal cross section crystallize into triangular lattices, they have different maximum packing densities. Rods with square and triangular cross sections crystallize into square and hexagonal lattices, respectively, already below their maximum packing densities of 1.

To investigate the effect of cross-section shape on PGBs, we studied rods of square cross section (Fig.~\ref{squares}). The photonic band structure of rods with square cross section are similar to those for rods with circular cross section, but the gap widths, are generally different at the same packing density and with the same dielectric constant, which suggests the possibility of engineering PBGs using particle shape \cite{Cersonsky2018}. Moreover, the existence of PBGs in the square rod system suggests that Voronoi particles, which also tile space at densest packing, can be candidates for three-dimensional photonic crystals \cite{Schultz2015}. Another interesting aspect of the square rod system is that at some intermediate packing densities, the TM band gap of the self-assembled system can be larger than that of rods in a perfect square lattice at the same value of $\epsilon$. The white dashed line in Fig.~\ref{squares}(d) indicates the region ($\phi=0.5$, $\epsilon \in [6,20]$) where both the self-assembled structure and the perfect square lattice have a non-vanishing PBG, and the gap size of the former is wider. For example, at $\phi=0.5$ and $\epsilon=20$, the gap size of the self-assembled system is about $10.72\%$ and that of the perfect system is about $6.94\%$.  To further demonstrate the TM band gap, we show the electric field distribution of some TM modes around the PBG edges in Fig.~\ref{modes}. Localized and extended modes around a PBG gap have also been observed in other 2D systems (e.g.~\cite{Florescu2009, Ricouvier2019}).
At $\phi=0.6$ and $\epsilon=20$, the gap size of the self-assembled system decreases to about $5.96\%$ (Fig.~\ref{squares}(d)) while that of the perfect system vanishes (at $\phi=0.6$, the perfect square lattice has no TM band gap for all $\epsilon$ values) (Fig.~\ref{squares}(b)).  Randomness due to fabrication errors in traditional lithography-based approaches is usually regarded as a bad aspect which deminishes PBGs, here instead we find a counterintuitive example where randomness actually helps to increase the PBG.

The PBGs for rods with circular or square cross section, however, are not large. In an attempt to increase the PBGs, we consider the design of the internal structure of the rods. Hollow or double-layer rods of various shapes and on various perfect lattice structures have been explored in previous studies \cite{Chang2006, Duque2012, Niu2013, Liu2013, Liu2015, Rezaei2006, Xiao2008}. Here we consider the simplest situation of hollow rods as indicated in Fig.~\ref{hollow}(a). We optimize the value of the inside radius $r_{in}$ to that which maximizes the TE(TM) band gap when the rods are in a perfect triangular lattice; the optimal $r_{in}$ for the TE mode is shown in Fig.~\ref{hollow}(b). Fig.~\ref{hollow}(c,d) show the band gap of the self-assembled structure at the optimal radius for the TE and TM modes, respectively. Compared with solid rods Fig.~\ref{disks}(e,f), the band gap of the TE mode is largely increased while that of the TM mode remains nearly the same. The value of $\phi$ at which the maximal PBG of the TE mode occurs changes from $\phi=0.6$ to $\phi=0.9$. Note that the band gap of the self-assembled structure with the optimal $r_{in}$ is smaller than that of the perfect lattice with the optimal $r_{in}$ (see the supplementary material). For rods with square cross section (Fig.~\ref{hollow}(e)), making the rods hollow also increases the PBGs of the TE mode while hardly changing the TM mode (Fig.~\ref{hollow}(f-h)). The optimal value of the inside side length $a_{in}$ maximizes the band gap when the rods are in a perfect square lattice.

\section{Conclusions}
We studied the PBGs of a self-assembled system consisting of parallel dielectric colloidal rods with circular or square cross section and interacting through a hard core potential. For the square rod system, we found that although the system has no PBGs when it is highly/perfectly ordered, i.e., close to or at the densest packing fraction, there is a wide range of intermediate packing fractions where a PBG exists in the self-assembled (noisy) structures.
Moreover, at some intermediate packing densities, the ``randomness" of the self-assembled system improves the TM band gap compared to that of rods in the corresponding perfect lattice (provided a PBG exists). The width of the PBG is packing-density dependent, suggesting that in experiments the PBG can be adjusted through, e.g., changing the concentration of colloidal particles in suspension by adding or decreasing solvent. A comparison of rods with circular and square cross section suggests that shape can be used as another control factor to engineer the PBG of a self-assembled system. Further, we showed that by making rods hollow and optimizing the internal radius, the PBG of the TE mode can be significantly increased. Internal rod structure is but one ``dimension" (examples of other dimensions are illustrated in Ref.~\cite{Glotzer2007,Gang2011}) available to engineer particles that can produce a PBG. In all, our study suggests that by suitably choosing the packing density, particle shape, and engineering other dimensions such as particle internal structure, self-assembly can indeed be a promising method to make photonic crystals with large band gaps despite the inherent thermal noise (randomness).

\acknowledgments
This work was supported by a grant from the Simons Foundation (256297, SCG) and used the Extreme Science and Engineering Discovery Environment (XSEDE) \cite{Towns2014}, which is supported by National Science Foundation grant number ACI-1548562; XSEDE award DMR 140129. Computational resources and services were also supported by Advanced Research Computing at the University of Michigan, Ann Arbor. D.W. also thanks the National Natural Science Foundation of China (Grant No. 11904265) for support.