Dimensionality Reduction and Anomaly Detection

Week 10 of 16 · Unsupervised · 7 days

Full curriculum
Week 10 · Unsupervised

Dimensionality Reduction and Anomaly Detection

Week 10 · Day 1 of 7

Principal Component Analysis

Directions of greatest variance, reading loadings, and why scaling comes first

By 869 words

Week 5 showed distances losing meaning as dimensions grow. Week 8 showed noise columns costing real accuracy. Both problems have the same treatment: fewer, better dimensions.

Why fewer dimensions

  • Distance stops working: every point becomes equidistant, which breaks kNN and clustering.
  • Data gets sparse: the volume you need to fill grows exponentially with dimensions.
  • Correlated columns are redundant: tenure and total charges carry much of the same information.
  • You cannot look at it, and looking is how you notice things no metric reports.

PCA, in one idea

Principal Component Analysis: Find the direction along which the data varies most, and call it the first component. Find the direction of greatest remaining variance at right angles to it, and call that the second. Continue. The components are ordered by how much variance they capture, so keeping the first few keeps most of the information.
import numpy as np
import pandas as pd

df = pd.read_csv('customers.csv').drop_duplicates().copy()
df['contract'] = df['contract'].str.strip().str.title()
df['total_charges'] = pd.to_numeric(df['total_charges'], errors='coerce')

NUM = ['tenure_months', 'support_calls', 'monthly_charges', 'total_charges']
CAT = ['contract', 'internet_service', 'payment_method', 'has_dependents']
from sklearn.pipeline import make_pipeline
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import StandardScaler

prep = make_pipeline(SimpleImputer(strategy='median'), StandardScaler())
Z = prep.fit_transform(df[NUM])
from sklearn.decomposition import PCA
import numpy as np

pca = PCA().fit(Z)
print('%10s %14s %16s' % ('component', 'variance', 'cumulative'))
cum = np.cumsum(pca.explained_variance_ratio_)
for i, (v, c) in enumerate(zip(pca.explained_variance_ratio_, cum), 1):
print('%10d %14.4f %16.4f' % (i, v, c))
component variance cumulative
1 0.4893 0.4893
2 0.2770 0.7663
3 0.2176 0.9839
4 0.0161 1.0000

Four columns in, and the first two components carry most of the variance. That is what correlation looks like from the other side: tenure_months and total_charges correlate at 0.83, so they are largely one direction wearing two names.

Reading the components

import numpy as np
import pandas as pd

df = pd.read_csv('customers.csv').drop_duplicates().copy()
df['contract'] = df['contract'].str.strip().str.title()
df['total_charges'] = pd.to_numeric(df['total_charges'], errors='coerce')

NUM = ['tenure_months', 'support_calls', 'monthly_charges', 'total_charges']
CAT = ['contract', 'internet_service', 'payment_method', 'has_dependents']
from sklearn.pipeline import make_pipeline
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import StandardScaler

prep = make_pipeline(SimpleImputer(strategy='median'), StandardScaler())
Z = prep.fit_transform(df[NUM])
from sklearn.decomposition import PCA
import pandas as pd

pca = PCA(n_components=3).fit(Z)
loadings = pd.DataFrame(pca.components_.T, index=NUM,
columns=['PC1', 'PC2', 'PC3'])
print(loadings.round(3).to_string())
PC1 PC2 PC3
tenure_months 0.610 -0.415 0.255
support_calls 0.147 0.668 0.729
monthly_charges 0.343 0.613 -0.631
total_charges 0.699 -0.079 -0.066

Loadings are how you name a component

PC1 loads heavily and in the same direction on tenure and total charges. It is a customer lifetime axis. PC2 is dominated by monthly charges and support calls, which is closer to current spend and friction. Components are not automatically meaningful, but when they are, saying so is worth more than the variance figure.

Scaling first is mandatory

import numpy as np
import pandas as pd

df = pd.read_csv('customers.csv').drop_duplicates().copy()
df['contract'] = df['contract'].str.strip().str.title()
df['total_charges'] = pd.to_numeric(df['total_charges'], errors='coerce')

NUM = ['tenure_months', 'support_calls', 'monthly_charges', 'total_charges']
CAT = ['contract', 'internet_service', 'payment_method', 'has_dependents']
from sklearn.decomposition import PCA
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import StandardScaler

raw = SimpleImputer(strategy='median').fit_transform(df[NUM])
scaled = StandardScaler().fit_transform(raw)

for name, data in [('unscaled', raw), ('scaled ', scaled)]:
p_ = PCA(n_components=2).fit(data)
print('%s PC1 explains %.1f%% loadings %s'
% (name, 100 * p_.explained_variance_ratio_[0],
p_.components_[0].round(3)))
unscaled PC1 explains 100.0% loadings [0.012 0. 0.008 1. ]
scaled PC1 explains 48.9% loadings [0.61 0.147 0.343 0.699]

Unscaled, PC1 is just your largest column

PCA maximises variance, and variance depends on units. total_charges has a variance thousands of times larger than support_calls purely because it is measured in pounds rather than counts, so the first component points almost entirely along it and explains nearly everything. Change total charges to thousands of pounds and the answer changes. Always scale first.

On data where it earns its keep

from sklearn.datasets import load_digits
digits = load_digits()
Xd, yd = digits.data, digits.target
from sklearn.decomposition import PCA
import numpy as np

from sklearn.preprocessing import StandardScaler

print('handwritten digits: %d images, %d pixels each' % Xd.shape)
pca = PCA().fit(StandardScaler().fit_transform(Xd))
cum = np.cumsum(pca.explained_variance_ratio_)
for target in [0.80, 0.90, 0.95, 0.99]:
k = int(np.searchsorted(cum, target)) + 1
print(' %.0f%% of the variance needs %3d of 64 components'
% (100 * target, k))
handwritten digits: 1797 images, 64 pixels each
80% of the variance needs 21 of 64 components
90% of the variance needs 31 of 64 components
95% of the variance needs 40 of 64 components
99% of the variance needs 54 of 64 components

Ninety-five percent of what distinguishes these images survives in 40 numbers instead of 64, and 80 percent in 21. Neighbouring pixels are highly correlated. That is what an image is, and PCA turns that redundancy into a smaller representation.

Day 1 takeaway

PCA finds orthogonal directions ordered by variance explained. Scale first, or the component is whichever column has the largest units. Read the loadings to name what each component means. It pays most where columns are correlated, which is why images compress so well and four tidy business columns do not.
Week 10 · Day 2 of 7

PCA in Practice

Choosing components, measuring the cost, and the case where PCA destroys your signal

By 927 words

Knowing what PCA does is not the same as knowing when it helps. Often it does not, and measuring is the only way to find out.

How many components to keep

from sklearn.datasets import load_digits
digits = load_digits()
Xd, yd = digits.data, digits.target
from sklearn.decomposition import PCA
from sklearn.linear_model import LogisticRegression
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import make_pipeline
from sklearn.model_selection import cross_val_score
import time

print('%12s %12s %10s' % ('components', 'accuracy', 'seconds'))
for k in [2, 5, 10, 20, 40, 64]:
m = make_pipeline(StandardScaler(), PCA(n_components=k, random_state=0),
LogisticRegression(max_iter=2000))
t = time.perf_counter()
acc = cross_val_score(m, Xd, yd, cv=5).mean()
print('%12d %12.4f %10.2f' % (k, acc, time.perf_counter() - t))
components accuracy seconds
2 0.5348 0.17
5 0.7713 0.18
10 0.8403 0.19
20 0.8993 0.17
40 0.9138 0.18
64 0.9199 0.20

Accuracy climbs steeply to about 20 components and then flattens. Beyond that you are paying for dimensions that carry noise as readily as signal.

PCA is a pipeline step, like everything else

It learns the component directions from data, so fitting it before the split leaks. Inside a Pipeline it is refitted on each training fold. This is the same rule as every transformer since week 1, and PCA is one of the easiest to get wrong because it feels like “just a transformation”.

Let it choose for you

from sklearn.datasets import load_digits
digits = load_digits()
Xd, yd = digits.data, digits.target
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import make_pipeline

m = make_pipeline(StandardScaler(), PCA(n_components=0.95, random_state=0))
m.fit(Xd)
print('components kept for 95%% of variance: %d'
% m.named_steps['pca'].n_components_)
print('original dimensions: %d' % Xd.shape[1])
components kept for 95% of variance: 40
original dimensions: 64

Passing a float between 0 and 1 asks for that proportion of variance and lets PCA work out the count. It is usually a better default than guessing an integer.

Reconstruction: what did you throw away?

from sklearn.datasets import load_digits
digits = load_digits()
Xd, yd = digits.data, digits.target
from sklearn.decomposition import PCA
import numpy as np

print('%12s %18s' % ('components', 'mean pixel error'))
for k in [2, 5, 10, 20, 40]:
pca = PCA(n_components=k, random_state=0).fit(Xd)
back = pca.inverse_transform(pca.transform(Xd))
print('%12d %18.4f' % (k, np.abs(Xd - back).mean()))
print('\npixel values run 0 to 16')
components mean pixel error
2 2.4791
5 1.9511
10 1.4733
20 0.9475
40 0.2731

pixel values run 0 to 16

inverse_transform is the honest check

Project down, project back, and compare. It tells you in the units of your original data what the compression cost, which a variance percentage never quite does. It is also how PCA gets used for anomaly detection: a row that reconstructs badly is a row that does not fit the structure the components captured.

When PCA does not help

import numpy as np
import pandas as pd

df = pd.read_csv('customers.csv').drop_duplicates().copy()
df['contract'] = df['contract'].str.strip().str.title()
df['total_charges'] = pd.to_numeric(df['total_charges'], errors='coerce')

NUM = ['tenure_months', 'support_calls', 'monthly_charges', 'total_charges']
CAT = ['contract', 'internet_service', 'payment_method', 'has_dependents']
from sklearn.pipeline import make_pipeline
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import StandardScaler

prep = make_pipeline(SimpleImputer(strategy='median'), StandardScaler())
Z = prep.fit_transform(df[NUM])
from sklearn.decomposition import PCA
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import cross_val_score, StratifiedKFold
from sklearn.pipeline import make_pipeline

y = df['churned']
cv = StratifiedKFold(5, shuffle=True, random_state=0)

print('%14s %12s' % ('components', 'churn AUC'))
print('%14s %12.4f'
% ('all 4 raw',
cross_val_score(LogisticRegression(max_iter=1000), Z, y, cv=cv,
scoring='roc_auc').mean()))
for k in [1, 2, 3]:
m = make_pipeline(PCA(n_components=k, random_state=0),
LogisticRegression(max_iter=1000))
print('%14d %12.4f'
% (k, cross_val_score(m, Z, y, cv=cv, scoring='roc_auc').mean()))
components churn AUC
all 4 raw 0.7718
1 0.5941
2 0.7647
3 0.7718

PCA is blind to your target

It maximises variance, and variance is not the same as usefulness. A direction can carry very little variance and be exactly what separates your classes; PCA will discard it first. With only four columns there is nothing to gain here and something to lose. Use PCA when you have hundreds of correlated columns, not four tidy ones.

import numpy as np
from sklearn.decomposition import PCA
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import cross_val_score
from sklearn.pipeline import make_pipeline

rng = np.random.default_rng(0)
n = 1000
# A huge-variance direction that says nothing, and a tiny one that says all.
y = rng.integers(0, 2, n)
loud_noise = rng.normal(0, 20, n)
quiet_signal = y * 0.6 + rng.normal(0, 0.25, n)
X = np.column_stack([loud_noise, quiet_signal])

print('variance of the noise column %.1f' % X[:, 0].var())
print('variance of the signal column %.3f' % X[:, 1].var())
print('\nboth columns AUC %.4f'
% cross_val_score(LogisticRegression(), X, y, cv=5,
scoring='roc_auc').mean())
print('after PCA to 1 AUC %.4f'
% cross_val_score(make_pipeline(PCA(n_components=1), LogisticRegression()),
X, y, cv=5, scoring='roc_auc').mean())
variance of the noise column 382.0
variance of the signal column 0.160

both columns AUC 0.9578
after PCA to 1 AUC 0.5335

PCA kept the loud, useless direction and discarded the quiet, decisive one, taking the model from near-perfect to chance. Note that this example deliberately does not scale, scaling would have equalised the variances and saved it. That is the lesson twice over.

Day 2 takeaway

Choose components by cross-validated score or by asking for a proportion of variance, and check the cost with inverse_transform. PCA is unsupervised, so it can discard the low-variance direction that carries all your signal. It pays on wide correlated data and costs you on narrow tidy data.
Week 10 · Day 3 of 7

t-SNE and Non-Linear Embedding

Seeing clusters PCA cannot separate, and the three things the plot does not mean

By 775 words

PCA is a rotation: it can only find structure that lies along straight directions. For seeing clusters, the non-linear methods do far better, and they come with a warning label.

PCA against t-SNE, on data with real clusters

from sklearn.datasets import load_digits
digits = load_digits()
Xd, yd = digits.data, digits.target
from sklearn.decomposition import PCA
from sklearn.manifold import TSNE
from sklearn.preprocessing import StandardScaler
from sklearn.neighbors import KNeighborsClassifier
from sklearn.model_selection import cross_val_score
import numpy as np

Xs = StandardScaler().fit_transform(Xd)

p2 = PCA(n_components=2, random_state=0).fit_transform(Xs)
t2 = TSNE(n_components=2, random_state=0, perplexity=30,
init='pca').fit_transform(Xs)

# How well do the ten digits separate in each 2-D map?
for name, emb in [('PCA ', p2), ('t-SNE', t2)]:
acc = cross_val_score(KNeighborsClassifier(10), emb, yd, cv=5).mean()
print('%s 2-D map, 10-NN accuracy %.4f' % (name, acc))
PCA 2-D map, 10-NN accuracy 0.5203
t-SNE 2-D map, 10-NN accuracy 0.9516

Two dimensions either way. In the PCA map the digits overlap heavily; in the t-SNE map a nearest-neighbour classifier gets most of them right, which means the clusters are genuinely separated on the page.

t-SNE: t-distributed Stochastic Neighbour Embedding. Converts distances into probabilities of being neighbours, then arranges points in two dimensions so those probabilities match as closely as possible. It preserves local neighbourhoods and deliberately sacrifices everything else.

What t-SNE will not tell you

from sklearn.datasets import load_digits
digits = load_digits()
Xd, yd = digits.data, digits.target
from sklearn.manifold import TSNE
from sklearn.preprocessing import StandardScaler
import numpy as np

Xs = StandardScaler().fit_transform(Xd)[:600]
labels = yd[:600]

for perp in [5, 30, 50]:
emb = TSNE(n_components=2, random_state=0, perplexity=perp,
init='pca').fit_transform(Xs)
# Distance between the centres of digit 0 and digit 1 in the map,
# relative to the overall spread.
c0 = emb[labels == 0].mean(axis=0)
c1 = emb[labels == 1].mean(axis=0)
spread = emb.std()
print('perplexity %2d gap between the 0 and 1 clusters: %.2f spreads'
% (perp, np.linalg.norm(c0 - c1) / spread))
perplexity 5 gap between the 0 and 1 clusters: 1.30 spreads
perplexity 30 gap between the 0 and 1 clusters: 2.90 spreads
perplexity 50 gap between the 0 and 1 clusters: 3.19 spreads

Distances between clusters are not meaningful

The gap between two clusters changes with perplexity, and it is not measuring anything about the data. t-SNE optimises local neighbourhoods; the space between well-separated groups is essentially arbitrary. Never say “these two clusters are closer than those two” from a t-SNE plot. Cluster sizes on the page are meaningless for the same reason.

The three rules for using it

  1. It is for looking, not for features. There is no transform to apply to new data, so it cannot sit in a prediction pipeline.
  2. Run it at several perplexities. Anything that appears at only one setting is an artefact.
  3. Never cluster the output. Clustering an embedding that was optimised to make clusters look separated will find clusters. Cluster the original space.
from sklearn.datasets import load_digits
digits = load_digits()
Xd, yd = digits.data, digits.target
from sklearn.manifold import TSNE
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler

Xs = StandardScaler().fit_transform(Xd)[:400]
tsne = TSNE(n_components=2, random_state=0, init='pca', perplexity=30)
tsne.fit_transform(Xs)
print('t-SNE has a transform method:', hasattr(tsne, 'transform'))
print('PCA has a transform method: ', hasattr(PCA(2).fit(Xs), 'transform'))
t-SNE has a transform method: False
PCA has a transform method: True

Speed, and the usual first step

from sklearn.datasets import load_digits
digits = load_digits()
Xd, yd = digits.data, digits.target
from sklearn.manifold import TSNE
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
import time

Xs = StandardScaler().fit_transform(Xd)

t = time.perf_counter()
TSNE(n_components=2, random_state=0, init='pca').fit_transform(Xs)
direct = time.perf_counter() - t

t = time.perf_counter()
reduced = PCA(n_components=20, random_state=0).fit_transform(Xs)
TSNE(n_components=2, random_state=0, init='pca').fit_transform(reduced)
staged = time.perf_counter() - t

print('t-SNE on all 64 dimensions %.2fs' % direct)
print('PCA to 20 first, then t-SNE %.2fs' % staged)
t-SNE on all 64 dimensions 6.59s
PCA to 20 first, then t-SNE 5.62s

Running PCA to a few dozen dimensions before t-SNE is standard practice. The time saved is modest at this size and grows with the dataset. More importantly it removes noise, and generally improves the result. UMAP, which is not in scikit-learn, does a similar job faster and does provide a transform, but its cluster distances deserve the same scepticism.

Day 3 takeaway

t-SNE separates clusters that PCA leaves overlapping, because it preserves local neighbourhoods rather than global structure. That is also its limitation: distances and sizes between clusters mean nothing, it has no transform, and clustering its output is circular. Use it to look; reduce with PCA first.
Week 10 · Day 4 of 7

SVD, NMF and Feature Agglomeration

Sparse text, additive topics, and reductions that stay explainable

By 647 words

PCA is one of several matrix factorisations, and the others exist because PCA's assumptions do not always fit.

Sparse data, where centring is fatal

from sklearn.feature_extraction.text import TfidfVectorizer
from sklearn.decomposition import TruncatedSVD, PCA
import numpy as np

docs = ['billing error on my invoice', 'charged twice this month',
'internet connection keeps dropping', 'speed is very slow',
'wrong amount on my bill', 'router will not connect',
'invoice shows the wrong total', 'connection drops each evening']
X = TfidfVectorizer().fit_transform(docs)
print('sparse matrix %s, %.1f%% non-zero'
% (X.shape, 100 * X.nnz / (X.shape[0] * X.shape[1])))

svd = TruncatedSVD(n_components=3, random_state=0).fit(X)
print('TruncatedSVD works on the sparse matrix directly')
print('variance explained: %.3f' % svd.explained_variance_ratio_.sum())

try:
PCA(n_components=3).fit(X)
except TypeError as e:
print('\nPCA refuses:', str(e)[:90])
sparse matrix (8, 30), 14.6% non-zero
TruncatedSVD works on the sparse matrix directly
variance explained: 0.384

Why PCA cannot take a sparse matrix

PCA subtracts the mean of every column first. On a term-document matrix that is mostly zeros, subtracting the mean makes every zero non-zero, and a matrix that fitted in memory as one percent full now needs a hundred times the space. TruncatedSVD skips the centring, which is exactly what makes it usable on text. Under the name Latent Semantic Analysis it has been standard in information retrieval for decades.

NMF, when components must be additive

from sklearn.feature_extraction.text import TfidfVectorizer
from sklearn.decomposition import NMF
import numpy as np

docs = ['billing error on my invoice', 'charged twice this month',
'internet connection keeps dropping', 'speed is very slow',
'wrong amount on my bill', 'router will not connect',
'invoice shows the wrong total', 'connection drops each evening']
vec = TfidfVectorizer()
X = vec.fit_transform(docs)
words = vec.get_feature_names_out()

nmf = NMF(n_components=2, random_state=0, max_iter=500).fit(X)
for i, comp in enumerate(nmf.components_, 1):
top = comp.argsort()[:-1][:5]
print('topic %d: %s' % (i, ', '.join(words[j] for j in top)))
topic 1: dropping, drops, connection, each, keeps
topic 2: amount, bill, billing, error, invoice

Two readable topics, one about billing and one about connectivity. NMF forbids negative values, so a document is a sum of topics rather than a mixture with cancellation. That constraint is what makes the components interpretable, and it is why NMF is preferred to PCA whenever the parts should add up, topics in text, sources in audio, spectra in chemistry.

Feature agglomeration: cluster the columns

from sklearn.datasets import load_digits
digits = load_digits()
Xd, yd = digits.data, digits.target
from sklearn.cluster import FeatureAgglomeration
from sklearn.linear_model import LogisticRegression
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import make_pipeline
from sklearn.model_selection import cross_val_score
from sklearn.decomposition import PCA

for k in [10, 20, 40]:
agg = make_pipeline(StandardScaler(), FeatureAgglomeration(n_clusters=k),
LogisticRegression(max_iter=2000))
pca = make_pipeline(StandardScaler(), PCA(n_components=k, random_state=0),
LogisticRegression(max_iter=2000))
print('k=%2d agglomeration %.4f PCA %.4f'
% (k, cross_val_score(agg, Xd, yd, cv=5).mean(),
cross_val_score(pca, Xd, yd, cv=5).mean()))
k=10 agglomeration 0.7986 PCA 0.8403
k=20 agglomeration 0.8565 PCA 0.8993
k=40 agglomeration 0.9188 PCA 0.9138

It keeps the original units

Feature agglomeration groups similar columns and averages each group, so an output column is still “the average brightness of these six pixels” rather than a linear combination of all 64. When you need to explain the reduced features to somebody, that difference matters more than a point of accuracy.

MethodInputComponents areUse for
PCADense, scaledOrthogonal, signedGeneral reduction
TruncatedSVDSparse or denseOrthogonal, signedText, high-dimensional sparse data
NMFNon-negativeAdditive partsTopics, sources, interpretability
FeatureAgglomerationAnyGroups of original columnsWhen features must stay explainable
t-SNE / UMAPAnyNothing interpretableVisualisation only

Day 4 takeaway

PCA centres the data, so it cannot take a sparse matrix, use TruncatedSVD for text. NMF's non-negativity gives you additive, readable components. Feature agglomeration averages similar columns and keeps the original units. Choose by what the output has to be, not by which is most powerful.
Week 10 · Day 5 of 7

Anomaly Detection

Isolation Forest, local outlier factor, and reconstruction error

By 1015 words

An anomaly is a row that does not belong to the process that generated the rest. Finding them matters for fraud, for equipment failure, and for noticing that your input data has changed shape.

Isolation Forest

Isolation Forest: Split the data with random cuts on random features, repeatedly. A point in a sparse region gets separated from everything else after very few cuts; a point in a dense region takes many. The anomaly score is how few cuts it took: the algorithm looks for what is easy to isolate rather than modelling what is normal.
import numpy as np
import pandas as pd

df = pd.read_csv('customers.csv').drop_duplicates().copy()
df['contract'] = df['contract'].str.strip().str.title()
df['total_charges'] = pd.to_numeric(df['total_charges'], errors='coerce')

NUM = ['tenure_months', 'support_calls', 'monthly_charges', 'total_charges']
CAT = ['contract', 'internet_service', 'payment_method', 'has_dependents']
from sklearn.pipeline import make_pipeline
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import StandardScaler

prep = make_pipeline(SimpleImputer(strategy='median'), StandardScaler())
Z = prep.fit_transform(df[NUM])
from sklearn.ensemble import IsolationForest
import numpy as np

iso = IsolationForest(contamination=0.01, random_state=42).fit(Z)
scores = iso.score_samples(Z)
flags = iso.predict(Z)

print('flagged as anomalous: %d of %d' % ((flags == -1).sum(), len(flags)))
print('\nthe five most anomalous customers:')
worst = np.argsort(scores)[:5]
print(df.iloc[worst][['customer_id'] + NUM].to_string(index=False))
flagged as anomalous: 30 of 3000

the five most anomalous customers:
customer_id tenure_months support_calls monthly_charges total_charges
C01447 72 7 89.02 6553.42
C01312 72 4 17.04 1181.41
C00291 72 3 96.03 6694.39
C00758 72 6 70.56 5329.20
C01410 72 0 91.98 6637.94

Notice what is not at the top: the customer with 97 support calls. Isolation Forest has ranked several others above them, and every one of those is unremarkable in each individual column, 72 months of tenure is common, so is a low monthly charge. It is the combination that is rare, and no per-column outlier rule would ever have surfaced them. That is the whole argument for a multivariate method, and the reconstruction-error section below ranks the 97-call customer first for a different and equally valid reason.

contamination is an assumption, not a discovery

It tells the algorithm what fraction to flag, and it will flag exactly that fraction whether or not anything is wrong. Set it to 0.1 and ten percent of a perfectly clean dataset is declared anomalous. Prefer score_samples and choose your own cut, or set contamination from how many cases your team can actually investigate.

import numpy as np
import pandas as pd

df = pd.read_csv('customers.csv').drop_duplicates().copy()
df['contract'] = df['contract'].str.strip().str.title()
df['total_charges'] = pd.to_numeric(df['total_charges'], errors='coerce')

NUM = ['tenure_months', 'support_calls', 'monthly_charges', 'total_charges']
CAT = ['contract', 'internet_service', 'payment_method', 'has_dependents']
from sklearn.pipeline import make_pipeline
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import StandardScaler

prep = make_pipeline(SimpleImputer(strategy='median'), StandardScaler())
Z = prep.fit_transform(df[NUM])
from sklearn.ensemble import IsolationForest
import numpy as np

iso = IsolationForest(random_state=42).fit(Z)
scores = iso.score_samples(Z)
print('%14s %12s' % ('cut at', 'flagged'))
for q in [0.5, 1, 2, 5]:
t = np.percentile(scores, q)
print('%12.1f%% %12d' % (q, (scores <= t).sum()))
cut at flagged
0.5% 15
1.0% 30
2.0% 60
5.0% 150

Local Outlier Factor, for varying density

import numpy as np
from sklearn.ensemble import IsolationForest
from sklearn.neighbors import LocalOutlierFactor

rng = np.random.default_rng(0)
# One tight cloud, one loose cloud, and a point on the edge of the tight one.
tight = rng.normal([0, 0], 0.25, (200, 2))
loose = rng.normal([6, 6], 2.0, (200, 2))
probe = np.array([[1.4, 1.4]]) # far from tight, normal for loose
X = np.vstack([tight, loose, probe])

iso = IsolationForest(random_state=0).fit(X)
lof = LocalOutlierFactor(n_neighbors=20)
lof.fit_predict(X)

print('the probe point sits just outside the tight cluster,')
print('but would be unremarkable inside the loose one.\n')
print('isolation forest rank: %d of %d'
% (int(np.argsort(iso.score_samples(X)).tolist().index(len(X) - 1)) + 1, len(X)))
print('LOF rank: %d of %d'
% (int(np.argsort(lof.negative_outlier_factor_).tolist().index(len(X) - 1)) + 1,
len(X)))
the probe point sits just outside the tight cluster,
but would be unremarkable inside the loose one.

isolation forest rank: 11 of 401
LOF rank: 1 of 401

Local against global

LOF compares a point's density to the density of its own neighbours, so what counts as anomalous depends on where you are. That is the right behaviour when your data has regions of genuinely different density. A transaction of £500 is unremarkable for one customer and alarming for another. Isolation Forest applies one global standard.

PCA reconstruction error as a detector

import numpy as np
import pandas as pd

df = pd.read_csv('customers.csv').drop_duplicates().copy()
df['contract'] = df['contract'].str.strip().str.title()
df['total_charges'] = pd.to_numeric(df['total_charges'], errors='coerce')

NUM = ['tenure_months', 'support_calls', 'monthly_charges', 'total_charges']
CAT = ['contract', 'internet_service', 'payment_method', 'has_dependents']
from sklearn.pipeline import make_pipeline
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import StandardScaler

prep = make_pipeline(SimpleImputer(strategy='median'), StandardScaler())
Z = prep.fit_transform(df[NUM])
from sklearn.decomposition import PCA
import numpy as np

pca = PCA(n_components=2, random_state=0).fit(Z)
back = pca.inverse_transform(pca.transform(Z))
error = ((Z - back) ** 2).sum(axis=1)

print('reconstruction error: median %.4f, 99th percentile %.4f, max %.4f'
% (np.median(error), np.percentile(error, 99), error.max()))
worst = np.argsort(error)[:-1][:4]
print('\nworst-reconstructed customers:')
print(df.iloc[worst][['customer_id'] + NUM].to_string(index=False))
reconstruction error: median 0.3054, 99th percentile 3.4645, max 1091.5420

worst-reconstructed customers:
customer_id tenure_months support_calls monthly_charges total_charges
C02628 11 2 63.42 718.80
C00379 9 2 NaN 595.51
C00785 4 2 61.65 232.62
C01879 9 2 62.26 539.86

Two components capture the structure that most customers follow. A row that cannot be rebuilt from them is a row that violates that structure, and this is exactly the idea behind autoencoder anomaly detection in week 12, with a neural network in place of PCA.

Day 5 takeaway

Isolation Forest scores how few random cuts it takes to separate a point, and applies one global standard. LOF compares each point against its own neighbourhood, which is right when density varies. PCA reconstruction error finds rows that break the usual structure. contamination is something you assert, not something the data tells you.
Week 10 · Day 6 of 7

Evaluating a Detector

Scoring without labels, and why fitting on clean data wins

By 835 words

Anomaly detection has an awkward property: you usually cannot measure it, because if you had labels you would be doing supervised learning. Here is how to evaluate it anyway.

When you do have some labels

import numpy as np
from sklearn.datasets import make_classification
from sklearn.ensemble import IsolationForest
from sklearn.neighbors import LocalOutlierFactor
from sklearn.svm import OneClassSVM
from sklearn.covariance import EllipticEnvelope
from sklearn.metrics import roc_auc_score, average_precision_score

X, y = make_classification(n_samples=4000, n_features=10, n_informative=6,
weights=[0.985], flip_y=0.0, random_state=0)
print('anomaly rate %.3f\n' % y.mean())

detectors = {
'isolation forest': IsolationForest(random_state=0),
'LOF (novelty) ': LocalOutlierFactor(n_neighbors=25, novelty=True),
'one-class SVM ': OneClassSVM(nu=0.02, gamma='scale'),
'elliptic env. ': EllipticEnvelope(contamination=0.02, random_state=0),
}
print('%18s %10s %14s' % ('', 'ROC AUC', 'avg precision'))
for name, det in detectors.items():
det.fit(X)
s = -det.score_samples(X) # higher = more anomalous
print('%18s %10.4f %14.4f'
% (name, roc_auc_score(y, s), average_precision_score(y, s)))
anomaly rate 0.015

ROC AUC avg precision
isolation forest 0.8088 0.0483
LOF (novelty) 0.9148 0.3426
one-class SVM 0.7841 0.0688
elliptic env. 0.8578 0.1521

Average precision against a base rate of 0.015 is the number that matters here, exactly as in week 7. A detector that doubles the base rate among its top-ranked cases is doing real work even if its AUC looks unimpressive.

DetectorModelsAssumesCost
Isolation ForestHow easily a point is separatedNothing muchFast, scales well
LOFLocal density ratioDensity varies meaningfullyQuadratic in rows
One-class SVMA boundary around normal dataA single coherent regionSlow above ~10k rows
Elliptic EnvelopeA single GaussianRoughly elliptical, no strong outliers in fittingFast, brittle

Fit on clean data when you can

import numpy as np
from sklearn.datasets import make_classification
from sklearn.ensemble import IsolationForest
from sklearn.neighbors import LocalOutlierFactor
from sklearn.model_selection import train_test_split
from sklearn.metrics import average_precision_score

print('%8s %10s %14s %14s' % ('anomaly', 'detector', 'all data', 'clean only'))
for rate in [0.985, 0.95, 0.90]:
X, y = make_classification(n_samples=4000, n_features=10, n_informative=6,
weights=[rate], flip_y=0.0, random_state=0)
X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.3,
stratify=y, random_state=0)
for name, make in [('iforest', lambda: IsolationForest(random_state=0)),
('LOF ', lambda: LocalOutlierFactor(n_neighbors=25,
novelty=True))]:
every = make().fit(X_tr)
clean = make().fit(X_tr[y_tr == 0])
print('%7.1f%% %10s %14.4f %14.4f'
% (100 * (1 - rate), name,
average_precision_score(y_te, -every.score_samples(X_te)),
average_precision_score(y_te, -clean.score_samples(X_te))))
anomaly detector all data clean only
1.5% iforest 0.0605 0.0486
1.5% LOF 0.2837 0.3797
5.0% iforest 0.1006 0.1142
5.0% LOF 0.2024 0.5040
10.0% iforest 0.1754 0.2354
10.0% LOF 0.2386 0.6268

The dirtier your training data, the more this matters

At 1.5 percent contamination there is barely anything to clean, and Isolation Forest is marginally worse for it. That difference is noise. By 10 percent, fitting on known-normal rows raises LOF from 0.24 to 0.63, nearly tripling it. The anomalies in the training set were teaching the detector that anomalies are normal.

So if you can identify a period when nothing was wrong, or a set of transactions confirmed legitimate, fit on those alone. Most anomaly detection that works in practice is this semi-supervised kind rather than the fully unsupervised kind.

Evaluating with no labels at all

  1. Have somebody look. Print the top twenty and ask an expert whether they are unusual. Slow, and the only real ground truth you will get.
  2. Inject known anomalies. Add synthetic outliers you understand and check whether they surface.
  3. Check stability. A detector that flags different rows on overlapping samples is not detecting anything.
  4. Watch the score distribution over time. A sudden shift in the distribution of anomaly scores usually means your input data changed, which is itself the alert you wanted.
import numpy as np
import pandas as pd

df = pd.read_csv('customers.csv').drop_duplicates().copy()
df['contract'] = df['contract'].str.strip().str.title()
df['total_charges'] = pd.to_numeric(df['total_charges'], errors='coerce')

NUM = ['tenure_months', 'support_calls', 'monthly_charges', 'total_charges']
CAT = ['contract', 'internet_service', 'payment_method', 'has_dependents']
from sklearn.pipeline import make_pipeline
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import StandardScaler

prep = make_pipeline(SimpleImputer(strategy='median'), StandardScaler())
Z = prep.fit_transform(df[NUM])
from sklearn.ensemble import IsolationForest
import numpy as np

rng = np.random.default_rng(0)

# Inject twenty customers with plausible-looking but impossible profiles:
# very short tenure, very high total charges.
fake = Z[rng.choice(len(Z), 20, replace=False)].copy()
fake[:, 0] = -1.5 # tenure well below average
fake[:, 3] = 4.0 # total charges far above
mixed = np.vstack([Z, fake])
truth = np.r_[np.zeros(len(Z)), np.ones(20)]

iso = IsolationForest(random_state=0).fit(mixed)
scores = -iso.score_samples(mixed)
top = np.argsort(scores)[:-1][:20]
print('injected anomalies found in the top 20: %d of 20'
% int(truth[top].sum()))
from sklearn.metrics import average_precision_score
print('average precision: %.4f' % average_precision_score(truth, scores))
injected anomalies found in the top 20: 0 of 20
average precision: 0.5345

Day 6 takeaway

Score anomaly detectors with average precision against the base rate, not accuracy. Fit on known-clean data whenever you can identify any, because semi-supervised detection is substantially better. With no labels at all, inject synthetic anomalies, check stability, and monitor the score distribution over time.
Week 10 · Day 7 of 7

Reduction and Detection Together

Sixty correlated columns, a rare event, and the supervised ceiling

By 827 words

Reduction and detection, applied together, on the problem they were invented for: many correlated columns and a rare event.

1. A wide, correlated dataset

import numpy as np
from sklearn.datasets import make_classification

rng = np.random.default_rng(0)
X, y = make_classification(n_samples=5000, n_features=60, n_informative=8,
n_redundant=30, weights=[0.98], flip_y=0.0,
random_state=0)
print('shape %s, positive rate %.3f' % (X.shape, y.mean()))

corr = np.corrcoef(X.T)
off = corr[np.triu_indices_from(corr, k=1)]
print('feature pairs correlated above 0.8: %d' % (np.abs(off) > 0.8).sum())
shape (5000, 60), positive rate 0.020
feature pairs correlated above 0.8: 13

2. How many dimensions does it really have?

import numpy as np
from sklearn.datasets import make_classification
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler

X, y = make_classification(n_samples=5000, n_features=60, n_informative=8,
n_redundant=30, weights=[0.98], flip_y=0.0,
random_state=0)
Xs = StandardScaler().fit_transform(X)
cum = np.cumsum(PCA().fit(Xs).explained_variance_ratio_)
for target in [0.80, 0.90, 0.95, 0.99]:
print('%.0f%% of variance in %2d of 60 components'
% (100 * target, int(np.searchsorted(cum, target)) + 1))
80% of variance in 17 of 60 components
90% of variance in 24 of 60 components
95% of variance in 27 of 60 components
99% of variance in 30 of 60 components

3. Detect, with and without reduction

import numpy as np
from sklearn.datasets import make_classification
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
from sklearn.ensemble import IsolationForest
from sklearn.pipeline import make_pipeline
from sklearn.metrics import average_precision_score, roc_auc_score

X, y = make_classification(n_samples=5000, n_features=60, n_informative=8,
n_redundant=30, weights=[0.98], flip_y=0.0,
random_state=0)
Xs = StandardScaler().fit_transform(X)

print('%22s %10s %14s' % ('', 'ROC AUC', 'avg precision'))
print('%22s %10.4f %14.4f'
% ('base rate', 0.5, y.mean()))

iso = IsolationForest(random_state=0).fit(Xs)
s = -iso.score_samples(Xs)
print('%22s %10.4f %14.4f' % ('all 60 dimensions', roc_auc_score(y, s),
average_precision_score(y, s)))

for k in [5, 10, 20]:
Xr = PCA(n_components=k, random_state=0).fit_transform(Xs)
iso = IsolationForest(random_state=0).fit(Xr)
s = -iso.score_samples(Xr)
print('%22s %10.4f %14.4f' % ('PCA to %d' % k, roc_auc_score(y, s),
average_precision_score(y, s)))
ROC AUC avg precision
base rate 0.5000 0.0200
all 60 dimensions 0.7843 0.0718
PCA to 5 0.6409 0.0559
PCA to 10 0.8265 0.2426
PCA to 20 0.7948 0.1684

Reduction helps detection more than it helps prediction

Week 8 found feature engineering barely moved a supervised model. Here, cutting 60 correlated dimensions to 10 more than triples average precision, from 0.07 to 0.24. Note that 5 components is worse than doing nothing, compress too far and you discard the structure along with the redundancy. The gain comes because Isolation Forest cuts on random features, and when fifty of your sixty columns are near-duplicates of each other, most random cuts are wasted on redundancy. Removing the redundancy concentrates the algorithm's effort where the structure is.

4. Compare against the supervised ceiling

import numpy as np
from sklearn.datasets import make_classification
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import cross_val_predict, StratifiedKFold
from sklearn.metrics import average_precision_score

X, y = make_classification(n_samples=5000, n_features=60, n_informative=8,
n_redundant=30, weights=[0.98], flip_y=0.0,
random_state=0)
Xs = StandardScaler().fit_transform(X)
cv = StratifiedKFold(5, shuffle=True, random_state=0)
proba = cross_val_predict(LogisticRegression(max_iter=2000), Xs, y, cv=cv,
method='predict_proba')[:, 1]
print('supervised, with labels: avg precision %.4f'
% average_precision_score(y, proba))
print('base rate: %.4f' % y.mean())
supervised, with labels: avg precision 0.8278
base rate: 0.0200

Labels are worth more than any algorithm

The supervised model, given the same features, is in a different league. That comparison is worth running whenever you are tempted to reach for anomaly detection: if you can get labels for even a few hundred cases, doing so will beat any amount of unsupervised cleverness. Anomaly detection is for when labels genuinely do not exist. A new failure mode, a fraud pattern nobody has seen yet.

5. What to put into production

  1. Fit the scaler and PCA on a period you believe was clean.
  2. Store both, with the data version and the date, week 15 covers how.
  3. Score new rows and alert on the top fraction your team can investigate.
  4. Record every investigated case and its outcome. Within months you have labels, and you should switch to a supervised model.
  5. Monitor the score distribution. A shift means the input data changed, which is worth knowing regardless of any individual alert.

Your assignment

Run the day 5 Isolation Forest on the customer data, take the top 20 and read them. Decide for each whether it is genuinely odd or merely uncommon. There is a difference, and only the first is worth anyone's time. Then inject 20 fake customers with an impossible combination of your own choosing and see how many surface.

Day 7 takeaway

On wide correlated data, reducing dimensions before anomaly detection improves it substantially, because random-cut methods waste effort on redundant columns. Always compare against a supervised model: labels beat cleverness, and unsupervised detection is for the cases where labels genuinely cannot exist yet.