Isolation Forest for Anomaly Detection in Python: Math, Algorithm, and Code Explained

Isolation Forest (or iForest) is an unsupervised machine learning algorithm designed specifically for anomaly detection. Unlike traditional outlier methods that construct computationally heavy density maps or profile normal behavior, Isolation Forest works on a fundamentally opposite premise: it directly exploits the two core properties of anomalies — they are few and different. Because outliers sit in sparse, peripheral regions of feature space, they can be separated from the rest of the dataset in far fewer random cuts than normal inliers.

1. Key Concepts & Notation Glossary

Before detailing the recursive tree partitioning, let us define the core variables and mathematical notation established by Liu, Ting, and Zhou (2008):

2. Mathematical Foundations & The Anomaly Score Derivation

An Isolation Tree recursively partitions a data subset $X = \{x_1, \dots, x_n\}$ of $m$ dimensions by randomly selecting a feature dimension $q \in \{1, \dots, m\}$ and a split point $p \sim \text{Uniform}(\min(X_{\cdot, q}), \max(X_{\cdot, q}))$. The dataset divides into two disjoint subsets:

$$X_{\text{left}} = \{x \in X \mid x_q < p\}, \quad X_{\text{right}} = \{x \in X \mid x_q \ge p\}$$

This recursive splitting continues until either: (1) $|X| \le 1$, (2) all feature values are identical, or (3) the tree hits the maximum depth limit $h_{\text{max}} = \lceil \log_2(\psi) \rceil$.

Why Binary Search Trees (BST) Provide the Normalization Factor c(n)

Because an Isolation Tree has the exact same structure as a random Binary Search Tree (BST), the average depth of an external terminating node in an iTree is mathematically equivalent to the average path length of an unsuccessful search in a BST constructed over $n$ random keys:

$$c(n) = 2H(n-1) - \frac{2(n-1)}{n}$$

Where $H(k)$ is the harmonic number, approximated by $H(k) = \ln(k) + \gamma$, with Euler-Mascheroni constant $\gamma \approx 0.5772156649$.

The Exponential Anomaly Score Equation

Using $c(n)$ to normalize the ensemble average path length $E[h(x)]$, the anomaly score $s(x, n)$ is defined as:

$$s(x, n) = 2^{-\frac{E[h(x)]}{c(n)}}$$

Because the exponent is $-\frac{E[h(x)]}{c(n)}$, the score exhibits three crucial operating regimes:

• When $E[h(x)] \to 0$: The point isolates near the root. $s(x, n) = 2^{-0} = 1.0$. The observation is definitely an anomaly.
• When $E[h(x)] \to c(n)$: The point isolates at the expected average depth of a random tree. $s(x, n) = 2^{-1} = 0.5$. The observation displays no distinct anomaly characteristics.
• When $E[h(x)] \to n - 1$: The point took maximum splits to isolate. $s(x, n) \to 2^{-\infty} \approx 0.0$. The observation is deeply buried in a dense cluster.

3. The Swamping & Masking Problem (Why Subsampling is Magical)

In most machine learning algorithms, larger datasets improve model accuracy. In Isolation Forest, the opposite is true: training trees on the full dataset severely degrades anomaly detection performance due to two fundamental phenomena:

Liu et al. proved that subsampling a small batch of observations per tree (standard $\psi = 256$) effectively separates anomaly clusters so they can be isolated in 1 or 2 cuts, while drastically reducing memory to $O(t \cdot \psi)$ and runtime to $O(t \cdot \psi \log \psi)$.

4. Step-by-Step Worked Numerical Trace

Let us trace the exact mathematics on a 1D dataset of 5 transaction values:

$$X = [10.0, 12.0, 11.0, 13.0, 500.0]$$

First, compute the normalization factor $c(n)$ for $n = 5$:

• $$H(5 - 1) = H(4) = \ln(4) + 0.5772156649 \approx 1.38629 + 0.57722 \approx 1.96351$$
• $$c(5) = 2(1.96351) - \frac{2(4)}{5} = 3.92702 - 1.60000 = 2.32702$$

Simulating Path Lengths Across Trees

When an isolation tree chooses a random split between $\min(X) = 10.0$ and $\max(X) = 500.0$, the random split threshold $p$ lands in the vast empty space between $13$ and $500$ with overwhelming probability ($487/490 \approx 99.4\%$ chance on the very first cut!).

• For point $500.0$: It is separated on cut 1 in virtually every tree: $E[h(500)] \approx 1.0$.
  $$s(500, 5) = 2^{-\frac{1.0}{2.32702}} = 2^{-0.4297} \approx 0.742 \quad (\text{High Anomaly Score})$$
• For point $11.0$: It is tightly flanked by $10, 12, 13$. It requires multiple binary cuts to strip away its neighbors: $E[h(11)] \approx 3.4$.
  $$s(11, 5) = 2^{-\frac{3.4}{2.32702}} = 2^{-1.4611} \approx 0.363 \quad (\text{Normal Inlier Score})$$

5. From-Scratch Python Implementation (Pure NumPy)

Here is the complete object-oriented from-scratch implementation of an Isolation Forest using pure NumPy, featuring recursive tree construction, depth-limit safeguards, and $c(n)$ normalization:

import numpy as np

class IsolationTreeNode:
    """A node in an isolation tree (iTree)."""
    def __init__(self, left=None, right=None, split_feature=None, split_value=None, size=0, is_leaf=False):
        self.left = left
        self.right = right
        self.split_feature = split_feature
        self.split_value = split_value
        self.size = size
        self.is_leaf = is_leaf


def c_factor(n):
    """
    Average path length of an unsuccessful search in a Binary Search Tree (BST)
    over n observations: c(n) = 2*H(n-1) - (2*(n-1)/n).
    """
    if n <= 1:
        return 0.0
    if n == 2:
        return 1.0
    euler_gamma = 0.5772156649
    harmonic_number = np.log(n - 1) + euler_gamma
    return 2.0 * harmonic_number - (2.0 * (n - 1) / n)


class IsolationTree:
    """Builds a single Isolation Tree by recursively partitioning feature space."""
    def __init__(self, max_depth=10):
        self.max_depth = max_depth
        self.root = None

    def fit(self, X, current_depth=0, rng=None):
        if rng is None:
            rng = np.random.default_rng()

        n_samples, n_features = X.shape

        # Base case: isolate 1 point, hit max depth, or all points identical
        if n_samples <= 1 or current_depth >= self.max_depth:
            return IsolationTreeNode(size=n_samples, is_leaf=True)

        # Randomly choose a feature
        feature_idx = rng.integers(0, n_features)
        feat_col = X[:, feature_idx]
        f_min, f_max = feat_col.min(), feat_col.max()

        if f_min == f_max:
            return IsolationTreeNode(size=n_samples, is_leaf=True)

        # Randomly choose split point between min and max
        split_val = rng.uniform(f_min, f_max)
        left_mask = feat_col < split_val
        right_mask = ~left_mask

        left_child = self.fit(X[left_mask], current_depth + 1, rng)
        right_child = self.fit(X[right_mask], current_depth + 1, rng)

        return IsolationTreeNode(
            left=left_child,
            right=right_child,
            split_feature=feature_idx,
            split_value=split_val,
            size=n_samples,
            is_leaf=False
        )

    def path_length(self, x, node=None, current_depth=0):
        """Measures the number of edges traversed to isolate observation x."""
        if node is None:
            node = self.root

        if node.is_leaf:
            # Add correction for un-split leaf clusters
            if node.size > 1:
                return current_depth + c_factor(node.size)
            return float(current_depth)

        if x[node.split_feature] < node.split_value:
            return self.path_length(x, node.left, current_depth + 1)
        else:
            return self.path_length(x, node.right, current_depth + 1)


class IsolationForestScratch:
    """Ensemble of Isolation Trees for unsupervised anomaly detection."""
    def __init__(self, n_estimators=100, max_samples=256, random_state=42):
        self.n_estimators = n_estimators
        self.max_samples = max_samples
        self.random_state = random_state
        self.trees = []
        self.subsample_size = None

    def fit(self, X):
        X = np.asarray(X, dtype=float)
        n_samples = len(X)
        self.subsample_size = min(self.max_samples, n_samples)
        # Standard maximum tree depth based on subsample size: ceil(log2(subsample_size))
        max_depth = int(np.ceil(np.log2(max(self.subsample_size, 2))))
        rng = np.random.default_rng(self.random_state)

        self.trees = []
        for _ in range(self.n_estimators):
            sub_indices = rng.choice(n_samples, size=self.subsample_size, replace=False)
            tree = IsolationTree(max_depth=max_depth)
            tree.root = tree.fit(X[sub_indices], current_depth=0, rng=rng)
            self.trees.append(tree)

        return self

    def anomaly_score(self, X):
        X = np.asarray(X, dtype=float)
        c_val = c_factor(self.subsample_size)
        if c_val == 0.0:
            return np.full(len(X), 0.5)

        # Average path length across all trees
        avg_paths = np.zeros(len(X))
        for x_idx, x in enumerate(X):
            avg_paths[x_idx] = np.mean([tree.path_length(x) for tree in self.trees])

        # Anomaly score s(x, n) = 2^(-E[h(x)] / c(n))
        return 2.0 ** (-avg_paths / c_val)

    def predict(self, X, threshold=0.6):
        """Predicts -1 for anomalies (score >= threshold) and 1 for normal."""
        scores = self.anomaly_score(X)
        return np.where(scores >= threshold, -1, 1)


# Verification test on the 5-point dataset
if __name__ == "__main__":
    X_demo = np.array([[10.0], [12.0], [11.0], [13.0], [500.0]])
    model = IsolationForestScratch(n_estimators=300, max_samples=5, random_state=42)
    model.fit(X_demo)
    scores = model.anomaly_score(X_demo)

    print("=== From-Scratch Isolation Forest Scores ===")
    for val, sc in zip(X_demo.ravel(), scores):
        status = "ANOMALY" if sc >= 0.6 else "Normal"
        print(f"Value: {val:5.1f} | Anomaly Score: {sc:.4f} | Prediction: {status}")

6. Production Scikit-Learn Pipeline & Sign Polarity Explained

In production, data scientists use sklearn.ensemble.IsolationForest. However, Scikit-Learn introduces a critical sign-convention twist that trips up many engineers:

import numpy as np
from sklearn.ensemble import IsolationForest

# 1. Real-world dataset: 5-point demo and synthetic transaction cluster
X_demo = np.array([[10.0], [12.0], [11.0], [13.0], [500.0]])

# 2. Fit Scikit-Learn IsolationForest
# contamination = expected fraction of anomalies (0.2 = 1 in 5 points)
clf = IsolationForest(
    n_estimators=100,
    max_samples=5,
    contamination=0.20,
    random_state=42
)
clf.fit(X_demo)

# 3. Model Predictions & Decision Function
# predict(): -1 indicates an anomaly, +1 indicates inlier
preds = clf.predict(X_demo)

# decision_function(): OPPOSITE sign convention!
# Negative values mean anomaly, positive values mean normal
raw_scores = clf.decision_function(X_demo)

# score_samples(): negative of anomaly score s(x, n)
offset_scores = clf.score_samples(X_demo)

print("=== Scikit-Learn Isolation Forest Evaluation ===")
for val, p, raw, off in zip(X_demo.ravel(), preds, raw_scores, offset_scores):
    label = "ANOMALY (-1)" if p == -1 else "Normal (+1)"
    print(f"Point: {val:5.1f} | Label: {label:14} | decision_function: {raw:+.4f} | score_samples: {off:.4f}")

7. Extended Isolation Forest (EIF) & Axis-Parallel Limitations

Standard Isolation Forest has one major geometric weakness: it strictly cuts feature space using axis-parallel hyperplanes ($x_q < p$). When continuous features exhibit strong diagonal correlations (for example, income correlated with home price):

• Axis-parallel splits create step-like rectangular partitions that leave artificial 'ghost' low-density corners.
• Real anomalies that violate the diagonal correlation but fall within individual coordinate ranges may require extra cuts to isolate.
• Extended Isolation Forest (EIF): Developed by Hariri et al. (2019), EIF selects random hyperplanes with arbitrary slope vectors ($w \cdot x < p$), cleanly isolating diagonal anomalies without grid-artifact distortion.

8. Anomaly Detection Comparison Matrix

9. Top 5 Pitfalls & Common Mistakes

1. Misinterpreting Scikit-Learn's Decision Function: Assuming negative scores are normal. In Scikit-Learn, negative is an anomaly (label $-1$); positive is an inlier (label $+1$).
2. Setting Subsample Size Equal to Dataset Size: Increasing `max_samples` to 100,000 causes swamping and masking, degrading accuracy. Keep `max_samples=256`.
3. Blindly Setting the Contamination Parameter: Setting `contamination` to an arbitrary number forces that exact percentage of data to be labeled $-1$. If ground truth is unknown, evaluate raw scores instead.
4. Expecting Local Anomaly Detection: An observation that is anomalous only within its small cluster will not be caught by iForest. Use Local Outlier Factor (LOF) instead.
5. Ignoring Categorical Preprocessing: Isolation Forest cannot perform numerical inequalities ($x_q < p$) on raw strings. Properly one-hot or target encode categorical columns first.

10. Summary & Key Takeaways

• Isolates Rather Than Profiles: Explores tree path length to isolate anomalies rather than modeling normal density.
• Short Path = Anomaly: Anomalies sit in sparse regions and get separated near the root of the tree in few cuts.
• BST Math Normalization: Scores are normalized against $c(n)$, yielding a probability-like metric where $s \to 1$ indicates strong anomalies.
• Subsampling is Essential: Small subsamples ($\psi = 256$) prevent swamping and masking while providing blazing-fast $O(n)$ inference.
• Polarity Awareness: Remember that Scikit-Learn's `decision_function` inverts signs (negative values indicate anomalies).