\section{Appendix}

\input{17_appendix_formulation}
\input{19_appendix_results}

\end{document}
\endinput

\section{Related Work}
\label{Section:RelatedWork}

Existing detection techniques span a wide variety of approaches. One approach converts last couple layers of a neural network into Bayesian layers~\cite{deng2021libre} while another does neural fingerprinting~\cite{dathathri2018detecting}. More common techniques include training a second-stage classifier~\cite{grosse2017statistical,gong2017adversarial,metzen2017detecting,aigrain2019detecting}, detecting statistical properties~\cite{bhagoji2017dimensionality,hendrycks2016early,li2017adversarial,crecchi2019detecting, pang2017towards}, and running statistical tests~\cite{grosse2017statistical,feinman2017detecting,roth2019odds}. We aim to build statistical AFVs and run both detection of attacks and classification of the type of attack.

To simplify the detection decision, detectors rely on generating an alternative
(often high dimensional) embedding space over which to project natural and adversarial examples
into manifolds where to then discern boundaries or nearest neighbors.
By comparing histograms from adversarial examples against that of natural examples, we compose
anomaly detection features to span such manifolds based on magnified discrepancies between those histograms.  Using an anomaly-oriented statistical profile to produce an embedded space,
we also leverage KNN distances~\cite{cohen2020detecting, papernot2018deep} to both neighboring natural and adversarial examples to train the detector neural network.
However, although these approaches claim high accuracy for several attack methods,
they require an iterative DkNN training phase which is computational expensive. Similarly, unlike the kernel density estimation in ~\cite{pang2018robust}, our approach approximates KDE features using inexpensive histogram comparison analysis while again requires no training on the pretrained model being defended -- only evaluation mode access.

 In ~\cite{metzen2017detecting},

 the authors also train a binary detector.  The detector is fed inputs from intermediate feature representations of a classifier, to discriminate between samples from the original data set and adversarial examples.

Their results on CIFAR10 assert that their detector achieves $91.3\%$ accuracy in detection of adversarial inputs.
Moreover, while the authors train the detector to detect perturbations of a specific adversary, the authors assert that the detector generalizes to similar and weaker adversaries.

Existing detection techniques span a wide variety of approaches. One approach converts last couple layers of a neural network into Bayesian layers~\cite{deng2021libre} while another does neural fingerprinting~\cite{dathathri2018detecting}. More common techniques include training a second-stage classifier~\cite{grosse2017statistical,gong2017adversarial,metzen2017detecting,aigrain2019detecting}, detecting statistical properties~\cite{bhagoji2017dimensionality,hendrycks2016early,li2017adversarial,crecchi2019detecting}, and running statistical tests~\cite{grosse2017statistical,feinman2017detecting,roth2019odds}. We aim to build statistical AFVs and run both detection of attacks and classification of the type of attack.

\js{
The authors in ~\cite{bhagoji2017dimensionality} improve the robustness of ML classification and training phases against evasion attacks by reducing data dimensionality. Although a few different distinct types of attacks have been selected from the literature in their threat model, the paper only focuses on the defense against attacks across multiple ML classifier, whereas DeClaw focuses on how to differentiate and detect different groups of attacks.

~\cite{dathathri2018detecting} encodes fingerprint patterns into the response of a DNN around real data, which characterize the DNN’s expected behavior around real data and can later be used to detect adversarial examples (similar to the idea of Dynamic Watermarking). However, fingerprint-based approach does require encoding certain features into the data distribution in a proactive way, and it does not differentiate different class of attacks like DeClaw.

~\cite{gong2017adversarial} also looks into disparity among adversarial samples. However, this work only studies 3 attacks (FGSM, TGSM and JSMA) and authors also admitted that "FGSM and JSMA generate adversarial datasets that are far away from each other. But it is too vague without defining precisely what is being far away". In contrast, DeClaw studies 11 different types of attacks. DeClaw also compares and differentiates them using a more quantitative distribution approach, which potentially leads to a more accurate and effective way for anomaly detection. One detection method proposed in

~\cite{li2017adversarial} is based the observation that correctly classified and many out-of-distribution examples tend to have different softmax distributions. However, this paper did not focus on the attack classification, nor study how exactly the distribution are different caused by different type of adversarial attacks.

In~\cite{aigrain2019detecting}, the authors augmented classification by branching off the main network at certain layer, and then this subnetwork detector is trained on the adversarial input and produce an output which is interpreted as the probability of the input being adversarial. However, the paper only compares the accuracy of the detector against 4 different adversaries with the focus more on their detector, rather than categorizing different type of attacks based on their distributional characterization.

The authors in~\cite{cohen2020detecting} distinguishes adversarial examples on six attack methods (FGSM, JSMA, Deepfool, CW, PGD,and EAD) based on the observation that the impact of training samples on the validation set data are different. By fitting a k-NN model on the DNN's activation layers, the
authors shows that the samples are highly correlated with the nearest neighbors of the normal inputs, while this correlation is much weaker for adversarial inputs. However, the work does not address how to identify different types of adversarial attacks based on these k-NN features.

The proposed method in~\cite{yang2019mlloo} is based on the observation that the feature attribution map of an adversarial example near the boundary always differs from that of the corresponding clean example. The authors also use statistical dispersion in feature attribution to detect 6 different adversarial examples. However, compared to DeClaw, this approach is focused on the confidence level rather than distributional characterization of different attacks. Therefore, it does not help classify and differentiate different attacks.
}

\section{Old Abstract}

\section{Introduction}
\label{Section:Introduction}

While deep neural networks (DNNs) help provide image classification, object detection, and speech recognition in autonomous driving, medical, and other domains, they remain vulnerable to adversarial examples.  In this paper, we focus on the detection problem -- determining whether an input is adversarial and even attempting to determine its type (e.g., PGD vs. Carlini-Wagner), so that potentially in the future, attack-specific countermeasures can be taken.

Existing detection techniques include training a second-stage classifier~\cite{grosse2017statistical,gong2017adversarial,metzen2017detecting,aigrain2019detecting}, detecting statistical properties~\cite{bhagoji2017dimensionality,hendrycks2016early,li2017adversarial,crecchi2019detecting,pang2017towards}, and running statistical tests~\cite{grosse2017statistical,feinman2017detecting,roth2019odds}. In contrast to our work, these approaches are focused on detection of an attack rather than determining the attack type. It is an open question as to whether attacks can be distinguished sufficiently. If the answer is yes, that may help provide more insight into how attacks differ, leading to more robust classifiers, or attack-specific mitigation strategies.

To explore the question, our hypothesis is that it is possible to do a fine-grained distributional characterization of latent variables at an intermediate hidden layer across natural input samples. Our hope is that adversarial samples are outliers with respect to the that characterization in different ways. Unfortunately, finding differences between natural and different types of adversarial inputs is not straightforward.  Fig.~\ref{distribution} shows histograms of (standarized) layer input values for a NATURAL/CLEAN sample as well as for two (first-stage) white-box attacks (Fast Gradient Sign (FGM) and $L_\infty$ Projected Gradient Descent (PGD) with 20 steps and $\epsilon = 8 / 255$. over a CIFAR10 pretrained model.
Each histogram shows frequency counts for z-scores of observed values.
While these distributions exhibit some differences, the differences are minute, making discrimination among them challenging.

We propose an approach called DeCLaW (Detecting, Classifying, and Warning of adversarial examples) that aims to overcome this challenge by constructing a set of anomaly  features, in the form of an {\em anomaly feature vector} (AFV), from an input image that better distinguishes among clean and adversarial images as well as several of the common attack types.
The key insight is to collect and look at statistical values in outlier areas in Figure~\ref{distribution} (shown with a magnifying glass) and try to characterize their distributions for different types of attacks, generating a collection of anomaly features that form the input image's AFV. We find that, given a model to be protected, each anomaly feature  in the vector induces some separation between different types of attacks from clean examples. While overlaps exist when looking at one anomaly feature, across multiple anomaly features, we find good success in distinguishing among attack types.

DeClaW uses a second stage neural network that is trained on AFV values rather than latent features. Given an image, its AFV is computed from the hidden layer values in the target classifier and then fed to the second stage neural network for analysis of the attack type. As a detector,  using 11 different popular attack methods from the IBM Adversarial Robustness Toolbox (ART), DeClaW\xspace achieves 96\% accuracy for CIFAR10 on adversarial examples and 93\% accuracy on clean examples.

Notably, AFV size is small compared to input pixel values or the latent variables, resulting in a simpler second stage classifier.  E.g., the current implementation of DeClaW  uses only 176 anomaly detection features in contrast to approximately 16K latent features used by Metzen et al.’s adversarial example detector [18].

As an attack-type classifier,  after grouping similar attacks into clusters, DeClaW\xspace achieves an overall $93.5\%$ classification accuracy for CIFAR10. To our knowledge, we are the first to produce good attack classification. For example, whereas the work in~\cite{metzen2017detecting} reports detection rate on three attack methods (i.e., FGM, DeepFool, BIM), it does not discriminate between those underlying attacks. DeClaW is able to discriminate among these attacks with a high success rate.

\begin{figure}[t]
\centerline{
\includegraphics[width=0.4\textwidth]{FIGURES/distribution_maskarading3.png}}
\caption{
Histograms of (standarized) layer input values for one NATURAL sample and a adversarial versions from FGM and PGD20 of the same CIFAR10 sample.  Attack perturbations tend to be small, but do induce potential differences when carefully analyzed. }
\label{distribution}
\vspace{-0.2in}
\end{figure}

\section{Motivating Goals}
\label{Section:MotivatingGoals}

Within the domain of image classification, adversarial examples are generally modeled as noise perturbations over natural examples (e.g., seemingly apparent random noise over an image) that force a classification error. This is referred to as an evasion attack. Another more subtle form of an adversarial example is a poisoning attack, where the attacker is able to poison the natural training samples with examples to distort the decision boundaries of the resulting model.
While developing robust models to handle adversarial inputs is difficult, new attack methods keep emerging. Therefore, a detection technique becomes a useful defense option. In particular, it would be desirable to obtain a defense that could identify:
\begin{itemize}
\item adversarial examples during classification (e.g., evasion)
\item adversarial examples during training (e.g., poisoning)
\item out of distribution, low confidence examples during classification (e.g., samples with very low or null trained classification support)
\item out of distribution, high distortion examples during training (e.g., samples with very high anomaly scores).
\end{itemize}

\subsection{Grand Vision}
\label{Section:GrandVision}

Fig.~\ref{vision_augmented_models} depicts our vision for augmenting classifier models with robust introspection features. In particular, these intelligent classifier qualification (shadow/second-stage) pipelines are foreseeable to be dominant in the automotive domain. Early detection and classification of adversarial inputs during early stages in a pipeline would allow the machine learning to reflect on the sample and speculate:
\begin{itemize}
    \item is the sample out-of-distribution? e.g., a sample may represent a statistical outlier with respect to a normative set of natural samples.
    \item is the sample anomalous? e.g., a sample could exhibit patterns of statistical departures from normative natural samples --- yet without being a statistical outlier.
    \item is the sample adversarial? e.g., a sample may align to existing statistical profiles, clusters, etc for known (or representative) attack samples.
    \item is there confidence that the sample aligns to concepts previously assimilated? e.g., a sample may align to an existing concept from the training data for which a number of similar samples reassert a support for the classification decision.
    \item is the collection of samples indicative of a process shift from the concepts previously assimilated? e.g., a stream of samples that drift away from concepts learned from the training data into a different concept may now raise a warning signal
\end{itemize}

A system (e.g., Declaw) providing answers to these questions would allow the defended machine learning intelligence, to adapt, retrain, warn and/or redirect samples to appropriate remedial measures such as human intervention, alternative models/weights, leverage of ensembles, etc.

\section{Our Approach}
\label{Section:OurApproach}

DeClaW\xspace is both an attack labeling classifier as well as a binary attack detector. The fundamental idea is to augment a pretrained model   with a sampling hook at a chosen layer to capture unperturbed data that can be analyzed and statistically characterized. That is then used to generate what we term as an anomaly feature vector when  presented a perturbed input of a given type (e.g., adversarial attack generated by PGD). DeClaW\xspace's AFV generation sampling hook as shown in Fig.~\ref{samplinghook}.  The AFVs from the perturbed inputs as well as natural inputs of different classes are used to train a second stage neural network (Fig.~\ref{declaw_network_model}) that can help distinguish among perturbed inputs and clean inputs and also help identify the type of perturbed input (e.g., a PGD attack versus a DeepFool attack). As we discuss later, some attack types may be difficult to distinguish reliably. When that is the case, those attack types are clustered together using a clustering algorithm and treated as one composite class.  Below, we describe the key concepts of the scheme.

\begin{figure}
\centerline{\includegraphics[width=0.4\textwidth]{FIGURES/DeClaWSamplingHook.png}}
\caption{DeClaW\xspace subclassing extends sampling hook into the target pretrained model at a batch normalization layer that collects latent feature statistics on the dataset. Another pass generates AFVs for natural and adversarial samples.
}
\label{samplinghook}
\end{figure}

\begin{figure}[t]
\centerline{\includegraphics[width=0.45\textwidth]{FIGURES/declaw_network_model.png}}
\caption{DeClaW\xspace's second stage classifier takes an image's AFV as input and outputs attack type (or clean). Dropout layers are used after each hidden layer for regularization.}

\label{declaw_network_model}
\end{figure}

One advantage of our approach is that the the size of AFV is designed to be independent of the latent feature vector size and is typically orders of magnitude smaller than the typical size of the feature vector from which it is derived. Thus, training and classification of an input using the second stage classifier of DeClaW\xspace is efficient. We elaborate on this in Section~\ref{Section:Experiments} (Experiments).

\subsection{Anomaly Feature Vectors}
An anomaly detection defense can be conceptualized as testing whether a sample is anomalous or not by checking whether statistics for the given sample are within the expected limits of the non-anomalous (natural) inputs.
However, because the goal of sophisticated attacks is to resemble as much as possible a natural input, a number of coarse-grain $A/B$ statistical tests simply lack the resolution to reliably discern anomalies.
We rely on various histogram comparative techniques, subsampling projections, aggregation predicates and distance metrics to account for shape and area differences between the histogram of expected layer inputs and the histogram for observed layer inputs associated with an evaluated input image.

DeClaW\xspace's neural network trains on a small and compact number of anomaly-detection distribution characterization features as opposed to on pixels or their derivative latent variables.
The basic DeClaW\xspace second stage neural network uses $176$ anomaly detection features -- a significantly smaller number than the number of features used by Metzen et al.'s adversarial example detector~\cite{metzen2017detecting}. In Metzen et al.'s detector, the preferred position reported is at stage $AD(2)$, the second residual block, where the inputs have a shape of $32 \times 32 \times 16 \approx 16K$ which is then fed to the detector convolution network.
Moreover, the number of AFV features is independent of the input size since AFV features are essentially distributional characterization features.
DeClaW\xspace trains fast, converging in few epochs.

\begin{figure*}[t]
\begin{subfigure}{0.32\textwidth}
\centering
\includegraphics[width=\linewidth]{FIGURES/LEFT_DISTRIBUTION.pdf}

    \label{fig:hierarchicalscene}
\end{subfigure}
\begin{subfigure}{0.32\textwidth}
    \centering
    \includegraphics[width=\linewidth]{FIGURES/ZSUM_distribution.png}

    \label{fig:lanedetection}
    \end{subfigure}

    \begin{subfigure}{0.32\textwidth}
    \centering
    \includegraphics[width=\linewidth]{FIGURES/OVERALL_DISTRIBUTION.png}

    \label{fig:sequential}

\end{subfigure}
\vspace{-0.2in}
\caption{Each anomaly feature helps provide partial discrimination among attacks and clean examples.  When used in combination, they can help towards classifying attack types and whether an example is clean.}
\label{fig:anomaly_feature_regions}
\end{figure*}

Fig.~\ref{fig:anomaly_feature_regions} depicts three histogram plots corresponding to different DeClaW\xspace anomaly features evaluated across the same set of $5000$ samples.
On each plot, the blue shade represents the distribution of feature values for natural samples
--- each other color represents the distribution of feature values for the application of a given attack method to the same $5000$ samples. Each anomaly feature is generated by applying a subsampling projection that zooms in a particular range of standardized scores (z-scores) and then by collecting aggregation metrics over the selected values (e.g., frequency counts of the selected z-scores). Each plot illustrates the discriminating power that such region-based aggregation feature induce. For example, the first plot shows a histogram for the values taken by an AFV feature that counts how many standardized scores are observed with value less than a minimum tolerance (i.e., the LEFT tail of the histogram).

We rely on four subsampling projections or regions: the LEFT tail (that is, low z-scores), the CENTER (that is, non-anomalous z-scores), the RIGHT tail (that is, high z-scores), and the OVERALL histogram -- and ratios of these.
Different aggregation predicates (e.g., count, sum, mean, ratio, etc) are used to aggregate over the subsampling projection.

A total of 12 different shades are shown corresponding to 11 attack methods that our second stage detector seeks to detect -- plus the natural/clean samples shown in a blue shade.

Intuitively, whereas the detection goal is to separate the blue shade from other shades, the classification goal is to find combinations of features that allow to separate each shade from the others.
Some features (e.g., Z.SUM()) induce linear separatability between natural samples and attack classes
--- hence, DeClaW\xspace's high detection accuracy and low FPR as well as FNR
--- whereas other features induce partial separatability between attack classes (hence DeClaW\xspace's high classification accuracy).
For example, the above-mentioned Z.SUM() feature measures the sum of absolute standardized scores of layer-input values observed for a given sample and empirically discriminates well among AFVs from CLEAN/NATURAL (i.e. blue shade) from AFVs from any adversarially perturbed samples (i.e., all others).

\section{Attack Classes, Datasets and Parameters}

\subsubsection{Distributional Characterization}\label{CHAR}
In~\cite{metzen2017detecting}, a binary classifier for detection attacks was built by feeding the latent variables at selected intermediate stage, to a separate CNN.
This approach was able to detect with $91\%$ accuracy, the presence of attacks on three attack classes:
PGD, BIM, DeepFool.
Rather than feeding the latent variables to the secondary neural net, we propose feeding a more compact intermediary representation centered around anomaly detection. The specific features of this intermediary representation are discussed next. Our approach achieves a detection rate of $95.3\%$ and a classification rate of $93\%$ with FPR and FNR rates below $5\%$ --- while classifying not three but rather eleven different attacks (three PGD variants, FGM, BIM. CWL2, DeepFool, Spatial Transform, Square Attack, Wasserstein, and Newtonfool).

Distributional A/B tests (e.g., difference in means, heterogeneity of variances, normality of distributions, etc) are typically weakened by virtue of relatively minute perturbations induced by modern attacks. Naive application of these techniques is usually ineffective. We augment the above with a number of fine-grained anomaly-oriented features that seek to magnify various types of observed discrepancies between the distribution layer values associated with natural samples and the distribution of layer values associated with attacked samples. These techniques are summarized in the subsections below.

\subsubsection{Region-based Aggregated Anomaly Counting Features}
Anomaly counting features divide the observed histogram (i.e., estimated probability density function) of the selected layer inputs into three regions:
\begin{itemize}
    \item (L, anomaly detection left tail): $\mu - k\sigma$,
    \item (R, anomaly detection right tail): $\mu + k\sigma$, and
    \item (C, non-anomalous center): $mu - \sigma \le x_i  \le \mu + k\sigma$.
\end{itemize}

Traditionally anomaly detection relies on $k=3$ for identifying single statistically significant outliers in a signal. However, our goal here is to determine whether the population of observed layer values as a whole is anomalous with regards to the expected population of values averaged across $10K$ natural samples. For this reason, our goal is to aggregate indicators across each of these (L, C, R) regions. Moreover, we use a value $k=1$ to magnify small anomalies -- statistically, we expect $65\%$ of all values to fall into C whereas approximately $17\%$ each will fall into regions L and R. For each region. we count the total number of values in the region, the average zscore for these values, the total sum of zscores, and various ratio between regions.
Fig.~\ref{}\rf{what is this supposed to point to?} illustrates that basic idea. Layer inputs are z-scored and considered as a single population. with a histogram built for the zscore values.

\subsubsection{Extreme Value Features}
In addition of these features, given the current layer values, we perform anomaly detection tests that compare these anomaly-detection signal values against for $10$ percentile and $90$ percentile extracted from the population of such signal values across all $10K$ natural samples. These per-region anomaly counting tests examine whether the current layer inputs when considered as a population exhibits extreme values.

\subsubsection{Histogram Bin Count Features}
Histogram binning features represent raw histogram binning counts that estimate the probability density function of the layer inputs $x_i$ across all channels.

\subsubsection{Histogram Tail Symmetry Features}
Symmetry features are used to quantify the symmetry between tail bins. Whereas the probability density function for natural samples is normally distributed, some attack notoriously generate layer inputs for which its probability density function although resembles the expected bell curve exhibits tail symmetry anomalies. These features quantifies differences on bin counts at corresponding positions on both tails of the sampled histogram for the layer inputs $x_i$.

\subsubsection{Statistical Test Features}
Statistical tests are used to generate p-values from the resulting tests (e.g., comparison of means, normality, same distribution, ratio of variances) between (A) a subsample from the population of layer inputs $x_i$ for the underlying image and (B) a subsample from the population of average valuesfor the same layer ---but averaged across the collection of $10K$ natural inputs. These tests examine whether statistical point estimators for the current layer inputs (when considered as a population) are statistically different from the observed estimators computed across the population of $10K$ natural samples.  In general, as the size of subsample being tested grows larger, these tests lose significance. For this reason, tests are conducted using a small size subsample from both populations being tested.

\subsubsection{Similar Shape Probability Distribution Features}
The Wasserstein distance is used to quantify how symmetric are similar are (regions of) binned distributions. In a nutshell, given two k-binned histograms, the Wasserstein distance estimates the work cost of making the skyline on one histogram resemble the other skyline by optimally removing/adding/moving from a mismatching bin into other adjacent bin(s). The larger the metric the larger the "work" (i.e., material over distance) needed to make the distributions look alike. Similar distributions have low distances. Moreover, the distance can be use to analyze symmetry in the tail regions of our distributions.

\subsubsection{Dimensionality Reduction Features}
We also leverage PCA and LDA dimensionality reduction components as features to test whether explicit the induced separability by these projections improve classification accuracy. Fig.~\ref{lda} shows a plot showing the manifold manifested by the first two LDA components for the space spanned by AFVs.

\subsubsection{Nearest Neighbor: Proximity to Known Example Features}
A Radius Nearest Neighbor classifier is used to generate radius neighbor probabilities as features. Let $f_{i}$ be the $i^{th}$ anomaly feature vector with corresponding label $y_i$. Then Radius Nearest Neighbor produces
   $f_i[k]= \sum_{k \in Q}{distance(f_i, f_j) \le r}$.

These features represent the ratio ($\frac{A}{B}$) of (A) nearby class members of some class $k$ over (B) the total of neighboring points with a given radius across all classes in $Q$. Notably, the classification is done using standard scaler inputs to preserve outlier anomalies. Using a grid search, a radius of $r=0.7$ (i.e., standard deviations) was found to deliver the best accuracy. Outliers having no neighbors are defaulted to the most frequent class.

As these features examine the peering radius neighborhood of every sample, they help our classifier identify out-of-distribution examples.  Given $Q+1$ classes in our attack classification problem, a Nearest Neighbor is trained with a subsample of the training set and predictions are computed for the entire training set. Given a large enough subsample, $Q+1$ class probabilities will be generated. These probabilities represent the nearest neighbor votes for a radius $r$ from the given sample. These help the classification support for in-distribution samples as well as help identify out-of-distribution outliers (i.e., those feature vectors having no neighbor counts).

\section{Experiments}
\label{Section:ExperimentalSetup}

We examine DeClaW\xspace classification performance with respect to one dozen different classes: one class $0$ for Natural samples and $11$ attack classes $A_1 \cdots A_{11}$ for attack-perturbed samples obtained from the corresponding application of $11$ different evasion methods. Each of attack class $A_1 \cdots A_{11}$ was generated by applying the attack function of the IBM Adversarial Toolkit over the same subset of natural samples. Default parameters were used in the attacks.  The Appendix lists the attack classes. Table~\ref{table_attack_classes}.
We also examine DeClaW\xspace detection performance with respect to the Class $0$ above for Natural samples and the agglomeration of the above $11$ classes aboves into one class representing attack-perturbed samples by any type of attack method.
For most of the results in this paper, unless noted otherwise as in the case of clustering of attack classes that are too similar to discriminate, we focused on $11$ different attack classes plus 1 class for the natural samples.

For the CIFAR-10 dataset, we used all the $60000$ natural samples to train and test the second stage classifier. Using a subset of natural CIFAR-10 samples $C_6$, we first inferred the reference temperature of the dataset of the natural samples. For each attack class $A_k$, we  collected the resulting anomaly-detection feature vector produced by our hook for the same $C_6$ natural samples. In total, we generated a total of $60000$ anomaly detection feature vectors from the CIFAR-10 natural samples and
a total of $102000$ attack samples across the $11$ attack classes. About $7.3\% = \frac{8000}{110000}$ of attacks failed to be successfully computed within our time bounds and we omitted those from consideration. We used a $70\%$ to $30\%$ split between train and test. More details of experimental setup are in the Appendix.

\label{Section:Experiments}

We use a 32-layer Wide Residual Network~\cite{zagoruyko2016wide} as the classifier for CIFAR-10, with a sampling hook placed before the batch normalization stage of the model's pipeline (details in the Appendix). As a detector, for CIFAR-10, this network’s accuracy on clean test data is $93\%$ and $96\%$ on adversarial data (see Fig.~\ref{fig:binary-cmat}).
This compares well to state-of-the-art baselines. For instance, Metzen et al. report 91.3\% accuracy on clean data and 77\% to 100\% on adversarial data that used same permitted distortion during attacks as the detector was trained on~\cite{metzen2017detecting}. Aigran et al. achieve only 55.7\% false-positive rate when evaluated against FGSM attacks with a 95\% true positive rate~\cite{aigrain2019detecting}.

Moreover, DeClaW\xspace (see Fig.~\ref{fig:cluster-cmat}) achieves high classification accuracy among  attack methods, with similar attacks clustered. For instance, on adversarial versions of test data, it had almost perfect ability to distinguish among PGD and Carlini-Wagner attack images and close to 90\% accuracy in distinguishing PGD attacks from non-PGD attacks. F1-scores, which can be computed from the confusion matrix, ranged from 90\%-99\% for the class groups shown.  To our knowledge, we are the first to report such classification accuracy levels across different attack methods.

 \begin{figure
        \centerline{
        \includegraphics[width=0.35\textwidth]{FIGURES/det-detc-reduced.png}}

        \caption
        Confusion matrix for the detection of adversarial examples where CLEAN bit is set only if no ATTACK class is detected and ATTACK bit is set otherwise.
        \label{fig:binary-cmat}
    \end{figure}

    \begin{figure
        \centerline{
        \includegraphics[width=1\linewidth]{FIGURES/det-class-reduced-plot.png}}

        \caption
        Declaw's classification accuracy among attacks after attack methods are clustered based on large FPR and FNR. Appendix (Fig.~\ref{fig:cluster-cmat2}) is an expanded figure.}
        \label{fig:cluster-cmat}
    \end{figure}

    %

    %

\section{Ablation Studies}
\label{Section:AblationStudies}

We performed ablation studies to isolate the contribution of various techniques used to increase classification performance.

These techniques are
\begin{itemize}
\item dimensionality reduction features
\item histogram comparison features
\item synthetic augmentation of feature vectors
\item range normalization of feature vectors
\item weakened strain boosting of poor-performing attack classes
\item clustering often confounded attack classes
\end{itemize}

To evaluate the performance of the detector and the classifier, we first will focus on F1 score for natural samples ($C0$-$F1$) and the macro-averaged F1 score across the various attacks ($AVG$-$F1$). We also compare the average and maximum accuracies ($MuAcc, MxAcc$) reported for classification as well as detection ($CLF$, $DTC$) modes.

The ablation studies below itemize the contribution of these techniques to increase in accuracy metrics of interest.

\subsection{Impact of Doubling The Data}
 Will our detector respond well to an increase of training data? Next, we doubled the initial 5000 CIFAR-10 samples data, resulting in a total of  $(N=10000)\dot (13+1)$ samples. Our initial baseline used 5000 CIFAR-10 samples and subjected them to the 13 attack classes, generating $(N=5000)\dot (13+1)$ samples for our second stage classifier. The following ablation shows the impact of doubling the available data over accuracy metrics. The remainder of experiments rely on N=10000.

\subsection{Earth Moving Distances}

Wasserstein distances or so-called earth moving distances techniques to estimate work needed to make make-alike two distributions.

In a nutshell, Wasserstein distance is a measure of the optimal effort to make two binned distributions become alike. Lower values represent similarly binned distributions, higher distributions for which significant shaping work is needed to make them alike. We use 23 bins.

The following ablation table compares detection and classification accuracy metrics between basic fine-grained features versus fine-grained features plus distributional shape features obtained using Wasserstein distances (also referred to as earth moving distances). The extra features yield a very modest $1\%$ increase on accuracy metrics.

\input{TABLES/table_extrafeatures}

\subsection{Dimensionality Reduction Features}
We use a fine-grained anomaly detection feature space --- some features being correlated and/or redundant. It is worth to examine if dimensionality reduction techniques will enhance our performance.

Principal Component Analysis (PCA) and Linear Discriminant Analysis (LDA) techniques can be used to introduce proxy features that capture most of the variability or reduce most of the cross entropy. While it is known that these proxy features are a gamble and can in general these features can hurt accuracy metrics, could it possible that they will contribute positively in differentiating between two distributions?

The following ablation table compares detection and classification accuracy metrics between normal, PCA-augmented, LDA-augmented and PCA and LDA augmented runs. Under the non-clustered mode, PCA adds again a modest $1\%$ increase in performance, whether LDA is used or not.

\input{TABLES/table_pcalda}

\subsection{Synthetic Data Augmentation}
Generative techniques to randomly sample a minute epsilon-ball around the point ascribed to an input feature vector. By choosing a very minute epsilon, we speculate the label of these samples to be the same as the input sample and generate additional samples to drive up robustness of the neural network with respect to minute perturbations on the underlying input feature vector.

We investigated whether by synthetically augmenting our dataset, we could increase accuracy metrics. Specifically, we introduced a very small noise perturbation around each feature value of a sample's profile. An adversarial attack essentially perturbs samples but with the optimization goal of switching the classifier's label. On the other hand, we introduce minute noise, but to expand the recognition range of the resulting classifier towards slightly perturbed samples. The perturbations of the (statistical profile) feature vector for the $i^{th}$ sample (with elements $x_{ij}$ was of the form
\begin{equation}
    x^*_{ij} = x_{ij} \pm (\sigma_{j} \* \rho_{j} \* \epsilon + \eta_{ij})
\end{equation}
where $\sigma_j$ corresponds to the standard deviation for the population of the given feature $x_j$, $\rho$ is a random magnitude in the range $[-1, 1]$ and $\epsilon$ the scaling factor of the resulting perturbation (e.g., $10^{-10}$).

The following ablation table compares detection and classification accuracy metrics between the synthetic augmentation at various $\epsilon$ values in the range $[10^{-6}, 10^{-15}]$.

\input{TABLES/table_aug}

\input{TABLES/table_augeps}

\subsection{Robustness via Weakened Strains}
In constructing the attack samples, each attack is applied over the set of natural samples using basic or default parameters. Each natural sample has two possible outcomes: successfully attacked or failed attack. The following table itemizes for each class, how many natural samples were processed, how many were successfully attacked, how many were failed attacks and what was the resulting accuracy metrics for the attack class during evaluation mode.

Each attack could be made harsher or be weakened. Using a metaphor of immunization, where a weakened strain is used to help the immune system ML to learn to detect and target such viruses, for attack classes were we obtained low detection/classification performance, we double the data with similarly generated attack samples but using weakened parameters. The following table shows the increase in accuracy metrics obtained in classification/detection for each of these attacks and the overall increase in classification and detection accuracy metrics.

\subsection{Handling of Failed Attacks}
What to do with failed attack samples?
A failed attack sample could lie anywhere within close proximity to the successful new class label or just the opposite, towards close proximity to the natural sample's label. In the case of the former, it is detrimental to the classification margin/boundary to include this sample as representative of the attack class -- as it pushes the attack membership towards the natural/clean cluster. This problem requires either dropping the failed attacks, including them into the attack class, or a mix of these two based on some clustering criteria.

\subsection{Best Classification Performance}
\input{TABLES/table_best}

\subsection{Experiment Control Parameters}
Declaw allows command line user-control of the application of these various hypothesis that expand the basic anomaly feature vector. These parameters are shown in Table~\ref{tab:controlparameters}.

\input{TABLES/table_controlparameters}

Declaw also provides user control of various basic hyper-parameters. These parameters are shown in Table~\ref{tab:hyperarameters}.

\input{TABLES/table_hyperparameters}

\section{Clustering of Similar Attacks}
\label{Section:ClusteringSimilarAttacks}

DeClaW\xspace learns a statistical envelope for each attack class. However, as some attacks are more sophisticated derivative variants of other attacks, the resulting envelope may be similar.

Attacks are essentially perturbations optimized to reduce visible artifacts against the natural inputs. It is therefore to be expected that given the set of selected distributional features selected to identify those perturbations, some of these attacks would produce distributional profiles too close to each other to clearly differentiate between them. The result is that macro-averaged metrics across such confounded classes introduces losses in recall, precision, and F1 values that distort the perceived quality of the detector and classifier.

Higher detection and classification accuracy can be obtained by treating such clustered classes as one meta attack class. For example, we found that a set of twenty attack classes clustered into ten different similarity clusters of different sizes. Similarly, in this paper, our twelve attack classification classes were found to cluster into just the six different classification clusters identified in Table~\ref{table_attack_classes}.

Our proposed solution is to cluster similar attacks. To do so, we use a transitive closure with respect to ranked FPR and FNR (being above a given tolerance (e.g., 10\%), currently favoring FNR when conflicts arise. For example, given a normalized confusion matrix, the decision whether to cluster to attack classes $A_i, A_j$ is solely decided on the $2x2$ pairwise confusion matrix between those two attack classes --- i.e.,
\[
M(A_i, A_j) =
\left( \begin{smallmatrix}
        acc(A_i) & fnr(A_i, A_j) \\
        fpr(A_i, A_j) & acc(A_j)
        \end{smallmatrix} \right).
\]

Now, we threshold the above matrix, for $fpr$ and $fnr$ values greater than some threshold $t\ge0.2$, we focus solely on large attack classification error. For example, given
$M(A_i, A_j)
= \left( \begin{smallmatrix}
        1 & a \\
        b & 1
        \end{smallmatrix} \right)$
where $a, b$ could either be $0$ or $1$. If any is $1$, attack classes $a, b$ are clustered as one. As a convention, we assign the parent class as the one with the lowest classification label.
The same principle applies for larger number of classes, except that a transitive closure over large classification error is used to drive the merge. As the number of classes is small, this process is done by hand, however, a union-join() or dfs() algorithm can be used to implement the transitive closure for the case when the number of classes to be defended is extremely large or it is desirable that such procedure is automated.
Table~\ref{table:cluster} compares detection and classification accuracy metrics when using the clustered (6 attack clusters) mode. After clustering attacks that were found to be difficult to distinguish based on AFVs, reducing the number of classification labels from 12 to 6, the average classification accuracy increased by $21\%$ while average detection accuracy remained the same.

\input{TABLES/table_clustered}

\section{Attacking the Defense}
\label{Section:AttackingTheDetector}

Finally, we analyze the performance of a white box iterative attack.

Using Carlini amd Warner's approach, we attacked the detector by incorporating into the gradient recomputation of the first state attack, a cross-entropy term against the second stage detector's logit prediction.

\section{Limitations}
Our approach shows promise on CIFAR-10 dataset in distinguishing among several types of attacks. But the approach should be tried on additional datasets to assess  generalizability. Preliminary results on CIFAR-100 dataset are promising and given in the Appendix under Additional Results.

The approach also needs to be evaluated against adversarial attacks where the adversary has knowledge of the DeClaW pipeline. DeClaw pipeline, as it stands, is also not adversarially trained.  We do note that in corporate environments, the DeClaW pipeline need not be made visible to an adversary, nor its results need to be directly visible, since it is primarily used for error detection and error classification rather than as the mainstream classifier (e.g., YouTube or Facebook could use an attack detector/classifier inspired by DeClaW on the backend for analysis without exposing that to end-users.)  Furthermore, DeClaW pipeline could be adversarially trained. We note that for an adversary  to recreate the  DeClaW pipeline, they require both the training data and whitebox access to the target model so as to create AFVs -- a higher threshold than just whitebox access.
\section{Future Work}
\label{Section:FutureWork}

\subsection{Adaptive attacks}
One limitation of this work is the lack of adaptive attack evaluation. Prior work shows that previous detectors can be attacked with adaptive attacks~\cite{carlini2017adversarial}. While we use a lot of statistical profiling in our AFVs, it remains to be seen how an adaptive adversary would perform.

\subsection{Handling of Out of Distribution Inputs}
The validation and/or certification of inputs is of importance to both training as well as evaluation tasks.

\begin{itemize}
    \item During training- it will be desirable to know if samples represent out-of-distribution instances or even if the classification support for these is too small. The former is typically associated with population-based anomaly detection limits while the latter is associated with whether the sample's anomaly-detection feature vector aligned to a known cluster, centroid or representative sample in the feature space.

    \item During evaluation- it will be desirable to reliably distinguish between samples suspected to be compromised by an attack and decide whether to reject the sample, log and trace the anomaly, ignore the attack or perhaps reroute the sample to a more robust classifier.
\end{itemize}

\subsection{Concept Shift}
Our approach relies on a single reference temperature (i.e., representative sample or centroid) for the entire natural samples dataset. However, this could be made more robust. Specifically, two simple but very important enhancements could be made.
\begin{itemize}
    \item removal of outliers: the dataset of natural samples may contain extreme outliers that could bias or even nullify the computation of the reference temperature for the overall dataset (i.e., a series of fine-grained centroid computations). We could first perform a first pass over the dataset of natural samples and compute an initial reference temperature, then identify any anomalous natural samples with respect to this initial reference temperature (e.g., plus/minus three sigma from the our first reference temperature estimators) and then recompute the final reference temperature on a second pass over the set of non-anomalous natural samples.

    \item multimodal distributions: the dataset of natural samples could manifest a multimodal distribution with respect to our features. Visually, this  represents the case where the anomaly detection feature vectors for natural samples do not cluster into one well-defined Gaussian multivariate cluster but rather into two or more clusters. A representative sample for both could be misleading or meaningless (e.g., cluster one represents daylight conditions and cluster two represents nighttime conditions). Trivially, one could use multiple detectors for expressing these ultimately different concepts. Alternatively, our second stage network architecture could be extended to explicitly learn such user-provided or domain-related features of importance to the classification.

\end{itemize}

\subsection{Large Number of Classes}
There is no reason for generating and partitioning anomaly feature vectors classification classes solely on attack method. Thanks to our clustering methodology, it is possible to train a large cardinality classifier where dozens classes to express variants of a given attack method are used (e.g., $PGD3, PGD10, PGD20, PGD40$, etc). Similarly, it might be desirable to express a class for failed versions of attacks, which are expected to have a different distribution to that of successful attacks. Same principles apply to addressing multimodal natural samples.

\subsection{Healing of damaged or attacked samples}
By leveraging the concepts of a representative sample and a normative subset of natural samples, we can explore the concept of automatically healing anomalous values in the observed layer values $l^*_i$ that exceed a strong anomaly $hi$ watermark (e.g.
$\mu + 3\sigma$). The autonomous healing comprises a weighted function of the actual pixel value plus the corresponding element wise value from the normative baseline.
\section{Conclusion}

\label{Section:Conclusion}

We present DeClaW\xspace, a system to detect, classify and warn of adversarial examples. We achieve high detection accuracy on CIFAR-10 (96\% accuracy on adversarial data and 93\% on clean data) while being the first to our knowledge to classify the type of attack. We are able to distinguish among many of the 11 different attacks on CIFAR-10 with F1-scores ranging from 90\%-99\%.

\section*{Acknowledgment}

We acknowledge research support from Ford and University of Michigan. This research is partially based upon work supported by the National Science Foundation under grant numbers 1646392 and 2939445 and DARPA award 885000.

\subsection{Attack Classes}
Table~\ref{table_attack_classes} shows the attack classes that we trained and evaluated DeClaw on.

\input{TABLES/table_split_11_classes}

\subsection{Experimental Setup}
For the CIFAR-10 dataset, we used all the $60000$ natural samples to train and test the second stage classifier. Using a subset of natural CIFAR-10 samples $C_6$, we first inferred the reference temperature of the dataset of the natural samples.\footnote{In the Future Work section, we identify issues with generating a reference temperature of a dataset such as the presence of concept drift or multimodality.} For each attack class $A_k$, we  collected the resulting anomaly-detection feature vector produced by our hook for the same $C_6$ natural samples. In total, we generated a total of $60000$ anomaly detection feature vectors from the CIFAR-10 natural samples and
a total of $102000$ attack samples across the $11$ attack classes.\footnote{About $7.3\% = \frac{8000}{110000}$ of attacks failed to be successfully computed within time bounds.} We used a $70\%$ to $30\%$ split between train and test.
-- resulting in second stage datasets as follows:
\begin{itemize}
\item training with $42000$ and testing with $18000$ anomaly detection feature vectors for natural class $C$ and
\item training with $72200$ and testing with $29790$ anomaly detection feature vectors across attacked classes $A_1 \dots A_{11}$.
\end{itemize}

We performed the experiments using $12$ classes:
\begin{itemize}
\item one class $C$ of anomaly detection feature vectors from natural samples
(divided into subsets $C_1 \dots C_6$, each of $10000$ samples).
\item
and $11$ classes $A_k$, each resulting from the application of a given attack class $A_k$
over the same ($10000$ sample size) subset $C_6$ of the natural samples.
\end{itemize}

DeClaW\xspace provides user control of various basic hyper-parameters. These parameters are shown in Table~\ref{tab:hyperarameter}.

\input{TABLES/table_hyperparameters}

We use a 32-layer Wide Residual Network~\cite{zagoruyko2016wide} as the classifier for CIFAR-10. The network has been trained for 80 epochs with stochastic gradient descent and momentum on 10000 data points from the train set. The momentum term was set to 0.9 and the initial learning rate was set to 0.1, reduced to 0.01 after 50 epochs, and further reduced to 0.1 after 70 epochs. After each epoch, the network’s performance on the validation data (the remaining
5000 data points from the train set) was determined. The network with maximal performance on
the validation data was used in the subsequent experiments (with all tunable weights being fixed).

We attach an adversary detection
subnetwork (called “detector” below) to the ResNet. The detector is a convolutional neural network
using batch normalization~\cite{ioffe2015batch} and rectified linear units. At this time, we focused on placing our sampling hook into the batch normalization stage of the original model's (e.g., CIFAR-10) pipeline.

\subsection{Clustering of similar attack types}
\input{11_clustering}

\subsection{Feature Generation Strategy}
Rather than training our second stage network with pixels as features,
we trained the second stage classifier using anomaly-detection features.
For each input sample $x_i$, we examined the resulting input tensor $l_{i}$
presented to the sampling hook layer.  That is, $l_i$ represents the layer
inputs to the final batch normalization layer shown in Table~\ref{samplinghook}.

Given that our CIFAR-10 input tensor $l_i$ to the batch normalization layer contains 640 channels of $8\times 8$ features, our sampling hook observes flattened vectors $v_i$ containing $40960$ values. We compare statistical features of these vectors derived from subsampling projections from it and comparing the outlier distribution against normative mean, variance, and count references. By examining the observed features location relative to channel-wise means $\mu_c$ and standard deviations $\sigma_c$, we identify outliers that contribute to an image being clean or adversarial, and if it's adversarial, which type of attack it is.

In total, we consider eight categories of features: region-based, extrema-based, histogram counts, histogram bin comparisons, probability tests, Wasserstein distances, PCA/LDA dimensionality reduction features, and nearest neighbor features.

\subsection{Basic Building Blocks for Anomaly Features}\label{sec:buildingblocks}

The input tensor values $l_{i}$ have a shape of $[1,640,8,8]$
--- meaning each sample has $640$ channels of $64$ values each.
By considering each tensor $l_i$ as a vector in $R^{40960}$ of random variables,
we estimate means and variances within each of the $640$ channels as well as across all the $640$ channels of $l_i$.
For convenience, let $L = len(l_i)$, -- that is, $L=640\times8\times8=40960$ for CIFAR10.

We first estimate the baseline channel-wise means and standard deviations, notated as $\mu_c$ and $\sigma_c$ for a channel $c$. For a given input tensor $l_i$ and its values in channel $c$ notated as $l_{i,c}$, we define its mean $\mu_c(l_i) = \frac{\sum{l_{i,c}}}{L}$ and standard deviation $\sigma_c(l_i) = \sqrt{\frac{\sum{(l_{i,c} - \mu_c(l_i))^2}} { L  }}$. Then, similar to the process in the original batch normalization paper~\cite{ioffe2015batch}, we compute running means and standard variances for each channel with an exponential smoother to compute our baseline means. For a channel $c$, let $\mu_{base,c}$ refer to the baseline mean and let $\sigma_{base,c}$ refer to the baseline standard deviation.

We hypothesize that using aggregation schemes to accumulate
outliers would yield stronger anomaly detection signals.
To this end, we set $lo$ and $hi$ watermarks
with respect to distribution of values in a channel $c$.
Trivially, these watermarks are (as is commonly done) implemented as $\sigma$-based control limits -- specifically,
$lo_c = \mu_{base,c} -  \sigma_{base,c}$ and
$hi_c = \mu_{base,c} +  \sigma_{base,c}$.

Given $lo$ and $hi$ control limits or watermarks for each channel,
we create three indicator vectors $S_{lo}, S_{mi}, S_{hi}$ for all $l_i$ that places each feature into one of those three disjoint sets. Let $v_i$ refer to a flattened vector of $l_i$ and $z_i$ refer to a z-scored version of $v_i$. For any given index $k$, and channel $c$ corresponding to the channel that $v_i[k]$ belongs to, the value of $S_{lo}[k] = 1$ if $v_i[k] \leq \mu_{base,c} - \sigma_{base,c}$ and is 0 otherwise. Likewise, the value of $S_{hi}[k] = 1$ if $v_i[k] \geq \mu_{base,c} + \sigma_{base,c}$ and is 0 otherwise. The value of $S_{mi}[k] = 1$ if $\mu_{base,c} - \sigma_{base,c} < v_i[k] < \mu_{base,c} + \sigma_{base,c}$ and is 0 otherwise. Finally, we let an overall outlier vector $S_{ov}$ be a vector where its value is 1 where either $S_{lo}$ or $S_{hi}$ is 1 and 0 otherwise.

Using the above $4$ indicator vectors, we generate $2$ basic aggregation-based anomaly features for each vector $S_x$:
\begin{itemize}
        \item $C(S_x) = \sum\limits_{\forall k \in [1, L ]}{S_x[k]}$
        \item $Z(S_x) = \sum\limits_{\forall k \in [1, L ]}{S_x[k] \cdot z_i[k]}$
\end{itemize}
For any input tensor $l_i$, these $8$ scores are computed and then used
as building blocks to generate features used by DeClaW\xspace.

Some features also compare against values from a reference dataset of representative samples. Let $D_N$ refer to such a dataset, formed from 10k randomly chosen values from the natural samples in the training set.

\subsection{Region-based Features}\label{sec:region}

With the building block and aggregatation features described in Section~\ref{sec:buildingblocks}, we specify region-based features. These region-based features count and measure the signal of outliers in the left and right regions of the distribution curve. There are 25 features derived from the $C(S_x)$ and $Z(S_x)$ building blocks as described in Table~\ref{basic_anomaly_features}.

\begin{table}[h]
    \centering
    \begin{tabular}{|l|l|}
    \hline
    Feature & Region-based \\
    name & feature specification \\
        \hline
        \hline
        $\text{error density}          $ &  $ f_{01} = \frac{C(S_{lo}) + C(S_{hi})} {L } $ \\
        \hline
        $\text{RIGHT error count}      $ & $ f_{02} = C(S_{hi}) $\\
        $\text{LEFT error count}       $ & $ f_{03} = C(S_{lo})$ \\
        $\text{OVERALL error count}    $ & $ f_{04} = C(S_{ov})$ \\
        \hline

        $\text{RIGHT/OVERALL}          $ & $ f_{05} = \frac{C(S_{hi})} {L } $ \\
        $\text{LEFT/OVERALL}           $ &  $ f_{06} = \frac{C(S_{lo})} {L } $ \\
        $\text{(RIGHT - LEFT)/overall}   $ & $ f_{07} = \frac{C(S_{hi}) - C(S_{lo})} {L } $ \\
        \hline

        $\text{RIGHT anomaly signal}   $ & $ f_{08} = Z(S_{hi})$ \\
        $\text{LEFT anomaly signal}    $ & $ f_{09} = Z(S_{lo})$ \\
        $\text{OVERALL anomaly signal} $ & $ f_{10} = Z(S_{ov})$ \\
        \hline

        $\text{RIGHT norm signal}      $ & $ f_{11} = \frac{Z(S_{hi})} {L}$ \\
        $\text{LEFT norm signal}       $ & $ f_{12} = \frac{Z(S_{lo})} {L} $\\
        $\text{OVERALL norm signal}    $ & $ f_{13} = \frac{Z(S_{ov})} {L} $\\
        \hline

        $\text{ANOMALY error}          $ & $ f_{14} = Z(S_{lo}) + Z(S_{hi})$ \\
        $\text{WITHIN anomaly error}   $ & $ f_{15} = Z(S_{mi}) $\\
        $\text{OVERALL error}          $ & $ f_{16} = Z(S_{ov})$ \\
        \hline

        $\text{ANOMALY/WITHIN ratio}   $ & $ f_{17} = \frac{Z(S_{lo}) + Z(S_{hi})} { Z(S_{mi}) }$ \\
        $\text{ANOMALY/num errors}     $ & $ f_{18} = \frac{Z(S_{lo}) + Z(S_{hi})} { C(S_{lo}) + C(S_{hi}) }$ \\
        $\text{WITHIN/(n - num errors)}$ & $ f_{19} = \frac{Z(S_{mi})} { C(S_{mi}) + 1 } $\\
        \hline

        $\text{normalized WITHIN area}    $ & $ f_{20} = \frac{Z(S_{mi})} { L }$ \\
        $\text{WITHIN/OVERALL ratio}$       & $ f_{21} = \frac{Z(S_{mi})} { Z(S_{ov} }$ \\
        $\text{norm. OVERALL error ratio} $ & $ f_{22} = \frac{Z(S_{ov})} { L }$ \\
        \hline

        $\text{average ANOMALY score}   $ & $ f_{23} = \frac{Z(S_{lo}) + Z(S_{hi})} { C(S_{lo}) + C(S_{hi}) } $\\
        $\text{average WITHIN score}    $ & $ f_{24} = \frac{Z(S_{mi}) } { C(S_{mi}) } $\\
        $\text{average OVEREALL score}  $ & $ f_{25} = \frac{Z(S_{ov}) } { C(S_{ov}) } $\\
        \hline

    \end{tabular}
    \caption{The basic anomaly detection building block features based on region-based aggregations.}
    \label{basic_anomaly_features}
\end{table}

\subsection{Extreme Value Features}
The above region-based features in Section~\ref{sec:region} deal with central and fat portions of the tails of the distribution of $l_i$ values.
Next, we develop features that isolate extrema tail values of the distribution of $l_i$ values.
Let $R=25$ be the number of region-based features above.
For each feature $f_r, \forall r \in [1,25]$, we note individual feature extrema values. For a given $f_r$, let $\phi_{lo}(f_r)$ and $\phi_{hi}(f_r)$ be extrema percentile thresholds (e.g., $lo = 10\%$ and $hi = 90\%$). That is, $10\%$ of feature values $f_r(l_i^N)$ from the population of samples $l_i^N$ from $D_N$ are less than $\phi_{lo=10\%}(f_r)$.
Similarly, $90\%$ of feature values $f_r(l_i^N)$ are less than $\phi_{hi=90\%}(f_r)$.
Finally, we define indicator flags $e_r(l_i)$ that identify, for $l_i$, whether or not the
feature value $f_r(l_i)$ belongs to the normative population of $f_r$ values.
\[ e_r(l_i) =
        \begin{cases}
              1,             & \text{if } f_r(l_i) \ge \phi_{hi=90\%}(f_r) \\
              1,             & \text{if } f_r(l_i) \le \phi_{lo=10\%}(f_r) \\
              0,             & \text{otherwise} \\
         \end{cases}
\]
$\forall r \in [1, R]$.

For each feature value $f_r(l_i)$ of every layer input $l_i$, values
within the range $[(\phi_{lo}(f_r), \phi_{hi}(f_r))$ will be considered normative values
and values outside will be flagged as outliers. After computing these $R$ boolean
flags for every layer input values $l_i$, we compute a normative score as simply:
\[ \text{score}(l_i) = \sum\limits_{\forall r \in [1,R]} {e_r(l_i)}.
\]

We also add a normalized score, which divides the normative score by the number of region-based features $R$. A total of $25$ extra features (see Table~\ref{tab:my_label}) are thus generated to flag unusually anomalous feature values.

\begin{table}[h]
    \centering
    \begin{tabular}{|l|l|l|}
    \hline
    Feature & Region-based \\
    name  & specification \\
        \hline
        \hline
        $\text{normative score:}          $ & $ f_{26} = \text{score}(l_i)$  \\
        $\text{normalized score:}         $ & $ f_{27} = \frac{\text{score}(l_i)}{R}$  \\
        \hline
        $\text{flag:} e_{01}(l_i)         $ & $ f_{28} = e_{01}(l_i) $ \\
        $\text{flag:} e_{02}(l_i)         $ & $ f_{29} = e_{02}(l_i) $ \\
        $\text{flag:} e_{03}(l_i)         $ & $ f_{30} = e_{03}(l_i) $ \\
        \hline
        $\cdots                           $ & $ f_{31} \cdots f_{50} $ \\
        \hline
    \end{tabular}
        \caption{Twenty five additional features are used to identify whether a given feature value $f_r(l_i)$ from
        Table~\ref{basic_anomaly_features} represent extrema events with respect to the overall distribution and percentiles of feature values
        obtained for normative natural samples.}
    \label{tab:my_label}
\end{table}

\subsection{Basic and Comparative Histogram Features}\label{sec:appendix_hist}

\begin{figure}[t]
\centerline{\includegraphics[width=0.5\textwidth]{FIGURES/histogram_binning_anomaly_detection.png}}
\caption{To compare the observed distribution $z_i$ values against the expected distribution of $z^N_i$ values
        histograms are used. We use 23 bins and construct features that measure discrepancies between observed
        and expected bincounts as well as symmetry between selected regions $R_x$ and $R_y$ of the compated distributions.}
\label{binning}
\end{figure}

In addition to the above $R$ region-based anomaly features,
to train the neural network to understand fine-grained differences
between the estimated $pdf$ of observed layer inputs $l_i$
and 10k randomly chosen normative layer inputs $l^N_i$ from the training set in $D_N$,
histogram bincount features were used over the standarized values $z_i$ of $l_i$, where $z_i$ is computed by subtracting the mean from $l_i$ and dividing by the standard deviation.

Until now, each of the features above represent features evaluated solely
with respect to the observed layer inputs $l_i$.
Next, we develop a set of comparative (logical/boolean, differential, ratio)-based features,
that seek to train the neural network to discern the observed distribution of $l_i$ values
against the distribution of normative $l^N_T \in R^{L}$ layer input values
--- extracted by element-wise averaging across the $D$ normative natural samples in $D_N$:
\[
        l^N_T = \frac{\sum\limits_{\forall i \in [1,D]}{l^N_i}}{D}.
\]
That is, the average $l^N_T$ input layer values for natural samples at the preselected batch normalization layer
is a vector of size $L$ which is the element-wise arithmetic mean of every $l^N_i$ from the normative set of $D_N$.

Given any vector of $l_i$ values of size $L=40960$,
a histogram of $B=23$ bins is computed using the range $[-3, 3]$, resulting in
intervals of size $q = \frac{6 \dot \sigma}{B}$.
The resulting $B=23$ bincounts ($h_{01}(l^*_i) \cdots h_{B}(l^*_i)$) represent $23$ additional features.\footnote{These
features provide a count of the number of standarized $z_i$ that fall into each of the $B=23$
bins that span the interval $[-3, 3]$ of the normative range of standarized values.}

Moreover, just as we generate $B=23$ histogram bincount features ($h_{01}(l_i) \cdots h_{B}(l_i)$)
for observed layer input values $l_i$,
we also compute the histogram bincounts ($h_{01}(l^N_i) \cdots h_{B}(l^N_i)$) for the normative layer input values.
Then, using these two sets of $B=23$ histogram bincounts (i.e., observed bincounts against normative bincounts),
we compute features that magnify dissimilarity (specifically, square difference as well as relative error) between
observed bincounts for a sample and the reference normative bincounts from natural samples.
In all, the following three sets of $B=23$ additional features are added:
\begin{itemize}
        \item ($B=23$ features  $f_{51} \cdots f_{74}$:) the OBSERVED bincounts for observed layer values $l_i$:
 $h_{01}(l_i) \cdots h_{B}(l_i)$,
        \item ($B=23$ features  $f_{75} \cdots f_{98}$ :) the SQUARE DIFFERENCE between observed and normative bincounts:
 ${(h_{01}(l_i) - h_{01}(l^N_T}))^2$ $\cdots$ $({h_{B}(l_i) - h_{B}(l^N_T))}^2$,
        \item ($B=23$ features  $f_{99} \cdots f_{122}$:) the RELATIVE ERROR between observed and normative bincounts:
 $\frac{h_{01}(l_i) - h_{01}(l^N_T)}{h_{01}(l^N_T)}$ $\cdots$ $\frac{h_{B}(l_i) - h_{B}(l^N_T)}{h_{B}(l^N_T)}$
        \item and ($3$ features $f_{123} \cdots f_{125}$:) the sum, mean, and variance of the RELATIVE ERROR above.
\end{itemize}

The fine-grain targeting scope of our binning features is illustrated in Fig.~\ref{binning} by juxtaposing along the three previously discussed coarse-grain aggregation regions.

\subsection{Comparative A/B Probability Test Features}
A couple of A/B tests are used to compare the size $L=40960$ populations of
observed $l_i$ layer values from normative layer values.
We examine whether observed values are normally distributed,
whether observed values have a similar mean to that of normative values,
whether channel means of observed values are similar to the means of channels in normative values,
whether observed values have the same distribution as normative values,
and whether observed values have a similar variance to that of normative values. These features are described in Table~\ref{ab}.
\begin{table}[t]
    \centering
    \begin{tabular}{|l|l|l|}
    \hline
        Feature  & Feature \\
        name  & description \\
    \hline
    \hline
        $\text{normally distributed?} $ & $ f_{124} = \text{KolmSmirnov-test}(l_i, l^N_T)$  \\
        \hline
        $\text{same channel means?}   $ & $ f_{125} = \text{mannwhitneyu-test}(l_i, l^N_T)$  \\
        \hline
        $\text{same pop. means?}      $ & $ f_{126}  = \text{t-test}(l_i, l^N_T)$  \\
        \hline
        $\text{same distribution?}    $ & $ f_{127} = \text{ks-test}(l_i, l^N_T)$  \\
        \hline
        $\text{same variance?}        $ & $ f_{128} = \text{bartlett-test}(l_i, l^N_T)$  \\
        \hline
    \end{tabular}
    \caption{P-values resulting from the application of several statistical tests are used as features.
        These tests compare a random small subsample of the population of $40960$ observed layer input values
        against a random small subsample of the population of $40960$ normative layer input values.}
    \label{ab}
\end{table}

\begin{figure}[h]
\centerline{\includegraphics[width=0.5\textwidth]{FIGURES/plot_lda_visualization_clustered.png}}
\caption{Plot of the first two standarized dimensions resulting from the LDA dimensionality reduction. Note how the $LDA$ first two dimensions provide some signal separation between natural samples (i.e., grey squares) and several attack classes? Can dimensionality reduction features help us discriminate anomalies between natural and attacked samples?}
\label{lda}
\end{figure}

\subsection{Wasserstein Distance Features}
Given the bin-counts from Section~\ref{sec:appendix_hist},
                the earth moving distances (i.e., Wasserstein distances)
                compute how much
                work is needed to make look alike the skyline of a histogram $H_1$ to the skyline
                of a reference histogram $H_0$. Here, $H_1$ represents the histogram for the distribution of observed
                layer input values for a given sample and $H_0$ represents the histogram for the distribution of
                normative layer input values -- that is, extracted from natural samples. Our Wassertein histogram
                bin-count comparison features, built using features $f_{99}$ through $f_{122}$, compare:
        \begin{itemize}
        \item ($H_1$) Left Tail from Observed Distribution for the layer inputs of a given sample against ($H_0$) Left Tail of the distribution across all natural samples.
        \item ($H_1$) Right Tail from Observed Distribution for the layer inputs of a given sample against ($H_0$) Right Tail of the distribution across all natural samples.
        \item ($H_1$) Center Tail from Observed Distribution for the layer inputs of a given sample against ($H_0$) Center Tail of the distribution across all natural samples.
        \end{itemize}
                This hypothesis tests whether the shape of the $pdf$ carries anomaly detection value.

\subsection{Dimensionality Reduction Features}
                We generate PCA(ndim=2) and LDA(ndim=2) features to test whether separability features improve classification accuracy. Fig.~\ref{lda} shows a xy-plot of AFVs using the first two LDA features. It illustrates the spatial clustering patterns across samples of different attack methods and whether additional separatability features ought to be introduced and for which attacks.

\subsection{Nearest Neighbor Classification Support Features}

Given a subsample of generated AFVs for both natural samples and samples perturbed by various attack methods,
     we trained a Radius Nearest Neighbor classifier (RNN).
     Then, for any sample to be evaluated, we estimated dataset support for classification by means of
     the RNN's class voting probabilities and embedded them as features to the DeClaW\xspace Neural Network.
 The RNN classifier
                computes the number of samples of each class in the $(r=\epsilon)$-ball
                neighborhood of the anomaly feature vector to be classified.
                For CIFAR-10, we used $r=3$ and $5K$ anomaly feature vectors.

\section{Clustering of Similar Attacks}

\subsection{Additional Results}

For the CIFAR-10 dataset, Fig.~\ref{fig:cluster-cmat2} reproduces Figure~\ref{fig:cluster-cmat} at a larger scale as well as incorporating data for the CLEAN and APGD classes.  Fig.~\ref{fig:roc} shows a xy-plot of TPR vs FPR values observed across hundreds of runs for CIFAR10 data. It shows that Declaw detection metrics consistently manifest desirable low FPR and high TPR -- across the space of our grid parameter search.

    \begin{figure*}[t]
        \centerline{
        \includegraphics[width=.9\textwidth]{FIGURES/det-clf-cnat-clustered-cropped-redone.png}}

        \caption{
        Confusion matrix for the classification results for the CIFAR-10 dataset produced by Declaw after the 12 classes above are clustered into six cluster groups based on large FPR and FNR.}
        \label{fig:cluster-cmat2}
    \end{figure*}

       \begin{figure
        \centerline{
        \includegraphics[width=0.45\textwidth]{FIGURES/ROC_c.png}}
        \caption{Plot of averaged TPR vs FPR values for the CIFAR-10 dataset across hundreds of different runs resulting from our grid parameter search. Blue diamonds correspond to results obtained with a learning rate of $0.01$ and black squares correspond to values obtained using a learning rate of $1.00$. Declaw detection consistently plots in the ROC sweetspot (i.e., low FPR, high TPR).}
        \label{fig:roc}
    \end{figure}

  Finally, DeClaW\xspace was preliminarily applied to the CIFAR-100 dataset. Without tuning, we achieved detection accuracy on CIFAR-100 (88\% accuracy on adversarial data and 79\% on clean data). As before, we were able to distinguish among many of the 11 different attacks on CIFAR-100 with F1-scores ranging from 74\%-96\%. Figs.~\ref{fig:binary-cmat-c100}-~\ref{fig:cluster-cmat-c100} show the corresponding detection and classification confusion matrices for the CIFAR-100 experiment.

\begin{figure}[p
        \centerline{
        \includegraphics[width=0.35\textwidth]{FIGURES/cmat-detect-c100.png}}
        \caption{
        Confusion matrix for the detection of adversarial examples for the CIFAR-100 dataset where CLEAN bit is set only if no ATTACK class is detected and ATTACK bit is set otherwise.
        \label{fig:binary-cmat-c100}
    \end{figure}

    \begin{figure*}[p
        \centerline{
        \includegraphics[width=0.8\textwidth]{FIGURES/cmat-classify-c100.png}}
        \caption{
        Confusion matrix for Declaw's classification accuracy for the CIFAR-100 dataset among attacks after attack methods are clustered into clusters based on large FPR and FNR.}
        \label{fig:cluster-cmat-c100}
    \end{figure*}


\section*{Methods}

\subsection*{One-dimensional Nonlinear Schr\"odinger equation, variables and parameters}
In the context of uni-directionnal deep water waves, 1DNLSE reads:
 \begin{equation}
i\frac{\partial \psi}{\partial z}= \frac{1}{g}\frac{\partial^2 \psi}{\partial t^2}+\gamma |\psi|^2\psi, \label{eq:NLSE}
\end{equation}
where $z$ is the propagation distance, $t$ is the time measured in the frame co-propagating a the group velocity $c_g=g/(2 \omega_0)$  evaluated at the carrier frequency $f_0=\omega_0/(2\pi)$. In the limit of infinite depth, the third order nonlinear coefficient is:
\begin{equation}
\gamma=k_0^3
\end{equation}
where $k_0=\omega_0^2/g$ is the modulus of the wavevector of carrier wave. The surface elevation $\eta(t,z)$ is related at the leading order to the  slowly varying complex envelope $\psi$ as:
\begin{equation}
\eta(t,z)=\frac{1}{2}\left( \psi(t,z)e^{i (k_0z -\omega_0 t)}+c.c.\right); \label{etaA}
\end{equation}

From (\ref{etaA}) the following relation can be derived :
\begin{equation}
\langle|\psi|^2\rangle=2\langle\eta^2\rangle=2 \sigma^2, \label{amplitudesurface}
\end{equation}
where $\langle...\rangle$ implies averages over time and $\sigma^2$ is the variance of the rapidly oscillating wave field.

In order to design the experiment, it is useful to introduce  a linear and a nonlinear propagation length  as follows:
\begin{equation}
z_{lin}=\frac{g}{\Delta \omega^2}\;\;\ {\rm and}\;\;\;\;
z_{nlin}=\frac{1}{\gamma P_0},
\label{eq:zlznl}
\end{equation}
where $\Delta \omega=2\pi \Delta f$ is a typical spectral bandwidth  and $P_0=\langle|\psi_0|^2\rangle$ is the average value of the envelope square both calculated at $z=0$.

The degree of nonlinearity of the wave propagation is given by the parameter:
\begin{equation}
\Gamma=\frac{z_{lin}}{z_{nlin}}=\frac{\gamma g P_0}{\Delta \omega^2}.
\end{equation}

In the context of ocean waves, $\Gamma=BFI^2$ where BFI is the Benjamin-Feir Index \cite{Janssen:03,Onorato:01}. It has been shown that initial conditions with a large value of $\Gamma$ will eventually lead to the formation of rogue waves characterized by heavy-tailed probability density function of the surface elevation~\cite{Walczak:15,Koussaifi:18,Suret:16,Tikan:18}.

BFI index can be expressed as follows:
\begin{equation}
BFI = \frac{\epsilon}{(\Delta f /f_0)},
\end{equation}
where $\epsilon= k_0 \sqrt{2} \sigma$ is the wave steepness, and $\Delta f$ and $f_0$ are average spectral width and central frequency of the initial wave packets.
\\

\subsection*{Resonances in the 1DNLSE}

The non trivial resonances are forbidden in integrable turbulence~\cite{Suret:11}. We consider the third order nonlinear interaction of four monochromatic waves $i=1,2,3,4$ of pulsation $\omega_i$ and wavenumber $k_i(\omega_i)$ in a unidirectional dispersive media. In this context, the resonances conditions of four wave mixing in the standard wave turbulence read $\omega_1+\omega_2=\omega_3+\omega_4$ and $k_1+k_2=k_3+k_4$ where $k_i=k(\omega_i)$. The linear dispersion relation of deep water waves and of the Eq.~({\ref{eq:NLSE}}) is:

\begin{equation}
    k(\omega)=\frac{\omega^2}{g}
    \label{eq:kL}
\end{equation}

Exact resonances thus lead to:

\begin{equation}
     \omega_1(\omega_3-\omega_2)=\omega_3(\omega_3-\omega_2)
\end{equation}

{\it i.e.} $\omega_1=\omega_3$ and $\omega_2=\omega_4$ or $\omega_2=\omega_3$ and $\omega_1=\omega_4$.

Exact resonances of non trivial interactions ($\omega_1\ne\omega_3$ and $\omega_1\ne\omega_4$) are thus forbidden.

\subsection*{Initial conditions in simulations and experiments: Partially coherent waves}

In the numerical simulations and in the experiments, the initial conditions are partially coherent waves~\cite{Randoux:16}. The slow varying amplitude of the initial condition at $z=0$ reads:

\begin{equation}
    \psi(t,z=0)=A_0 \sum_{l=-N/2}^{+N/2} e^{-\frac{1}{2}(f_l/\Delta f)^2} e^{i \,2 \pi\,f_l\,t} e^{i\phi_l}
\end{equation}

where $f_l=l/T_{max}$, $T_{max}$ is the temporal duration of the experiments and $\phi_l$ are independently and randomly distributed over $[0,2\pi]$.

 Note that the statistical characteristics of partially coherent waves are very different than plane waves initially perturbed by noise which is also investigated in \cite{Leisman:EDR:19}. For a comparison between the two cases, refer {\it e.g.} to \cite{Copie2020}.

\subsection*{Numerical Simulations}
Numerical simulations of Eq.\ref{eq:NLSE} are realized using step-adaptive high order Runge-Kutta method. We construct different initial conditions using the random phase approach where a uniformly distributed phase is added to every Fourier component of a Gaussian spectrum with $\Delta f=$ 0.2 Hz. Typical temporal windows used in the numerical simulations correspond to 100 seconds.  Three parameters of simulations depend on the value of the steepness : the number of points $N$, the length of propagation $L_{max}$ and the number of realizations $N_{sample}$. We separate the numerical studies into three ranges of steepness $\epsilon$:

\begin{tabular}{llll}
    & $\;\;\;N\;\;\;$ & $L_{max}(m)$ & $N_{sample}$ \\
  $\;\;\;\;\;\;\;\;\;\;\;\epsilon\le 0.05$ & $\;\;2048\;\;$ & $\;\;\;500$ & $\;10000$ \\
  $0.05<\epsilon\le 0.09$ & $\;\;2048\;\;$ & $\;\;\;500$ & $\;\;\;\;500$ \\
  $0.09<\epsilon\le 0.19$ & $\;\;1024\;\;$ & $\;\;2000$ & $\;\;\;\;\;100$ \\
\end{tabular}\\

In order to reconstruct numerically NDR, we multiply the spatiotemporal diagram by a Super-Gaussian window with power 15 along z direction avoiding thereby undesirable effects related to the Fourier analysis of non-periodic signals.

\subsection*{Resolution of the measurement of $k$}

 The key point in our approach is to remove the carrier wave before computing the NDR. This allows us to reveal the details of the NDR of the slowly varying envelop. Note that we only use 20 gauges in the water tank. They are separated by $6$~m and the maximum measurable wavevector is around $1.05$~m$^{-1}$. As a consequence, contrary to Taklo {\it et al.} who use 384 probes, we do not resolve the wavevectors of the carrier wave $2\pi/\lambda_0\simeq5.3$~m$^{-1}$ and of the harmonics. Our strategy enables the accurate  measurement of the NLDR of the slow varying envelop of the wave by using only 20 probes. This provides an effective measuring range of $120$~m with a resolution of the measurement of the $k-$spectrum  of $\Delta k_{min}=2\pi/120$~rad~m$^{-1}$. Note finally that, in~\cite{Taklo2015}, the accuracy of the  measurement of $k$ given by the length of the water tank is $\Delta k_{min}/k_0=0.014$ while our setup enables an accuracy of $\Delta k_{min}/k_0 < 0.01$.

In order to measure the averaged spectra and NDR, we use 3,6 and, 3  experimental runs with a duration of 512~s for $\Gamma=0.12$, $0.33$, $0.65$, respectively. One run of 128~s have been used for $\Gamma=6.18$.

Note that in NLS and HOS simulations, the chosen lengths of propagation depend on the parameters and vary typically from 300 to 500~m. The uncertainty of measurement of $k_M$ and $\Delta k$ is therefore significantly lower in simulations than in experiments.

\subsection*{Evaluation of the full width at half maximum $\Delta k$ of the position of the maximum $k_M$}

In Fig. 4, we report the evaluation of the full width at half maximum $\Delta k$ and of the position of the maximum $k_M$ of the function $f(k)=|\widetilde{\psi}(k,\omega=0)|^2$ in 1DNLSE, HOS simulations and in experiments. The accuracy of the measurement of $\Delta k$ and $k_M$ is limited both by  the discretization of $k$ (see above the uncertainty $\Delta k_{min}$) and by the random fluctuations of $f(k)$. In order to overcome these difficulties, when it is appropriate, we evaluate $\Delta k$ and $k_M$ by using best fitting procedure with Lorentzian function.

\newpage
\section*{Supplementary Material}

The purpose of this Supplementary Information is to provide some mathematical and numerical details on the High-Order Spectral (HOS) method that is utilized as a reference in the manuscript. Simulation results (spatio-temporal diagrams and NDR) obtained with the model are also presented.

\subsection*{High-Order Spectral method}

One-dimensional Nonlinear Schr\"odinger equation (Eq. (3) of manuscript) is established under the assumptions of weak nonlinearity of the wave field as well as the narrowbandedness of its energy content. Different degrees of nonlinearity ($\Gamma$ parameter) have been tested during the wave tank experiments reported in the manuscript. Then, it appears of great interest to have a high-fidelity numerical model able to replicate the experiments.

This is achieved thanks to the direct numerical simulation of the Euler's equation solved using the High-Order Spectral (HOS) method \cite{dommermuth1987high,west1987new}. The present numerical results have been obtained with the open-source solver HOS-NWT \cite{HOS-NWT}. This model solves the spatio-temporal evolution of the free surface elevation $\eta(z,t)$ in a so-called numerical wave tank. It reproduces all the physical features of the experimental (physical) facility including i) the generation of waves thanks to a wave maker and ii) the absorption of those waves when they reach the opposite wall. More details on the HOS-NWT numerical model as well as different validations performed against wave tank experiments can be found in \cite{bonnefoy2010time,ducrozet2012modified}.

This digital twin is utilized as a complement to the experiments. In the present study, the main advantage of the numerical solution is that it allows to overcome the limitation associated with the finite extent of the physical wave tank (and the associated measuring range $L_m=120~m$). As long as the spatial discretization is kept identical, the numerical wave tank can be of arbitrary size. Increasing the size results in a more accurate resolution $\Delta k_{min}$ for the estimation of wave number $k$, essential to the evaluation of the NDR and its properties (Fig.  4 of the manuscript). \\

\begin{figure*}[!t]
  \includegraphics[width=0.73\linewidth]{File_C12_C10_series_04_run_24_HOS_eps_0p04.png}
  \caption{\textbf{Numerical simulation using high-order spectral methods $\Gamma$ = 0.05,
 $\epsilon$ = 0.04.} \textbf{(top)}  Spatiotemporal diagram of the wave envelope amplitude $|\Psi(t,z)|$ evolution over the 531 m of the numerical water tank having 177 probes. \textbf{(bottom)} Nonlinear dispersion relation $|\widetilde{\psi}(\omega,k)|^2$ reconstructed from the evolution of the complex wave envelope.}
  \label{fig1si}
\end{figure*}

\subsection*{Numerical Simulations}

In a wave tank environment, numerical and physical experiments start from temporal initial conditions $\eta(z,t=0)$ at rest. Waves are generated thanks to a wave maker that imposes the spatial initial condition $\eta(z=0,t)$. The HOS digital twin of the \'Ecole Centrale de Nantes water tank facility (Fig. 2.a of original manuscript) uses the exact same wave maker motions than the experimental ones. This allows the direct deterministic comparison between the experimental and numerical wave probes measurements for validation purposes.

Numerical parameters are chosen after a careful convergence study that ensures the accuracy of the numerical solution. The time integration, achieved thanks to a step-adaptive Runge-Kutta method, is controlled by a tolerance parameter chosen as $10^{-8}$ for the present long-time integration. An HOS order of nonlinearity set to $5$ ensures an accurate numerical solution for all configurations tested. Regarding the spatial discretization, the total length of the numerical domain is set to $L_x=560~m$, discretized with $N_x=12288$ points/modes free of aliasing errors.\\

The wave conditions simulated cover a wide range of nonlinearity with $\Gamma=\left[ 0.05 ; 0.12 ; 0.33 ; 0.48 ; 0.65\right]$. The spatiotemporal diagram of the wave enveloped amplitude $|\Psi(t,z)|$ as well as the corresponding NDR $|\widetilde{\psi}(k,\omega)|^2$ are presented in Figures \ref{fig1si} to \ref{fig5si} with an increasing level of nonlinearity.

Taking into account the evolution of the nonlinear scale $z_{nlin}$ with $\Gamma$ (Eqs. (7) and (8) of manuscript), the analysis is performed on a length of $531~m$ for $\Gamma=0.05$ and $\Gamma=0.12$ in Figs. \ref{fig1si} and \ref{fig2si} respectively. With an increased nonlinearity $\Gamma=\left[ 0.33 ; 0.48 ; 0.65\right]$, the analysis is conducted on a shorter length of $378~m$. Similarly to the experimental facility, the numerical results use wave gauges that are equally spaced in the computational domain every $6~m$.\\

Overall, the HOS numerical simulations confirm observations in the physical wave tank. As expected, comparing left column of Fig. 1 in the manuscript and Fig. SI \ref{fig2si}, for small values of $\Gamma$ 1DNLSE and HOS results are very similar with a small departure from linear dispersion relation. Then, an increased nonlinearity is associated with coherent structures that are more visible in the $(t,z)$ diagram, see Fig. SI \ref{fig5si} for instance. However, compared to the 1DNLSE results (Fig.1 right column of manuscript), the corresponding straight lines are less frequent in the Euler's simulations. In addition, the evolution of those structures is clearly subjected to a negative speed when taking into account high-order nonlinear and dispersion effects. The corresponding frequency down-shift is also clearly visible in the NDR.

\begin{figure*}[!h]
  \includegraphics[width=0.73\linewidth]{File_C12_C10_series_04_run_10_HOS_eps_0p06.png}
  \caption{{Numerical simulation using high-order spectral methods $\Gamma$ = 0.12,
 $\epsilon$ = 0.06.} \textbf{(top)}  Spatiotemporal diagram of the wave envelope amplitude $|\Psi(t,z)|$ evolution over the 531 m of the numerical water tank having 177 probes. \textbf{(bottom)} Nonlinear dispersion relation $|\widetilde{\psi}(\omega,k)|^2$ reconstructed from the evolution of the complex wave envelope.}
  \label{fig2si}
\end{figure*}

\begin{figure*}[!h]
  \includegraphics[width=0.73\linewidth]{File_C12_C10_series_04_run_09_HOS_eps_0p1.png}
  \caption{{Numerical simulation using high-order spectral methods $\Gamma$ = 0.33,
 $\epsilon$ = 0.1.} \textbf{(top)}  Spatiotemporal diagram of the wave envelope amplitude $|\Psi(t,z)|$ evolution over the 378 m of the numerical water tank having 126 probes. \textbf{(bottom)} Nonlinear dispersion relation $|\widetilde{\psi}(\omega,k)|^2$ reconstructed from the evolution of the complex wave envelope.}
  \label{fig3si}
\end{figure*}

\begin{figure*}[!h]
  \includegraphics[width=0.73\linewidth]{File_C12_C10_series_04_run_22_HOS_eps_0p12.png}
  \caption{{Numerical simulation using high-order spectral methods $\Gamma$ = 0.48,
 $\epsilon$ = 0.12.} \textbf{(top)}  Spatiotemporal diagram of the wave envelope amplitude $|\Psi(t,z)|$ evolution over the 378 m of the numerical water tank having 126 probes. \textbf{(bottom)} Nonlinear dispersion relation $|\widetilde{\psi}(\omega,k)|^2$ reconstructed from the evolution of the complex wave envelope.}
  \label{fig4si}
\end{figure*}

\begin{figure*}[!h]
  \includegraphics[width=0.73\linewidth]{File_C12_C10_series_04_run_20_HOS_eps_0p14.png}
  \caption{{Numerical simulation using high-order spectral methods $\Gamma$ = 0.65,
 $\epsilon$ = 0.14.} \textbf{(top)}  Spatiotemporal diagram of the wave envelope amplitude $|\Psi(t,z)|$ evolution over the 378 m of the numerical water tank having 126 probes. \textbf{(bottom)} Nonlinear dispersion relation $|\widetilde{\psi}(\omega,k)|^2$ reconstructed from the evolution of the complex wave envelope.}
  \label{fig5si}
\end{figure*}

 \clearpage

\begin{acknowledgments}
This work has been partially supported  by the Agence Nationale de la
Recherche  through the LABEX CEMPI project (ANR-11-LABX-0007) and the ANR DYSTURB Project (ANR-17-CE30-0004),   the
Ministry of Higher Education and Research, Hauts de France council and
European Regional Development  Fund (ERDF) through the the Nord-Pas de
Calais Regional Research Council and the European Regional Development
Fund (ERDF) through the Contrat de Projets Etat-R\'egion (CPER
Photonics for Society P4S).  E.F. is thankful for partial support from the Simons Foundation/MPS No. 651463
\end{acknowledgments}