III-28 · Hierarchical and Density-Based Clustering

Decide whether a cluster is a merge or a neighborhood

Inspect actual linkage merges, then identify the core observations that can propagate a density cluster and the border points that cannot.

What makes a cluster?

Hand-authored, unitless toy observations. Labels identify groups; they have no numerical order.

Density roles: filled circle = core, hollow diamond = border, cross = noise. Selected ID and role appear above the plane; all coordinates are in the table. Drag snaps to 0.1; arrow keys preserve full input precision. Only handles capture touch. The selector brings an overlapping handle forward; Enter or Space does the same for a focused handle.

Edit handlesSelected handle

Handle x-5 to 5 · Enter to applyHandle y-5 to 5 · Enter to applyScale x0.1 to 10 · Enter to applyScale y0.1 to 10 · Enter to apply

Linkage

Merges displayed0 to 5 · Enter to applyEpsilon0 to 2 · Enter to applyMinimum points1 to 6 · Enter to apply

Frame 1 / 1 · edit begins at the displayed state.

X · samples in rowsEach mathematical feature vector is a 2 × 1 column.: 6 × 2
Centroids / squared distances / assignments: 2 × 2 / 6 × 2 / 6 × 1
Held-assignment WCSSNearest-label preview; press Assign to commit.: ≈ 4.1
Next nearest-assignment WCSSAllows labels to change; same current centroids.: ≈ 4.1
Mean update · current → proposed WCSSHeld assignments; empty centroids retained. Inspect before applying Update means.: ≈ 4.1 → ≈ 2.79333
Phase and stopNot stopped; replay is capped at 40 half-steps.: ready · 0 mean updates
Empty-cluster policyEmpty mean is undefined. Retained coordinates are a policy, not a computed empty mean.: Retain previous centroid
Input provenance / initializationSeed affects partial Fisher–Yates observation selection only. Current centroids may be edited.: Hand-authored observations / seed 7

d_{ij}^2=\sum_{r=1}^2((x_{ir}-c_{jr})/s_r)^2,\quad a_i=\operatorname*{arg\,min}_j d_{ij}^2

\mu_j=\frac{1}{|C_j|}\sum_{i\in C_j}x_i,\quad J=\sum_i d_{i,a_i}^2

Actual merge tree · height is linkage distance

[s1] · [s2] · [s5] · [s6] · [s3, s4]

single: exact ties choose the first pair by ascending persistent group IDs.
Merge	Left IDs	Right IDs	Height
1	s3	s4	≈ 0.1
2	s1	s2	≈ 0.1
3	s3, s4	s1, s2	≈ 0.1
4	s5	s1, s2, s3, s4	≈ 0.3
5	s6	s1, s2, s3, s4, s5	≈ 1.4

Closed ε-neighborhood includes self. Core iff count ≥ minimum points. Noise uses −1; a shared border belongs to the first reached core component in input order.
Observation	x₁	Neighbor IDs	Count	Role	Cluster
s1	0	s1, s2, s3, s4	4	core	0
s2	0.1	s1, s2, s3, s4	4	core	0
s3	0.2	s1, s2, s3, s4	4	core	0
s4	0.3	s1, s2, s3, s4, s5	5	core	0
s5	0.6	s4, s5	2	border	0
s6	2	s6	1	noise	-1

The feature plane keeps actual raw positions. Filled circles are core, hollow diamonds border, crosses noise. A selected-point label avoids dense annotation overlap. Inspect every coordinate and neighborhood in the table; the second SVG projects the actual linkage merge tree.

Change the definition of a group

At reset, the six x coordinates are 0, 0.1, 0.2, 0.3, 0.6 and 2; all y coordinates are zero. With ε=0.31 and minimum points 4, the first four observations are core, the fifth is border and the last is noise. Select a handle or use its numeric coordinates to inspect its closed neighborhood. Every neighborhood count is recomputed from these actual positions and the same explicit feature scales.

Now change Linkage from single to complete. The observations stay fixed, but the merge heights change because “distance between groups” changes. Replay actual merges selects successive cuts of that calculated tree. The immutable replay contains complete numerical input states; it does not manufacture a drawn tree or change DBSCAN parameters behind the scenes.

These methods answer different questions from k-means. Hierarchical clustering asks how to build nested groups. Density clustering asks which regions have enough nearby observations to propagate connectivity. Neither needs labels supplied as ground truth.

Start with one observation per group

Agglomerative clustering begins with n singleton groups. At each step it merges the closest active pair, creating a new group. The full sequence is a hierarchy; selecting a merge count creates a cut with n−count groups. A dendrogram’s vertical merge height encodes the declared linkage distance, not time, a centroid objective, or a probability.

Let $d(x,y)$ be the square root of the scaled squared Euclidean distance. For disjoint nonempty groups A and B, our three linkages are

\begin{aligned} D_{\mathrm{single}}(A,B)&=\min_{x\in A,y\in B}d(x,y),\\ D_{\mathrm{complete}}(A,B)&=\max_{x\in A,y\in B}d(x,y),\\ D_{\mathrm{average}}(A,B)&=\frac{1}{|A||B|}\sum_{x\in A,y\in B}d(x,y). \end{aligned}

Average linkage averages all cross-observation distances. It is neither distance between centroids nor the unweighted average of earlier group scores. Single linkage can connect long chains through short local links. Complete linkage waits for all cross distances to be small. These choices encode different group preferences; no one is universally correct.

Initial group IDs are input indices 0…n−1. New IDs are n, n+1, and so on. Active groups are sorted by these persistent IDs. Exact computed-distance ties choose the first enumerated pair in that order. Input row order therefore resolves ties; sample IDs remain visible identifiers. A tie does not supply evidence that one merge is naturally preferred.

Work a three-observation hierarchy

Use the tiny canonical fixture with points at x=0, 1 and 3, y=0 and unit scales. The initial distances are 1, 3 and 2. All three methods first merge rows 0 and 1 at height 1.

The remaining active groups have IDs 2 and 3, where ID 3 contains [0,1]. Their cross distances are 3 and 2:

Linkage	Calculation for [0,1] vs [2]	Final merge height
Single	min(3,2)	2
Complete	max(3,2)	3
Average	(3+2)/2	2.5

The final merge table lists left [2], right [0,1] because the persistent IDs are ordered. That left/right drawing orientation is a reproducibility convention, not a meaningful ordering of clusters. A cut after the first merge has two groups; after the second, every observation is in one group.

The Stage’s dendrogram paths use these actual merge heights and child memberships. Wide tree geometry retains readable labels in a bounded scroll region. Dense observations keep their true raw coordinates; the density plane uses one selected annotation and a complete numeric table instead of stacking unreadable labels.

Define density before propagating it

For ε≥0, the closed neighborhood includes the observation itself:

N_\varepsilon(x_i)=\{j:d(x_i,x_j)\le\varepsilon\}, \qquad \text{core}(i)\iff |N_\varepsilon(x_i)|\ge m.

Here m is Minimum points, including self. In the implementation the equivalent comparison is squared distance ≤ ε², with checked finite arithmetic. A row at exactly ε is included. Coincident rows with distinct IDs count separately, including when ε=0.

A core point can expand the search to all of its neighbors. A border point is noncore but reached from a core component; it receives that cluster label and does not expand the component. Noise is not reached from any core component. The displayed label −1 is an explicit noise marker, not a learned negative class or a prediction of an observed target.

Operationally, scan input indices in ascending order. Visit a new point; if it is core, start a cluster and expand its neighborhood with a queue. Neighbors enter that queue in ascending input order; each enters at most once per expansion. Previously visited noise can become border when a core expansion reaches it. A point already assigned to a cluster is never stolen by a later expansion.

Count a core, a border and noise literally

At the six-row reset, s1=(0,0) sees s1, s2, s3 and s4: four observations within ε=0.31. Rows s2/s3 also see that same four. Row s4=(0.3,0) additionally sees s5=(0.6,0), since their distance is 0.3, so it sees five. All four are core for m=4.

s5 sees itself and s4 only. It is noncore, yet s4’s expansion reaches it, so it is border with label 0. s6=(2,0) sees only itself. No core expansion reaches it, so it is noise with label −1. The labels are [0,0,0,0,0,−1]. Lowering m to 1 makes even an isolated point core, illustrating how parameter choice changes the definition rather than discovering a uniquely true answer.

A separate exact boundary fixture uses dyadic x coordinates 0, 1/8, 1/4, 3/8, 3/4 and 2 with ε=3/8. The first point includes the fourth at equality. The neighborhoods are [0,1,2,3], [0,1,2,3], [0,1,2,3], [0,1,2,3,4], [3,4], [5]. This rational example checks both the closed boundary and self-count rule independently of decimal rounding.

A shared border must not merge two core components

Load Shared border. Four left-core observations and four right-core observations are internally connected. The central point s9=(0,0) sees the two inner cores s4=(−0.9,0), s8=(0.9,0), plus itself. With ε=1 and m=4 it has only three neighbors, so it is border rather than core.

It touches both core components but cannot expand through itself to join them. The left component is encountered first in input order, so s9 receives its label 0. The right component retains label 1. Swapping component order could assign this ambiguous border to the right without changing core connectivity. Our convention is explicit; other DBSCAN implementations may resolve such border ownership differently.

The independent oracle builds the graph of core-to-core neighborhood edges, computes its transitive closure, and finds components. It then assigns each noncore point to the first adjacent component by smallest core input index. That is mathematically different from translating the canonical queue traversal and exposes the mistaken claim that border vertices propagate connectivity.

Limits and numerical assumptions

The batch is n×2, with each feature vector a 2×1 mathematical column. Pairwise distances and the Boolean neighborhood relation are n×n; labels are n×1. The stored neighborhoods are ragged per-row lists of neighbor indices, not an n×n matrix. At reset their lengths are [4,4,4,5,2,1]. Core flags are Boolean decisions, not confidence probabilities. Single, complete and average linkage report distances; the density method does not optimize k-means WCSS. The comparison WCSS inspector refers to the separately defined centroid model and is not a DBSCAN objective.

The educational model supports 1–32 rows, ε≥0 and 1≤m≤n; the Stage bounds ε to 0…2 and coordinates to −5…5. No observations is an invalid state, with no invented objective or hierarchy. The implementation is a transparent tiny-data mechanism, not a spatially indexed large-data implementation. It does not handle missing values or varying local density automatically. A dense bridge of actual core points may legitimately connect two regions, while a border bridge cannot.

Feature units still matter: the ε neighborhood is circular in normalized coordinates and elliptical on raw axes when reference scales differ. The visible selected neighborhood ellipse uses those same scales. Coincidence does not erase observation identity. Distinct cluster labels can be permuted; noise treatment must be declared separately in any external evaluation.

Read the actual code and connect the definitions

The imported reference implements hierarchical and dbscan from scratch, using NumPy distance arrays and explicit loops. Send an operation of hierarchical or dbscan with the fixture’s JSON input on stdin to inspect that operation. Running without stdin prints the Part’s original Lloyd mean-update example:

.venv/bin/python -B code/python/p05/clustering.py

code/python/p05/clustering.py

"""Canonical from-scratch NumPy clustering. X is n×2; centers are k×2.
No target column, fitting library, hidden reseeding, or optimum guarantee.
"""
import json
import math
import sys
from copy import deepcopy
import numpy as np


def number(v, name, lo=None, hi=None, integer=False):
    try:
        valid = isinstance(v, (int, float)) and not isinstance(v, bool) and math.isfinite(float(v))
    except (OverflowError, TypeError):
        valid = False
    if not valid:
        raise ValueError(f"{name} must be a finite number.")
    v = float(v)
    if (lo is not None and v < lo) or (hi is not None and v > hi) or (integer and v != math.floor(v)):
        raise ValueError(f"{name} is outside its supported range.")
    return int(v) if integer else v


def vector(v, name):
    if not isinstance(v, (list, tuple)) or len(v) != 2:
        raise ValueError(f"{name} needs a dense two-feature row.")
    return [number(x, name) for x in v]


def observations(points):
    if not isinstance(points, list) or not 1 <= len(points) <= 32:
        raise ValueError("Use 1–32 observations; empty data have no objective.")
    ids, rows = set(), []
    for p in points:
        if not isinstance(p, dict) or not isinstance(p.get("id"), str) or not p["id"].strip() or p["id"] in ids:
            raise ValueError("Observation IDs must be unique nonempty strings.")
        ids.add(p["id"])
        rows.append(vector(p.get("x"), "Observation"))
    return np.asarray(rows, dtype=np.float64)


def scales_vector(scales):
    values = vector(scales, "Scales")
    if any(x <= 0 for x in values):
        raise ValueError("Reference scales must be positive.")
    return values


def squared_distance(a, b, scales=(1, 1)):
    a, b, scales = vector(a, "First vector"), vector(b, "Second vector"), scales_vector(scales)
    total = 0.0
    for j in range(2):
        difference = number(a[j]-b[j], "Coordinate difference")
        delta = number(difference/scales[j], "Scaled difference")
        square = number(delta*delta, "Squared difference")
        if difference != 0 and (delta == 0 or square == 0):
            raise ValueError("Distance underflow: rescale coordinates or reference scales so distinct points remain distinguishable.")
        total = number(total+square, "Squared distance")
    return total


def distance_matrix(X, centers, scales):
    # Explicit scalar boundary checks keep NumPy warning settings from changing semantics.
    return np.asarray([[squared_distance(x.tolist(), c.tolist(), scales) for c in centers] for x in X], dtype=np.float64)


def inspect_lloyd(state):
    X = observations(state["points"])
    scales = scales_vector(state["scales"])
    centers = state["centers"]
    if not isinstance(centers, list):
        raise ValueError("Centroids must be a dense array.")
    number(len(centers), "k", 1, len(X), True)
    C = np.asarray([vector(c, "Centroid") for c in centers], dtype=np.float64)
    phase, labels = state["phase"], state["labels"]
    if phase not in ("ready", "assign", "update") or type(state["stopped"]) is not bool:
        raise ValueError("Invalid phase or stop flag.")
    number(state["iteration"], "Iteration", 0, 1000, True)
    if phase == "ready":
        if labels is not None or state["stopped"] or state["iteration"] != 0:
            raise ValueError("Invalid ready state.")
    elif not isinstance(labels, list) or len(labels) != len(X):
        raise ValueError("Assignments must match rows.")
    if state["stopped"] and phase != "update":
        raise ValueError("Stop is reported only after mean update.")
    if labels is not None:
        labels = [number(v, "Assignment", 0, len(C)-1, True) for v in labels]
    D = distance_matrix(X, C, scales)
    nearest = np.argmin(D, axis=1).tolist()  # first index on computed equality
    held = nearest if labels is None else labels
    clusters = []
    for j, c in enumerate(C):
        indices = [i for i, label in enumerate(held) if label == j]
        sums = [0.0, 0.0]
        for i in indices:
            for d in range(2):
                sums[d] = number(sums[d]+float(X[i, d]), "Cluster sum")
        mean = [number(v/len(indices), "Mean") for v in sums] if indices else c.tolist()
        clusters.append({"index":j, "ids":[state["points"][i]["id"] for i in indices], "count":len(indices), "sums":sums, "mean":mean, "empty":not indices})
    def objective(assignment, matrix=D):
        total = 0.0
        for i, label in enumerate(assignment):
            total = number(total+float(matrix[i, label]), "WCSS")
        return total
    return {"distances":D.tolist(), "nearest":nearest, "labels":held, "clusters":clusters,
            "objective":objective(held), "nearestObjective":objective(nearest),
            "proposedObjective":objective(held, distance_matrix(X, np.asarray([c["mean"] for c in clusters]), scales)),
            "ties":[sum(v == row[nearest[i]] for v in row) > 1 for i, row in enumerate(D.tolist())]}


def lloyd_step(state):
    active = inspect_lloyd(state)
    result = deepcopy(state)
    if state["stopped"]:
        return result
    if state["phase"] != "assign":
        result.update(labels=active["nearest"], phase="assign")
    else:
        centers = [c["mean"] for c in active["clusters"]]
        result.update(centers=centers, phase="update", iteration=state["iteration"]+1, stopped=centers == state["centers"])
    inspect_lloyd(result)
    return result


def lloyd_trace(state, half_steps=40):
    number(half_steps, "Half steps", 1, 100, True)
    frames = [deepcopy(state)]
    inspect_lloyd(state)
    for _ in range(half_steps):
        if frames[-1]["stopped"]:
            break
        frames.append(lloyd_step(frames[-1]))
    return frames


def initialize(points, k, seed):
    observations(points)
    k = number(k, "k", 1, len(points), True)
    state = number(seed, "Seed", 0, 2**32-1, True)
    order = list(range(len(points)))
    for i in range(k):
        state = (1664525*state+1013904223) % 2**32
        j = i+math.floor((state/2**32)*(len(points)-i))
        order[i], order[j] = order[j], order[i]
    return [vector(points[i]["x"], "Selected centroid") for i in order[:k]]


def same_partition(a, b):
    if not isinstance(a, list) or not isinstance(b, list) or len(a) != len(b):
        raise ValueError("Partitions need equal row counts.")
    a = [number(x, "Label", 0, 31, True) for x in a]
    b = [number(x, "Label", 0, 31, True) for x in b]
    return all((a[i] == a[j]) == (b[i] == b[j]) for i in range(len(a)) for j in range(i))


def hierarchical(points, linkage, scales=(1, 1)):
    X = observations(points)
    scales_vector(scales)
    if linkage not in ("single", "complete", "average"):
        raise ValueError("Choose single, complete or average linkage.")
    D = np.sqrt(distance_matrix(X, X, scales)).tolist()
    groups = [{"id":i, "members":[i]} for i in range(len(X))]
    merges = []
    while len(groups) > 1:
        best = None
        for i in range(len(groups)):
            for j in range(i+1, len(groups)):
                values = [D[a][b] for a in groups[i]["members"] for b in groups[j]["members"]]
                if linkage == "single":
                    height = min(values)
                elif linkage == "complete":
                    height = max(values)
                else:
                    total = 0.0
                    for d in values:
                        total = number(total+d, "Linkage sum")
                    height = number(total/len(values), "Average linkage")
                if best is None or height < best[0]:
                    best = height, i, j
        height, i, j = best
        left, right = groups[i], groups[j]
        group = {"id":len(X)+len(merges), "members":sorted(left["members"]+right["members"])}
        merges.append({"left":deepcopy(left), "right":deepcopy(right), "group":deepcopy(group), "height":height})
        groups = sorted([g for g in groups if g["id"] not in (left["id"], right["id"])]+[group], key=lambda g:g["id"])
    return {"distances":D, "merges":merges}


def dbscan(points, epsilon, min_points, scales=(1, 1)):
    X = observations(points)
    epsilon = number(epsilon, "Epsilon", 0)
    min_points = number(min_points, "Minimum points", 1, len(X), True)
    limit = number(epsilon*epsilon, "Squared epsilon")
    D = distance_matrix(X, X, scales).tolist()
    neighborhoods = [[i for i, d in enumerate(row) if d <= limit] for row in D]
    core = [len(ns) >= min_points for ns in neighborhoods]
    labels, visited, count = [-1]*len(X), [False]*len(X), 0
    for i in range(len(X)):
        if visited[i]:
            continue
        visited[i] = True
        if not core[i]:
            continue
        labels[i] = count
        queue, queued = list(neighborhoods[i]), set(neighborhoods[i])
        q = 0
        while q < len(queue):
            j = queue[q]
            if not visited[j]:
                visited[j] = True
                if core[j]:
                    for v in neighborhoods[j]:
                        if v not in queued:
                            queue.append(v); queued.add(v)
            if labels[j] == -1:
                labels[j] = count
            q += 1
        count += 1
    return {"distances":D, "neighborhoods":neighborhoods, "core":core, "labels":labels, "count":count,
            "kinds":["core" if core[i] else "noise" if label < 0 else "border" for i, label in enumerate(labels)]}


def dispatch(request):
    operation = request["operation"]
    if operation == "lloyd":
        state = request["state"]
        before = inspect_lloyd(state)
        after_state = lloyd_step(state)
        after = inspect_lloyd(after_state)
        return {"distances":before["distances"], "labels":before["labels"], "means":[c["mean"] for c in before["clusters"]],
                "empty":[c["empty"] for c in before["clusters"]], "objective":before["objective"],
                "proposedObjective":before["proposedObjective"], "afterCenters":after_state["centers"], "afterObjective":after["objective"]}
    if operation == "hierarchical":
        result = hierarchical(request["points"], request["linkage"], request.get("scales", [1, 1]))
        return {"merges":[{"left":m["left"]["members"], "right":m["right"]["members"], "height":m["height"]} for m in result["merges"]]}
    if operation == "dbscan":
        result = dbscan(request["points"], request["epsilon"], request["minPoints"], request.get("scales", [1, 1]))
        return {k:result[k] for k in ("neighborhoods", "core", "labels", "kinds", "count")}
    if operation == "initialize":
        return {"centers":initialize(request["points"], request["k"], request["seed"])}
    raise ValueError("Unknown clustering operation.")


if __name__ == "__main__":
    raw = "" if sys.stdin.isatty() else sys.stdin.read()
    default = {"operation":"lloyd", "state":{"points":[{"id":f"s{i+1}", "x":[x,0]} for i, x in enumerate([-3,-1,1,3])], "centers":[[-3,0],[1,0]], "scales":[1,1], "labels":[0,0,1,1], "phase":"assign", "iteration":0, "stopped":False}}
    try:
        print(json.dumps(dispatch(json.loads(raw) if raw.strip() else default), allow_nan=False, indent=2))
    except (ValueError, KeyError, TypeError, OverflowError) as error:
        print(f"Clustering input error: {error}", file=sys.stderr)
        raise SystemExit(2)

III-26 optimized squared distance to means; III-27 showed its initialization and scale dependence. Hierarchical linkage now supplies nested relations; DBSCAN supplies core connectivity and an explicit noise category. Select the definition that matches the use of a grouping before comparing algorithm names.

Five graded investigations

Checkpoints

You can calculate literal linkage heights, read a genuine merge tree and its cut, enumerate closed neighborhoods, and distinguish core propagation from border ownership and noise. The choice of cluster definition, input traversal order, feature scale and noise convention is part of the result. It is not evidence that an algorithm’s numeric labels are true categories.