When I see a bird that walks like a duck and swims like a duck and quacks like a duck, I call that bird a duck.
–James Whitcomb Riley (1849 - 1916)
Like any classifier, \(k\)-Nearest Neighbors is trained by providing it a set of labeled data.
However, at training time, the classifier does very little. It just stores away the training data.
plt.figure(figsize=(5, 5))
demo_y = np.array([1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0])
demo_X = np.array([[-3,1], [-2, 4], [-2, 2], [-1.5, 1], [-1, 3], [0, 0], [1, 1.5], [2, 0.5], [2, 3], [2, 0], [3, 1], [4, 4], [0, 1]])
test_X = np.array([-0.3, 0.7])
#
plt.scatter(demo_X[:,0], demo_X[:,1], c=demo_y, cmap=cmap2_bold, s=80)
plt.axis('equal')
plt.axis('off')
plt.show()
Label 0 is drawn in orange, label 1 in sky blue.
plt.figure(figsize=(5, 5))
plt.scatter(demo_X[:, 0], demo_X[:, 1], c=demo_y, cmap=cmap2_bold, s=80)
plt.plot(test_X[0], test_X[1], 'ok', markersize=10)
plt.annotate('Test Point', test_X, [75, 25],
textcoords = 'offset points', fontsize = 14,
arrowprops = {'arrowstyle': '->'})
plt.axis('equal')
plt.axis('off')
plt.show()
At test time, simply
plt.figure(figsize=(5, 5))
plt.scatter(demo_X[:, 0], demo_X[:, 1], c=demo_y, cmap=cmap2_bold, s=80)
plt.plot(test_X[0], test_X[1], 'ok', markersize=10)
ax=plt.gcf().gca()
circle = mp.patches.Circle(test_X, 0.5, facecolor = OKABE_ITO[0], alpha = 0.3)
plt.axis('equal')
plt.axis('off')
ax.add_artist(circle)
plt.show()
The nearest point is orange, so 1-NN predicts orange (label 0).
plt.figure(figsize=(5, 5))
plt.scatter(demo_X[:,0], demo_X[:,1], c=demo_y, cmap=cmap2_bold, s=80)
plt.plot(test_X[0], test_X[1], 'ok', markersize=10)
ax=plt.gcf().gca()
circle = mp.patches.Circle(test_X, 0.9, facecolor = 'gray', alpha = 0.3)
plt.axis('equal')
plt.axis('off')
ax.add_artist(circle)
plt.show()
One orange and one sky-blue neighbour: a tie. The majority vote is undefined.
plt.figure(figsize=(4, 4))
ax=plt.gcf().gca()
circle = mp.patches.Circle(test_X, 1.4, facecolor = OKABE_ITO[1], alpha = 0.3)
ax.add_artist(circle)
plt.scatter(demo_X[:,0], demo_X[:,1], c=demo_y, cmap=cmap2_bold, s=80)
plt.plot(test_X[0], test_X[1], 'ok', markersize=10)
plt.axis('equal')
plt.axis('off')
plt.tight_layout()
plt.show()
Two sky-blue neighbours against one orange: 3-NN predicts sky blue (label 1).
Consider this dataset where items fall into three classes:
import sklearn.datasets as sk_data
X, y = sk_data.make_blobs(n_samples=150,
centers=[[-2, 0],[1, 5], [2.5, 1.5]],
cluster_std = [2, 2, 3],
n_features=2,
center_box=(-10.0, 10.0),random_state=0)
plt.figure(figsize = (5, 5))
plt.axis('equal')
plt.axis('off')
plt.scatter(X[:, 0], X[:, 1], c = y, cmap = cmap3_bold, s = 80)
plt.show()
Let’s observe how the complexity of the resulting model changes as we vary \(k\).
We’ll do this by plotting the decision regions. These show how the method would classify each potential test point in the space.
# Plot the decision boundary. For that, we will assign a color to each
# point in the mesh [x_min, x_max]x[y_min, y_max].
h = .1 # step size in the mesh
x_min, x_max = X[:, 0].min() - 1, X[:, 0].max() + 1
y_min, y_max = X[:, 1].min() - 1, X[:, 1].max() + 1
xx, yy = np.meshgrid(np.arange(x_min, x_max, h),
np.arange(y_min, y_max, h))f, axs = plt.subplots(1, 3, figsize=(15, 4))
for i, k in enumerate([1, 5, 25]):
knn = KNeighborsClassifier(n_neighbors = k)
knn.fit(X, y)
Z = knn.predict(np.c_[xx.ravel(), yy.ravel()])
Z = Z.reshape(xx.shape)
axs[i].pcolormesh(xx, yy, Z, cmap = cmap3_light, shading = 'auto')
# overlay the training points that produced these regions
axs[i].scatter(X[:, 0], X[:, 1], c = y, cmap = cmap3_bold, s = 25, edgecolors = 'k', linewidths = 0.5)
axs[i].axis('equal')
axs[i].axis('off')
axs[i].set_title(f'Decision Regions for $k$ = {k}')
plt.show()
With \(k = 1\) every training point owns a little island, so the boundary wraps around each stray point. Increasing \(k\) smooths the boundary: the stray points are out-voted by their neighbours.
Smoother boundaries are more likely to generalize well – but not always, as we will see.
The curse of dimensionality is a somewhat tongue-in-cheek term for serious problems that arise when we use geometric algorithms in high dimensions.
There are various aspects of the curse that affect \(k\)-NN.
plt.figure(figsize=(7, 4))
for name, colour in zip(ratio, OKABE_ITO):
plt.plot(dims, ratio[name], 'o-', color=colour, label=name)
plt.axhline(1, linestyle='--', color='gray')
plt.xscale('log')
plt.yscale('log')
plt.xlabel('Dimension $d$')
plt.ylabel('max distance / min distance')
plt.title('1000 random points: ratio of farthest to closest pair')
plt.legend()
plt.show()
A ratio of 1 (dashed) would mean every pair is at the same distance. All three metrics head toward it. Manhattan and cosine distance stay a little further from 1 than Euclidean, which is one reason they are sometimes preferred in high dimensions – but none escapes the trend.
Here is a geometric way to see the same concentration. Compare a ball of radius \(1\) with the slightly smaller ball of radius \(1-\epsilon\) inside it:
ax = plt.figure(figsize = (4, 4)).add_subplot(projection = '3d')
# coordinates of sphere surface
u, v = np.mgrid[0:2*np.pi:50j, 0:np.pi:50j]
x = np.cos(u)*np.sin(v)
y = np.sin(u)*np.sin(v)
z = np.cos(v)
#
ax.plot_surface(x, y, z, color=OKABE_ITO[0], alpha = 0.2)
s3 = 1/np.sqrt(3)
ax.quiver(0, 0, 0, s3, s3, s3, color = OKABE_ITO[5])
ax.text(s3/2, s3/2, s3/2-0.2, '1', size = 14)
#
eps = 0.9
#
ax.plot_surface(eps * x, eps * y, eps * z, color=OKABE_ITO[5], alpha = 0.2)
ax.quiver(0, 0, 0, eps, 0, 0, color = 'k')
ax.text(1/2-0.2, 0, -0.4, r'$1-\epsilon$', size = 14)
ax.set_axis_off()
plt.title('Inner and Outer Hyperspheres')
plt.show()
In 3-D the shell between them is thin. In high dimensions it is almost the whole ball.
d_vals = np.arange(1, 301)
plt.figure(figsize = (6, 4))
for eps, colour in zip([0.01, 0.1], [OKABE_ITO[5], OKABE_ITO[3]]):
plt.plot(d_vals, 1 - (1 - eps) ** d_vals, color=colour, label=fr'$\epsilon$ = {eps}')
plt.xlabel('Dimensions $d$')
plt.ylabel(r'Fraction of volume in outer $\epsilon$ shell')
plt.title(r'$f_d = 1 - (1-\epsilon)^d$')
plt.legend(loc = 'lower right')
plt.show()
With \(\epsilon = 0.1\) the shell holds 96% of the ball by \(d = 30\); with \(\epsilon = 0.01\) it holds 95% by \(d = 300\).
rng = np.random.default_rng(0) # seeded, so the figure is reproducible
nsamples = 1000
unif_X = rng.uniform(0, 1, nsamples).reshape(-1, 1)
dists = pdist(unif_X)
mean_dists = [np.mean(dists)]
min_dists = [np.min(dists)]
for d in range(2, 101):
unif_X = np.column_stack([unif_X, rng.uniform(0, 1, nsamples)])
dists = pdist(unif_X)
mean_dists.append(np.mean(dists))
min_dists.append(np.min(dists))plt.figure(figsize = (6, 3))
plt.plot(range(1, 101), min_dists, label = "Minimum Distance")
plt.plot(range(1, 101), mean_dists, label = "Average Distance")
plt.xlabel(r'Number of dimensions ($d$)')
plt.ylabel('Distance')
plt.legend(loc = 'best')
plt.title(f'Comparison of Minimum Versus Average Distance Between {nsamples} Points\nAs Dimension Grows')
plt.show()
The average distance between points grows. However we also observe that the minimum distance between points grows at a similar rate.
Let’s look at the ratio between the minimum distance between points and the average distance between points.
This shows that, for any test point \(x\), the distance to the closest point to \(x\), gets closer and closer to the average distance between points.
If we used a point at the average distance for classifying \(x\) we’d get a very poor classifier.
A better way to measure the classifier’s performance is using a Confusion Matrix.
Diagonal elements represent successes and off diagonals represent errors.
Using the confusion matrix we can define some more useful measures:
First we’ll generate some synthetic data to work with.
X, y = datasets.make_circles(noise=.1, factor=.5, random_state=1)
print('Shape of data: {}'.format(X.shape))
print('Unique labels: {}'.format(np.unique(y)))
def plot_points(X, y, title='', ax=None):
"""Scatter two-class points: label 0 orange, label 1 sky blue."""
ax = ax or plt.gca()
ax.scatter(X[:, 0], X[:, 1], c=y, cmap=cmap2_bold, vmin=0, vmax=1, s=80)
ax.set_title(title)
ax.axis('equal')
ax.axis('off')
return axShape of data: (100, 2)
Unique labels: [0 1]Recall that we always want to test on data separate from our training data.
For now, we will do something very simple: take the first 50 examples for training and the rest for testing.
Indeed, the Positive (orange) points in the upper half of the test data are all classified incorrectly – they come out sky blue.
Let’s look at one of the points that the classifier got wrong.
def plot_neighborhood(knn, k, X_train, y_train, X_test, y_pred_test, test_point):
"""Left: training data with the k neighbours of one test point circled.
Right: test predictions with that test point circled."""
neighbors = knn.kneighbors([X_test[test_point]])[1][0]
radius = np.max(np.linalg.norm(X_train[neighbors] - X_test[test_point], axis=1))
fig, axs = plt.subplots(1, 2, figsize = (10, 4))
plot_points(X_train, y_train, 'Training', ax=axs[0])
axs[0].scatter(X_train[neighbors, 0], X_train[neighbors, 1],
facecolors='none', edgecolors='k', s=160)
axs[0].add_artist(mp.patches.Circle(X_test[test_point], radius, facecolor='gray', alpha=0.25))
plot_points(X_test, y_pred_test,
f'Testing $k$={k}\nAccuracy: {knn.score(X_test, y_test)}', ax=axs[1])
axs[1].scatter(X_test[test_point, 0], X_test[test_point, 1],
facecolors='none', edgecolors='k', s=160)
plt.show()
test_point = np.argmax(X_test[:, 1]) # the top-most test point
plot_neighborhood(knn5, 5, X_train, y_train, X_test, y_pred_test, test_point)
Three of its five nearest training points are sky blue, so the vote goes against it.
For comparison purposes, let’s try \(k\) = 3.
k = 3
knn3 = KNeighborsClassifier(n_neighbors=k)
knn3.fit(X_train,y_train)
y_pred_test = knn3.predict(X_test)
fig, axs = plt.subplots(1, 2, figsize = (10, 4))
plot_points(X_train, y_train, 'Training', ax=axs[0])
plot_points(X_test, y_pred_test, f'Testing $k$={k}\nAccuracy: {knn3.score(X_test, y_test)}', ax=axs[1])
plt.show()
And let’s look at the same individual point as before.
import sklearn.model_selection as model_selection
nreps = 50
kvals = range(1, 16, 2)
acc = []
se = []
np.random.seed(4)
for k in kvals:
val_rep = []
train_rep = []
for i in range(nreps):
X_train, X_val, y_train, y_val = model_selection.train_test_split(X,
y,
test_size = 0.5)
knn = KNeighborsClassifier(n_neighbors = k)
knn.fit(X_train, y_train)
train_rep.append(knn.score(X_train, y_train))
val_rep.append(knn.score(X_val, y_val))
acc.append([np.mean(val_rep), np.mean(train_rep)])
se.append(np.std(val_rep) / np.sqrt(nreps)) # standard error of the mean
accy = np.array(acc)plt.errorbar(kvals, accy[:, 0], se, fmt = '.-', capsize = 3, label = 'Accuracy on Validation Data')
plt.plot(kvals, accy[:, 1], '.-', label = 'Accuracy on Training Data')
plt.xlabel(r'$k$')
plt.ylabel('Accuracy')
plt.xticks(kvals)
plt.title(r'Train/Validation Comparison of $k$-NN ($\pm 1$ standard error)')
plt.legend(loc = 'best')
plt.show()
print(f'Best validation accuracy at k = {kvals[np.argmax(accy[:, 0])]}: {np.max(accy[:, 0]):.3f}')
Best validation accuracy at k = 1: 0.978Here is the decision boundary for \(k\)-NN with \(k = 1\).
np.random.seed(1)
X_train, X_val, y_train, y_val = model_selection.train_test_split(X, y, test_size = 0.5)
k = 1
knn = KNeighborsClassifier(n_neighbors = k)
knn.fit(X_train, y_train)
Z = knn.predict(np.c_[xx.ravel(), yy.ravel()])
Z = Z.reshape(xx.shape)
fig, axs = plt.subplots(1, 2, figsize = (12, 5))
for ax, (XX, yy_, name) in zip(axs, [(X_train, y_train, 'Training'), (X_val, y_val, 'Validation')]):
ax.contourf(xx, yy, Z, cmap=cmap2_light, vmin=0, vmax=1)
ax.scatter(XX[:, 0], XX[:, 1], c=yy_, cmap=cmap2_bold, vmin=0, vmax=1, s=30)
ax.axis('equal')
ax.axis('off')
ax.set_xlim((x_min, x_max))
ax.set_ylim((y_min, y_max))
ax.set_title(f'{k}-NN - {name} Data\nAccuracy: {knn.score(XX, yy_)}')
plt.show()Ignoring fixed x limits to fulfill fixed data aspect with adjustable data limits.
Ignoring fixed x limits to fulfill fixed data aspect with adjustable data limits.
The Iris dataset is a famous dataset used by Ronald Fisher in a classic 1936 paper on classification.
To set the hyperparameter \(k\) we follow the protocol above: score each \(k\) on the validation split of many random train/validation splits.
X = iris.data
y = iris.target
kvals = range(2, 20)
nreps = 50
acc = []
std = []
np.random.seed(0)
for k in kvals:
val_rep = []
train_rep = []
for i in range(nreps):
X_train, X_val, y_train, y_val = model_selection.train_test_split(
X, y, test_size = 0.33)
knn = KNeighborsClassifier(n_neighbors = k)
knn.fit(X_train, y_train)
train_rep.append(knn.score(X_train, y_train))
val_rep.append(knn.score(X_val, y_val))
acc.append([np.mean(val_rep), np.mean(train_rep)])
std.append([np.std(val_rep), np.std(train_rep)])plt.figure(figsize= (6, 4))
accy = np.array(acc)
stds = np.array(std)/np.sqrt(nreps)
print(f'Max validation accuracy at k = {kvals[np.argmax(accy[:, 0])]} with accuracy {np.max(accy[:, 0]):.03f}')
plt.plot(kvals, accy[:, 0], '.-', label = 'Accuracy on Validation Data')
plt.plot(kvals, accy[:, 1], '.-', label = 'Accuracy on Training Data')
plt.xlabel('k')
plt.ylabel('Accuracy')
plt.ylim(0.9, None)
plt.xticks(kvals)
plt.legend(loc = 'best')
plt.show()Max validation accuracy at k = 13 with accuracy 0.969
# we will use only the first two (of four) features, so we can visualize
X2 = iris.data[:, :2]
h = .02 # step size in the mesh
k = 13
knn = KNeighborsClassifier(n_neighbors=k)
knn.fit(X2, y)
# Plot the decision boundary. For that, we will assign a color to each
# point in the mesh [x_min, x_max]x[y_min, y_max].
x_min, x_max = X2[:, 0].min() - 1, X2[:, 0].max() + 1
y_min, y_max = X2[:, 1].min() - 1, X2[:, 1].max() + 1
xx, yy = np.meshgrid(np.arange(x_min, x_max, h),
np.arange(y_min, y_max, h))
Z = knn.predict(np.c_[xx.ravel(), yy.ravel()])
# Put the result into a color plot
Z = Z.reshape(xx.shape)
plt.figure(figsize=(6, 4))
plt.pcolormesh(xx, yy, Z, cmap=cmap3_light, shading='auto')
# Plot also the training points
plt.scatter(X2[:, 0], X2[:, 1], c=y, cmap=cmap3_bold, edgecolors='k', linewidths=0.5)
plt.xlim(xx.min(), xx.max())
plt.ylim(yy.min(), yy.max())
plt.xlabel(iris.feature_names[0])
plt.ylabel(iris.feature_names[1])
plt.title(f"3-Class $k$-NN classification ($k$ = {k})")
plt.show()
There are a few artifacts, though overall this looks like a reasonably smooth set of decision boundaries.
Here we show the matrix as an image.
It is easier to visualize if we blur the pixels a little bit.
Here are some more samples from the dataset.
def sweep_k(kvals, val_fraction, nreps):
acc = []
std = []
np.random.seed(0)
#
for k in kvals:
val_rep = []
for i in range(nreps):
X_train, X_val, y_train, y_val = model_selection.train_test_split(
X, y, test_size = val_fraction)
knn = KNeighborsClassifier(n_neighbors = k)
knn.fit(X_train, y_train)
val_rep.append(knn.score(X_val, y_val))
acc.append(np.mean(val_rep))
std.append(np.std(val_rep))
return(np.array(acc), np.array(std)/np.sqrt(nreps))
val_fraction1 = 0.33
accy1, stds1 = sweep_k(kvals, val_fraction1, 50)
val_fraction2 = 0.10
accy2, stds2 = sweep_k(kvals, val_fraction2, 50)plt.figure(figsize = (6, 4))
plt.errorbar(kvals, accy1, stds1,
label = f'{val_fraction1:.0%} held out for validation; best accuracy {np.max(accy1):.03f}')
plt.errorbar(kvals, accy2, stds2,
label = f'{val_fraction2:.0%} held out for validation; best accuracy {np.max(accy2):.03f}')
plt.xlabel('k')
plt.ylabel('Accuracy')
plt.legend(loc = 'best')
plt.xticks(kvals)
plt.title(f'Validation Accuracy with $\\pm 1\\sigma$ Error Bars')
plt.show()
We can get a sense of why \(k\)-NN can succeed at this task by looking at the nearest neighbors of some points.
We hold the first three images out as queries and fit on the rest, so that a query cannot be its own nearest neighbour:
Query labels: [1 5 0]
Neighbour labels: [1 1 1] [5 5 5] [0 0 0]plt.rc("image", cmap="binary") # this sets a black on white colormap
for t in range(3):
plt.figure(figsize=(8, 2))
plt.subplot(1, 4, 1)
plt.imshow(X[t].reshape(8, 8))
plt.axis('off')
plt.title(f"Query (label {y[t]})")
# plot three nearest neighbors from the training set
for i in [0, 1, 2]:
plt.subplot(1, 4, 2 + i)
plt.title(f"neighbor {i + 1} (label {y[neighbors[t, i]]})")
plt.imshow(X[neighbors[t, i]].reshape(8, 8))
plt.axis('off')
plt.show()
Classify a new taxpayer with Refund = No, Married (ignoring income).
\[ \begin{align*} Evade &= No: P(\text{Refund=No} \vert \text{Evade=No}) \cdot P(\text{Married} \vert \text{Evade=No}) \cdot P(\text{Evade=No}) \\ Evade &=No: \tfrac{4}{7} \cdot \tfrac{4}{7} \cdot 0.7 \approx 0.23 \end{align*} \]
\[ \begin{align*} Evade &=Yes: P(\text{Refund=No} \vert \text{Evade=Yes}) \cdot P(\text{Married} \vert \text{Evade=Yes}) \cdot P(\text{Evade=Yes}) \\ Evade &=Yes: \tfrac{3}{3} \cdot \tfrac{0}{3} \cdot 0.3 = 0 \end{align*} \]
so we predict Evade = No.
Note the zero: one empty histogram bin wipes out every other factor. A single value never seen with a class vetoes that class, however strong the rest of the evidence.