Dimensionality Reduction¶

In this exercise, we will use principal component analysis (PCA) to perform dimensionality reduction on a network traffic dataset. Recall that the netml library has mechanisms to extract a wide range of statistics from traffic flows. We will start with those (an N-dimensional feature set), and then use PCA to reduce the dimensionality of the data.

The notebook proceeds as follows:

  1. Load two packet captures (benign HTTP traffic and a Log4j scan) and convert them into flows.
  2. Extract per-flow statistical features and build a labeled dataset X, y.
  3. Apply PCA: project the data into two dimensions, examine explained variance, interpret the top components, and evaluate a classifier trained on the lower-dimensional data.
  4. (Optional) Visualize the data with t-SNE.
  5. (Optional) Train an autoencoder for dimensionality reduction and anomaly detection.
In [ ]:
import logging
logging.getLogger("scapy.runtime").setLevel(logging.ERROR)

from netml.pparser.parser import PCAP

import numpy as np
import pandas as pd

import matplotlib as mpl
import matplotlib.pyplot as plt
%matplotlib inline

Load the Packet Capture Files¶

We will use two packet captures from the data/ directory:

File Traffic Label
data/http.pcap Benign web (HTTP) traffic 0 ("http")
data/log4j.pcap Log4j vulnerability scan 1 ("log4j")

Each capture is loaded with netml, converted into flows, and then summarized with the STATS feature set (one row of statistics per flow). The label tells us which capture each flow came from, which we will use later to color-code plots and to evaluate classifiers trained on the reduced-dimension data.

In [ ]:
hpcap = PCAP('data/http.pcap', flow_ptks_thres=2, verbose=10)
lpcap = PCAP('data/log4j.pcap', flow_ptks_thres=2, verbose=10)

Convert the Packet Capture Into Flows¶

The pcap2flows function in netml converts each pcap file into a list of flows. Examine the resulting data structure (hpcap.flows, lpcap.flows). What does it contain?

In [ ]:
# extract flows from each pcap
hpcap.pcap2flows()
lpcap.pcap2flows()

Explore the Flows¶

How many flows are in each of these pcaps? (Use the netml library output to determine the size of each data structure.)

In [ ]:
# TODO: print the number of flows in the HTTP capture and in the Log4j capture

Extract Features and Build the Dataset¶

Use the STATS feature set in netml (flow2features('STATS', ...)) to compute twelve statistics for each flow. The features are, in order: flow duration, packets per second, bytes per second, and the mean, standard deviation, 25th percentile, median, 75th percentile, minimum, and maximum of the packet sizes, followed by the number of packets and the number of bytes in the flow.

The code below extracts the features from both captures, attaches the label (0 for HTTP, 1 for Log4j), and concatenates the two into a single dataset. It then defines:

  • X: the feature matrix (one row per flow, twelve columns), and
  • y: the label for each row.

These two variables are used throughout the rest of the notebook.

In [ ]:
# extract features from each flow via STATS
hpcap.flow2features('STATS', fft=False, header=False)
hd = pd.DataFrame(hpcap.features)

lpcap.flow2features('STATS', fft=False, header=False)
ld = pd.DataFrame(lpcap.features)

# human-readable names for the twelve STATS features
labels = ['duration', 'packets per second',
          'bytes per second',
          'mean', 'standard deviation',
          '25', 'median', '75',
          'minimum', 'maximum',
          'packets', 'bytes']

# attach labels: 0 = http (benign), 1 = log4j (scan)
pd.set_option('mode.chained_assignment', None)
hd['label'] = 0
ld['label'] = 1

data = pd.concat([hd, ld], ignore_index=True)

X = data.loc[:, :11]     # features
y = data['label']        # labels

data

Dimensionality Reduction¶

One way to reduce the dimensionality of a dataset is with a technique called principal components analysis (PCA). PCA is a linear transformation that maps the points into a space where the lower dimensions are orthogonal and also capture the highest variance in the dataset.

So, the first principal component (PC1) is a vector that is oriented in the direction that captures the highest variance in the dataset. You can think of this as the single dimension that has the most information in the dataset. PC2 captures the next highest variance, in the direction that is orthogonal to PC1, and so forth.

There are many applications of PCA, but one application is the visualization of a high-dimension dataset, since PCA is just a transformation of the data that does not inherently lose information. When we only project into the top two dimensions, some information is lost (whatever is in the lower principal components), but we can visualize the data in terms of the two dimensions that capture the most "information" (i.e., variance) in the dataset.

Project the Data Onto the Top Two Principal Components¶

Use sklearn.decomposition.PCA with n_components=2 to fit and transform X. Store the result in X_2D and make a scatter plot of the two components, colored by the label y (HTTP vs. Log4j). Does the two-dimensional projection separate the two kinds of traffic?

In [ ]:
from sklearn.decomposition import PCA

colors = ['red', 'blue']   # 0 = http, 1 = log4j

# TODO: fit a PCA with two components on X and transform X into X_2D
pca = PCA(n_components=2)
X_2D = None

# TODO: scatter plot of X_2D, colored by y
#   hint: plt.scatter(*zip(*X_2D), c=y, cmap=mpl.colors.ListedColormap(colors), alpha=0.5)
plt.xlabel('PC 1')
plt.ylabel('PC 2')
plt.show()

Understanding Explained Variance¶

The first principal components will capture most of the variance. We can look at the explained variance ratio to understand how much variance each of those principal components captures, which will give us a sense of where we might be able to cut off this data, particularly if we look at cumulative explained variance.

In [ ]:
pca.explained_variance_ratio_

Below, we fit a PCA with all twelve components and plot the cumulative explained variance (a "scree" plot) alongside the per-component explained variance ratio. Note that the component count on the x-axis starts at 1: the first point/bar is PC1, not "zero components".

In [ ]:
pca = PCA(n_components=12)
pca.fit(X)

ratio = pca.explained_variance_ratio_
components = range(1, len(ratio) + 1)   # components are numbered starting at 1

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4))

# cumulative explained variance (scree plot)
ax1.plot(components, np.cumsum(ratio), marker='o')
ax1.set_ylim([0.9, 1])
ax1.set_xticks(list(components))
ax1.set_xlabel('number of components')
ax1.set_ylabel('cumulative explained variance')

# per-component explained variance ratio
ax2.bar(components, ratio)
ax2.set_xticks(list(components))
ax2.set_xlabel('principal component')
ax2.set_ylabel('explained variance ratio')

plt.tight_layout()
plt.show()

Understanding the Top Principal Components¶

If we look at the top principal component, we can also see that it is a linear combination of our original features, which has a disproportionate "amount" of the 25th percentile of packet size.

In [ ]:
pca.components_[0]
In [ ]:
# Sort in descending order
indices = np.argsort(pca.components_[0])[::-1]

# Sort the labels in a corresponding fashion
names = [labels[i] for i in indices]
names

Evaluate the Lower Dimension Model¶

Train a classifier (for example, a RandomForestClassifier) and evaluate it with 5-fold cross-validation, first on the full twelve-feature dataset X and then on the two-dimensional projection X_2D. What happens to the accuracy? Why?

What might be a good reason to reduce the dimensionality of the dataset?

In [ ]:
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import KFold, cross_val_score

kf = KFold(n_splits=5, shuffle=True, random_state=23)

rf = RandomForestClassifier(random_state=0,
                            n_jobs=-1,
                            n_estimators=100,
                            class_weight='balanced')

# TODO: cross-validated accuracy on the full feature set X
#   hint: cross_val_score(rf, X, y, cv=kf, scoring='accuracy').mean()

# TODO: cross-validated accuracy on the two-dimensional projection X_2D

Optional: t-SNE for Visualization¶

t-distributed Stochastic Neighbor Embedding (t-SNE) is a nonlinear dimensionality reduction technique that tries to place points that are close together in the original high-dimensional space close together in two dimensions. Unlike PCA, it is stochastic (different random seeds give different pictures), it is intended only for visualization, and it has no transform method: it cannot map new data points into an existing embedding, so it cannot be used as a preprocessing step for a classifier.

Standardize X, subsample to at most 2,000 rows (t-SNE is slow on large datasets), fit TSNE(n_components=2, perplexity=30, random_state=0), and plot the result colored by y. How does the picture compare to the PCA projection above?

In [ ]:
from sklearn.manifold import TSNE
from sklearn.preprocessing import StandardScaler

X_scaled = StandardScaler().fit_transform(X)

# subsample to at most 2000 rows
rng = np.random.default_rng(0)
n = min(2000, X_scaled.shape[0])
idx = rng.choice(X_scaled.shape[0], size=n, replace=False)
X_sub = X_scaled[idx]
y_sub = np.asarray(y)[idx]

# TODO: fit t-SNE on X_sub (n_components=2, perplexity=30, random_state=0)
tsne = None
X_tsne = None

# TODO: scatter plot of X_tsne, colored by y_sub (0 = http, 1 = log4j)
plt.xlabel('t-SNE 1')
plt.ylabel('t-SNE 2')
plt.show()

Optional: Autoencoder for Dimensionality Reduction and Anomaly Detection¶

An autoencoder is a neural network trained to reproduce its input: an encoder compresses the input down to a small "bottleneck" representation (here, two dimensions, analogous to the top two principal components), and a decoder reconstructs the original features from it. Because the bottleneck is much smaller than the input, the network must learn the structure of the data it is trained on.

If we train the autoencoder only on benign HTTP flows, it learns to reconstruct HTTP traffic well. Log4j scan flows have different statistics (for example, packet sizes and rates), which the network has never seen, so the decoder reconstructs them poorly and the reconstruction error is large. Reconstruction error therefore acts as an anomaly score, and we never needed the attack labels to train the model.

Build a small dense autoencoder with keras (input dimension = number of features, bottleneck of 2), train it on the HTTP rows only, then compute the per-row reconstruction error on all rows and plot the error distribution for HTTP vs. Log4j flows. This section requires TensorFlow/Keras; the cell below prints a message and skips if it is not installed.

In [ ]:
from sklearn.preprocessing import StandardScaler

try:
    import keras
    from keras import layers
except Exception:
    keras = None
    print("keras/tensorflow is not installed. Run `pip install tensorflow` to run this optional section.")

if keras is not None:
    n_features = X.shape[1]

    # standardize using statistics from the benign (HTTP) rows only
    X_http = X[y == 0]
    scaler = StandardScaler().fit(X_http)
    X_http_scaled = scaler.transform(X_http)
    X_all_scaled = scaler.transform(X)

    # TODO: build a small dense autoencoder, e.g. n_features -> 8 -> 2 (bottleneck) -> 8 -> n_features
    autoencoder = keras.Sequential([
        layers.Input(shape=(n_features,)),
        # TODO: encoder layers (ending in a Dense(2) bottleneck)
        # TODO: decoder layers (ending in Dense(n_features))
    ])
    autoencoder.compile(optimizer='adam', loss='mse')

    # TODO: train on the HTTP rows only (the input is also the target)
    #   hint: autoencoder.fit(X_http_scaled, X_http_scaled, epochs=50, batch_size=32, verbose=0)

    # TODO: reconstruct all rows and compute the per-row mean squared reconstruction error
    #   hint: reconstructed = autoencoder.predict(X_all_scaled, verbose=0)
    #         error = np.mean((X_all_scaled - reconstructed) ** 2, axis=1)
    error = None

    # TODO: plot histograms of the reconstruction error for http (y == 0) and log4j (y == 1)
    #   hint: plt.hist(error[np.asarray(y) == 0], bins=50, alpha=0.5, label='http')
    plt.xlabel('reconstruction error')
    plt.ylabel('number of flows')
    plt.legend()
    plt.show()