<jupyter_start><jupyter_text>Autoencoders for Anomaly Detection  High Performance Computing**High Performance Computing**HPC refers to HW/SW infrastructures for particularly intensive workloads  High Performance Computing**HPC is (somewhat) distinct from cloud computing*** Cloud computing is mostly about running (and scaling services)* ...HPC is all about _performance_Typical applications: simulation, massive data analysis, training large ML models**HPC systems follow a batch computation paradigm*** Users send _jobs_ to the systems (i.e. configuration for running a program)* Jobs end in one of several _queues_* A _job scheduler_ draws from the queue* ...And dispatches jobs to computational _nodes_ for execution  High Performance Computing**HPC systems can be large and complex**E.g. Marconi-100 at CINECA, which was the 9-th most powerful supercompuer [(as of June 2020)](https://www.top500.org/lists/top500/2020/06/)* The system has 31,360 cores overall!**Configuring (and maintaining the configuration) of these systems*** ...Is of paramount important, as it has an impact on the performance* ...Is very challenging, due to their _large scale_ and the presence of _node heterogeneity_ Hence the interest in _detecting anomalous conditions_  The Dataset**As an example, we will consider the DAVIDE system**Small scale, energy-aware architecture:* Top of the line components (at the time), liquid cooled* An advanced monitoring and control infrastructure (ExaMon)* ...Developed together with UniBoThe system went out of production in January 2020**The monitoring system enables anomaly detection*** Data is collected from a number of samples with high-frequency* Long term storage only for averages over 5 minute intervals* Anomalies correspond to unwanted configurations of the _frequency governor_* ...Which can throttle performance to save power or prevent overheating  A Look at the Dataset**Our dataset refers to the non-idle periods of a single node**<jupyter_code>print(f'#examples: {hpc.shape[0]}, #columns: {hpc.shape[1]}')
hpc.iloc[:3]<jupyter_output>#examples: 6667, #columns: 161<jupyter_text>* This still a time series, but a _multivariate_ one  A Look at the Dataset**How to display a multivariate series? Approach 1: showing _individual columns_**<jupyter_code>tmp = pd.Series(index=hpc['timestamp'], data=hpc[inputs[0]].values)
nn.plot_series(tmp, figsize=figsize)<jupyter_output><empty_output><jupyter_text>* The series contains significant gaps (i.e. the idle periods)  A Look at the Dataset**Approach 2: obtaining _statistics_**<jupyter_code>hpc[inputs].describe()<jupyter_output><empty_output><jupyter_text>* No missing value, _normalized_ data  A Look at the Dataset**Approach 3: standardize, then use a _heatmap_**<jupyter_code>hpcsv = hpc.copy()
hpcsv[inputs] = (hpcsv[inputs] - hpcsv[inputs].mean()) / hpcsv[inputs].std()
nn.plot_dataframe(hpcsv[inputs], figsize=figsize)<jupyter_output><empty_output><jupyter_text>* White = mean, red = below mean, blue = above mean  Anomalies**There are three possible configurations of the frequency governor:*** Mode 0 or "normal": frequency proportional to the workload* Mode 1 or "power saving": frequency always at the minimum value* Mode 2 or "performance": frequency always at the maximum value**On this dataset, this information is known**...And it will serve as our ground truth* We will focus on discriminating normal from non-normal behavior* I.e. we will treat both "power saving" and "performance" configurations as anomalous**Detecting them will be _challenging_*** Since the signals vary so much when the running job changes  Anomalies**We can plot the location of the anomalies:**<jupyter_code>labels = pd.Series(index=hpcsv.index, data=(hpcsv['anomaly'] != 0), dtype=int)
nn.plot_dataframe(hpcsv[inputs], labels, figsize=figsize)<jupyter_output><empty_output><jupyter_text>* On the top, blue = normal, orange = anomaly  A KDE Approach   KDE Approach**Let's try first a density estimation approach (once again using KDE)**First, we need to _standardize_ the data again, based on _training information alone_<jupyter_code>tr_end, val_end = 3000, 4500

hpcs = hpc.copy()
tmp = hpcs.iloc[:tr_end]
hpcs[inputs] = (hpcs[inputs] - tmp[inputs].mean()) / tmp[inputs].std()<jupyter_output><empty_output><jupyter_text>* This is needed so that we do not accidentally exploit test set information* The training set separator was chosen so as not to include anomalies**Then we can separate training, validation, and test data:**<jupyter_code>trdata = hpcs.iloc[:tr_end]
valdata = hpcs.iloc[tr_end:val_end]
tsdata = hpcs.iloc[val_end:]<jupyter_output><empty_output><jupyter_text>A KDE Approach**Then we estimate the optimal bandwidth:**<jupyter_code>params = {'bandwidth': np.linspace(0.1, 1, 10)}
opt = GridSearchCV(KernelDensity(kernel='gaussian'), params, cv=5)
opt.fit(trdata[inputs])
opt.best_params_<jupyter_output><empty_output><jupyter_text>**...Ad we can train the estimator and generate the anomaly signal:**<jupyter_code>h = opt.best_params_['bandwidth']
kde = KernelDensity(bandwidth=h)
kde.fit(trdata[inputs])
ldens = kde.score_samples(hpcs[inputs])
signal_kde = pd.Series(index=hpcs.index, data=-ldens)<jupyter_output><empty_output><jupyter_text>* Tuning the bandwidth and obtaining densities are relatively expensive operations  KDE Approach**There is _a good match_ with the anomalies, but also many _spurious peaks_**<jupyter_code>nn.plot_signal(signal_kde, labels, figsize=figsize)<jupyter_output><empty_output><jupyter_text>* This is mostly due to the large variations due to job changes  KDE Approach**We then need to define the threshold, but for that we need a cost model**Our main goal is to detect anomalies, not anticipating them* Misconfigurations in HPC are usually not critical * ...And cause little issue, unless they stay unchecked for very long**We will use a simple _cost model_:*** $c_{\mathit{alarm}}$ for false positive (erroneous detections)* $c_{\mathit{missed}}$ for false negatives (undetected anomalies)* Detections are fine as long as they are within $\mathit{tolerance}$ units from the anomaly<jupyter_code>c_alarm, c_missed, tolerance = 1, 5, 12
cmodel =  nn.HPCMetrics(c_alarm, c_missed, tolerance)<jupyter_output><empty_output><jupyter_text>The implementation details can be found in the `nn` utility module  KDE Approach**We can now optimize the threshold over the validation set**<jupyter_code>th_range = np.linspace(1e4, 1e5, 100)
th_kde, val_cost_kde = nn.opt_threshold(signal_kde[tr_end:val_end],
                                        valdata['anomaly'],
                                        th_range, cmodel)
print(f'Best threshold: {th_kde:.3f}')
tr_cost_kde = cmodel.cost(signal_kde[:tr_end], hpcs['anomaly'][:tr_end], th_kde)
print(f'Cost on the training set: {tr_cost_kde}')
print(f'Cost on the validation set: {val_cost_kde}')
ts_cost_kde = cmodel.cost(signal_kde[val_end:], hpcs['anomaly'][val_end:], th_kde)
print(f'Cost on the test set: {ts_cost_kde}')<jupyter_output>Best threshold: 82727.273
Cost on the training set: 0
Cost on the validation set: 266
Cost on the test set: 265<jupyter_text>* The `opt_threshold` function runs the usual line search process* In this case the training and validation set are completely separated  The Trouble with KDE**KDE-based approach works well, but have some _issues_**First, KDE itself runs into trouble with high-dimensional data:* With a larger dimensionality, _prediction times_ grows...* ...And _more data_ is needed to obtain reliable results* We are not seeing that too much, but eventually it will become a problem Second, KDE gives you _nothing more_ than an anomaly signal* Determining the cause of the anomaly is up to a domain expert* ...Who needs to take a look at all the supposed anomalous data* This is doable in low-dimensional spaces, but _way harder on high-dimensional one_  Autoencoders for Anomaly Detection  Autoencoders**An autoencoder is _a type of neural network_**The network is designed to _reconstruct its input vector_* The input is some tensor $x$ and the output _should be_ the same tensor $x$**Autoencoders can be broken down in two halves*** An encoding part, i.e. $\mathit{encode}(x, \theta_e)$, mapping $x$ into a vector of _latent variables $z$_* A decoding part, i.e. $\mathit{decode}(z, \theta_d)$, mapping $z$ into reconstructed input tensor**Autoencoders are trained so as to satisfy:**$$\mathit{decode}(\mathit{encode}(\hat{x}_i, \theta_e), \theta_d) \simeq \hat{x}_i$$* I.e. $\mathit{decode}$, when applied to the output of $\mathit{encode}$* ...Should approximately return the input vector itselfA nice introduction and tutorial about autoencoders can be found [on the Keras blog](https://blog.keras.io/building-autoencoders-in-keras.html)  Autoencoders**Formally, we typically employ an MSE loss**$$L(\theta_e, \theta_d) = \sum_{i=1}^n \left\|\hat{x}_i - \mathit{decode}(\mathit{encode}(\hat{x}_i, \theta_e), \theta_d)\right\|_2^2$$* This is trivial to satisfy if both $\mathit{encode}$ and $\mathit{decode}$ learn an identity relation* ...So we need to prevent that**There are _two main approaches_ to avoid learning a trivial mapping*** Using an _information bottleneck_, i.e. making sure that $z$ has fewer dimensions that $x$* Use a regularization to enforce _sparse encodings_, e.g.:$$L(\theta_e, \theta_d) = \sum_{i=1}^n \left\|\hat{x}_i - \mathit{decode}(\mathit{encode}(\hat{x}_i, \theta_e), \theta_d)\right\|_2^2 + \alpha \| \mathit{encode}(x, \theta_e) \|_1$$  Autoencoders for Anomaly Detection**Autoencoders can be used for anomaly detection**...By using the _reconstruction error as an anomaly signal_, e.g.:$$\left\|x - \mathit{decode}(\mathit{encode}(x, \theta_e), \theta_d)\right\|_2^2 > \theta$$**This approach has some PROs and CONs:*** Compared to KDE  - Neural Networks have good _support for high dimensional data_  - ...Plus _limited overfitting_ and _fast prediction/detection time_  - However, error reconstruction can be _harder than density estimation_* Compared to autoregressors  - Reconstructing an input is _easier than predicting the future_  - ...So, we tend to get higher reliability  Autoencoders in Keras**Let's build an autoencoder in practice (with tensorflow 2.0 and keras)**First, we build the model using (e.g.) the functional API<jupyter_code>from tensorflow import keras
from tensorflow.keras import layers, callbacks

input_shape = (len(inputs), )
ae_x = keras.Input(shape=input_shape, dtype='float32')
ae_z = layers.Dense(64, activation='relu')(ae_x)
ae_y = layers.Dense(len(inputs), activation='linear')(ae_z)
ae = keras.Model(ae_x, ae_y)<jupyter_output><empty_output><jupyter_text>* `Input` builds the entry point for the input data* `Dense` builds a fully connected layer* "Calling" layer A with parameter B attaches B to A* `Model` builds a model object with the specified input and output  Autoencoders in Keras**Then we can prepare our model for training**In keras terms, we _compile_ it:<jupyter_code>ae.compile(optimizer='RMSProp', loss='mse')<jupyter_output><empty_output><jupyter_text>* We are using the `RMSProp` optimizer (a variant of Stochastic Gradient Descent)**Then we can start training:**<jupyter_code>cb = [callbacks.EarlyStopping(patience=3, restore_best_weights=True)]
history = ae.fit(trdata[inputs], trdata[inputs], validation_split=0.1,
                 callbacks=cb,
                 batch_size=32, epochs=20, verbose=0)<jupyter_output><empty_output><jupyter_text>* We are using a callback to stop training early* ...If no improvement on the validation set is observed for 3 epochs  Autoencoders in Keras**Let's have a look at the loss evolution over different epochs**<jupyter_code>nn.plot_training_history(history, figsize=figsize)<jupyter_output><empty_output><jupyter_text>Autoencoders in Keras**Finally, we can obtain the predictions**<jupyter_code>preds = pd.DataFrame(index=hpcs.index, columns=inputs, data=ae.predict(hpcs[inputs]))
preds.head()<jupyter_output><empty_output><jupyter_text>Alarm Signal**We can finally obtain our alarm signal, i.e. the sum of squared errors**<jupyter_code>sse = np.sum(np.square(preds - hpcs[inputs]), axis=1)
signal_ae = pd.Series(index=hpcs.index, data=sse)
nn.plot_signal(signal_ae, labels, figsize=figsize)<jupyter_output><empty_output><jupyter_text>* It is actually quite similar to the KDE signal  Threshold Optimization**Then we can optimize the threshold as usual**<jupyter_code>th_ae, val_cost_ae = nn.opt_threshold(signal_ae[tr_end:val_end],
                                        hpcs['anomaly'][tr_end:val_end],
                                        th_range, cmodel)
print(f'Best threshold: {th_ae:.3f}')
tr_cost_ae = cmodel.cost(signal_ae[:tr_end], hpcs['anomaly'][:tr_end], th_ae)
print(f'Cost on the training set: {tr_cost_ae}')
print(f'Cost on the validation set: {val_cost_ae}')
ts_cost_ae = cmodel.cost(signal_ae[val_end:], hpcs['anomaly'][val_end:], th_ae)
print(f'Cost on the test set: {ts_cost_ae}')<jupyter_output>Best threshold: 93636.364
Cost on the training set: 0
Cost on the validation set: 264
Cost on the test set: 265<jupyter_text>* We have more or less the same performance as KDE  Mutiple Signal Analysis**But autoencoders do _more than just anomaly detection_!*** Instead of having a single signal we have _many_* So we can look at the _individual_ reconstruction errors<jupyter_code>se = np.square(preds - hpcs[inputs])
signals_ae = pd.DataFrame(index=hpcs.index, columns=inputs, data=se)
nn.plot_dataframe(signals_ae, labels, vmin=-5e4, vmax=5e4, figsize=figsize)<jupyter_output><empty_output><jupyter_text>Mutiple Signal Analysis**Reconstruction errors are often concentrated on a few signals*** These correspond to the properties of the input vector that were harder to reconstruct* ...And often they are useful clues about the _nature of the anomaly_<jupyter_code>se = np.square(preds - hpcs[inputs])
signals_ae = pd.DataFrame(index=hpcs.index, columns=inputs, data=se)
nn.plot_dataframe(signals_ae, labels, vmin=-5e4, vmax=5e4, figsize=figsize)<jupyter_output><empty_output><jupyter_text>Multiple Signal Analysis**Let's focus on the last _mode 1_ anomaly ("power saving" mode)**Here are the 8 largest errors in descending order<jupyter_code>last_mode_1 = hpcs.index[hpcs['anomaly']==1][-1]
se.iloc[last_mode_1].sort_values(ascending=False)[:8]<jupyter_output><empty_output><jupyter_text>* They are mostly related to performance (e.g. "ips" - Instructions Per Second)* ...As it should be!  Multiple Signal Analysis**Now, let's move to the last _mode 2_ anomaly ("performance" mode)**Here are the 8 largest errors in descending order<jupyter_code>last_mode_2 = hpcs.index[hpcs['anomaly']==2][-1]
se.iloc[last_mode_2].sort_values(ascending=False)[:8]<jupyter_output><empty_output><jupyter_text>* Again, they are performance related  Multiple Signal Analysis**Here are the _average errors_ for mode 1 anomalies**<jupyter_code>mode_1 = hpcs.index[hpcs['anomaly']==1]
tmp =  se.iloc[mode_1].mean().sort_values(ascending=False)
nn.plot_bars(tmp, tick_gap=20, figsize=figsize)<jupyter_output><empty_output><jupyter_text>* Errors are concentrated on 10-20 features  Multiple Signal Analysis**These are the 20 _largest_ average errors for _mode 1_ anomalies**<jupyter_code>mode_1 = hpcs.index[hpcs['anomaly']==1]
tmp =  se.iloc[mode_1].mean().sort_values(ascending=False)
nn.plot_bars(tmp.iloc[:20], figsize=figsize)<jupyter_output><empty_output><jupyter_text>* The largest errors are on "ips", then on "util" (utilization)  Multiple Signal Analysis**Let's repeat the analysis for _mode 2_. Here are the _average errors_**<jupyter_code>mode_2 = hpcs.index[hpcs['anomaly']==2]
tmp =  se.iloc[mode_2].mean().sort_values(ascending=False)
nn.plot_bars(tmp, tick_gap=20, figsize=figsize)<jupyter_output><empty_output><jupyter_text>* The situation is similar to mode 1  Multiple Signal Analysis**The 20 _largest_ average errors for _mode 2_**<jupyter_code>mode_2 = hpcs.index[hpcs['anomaly']==2]
tmp =  se.iloc[mode_2].mean().sort_values(ascending=False)
nn.plot_bars(tmp.iloc[:20], figsize=figsize)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10805097917314269, 'fraction_numerical': 0.011501398818775257, 'mean_word_length': 6.09572121746802, 'pattern_counts': {'":': 3, '<': 106, '<?xml version=': 0, '>': 107, 'https://': 2, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '27', 'dir': 'jupyter-structured-clean-dedup', 'id': '82789', 'license': 'MIT', 'n_tokens_mistral': 5257, 'n_tokens_neox': 4947, 'n_words': 1655, 'path': 'notebooks/1. Autoencoders for Anomaly Detection.ipynb', 'repo_name': 'lompabo/aiiti-course-2021-04', 'size': '16085'}