Isolation forest

Deep dive

How it works

In: n×d numeric matrix, t, ψ · Out: score in (0,1) per row · State: t small trees (≤ 2ψ−1 nodes each) · Cost: train O(t·ψ·log ψ), score O(t·log ψ) per point

1. For each of t trees, draw ψ rows without replacement (ψ = min(256, n)). Height limit l = ⌈log₂ ψ⌉.
2. Grow: at a node with more than one row and depth < l, pick a random feature and a split value uniform between its min and max in the node; x[f] < v goes left. A constant feature makes a leaf.
3. A leaf keeps only its row count size.
4. Path length in one tree: edges walked plus c(size), the depth the unbuilt subtree would have added.
5. c(n) = 2H(n−1) − 2(n−1)/n (average unsuccessful BST search), H(i) ≈ ln i + 0.5772156649; c(2) = 1, c(1) = c(0) = 0.
6. Score s(x) = 2^(−E[h(x)]/c(ψ)), E over the t trees. Short mean path gives s near 1.
7. Threshold: sklearn's contamination='auto' uses the paper's cut (s = 0.5, offset −0.5); a float (0, 0.5] sets the cut at that fraction of the training scores.

Why subsample: a big sample lets anomalies clump and hide each other (masking) or sit among normals (swamping). Small samples keep them sparse and trees tiny.

Code

2-D blob from summed uniforms plus 5 planted outliers; own xorshift RNG.

use std::time::Instant;

struct Rng(u64);
impl Rng {
    fn f64(&mut self) -> f64 { // xorshift64, uniform in [0,1)
        self.0 ^= self.0 << 13; self.0 ^= self.0 >> 7; self.0 ^= self.0 << 17;
        (self.0 >> 11) as f64 / (1u64 << 53) as f64
    }
    fn below(&mut self, n: usize) -> usize { (self.f64() * n as f64) as usize }
}

const D: usize = 2;
type P = [f64; D];

enum Node { Leaf(usize), Split { f: usize, v: f64, l: usize, r: usize } }

/// Average path length of an unsuccessful BST search over n points: c(n) = 2H(n-1) - 2(n-1)/n.
fn c(n: usize) -> f64 {
    match n {
        0 | 1 => 0.0,
        2 => 1.0,
        _ => 2.0 * ((n as f64 - 1.0).ln() + 0.577_215_664_901_532_9) - 2.0 * (n as f64 - 1.0) / n as f64,
    }
}

/// Build one tree over `pts` (reordered in place); returns the node index (root = last pushed).
fn build(nodes: &mut Vec<Node>, pts: &mut [P], depth: usize, limit: usize, rng: &mut Rng) -> usize {
    if depth >= limit || pts.len() <= 1 { nodes.push(Node::Leaf(pts.len())); return nodes.len() - 1; }
    let f = rng.below(D);
    let (lo, hi) = pts.iter().fold((f64::MAX, f64::MIN), |(a, b), p| (a.min(p[f]), b.max(p[f])));
    if lo == hi { nodes.push(Node::Leaf(pts.len())); return nodes.len() - 1; }
    let v = lo + rng.f64() * (hi - lo); // uniform split value between min and max
    let mut i = 0;                      // partition: < v on the left
    for j in 0..pts.len() { if pts[j][f] < v { pts.swap(i, j); i += 1; } }
    let (a, b) = pts.split_at_mut(i);
    let l = build(nodes, a, depth + 1, limit, rng);
    let r = build(nodes, b, depth + 1, limit, rng);
    nodes.push(Node::Split { f, v, l, r });
    nodes.len() - 1
}

struct Forest { trees: Vec<Vec<Node>>, psi: usize }

impl Forest {
    fn fit(data: &[P], t: usize, psi: usize, rng: &mut Rng) -> Forest {
        let limit = (psi as f64).log2().ceil() as usize;
        let trees = (0..t).map(|_| {
            let mut idx: Vec<usize> = Vec::with_capacity(psi); // sample psi distinct rows (psi << n)
            while idx.len() < psi { let i = rng.below(data.len()); if !idx.contains(&i) { idx.push(i); } }
            let mut sample: Vec<P> = idx.iter().map(|&i| data[i]).collect();
            let mut nodes = Vec::new();
            build(&mut nodes, &mut sample, 0, limit, rng);
            nodes
        }).collect();
        Forest { trees, psi }
    }
    fn path(nodes: &[Node], x: &P) -> f64 {
        let (mut n, mut h) = (nodes.len() - 1, 0.0);
        loop {
            match nodes[n] {
                Node::Leaf(size) => return h + c(size), // c(size) credits the unbuilt subtree
                Node::Split { f, v, l, r } => { n = if x[f] < v { l } else { r }; h += 1.0; }
            }
        }
    }
    /// s(x) = 2^(-E[h(x)] / c(psi)); near 1 = anomaly, below 0.5 = normal.
    fn score(&self, x: &P) -> f64 {
        let eh = self.trees.iter().map(|t| Self::path(t, x)).sum::<f64>() / self.trees.len() as f64;
        2f64.powf(-eh / c(self.psi))
    }
}

fn main() {
    let mut rng = Rng(0x9E37_79B9_7F4A_7C15);
    let blob = |n: usize, rng: &mut Rng| -> Vec<P> { // sum of 4 uniforms: bell-shaped, sd ~0.58
        (0..n).map(|_| { let mut p = [0.0; D]; for v in p.iter_mut() { *v = (0..4).map(|_| rng.f64()).sum::<f64>() - 2.0; } p }).collect()
    };
    let mut data = blob(10_000, &mut rng);
    let outliers: [P; 5] = [[6.0, 5.0], [-7.0, 4.0], [5.0, -6.0], [-5.0, -5.0], [8.0, 0.0]];
    data.extend_from_slice(&outliers);

    let t0 = Instant::now();
    let forest = Forest::fit(&data, 100, 256, &mut rng);
    let train = t0.elapsed();
    let t1 = Instant::now();
    let scores: Vec<f64> = data.iter().map(|p| forest.score(p)).collect();
    let score_t = t1.elapsed();

    let mut inl: Vec<f64> = scores[..10_000].to_vec();
    inl.sort_by(|a, b| a.partial_cmp(b).unwrap());
    let p99 = inl[(0.99 * inl.len() as f64) as usize];
    let omin = scores[10_000..].iter().cloned().fold(f64::MAX, f64::min);
    println!("train {:?}, score {:?} for {} points, inlier p99 {:.3}, weakest outlier {:.3}", train, score_t, data.len(), p99, omin);
    assert!(omin > p99);
    assert!(scores.iter().all(|&s| s > 0.0 && s < 1.0));
    assert!((c(256) - 10.2448).abs() < 1e-3); // 2(ln 255 + gamma) - 2*255/256
    let exact_h: f64 = (1..256).map(|i| 1.0 / i as f64).sum();
    assert!((c(256) - (2.0 * exact_h - 2.0 * 255.0 / 256.0)).abs() < 0.01);
}

Run on an Apple M4 Pro (rustc 1.99.1, --release, one thread): train 1.0 ms, scoring 10 005 points 30 ms. Asserts pass: weakest outlier 0.668, inlier 99th percentile 0.639. With 10⁶ inliers the outlier is 0.645 and the percentile 0.653: the check fails (see Pitfalls).

Complexity and constants

• Train: t·ψ·log ψ node visits, order of 2·10⁵ for the defaults. Measured here: 1.0 ms at n = 10⁴, 1.35 ms at n = 10⁶.
• Score: t·(≤ 9) hops. Measured 3.0 µs per point (n = 10⁴), 3.2 µs (n = 10⁶): about 30 ns per tree.
• Memory: 40 B per node here, ≤ 511 nodes per tree, so ≤ 2 MB for 100 trees (computed).

Efficiency

• Cost: score = t trees × ≤ 9 dependent loads and compares; measured 3 µs per point in 2-D on one M4 Pro core, about 330k points/s/core. Bound: dependent-load latency; the forest (≤ 2 MB, computed) stays in cache. Training is milliseconds; scoring dominates.
• Levers:
  • Trees and points are independent — near-linear speedup with cores (n_jobs in sklearn) — large scoring batches.
  • Fewer trees — cost linear in t (this page) — latency-bound serving.
  • Compact nodes (u16 children, f32 split, u8 feature: about 12 B instead of 40 B, computed) — more trees per cache — many models, edge devices.
• Naive → standard → state of the art:
  • Full data per tree, no depth limit · train O(t·n log n), memory O(t·n)
  • ψ = 256, depth ⌈log₂ ψ⌉ (sklearn defaults) · train about 1 ms, score about 3 µs/point (measured here)
  • Parallel, compact layout; random cut forest for streams, extended trees for artefacts · same order per point, adds incremental updates

Implementation notes

• Sample ψ rows without replacement; the code rejects duplicates, fine while ψ ≪ n (else Floyd's algorithm).
• Scores depend on ψ via c(ψ): compare only across equal ψ.
• Encode categoricals first. Splits are uniform per feature, so scaling is unnecessary, but a heavy-tailed feature's range dominates its cuts: clip or log-transform.

Pitfalls

• Score saturation: a point far outside the data range isolates at about the depth of one just outside, so the score hardly tells "far" from "very far". With 10⁶ inliers the planted outliers did not beat the inlier 99th percentile here (measured).
• A contamination fraction is a setting, not a finding. Choose the cut on a labelled sample or an alert budget (evaluation metrics, calibration).
• Many irrelevant features dilute the random feature choice; local anomalies beside a dense cluster look normal.
• Repeated anomalies (a fraud ring, a stuck sensor) can fill the subsample and mask each other.
• A fitted forest ages: refit on a schedule; watch drift.

Variants and alternatives

In production

scikit-learn IsolationForest, PyOD IForest, SageMaker Random Cut Forest (streaming cousin), the eif package.

Implementations

Sources

• Liu, Ting, Zhou, Isolation Forest, ICDM 2008. paper
• Liu, Ting, Zhou, Isolation-Based Anomaly Detection, ACM TKDD 2012. paper
• Hariri, Carrasco Kind, Brunner, Extended Isolation Forest, IEEE TKDE 2021. paper · code
• Guha, Mishra, Roy, Schrijvers, Robust Random Cut Forest Based Anomaly Detection on Streams, ICML 2016. paper · code
• scikit-learn, IsolationForest documentation (defaults, depth limit, contamination).
• Amazon SageMaker, Random Cut Forest and its hyperparameters.