Model Families

Recall that we expressed our fraud dataset in terms of seven numerical features: four features from the raw data and three from the one-hot encoding of the payment method. We can thus think of each transaction as a point in a seven-dimensional real vector space, and our goal is to divide up the space into areas of fraud and nonfraud transactions. The “model” output by our machine learning algorithm is a description of this division of the vector space.

In theory, the division of our vector space into fraud and nonfraud areas can be infinitely complex; in practice, most algorithms produce a decision boundary, which is a surface in the vector space.⁷ One side of the decision boundary consists of the points labeled as fraud, and the other side consists of the points labeled as nonfraud. The boundary can be as simple as a line (or hyperplane in higher dimensions) or as complex as a union of nonlinear disconnected regions. Figure 2-1 presents some examples.

Figure 2-1. Examples of two-dimensional spaces divided by a decision boundary

If we want to be more granular, instead of mapping each point in the vector space to either “fraud” or “nonfraud,” we can map each point to a probability of fraud. In this case our machine learning algorithm outputs a function that assigns each point in the vector space a value between 0 and 1, to be interpreted in our example as the probability of fraud.

Any given machine learning algorithm restricts itself to finding a certain type of decision boundary or probability function that can be described by a finite number of model parameters. The simplest decision boundary is a linear decision boundary—that is, a hyperplane in the vector space. An oriented hyperplane H in an n-dimensional vector space can be described by an n-dimensional vector $\mathit{\theta }$ orthogonal to the hyperplane, plus another vector $\mathit{\beta }$ indicating how far the hyperplane is from the origin:

This description allows us to divide the vector space in two; to assign probabilities we want to look at the distance of the point $\mathit{x}$ from the hyperplane H. We can thus compute a real-valued “score”:

where we have let $b=-\mathit{\theta }·\mathit{\beta }$. Our model to compute the score can thus be described by n + 1 model parameters: n parameters to describe the vector $\mathit{\theta }$, and one “offset” parameter b. To turn the score into a classification, we simply choose a threshold t above which all scores indicate fraud, and below which all scores indicate nonfraud.

If we want to map the real-valued score $s\left(\mathit{x}\right)$ to a probability, we must apply a function that maps the real numbers to the interval [0,1]. The standard function to apply is known as the logistic function or sigmoid function,⁸ as illustrated in Figure 2-2. It is formulated as:

Figure 2-2. The sigmoid function

The output of the logistic function can be interpreted as a probability, allowing us to define the likelihood of the dependent variable taking on a particular value given some input feature vector $\mathit{x}$.

Loss Functions

Now that we have restricted our choice of prediction algorithms to a certain parametrized family, we must choose the best one for the given training data. How do we know when we have found the best algorithm? We define the best algorithm to be one that optimizes some quantity computed from the data. This quantity is called an objective function. In the case of machine learning, the objective function is also known as a cost function or loss function, because it measures the “cost” of wrong predictions or the “loss” associated with them.

Mathematically, a loss function is a function that maps a set of pairs of (predicted label, truth label) to a real number. The goal of a machine learning algorithm is to find the model parameters that produce predicted labels for the training set that minimize the loss function.

In regression problems, for which the prediction algorithm outputs a real number instead of a label, the standard loss function is the sum of squared errors. If $y_i$ is the true value and ${\widehat{y}}_i$ is the predicted value, the loss function is as follows:

We can use this loss function for classification problems as well, where $y_i$ is either 0 or 1, and ${\widehat{y}}_i$ is the probability estimate output by the algorithm.

For logistic regression, we use negative log likelihood as the loss function. The likelihood of a set of probability predictions $\left\{p_i\right\}$ for a given set of ground truth labels $\left\{y_i\right\}$ is defined to be the probability that these truth labels would have arisen if sampled from a set of binomial distributions according to the probabilities $\left\{p_i\right\}$. Concretely, if the truth label $y_i$ is 0, the likelihood of 0 is $1-p_i$—the probability that 0 would have been sampled from a binomial distribution with mean $p_i$. If the truth label $y_i$ is 1, the likelihood is $p_i$.

The likelihood of the entire set of predictions is the product of the individual likelihoods:

The goal of logistic regression is to find parameters that produce probabilities $\left\{p_i\right\}$ that maximize the likelihood.

To make computations easier, and, in particular, because most optimization methods require computing derivatives of the loss function, we take the negative log of the likelihood. As maximizing the likelihood is equivalent to minimizing the negative log likelihood, we call negative log likelihood the loss function:

Here we have used the fact that $y_i$ is always 0 or 1 to combine the two products into a single term.

Optimization

The last step in the machine learning procedure is to search for the optimal set of parameters that minimizes the loss function. To carry out this search we use an optimization algorithm. There may be many different optimization algorithms available to you when fitting your machine learning model.⁹ Most scikit-learn estimators (e.g., LogisticRegression) allow you to specify the numerical solver to use, but what are the differences between the different options, and how do you go about selecting one?

The job of an optimization algorithm is to minimize (or maximize) an objective function. In the case of machine learning, the objective function is expressed in terms of the model’s learnable parameters ($\mathit{\theta }$ and b in the previous example), and the goal is to find the values of $\mathit{\theta }$ and b that optimize the objective function.

Optimization algorithms mainly come in two different flavors:

First-order algorithms

These algorithms optimize the objective function using the first derivatives of the function with respect to the learnable parameters. Gradient descent methods are the most popular types of first-order optimization algorithms; we can use them to find the inputs to a function that give the minimum (or maximum) value. Computing the gradient of a function (i.e., the partial derivatives with respect to each variable) allows us to determine the instantaneous direction that the parameters need to move in order to achieve a more optimal outcome.

Second-order algorithms

As the name suggests, these algorithms use the second derivatives to optimize the objective function. Second-order algorithms will not fall victim to paths of slow convergence. For example, second-order algorithms are good at detecting saddle points, whereas first-order algorithms are likely to become stuck at these points. However, second-order methods are often slower and more expensive to compute.

First-order methods tend to be much more frequently used because of their relative efficiency. Picking a suitable optimization algorithm depends on the size of the dataset, the nature of the cost function, the type of learning problem, and speed/resource requirements for the operation. In addition, some regularization techniques can also have compatibility issues with certain types of optimizers. First-order algorithms include the following:

• LIBLINEAR¹⁰ is the default solver for the linear estimators in scikit-learn. This algorithm tends to not do well on larger datasets; as suggested by the scikit-learn documentation, the Stochastic Average Gradient (SAG) or SAGA (improvement to SAG) methods work better for large datasets.¹¹

  Different optimization algorithms deal with multiclass classification differently. LIBLINEAR works only on binary classification. For it to work in a multiclass scenario, it has to use the one-versus-rest scheme; we discuss this scheme more fully in Chapter 5.

• Stochastic Gradient Descent (SGD) is a very simple and efficient algorithm for optimization that performs a parameter update for each separate training example. The stochastic nature of the gradient descent means that the algorithm is more likely to discover new and possibly better local minima as compared to standard gradient descent. However, it typically results in high-variance oscillations, which can result in a delay in convergence. This can be solved with a decreasing learning rate (i.e., exponentially decrease the learning rate) that results in smaller fluctuations as the algorithm approaches convergence.

  The technique called momentum also helps accelerate SGD convergence by navigating the optimization movement only in the relevant directions and softening any movement in irrelevant directions, which stabilizes SGD.

• Optimization algorithms such as AdaGrad, AdaDelta, and Adam (Adaptive Moment Estimation) allow for separate and adaptive learning rates for each parameter that solve some problems in the other simpler gradient descent algorithms.

• When your training dataset is large, you will need to use a distributed optimization algorithm. One popular algorithm is the Alternating Direction Method of Multipliers (ADMM).¹²

Example: Gradient descent

To conclude this section, we go briefly into the details of gradient descent, a powerful optimization algorithm that has been applied to many different machine learning problems.

The standard algorithm for gradient descent is as follows:

1. Select random starting parameters for the machine learning model. In the case of a linear model, this means selecting a random normal vector $\mathit{\theta }$ and offset $\mathit{\beta }$, which results in a random hyperplane in n-dimensional space.

2. Compute the value of the gradient of the loss function for this model at the point described by these parameters.

3. Change the model parameters in the direction of greatest gradient decrease by a certain small magnitude, typically referred to as the α or learning rate.

4. Iterate: repeat steps 2 and 3 until convergence or a satisfactory optimization result is attained.

Figure 2-3 illustrates the intermediate results of a gradient descent optimization process of a linear regression. At zero iterations, observe that the regression line, formed with the randomly chosen parameters, does not fit the dataset at all. As you can imagine, the value of the sum-of-squares cost function is quite large at this point. At three iterations, notice that the regression line has very quickly moved to a more sensible position. Between 5 and 20 iterations the regression line slowly adjusts itself to more optimal positions where the cost function is minimized. If performing any more iterations doesn’t decrease the cost function significantly, we can say that the optimization has converged, and we have the final learned parameters for the trained model.

Figure 2-3. Regression line after a progressive number of iterations of gradient descent optimization (0, 3, 5, and 20 iterations of gradient descent, shown at the upper left, upper right, lower left, and lower right, respectively)

Logistic Regression

Although we discussed logistic regression in some detail earlier, we go over its key properties here. Logistic regression takes as input numerical feature vectors and attempts to predict the log odds¹³ of each data point occurring; we can convert the log odds to probabilities by using the sigmoid function discussed earlier. In the log odds space, the decision boundary is linear, so increasing the value of a feature monotonically increases or decreases (depending on the sign of the coefficient) the score output by the model.

Logistic regression is one of the most popular algorithms in practice due to a number of properties: it can be trained very efficiently and in a distributed manner, it scales well to millions of features, it admits a concise description and fast scoring algorithm (a simple dot product), and it is explainable—each feature’s contribution to the final score can be computed.

However, there are a few important things to be aware of when considering logistic regression as a supervised learning technique:

• Logistic regression assumes linearity of features (independent variables) and log odds, requiring that features are linearly related to the log odds. If this assumption is broken the model will perform poorly.

• Features should have little to no multicollinearity;¹⁴ that is, independent variables should be truly independent from one another.

• Logistic regression typically requires a larger sample size compared to other machine learning algorithms like linear regression. Maximum likelihood estimates (used in logistic regression) are less powerful than ordinary least squares (used in linear regression), which results in requiring more training samples to achieve the same statistical learning power.¹⁵

Decision Trees

Decision trees are very versatile supervised learning models that have the important property of being easy to interpret. A decision tree is, as its name suggests, a binary tree data structure that is used to make a decision. Trees are a very intuitive way of displaying and analyzing data and are popularly used even outside of the machine learning field. With the ability to predict both categorical values (classification trees) and real values (regression trees) as well as being able to take in numerical and categorical data without any normalization or dummy variable creation,¹⁶ it’s not difficult to see why they are a popular choice for machine learning.

Let’s see how a typical (top-down) learning decision tree is constructed:

1. Starting at the root of the tree, the full dataset is split based on a binary condition into two child subsets. For example, if the condition is “age ≥ 18,” all data points for which this condition is true go to the left child and all data points for which this condition is false go to the right child.

2. The child subsets are further recursively partitioned into smaller subsets based on other conditions. Splitting conditions are automatically selected at each step based on what condition best splits the set of items. There are a few common metrics by which the quality of a split is measured:

   Gini impurity

   If samples in a subset were randomly labeled according to the distribution of labels in the set, the proportion of samples incorrectly labeled would be the Gini impurity. For example, if a subset were made up of 25% samples with label 0 (and 75% with label 1), assigning label 0 to a random 25% of all samples (and label 1 to the rest) would give 37.5% incorrect labels: 75% of the label-0 samples and 25% of the label-1 samples would be incorrect. A higher-quality decision tree split would split the set into subsets cleanly separated by their label, hence resulting in a lower Gini impurity; that is, the rate of misclassification would be low if most points in a set belong to the same class.

   Variance reduction

   Often used in regression trees, where the dependent variable is continuous. Variance reduction is defined as the total reduction in a set’s variance as a result of the split into two subsets. The best split at a node in a decision tree would be the split that results in the greatest variance reduction.

   Information gain

   Information gain is a measure of the purity of the subsets resulting from a split. It is calculated by subtracting the weighted sum of each decision tree child node’s entropy from the parent node’s entropy. The smaller the entropy of the children, the greater the information gain, hence the better the split.

3. There are a few different methods for determining when to stop splitting nodes:

   • When all leaves of the tree are pure—that is, all leaf nodes each only contain samples belonging to the same class—stop splitting.

   • When a branch of the tree has reached a certain predefined maximum depth, the branch stops being split.

   • When either of the child nodes will contain fewer than the minimum number of samples, the node will not be partitioned.

4. Ultimately the algorithm outputs a tree structure where each node represents a binary decision, the children of each node represent the two possible outcomes of that decision, and each leaf represents the classification of data points following the path from the root to that leaf. (For impure leaves, the decision is determined by majority vote of the training data samples at that leaf.)

An important quality of decision trees is the relative ease of explaining classification or regression results, since every prediction can be expressed in a series of Boolean conditions that trace a path from the root of the tree to a leaf node. For example, if a decision tree model predicted that a malware sample belongs to malware family A, we know it is because the binary was signed before 2015, does not hook into the window manager framework, does make multiple network calls out to Russian IP addresses, etc. Because each sample traverses at most the height of the binary tree (time complexity $O\left(\log n\right)$), decision trees are also efficient to train and make predictions on. As a result, they perform favorably for large datasets.

Nevertheless, decision trees have some limitations:

• Decision trees often suffer from the problem of overfitting, wherein trees are overly complex and don’t generalize well beyond the training set. Pruning is introduced as a regularization method to reduce the complexity of trees.

• Decision trees are more inefficient at expressing some kinds of relationships than others. For example, Figures 2-5 and 2-6 present the minimal decision tree required to represent the AND, OR, and XOR relationships. Notice how XOR requires one more intermediate node and split to be appropriately represented, even in this simple example. For realistic datasets, this can quickly result in exploding model complexity.

  

  Figure 2-5. Decision tree for A-AND-B → y=1 (left), A-OR-B → y=1 (right)

  

  Figure 2-6. Decision tree for A-XOR-B → y=1 (bottom)

• Decision trees tend to be less accurate and robust than other supervised learning techniques. Small changes to the training dataset can result in large changes to the tree, which in turn result in changes to model predictions. This means that decision trees (and most other related models) are unsuitable for use in online learning or incremental learning.

• Split-quality metrics for categorical variables in decision trees are biased toward variables with more possible values; that is, splits on continuous variables or categorical variables with three or more categories will be chosen with a greater probability than binary variables.

• Greedy training of decision trees (as it is almost always done) does not guarantee an optimal decision tree because locally optimal and not globally optimal decisions are made at each split point. In fact, the training of a globally optimal decision tree is an NP-complete problem.¹⁷

Decision Forests

An ensemble refers to a combination of multiple classifiers that creates a more complex, and often better performing, classifier. Combining decision trees into ensembles is a proved technique for creating high-quality classifiers. These ensembles are aptly named decision forests. The two most common types of forests used in practice are decision forests and gradient-boosted decision trees:

• Random forests are formed by simple ensembling of multiple decision trees, typically ranging from tens to thousands of trees. After training each individual decision tree, overall random forest predictions are made by taking the statistical mode of individual tree predictions for classification trees (i.e., each tree “votes”), and the statistical mean of individual tree predictions for regression trees.

  You might notice that simply having many decision trees in the forest will result in highly similar trees and a lot of repeated splits across different trees, especially for features that are strong predictors of the dependent variable. The random forest algorithm addresses this issue using the following training algorithm:

  1. For the training of each individual tree, randomly draw a subset of N samples from the training dataset.

  2. At each split point, we randomly select m features from the p available features, where $m\le p$,¹⁸ and pick the optimal split point from these m features.¹⁹

  3. Repeat step 2 until the individual tree is trained.

  4. Repeat steps 1, 2, and 3 until all trees in the forest are trained.

  Single decision trees tend to overfit to their training sets, and random forests mitigate this effect by taking the average of multiple decision trees, which usually improves model performance. In addition, because each tree in the random forest can be trained independently of all other trees, it is straightforward to parallelize the training algorithm and therefore random forests are very efficient to train. However, the increased complexity of random forests can make them much more storage intensive, and it is much harder to explain predictions than with single decision trees.

• Gradient-boosted decision trees (GBDTs) make use of smarter combinations of individual decision tree predictions to result in better overall predictions. In gradient boosting, multiple weak learners are selectively combined by performing gradient descent optimization on the loss function to result in a much stronger learning model.

  The basic technique of gradient boosting is to add individual trees to the forest one at a time, using a gradient descent procedure to minimize the loss when adding trees. Addition of more trees to the forest stops either when a fixed limit is hit, when validation set loss reaches an acceptable level, or when adding more trees no longer improves this loss.

  Several improvements to basic GBDTs have been made to result in better performing, better generalizing, and more efficient models. Let’s look at a handful of them:

  1. Gradient boosting requires weak learners. Placing artificial constraints on trees, such as limits on tree depth, number of nodes per tree, or minimum number of samples per node, can help constrain these trees without overly diminishing their learning ability.

  2. It can happen that the decision trees added early on in the additive training of gradient-boosted ensembles contribute much more to the overall prediction than the trees added later in the process. This situation results in an imbalanced model that limits the benefits of ensembling. To solve this problem, the contribution of each tree is weighted to slow down the learning process, using a technique called shrinkage, to reduce the influence of individual trees and allow future trees to further improve the model.

  3. We can combine the stochasticity of random forests with gradient boosting by subsampling the dataset before creating a tree and subsampling the features before creating a split.

  4. We can use standard and popular regularization techniques such as $L_1$ and $L_2$ regularization to smooth final learned weights to further avoid overfitting.

  XGBoost²⁰ is a popular GBDT flavor that achieves state-of-the-art results while scaling well to large datasets. As the algorithm that was responsible for many winning submissions to machine learning competitions, it garnered the attention of the machine learning community and has become the decision forest algorithm of choice for many practitioners. Nevertheless, GBDTs are more prone to overfitting than random forests, and also more difficult to parallelize because they use additive training, which relies on the results of a given tree to update gradients for the subsequent tree. We can mitigate overfitting of GBDTs by using shrinkage, and we can parallelize training within a single tree instead of across multiple trees.

Support Vector Machines

Like logistic regression, a support vector machine (SVM) is (in its simplest form) a linear classifier, which means that it produces a hyperplane in a vector space that attempts to separate the two classes in the dataset. The difference between logistic regression and SVMs is the loss function. Logistic regression uses a log-likelihood function that penalizes all points proportionally to the error in the probability estimate, even those on the correct side of the hyperplane. An SVM, on the other hand, uses a hinge loss, which penalizes only those points on the wrong side of the hyperplane or very near it on the correct side.

More specifically, the SVM classifier attempts to find the maximum-margin hyperplane separating the two classes, where “margin” indicates the distance from the separating plane to the closest data points on each side. For the case in which the data is not linearly separable, points within the margin are penalized proportionately to their distance from the margin. Figure 2-7 shows a concrete example: the two classes are represented by white and black points, respectively. The solid line is the separating plane and the dashed lines are the margins. The square points are the support vectors; that is, those that provide nonzero contribution to the loss function. This loss function is expressed mathematically as:

where $\beta$ is the margin, ${\xi }_i$ is the distance from the ith support vector to the margin, and C is a model hyperparameter that determines the relative contribution of the two terms.

Figure 2-7. Classification boundary (dark line) and margins (dashed lines) for linear SVM separating two classes (black and white points); squares represent support vectors

To classify a new data point x, we simply determine which side of the plane x falls on. If we want to get a real-valued score we can compute the distance from x to the separating plane and then apply a sigmoid to map to [0,1].

The real power of SVMs comes from the kernel trick, which is a mathematical transformation that takes a linear decision boundary and produces a nonlinear boundary. At a high level, a kernel transforms one vector space, $V_1$, to another space, $V_2$. Mathematically, the kernel is a function on $V_1\times V_1$ defined by $K(x,y)$, and each $x\in V_1$ is mapped to the function $K(x,·)$; $V_2$ is the space spanned by all such functions. For example, we can recover the linear SVM by defining K to be the usual dot product $K(x,y)=x·y$.

If the kernel is nonlinear, a linear classifier in $V_1$ will produce a nonlinear classifier in $V_2$. The most popular choice is the radial basis function $K(x,y)=e^{-\gamma |x-y|}$. Even though we omit the mathematical details here,²¹ you can think of an SVM with the RBF kernel as producing a sort of smoothed linear combination of spheres around each point x, where the inside of each sphere is classified the same as x, and the outside is assigned the opposite class. The parameter $\gamma$ determines the radius of these spheres; i.e., how close to each point you need to be in order to be classified like that point. The parameter C determines the “smoothing”; a large value of C will lead to the classifier consisting of a union of spheres, while a small value will yield a wigglier boundary influenced somewhat by each sphere. We can therefore see that too large a value of C will produce a model that overfits the training data, whereas too small a value will give a poor classifier in terms of accuracy. (Note that $\gamma$ is unique to the RBF kernel, whereas C exhibits the same properties for any kernel, including the linear SVM. Optimal values of these parameters are usually found using grid search.)

SVMs have shown very good performance in practice, especially in high-dimensional spaces, and the fact that they can be described in terms of support vectors leads to efficient implementations for scoring new data points. However, the complexity of training a kernelized SVM grows quadratically with the number of training samples, so that for training set sizes beyond a few million, kernels are rarely used and the decision boundary is linear. Another disadvantage is that the scores output by SVMs are not interpretable as probabilities; converting scores to probabilities requires additional computation and cross-validation, for example using Platt scaling or isotonic regression. The scikit-learn documentation has further details.

Naive Bayes

The Naive Bayes classifier is one of the oldest statistical classifiers. The classifier is called “naive” because it makes a very strong statistical assumption, namely that features are chosen independently from some (unknown) distribution. This assumption never actually holds in real life. For example, consider a spam classifier for which the features are the words in the message. The Naive Bayes assumption posits that a spam message is composed by sampling words independently, where each word w has a probability $p_{w,spam}$ of being sampled, and similarly for good messages. This assumption is clearly ludicrous; for one thing, it completely ignores word ordering. Yet despite the fact that the assumption doesn’t hold, Naive Bayes classifiers have been shown to be quite effective for problems such as spam classification.

The main idea behind Naive Bayes is as follows: given a data point with feature set $X=x_1,...,x_n$, we want to determine the probability that the label Y for this point is the class C. In Equation 2-1, this concept is expressed as a conditional probability.

Now using Bayes’ Theorem, this probability can be reexpressed as in Equation 2-2.

If we make the very strong assumption that for samples in each class the features are chosen independently of one another, we get Equation 2-3.

we can estimate the numerator of Equation 2-2 from labeled data: $Pr[X_i=x_i\mid Y=C]$ is simply the fraction of samples with the ith feature equal to $x_i$ out of all the samples in class C, whereas $Pr[Y=C]$ is the fraction of samples in class C out of all the labeled samples.

What about the denominator of Equation 2-2? It turns out that we don’t need to compute this, because in two-class classification it is sufficient to compute the ratio of the probability estimates for the two classes $C_1$ and $C_2$. This ratio (Equation 2-4) gives a positive real number score.

A score of $\theta$ > 1 indicates that $C_1$ is the more likely class, whereas $\theta$ < 1 indicates that $C_2$ is more likely. (If optimizing for one of precision or recall you might want to choose a different threshold for the classification boundary.)

The astute observer will notice that we obtained our score $\theta$ without reference to a loss function or an optimization algorithm. The optimization algorithm is actually hidden in the estimate of $Pr[X_i=x_i\mid Y=C]$; using the fraction of samples observed in the training data gives the maximum likelihood estimate, the same loss function used for logistic regression. The similarity to logistic regression doesn’t stop there: if we take logarithms of Equation 2-4 the righthand side becomes a linear function of the features, so we can view Naive Bayes as a linear classifier, as well.

A few subtleties arise when trying to use Naive Bayes in practice:

• What happens when all of the examples of some feature are in the same class (e.g. brand names of common sex enhancement drugs appear only in spam messages)? Then, one of the terms of Equation 2-4 will be zero, leading to a zero or infinity estimate for $\theta$, which doesn’t make sense. To get around this problem we use smoothing, which means adding “phantom” samples to the labeled data for each feature. For example, if we are smoothing by a factor of $\alpha$, we would calculate the value for a feature as follows:

  $$Pr[X_i=x_i\mid Y=C]=\frac{\left(\text{\# samples in class}C\text{with feature}{x}_i\right)+\alpha }{\left(\text{\# samples in class}C\right)+\alpha ·\left(\text{\# of features}\right)}$$

  The choice $\alpha$ = 1 is called Laplace smoothing, whereas $\alpha$ < 1 is called Lidstone smoothing.

• What happens if a feature $x_i$ appears in our validation set (or worse, in real-life scoring) that did not appear in the training set? In this case we have no estimate for $Pr[X_i=x_i\mid Y=C]$ at all. A naive estimate would be to set the probability to $Pr[Y=C]$; for more sophisticated approaches see the work of Freeman.²²

As a final note, we can map the score $\theta \in (0,\infty )$ to a probability in (0,1) using the mapping $\theta \to \frac{\theta }{1+\theta }$; however, this probability estimate will not be properly calibrated. As with SVMs, to obtain better probability estimates we recommend techniques such as Platt scaling or isotonic regression.

k-Nearest Neighbors

The k-nearest neighbors (k-NN) algorithm is the most well-known example of a lazy learning algorithm. This type of machine learning technique puts off most computations to classification time instead of doing the work at training time. Lazy learning models don’t learn generalizations of the data during the training phase. Instead, they record all of the training data points they are passed and use this information to make the local generalizations around the test sample during classification. k-NN is one of the simplest machine learning algorithms:

• The training phase simply consists of storing all the feature vectors and corresponding sample labels in the model.

• The classification prediction²³ is simply the most common label out of the test sample’s k nearest neighbors (hence the name).

The distance metrics for determining how “near” points are to each other in an n-dimensional feature space (where n is the size of the feature vectors) are typically the Euclidean distance for continuous variables and the Hamming distance for discrete variables.

As you might imagine, with such a simple algorithm the training phase of k-NN is typically very fast compared to other learning algorithms, at the cost of classification processing time. Also, the fact that all feature vectors and labels need to be stored within the model results in a very space-inefficient model. (A k-NN model that takes in 1 GB of training feature vectors will at least be 1 GB in size.)

The simplicity of k-NN makes it a popular example for teaching the concept of machine learning to novices, but it is rarely seen in practical scenarios because of the serious drawbacks it has. These include:

• Large model sizes, because models must store (at least) all training data feature vectors and labels.

• Slow classification speeds, because all generalization work is pushed off until classification time. Searching for the nearest neighbors can be time consuming, especially if the model stores training data points in a manner that is not optimized for spatial search. k-d trees (explained in “k-d trees”) are often used as an optimized data structure to speed up neighbor searches.²⁴

• High sensitivity to class imbalance²⁵ in the dataset. Classifications will be skewed towards the classes with more samples in the training data since there is a greater likelihood that samples of these classes will make it into the k-NN set of any given test sample.

• Diminished classification accuracy due to noisy, redundant, or unscaled features. (Choosing a larger k reduces the effect of noise in the training data, but also can result in a weaker learner.)

• Difficulty in choosing the parameter k. Classification results are highly dependent on this parameter, and it can be difficult to choose a k that works well across all parts of the feature space because of differing densities within the dataset.

• Breaks down in high dimensions due to the “curse of dimensionality.” In addition, with more dimensions in the feature space, the “neighborhood” of any arbitrary point becomes larger, which results in noisier neighbor selection.

Neural Networks

Artificial neural networks (ANNs) are a class of machine learning techniques that have seen a resurgence in popularity recently. One can trace the origins of neural networks all the way back to 1942, when McCulloch and Pitts published a groundbreaking paper postulating how neurons in the human nervous system might work.²⁶ Between then and the 1970s, neural network research advanced at a slow pace, in large part due to von Neumann computing architectures (which are quite in opposition to the idea of ANNs) being in vogue. Even after interest in the field was renewed in the 1980s, research was still slow because the computational requirements of training these networks meant that researchers often had to wait days or weeks for the results of their experiments. What triggered the recent popularity of neural networks was a combination of hardware advancements—namely graphics processing units (GPUs) for “almost magically” parallelizing and speeding up ANN training—and the availability of the huge amounts of data that ANNs need to get good at complex tasks like image and speech recognition.

The human brain is composed of a humongous number of neurons (on the order of 10 billion), each with connections to tens of thousands of other neurons. Each neuron receives electrochemical inputs from other neurons, and if the sum of these electrical inputs exceeds a certain level, the neuron then triggers an output transmission of another electrochemical signal to its attached neurons. If the input does not exceed this level, the neuron does not trigger any output. Each neuron is a very simple processing unit capable of only very limited functionality, but in combination with a large number of other neurons connected in various patterns and layers, the brain is capable of performing extremely complex tasks ranging from telling apart a cat from a dog to grasping profound philosophical concepts.

ANNs were originally attempts at modeling neurons in the brain to achieve human-like learning. Individual neurons were modeled with simple mathematical step functions (called activation functions), taking in weighted input from some neurons and emitting output to some other neurons if triggered. This mathematical model of a biological neuron is also called a perceptron. Armed with perceptrons, we then can form a plethora of different neural networks by varying the topology, activation functions, learning objectives, or training methods of the model.

Typically, ANNs are made up of neurons arranged in layers. Each neuron in a layer receives input from the previous layer and, if activated, emits output to one or more neurons in the next layer. Each connection of two neurons is associated with a weight, and each neuron or layer might also have an associated bias. These are the parameters to be trained by the process of backpropagation, which we describe simply and briefly. Before starting, all of the weights and biases are randomly initialized. For each sample in the training set, we perform two steps:²⁷

1. Forward pass. Feed the input through the ANN and get the current prediction.

2. Backward pass. If the prediction is correct, reward the connections that produced this result by increasing their weights in proportion to the confidence of the prediction. If the prediction is wrong, penalize the connections that contributed to this wrong result.

Neural networks have been extensively used in industry for many years and are backed by large bodies of academic research. There are hundreds of variations of neural network infrastructures, and we can use them for both supervised and unsupervised learning. An important quality of some ANNs is that they can perform unsupervised feature learning (which is different from unsupervised learning), which means that minimal or no feature engineering is required. For example, building a malware classifier with an SVM requires domain experts to generate features (as we discuss in Chapter 4). If we use ANNs, it is possible to feed the raw processor instructions or call graphs into the network and have the network itself figure out which features are relevant for the classification task.

Even though there are a lot of hardware and software optimizations available for the training of neural networks, they are still significantly more computationally expensive than training a decision tree, for instance. (Predictions, on the other hand, can be made rather efficiently because of the layered structure.) The number of model hyperparameters to tune in the training of ANNs can be enormous. For instance, you must choose the configuration of network architecture, the type of activation function, whether to fully connect neurons or leave layers sparsely connected, and so on. Lastly, although ANNs—or deep learning networks, as the deep (many-layered) variants of neural networks are popularly called—are frequently thought to be a silver bullet in machine learning, they can be very complex and difficult to reason about. There are often alternatives (such as the other algorithms that we discussed in this section) that are faster, simpler, and easier to train, understand, and explain.

Selecting a Model Family

In the previous section, we discussed a number of different supervised classification algorithms. How are you supposed to decide which one to pick for a given task? The best answer is to let the data decide: try a number of different approaches and see what works best. If you don’t have the time or infrastructure to do this, here are some factors to consider:

Computational complexity

You want to be able to train your model in a reasonable amount of time, using all of the data you have. Logistic regression models and linear SVMs can be trained very quickly. In addition, logistic regression and decision forests have very efficient parallel implementations that can handle very large amounts of data. Kernel SVMs and neural networks can take a long time to train.

Mathematical complexity

Your data might be such that a linear decision boundary will provide good classification; on the other hand, many datasets have nonlinear boundaries. Logistic regression, linear SVM, and Naive Bayes are all linear algorithms; for nonlinear boundaries you can use decision forests, kernel SVMs, or neural networks.

Explainability

A human might need to read the output of your model and determine why it made the decision it did; for example, if you block a legitimate user, you should educate them as to what they did that was suspicious. The best model family for explainability is a decision tree, because it says exactly why the classification was made. We can use relative feature weights to explain logistic regression and Naive Bayes. Decision forests, SVMs, and neural networks are very difficult to explain.

There are many opinions on what model to use in which situation; an internet search for “machine learning cheat sheet” will produce dozens. Our basic recommendation is to use decision forests for high accuracy on a moderate number of features (up to thousands), and logistic regression for fast training on a very large number of features (tens of thousands or more). But ultimately, you should experiment to find the best choice for your data.

Training Data Construction

In a supervised learning problem, you will (somehow) have assembled labeled examples of the thing that you want to classify—account registrations, user logins, email messages, or whatever. Because your goal is to come up with a model that can predict the future based on past data, you will need to reserve some of your labeled data to act as this “future” data so that you can evaluate your model. There are a number of different ways you can do this:

Cross-validation

This technique is a standard method for evaluating models when there is not a lot of training data, so every labeled example makes a significant contribution to the learned model. In this method, you divide the labeled data into k equal parts (usually 5 or 10) and train k different models: each model “holds out” a different one of the k parts and trains on the remaining k–1 parts. The held-out part is then used for validation. Finally, the k different models are combined by averaging both the performance statistics and the model parameters.

Train/validate/test

This approach is most common when you have enough training data so that removing up to half of it doesn’t affect your model very much. In this approach you randomly divide your labeled data into three parts. (See the left side of Figure 2-8.) The majority (say, 60%) is the training set and is input to the learning algorithm. The second part (say, 20%) is the validation set and is used to evaluate and iterate on the model output by the learning algorithm. For example, you can tune model parameters based on performance on the validation set. Finally, after you have settled on an optimal model, you use the remaining part (say, 20%) as the test set to estimate real-world performance on unseen data.

Out-of-time validation

This approach addresses the problem of validating models when the data distribution changes over time. In this case, sampling training and validation sets at random from the same labeled set is “cheating” in some sense—the distributions of the training and validation sets will match closely and not reflect the time aspect. A better approach is to split the training and validation sets based on a time cutoff: samples before time t are in the training set and samples after time t are in the validation set. (See the right side of Figure 2-8.) You might want to have the test set come from the same time period as the validation set, or a subsequent period entirely.

Figure 2-8. In-time validation (left) and out-of-time validation (right)

The composition of the training set might also require some attention. Let’s take a look at a few common issues.

Feature Selection

If you have a reasonably efficient machine learning infrastructure, most of your time and energy will be spent on feature engineering—figuring out signals that you can use to identify attacks, and then building them into your training and scoring pipeline. To make best use of your effort, you want to use only features that provide high discriminatory power; the addition of each feature should noticeably improve your model.

In addition to requiring extra effort to build and maintain, redundant features can hurt the quality of your model. If the number of features is greater than the number of data points, your model will be overfit: there are enough model parameters to draw a curve through all the training data. In addition, highly correlated features can lead to instability in model decisions. For example, if you have a feature that is “number of logins yesterday” and one that is “number of logins in the last two days,” the information you are trying to collect will be split between the two features essentially arbitrarily, and the model might not learn that either of these features is important.

You can solve the feature correlation problem by computing covariance matrices between your features and combining highly correlated features (or projecting them into orthogonal spaces; in the previous example, “number of logins in the day before yesterday” would be a better choice than “number of logins in the last two days”).

There are number of techniques to address the feature selection problem:

• Logistic regression, SVMs, and decision trees/forests have methods of determining relative feature importance; you can run these and keep only the features with highest importance.

• You can use $L_1$ regularization (see the next section) for feature selection in logistic regression and SVM classifiers.

• If the number n of features is reasonably small (say, n < 100) you can use a “build it up” approach: build n one-feature models and determine which is best on your validation set; then build n–1 two-feature models, and so on, until the gain of adding an additional feature is below a certain threshold.

• Similarly, you can use a “leave one out” approach: build a model on n features, then n models on n–1 features and keep the best, and so on until the loss of removing an additional feature is too great.

scikit-learn implements the sklearn.feature_selection.SelectFromModel helper utility that assists operators in selecting features based on importance weights. As long as a trained estimator has the feature_importances_ or coef_ attribute after fitting,²⁸ it can be passed into SelectFromModel for feature selection importance ranking. Assuming that we have a pretrained DecisionTreeClassifier model (variable name clf) and an original training dataset (variable name train_x) with 119 features, here is a short code snippet showing how to use SelectFromModel to keep only features with a feature_importance that lies above the mean:^29,30

from sklearn.feature_selection import SelectFromModel

sfm = SelectFromModel(clf, prefit=True)

# Generate new training set, keeping only the selected features
train_x_new = sfm.transform(train_x)

print("Original num features: {}, selected num features: {}"
      .format(train_x.shape[1], train_x_new.shape[1]))

> Original num features: 119, selected num features: 7

Choosing Thresholds and Comparing Models

The supervised classification algorithm you choose will typically output a real-valued score, and you will need to choose a threshold or thresholds³¹ above which to block the activity or show additional friction (e.g., require a phone number). How do you choose this threshold? This choice is ultimately a business decision, based on the trade-off between security and user friction. Ideally you can come up with some cost function—for example, 1 false positive is equal to 10 false negatives—and minimize the total cost on a representative sample of the data. Another option is to fix a precision or recall target—for example, 98%—and choose a threshold that achieves that target.

Now suppose that you have two versions of your model with different parameters (e.g., different regularization) or even different model families (e.g., logistic regression versus random forest). Which one is better? If you have a cost function this is easy: compute the cost of the two versions on the same dataset and choose the lower-cost option. If fixing a precision target, choose the version that optimizes recall, and vice versa.

Another common method for model comparison is to plot the receiver operating characteristic (ROC) curve and compute the area under the curve (AUC). The ROC curve plots false positive rate (FP / (FP + TN)) on the x-axis and true positive rate (TP / (TP + FN), also known as recall) on the y-axis. Each point on the curve corresponds to a score threshold and represents the (FPR, TPR) pair at that threshold. The AUC can be interpreted as the probability that a randomly chosen positive example has a higher score than a randomly chosen negative example; under this interpretation it’s easy to see that the worst case is AUC 0.5, which is equivalent to a random ordering of samples.³²

Figure 2-10 shows an example ROC curve, with the line y = x plotted for comparison. Because the AUC is very high, we have used a log scale to zoom in on the lefthand side, which is where differences between high-performance models will appear; if you are operating only in the high-precision region of the curve, you may want to calculate up to a threshold false positive rate, such as 1%.

One nice property of AUC is that it is unaffected by sampling bias. Thus, if you sample two classes with different weights, the AUC you get on the resulting dataset will be representative of the AUC on the unsampled dataset.

Figure 2-10. ROC curves

One common metric that is limited in real-world usefulness is the F-score, which is defined as follows:

The F-score combines precision and recall and harshly penalizes extremes; however, it requires choosing a threshold and a relative weighting of precision and recall (parametrized by $\alpha$).

Clustering Algorithms

The geometric intuition behind clustering is straightforward: you want to group together data points that are “close together” in some sense. Thus, for any algorithm to work you need to have some concrete way to measure “closeness”; such a measurement is called a metric. The metric and clustering algorithm you use will depend on the form your data is in; for example, your data might consist of real-valued vectors, lists of items, or sequences of bits. We now consider the most popular algorithms.

k-means

k-means is usually the first algorithm that comes to mind when you think of clustering. k-means applies to real-valued vectors, when you know how many clusters you expect; the number of clusters is denoted by k. The goal of the algorithm is to assign each data point to a cluster such that the sum of the distances from each point to its cluster centroid is minimized. Here the notion of distance is the usual Euclidean distance in a vector space:

$$d(x,y)=\sqrt{{\sum }_i{(x_i-y_i)}^2}$$

In mathematical terms, the k-means algorithm computes a cluster assignment f: X → {1,…k} that minimizes the loss function:

$$L(X)=\sum _{i}d(x_i,c_{f(x_i)})$$

where $X=\{x_1,...,x_n\}$ is your dataset, $c_j$ is the jth centroid, and d is the distance between two points. The value $L\left(X\right)$ is called “inertia.”

The standard algorithm for computing k-means clusters is as follows:

1. Choose k centroids $c_1,...,c_k$ at random.

2. Assign each data point $x_i$ to its nearest centroid.

3. Recompute the centroids $c_j$ by taking the average of all the data points assigned to the jth cluster.

4. Repeat (2) and (3) until the algorithm converges; that is, the difference between $L\left(X\right)$ on successive iterations is below a predetermined threshold.

k-means is a simple and effective clustering algorithm that scales well to very large datasets. However, there are some things for which you need to be on the lookout:

• Since k is a fixed parameter of the algorithm, you must choose it appropriately. If you know how many clusters you are looking for (e.g., if you are trying to cluster different families of malware), you can simply choose k to be that number. Otherwise, you will need to experiment with different values of k. It is also common to choose values of k that are between one to three times the number of classes (labels) in your data, in case some categories are discontinuous. Warning: loss functions computed using different values of k are not comparable to each other!

• You must normalize your data before using k-means. A typical normalization is to map the jth coordinate x_ij to (x_ij-µ_j)/σ_j, where µ_j is the mean of the jth coordinates, and σ_j is the standard deviation.

  To see why normalization is necessary, consider a two-dimensional dataset whose first coordinate ranges between 0 and 1, and whose second coordinate ranges between 0 and 100. Clearly the second coordinate will have a much greater impact on the loss function, so you will lose information about how close together points are in the first coordinate.

• Do not use k-means with categorical features, even if you can represent them as a number. For example, you could encode “red,” “green,” and “blue” as 0, 1, and 2, respectively, but these numbers don’t make sense in a vector space—there is no reason that blue should be twice as far from red as green is. This problem can be addressed by one-hot encoding the categorical features as multiple binary features (as discussed in the worked example earlier in this chapter), but…

• Beware when using k-means with binary features. k-means can sometimes be used with binary features, encoding the two responses as 0 and 1, or −1 and 1, but results here can be unpredictable; the binary feature might become the dominant feature determining the cluster, or its information might be lost entirely.

• k-means loses effectiveness in high dimensions, due to the “curse of dimensionality”—all points are roughly equally distant from each other. For best results use k-means in low dimensions or after applying a dimensionality reduction algorithm such as principal component analysis (PCA). Another option is to use the L-infinity distance, where the distance between two points is taken to be the maximum of the difference of any coordinate:

  $$d(x,y)=ma{x}_i\left(|x_i-y_i|\right)$$

• k-means works best when the initial centroids are chosen at random; however, this choice can make reproducing results difficult. Try different choices of initial centroids to see how the results depend on the initialization.

• k-means assumes that the clusters are spherical (globular) in nature. As you can imagine, it does not work well on non-spherical distributions, such as the one illustrated in Figure 2-11.

Figure 2-11. Nonspherical data distribution

Hierarchical clustering

Unlike the k-means algorithm, hierarchical clustering methods are not parametrized by an operator-selected value k (the number of clusters you want to create). Choosing an appropriate k is a nontrivial task, and can significantly affect clustering results. Agglomerative (bottom-up) hierarchical clustering builds clusters as follows (illustrated in Figure 2-12):

1. Assign each data point to its own cluster (Figure 2-12, bottom layer).

2. Merge the two clusters that are the most similar, where “most similar” is determined by a distance metric such as the Euclidean distance or Mahalanobis distance.

3. Repeat step 2 until there is only one cluster remaining (Figure 2-12, top layer).

4. Navigate the layers of this tree (dendrogram) and select the layer that gives you the most appropriate clustering result.

Divisive (top-down) hierarchical clustering is another form of hierarchical clustering that works in the opposite direction. Instead of starting with as many clusters as there are data points, we begin with a single cluster consisting of all data points and start dividing clusters based on the distance metric, stopping when each data point is in its own separate cluster.

Figure 2-12. Agglomerative hierarchical clustering dendrogram

There are some important points to take note of when considering whether to use hierarchical clustering:

• Hierarchical clustering produces a dendrogram tree model, as illustrated in Figure 2-12. This model can be more complex to analyze and takes up more storage space than the centroids produced by k-means, but also conveys more information about the underlying structure of the data. If model compactness or ease of analysis is a priority, hierarchical clustering might not be your best option.

• k-means works with only a small selection of distance metrics (mostly Euclidean distance) and requires numerical data to work. In contrast, hierarchical clustering works with almost any kind of distance metric or similarity function, as long as it produces a result that can be numerically compared (e.g., C is more similar to A than to B). You can use it with categorical data, mixed type data, strings, images, and so on as long as an appropriate distance function is provided.

• Hierarchical clustering has high time complexity, which makes it unsuitable for large datasets. Taking n to be the number of data points, agglomerative hierarchical clustering has a time complexity of $O(n^{2}log(n))$, and naive divisive clustering has a time complexity of $O(2^n)$.

Evaluating Clustering Results

It can sometimes be difficult to make sense of the results of clustering operations. Evaluating supervised learning algorithms is a much more straightforward task because we have access to the ground truth labels: we can simply count the number of samples to which the algorithm correctly and wrongly assigns labels. In the case of unsupervised learning, it is unlikely we have access to labels, though if we do evaluation becomes much easier.^37,38 For example, commonly used metrics for evaluating clustering results if the ground truth labels are known are:

Homogeneity

The degree to which each cluster contains only members of a single class.

Completeness

The degree to which all members of a certain class are assigned to the same cluster.

The harmonic mean of these two metrics is known as the V-measure,³⁹ an entropy-based score representing the clustering operation’s accuracy, which has the formula:

where h is homogeneity and c is completeness.

We will use the scikit-learn implementations of these clustering evaluation metrics in Chapter 5.

On the other hand, the fact that we are considering clustering at all means it is likely that we do not have any ground truth labels for our dataset. Instead, we can no longer use the V-measure because both homogeneity and completeness can be measured only by comparing prediction labels to ground truth labels. In this case, we have to rely on signals derived directly from the trained model itself. Essentially, we consider the clustering operation to be successful if samples assigned to the same cluster are all similar to one another, and samples assigned to different clusters are all completely dissimilar from one another. There are two popular ways to measure this:

Silhouette coefficient

This score is calculated separately for each sample in the dataset. Using some distance metric (e.g., Euclidean distance), we find the following two mean distances for some sample x:

• a: the mean distance between sample x and all other samples in the same cluster

• b: the mean distance between sample x and all other samples in the next nearest cluster

The Silhouette coefficient s is then defined as follows:⁴⁰

$$s=\frac{b-a}{max(a,b)}$$

For a clustering result that is hugely incorrect, the a will be much larger than b, resulting in the numerator being negative. Because the denominator will never be negative (distance metrics do not take on negative values), s will itself be negative. In fact, s is bounded between −1 and +1. The closer s is to +1, the better the clustering result. When s is very close to 0, we have highly overlapping clusters. Nevertheless, it is found that the Silhouette coefficient is biased toward distance-based (e.g., k-means), grid-based (e.g., STING⁴¹), and hierarchical clustering algorithms that produce convex clusters. It does not work well for clusters produced by density-based algorithms such as DBSCAN and OPTICS, because neither a nor b take cluster density into consideration. Also, because the Silhouette coefficient must be independently calculated for each sample in the dataset, it can be quite slow.

Calinski-Harabaz index

The Calinski-Harabaz (C-H) index⁴² is higher (better) when clusters are dense and well separated. By most measures, this measure of clustering effectiveness is closer to how humans evaluate clustering results. We consider a clustering result to be good if there are visibly separated and densely gathered groups of samples. The C-H index uses two measures:

• $W_k$: the within-cluster dispersion, a matrix of distances between samples in a cluster and the geometric center of the cluster

• $B_k$: the between-group dispersion, a matrix of distances between the center of a cluster and the centers of all other clusters

(Here, k is the number of clusters in the model.)

The C-H score (s) takes the ratio of $W_k$ and $B_k$ (where N is the number of samples in the dataset and tr is the trace of the matrix⁴³):

$$s=\frac{\mathrm{tr}\left(B_k\right)}{\mathrm{tr}\left(W_k\right)}\times \frac{N-k}{k-1}$$

The C-H score can be computed much more efficiently than the Silhouette index. Even though it is still biased against density-based clusters, it is comparatively more reliable than the Silhouette index.

No matter which method you choose, there will be some inherent limitations in evaluating clustering results without any human labeling. If you suspect that clustering models are not performing well in practice in spite of good evaluation metrics, it may be necessary to invest resources in manually labeling a subset of samples and use semi-supervised approaches to evaluate the model.

Host Intrusion Detection

When developing an intrusion detection agent for hosts (e.g., servers, desktops, laptops, embedded systems), you will likely need to generate your own metrics and might even want to perform correlations of signals collected from different sources. The relevance of different metrics varies widely depending on the threat model, but basic system- and network-level statistics make for a good starting point. You can carry out the collection of these system metrics in a variety of ways, and there is a diversity of tools and frameworks to help you with the task. We’ll take a look at osquery, a popular operating system (OS) instrumentation framework that collects and exposes low-level OS metrics, making them available for querying through a SQL-based interface. Making scheduled queries through osquery can allow you to establish a baseline of host and application behavior, thereby allowing the intrusion detector to identify suspicious events that occur unexpectedly.

Malware is the most dominant threat vector for hosts in many environments. Of course, malware detection and analysis warrants its own full discussion, which we provide in Chapter 4. For now, we base our analysis on the assumption that most malware affects system-level actions, and therefore we can detect the malware by collecting system-level activity signals and looking for indicators of compromise (IoCs) in the data. Here are examples of some common signals that you can collect:

• Running processes

• Active/new user accounts

• Kernel modules loaded

• DNS lookups

• Network connections

• System scheduler changes

• Daemon/background/persistent processes

• Startup operations, launchd entries

• OS registry databases, .plist files

• Temporary file directories

• Browser extensions

This list is far from exhaustive; different types of malware naturally generate different sets of behavior, but collecting a wide range of signals will ensure that you have visibility into the parts of the system for which the risk of compromise by malware is the highest.

SELECT * FROM users;

SELECT pid, name, resident_size FROM processes
ORDER BY resident_size DESC LIMIT 5

SELECT * FROM processes WHERE on_disk = 0;

2017-06-04T18:24:17+00:00        []
2017-06-04T18:54:17+00:00        []
2017-06-04T19:24:17+00:00        ["/tmp/YBBHNCA8J0"]
2017-06-04T19:54:17+00:00        []

{
 ...
  // Define a schedule of queries to run periodically:
  "deleted_running_binary": {
    "query": "SELECT * FROM processes WHERE on_disk = 0;",
    "interval": 1800
  }
  ...
}

Network Intrusion Detection

Almost all forms of host intrusion instigate communication with the outside world. Most breaches are carried out with the objective of stealing some valuable data from the target, so it makes sense to detect intrusions by focusing on the network. For botnets, remote command-and-control servers communicate with the compromised “zombie” machines to give instructions on operations to execute. For APTs, hackers can remotely access the machines through a vulnerable or misconfigured service, allowing them shell and/or root access. For adware, communication with external servers is required for downloading unsolicited ad content. For spyware, results of the covert monitoring are often transmitted over the network to an external receiving server.

From simple protocol-tapping utilities like tcpdump to some more complex sniffing tools like Bro, the network intrusion detection software ecosystem has many utilities and application suites that can help you collect signals from network traffic of all sorts. Network intrusion detection tools operate on the basic concept of inspecting traffic that passes between hosts. Just like with host intrusion detection, attacks can be identified either by matching traffic to a known signature of malicious traffic or by anomaly detection, comparing traffic to previously established baselines. In this section, we focus on anomaly detection rather than signature matching; however, we do examine the latter in Chapter 4, which discusses malware analysis in depth.

Snort is a popular open source IDS that sniffs packets and network traffic for real-time anomaly detection. It is the de facto choice for intrusion-detection monitoring, providing a good balance of usability and functionality. Furthermore, it is backed by a vibrant open source community of users and contributors who have created add-ons and GUIs for it. Snort has a relatively simple architecture, allowing users to perform real-time traffic analysis on IP networks, write rules that can be triggered by detected conditions, and compare traffic to an established baseline of the normal network communication profile.

In extracting features for network intrusion detection, there is a noteworthy difference between extracting network traffic metadata and inspecting network traffic content. The former is used in stateful packet inspection (SPI), working at the network and transport layers—OSI layers 3 and 4—and examining each network packet’s header and footer without touching the packet context. This approach maintains state on previous packets received, and hence is able to associate newly received packets with previously seen packets. SPI systems are able to know whether a packet is part of a handshake to establish a new connection, a section of an existing network connection, or an unexpected rogue packet. These systems are useful in enforcing access control—the traditional function of network firewalls—because they have a clear picture of the IP addresses and ports involved in correspondence. They can also be useful in detecting slightly more complex layer 3/4 attacks⁶ such as IP spoofing, TCP/IP attacks (such as ARP cache poisoning or SYN flooding), and denial-of-service (DOS) attacks. However, there are obvious limitations to restricting analysis to just packet headers and footers. For example, SPI cannot detect signs of breaches or intrusions on the application level, because doing so would require a deeper level of inspection.

Web Application Intrusion Detection

We saw earlier that we can detect web application attacks like XSS and SQL injections by using deep network packet inspection tools such as Bro. Inspecting HTTP server logs can provide you with a similar level of information and is a more direct way of obtaining features derived from web application user interactions. Standard web servers like Apache, IIS, and Nginx generate logs in the NCSA Common Log Format, also called access logs. NCSA combined logs and error logs also record information about the client’s user agent, referral URL, and any server errors generated by requests. In these logs, each line represents a separate HTTP request made to the server, and each line is made up of tokens in a well-defined format. Here is an example of a record in the combined log format that includes the requestor’s user agent and referral URL:

123.123.123.123 - jsmith [17/Dec/2016:18:55:05 +0800] "GET /index.html HTTP/1.0"
200 2046 "http://referer.com/" "Mozilla/5.0 (Macintosh; Intel Mac OS X 10.17.3)
AppleWebKit/536.27.14 (KHTML, like Gecko) Chrome/55.0.2734.24 Safari/536.27.14"

Unlike DPI, the standard web access logs do not log POST body data out of the box. This means that attack vectors embedded in the user input cannot be detected by inspecting standard access logs.

Even with the comparatively limited scope of visibility provided in standard HTTP server log files, there are still some interesting features that you can extract:

IP-level access statistics

High frequency, periodicity, or volume by a single IP address or subnet is suspicious.

URL string aberrations

Self-referencing paths (/./) or backreferences (/../) are frequently used in path-traversal attacks.

Decoded URL and HTML entities, escaped characters, null-byte string termination

These are frequently used by simple signature/rule engines to avoid detection.

Unusual referrer patterns

Page accesses with an abnormal referrer URL are often a signal of an unwelcome access to an HTTP endpoint.

Sequence of accesses to endpoints

Out-of-order access to HTTP endpoints that do not correspond to the website’s logical flow is indicative of fuzzing or malicious explorations.

For instance, if a user’s typical access to a website is a POST to /login followed by three successive GETs to /a, /b, and /c, but a particular IP address is repeatedly making GET requests to /b and /c without a corresponding /login or /a request, that could be a sign of bot automation or manual reconnaissance activity.

User agent patterns

You can perform frequency analysis on user agent strings to alert on never-before-seen user agent strings or extremely old clients (e.g., a “Mosaic/0.9” user agent from 1993) which are likely spoofed.

Web logs provide enough information to detect different kinds of attacks on web applications,¹⁰ including, but not limited to, the OWASP Top Ten—XSS, Injection, CSRF, Insecure Direct Object References, etc.

In Summary

Generating a reliable and comprehensive set of features is critical for the anomaly detection process. The goal of feature engineering is to distill complex information into a compact form that removes unnecessary information, but does not sacrifice any important characteristics of the data. These generated features will then be fed into algorithms, which will consume the data and use it to train machine learning models. In the next section, we will see how you can convert feature sets into valuable insights that drive anomaly detection systems.

Forecasting (Supervised Machine Learning)

Forecasting is a highly intuitive way of performing anomaly detection: we learn from prior data and make a prediction about the future. We can consider any substantial deviations between the forecasts and observations as anomalous. Taking the weather as an example, if it had not been raining for weeks, and there was no visible sign of upcoming rain, the forecast would predict a low chance of rain in the coming days. If it did rain in the coming days, it would be a deviation from the forecast.

This class of anomaly detection algorithms uses past data to predict current data, and measures how different the currently observed data is from the prediction. By this definition, forecasting lies in the realm of supervised machine learning because it trains a regression model of data values versus time. Because these algorithms also operate strictly within the notion of past and present, they are suitable only for single-dimension time series datasets. Predictions made by a forecasting model will correspond to the expected value that this time series will have in the next time step, so applying forecasting to datasets other than time series data does not make sense.

Time series data is naturally suited for representation in a line chart. Humans are adept at studying line charts, recognizing trends, and identifying anomalies, but machines have a more difficult time of it. A major reason for this difficulty is the noise embedded within time series data, caused either by measurement inaccuracies, sampling frequency, or other external factors associated with the nature of the data. Noise results in a choppy and volatile series, which can camouflage outbreaks or spikes that we are interested in identifying. In combination with seasonality and cyclic patterns that can sometimes be complex, attempting to use naive linear-fit methods to detect anomalies would likely not give you great results.

In forecasting, it is important to define the following descriptors of time series:

Trends

Long-term direction of changes in the data, undisturbed by relatively small-scale volatility and perturbations. Trends are sometimes nonlinear, but can typically be fit to a low-order polynomial curve.

Seasons

Periodic repetitions of patterns in the data, typically coinciding with factors closely related to the nature of the data; for example, day-night patterns, summer-winter differences, moon phases.

Cycles

General changes in the data that have pattern similarities but vary in periodicity, e.g., long-term stock market cycles.

Figure 3-1 depicts a diurnal-patterned seasonality, with a gentle upward trend illustrated by a regression line fitted to the data.

Figure 3-1. A diurnal season and upward trend

ARIMA

Using the ARIMA (autoregressive integrated moving average) family of functions is a powerful and flexible way to perform forecasting on time series. Autoregressive models are a class of statistical models that have outputs that are linearly dependent on their own previous values in combination with a stochastic factor.¹² You might have heard of exponential smoothing, which can often be equated/approximated to special cases of ARIMA (e.g., Holt-Winters exponential smoothing). These operations smooth jagged line charts, using different variants of weighted moving averages to normalize the data. Seasonal variants of these operations can take periodic patterns into account, helping make more accurate forecasts. For instance, seasonal ARIMA (SARIMA) defines both a seasonal and a nonseasonal component of the ARIMA model, allowing periodic characteristics to be captured.¹³

In choosing an appropriate forecasting model, always visualize your data to identify trends, seasonalities, and cycles. If seasonality is a strong characteristic of the series, consider models with seasonal adjustments such as SARIMA and seasonal Holt-Winters methods. Forecasting methods learn characteristics of the time series by looking at previous points and making predictions about the future. In exploring the data, a useful metric to learn is the autocorrelation, which is the correlation between the series and itself at a previous point in time. A good forecast of the series can be thought of as the future points having high autocorrelation with the previous points. ARIMA uses a distributed lag model in which regressions are used to predict future values based on lagged values (an autoregressive process). Autoregressive and moving average parameters are used to tune the model, along with polynomial factor differencing—a process used to make the series stationary (i.e., having constant statistical properties over time, such as mean and variance), a condition that ARIMA requires the input series to have.

In this example, we attempt to perform anomaly detection on per-minute metrics of a host’s CPU utilization.¹⁴ The y-axis of Figure 3-2 shows the percentage CPU utilization, and the x-axis shows time.

Figure 3-2. CPU utilization over time

We can observe a clear periodic pattern in this series, with peaks in CPU utilization roughly every 2.5 hours. Using a convenient time series library for Python, PyFlux, we apply the ARIMA forecasting algorithm with autoregressive (AR) order 11, moving average (MA) order 11, and a differencing order of 0 (because the series looks stationary).¹⁵ There are some tricks to determining the AR and MA orders and the differencing order, which we will not elaborate on here. To oversimplify matters, AR and MA orders are needed to correct any residual autocorrelations that remain in the differenced series (i.e., between the time-shifted series and itself). The differencing order is a term used to make the series stationary—an already stationary series should have a differencing order of 0, a series with a constant average trend (steadily trending upward or downward) should have a differencing order of 1, and a series with a time-varying trend (a trend that changes in velocity and direction over the series) should have a differencing order of 2. Let’s plot the in-sample fit to get an idea of how the algorithm does:¹⁶

import pandas as pd
import pyflux as pf
from datetime import datetime

# Read in the training and testing dataset files
data_train_a = pd.read_csv('cpu-train-a.csv',
    parse_dates=[0], infer_datetime_format=True)
data_test_a = pd.read_csv('cpu-test-a.csv',
    parse_dates=[0], infer_datetime_format=True)

# Define the model
model_a = pf.ARIMA(data=data_train_a,
                   ar=11, ma=11, integ=0, target='cpu')

# Estimate latent variables for the model using the
# Metropolis-Hastings algorithm as the inference method
x = model_a.fit("M-H")

# Plot the fit of the ARIMA model against the data
model_a.plot_fit()

Figure 3-3 presents the result of the plot.

Figure 3-3. CPU utilization over time fitted with ARIMA model prediction

As we can observe in Figure 3-3, the results fit the observed data quite well. Next, we can do an in-sample test on the last 60 data points of the training data. The in-sample test is a validation step that treats the last subsection of the series as unknown and performs forecasting for those time steps. This process allows us to evaluate performance of the model without running tests on future/test data:

> model_a.plot_predict_is(h=60)

The in-sample prediction test (depicted in Figure 3-4) looks pretty good because it does not deviate from the original series significantly in phase and amplitude.

Figure 3-4. In-sample (training set) ARIMA prediction

Now, let’s run the actual forecasting, plotting the most recent 100 observed data points followed by the model’s 60 predicted values along with their confidence intervals:

> model_a.plot_predict(h=60, past_values=100)

Bands with a darker shade imply a higher confidence; see Figure 3-5.

Figure 3-5. Out-of-sample (test-set) ARIMA prediction

Comparing the prediction illustrated in Figure 3-5 to the actual observed points illustrated in Figure 3-6, we see that the prediction is spot-on.

Figure 3-6. Actual observed data points

To perform anomaly detection using forecasting, we compare the observed data points with a rolling prediction made periodically. For example, an arbitrary but sensible system might make a new 60-minute forecast every 30 minutes, training a new ARIMA model using the previous 24 hours of data. Comparisons between the forecast and observations can be made much more frequently (e.g., every three minutes). We can apply this method of incremental learning to almost all the algorithms that we will discuss, which allows us to approximate streaming behavior from algorithms originally designed for batch processing.

Let’s perform the same forecasting operations on another segment of the CPU utilization dataset captured at a different time:

data_train_b = pd.read_csv('cpu-train-b.csv',
    parse_dates=[0], infer_datetime_format=True)
data_test_b = pd.read_csv('cpu-test-b.csv',
    parse_dates=[0], infer_datetime_format=True)

Forecasting using the same ARIMAX model¹⁷ trained on data_train_b, the prediction is illustrated in Figure 3-7.

Figure 3-7. Out-of-sample (test-set, data_train_b) ARIMAX prediction

The observed values are, however, very different from the predictions illustrated in Figure 3-8.

Figure 3-8. Actual observed data points (data_train_b)

We see a visible anomaly that occurs a short time after our training period. Because the observed values fall within the low-confidence bands, we will raise an anomaly alert. The specific threshold conditions for how different the forecasted and observed series must be to raise an anomaly alert is something that is highly application specific but should be simple enough to implement on your own.

Artificial neural networks

Another way to perform forecasting on time series data is to use artificial neural networks. In particular, long short-term memory (LSTM) networks^18,19 are suitable for this application. LSTMs are a variant of recurrent neural networks (RNNs) that are uniquely architected to learn trends and patterns in time series input for the purposes of classification or prediction. We will not go into the theory or implementation details of neural networks here; instead, we will approach them as black boxes that can learn information from time series containing patterns that occur at unknown or irregular periodicities. We will use the Keras LSTM API, backed by TensorFlow, to perform forecasting on the same CPU utilization dataset that we used earlier.

The training methodology for our LSTM network is fairly straightforward. We first extract all continuous length-n subsequences of data from the training input, treating the last point in each subsequence as the label for the sample. In other words, we are generating n-grams from the input. For example, taking n = 3, for this raw data:

raw: [0.51, 0.29, 0.14, 1.00, 0.00, 0.13, 0.56]

we get the following n-grams:

n-grams: [[0.51, 0.29, 0.14],
          [0.29, 0.14, 1.00],
          [0.14, 1.00, 0.00],
          [1.00, 0.00, 0.13],
          [0.00, 0.13, 0.56]]

and the resulting training set is:

sample	label
(0.51, 0.29) (0.29, 0.14) (0.14, 1.00) (1.00, 0.00) (0.00, 0.13)	0.14 1.00 0.00 0.13 0.56

The model is learning to predict the third value in the sequence following the two already seen values. LSTM networks have a little more complexity built in that deals with remembering patterns and information from previous sequences, but as mentioned before, we will leave the details out. Let’s define a four-layer²⁰ LSTM network:²¹

from keras.models import Sequential
from keras.layers.recurrent import LSTM
from keras.layers.core import Dense, Activation, Dropout

# Each training data point will be length 100-1,
# since the last value in each sequence is the label
sequence_length = 100

model = Sequential()

# First LSTM layer defining the input sequence length
model.add(LSTM(input_shape=(sequence_length-1, 1),
               units=32,
               return_sequences=True))
model.add(Dropout(0.2))

# Second LSTM layer with 128 units
model.add(LSTM(units=128,
               return_sequences=True))
model.add(Dropout(0.2))

# Third LSTM layer with 100 units
model.add(LSTM(units=100,
               return_sequences=False))
model.add(Dropout(0.2))

# Densely connected output layer with the linear activation function
model.add(Dense(units=1))
model.add(Activation('linear'))

model.compile(loss='mean_squared_error', optimizer='rmsprop')

The precise architecture of the network (number of layers, size of each layer, type of layer, etc.) is arbitrarily chosen, roughly based on other LSTM networks that work well for similar problems. Notice that we are adding a Dropout(0.2) term after each hidden layer—dropout²² is a regularization technique that is commonly used to prevent neural networks from overfitting. At the end of the model definition, we make a call to the model.compile() method, which configures the learning process. We choose the rmsprop optimizer because the documentation states that it is usually a good choice for RNNs. The model fitting process will use the rmsprop optimization algorithm to minimize the loss function, which we have defined to be the mean_squared_error. There are many other tunable knobs and different architectures that will contribute to model performance, but, as usual, we opt for simplicity over accuracy.

Let’s prepare our input:

...
# Generate n-grams from the raw training data series
n_grams = []
for ix in range(len(training_data)-sequence_length):
n_grams.append(training_data[ix:ix+sequence_length])

# Normalize and shuffle the values
n_grams_arr = normalize(np.array(n_grams))
np.random.shuffle(n_grams_arr)

# Separate each sample from its label
x = n_grams_arr[:, :-1]
labels = n_grams_arr[:, −1]
...

Then, we can proceed to run the data through the model and make predictions:

...
model.fit(x,
   labels,
   batch_size=50,
   nb_epochs=3,
   validation_split=0.05)

y_pred = model.predict(x_test)
...

Figure 3-9 shows the results alongside the root-mean-squared (RMS) deviation.

We see that the prediction follows the nonanomalous observed series closely (both normalized), which hints to us that the LSTM network is indeed working well. When the anomalous observations occur, we see a large deviation between the predicted and observed series, evident in the RMS plot. Similar to the ARIMA case, such measures of deviations between predictions and observations can be used to signal when anomalies are detected. Thresholding on the observed versus predicted series divergence is a good way to abstract out the quirks of the data into a simple measure of unexpected deviations.

Figure 3-9. Observed, predicted, and RMS deviation plots of LSTM anomaly detection applied on the CPU time series

Statistical Metrics

We can use statistical tests to determine whether a single new data point is similar to the previously seen data. Our example at the beginning of the chapter, in which we made a threshold-based anomaly detector adapt to changing data by maintaining an aggregate moving average of the series, falls into this category. We can use moving averages of time series data as an adaptive metric that indicates how well data points conform to a long-term trend. Specifically, the moving average (also known as a low-pass filter in signal processing terminology) is the reference point for statistical comparisons, and significant deviations from the average will be considered anomalies. Here we briefly discuss a few noteworthy metrics, but we do not dwell too long on each, because they are fairly straightforward to use.

import numpy as np

# Input data series
x = [1, 2, 3, 4, 5, 6]

# Calculate median absolute deviation
mad = np.median(np.abs(x - np.median(x)))

# MAD of x is 1.5

Goodness-of-Fit

In building an anomaly detection system, it is important to consider whether the data used to train the initial model is contaminated with anomalies. As discussed earlier, this question can be difficult to answer, but you can often make an informed guess given a proper understanding of the nature of the data source and threat model. In a perfect world, the expected distributions of a dataset can be accurately modeled with known distributions. For instance, the distribution of API calls to an application server per day (over time) might closely fit a normal distribution, and the number of hits to a website in the hours after a promotion is launched might be accurately described by an exponential decay. However, because we do not live in a perfect world, it is rare to find real datasets that conform perfectly to simple distributions. Even if a dataset can be fitted to some hypothetical analytical distribution, accurately determining what this distribution is can be a challenge. Nevertheless, this approach can be feasible in some cases, especially when dealing with a large dataset for which the expected distribution is well known.²⁵ In such cases, comparing the divergence between the expected and observed distributions can be a method of anomaly detection.

Goodness-of-fit tests such as the chi-squared test, the Kolmogorov–Smirnov test, and the Cramér–von Mises criterion can be used to quantify how similar two continuous distributions are. These tests are mostly only valid for one-dimensional datasets, however, which can largely limit their usefulness. We will not dive too deeply into traditional goodness-of-fit tests because of their limited usefulness in real-world anomaly detection. Instead, we will take a closer look at more versatile methods such as the elliptic envelope fitting method provided in scikit-learn.

Elliptic envelope fitting (covariance estimate fitting)

For normally distributed datasets, elliptic envelope fitting can be a simple and elegant way to perform anomaly detection. Because anomalies are, by definition, points that do not conform to the expected distribution, it is easy for these algorithms to exclude such outliers in the training data. Thus, this method is only minimally affected by the presence of anomalies in the dataset.

The use of this method requires that you make a rather strong assumption about your data—that the inliers come from a known analytical distribution. Let’s take an example of a hypothetical dataset containing two appropriately scaled and normalized features (e.g., peak CPU utilization and start time of user-invoked processes on a host in 24 hours). Note that it is rare to find time series datasets that correspond to simple and known analytical distributions. More likely than not, this method will be suitable in anomaly detection problems with the time dimension excluded. We will synthesize this dataset by sampling a Gaussian distribution and then including a 0.01 ratio of outliers in the mixture:

import numpy as np

num_dimensions = 2
num_samples = 1000
outlier_ratio = 0.01
num_inliers = int(num_samples * (1-outlier_ratio))
num_outliers = num_samples - num_inliers

# Generate the normally distributed inliers
x = np.random.randn(num_inliers, num_dimensions)

# Add outliers sampled from a random uniform distribution
x_rand = np.random.uniform(low=-10, high=10, size=(num_outliers, num_dimensions))
x = np.r_[x, x_rand]

# Generate labels, 1 for inliers and −1 for outliers
labels = np.ones(num_samples, dtype=int)
labels[-num_outliers:] = −1

Plotting this dataset in a scatter plot (see Figure 3-10), we see that the outliers are visibly separated from the central mode cluster:

import matplotlib.pyplot as plt

plt.plot(x[:num_inliers,0], x[:num_inliers,1], 'wo', label='inliers')
plt.plot(x[-num_outliers:,0], x[-num_outliers:,1], 'ko', label='outliers')
plt.xlim(-11,11)
plt.ylim(-11,11)
plt.legend(numpoints=1)
plt.show()

Figure 3-10. Scatter plot of synthetic dataset with inlier/outlier ground truth labels

Elliptical envelope fitting does seem like a suitable option for anomaly detection given that the data looks normally distributed (as illustrated in Figure 3-10). We use the convenient sklearn.covariance.EllipticEnvelope class in the following analysis:²⁶

from sklearn.covariance import EllipticEnvelope

classifier = EllipticEnvelope(contamination=outlier_ratio)
classifier.fit(x)
y_pred = classifier.predict(x)
num_errors = sum(y_pred != labels)
print('Number of errors: {}'.format(num_errors))

> Number of errors: 0

This method performs very well on this dataset, but that is not surprising at all given the regularity of the distribution. In this example, we know the accurate value for outlier_ratio to be 0.01 because we created the dataset synthetically. This is an important parameter because it informs the classifier of the proportion of outliers it should look for. For realistic scenarios in which the outlier ratio is not known, you should make your best guess for the initial value based on your knowledge of the problem. Thereafter, you can iteratively tune the outlier_ratio upward if you are not detecting some outliers that the algorithm should have found, or tune it downward if there is a problem with false positives.

Let’s take a closer look at the decision boundary formed by this classifier, which is illustrated in Figure 3-11.

Figure 3-11. Decision boundary for elliptic envelope fitting on Gaussian synthetic data

The center mode is shaded in gray, demarcated by an elliptical decision boundary. Any points lying beyond the decision boundary of the ellipse are considered to be outliers.

We need to keep in mind that this method’s effectiveness varies across different data distributions. Let’s consider at a dataset that does not fit a regular Gaussian distribution (see Figure 3-12).

Figure 3-12. Decision boundary for elliptic envelope fitting on non-Gaussian synthetic data

There are now eight misclassifications: four outliers are now classified as inliers, and four inliers that fall just outside the decision boundary are flagged as outliers.

Applying this method in a streaming anomaly detection system is straightforward. By periodically fitting the elliptical envelope to new data, you will have a constantly updating decision boundary with which to classify incoming data points. To remove effects of drift and a continually expanding decision boundary over time, it is a good idea to retire data points after a certain amount of time. However, to ensure that seasonal and cyclical effects are covered, this sliding window of fresh data points needs to be wide enough to encapsulate information about daily or weekly patterns.

The EllipticEnvelope() function in sklearn is located in the sklearn.covariance module. The covariance of features in a dataset refers to the joint variability of the features. In other words, it is a measure of the magnitude and direction of the effect that a change in one feature has on another. The covariance is a characteristic of a dataset that we can use to describe distributions, and in turn to detect outliers that do not fit within the described distributions. Covariance estimators can be used to empirically estimate the covariance of a dataset given some training data, which is exactly how covariance-based fitting for anomaly detection works.

Robust covariance estimates²⁷ such as the Minimum Covariance Determinant (MCD) estimator will minimize the influence of training data outliers on the fitted model. We can measure the quality of a fitted model by the distance between outliers and the model’s distribution, using a distance function such as Mahalanobis distance. Compared with nonrobust estimates such as the Maximum Likelihood Estimator (MLE), MCD is able to discriminate between outliers and inliers, generating a better fit that results in inliers having small distances and outliers having large distances to the central mode of the fitted model.

The elliptic envelope fitting method makes use of robust covariance estimators to attain covariance estimates that model the distribution of the regular training data, and then classifies points that do not meet these estimates as anomalies. We’ve seen that elliptic envelope fitting works reasonably well for a two-dimensional contaminated dataset with a known Gaussian distribution, but not so well on a non-Gaussian dataset. You can apply this technique to higher-dimensional datasets as well, but your mileage may vary—elliptic envelopes work better on datasets with low dimensionality. When using it on time series data, you might find it useful in some scenarios to remove time from the feature set and just fit the model to a subset of other features. In this case, however, note that you will not be able to capture an anomaly that is statistically regular relative to the aggregate distribution, but in fact is anomalous relative to the time it appeared. For example, if some anomalous data points from the middle of the night have features that exhibit values that are not out of the ordinary for a midday data point, but are highly anomalous for nighttime measurements, the outliers might not be detected if you omit the time dimension.

Unsupervised Machine Learning Algorithms

We now turn to a class of solutions to the anomaly detection problem that arise from modifications of typical supervised machine learning models. Supervised machine learning classifiers are typically used to solve problems that involve two or more classes. However, when used for anomaly detection, the modifications of these algorithms give them characteristics of unsupervised learning. In this section we look at a couple such algorithms.

One-class support vector machines

We can use a one-class SVM to detect anomalies by fitting the SVM with data belonging to only a single class. This data (which is assumed to contain no anomalies) is used to train the model, creating a decision boundary that can be used to classify future incoming data points. There is no robustness mechanism built into standard one-class SVM implementations, which means that the model training is less resilient to outliers in the dataset. As such, this method is more suitable for novelty detection than outlier detection; that is, the training data should ideally be thoroughly cleaned and contain no anomalies.

Where the one-class SVM method pulls apart from the pack is in dealing with non-Gaussian or multimodal distributions (i.e., when there is more than one “center” of regular inliers), as well as high-dimensional datasets. We will apply the one-class SVM classifier to the second dataset we used in the preceding section. Note that this dataset is not ideal for this method, because outliers comprise one percent of the data, but let’s see how much the resulting model is affected by the presence of contaminants:²⁸

from sklearn import svm

classifier = svm.OneClassSVM(nu=0.99 * outlier_ratio + 0.01,
                             kernel="rbf",
                             gamma=0.1)
classifier.fit(x)
y_pred = classifier.predict(x)
num_errors = sum(y_pred != labels)
print('Number of errors: {}'.format(num_errors))

Let’s examine the custom parameters that we specified in the creation of the svm.OneClassSVM object. Note that these parameters are dependent on datasets and usage scenarios; in general, you should always have a good understanding of all tunable parameters offered by a classifier before you use it. To deal with a small proportion of outliers in the data, we set the nu parameter to be roughly equivalent to the outlier ratio. According to the sklearn documentation, this parameter controls the “upper bound on the fraction of training errors and the lower bound of the fraction of support vectors.” In other words, it represents the acceptable range of errors generated by the model that can be caused by stray outliers, allowing the model some flexibility to prevent overfitting the model to outliers in the training set.

The kernel is selected by visually inspecting the dataset’s distribution. Each cluster in the bimodal distribution has Gaussian characteristics, which suggests that the radial basis function (RBF) kernel would be a good fit given that the values of both the Gaussian function and the RBF decrease exponentially as points move radially further away from the center.

The gamma parameter is used to tune the RBF kernel. This parameter defines how much influence any one training sample has on the resulting model. Its default value is 0.5. Smaller values of gamma would result in a “smoother” decision boundary, which might not be able to adequately capture the shape of the dataset. Larger values might result in overfitting. We chose a smaller value of gamma in this case to prevent overfitting to outliers that are close to the decision boundary.

Inspecting the resulting model, we see that the one-class SVM is able to fit this strongly bimodal dataset quite well, generating two mode clusters of inliers, as demonstrated in Figure 3-13. There are 16 misclassifications, so the presence of outliers in the training data did have some effect on the resulting model.

Figure 3-13. Decision boundary for one-class SVM on bimodal synthetic data—trained using both outliers and inliers

Let’s retrain the model on purely the inliers and see if it does any better. Figure 3-14 presents the result.

Figure 3-14. Decision boundary for one-class SVM on bimodal synthetic data—trained using only inliers

Indeed, as can be observed from Figure 3-14, there now are only three classification errors.

One-class SVMs offer a more flexible method for fitting a learned distribution to your dataset than robust covariance estimation. If you are thinking of using one as the engine for your anomaly detection system, however, you need to pay special attention to potential outliers that might slip past detection and cause the gradual degradation of the model’s accuracy.

Isolation forests

Random forest classifiers have a reputation for working well as anomaly detection engines in high-dimensional datasets. Random forests are algorithmic trees, and stream classification on tree data structures is much more efficient compared to models that involve many cluster or distance function computations. The number of feature value comparisons required to classify an incoming data point is the height of the tree (vertical distance between the root node and the terminating leaf node). This makes it very suitable for real-time anomaly detection on time series data.

The sklearn.ensemble.IsolationForest class helps determine the anomaly score of a sample using the Isolation Forest algorithm. This algorithm trains a model by iterating through data points in the training set, randomly selecting a feature and randomly selecting a split value between the maximum and minimum values of that feature (across the entire dataset). The algorithm operates in the context of anomaly detection by computing the number of splits required to isolate a single sample; that is, how many times we need to perform splits on features in the dataset before we end up with a region that contains only the single target sample. The intuition behind this method is that inliers have more feature value similarities, which requires them to go through more splits to be isolated. Outliers, on the other hand, should be easier to isolate with a small number of splits because they will likely have some feature value differences that distinguish them from inliers. By measuring the “path length” of recursive splits from the root of the tree, we have a metric with which we can attribute an anomaly score to data points. Anomalous data points should have shorter path lengths than regular data points. In the sklearn implementation, the threshold for points to be considered anomalous is defined by the contamination ratio. With a contamination ratio of 0.01, the shortest 1% of paths will be considered anomalies.

Let’s see this method in action by applying the Isolation Forest algorithm on the non-Gaussian contaminated dataset we saw in earlier sections (see Figure 3-15):²⁹

from sklearn.ensemble import IsolationForest

rng = np.random.RandomState(99)

classifier = IsolationForest(max_samples=num_samples,
                             contamination=outlier_ratio,
                             random_state=rng)
classifier.fit(x)
y_pred = classifier.predict(x)
num_errors = sum(y_pred != labels)
print('Number of errors: {}'.format(num_errors))

> Number of errors: 8

Figure 3-15. Decision boundary for isolation forest on synthetic non-Gaussian data

Using isolation forests in streaming time series anomaly detection is very similar to using one-class SVMs or robust covariance estimations. The anomaly detector simply maintains a tree of isolation forest splits and updates the model with new incoming points (as long as they are not deemed anomalies) in newly isolated segments of the feature space.

It is important to note that even though testing/classification is efficient, initially training the model is often more resource and time intensive than other methods of anomaly detection discussed earlier. On very low-dimensional data, using isolation forests for anomaly detection might not be suitable, because the small number of features on which we can perform splits can limit the effectiveness of the algorithm.

Density-Based Methods

Clustering methods such as the k-means algorithm are known for their use in unsupervised classification and regression. We can use similar density-based methods in the context of anomaly detection to identify outliers. Density-based methods are well suited for high-dimensional datasets, which can be difficult to deal with using the other classes of anomaly detection methods. Several different density-based methods have been adapted for use in anomaly detection. The main idea behind all of them is to form a cluster representation of the training data, under the hypothesis that outliers or anomalies will be located in low-density regions of this cluster representation. This approach has the convenient property of being resilient to outliers in the training data because such instances will likely also be found in low-density regions.

Even though the k-nearest neighbors (k-NN) algorithm is not a clustering algorithm, it is commonly considered a density-based method and is actually quite a popular way to measure the probability that a data point is an outlier. In essence, the algorithm can estimate the local sample density of a point by measuring its distance to the k^th nearest neighbor. You can also use k-means clustering for anomaly detection in a similar way, using distances between the point and centroids as a measure of sample density. k-NN has the potential to scale well to large datasets by using k-d trees (k-dimensional trees), which can greatly improve computation times for smaller-dimensional datasets.³⁰ In this section, we will focus on a method called the local outlier factor (LOF), which is a classic density-based machine learning method for anomaly detection.

Local outlier factor

The LOF is an anomaly score that you can generate using the scikit-learn class sklearn.neighbors.LocalOutlierFactor. Similar to the aforementioned k-NN and k-means anomaly detection methods, LOF classifies anomalies using local density around a sample. The local density of a data point refers to the concentration of other points in the immediate surrounding region, where the size of this region can be defined either by a fixed distance threshold or by the closest n neighboring points. LOF measures the isolation of a single data point with respect to its closest n neighbors. Data points with a significantly lower local density than that of their closest n neighbors are considered to be anomalies. Let’s run an example on a similar non-Gaussian, contaminated dataset once again:³¹

from sklearn.neighbors import LocalOutlierFactor

classifier = LocalOutlierFactor(n_neighbors=100)
y_pred = classifier.fit_predict(x)

Z = classifier._decision_function(np.c_[xx.ravel(), yy.ravel()])
Z = Z.reshape(xx.shape)

num_errors = sum(y_pred != labels)
print('Number of errors: {}'.format(num_errors))

> Number of errors: 9

Figure 3-16 presents the result.

Figure 3-16. Decision boundary for local outlier factor on bimodal synthetic distribution

As we can observe from Figure 3-16, LOF works very well even when there is contamination in the training set, and it is not very strongly affected by dimensionality of the data. As long as the dataset maintains the property that outliers have a weaker local density than their neighbors in a majority of the training features, LOF can find clusters of inliers well. Because of the algorithm’s approach, it is able to distinguish between outliers in datasets with varying cluster densities. For instance, a point in a sparse cluster might have a higher distance to its nearest neighbors than another point in a denser cluster (in another area of the same dataset), but because density comparisons are made only with local neighbors, each cluster will have different distance conditions for what constitutes an outlier. Lastly, LOF’s nonparametric nature means that it can easily be generalized across multiple different dimensions as long as the data is numerical and continuous.

In Summary

Having analyzed the five categories of anomaly detection algorithms, it should be clear that there is no lack of machine learning methods applicable to this classic data mining problem. Selecting which algorithm to use can sometimes be daunting and can take a few iterations of trial and error. However, using our guidelines and hints for which classes of algorithms work better for the nature of the data you have and for the problem you are solving, you will be in a much better position to take advantage of the power of machine learning to detect anomalies.

Optimizing for Explainability

As mentioned earlier, the semantic gap of alert explainability is one of the biggest stumbling blocks of anomaly detectors using machine learning. Many practical machine learning applications value explanations of results. However, true explainability of machine learning is an area of research that hasn’t yet seen many definitive answers.

Simple machine learning classifiers, and even non–machine learning classification engines, are quite transparent in their predictions. For example, a linear regression model on a two-dimensional dataset generates very explainable results, but lacks the ability to learn more complex and nuanced features. More complex machine learning models, such as neural networks, random forest classifiers, and ensemble techniques, can fit real-world data better, but they are very black-box—the decision-making processes are completely opaque to an external observer. However, there are ways to approach the problem that can alleviate the concern that machine learning predictions are difficult to explain, proving that explainability is not in fact at odds with accuracy.³² Having an external system generate simple, human-readable explanations for the decisions made by a black-box classifier satisfies the conditions of result explainability,³³ even if the explanations do not describe the actual decision-making conditions of the machine learning system. This external system analyzes any output from the machine learning system and performs context-aware data analysis to generate the most probable reasons for why the alert was raised.

Maintainability of Anomaly Detection Systems

The longevity and usefulness of machine learning systems is dictated not by accuracy or efficacy, but by the understandability, maintainability, and ease of configuration of the software. Designing a modular system that allows for swapping out, removing, and reimplementing subcomponents is crucial in environments that are in constant flux. The nature of data constantly changes, and a well-performing machine learning model today might no longer be suitable half a year down the road. If the anomaly detection system is designed and implemented on the assumption that elliptic envelope fitting is to be used, it will be difficult to swap the algorithm out for, say, isolation forests in the future. Flexible configuration of both system and algorithm parameters is important for the same reason. If tuning model parameters requires recompiling binaries, the system is not configurable enough.

Integrating Human Feedback

Having a feedback loop in your anomaly detection system can make for a formidable adaptive system. If security analysts are able to report false positives and false negatives directly to a system that adjusts model parameters based on this feedback, the maintainability and flexibility of the system can be vastly elevated. In untrusted environments, however, directly integrating human feedback into the model training can have negative effects.

Mitigating Adversarial Effects

As mentioned earlier, in a hostile environment your machine learning security systems almost certainly will be attacked. Attackers of machine learning systems generally use one of two classes of methods to achieve their goals. If the system continually learns from input data and instantaneous feedback labels provided by users (online learning model), attackers can poison the model by injecting intentionally misleading chaff traffic to skew the decision boundaries of classifiers. Attackers can also evade classifiers with adversarial examples that are specially crafted to trick specific models and implementations. It is important to put specific processes in place to explicitly prevent these threat vectors from penetrating your system. In particular, designing a system that blindly takes user input to update the model is risky. In an online learning model, inspecting any input that will be converted to model training data is important for detecting attempts at poisoning the system. Using robust statistics that are resilient to poisoning and probing attempts is another way of slowing down the attacker. Maintaining test sets and heuristics that periodically test for abnormalities in the input data, model decision boundary, or classification results can also be useful. We further explore adversarial problems and their solutions in Chapter 8.

Defining Malware Classification

Before we begin tearing apart binaries, let’s ground the discussion with some definitions. Malware classification groups distinct malware samples together based on common properties. We can classify malware in many different ways, depending on the purpose of the task. For instance, a security operations team might group malware by severity and function in order to effectively triage the risk that it poses to an organization. Security response teams might group malware by potential scope of damage and entry vector in order to devise remediation and mitigation strategies. Malware researchers might categorize malware by origin and authorship in order to understand its genealogy and purpose.

For general-purpose malware analysis, industry practice is to group samples by family—a term used by malware analysts that allows for tracking authorship, correlating information, and identifying new variants of newly found malware. Malware samples of the same family can have similar code, capabilities, authorship, functions, purposes, and/or origins. A famous example of a malware family is Conficker, a worm targeting the Microsoft Windows operating system. Even though there are many variations of the Conficker worm, each with different code, authors, and behavior, certain characteristics of the worms cause them to be attributed to the same malware family, indicating that they have likely evolved from a previously known ancestor. For example, all of the Conficker worms exploit Windows OS vulnerabilities and engage in dictionary attacks to crack the password of the administrator account, thereafter installing covert software on the exploited host to engage in botnet activity.

Differences between malware samples within the same family can originate from different compilers used to compile the source code, or from sections of code added and/or removed to modify the functionality of the malware itself. Malware samples that evolve over time in response to changing detection or mitigation strategies often also exhibit similarities between the older and newer versions, allowing analysts to trace the evolution of a family of malware. Nevertheless, malware family attribution is a notoriously difficult task that can have different results depending on the classification definitions and methods used by the analyst.

Malware classification can also be generalized to include the classification of nonmalicious binaries. This type of classification is used to determine whether a piece of software is malicious. Given an arbitrary binary, we want to know the likelihood that we are able to trust it and execute it in a trusted environment. This is a core objective of antivirus software and is an especially critical task for computer security practitioners, because this knowledge can help to prevent the spread of malware within an organization. Traditionally, this task is driven by signature matching: given a trove of properties and behavior of previously seen malware, new incoming binaries can be compared against this dataset to determine whether it matches something seen before.

The signature-matching method performs well so long as malware authors fail to significantly change properties and behavior of the malware to avoid detection, and the selected properties and behavior have a good balance of signal stability (so all malware samples belonging to this family exhibit this signal) and distinctiveness (so benign binaries will not exhibit properties or behaviors that cause them to be wrongly classified as malware). However, malware authors have a strong incentive to continuously alter the properties and behavior of their software to avoid detection.

Metamorphic or polymorphic⁶ viruses and worms employ static and dynamic obfuscation techniques to change characteristics of their code, behavior, and properties used in the signature generation algorithms of malware identification engines. This level of sophistication in malware used to be rare but has become more common due to its continued success in thwarting syntactic signature malware engines. Syntactic signature engines continue to chase the ever-narrowing set of static signals that malware authors neglect to obfuscate or fundamentally cannot change.

Malware: Behind the Scenes

To generate a descriptive dataset for classifying malware, we need to understand how malware works. This in turn requires some discussion of the malware economy, common types of malware, and general software execution processes in modern computing environments.

Modern code execution processes

We now examine how general classes of modern programs are written and executed, and consider how one might inspect binaries and executing programs to understand their inner workings without any access to the written code.

Note

The following discussion describes the code execution process for a large class of common computer programs, which applies to many modern programming languages and execution platforms. It is by no means a comprehensive or representative depiction of how all kinds of programs are executed. The vast and diverse ecosystem of programming environments and system runtimes results in a range of subtle to obtuse differences in how code executes in different environments. Nevertheless, many of the concepts we discuss are generalizable and parallels can often be drawn with other types of code execution processes.

In general, there are two types of code execution: compiled execution and interpreted execution. In compiled execution, the written code is translated into native machine instructions by a series of conversion steps⁷ (often referred to as the software build process). These machine instructions are packaged into binaries, which can then be executed directly by hardware. In interpreted execution implementations, the written code (sometimes referred to as a script) is translated into an intermediate format which is then fed into an interpreter for program execution. The interpreter is in charge of enacting the program’s instructions on the hardware it is running on. The intermediate format varies between different implementations, but is most commonly a form of bytecode (binary machine instructions) that will be executed on a virtual machine.⁸ Some implementations are a hybrid of compiled and interpreted executions, often using a process called just-in-time (JIT) compilation, in which interpreter bytecode is compiled into native machine instructions in real time.⁹

Figure 4-1 depicts common code execution processes for some modern software implementations.

Figure 4-1. Code execution and program analysis flowchart

Let’s take a closer look at the elements in Figure 4-1:

• The rectangular boxes represent the program in its various states of existence.

• The ellipses represent software conversion steps that translate the program from one state to another.

• The solid arrows between nodes represent the progression of the code from its human-written state to its eventual execution on the hardware.

• The gray box contains some tools that reverse engineers can use to inspect the static or dynamic state of a binary or running program (as indicated by the dashed arrows), providing valuable points of visibility into the code execution process.

Compiled code execution

As an example, we’ll take a piece of C code that performs some simple arithmetic and go step-by-step through the build process for a compiled implementation. Referring to Figure 4-1, we follow the path of a program from the initial “Source Code (Compiled Execution)” state. Here is the code we want to build, saved in a file named add.c:¹⁰

#include <stdio.h>

int main() {
    // Adds 1 to variable x
    int x = 3;
    printf("x + 1 = %d", x + 1);
    return 0;
}

1. The first step of the build process is a small but important one: preprocessing (omitted in Figure 4-1). In C, lines starting with the # character are interpreted by the preprocessor as preprocessor directives. The preprocessor simply iterates through the code and treats these directives as macros, preparing the code for compilation by inserting contents of included libraries and removing code comments, amongst other similar actions it performs. To inspect the results of the preprocessing stage, you can run the following command:

   > cc -E add.c
   
   [above lines omitted for brevity]
   extern void funlockfile (FILE *__stream)
                             __attribute__ ((__nothrow__ , __leaf__));
   # 942 "/usr/include/stdio.h" 3 4

   # 2 "add.c" 2
   
   # 3 "add.c"
   int main() {
   
       int x = 3;
       printf("x + 1 = %d", x + 1);
       return 0;
   }

   Note that the output of the preprocessing stage contains many lines of code that weren’t in the original add.c file. The preprocessor has replaced the #include <stdio.h> line with some contents from the standard C library stdio.h. Also note that the inline comment in the original code no longer shows up in this output.

2. The next step of the build process is compilation. Here, the compiler translates the preprocessed code into assembly code. The assembly code generated is specific to the target processor architecture, since it contains instructions that the CPU has to understand and execute. The assembly instructions generated by the compiler must be part of the instruction set understood by the underlying processor. To inspect the output of the C compiler by saving the assembly code to a file, you can run the following:

   > cc -S add.c

   This is the assembly code generated on a specific version of the GCC¹¹ (GNU Compiler Collection) C compiler, targeted at a 64-bit Linux system (x86_64-linux-gnu):

   > cat add.s
   
       .file   "add.c"
       .section    .rodata
   .LC0:
       .string "x + 1 = %d"
       .text
       .globl  main
       .type   main, @function
   main:
   .LFB0:
       .cfi_startproc
       pushq   %rbp
       .cfi_def_cfa_offset 16
       .cfi_offset 6, −16
       movq    %rsp, %rbp
       .cfi_def_cfa_register 6
       subq    $16, %rsp
       movl    $3, −4(%rbp)
       movl    −4(%rbp), %eax
       addl    $1, %eax
       movl    %eax, %esi
       movl    $.LC0, %edi
       movl    $0, %eax
       call    printf
       movl    $0, %eax
       leave
       .cfi_def_cfa 7, 8
       ret
       .cfi_endproc
   .LFE0:
       .size   main, .-main
       .ident  "GCC: (Ubuntu 5.4.0-6ubuntu1~16.04.4) 5.4.0 20160609"
       .section    .note.GNU-stack,"",@progbits

   This output will seem unintelligible unless you are familiar with assembly code (in this case, x64 assembly code). However, with some knowledge of assembly it is possible to gather quite a complete picture of what the program is doing solely based on this code. Looking at the two lines in bold in example output, addl ... and call printf, it is pretty easy to guess that the program is doing an addition and then invoking the print function. Most of the other lines just make up the plumbing—moving values in and out of CPU registers and memory locations where other functions can access them. Nevertheless, analyzing assembly code is an involved topic, and we will not go into further detail here.¹²

3. After the assembly code is generated, it is then up to the assembler to translate this into object code (machine code). The output of the assembler is a set of machine instructions that the target processor will directly execute:

   > cc -c add.c

   This command creates the object file add.o. The contents of this file are in binary format and are difficult to decipher, but let’s inspect it anyway. We can do this using tools such as hexdump and od. The hexdump utility, by default, displays the contents of the target file in hexadecimal format. The first column of the output indicates the offset of the file (in hexadecimal) where you can find the corresponding content:

   > hexdump add.o
   
   0000000 457f 464c 0102 0001 0000 0000 0000 0000
   0000010 0001 003e 0001 0000 0000 0000 0000 0000
   0000020 0000 0000 0000 0000 02b8 0000 0000 0000
   0000030 0000 0000 0040 0000 0000 0040 000d 000a
   0000040 4855 e589 8348 10ec 45c7 03fc 0000 8b00
   0000050 fc45 c083 8901 bfc6 0000 0000 00b8 0000
   0000060 e800 0000 0000 00b8 0000 c900 78c3 2b20
                [omitted for brevity]
   00005c0 0000 0000 0000 0000 0000 0000 0000 0000
   00005d0 01f0 0000 0000 0000 0013 0000 0000 0000
   00005e0 0000 0000 0000 0000 0001 0000 0000 0000
   *
   00005f8

   The od (which stands for octal dump) utility dumps contents of files in octal and other formats. Its output might be slightly more readable, unless you are a hex-reading wizard:

   > od -c add.o
   
   ...0000 177   E   L   F 002 001 001  \0  \0  \0  \0  \0  \0  \0  \0  \0
   ...0020 001  \0   >  \0 001  \0  \0  \0  \0  \0  \0  \0  \0  \0  \0  \0
   ...0040  \0  \0  \0  \0  \0  \0  \0  \0 270 002  \0  \0  \0  \0  \0  \0
   ...0060  \0  \0  \0  \0   @  \0  \0  \0  \0  \0   @  \0  \r  \0  \n  \0
   ...0100   U   H 211 345   H 203 354 020 307   E 374 003  \0  \0  \0 213
                            [omitted for brevity]
   ...2700  \0  \0  \0  \0  \0  \0  \0  \0  \0  \0  \0  \0  \0  \0  \0  \0
   ...2720 360 001  \0  \0  \0  \0  \0  \0 023  \0  \0  \0  \0  \0  \0  \0
   ...2740  \0  \0  \0  \0  \0  \0  \0  \0 001  \0  \0  \0  \0  \0  \0  \0
   ...2760  \0  \0  \0  \0  \0  \0  \0  \0
   ...2770

   This allows us to directly make out some of the structure of the binary file. For instance, notice that around the beginning of the file lie the characters E, L, and F. The assembler produced an ELF file (Executable and Linkable Format, specifically ELF64), and every ELF file begins with a header indicating some properties of the file, including what type of file it is. A utility such as readelf can help us to parse out all of the information embedded within this header:

   > readelf -h add.o
   
   ELF Header:
     Magic:   7f 45 4c 46 02 01 01 00 00 00 00 00 00 00 00 00
     Class:                             ELF64
     Data:                              2's complement, little endian
     Version:                           1 (current)
     OS/ABI:                            UNIX - System V
     ABI Version:                       0
     Type:                              REL (Relocatable file)
     Machine:                           Advanced Micro Devices X86-64
     Version:                           0x1
     Entry point address:               0x0
     Start of program headers:          0 (bytes into file)
     Start of section headers:          696 (bytes into file)
     Flags:                             0x0
     Size of this header:               64 (bytes)
     Size of program headers:           0 (bytes)
     Number of program headers:         0
     Size of section headers:           64 (bytes)
     Number of section headers:         13
     Section header string table index: 10

4. At this stage, let’s try to execute the object file generated by the assembler:¹³

   > chmod u+x add.o
   > ./add.o
   bash: ./add.o: cannot execute binary file: Exec format error

   Why does the Exec format error show up? The object code generated by the assembler is missing some crucial pieces of the program that are required for execution. Furthermore, sections of the program are not arranged properly, so library and program functions cannot be successfully invoked. Linking, the final stage of the build process, will fix these issues. In this case, the linker will insert the object code for the printf library function into the binary. Let’s invoke cc to generate the final executable binary (specifying the name of the output as add; otherwise, cc will use the default name of a.out), and then run the program:

   > cc -o add add.c
   > chmod u+x add
   > ./add
   
   x + 1 = 4
   

   This concludes the build process for a simple program in C, from code to execution.

In the preceding example, the stdio.h external library was statically linked into the binary, which means that it was compiled together with the rest of the code in a single package. Some languages and implementations allow for the dynamic inclusion of external libraries, which means that library components referenced in the code are not included in the compiled binary. Upon execution, the loader is invoked, scanning the program for references to dynamically linked libraries (or shared libraries,¹⁴ with extensions .so, .dll, etc.) and then resolving these references by locating the libraries on the system. We do not go into further detail on dynamic library loading mechanisms here.

Interpreted code execution

As an example of interpreted language implementations, we will dissect the typical Python script execution process. Note that there are several different implementations of Python with different code execution processes. In this example, we look at CPython,¹⁵ the standard and original implementation of Python written in C. In particular, we are using Python 3.5.2 on Ubuntu 16.04. Again referring to Figure 4-1, we follow the path of a program from the initial “Source Code (Interpreted Execution)” state:¹⁶

class AddOne():
    def __init__(self, start):
        self.val = start
    def res(self):
        return self.val + 1

def main():
    x = AddOne(3)
    print('3 + 1 = {} '.format(x.res()))

if __name__ == '__main__':
    main()

1. We begin with this Python source code saved in a file, add.py. Running the script by passing it as an argument to the Python interpreter yields the expected result:

   > python add.py
   3 + 1 = 4

   Admittedly, this is quite a convoluted way to add two numbers, but this example gives us a chance to explore the Python build mechanism. Internally, this human-written Python code is compiled into an intermediate format known as bytecode, a platform-independent representation of the program. We can see compiled Python modules (.pyc files¹⁷) created if the script imports external modules and is able to write to the target directory.¹⁸ In this case, no external modules were imported, so no .pyc files were created. For the sake of inspecting the build process, we can force the creation of this file by using the py_compile module:

   > python -m py_compile add.py

   This creates the .pyc file, which contains the compiled bytecode for our program. In Python 3.5.2, the compiled Python file is created as pycache/add.cpython-35.pyc. We then can inspect the contents of this binary file by removing the header and unmarshaling the file into a types.CodeType structure:

   import marshal
   import types
   
   # Convert a big-endian 32-bit byte array to a long
   def to_long(s):
       return s[0] + (s[1] << 8) + (s[2] << 16) + (s[3] << 24)
   
   # Print out hierarchy of code names and line numbers
   def inspect_code(code, indent='    '):
       print('{}{}(line:{})'.format(indent,
           code.co_name, code.co_firstlineno))
       for c in code.co_consts:
           if isinstance(c, types.CodeType):
               inspect_code(c, indent + '    ')
   
   f = open('__pycache__/add.cpython-35.pyc', 'rb')
   
   # Read .pyc file header
   magic = f.read(4)
   print('magic: {}'.format(magic.hex()))
   mod_time = to_long(f.read(4))
   print('mod_time: {}'.format(mod_time))
   
   # Only Python >=3.3 .pyc files contain the source_size header 
   source_size = to_long(f.read(4))
   print('source_size: {}'.format(source_size))
   
   print('\ncode:')
   code = marshal.load(f)
   inspect_code(code)
   
   f.close()

   

   Python .pyc files from version 3.2 and below have a header containing two 32-bit big-endian numbers followed by the marshaled code object. In versions 3.3 and above, a new 32-bit field that encodes the size of the source file is included in the header, as well (increasing the size of the header from 8 bytes to 12 bytes in Python 3.3 compared to earlier version).

   Executing this script yields the following results:

   magic: 160d0d0a
   mod_time: 1493965574
   source_size: 231
   
   code:
       <module>(line:1)
           AddOne(line:1)
               __init__(line:2)
               res(line:4)
           main(line:7)

   There is more information encoded in the CodeType object that we are not displaying, but this shows the general structure of the bytecode binary.

2. This bytecode is executed by the Python virtual machine runtime. Note that this bytecode is not binary machine code, but rather Python-specific opcodes that are interpreted by this virtual machine, which then translates the code into machine instructions. Using the dis.disassemble() function to disassemble the code object we created previously, we get the following:

   > import dis
   > dis.disassemble(code)
   
     1           0 LOAD_BUILD_CLASS
                 1 LOAD_CONST               0 (< code object AddOne at
                                               0x7f78741f7930, file
                                               "add.py", line 1> )
                 4 LOAD_CONST               1 ('AddOne')
                 7 MAKE_FUNCTION            0
                10 LOAD_CONST               1 ('AddOne')
                13 CALL_FUNCTION            2 (2 positional,
                                               0 keyword pair)
                16 STORE_NAME               0 (AddOne)
   
     7          19 LOAD_CONST               2 (< code object main at
                                               0x7f78741f79c0, file
                                               "add.py", line 7> )
                22 LOAD_CONST               3 ('main')
                25 MAKE_FUNCTION            0
                28 STORE_NAME               1 (main)
   
    11          31 LOAD_NAME                2 (__name__)
                34 LOAD_CONST               4 ('__main__')
                37 COMPARE_OP               2 (==)
                40 POP_JUMP_IF_FALSE       50
   
    12          43 LOAD_NAME                1 (main)
                46 CALL_FUNCTION            0 (0 positional,
                                               0 keyword pair)
                49 POP_TOP
           > >    50 LOAD_CONST               5 (None)
                53 RETURN_VALUE

   You can also obtain the output shown in the previous two steps by invoking the Python trace module on the command line via python -m trace add.py.

   You can immediately see the similarities between this output and the x86 assembly code we discussed earlier. The Python virtual machine reads in this bytecode and converts it to machine code,¹⁹ which executes on the target architecture, thereby completing the code execution process.

The code execution process for interpreted languages is shorter because there is no required build or compilation step: you can run code immediately after you write it. You cannot run Python bytecode directly on target hardware, as it relies on interpretation by the Python virtual machine and further translation to machine code. This process results in some inefficiencies and performance losses compared to “lower-level” languages like C. Nevertheless, note that the Python file code execution described earlier involves some degree of compilation, so the Python virtual machine doesn’t need to reanalyze and reparse each source statement repeatedly through the course of the program. Running Python in interactive shell mode is closer to the model of a pure interpreted language implementation, because each line is analyzed and parsed at execution time.

With access to human-written source code, we can easily parse specific properties and intentions of a piece of software that allow us to accurately classify it by family and function. However, because we don’t often have access to the code, we must resort to more indirect means to extract information about the program. With an understanding of modern code execution processes, we can now begin to look at some different ways to approach static and runtime analysis of malware. Code traverses a well-defined path in its journey from authorship to execution. Intercepting it at any point along the path can reveal a great deal of information about the program.