What it actually does
Every other detector in this section models what normal looks like and flags whatever falls outside. This one models the anomalies directly, from a single observation: anomalies are few and different, so a random split is more likely to separate one from the pack than to separate two ordinary points from each other. Cut the data with random splits until every point is alone, and record how many cuts each one needed. Outliers fall out early, in a handful of cuts; points inside a dense cluster survive many. The depth at which a point becomes alone IS the anomaly score, and no density was ever estimated to get it.
- Anomalies are few and different — the entire premise, and it fails when anomalies form a substantial or clustered fraction of the data
- Anomalies are distinguishable along individual axes, since splits are axis-aligned and cannot see a correlated combination
- Anomalies are global rather than local: a point sparse relative to its own cluster but dense relative to the whole is not detected
- Numeric features carry the signal; a random threshold on an arbitrary categorical encoding splits on nothing meaningful
Loss function
There is no training objective — nothing here is minimized, and the expression defines a scoring rule rather than something the fit pursues. Read it as a comparison: how deep did this point sit, against how deep a point in a random tree of this size would be expected to sit. Dividing by c(psi) is what makes the number interpretable rather than a raw depth, and it is why scores from forests built on different subsample sizes can be compared at all. The exponential base of two turns the ratio into a bounded score, which is convenience rather than theory.
- path length — the number of splits needed to isolate x in one tree
- average path length across the forest, which is what the ensemble is for
- expected path length in a random binary tree over psi points — the normalizer that makes scores comparable across subsample sizes
- subsample size per tree, conventionally 256 and deliberately small
- score in (0, 1): near 1 is an anomaly, near 0.5 is unremarkable, well below 0.5 means the point sits deep inside a cluster
How it is minimized
None — an ensemble of trees grown by uniformly random splits on small subsamples
- the split feature, drawn uniformly at random — no criterion is evaluated
- the threshold, drawn uniformly between the min and max of that feature within the node
- the range within the current node, not the global range — which is what makes splits adapt to local extent
- the path length, with a correction term added when a node is truncated at the depth limit
Nothing is optimized, and that is the design rather than a shortcut. Evaluating a criterion at every split is what makes tree building expensive; drawing the split uniformly makes it O(1), so a tree costs a single pass and the whole forest is linear in the data. Two choices look like limitations and are not. Trees are grown on small subsamples — 256 points is the standard, and using more usually makes the detector worse, because a large subsample lets normal points crowd around anomalies until they no longer isolate early. And trees are truncated at roughly log2(psi) depth, since a point still unresolved by then is already deep inside normal and its exact depth carries no information; the truncation is corrected analytically rather than by growing further.
max_samples (psi)
Subsample size per tree. The counterintuitive one: 256 is the standard default and increasing it typically degrades detection, because swamping and masking both worsen with more points per tree (256, and rarely worth changing)
n_estimators
Number of trees. Averages the path length, so more is monotonically better and only costs time (100 to 300; the score curve is flat well before 200)
contamination
Not a model parameter at all — it converts scores into labels by picking a quantile. Treat it as the alert budget it is (set from the operational alert rate, not estimated)
max_features
Features considered per tree. Below 1.0 it helps on wide data, where most features are irrelevant to any given anomaly
Nothing iterates; the score converges statistically as trees are added, and adding trees cannot make it worse. The failure modes are structural and worth naming precisely. Splits are axis-aligned, so an anomaly that is unusual only in a correlated combination of features — off the diagonal of an elongated cloud, but within range on every axis individually — is invisible, and the extended isolation forest exists to fix exactly this by cutting with random hyperplanes instead. Masking: a tight cluster of anomalies isolates no faster than a cluster of normal points, so a coordinated attack scores as normal. Swamping: normal points sitting near a dense anomaly region get dragged into early isolation and produce false positives. Both are why the subsample is small — fewer points per tree means less crowding, which is the opposite of the usual more-data intuition.
Training: O(t · psi · log psi) — independent of n entirely, because each tree only ever sees psi points. That is the headline property: the fit does not get slower as the dataset grows. Memory is O(t · psi), a few hundred tiny trees, so the model is kilobytes. Scoring: O(t · log psi) per point, a few hundred short tree walks with no arithmetic beyond a comparison, which makes it viable at streaming volumes where a distance- or density-based detector is not.
How, where, and why it is used
Time-Series Forecasting
Not applicableIt produces an outlier score, never a value for a future period; applying it to a series means finding anomalies within it, which is the anomaly-detection domain rather than this one.
Anomaly Detection
Primary fitFit an ensemble of random trees on small subsamples of the data and score each point by its average path length, normalized against the expected depth for that subsample size. There is no reference distribution to specify, no metric to choose, and no scaling step — the score comes from how quickly random cuts separate the point from everything else.
- High-volume fraud and abuse screening, where linear cost and a tiny model matter more than the last point of recall
- Infrastructure and application telemetry, scoring metrics that arrive faster than any density model could be refitted
- Data-quality screening across wide tables of mixed-scale numeric columns, where scaling every column correctly is not realistic
- The first detector to try on a new unlabelled problem, because it needs almost no decisions to be made first
It is the default for a specific and defensible reason: it is the only detector here whose cost does not grow with the dataset, and it demands almost nothing of the analyst — no metric, no scaling, no distributional assumption, and one hyperparameter that is fine at its default. Where it should not be used is equally clear. It finds global anomalies, so a point that is only unusual relative to its own local cluster is missed, and that is the case Local Outlier Factor exists for. It is blind to correlation, so an anomaly that is off-diagonal but in-range on every axis scores as normal — the extended variant with oblique splits is the honest fix. And a coordinated cluster of anomalies masks itself, which matters in exactly the adversarial settings where detection is most valuable.
- No scaling required, which is a genuine differentiator: the split threshold is drawn within each node’s own range, so features in different units cost nothing
- Drop obviously irrelevant columns or lower max_features — random feature choice means every irrelevant column dilutes the signal
- Encode categoricals as ordinals only if the ordering is meaningful; a random threshold on an arbitrary integer code splits on nothing
- Fit on a confirmed-clean window where one exists, though it is far more tolerant of contamination than a density-based detector
Precision@k against confirmed incidents and PR-AUC rather than ROC-AUC, at the class imbalance this domain always has. Hold out whole anomaly episodes rather than points from inside one, or a masked cluster will look like a set of individually detected anomalies. Set the threshold from the score distribution on clean data, and treat contamination as the alert budget it actually is.
- Raising max_samples in the belief that more data per tree helps — it usually makes detection worse, which is the least intuitive fact about the method
- Correlated anomalies: a point in-range on every axis but far off the joint diagonal scores as perfectly normal under axis-aligned cuts
- Local anomalies in a multi-density dataset, where a point sparse relative to its own cluster is dense relative to the whole and is missed
- Reading the score as a probability; it is a normalized depth ratio, and the 0.5 midpoint is a convention rather than a decision boundary
Optimization
Not applicableNothing is optimized — splits are drawn uniformly at random and no criterion is ever evaluated — and the output is an outlier score rather than a decision under constraints.
Risk & Fraud
Primary fitA first-line unsupervised screen over transaction, account, or claim features, running ahead of the supervised model so that a cheap linear-cost pass reduces what an expensive one has to see. It also covers the case supervised models structurally cannot: novel fraud patterns that have no labelled examples yet, because labels in this domain arrive weeks late if at all.
- Real-time transaction screening at volumes where a distance-based detector cannot meet latency
- Anti-money-laundering triage over account behaviour, where labels are scarce and delayed
- New-account and bot detection, where the attack pattern changes faster than a supervised model can be retrained
- Monitoring for population shift in the inputs to a deployed supervised model
The economics fit the domain exactly: constant-size model, linear scoring, no scaling to maintain, and no labels required — which means it keeps working during the window when a new fraud pattern exists and its labels do not. That window is precisely when a supervised model is blind. Against it, and this matters more here than anywhere else: sophisticated fraud is not a scattered set of odd points but a coordinated cluster, and a cluster masks itself under isolation. Treat it as a screen and a novelty tripwire, never as the fraud model.
- Aggregate to entity-level behavioural features over a window; raw transactions are too heavy-tailed for path length to say much
- Include velocity and ratio features, which is where genuinely anomalous behaviour usually shows up as an extreme value
- Keep the feature count moderate — random feature selection means each irrelevant column costs signal directly
- Refit on a rolling window so the notion of normal tracks the population, which is cheap because the fit does not depend on n
Precision at a fixed alert budget on an out-of-time split, and separately the fraction of confirmed fraud it discards — that second number is the only one that matters for a screen. Track it alongside the supervised model to see what each catches that the other does not.
- Coordinated fraud rings masking each other, which is the failure mode this domain is most likely to hit
- Legitimate high-value outliers — a genuine large transaction is anomalous and not fraud, so the score needs business context before it becomes a decision
- Treating the contamination parameter as an estimate of the fraud rate rather than as the alert budget it sets
Control & Operations
ViableScore sensor and equipment telemetry continuously against a forest refitted on a rolling window of recent normal operation. The model is small enough to sit on an edge device or a gateway, and the scoring cost per reading is a few hundred short tree walks.
- Predictive maintenance screening on vibration, temperature, and current signatures
- Process monitoring across many sensors of wildly different scales, where a distance metric would need per-sensor calibration
- Edge deployment on constrained hardware, where a kilobyte-scale model and comparison-only scoring are the binding constraints
Sensors arrive in incompatible units and drift constantly, and this is the detector that does not care — thresholds are drawn within each node’s own range, so no scaling has to be maintained as the units and ranges shift. Combined with a tiny model and a fit that ignores dataset size, it deploys where most detectors cannot. The reservation is the same structural one: equipment faults often appear as an unusual combination of individually normal readings, which is exactly what axis-aligned cuts cannot see, so a correlation-aware detector belongs alongside it.
- Window the signal into statistics — mean, variance, peak, spectral band energy — since a single raw reading rarely isolates
- Refit on a rolling recent window so gradual drift becomes the new normal instead of a growing stream of alerts
- Keep the feature set tight; on wide sensor arrays lower max_features rather than adding trees
Time-to-detection ahead of a confirmed fault, not just whether the fault was flagged eventually, and false alerts per machine per week against what the maintenance team can actually action.
- Gradual drift being absorbed by the rolling refit until a genuine slow degradation has become the baseline
- Multi-sensor faults that stay in range on each axis and are therefore invisible to axis-aligned splits
- Alert volume tuned by the contamination parameter rather than by what an operator can respond to
What it costs to run
The distinguishing property: cost is independent of dataset size. Each tree sees 256 points regardless of whether the dataset has ten thousand rows or a billion, so a fit is milliseconds and a rolling refit is essentially free. No other detector in this section can say that.
A few hundred short tree walks per point, comparisons only, no arithmetic and no distance computation. Microseconds, and the model is kilobytes — a few hundred trees over 256 points each — which is what makes edge and streaming deployment realistic.
As often as you like, and frequently the right answer is continuously on a rolling window. Because the fit ignores dataset size, keeping the notion of normal current costs almost nothing, which inverts the usual constraint where refresh frequency is limited by training cost.
- Track the score distribution rather than the alert count, since the alert count is fixed by the contamination threshold and will therefore look stable while the underlying scores move
- Watch the share of points scoring near 0.5, which is where a forest that has lost discriminating power piles up
- Compare a rolling-window model against a fixed-baseline model — divergence between them separates genuine drift from a rolling fit absorbing an ongoing problem
- Alert on features whose in-node ranges have shifted substantially, because the split thresholds are drawn from those ranges
- The score is a normalized depth ratio, not a probability; 0.5 is a convention rather than a decision boundary, and any expected-cost calculation on it is unfounded
- contamination sets a quantile threshold, so alert volume is fixed by configuration rather than discovered — raising it does not make the model more sensitive, it just alerts more
- A rolling refit will quietly absorb a slowly worsening condition into the baseline; keep a fixed-baseline model alongside it if slow degradation matters
- Scores are only comparable across models with the same subsample size, since c(psi) is the normalizer — changing max_samples silently rescales every threshold
- Determinism requires pinning the seed; with random splits, two fits on identical data give different scores, and a threshold tuned on one does not transfer exactly to the other
Pros and cons, in context
Training cost is independent of dataset size
Each tree sees a fixed 256-point subsample, so the fit does not slow down as data accumulates. Decisive at streaming and web scale, and irrelevant on a dataset small enough that a Gaussian mixture would fit instantly anyway.
No distance metric, no scaling, no distributional assumption
Thresholds are drawn within each node’s own range, so mixed units cost nothing and there is no scaler to persist or drift. This is the practical reason it is so often the first thing tried on a new unlabelled problem.
The model is tiny and scoring is comparison-only
A few hundred trees over 256 points each is kilobytes, and a score is a few hundred short walks with no arithmetic — which is what makes edge and in-database deployment realistic where a density model is not.
Models anomalies directly rather than modelling normality and inverting it
Nothing has to be assumed about the shape of normal, which is exactly where parametric detectors fail on multi-modal data. The cost is that it only sees the kind of anomaly its premise describes: few, and separable on an axis.
Blind to anomalies defined by correlation between features
Axis-aligned cuts cannot isolate a point that is in range on every axis but far off the joint structure — a very common shape for a real anomaly. The extended isolation forest with oblique splits is the fix, and it is not the default anywhere.
Finds global anomalies, not local ones
In a dataset with clusters of different densities, a point that is sparse relative to its own cluster still isolates slowly relative to the whole. That case is precisely what Local Outlier Factor is for, and the two are complements rather than alternatives.
Clustered anomalies mask each other
A group of similar anomalies isolates no faster than a group of normal points, so coordinated behaviour scores as normal. Most damaging in adversarial settings — which is where detection matters most.
The score is uncalibrated and the threshold is a configuration choice
contamination picks a quantile; it is an alert budget, not an estimate of how many anomalies exist. Teams routinely read it as the latter and then conclude the model is over- or under-detecting when they have simply set a dial.
Intuitive, idiomatic, optimized
The same algorithm in three languages, each progressing through Kent Beck's three stages. Every optimization names the lever it pulls and the trade-off it accepts — cost profiles are illustrative, not measured benchmarks.
The math, transcribed literally. No dependencies.
"""Isolation forest - random splits, path lengths, the score. Transcribed.
No criterion is ever evaluated. The split feature and threshold are drawn
uniformly, which is what makes tree construction O(1) per node and the whole
method linear. What is measured is only how deep a point had to go before it
was alone.
"""
import math
import random
def average_path_length(n):
"""c(psi) = 2H(n-1) - 2(n-1)/n, the expected depth in a random binary tree.
Dividing by this is what makes a raw depth into a comparable score, and
what lets forests built on different subsample sizes be compared at all.
"""
if n <= 1:
return 0.0
if n == 2:
return 1.0
harmonic = math.log(n - 1.0) + 0.5772156649015329 # Euler-Mascheroni
return 2.0 * harmonic - 2.0 * (n - 1.0) / n
def build_tree(X, rows, depth, max_depth, rng):
"""Split uniformly at random until isolated or the depth limit is hit."""
if depth >= max_depth or len(rows) <= 1:
# A node truncated at the limit still holds several points; the
# analytic correction estimates how much deeper they would have gone.
return {"leaf": True, "size": len(rows)}
d = len(X[0])
feature = rng.randrange(d)
# Range within THIS node, not the global range - which is what makes the
# splits adapt to local extent as the tree descends.
lowest = min(X[i][feature] for i in rows)
highest = max(X[i][feature] for i in rows)
if lowest == highest:
return {"leaf": True, "size": len(rows)}
threshold = rng.uniform(lowest, highest)
left = [i for i in rows if X[i][feature] <= threshold]
right = [i for i in rows if X[i][feature] > threshold]
return {
"leaf": False,
"feature": feature,
"threshold": threshold,
"left": build_tree(X, left, depth + 1, max_depth, rng),
"right": build_tree(X, right, depth + 1, max_depth, rng),
}
def fit(X, n_trees=100, subsample=256, seed=0):
"""Each tree sees only `subsample` points - and more is usually worse."""
rng = random.Random(seed)
n = len(X)
psi = min(subsample, n)
max_depth = int(math.ceil(math.log2(psi))) if psi > 1 else 1
forest = []
for _ in range(n_trees):
rows = rng.sample(range(n), psi)
forest.append(build_tree(X, rows, 0, max_depth, rng))
return forest, psi
def path_length(tree, x, depth=0):
if tree["leaf"]:
# Add the expected remaining depth for the points still in this node.
return depth + average_path_length(tree["size"])
if x[tree["feature"]] <= tree["threshold"]:
return path_length(tree["left"], x, depth + 1)
return path_length(tree["right"], x, depth + 1)
def score(forest, psi, x):
"""s(x, psi) = 2^(-E[h(x)] / c(psi)); near 1 is anomalous."""
total = 0.0
for tree in forest:
total += path_length(tree, x)
expected = total / len(forest)
return 2.0 ** (-expected / average_path_length(psi))O(t * psi * log psi) to fit and O(t * log psi) per query, in interpreter loops with a Python list rebuilt at every node.
How this language actually wants it written.
"""Isolation forest - typed, flat trees, vectorized subsampling."""
from dataclasses import dataclass
import numpy as np
from numpy.typing import NDArray
Matrix = NDArray[np.float64]
Vector = NDArray[np.float64]
EULER_MASCHERONI = 0.5772156649015329
def average_path_length(n: NDArray | int) -> NDArray | float:
"""c(psi) = 2H(n-1) - 2(n-1)/n, vectorized over node sizes.
This is the normalizer that turns a raw depth into a comparable score, and
also the analytic correction applied at a truncated leaf.
"""
sizes = np.asarray(n, dtype=np.float64)
harmonic = np.log(np.maximum(sizes - 1.0, 1e-12)) + EULER_MASCHERONI
corrected = 2.0 * harmonic - 2.0 * (sizes - 1.0) / np.maximum(sizes, 1.0)
return np.where(sizes <= 1.0, 0.0, corrected)
@dataclass(frozen=True)
class IsolationTree:
"""Parallel flat arrays rather than nested dicts.
A tree is at most 2*psi - 1 nodes and psi is 256, so the whole tree fits in
a few kilobytes of contiguous memory and a traversal never chases a pointer.
"""
feature: NDArray # int32; negative marks a leaf
threshold: Vector
left: NDArray # int32
right: NDArray # int32
leaf_size: NDArray # int32; drives the truncation correction
def path_length(self, x: Vector) -> float:
cursor = 0
depth = 0
while self.feature[cursor] >= 0:
cursor = (
self.left[cursor]
if x[self.feature[cursor]] <= self.threshold[cursor]
else self.right[cursor]
)
depth += 1
return depth + float(average_path_length(self.leaf_size[cursor]))
@dataclass(frozen=True)
class IsolationForest:
trees: list[IsolationTree]
subsample: int
def score_samples(self, X: Matrix) -> Vector:
"""Near 1 is anomalous, near 0.5 is unremarkable. NOT a probability."""
if X.ndim != 2:
raise ValueError(f"X must be 2-D, got shape {X.shape}")
expected = np.empty(X.shape[0], dtype=np.float64)
for index, row in enumerate(X):
expected[index] = np.mean([tree.path_length(row) for tree in self.trees])
return 2.0 ** (-expected / float(average_path_length(self.subsample)))
def _grow(X: Matrix, rows: NDArray, max_depth: int, rng: np.random.Generator) -> IsolationTree:
feature: list[int] = []
threshold: list[float] = []
left: list[int] = []
right: list[int] = []
leaf_size: list[int] = []
def build(node_rows: NDArray, depth: int) -> int:
index = len(feature)
feature.append(-1)
threshold.append(0.0)
left.append(-1)
right.append(-1)
leaf_size.append(node_rows.size)
if depth >= max_depth or node_rows.size <= 1:
return index
chosen = int(rng.integers(X.shape[1]))
column = X[node_rows, chosen]
lowest, highest = float(column.min()), float(column.max())
if lowest == highest:
return index
# Range within THIS node, so splits adapt to local extent.
cut = float(rng.uniform(lowest, highest))
mask = column <= cut
feature[index] = chosen
threshold[index] = cut
leaf_size[index] = 0
left[index] = build(node_rows[mask], depth + 1)
right[index] = build(node_rows[~mask], depth + 1)
return index
build(rows, 0)
return IsolationTree(
feature=np.asarray(feature, dtype=np.int32),
threshold=np.asarray(threshold, dtype=np.float64),
left=np.asarray(left, dtype=np.int32),
right=np.asarray(right, dtype=np.int32),
leaf_size=np.asarray(leaf_size, dtype=np.int32),
)
def fit(X: Matrix, n_trees: int = 100, subsample: int = 256, seed: int = 0) -> IsolationForest:
"""Fit a forest. Raises ValueError on malformed input.
subsample defaults to 256 and raising it usually makes detection WORSE:
more points per tree lets normal observations crowd around anomalies until
they no longer isolate early. That is swamping, and it is the single least
intuitive property of the method.
"""
if X.ndim != 2:
raise ValueError(f"X must be 2-D, got shape {X.shape}")
if n_trees < 1:
raise ValueError(f"n_trees must be at least 1, got {n_trees}")
if subsample < 2:
raise ValueError(f"subsample must be at least 2, got {subsample}")
rng = np.random.default_rng(seed)
psi = min(subsample, X.shape[0])
max_depth = max(1, int(np.ceil(np.log2(psi))))
trees = [
_grow(X, rng.choice(X.shape[0], size=psi, replace=False), max_depth, rng)
for _ in range(n_trees)
]
return IsolationForest(trees=trees, subsample=psi)What changed: Trees become parallel flat arrays instead of nested dicts — a tree is at most 2·psi−1 nodes with psi = 256, so the whole structure is a few kilobytes of contiguous memory and a traversal never chases a pointer or performs a dictionary lookup. The path-length normalizer becomes vectorized so it can be applied to whole arrays of leaf sizes rather than one at a time. Around that, the idiomatic changes: frozen dataclasses, validated input, and the subsample contract documented where it will actually be read, since raising max_samples is the most common way this model is made worse by someone trying to improve it.
O(t * psi * log psi) to fit, independent of n; O(t * log psi) per query over contiguous arrays.
Faster, and the justification for each change.
"""Isolation forest - one packed forest, level-synchronous batch scoring."""
import numpy as np
from numpy.typing import NDArray
Matrix = NDArray[np.float64]
Vector = NDArray[np.float64]
EULER_MASCHERONI = 0.5772156649015329
class PackedIsolationForest:
"""Every tree concatenated into one set of arrays.
A single tree is tiny, so the cost of scoring is dominated by interpreter
overhead rather than arithmetic: one Python-level call per tree per point.
Packing all t trees into one contiguous int32 block and traversing a whole
batch level by level removes that entirely - the interpreter runs once per
tree per LEVEL, not once per tree per point.
"""
def __init__(self, feature: NDArray, threshold: Vector, left: NDArray,
right: NDArray, leaf_size: NDArray, roots: NDArray, subsample: int) -> None:
# One contiguous block per field, single dtype: the traversal gathers
# from these arrays on every level and they stay cache-resident.
self._feature = np.ascontiguousarray(feature, dtype=np.int32)
self._threshold = np.ascontiguousarray(threshold, dtype=np.float64)
self._left = np.ascontiguousarray(left, dtype=np.int32)
self._right = np.ascontiguousarray(right, dtype=np.int32)
self._roots = np.ascontiguousarray(roots, dtype=np.int32)
self._subsample = subsample
# Precompute the truncation correction per node, so scoring never
# evaluates the harmonic expansion.
sizes = np.asarray(leaf_size, dtype=np.float64)
harmonic = np.log(np.maximum(sizes - 1.0, 1e-12)) + EULER_MASCHERONI
self._correction = np.where(
sizes <= 1.0, 0.0, 2.0 * harmonic - 2.0 * (sizes - 1.0) / np.maximum(sizes, 1.0)
)
self._normalizer = float(
2.0 * (np.log(subsample - 1.0) + EULER_MASCHERONI)
- 2.0 * (subsample - 1.0) / subsample
)
def score_samples(self, X: Matrix) -> Vector:
design = np.ascontiguousarray(X, dtype=np.float64)
n = design.shape[0]
total_depth = np.zeros(n, dtype=np.float64) # allocated once
cursor = np.empty(n, dtype=np.int32)
rows = np.arange(n)
for root in self._roots:
cursor.fill(root)
depth = np.zeros(n, dtype=np.float64)
active = rows
# One iteration per LEVEL. Every point still at an internal node
# advances together; trees are at most log2(256) = 8 deep, so this
# loop runs eight times for the whole batch instead of n times.
while active.size:
at = cursor[active]
internal = self._feature[at] >= 0
if not internal.any():
break
moving = active[internal]
nodes = cursor[moving]
goes_left = (
design[moving, self._feature[nodes]] <= self._threshold[nodes]
)
cursor[moving] = np.where(goes_left, self._left[nodes], self._right[nodes])
depth[moving] += 1.0
active = moving
# Truncation correction is a gather from the precomputed table.
np.add(total_depth, depth + self._correction[cursor], out=total_depth)
expected = total_depth / self._roots.size
return 2.0 ** (-expected / self._normalizer)What changed: A single isolation tree is so small — eight levels over 256 points — that scoring is dominated by interpreter overhead rather than arithmetic: the previous stage paid one Python-level call per tree per point. Packing every tree into one contiguous int32 block and traversing level-synchronously inverts that: all points still at an internal node advance together, so the interpreter runs once per tree per level rather than once per tree per point, which for a 100-tree forest is 800 iterations for an entire batch. The truncation correction moves into a precomputed per-node table, so the harmonic expansion is evaluated once at construction rather than at every leaf arrival.
Eliminate Python-level loops over samples
Level-synchronous traversal advances the whole batch one level at a time, replacing an interpreter call per point per tree with one per level per tree.
Trade-off: Every point pays for the deepest path in the tree rather than its own, and the traversal is considerably harder to read than a recursive walk — a real maintenance cost for a real speedup.
Ensure contiguous memory layout and a single dtype
All trees are packed into one int32 block per field, so the gathers performed at every level hit a small contiguous region that stays cache-resident across the whole batch.
Trade-off: Trees can no longer be added or removed individually — the forest must be repacked — and the int32 node indices cap total forest size, which is ample here only because each tree is tiny.
Pre-allocate output arrays and use in-place operations
The depth accumulator and the cursor array are allocated once for the batch and reused across all t trees, and the truncation correction is a table lookup rather than a computation.
Trade-off: Both buffers are mutable state shared across the tree loop, so the method is not reentrant and two threads scoring through one object would corrupt each other.
O(t * log psi) interpreter steps for a whole batch. Illustrative, not a measured benchmark.