Nonlinear Embedding

In the two lessons the projection was always of the same form: every new number was a weighted sum of the old columns. This constraint made the error computable and let it be read from the singular values. It had a cost too — a weighted sum projects onto a flat subspace, and if the points do not stand on a flat subspace, that plane cannot follow them. Two points standing close to each other on a curved surface may land far from each other in a projection that presses the surface onto a plane.

The approaches this lesson covers remove that constraint. What they share is this: they target neighborhoods, not distances. Whom each point is close to is determined once in the original space, and then a layout is sought in a low-dimensional space that rebuilds that neighborhood relation as well as it can. There are dozens of distinct procedures in the family and they differ in the details; here the common core is written and measured. The measure has to change too, because the reconstruction error is undefined here: there is no path from the embedding back to the original space. The question to ask is this — of the k nearest neighbors in the original space, how many remain neighbors in the embedding.

• DR17. The set measured is the test set: 252 points, 12 scaled columns. The scale is again read from the training set only.
• DR18. The neighborhood preservation rate is the average, over all points, of how many of each point’s $k$ nearest neighbors in the original space are also among its $k$ nearest neighbors in the embedding. $k = 10$ is taken.
• DR19. The baseline is the random layout: 252 points are scattered randomly on the plane. Its expected value is $k / (N-1)$ and three draws are averaged. This baseline was chosen because it directly gives the neighborhood preservation of a layout that does not learn.
• DR20. The second comparison is the linear projection: the first two and the first three of the component directions of the first lesson. Two dimensions is the basis, so that a comparison can be made at the same number of dimensions.
• DR21. The layout is written within the lesson and is the model: it is started randomly, at every step neighbors are pulled toward each other and randomly chosen non-neighbor points are pushed away, and the step size falls linearly to zero. No library is called.
• DR22. For the pushing, 10 random samples are drawn per point. If all pairs were scanned, the cost would grow with the square of the number of points.
• DR23. The tuning budget is six candidates: two step counts with three step sizes. All the candidates tried are printed.
• DR24. The selection is made with the measured number itself and no separate set is used. This means the number takes part of its value from the criterion itself; the range of the six candidates is printed for this reason.
• DR25. The embedding has no unit. In the distance comparisons the embedding is scaled so that its mean pairwise distance equals that of the original space; otherwise the comparison would measure the choice of scale.
• DR26. Two runs are made with the same setting from two separate starting points. The sensitivity to the starting point is part of the measure.

The Neighborhood Preservation Rate and Its Baseline

The measure is the intersection of two sets. For a point, the ten nearest neighbors in the original space are one set and the ten nearest neighbors in the embedding are a second; the size of the intersection is divided by ten and averaged over all points. In a random layout the expected value is as much as ten neighbors being chosen at random among 251 candidates.

# COMMON — MODEL: the FICTIONAL subscriber table of the Introduction to Machine Learning course is rebuilt with the same seed.
# The label column is never generated in this course.
import math

SEED, RAW, M32 = 20260218, 1400, 0xFFFFFFFF
REGION = [("north", 0.28, 21), ("south", 0.22, 17), ("east", 0.18, 26),
          ("west", 0.14, 14), ("central", 0.18, 23)]


def rng(t):
    x = ((t ^ (t >> 16)) * 2246822507) & M32
    x = ((x ^ (x >> 13)) * 3266489909) & M32
    s = [(x ^ (x >> 16)) & M32]

    def next_value():
        s[0] = (s[0] * 1664525 + 1013904223) & M32
        return s[0] / 4294967296
    return next_value


def pick(u, w):
    t = 0.0
    for i, x in enumerate(w):
        t += x
        if u < t:
            return i
    return len(w) - 1


SUBSCRIBERS, DATA = [], []
for i in range(RAW):
    r = rng(SEED + i)
    b = REGION[pick(r(), [x[1] for x in REGION])]
    household = 1 + pick(r(), [0.06, 0.24, 0.30, 0.24, 0.11, 0.05])
    rating = 1 + pick(r(), [0.08, 0.14, 0.27, 0.34, 0.17])
    if r() >= 0.046:
        SUBSCRIBERS.append({"id": 10001 + i, "region": b, "household": household,
                            "satisfaction": rating})
for k in SUBSCRIBERS:
    r, v = rng(SEED + 7000 + k["id"]), []
    for n in range(pick(r(), [0.05, 0.12, 0.21, 0.62])):
        v.append(0.0 if r() < 0.038 else math.floor(
            k["region"][2] * math.exp((r() + r() + r() - 1.5) * 0.62)
            * (1 - n * 0.05) * 100 + 0.5) / 100)
        r()
    if not v:
        continue
    mean = sum(v) / len(v)
    x = {"household": k["household"], "satisfaction": k["satisfaction"], "period": len(v),
         "mean_consumption": round(mean, 2), "volatility": round(max(v) - min(v), 2),
         "zero_readings": int(min(v) == 0.0), "per_person": round(mean / k["household"], 3)}
    for b in REGION:
        x["area_" + b[0]] = int(k["region"][0] == b[0])
    DATA.append(x)

FIELDS = ["mean_consumption", "volatility", "per_person", "household", "satisfaction", "period",
          "zero_readings", "area_north", "area_south", "area_east", "area_west", "area_central"]
r, s = rng(SEED + 90000), list(range(len(DATA)))
for i in range(len(s) - 1, 0, -1):
    j = int(r() * (i + 1))
    s[i], s[j] = s[j], s[i]
K = [DATA[i] for i in s]
TRAIN, TEST = K[:756], K[1008:]

n = len(FIELDS)
MEAN = {a: sum(x[a] for x in TRAIN) / len(TRAIN) for a in FIELDS}
STD_DEV = {a: max(1e-9, (sum((x[a] - MEAN[a]) ** 2 for x in TRAIN) / len(TRAIN)) ** 0.5)
           for a in FIELDS}
Z_TRAIN = [[(x[a] - MEAN[a]) / STD_DEV[a] for a in FIELDS] for x in TRAIN]
Z = [[(x[a] - MEAN[a]) / STD_DEV[a] for a in FIELDS] for x in TEST]
N = len(Z)

def multiply(m, v):
    return [sum(m[i][j] * v[j] for j in range(len(v))) for i in range(len(m))]


def orthogonalize(v, B):
    for u in B:
        c = sum(v[i] * u[i] for i in range(len(v)))
        v = [v[i] - c * u[i] for i in range(len(v))]
    return v


def find_directions(Y):
    C = [[sum(z[i] * z[j] for z in Y) / len(Y) for j in range(n)] for i in range(n)]
    B = []
    for _ in range(n):
        v, lam = [1.0 / (i + 1) for i in range(n)], 0.0
        for _ in range(200):
            w = orthogonalize(multiply(C, v), B)
            u = sum(q * q for q in w) ** 0.5
            if u < 1e-9:
                break
            v, lam = [q / u for q in w], u
        if lam < 1e-9:
            break
        B.append(v)
    return B


COMPONENTS, NEIGHBOR_COUNT = find_directions(Z_TRAIN), 10


def distance(P, a, b):
    return sum((P[a][i] - P[b][i]) ** 2 for i in range(len(P[a]))) ** 0.5


def nearest(P, k):
    return [{b for _, b in sorted((distance(P, a, b), b)
                                  for b in range(len(P)) if b != a)[:k]}
            for a in range(len(P))]


NEIGHBORS = nearest(Z, NEIGHBOR_COUNT)


def preservation(P, k=NEIGHBOR_COUNT):
    q = nearest(P, k)
    return sum(len(NEIGHBORS[a] & q[a]) for a in range(N)) / (N * k)


def project(d):
    return [[sum(z[j] * COMPONENTS[i][j] for j in range(n)) for i in range(d)] for z in Z]


def random_layout(seed, d=2):
    r = rng(seed)
    return [[r() * 2 - 1 for _ in range(d)] for _ in range(N)]


rs = sum(preservation(random_layout(SEED + 7 * (t + 1))) for t in range(3)) / 3
print(f"points {N}, columns {n}, neighbor count {NEIGHBOR_COUNT}")
print(f"baseline: random layout {rs:.4f} (expected {NEIGHBOR_COUNT / (N - 1):.4f})")
print(f"linear projection 2 dimensions {preservation(project(2)):.4f}")
print(f"linear projection 3 dimensions {preservation(project(3)):.4f}")

points 252, columns 12, neighbor count 10
baseline: random layout 0.0390 (expected 0.0398)
linear projection 2 dimensions 0.2639
linear projection 3 dimensions 0.3849

The baseline is 0.0390 and it sits on the theoretical expectation, that is, 0.0398. The linear projection gives 0.2639 in two dimensions: 0.2249 above the baseline, but only about a quarter of the ten nearest neighbors. In the distance table of the first lesson, two components had brought the mean of the pairwise distances down to 0.4965 times; this row is the same loss in terms of neighborhoods. Distances do not shorten proportionally, and their ordering is disturbed too.

A Layout That Targets Neighborhoods

The procedure has three parts. For each point the ten nearest neighbors in the original space are determined once. A random layout is started on the plane. Then at each step two forces are applied: neighbors are pulled toward each other, and randomly chosen non-neighbor points are pushed away. The pushing force falls with distance, so that pairs that already stand far apart do not push each other to infinity. The step size falls linearly to zero.

def embed(seed, steps, eta, k=NEIGHBOR_COUNT, neg=10):
    # MODEL: the neighborhood-preserving layout is written at small scale, using only the standard library.
    r = rng(seed)
    Y = [[(r() - 0.5) * 0.1 for _ in range(2)] for _ in range(N)]
    nl = [sorted(NEIGHBORS[a]) for a in range(N)]
    for it in range(steps):
        h = eta * (1 - it / steps)
        for a in range(N):
            gx = gy = 0.0
            for b in nl[a]:
                gx += Y[a][0] - Y[b][0]
                gy += Y[a][1] - Y[b][1]
            for _ in range(neg):
                c = int(r() * N)
                if c == a or c in NEIGHBORS[a]:
                    continue
                dx, dy = Y[a][0] - Y[c][0], Y[a][1] - Y[c][1]
                w = 1.0 / (1.0 + dx * dx + dy * dy)
                gx -= w * dx
                gy -= w * dy
            Y[a][0] -= h * gx / k
            Y[a][1] -= h * gy / k
    return Y


print(f"{'steps':>5} {'step size':>10} {'preserved':>9}")
best = None
for steps in (200, 400):
    for eta in (0.05, 0.10, 0.20):
        d = preservation(embed(SEED + 11, steps, eta))
        print(f"{steps:>5} {eta:>10.2f} {d:>9.4f}")
        if best is None or d > best[0]:
            best = (d, steps, eta)
print(f"candidates tried 6, chosen steps {best[1]}, step size {best[2]:.2f}")
EMBEDDING = embed(SEED + 11, best[1], best[2])
SECOND = embed(SEED + 22, best[1], best[2])
print(f"embedding {preservation(EMBEDDING):.4f}, second run with the same setting {preservation(SECOND):.4f}")

steps  step size preserved
  200       0.05    0.4099
  200       0.10    0.4869
  200       0.20    0.4286
  400       0.05    0.5032
  400       0.10    0.4377
  400       0.20    0.3940
candidates tried 6, chosen steps 400, step size 0.05
embedding 0.5032, second run with the same setting 0.4806

Three numbers stand side by side. The baseline 0.0390, the linear projection 0.2639, the embedding 0.5032. The embedding’s difference from the linear projection is 0.2393, and this is about the same size as the difference from the baseline to the linear projection (0.2249). In the same two dimensions, targeting neighborhoods almost doubles the neighborhood preservation.

Where the difference comes from is read in the tuning table. The best of the six candidates is 0.5032, the worst 0.3940. Even the worst candidate exceeds the linear projection by 0.1301, so the advantage does not come from the tuning choice. But because the chosen number was chosen with the measured number itself, 0.5032 is an upper bound; the honest range is 0.3940 to 0.5032. The sensitivity to the starting point points the same way: the second run with the same setting gives 0.4806. As with clustering, an embedding produces not a single number but a range of numbers.

What Distance in the Embedding Says

A high neighborhood preservation rate does not mean the embedding preserves distances. They are separate questions, and the embedding answers only the first. There are two ways to measure this: does the embedding preserve the farthest pairs too, and how scattered are the original distances of pairs that stand at the same distance in the embedding.

def farthest(P, k):
    return [{b for _, b in sorted((-distance(P, a, b), b)
                                  for b in range(len(P)) if b != a)[:k]}
            for a in range(len(P))]


FARTHEST = farthest(Z, NEIGHBOR_COUNT)


def mean_distance(P):
    return sum(distance(P, a, b) for a in range(N)
               for b in range(a + 1, N)) / (N * (N - 1) / 2)


ORIG = mean_distance(Z)
print(f"original mean pairwise distance {ORIG:.3f}")
for name, P in (("linear 2", project(2)), ("embedding", EMBEDDING), ("second run", SECOND)):
    q = farthest(P, NEIGHBOR_COUNT)
    far = sum(len(FARTHEST[a] & q[a]) for a in range(N)) / (N * NEIGHBOR_COUNT)
    scale = ORIG / mean_distance(P)
    Q = [[y * scale for y in p] for p in P]
    band = sorted(distance(Z, a, c) for a in range(N) for c in range(a + 1, N)
                  if 0.95 * ORIG <= distance(Q, a, c) <= 1.05 * ORIG)
    print(f"{name:<13} far preservation {far:.4f}  band {len(band):>4} pairs  original distance "
          f"{band[0]:.2f}-{band[-1]:.2f}  percentile 10-90 {band[len(band) // 10]:.2f}-"
          f"{band[len(band) * 9 // 10]:.2f}")
f1, f2 = ORIG / mean_distance(EMBEDDING), ORIG / mean_distance(SECOND)
Q1 = [[y * f1 for y in p] for p in EMBEDDING]
Q2 = [[y * f2 for y in p] for p in SECOND]
s, c, large = 0.0, 0, 0
for a in range(N):
    for b in range(a + 1, N):
        d1, d2 = distance(Q1, a, b), distance(Q2, a, b)
        o = abs(d1 - d2) / max(d1, d2, 1e-9)
        s, c, large = s + o, c + 1, large + (o > 0.5)
print(f"relative difference of the same pair between the two runs: mean {s / c:.4f}, "
      f"fraction of pairs deviating by more than half {large / c:.4f}")

original mean pairwise distance 4.695
linear 2      far preservation 0.6012  band 1899 pairs  original distance 2.87-7.20  percentile 10-90 3.93-5.74
embedding     far preservation 0.0071  band 2085 pairs  original distance 3.47-9.46  percentile 10-90 4.04-6.79
second run    far preservation 0.0000  band 5439 pairs  original distance 3.39-9.21  percentile 10-90 3.97-6.01
relative difference of the same pair between the two runs: mean 0.2906, fraction of pairs deviating by more than half 0.2088

The first column is unexpectedly sharp. The linear projection preserves 0.6012 of the ten farthest pairs; the embedding 0.0071, and the second run none. The two approaches preserve the opposite of each other. The linear projection carries large distances because it has chosen the directions that carry the most variability, and that is exactly what separates far pairs; it crushes close neighbors. The embedding rescues close neighborhoods and in return throws away the large scale of distance entirely. If two points stand far apart in the embedding, that does not mean they are far apart in the original space — it only means they are not neighbors.

The second measurement says the same thing with a single range. After the mean pairwise distance is equalized, the original distance of the 2,085 pairs that stand within five percent of 4.695 in the embedding varies between 3.47 and 9.46; even the range between the 10th and 90th percentiles is 4.04–6.79. In the same band the range of the linear projection is 2.87–7.20 and its 10th–90th percentile range is 3.93–5.74, that is, markedly narrower. To look at the distance between two points in an embedding and say “they differ this much” is to read a random number from a range whose largest is 2.7 times its smallest.

The last line closes the decision. Changing only the starting point with the same setting moves the embedded distance of the same pair by a mean relative difference of 0.2906, and in 0.2088 of the pairs the difference is larger than half. The coordinates of an embedding, and the distances between the coordinates, depend on the run. The only number that can be reported for an embedding is the neighborhood preservation rate; axes, directions and distances cannot be reported.

Summary

• The baseline of the neighborhood preservation rate is the random layout and it is 0.0390; with the theoretical expectation 0.0398 it sits on it.
• The linear projection gives 0.2639 in two dimensions and the embedding that targets neighborhoods 0.5032; the difference is 0.2393, and since even the worst of the six candidates, 0.3940, exceeds the linear projection, it does not come from the tuning choice.
• The number is a range: the tuning range is 0.3940–0.5032, and the second run with the same setting is 0.4806.
• The two approaches preserve the opposite of each other: the linear projection preserves 0.6012 of the farthest pairs, the embedding 0.0071.
• The embedding makes no claim of preserving distances: after the mean distance is equalized, the original distance of pairs standing in the same band in the embedding varies between 3.47 and 9.46, and in the second run the distance of the same pair wanders by a mean relative difference of 0.2906.

Next Step

In three lessons the points were reduced to fewer numbers, and each time what was lost was measured: the reconstruction error, the storage cost, the neighborhood preservation rate. The three measures had one thing in common — all of them were averaged over all points. A projection was counted good if it was good on the average of 252 points, and an embedding was counted as preserving half if it preserved half of the neighborhoods. The single points that stay below the average — a subscriber, for example, whose reconstruction error is ten times that of the others — were never looked at. Yet such a point may itself be a finding: a record that does not fit the table and stands alone. The next lesson turns this question into a measure — how a point that stands alone is flagged, and how many of those flagged really are so.