Isolation Forest: Finding Outliers by Isolating Them (Part 2)

This is the second post in the series on model-based outlier detection. The first post discussed what outliers are, why z-scores fail through masking, and how robust statistics help. In this post I shall discuss Isolation Forest, one of the most widely used methods for detecting anomalies, and build it from scratch before using scikit-learn.

The key idea: anomalies are easy to isolate

Most methods profile what "normal" looks like and then flag whatever does not fit. Isolation Forest (Liu, Ting and Zhou, 2008) turns this around. Anomalies are few and different, so if we keep cutting the feature space at random, an anomaly is separated from the rest after only a few cuts, while a point deep inside a dense cluster needs many more.

The number of cuts needed to isolate a point is its path length. Short paths suggest anomalies and long paths suggest regular points. No distance, density or distributional assumption is needed.

Building an isolation tree

An isolation tree (iTree) is built on a sample of \(\psi\) points as follows.

1. Pick a feature \(q\) at random.
2. Pick a split value \(p\) uniformly at random between the minimum and maximum of that feature in the current node.
3. Send points with \(x_q < p\) to the left child and the rest to the right child.
4. Repeat on each child until a node has one point, all its points are identical, or the tree reaches the height limit \(l = \lceil \log_2 \psi \rceil\).

The height limit is deliberate. We only care about points that are isolated early, so there is no need to grow the tree fully.

From path length to anomaly score

For a point \(x\), let \(h(x)\) be the number of edges from the root to the leaf where \(x\) ends up. If that leaf still holds \(n\) unseparated points (because the height limit was reached), we add an adjustment \(c(n)\) for the splits the tree would still have needed. Averaging over all trees gives \(E[h(x)]\).

The adjustment is the average path length of an unsuccessful search in a binary search tree built on \(n\) points:

$$c(n) = 2H(n-1) - \frac{2(n-1)}{n}, \qquad H(i) \approx \ln(i) + \gamma$$

where \(\gamma \approx 0.5772\) is Euler's constant. It also serves as a normaliser, so that scores are comparable across sample sizes. The anomaly score is

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

The score behaves as follows.

• As \(E[h(x)] \to 0\), \(s \to 1\). The point is almost certainly an anomaly.
• When \(E[h(x)] = c(\psi)\), \(s = 0.5\). The point is average.
• As \(E[h(x)] \to \psi - 1\), \(s \to 0\). The point is very safely regular.

If every point scores close to 0.5, the data has no distinct anomalies. For the default \(\psi = 256\), \(c(\psi) \approx 10.24\).

Why subsampling helps

Each tree is built on a small random subsample (256 points by default in scikit-learn), not on the whole dataset. This sounds like it should hurt, but it helps. With the full data, anomalies sitting close to normal points (swamping) or clustered together (masking) become hard to isolate. A small sample thins out both effects. It also makes the method fast, since training cost depends on \(\psi\) and the number of trees rather than the dataset size, which is why Isolation Forest scales to large datasets.

Implementing it from scratch

The code below follows the description above. The three classes together are about 50 lines.

import numpy as np

def c(n):
    """Average path length of an unsuccessful search in a binary search tree."""
    if n <= 1:
        return 0.0
    if n == 2:
        return 1.0
    return 2.0 * (np.log(n - 1) + np.euler_gamma) - 2.0 * (n - 1) / n

class IsoTree:
    def __init__(self, X, height, max_height, rng):
        n = X.shape[0]
        self.size = n
        self.leaf = (height >= max_height) or (n <= 1)
        if not self.leaf:
            q = rng.randint(X.shape[1])                  # random feature
            lo, hi = X[:, q].min(), X[:, q].max()
            if lo == hi:                                 # cannot split further
                self.leaf = True
            else:
                self.q, self.p = q, rng.uniform(lo, hi)  # random split value
                mask = X[:, q] < self.p
                self.left = IsoTree(X[mask], height + 1, max_height, rng)
                self.right = IsoTree(X[~mask], height + 1, max_height, rng)

    def path_length(self, x, height=0):
        if self.leaf:
            return height + c(self.size)                 # adjustment for unexpanded leaf
        child = self.left if x[self.q] < self.p else self.right
        return child.path_length(x, height + 1)

class MyIsolationForest:
    def __init__(self, n_trees=100, psi=256, seed=0):
        self.n_trees, self.psi, self.seed = n_trees, psi, seed

    def fit(self, X):
        rng = np.random.RandomState(self.seed)
        self.psi_ = min(self.psi, len(X))
        max_h = int(np.ceil(np.log2(self.psi_)))         # height limit
        self.trees = [
            IsoTree(X[rng.choice(len(X), self.psi_, replace=False)], 0, max_h, rng)
            for _ in range(self.n_trees)
        ]
        return self

    def avg_path_length(self, X):
        return np.array([[t.path_length(x) for t in self.trees] for x in X]).mean(axis=1)

    def anomaly_score(self, X):
        return 2.0 ** (-self.avg_path_length(X) / c(self.psi_))

Comparing with scikit-learn

Now generate a dataset with two clusters of different densities plus 30 uniformly scattered outliers, and run both our forest and scikit-learn's IsolationForest on it. (Run the code above first, since this block reuses c and MyIsolationForest.)

from scipy.stats import spearmanr
from sklearn.datasets import make_blobs
from sklearn.ensemble import IsolationForest
from sklearn.preprocessing import StandardScaler

rng = np.random.RandomState(42)
X_in, _ = make_blobs(n_samples=[300, 300], centers=[[-2, -2], [3, 3]],
                     cluster_std=[0.5, 1.2], random_state=42)
X_out = rng.uniform(-6, 8, size=(30, 2))
X = StandardScaler().fit_transform(np.vstack([X_in, X_out]))
is_out = np.r_[np.zeros(len(X_in), bool), np.ones(len(X_out), bool)]

# 1. Our own forest
mine = MyIsolationForest(n_trees=100, psi=256, seed=0).fit(X)
eh = mine.avg_path_length(X)          # E[h(x)]
s_mine = mine.anomaly_score(X)        # s(x)
print("c(256) = %.2f" % c(256))
print("mean E[h]: regular = %.2f, planted outliers = %.2f"
      % (eh[~is_out].mean(), eh[is_out].mean()))
print("mean s(x): regular = %.3f, planted outliers = %.3f"
      % (s_mine[~is_out].mean(), s_mine[is_out].mean()))

# 2. scikit-learn
iso = IsolationForest(n_estimators=100, random_state=0).fit(X)
s_sk = -iso.score_samples(X)          # score_samples returns the NEGATIVE of s(x)
print("max_samples_ = %d, offset_ = %.1f" % (iso.max_samples_, iso.offset_))
print("Spearman correlation (ours vs sklearn): %.3f" % spearmanr(s_mine, s_sk).correlation)

# 3. Effect of the contamination parameter
def report(name, model):
    flagged = model.predict(X) == -1
    tp = (flagged & is_out).sum()
    fp = (flagged & ~is_out).sum()
    print(f"{name}: flagged={flagged.sum()}, caught={tp}/{is_out.sum()}, false alarms={fp}")

report("contamination='auto'", iso)
for cont in [0.02, 0.05, 0.10, 0.20]:
    report(f"contamination={cont}",
           IsolationForest(n_estimators=100, contamination=cont, random_state=0).fit(X))

# 4. Sensitivity to the random seed (contamination fixed at 0.05)
for seed in range(3):
    report(f"seed={seed}, contamination=0.05",
           IsolationForest(n_estimators=100, contamination=0.05, random_state=seed).fit(X))

The output is as follows.

• The average path length \(E[h(x)]\) is 12.16 for regular points and 6.72 for the planted outliers. The mean anomaly score \(s(x)\) is 0.441 for regular points and 0.641 for the outliers.
• scikit-learn uses 256 samples per tree by default, and with contamination='auto' its decision threshold corresponds to \(s = 0.5\) (offset_ is −0.5).
• The Spearman rank correlation between our scores and scikit-learn's is 0.975. The two forests are built from different random trees and differ in some implementation details, yet they rank the points almost identically.

One detail to remember: scikit-learn's score_samples returns the negative of \(s(x)\), so higher means more normal. This is why the code above uses -iso.score_samples(X). Likewise decision_function is positive for inliers and negative for outliers, and predict returns +1 or −1.

The contamination parameter matters a lot

The last loop in the code above varies the contamination parameter, which is the fraction of points you expect to be anomalies. The table shows what happened with 600 regular points and 30 planted outliers.

contamination	flagged	outliers caught (of 30)	false alarms
'auto' (threshold s = 0.5)	97	27	70
0.02	13	13	0
0.05	32	22	10
0.10	63	27	36
0.20	126	27	99

Three lessons follow from this.

• The default is not a safe choice. With 'auto', the broad cluster's fringe scores above 0.5, so 70 regular points are flagged. A threshold of 0.5 is not a universal rule for "anomalous."
• Raising contamination trades precision for recall. Going from 0.02 to 0.10 catches more outliers, but false alarms jump from 0 to 36.
• The seed matters only slightly. With contamination=0.05, three different seeds all flag 32 points, and catch 22, 23 and 22 of the outliers.

In practice, treat the score as a ranking, look at the distribution of scores, and choose the threshold from the cost of a missed anomaly versus a false alarm. Reviewing the top-ranked cases by hand is often the best calibration you can get without labels.

Visualising the result

import matplotlib.pyplot as plt

m5 = IsolationForest(n_estimators=100, contamination=0.05, random_state=0).fit(X)
flagged = m5.predict(X) == -1

xx, yy = np.meshgrid(np.linspace(X[:, 0].min() - 0.5, X[:, 0].max() + 0.5, 300),
                     np.linspace(X[:, 1].min() - 0.5, X[:, 1].max() + 0.5, 300))
Z = -m5.score_samples(np.c_[xx.ravel(), yy.ravel()]).reshape(xx.shape)

fig, axes = plt.subplots(1, 2, figsize=(12, 4.8))

ax = axes[0]
cs = ax.contourf(xx, yy, Z, levels=20, cmap="viridis_r", alpha=0.85)
plt.colorbar(cs, ax=ax, label="anomaly score s(x)")
ax.scatter(X[~flagged, 0], X[~flagged, 1], s=12, c="white",
           edgecolors="k", linewidths=0.4, label="not flagged")
ax.scatter(X[flagged, 0], X[flagged, 1], s=30, c="red",
           edgecolors="k", linewidths=0.5, label="flagged")
ax.set_title("Isolation Forest anomaly-score map")
ax.set_xlabel("x1 (standardised)")
ax.set_ylabel("x2 (standardised)")
ax.legend(loc="lower right", fontsize=8, framealpha=0.9)

ax = axes[1]
bins = np.linspace(eh.min(), eh.max(), 30)
ax.hist(eh[~is_out], bins=bins, alpha=0.7, density=True, label="regular points")
ax.hist(eh[is_out], bins=bins, alpha=0.7, density=True, color="tab:red",
        label="planted outliers")
ax.axvline(c(mine.psi_), color="k", ls="--", lw=1, label=r"$c(\psi)$ (score = 0.5)")
ax.set_title("Average path length E[h(x)] (our own forest)")
ax.set_xlabel("E[h(x)]")
ax.set_ylabel("density")
ax.legend(fontsize=8)

plt.tight_layout()
plt.show()

 

Figure 1: Left: anomaly-score map from scikit-learn (contamination = 0.05), with flagged points in red. Right: average path lengths from our own forest.

In the left panel, yellow-green regions have a low score (regular) and purple regions a high score (anomalous). In the right panel, most planted outliers have an average path length between about 4 and 9, while most regular points are above 10. The dashed line marks \(c(\psi)\), where the score is 0.5.

Notice the blocky, axis-aligned bands in the score map, especially in the corners. They come from the axis-parallel splits, which can give regions along the axes slightly different scores from regions in other directions. Variants such as the Extended Isolation Forest use randomly oriented splits to reduce this artefact.

Why were three outliers never caught?

In the table, no setting catches more than 27 of the 30 planted outliers. The reason is in the data. The outliers were drawn uniformly from a wide box, so a few of them landed inside the clusters, where they are not really different from the regular points. The code below checks this.

from sklearn.neighbors import NearestNeighbors

m20 = IsolationForest(n_estimators=100, contamination=0.20, random_state=0).fit(X)
missed = np.where(is_out & (m20.predict(X) == 1))[0]       # planted outliers not flagged

nn = NearestNeighbors(n_neighbors=11).fit(X[~is_out])
d_missed = nn.kneighbors(X[missed])[0][:, 10]              # distance to 10th regular neighbour
d_regular = nn.kneighbors(X[~is_out])[0][:, 10]
print("missed outliers:", missed)
print("their 10th-neighbour distances:", np.round(d_missed, 3))
print("regular points: median = %.3f, 90th percentile = %.3f"
      % (np.median(d_regular), np.percentile(d_regular, 90)))

The three missed points have a 10th-nearest-neighbour distance of 0.14 to 0.17, which is within the range of regular points (median 0.115, 90th percentile 0.248). They are not outliers in any useful sense, and no detector should flag them without producing many false alarms. This is a good reminder that "planted outlier" and "detectable outlier" are not the same thing.

Practical notes

• Feature scaling is not required. Each split is drawn within a feature's own range, so rescaling a feature does not change the trees. I standardised the data here only to make plotting convenient.
• Irrelevant features hurt. Random features dilute the splits, so drop noisy features or select them first.
• Local anomalies are a weakness. Isolation Forest is strongest at global anomalies. A point that is unusual only relative to its immediate neighbourhood, especially next to a dense cluster, may score as normal. That is where density-based methods do better.
• Discrete and categorical features need care. Splits are on numeric ranges, so encode categories thoughtfully.
• Novelty detection: you can fit on clean data and call predict or decision_function on new observations.

Conclusion

In this post I discussed how Isolation Forest isolates anomalies with random splits, how path length becomes an anomaly score, and how to build the method from scratch and compare it with scikit-learn. I also showed that the choice of threshold matters at least as much as the algorithm. In the next post I shall discuss the Local Outlier Factor, which compares each point's density with that of its neighbours and so handles clusters of different densities.

Reference: F. T. Liu, K. M. Ting and Z.-H. Zhou, "Isolation Forest", Proceedings of the 8th IEEE International Conference on Data Mining (ICDM), 2008.

The code in this post was tested with Python 3 and scikit-learn 1.8.0.

Outliers, Anomalies and Novelties: Where to Begin (Part 1)

This is the first post in a series on model-based outlier detection. In the coming posts I shall discuss Isolation Forest, Local Outlier Factor and One-Class SVM. Before that, I want to build the foundation: what exactly an outlier is, why the simple rules we learned first fail, and how scikit-learn thinks about the problem.

What is an outlier?

An outlier is an observation that deviates so much from the rest of the data that it seems to come from a different mechanism. Sometimes it is an error (a mistyped price). Sometimes it is the most valuable record in the dataset (a fraudulent transaction, a failing machine, a market crash day). The method does not decide which, so a domain expert should review whatever is flagged.

Three types are worth separating:

• Point anomalies: a single observation is unusual on its own, such as a transaction of ₹50 lakh in an account that normally spends ₹5,000.
• Contextual anomalies: a value is normal in general but unusual in context. A temperature of 35°C is fine in May but strange in January.
• Collective anomalies: each point looks fine, but together they are abnormal, such as a flat line in a sensor signal that should oscillate.

Why the z-score rule fails

The usual rule flags a point when \(|z_i| > 3\), where

$$z_i = \frac{x_i - \bar{x}}{s}$$

There are two problems.

Masking. The mean \(\bar{x}\) and standard deviation \(s\) are themselves computed from the data, including the outliers. Outliers pull the mean toward them and inflate \(s\), so they hide each other.

A hard ceiling. For a sample of size \(n\), no point can have a z-score larger than

$$|z_i| \le \frac{n-1}{\sqrt{n}}$$

(using \(s\) with the \(n-1\) divisor). For \(n = 9\) this is about 2.67, so the "3-sigma rule" can never flag anything in such a sample.

Consider this tiny dataset: 10, 11, 9, 10, 12, 11, 10, 45, 50. Two values are clearly far from the rest, yet their z-scores are only about 1.6 and 1.9.

A robust alternative: median and MAD

Replace the mean with the median and the standard deviation with the median absolute deviation, \(\text{MAD} = \text{median}(|x_i - \tilde{x}|)\). The modified z-score is

$$M_i = \frac{0.6745\,(x_i - \tilde{x})}{\text{MAD}}$$

A common rule (Iglewicz and Hoaglin) flags points with \(|M_i| > 3.5\). The median and MAD barely move when a few outliers are present, so masking largely disappears.

import numpy as np

x = np.array([10, 11, 9, 10, 12, 11, 10, 45, 50])

z = (x - x.mean()) / x.std(ddof=1)
med = np.median(x)
mad = np.median(np.abs(x - med))
mz = 0.6745 * (x - med) / mad

print("z-score :", np.round(z, 2))
print("modified:", np.round(mz, 2))

The z-scores of 45 and 50 come out at 1.60 and 1.91, while their modified z-scores are 22.93 and 26.31. The robust version separates them from the rest without any difficulty.

The multivariate case: Mahalanobis distance

With several features, an observation can be unremarkable on every single axis yet unusual in combination (for example, a very tall and very light person). The Mahalanobis distance accounts for the correlation between features:

$$D^2(\mathbf{x}) = (\mathbf{x}-{\mu})^\top {\Sigma}^{-1} (\mathbf{x}-{\mu})$$

If the data are approximately multivariate normal in \(p\) dimensions, \(D^2\) follows a \(\chi^2_p\) distribution, so a natural cutoff is the 97.5th percentile \(\chi^2_{p,\,0.975}\).

The same masking problem returns, because \({\mu}\) and \({\Sigma}\) are estimated from contaminated data. The fix is the Minimum Covariance Determinant (MCD) estimator. It searches for the subset of \(h\) observations whose sample covariance matrix has the smallest determinant, and computes the mean and covariance from that subset alone. It can tolerate up to nearly 50% contamination, and scikit-learn implements it in MinCovDet.

The code below plants a small cluster of outliers in a correlated dataset and compares the classical and robust estimates.

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import chi2
from sklearn.covariance import EmpiricalCovariance, MinCovDet

rng = np.random.RandomState(42)
n_in, n_out = 200, 30
cov = [[1.0, 0.8], [0.8, 1.0]]
X_in = rng.multivariate_normal([0, 0], cov, size=n_in)
X_out = rng.multivariate_normal([4, -3], [[0.3, 0], [0, 0.3]], size=n_out)
X = np.vstack([X_in, X_out])
is_out = np.r_[np.zeros(n_in, bool), np.ones(n_out, bool)]

classical = EmpiricalCovariance().fit(X)
robust = MinCovDet(random_state=0).fit(X)

d2_classical = classical.mahalanobis(X)   # squared distances
d2_robust = robust.mahalanobis(X)
cutoff = chi2.ppf(0.975, df=X.shape[1])

fig, axes = plt.subplots(1, 2, figsize=(11, 4.5), sharex=True, sharey=True)
for ax, (name, d2) in zip(axes, [("Classical", d2_classical),
                                 ("Robust (MCD)", d2_robust)]):
    flagged = d2 > cutoff
    tp = (flagged & is_out).sum()
    fp = (flagged & ~is_out).sum()
    print(f"{name}: flagged={flagged.sum()}, "
          f"true outliers caught={tp}/{n_out}, false alarms={fp}")
    ax.scatter(X[~flagged, 0], X[~flagged, 1], s=18, alpha=0.7, label="not flagged")
    ax.scatter(X[flagged, 0], X[flagged, 1], s=28, c="red", label="flagged")
    ax.set_title(name)
    ax.set_xlabel("x1")
axes[0].set_ylabel("x2")
axes[0].legend()
plt.tight_layout()
plt.show()

Running this gives the following counts.

• Classical: 13 points flagged, 11 of the 30 planted outliers caught, 2 false alarms.
• Robust (MCD): 36 points flagged, all 30 planted outliers caught, 6 false alarms.

The classical estimate misses most of the outliers, because they drag the mean and covariance toward themselves. The robust version catches all of them, at the cost of a few false alarms. This is the usual trade-off with a 97.5% cutoff, and it is worth keeping in mind when choosing a threshold. The exact counts can change slightly with the random seed and library version.

Figure 1: Points flagged (red) by the classical Mahalanobis distance (left) and the robust MCD-based distance (right).

 

Outlier detection versus novelty detection

scikit-learn distinguishes two tasks, and it is important to choose the right one:

• Outlier detection: the training data itself contains outliers, and the goal is to find them. The methods are unsupervised, and you typically call fit_predict.
• Novelty detection: the training data is assumed clean, and the goal is to decide whether a new observation looks like the training data. You call fit on clean data and then predict on new data.

In scikit-learn, predictions use the convention +1 for inliers and −1 for outliers. The estimators I shall cover in this series are IsolationForest, LocalOutlierFactor (which needs novelty=True for the novelty setting) and OneClassSVM (designed for novelty detection and sensitive to contamination in the training set).

Conclusion

In this post I discussed what outliers are, why mean-and-standard-deviation rules fail through masking, how median/MAD and the MCD-based Mahalanobis distance give robust alternatives, and how outlier detection differs from novelty detection. In the next post I shall discuss Isolation Forest, an elegant method that isolates anomalies with random splits and works well without any distributional assumption.

The code in this post was tested with Python 3 and scikit-learn 1.8.0.