This is the third post in the series on model-based outlier detection. The first post covered the basics and robust statistics, and the second introduced Isolation Forest, which is strongest at global anomalies. In this post I shall discuss the Local Outlier Factor (LOF), which judges each point against the density of its own neighbourhood, and build it from scratch before using scikit-learn.
Why "local" matters
Suppose a dataset has a tight cluster and a loose cluster. Points in the loose cluster are naturally far from their neighbours, so a method that only looks at absolute distances (or at how easily a point is separated from the rest) treats the loose cluster's fringe as suspicious. Meanwhile, a point sitting a short distance outside the tight cluster looks harmless in absolute terms, even though it is clearly out of place relative to its surroundings.
LOF (Breunig, Kriegel, Ng and Sander, 2000) fixes this by asking a relative question: is this point's neighbourhood much sparser than the neighbourhoods of its own neighbours?
The definitions
Fix an integer \(k\) (called MinPts in the original paper) and a distance \(d\).
1. k-distance. The distance from a point \(A\) to its \(k\)-th nearest neighbour:
$$k\text{-dist}(A) = d(A, o_k)$$where \(o_k\) is the \(k\)-th nearest neighbour of \(A\). The k-neighbourhood \(N_k(A)\) contains every other point within that distance, which is \(k\) points unless there are ties.
2. Reachability distance. For two points \(A\) and \(B\):
$$\text{reach-dist}_k(A, B) = \max\{\, k\text{-dist}(B),\; d(A, B) \,\}$$If \(A\) is far from \(B\), this is just the actual distance. If \(A\) lies inside \(B\)'s own neighbourhood, the distance is rounded up to \(k\text{-dist}(B)\). This smooths out the random fluctuations of very small distances and makes the density estimates more stable.
3. Local reachability density (lrd). The inverse of the average reachability distance from \(A\) to its neighbours:
$$\text{lrd}_k(A) = \left( \frac{1}{|N_k(A)|} \sum_{B \in N_k(A)} \text{reach-dist}_k(A, B) \right)^{-1}$$A point surrounded by close neighbours has a high density, and an isolated point has a low one.
4. Local Outlier Factor. The average ratio of the neighbours' density to the point's own density:
$$\text{LOF}_k(A) = \frac{1}{|N_k(A)|} \sum_{B \in N_k(A)} \frac{\text{lrd}_k(B)}{\text{lrd}_k(A)}$$The interpretation is simple.
- \(\text{LOF} \approx 1\): the point is about as dense as its neighbours, so it is a regular point.
- \(\text{LOF} > 1\): the point is sparser than its neighbours, which makes it an outlier candidate. The larger the value, the stronger the evidence.
- \(\text{LOF} < 1\): the point is denser than its neighbours, which is typical of a cluster's core.
Implementing it from scratch
The code below is a direct translation of the four definitions. It uses a full distance matrix, so it is only suitable for small datasets. scikit-learn uses tree-based neighbour search instead.
import numpy as np
from scipy.spatial.distance import cdist
def lof_scores(X, k=20):
"""Local Outlier Factor from scratch (O(n^2) memory, fine for small data)."""
n = len(X)
D = cdist(X, X)
np.fill_diagonal(D, np.inf) # a point is not its own neighbour
k_dist = np.sort(D, axis=1)[:, k - 1] # k-distance: distance to the k-th neighbour
nbrs = [np.where(D[i] <= k_dist[i])[0] for i in range(n)] # N_k(i), ties included
lrd = np.empty(n) # local reachability density
for i in range(n):
reach = np.maximum(k_dist[nbrs[i]], D[i, nbrs[i]]) # reach-dist_k(i, o)
lrd[i] = 1.0 / reach.mean()
return np.array([lrd[nbrs[i]].mean() / lrd[i] for i in range(n)])
A dataset where "local" matters
Now build a dataset to show the difference. It has a tight cluster (300 points), a loose cluster (100 points), 4 planted local outliers close to the tight cluster, and 3 planted global outliers far from everything. I constructed this dataset deliberately to illustrate the point, so on real data the gap between methods will usually be smaller. The code also compares our implementation with scikit-learn's LocalOutlierFactor.
from scipy.stats import spearmanr
from sklearn.neighbors import LocalOutlierFactor
rng = np.random.RandomState(7)
dense = rng.normal([0, 0], 0.4, size=(300, 2)) # tight cluster
sparse = rng.normal([6, 6], 1.6, size=(100, 2)) # loose cluster
local_out = np.array([[1.8, 0.2], [-0.3, 1.9], [-1.6, -1.2], [0.6, -1.9]]) # near the dense cluster
global_out = np.array([[-5.0, 8.0], [12.0, -3.0], [10.0, 12.0]]) # far from everything
X = np.vstack([dense, sparse, local_out, global_out])
n_reg = len(dense) + len(sparse)
y = np.r_[np.zeros(n_reg), np.ones(len(local_out) + len(global_out))] # 1 = planted outlier
is_local = np.r_[np.zeros(n_reg), np.ones(len(local_out)), np.zeros(len(global_out))].astype(bool)
is_global = np.r_[np.zeros(n_reg + len(local_out)), np.ones(len(global_out))].astype(bool)
k = 20
lof_mine = lof_scores(X, k)
lof_sk = LocalOutlierFactor(n_neighbors=k).fit(X)
lof_sklearn = -lof_sk.negative_outlier_factor_ # sklearn stores the NEGATIVE of LOF
print("max |ours - sklearn| = %.2e" % np.abs(lof_mine - lof_sklearn).max())
print("Spearman correlation = %.4f" % spearmanr(lof_mine, lof_sklearn).correlation)
print("LOF, regular points: median = %.2f, 99th percentile = %.2f, max = %.2f"
% (np.median(lof_mine[y == 0]), np.percentile(lof_mine[y == 0], 99), lof_mine[y == 0].max()))
print("LOF, local outliers :", np.round(lof_mine[is_local], 2))
print("LOF, global outliers:", np.round(lof_mine[is_global], 2))
flagged = LocalOutlierFactor(n_neighbors=k).fit_predict(X) == -1
print("offset_ ('auto') = %.1f; flagged = %d, planted caught = %d/%d, false alarms = %d"
% (lof_sk.offset_, flagged.sum(), (flagged & (y == 1)).sum(), y.sum(), (flagged & (y == 0)).sum()))
The output shows the following.
- Our implementation matches scikit-learn. The largest difference between the two sets of scores is about \(6\times10^{-9}\), and the rank correlation is 1.0.
- Regular points score near 1, with a median LOF of 1.02, though a few reach 2.73 (the 99th percentile is 2.04).
- Planted outliers score much higher. The four local outliers have LOF values of 3.70, 4.87, 4.26 and 3.95, and the global outliers have 21.32, 5.84 and 3.56.
- The default threshold is crude. With
contamination='auto', scikit-learn flags points with LOF above 1.5 (offset_is −1.5). Here that catches all 7 planted outliers but also flags 23 regular points. Treat the score as a ranking and choose the cutoff from the cost of false alarms, as in the Isolation Forest post.
Note also that the global outlier at (−5, 8) has a lower LOF (3.56) than the local outliers. LOF measures relative density, not distance, so it should not be read as "how far away a point is."
One scikit-learn detail: negative_outlier_factor_ stores the negative of LOF, so higher means more normal. This is why the code uses -lof_sk.negative_outlier_factor_.
LOF against a global method and Isolation Forest
Now compare LOF with two methods that look at the data more globally: the distance to the \(k\)-th nearest neighbour (the raw quantity that LOF normalises) and Isolation Forest. The code ranks all points by each method's score, checks how many of the 7 planted outliers appear in the top 7, and computes the ROC AUC (the probability that a random planted outlier scores higher than a random regular point, so 1.0 is perfect). It then varies \(k\) for LOF.
from sklearn.ensemble import IsolationForest
from sklearn.metrics import roc_auc_score
from sklearn.neighbors import NearestNeighbors
knn_dist = NearestNeighbors(n_neighbors=k + 1).fit(X).kneighbors(X)[0][:, k] # column 0 is the point itself
iso_score = -IsolationForest(random_state=0).fit(X).score_samples(X)
scores = {"k-NN distance": knn_dist, "Isolation Forest": iso_score, "LOF": lof_mine}
n_planted = int(y.sum())
for name, s in scores.items():
top = np.argsort(s)[::-1][:n_planted] # the n_planted highest-scoring points
print(f"{name:17s} ROC AUC = {roc_auc_score(y, s):.3f} | top {n_planted}: "
f"{int(is_local[top].sum())}/{is_local.sum()} local, "
f"{int(is_global[top].sum())}/{is_global.sum()} global")
print()
for kk in [5, 10, 20, 40, 80]:
s = lof_scores(X, kk)
top = np.argsort(s)[::-1][:n_planted]
print(f"LOF k={kk:2d}: ROC AUC = {roc_auc_score(y, s):.3f}, planted in top {n_planted} = {int(y[top].sum())}")
| Method | ROC AUC | Local outliers in top 7 (of 4) | Global outliers in top 7 (of 3) |
|---|---|---|---|
| k-NN distance | 0.925 | 0 | 3 |
| Isolation Forest | 0.978 | 1 | 3 |
| LOF (k = 20) | 1.000 | 4 | 3 |
All three methods find the global outliers. Only LOF finds the local ones. The figure shows why.
import matplotlib.pyplot as plt
fig, axes = plt.subplots(2, 2, figsize=(11, 9), sharex=True, sharey=True)
def draw(ax, title, top=None):
ax.scatter(X[y == 0, 0], X[y == 0, 1], s=10, c="lightgrey", edgecolors="none")
if top is None: # panel 1: show the planted outliers
ax.scatter(X[is_local, 0], X[is_local, 1], s=60, c="tab:orange", edgecolors="k", label="planted local outliers")
ax.scatter(X[is_global, 0], X[is_global, 1], s=60, c="tab:red", edgecolors="k", label="planted global outliers")
else: # other panels: show the top-scoring points
hit = np.isin(np.arange(len(X)), top) & (y == 1)
miss = np.isin(np.arange(len(X)), top) & (y == 0)
ax.scatter(X[hit, 0], X[hit, 1], s=60, c="tab:green", edgecolors="k", label="top-scored, planted outlier")
ax.scatter(X[miss, 0], X[miss, 1], s=60, c="tab:blue", marker="X", edgecolors="k", label="top-scored, regular point")
ax.set_title(title)
ax.legend(loc="upper left", fontsize=8)
draw(axes[0, 0], "The data and the planted outliers")
for ax, (name, s) in zip([axes[0, 1], axes[1, 0], axes[1, 1]], scores.items()):
top = np.argsort(s)[::-1][:n_planted]
draw(ax, f"{name}: top {n_planted} points", top)
for ax in axes[1]:
ax.set_xlabel("x1")
for ax in axes[:, 0]:
ax.set_ylabel("x2")
plt.tight_layout()
plt.show()
Both the k-NN distance and Isolation Forest spend part of their top 7 on regular points at the fringe of the loose cluster (the blue crosses), because those points really are far from their neighbours or easy to isolate in absolute terms. They miss the planted points beside the tight cluster for the opposite reason: those points are close to something, so they look unremarkable. LOF compares each point with the density around it, so it ranks all 7 planted points on top.
Choosing k
The output of the last loop in the previous block shows how LOF depends on \(k\).
| k | ROC AUC | Planted outliers in top 7 (of 7) |
|---|---|---|
| 5 | 0.999 | 5 |
| 10 | 1.000 | 7 |
| 20 | 1.000 | 7 |
| 40 | 1.000 | 7 |
| 80 | 0.995 | 5 |
A very small \(k\) makes the density estimates noisy, and a very large \(k\) blurs the neighbourhood across clusters, so the sweet spot is in the middle. The original paper recommends using at least 10 and, if you can, evaluating a range of values and taking each point's maximum LOF over that range. In practice, try several values of \(k\) and favour points that stay high across the range. Your own data's cluster sizes matter here: \(k\) should be smaller than the smallest group you want to treat as a cluster.
Novelty detection with LOF
By default, LocalOutlierFactor is for outlier detection on the training data itself, and you use fit_predict or negative_outlier_factor_. To score new observations against clean training data, set novelty=True. Then predict, decision_function and score_samples become available, and fit_predict no longer is.
X_clean = np.vstack([dense, sparse]) # training data assumed clean
lof_nov = LocalOutlierFactor(n_neighbors=20, novelty=True).fit(X_clean)
X_new = np.array([[0.1, -0.1], # typical point in the dense cluster
[6.5, 5.0], # typical point in the sparse cluster
[1.7, 0.1], # just outside the dense cluster
[9.0, -4.0]]) # far from everything
print("predict:", lof_nov.predict(X_new)) # +1 inlier, -1 outlier
print("LOF :", np.round(-lof_nov.score_samples(X_new), 2)) # score_samples = -LOF
The typical points get +1 (LOF about 1.00 and 0.97), the point just outside the tight cluster gets −1 (LOF 3.31), and the distant point gets −1 (LOF 5.87). Remember that score_samples again returns the negative of LOF.
Practical notes
- Scale your features. LOF is built on distances, so a feature with a large range dominates. Standardise first (I skipped it here only because both coordinates are on the same scale).
- Distances degrade in high dimensions, so consider PCA or feature selection before LOF on wide data.
- Duplicates can cause trouble. If more than \(k\) points are identical, the reachability distances become zero and the density blows up. Deduplicate or add tiny jitter.
- LOF is not a probability. A value of 1.5 means something different in every dataset, as the 2.73 maximum among our regular points shows.
- Cost: the neighbour search dominates. scikit-learn's tree-based search handles moderate datasets comfortably, but for very large data Isolation Forest is much cheaper.
- No single winner. Isolation Forest is fast and good at global anomalies, and LOF is better when density varies across the data. Using both and looking at where they disagree is often informative.
Conclusion
In this post I discussed how LOF compares each point's density with that of its neighbours, built it from scratch, checked it against scikit-learn, and showed on a deliberately constructed dataset how it catches local outliers that a global method misses. In the next post I shall discuss the One-Class SVM, which learns a boundary around the normal data in a kernel-induced feature space.
Reference: M. M. Breunig, H.-P. Kriegel, R. T. Ng and J. Sander, "LOF: Identifying Density-Based Local Outliers", Proceedings of the ACM SIGMOD International Conference on Management of Data, 2000.
The code in this post was tested with Python 3 and scikit-learn 1.8.0.
No comments:
Post a Comment