In [1]:
__AUTHORS__ = {'am': ("Andrea Marino",
"andrea.marino@unifi.it",),
'mn': ("Massimo Nocentini",
"massimo.nocentini@unifi.it",
"https://github.com/massimo-nocentini/",)}
__KEYWORDS__ = ['Python', 'numpy', 'numerical', 'data',]
The topic is very broad: datasets can come from a wide range of sources and a wide range of formats, including be collections of documents, collections of images, collections of sound clips, collections of numerical measurements, or nearly anything else. Despite this apparent heterogeneity, it will help us to think of all data fundamentally as arrays of numbers.
For this reason, efficient storage and manipulation of numerical arrays is absolutely fundamental to the process of doing data science.
NumPy (short for Numerical Python) provides an efficient interface to store and operate on dense data buffers.
In some ways, NumPy arrays are like Python's built-in list type, but NumPy arrays provide much more efficient storage and data operations as the arrays grow larger in size.
NumPy arrays form the core of nearly the entire ecosystem of data science tools in Python, so time spent learning to use NumPy effectively will be valuable no matter what aspect of data science interests you.
In [2]:
import numpy
numpy.__version__
Out[2]:
By convention, you'll find that most people in the SciPy/PyData world will import NumPy using np as an alias:
In [1]:
import numpy as np
Throughout this chapter, and indeed the rest of the book, you'll find that this is the way we will import and use NumPy.
Effective data-driven science and computation requires understanding how data is stored and manipulated.
Here we outlines and contrasts how arrays of data are handled in the Python language itself, and how NumPy improves on this.
Python offers several different options for storing data in efficient, fixed-type data buffers.
The built-in array module (available since Python 3.3) can be used to create dense arrays of a uniform type:
In [5]:
import array
L = list(range(10))
A = array.array('i', L)
A
Out[5]:
Here 'i' is a type code indicating the contents are integers.
Much more useful, however, is the ndarray object of the NumPy package.
While Python's array object provides efficient storage of array-based data, NumPy adds to this efficient operations on that data.
In [10]:
np.array([1, 4, 2, 5, 3])
Out[10]:
Remember that unlike Python lists, NumPy is constrained to arrays that all contain the same type. If types do not match, NumPy will upcast if possible (here, integers are up-cast to floating point):
In [8]:
np.array([3.14, 4, 2, 3])
Out[8]:
If we want to explicitly set the data type of the resulting array, we can use the dtype keyword:
In [9]:
np.array([1, 2, 3, 4], dtype='float32')
Out[9]:
In [11]:
np.zeros(10, dtype=int)
Out[11]:
In [12]:
np.ones((3, 5), dtype=float)
Out[12]:
In [13]:
np.full((3, 5), 3.14)
Out[13]:
In [14]:
np.arange(0, 20, 2)
Out[14]:
In [15]:
np.linspace(0, 1, 5)
Out[15]:
In [16]:
np.random.random((3, 3))
Out[16]:
In [17]:
np.random.normal(0, 1, (3, 3))
Out[17]:
In [18]:
np.eye(3)
Out[18]:
NumPy arrays contain values of a single type, so have a look at those types and their bounds:
| Data type | Description |
|---|---|
bool_ |
Boolean (True or False) stored as a byte |
int_ |
Default integer type (same as C long; normally either int64 or int32) |
intc |
Identical to C int (normally int32 or int64) |
intp |
Integer used for indexing (same as C ssize_t; normally either int32 or int64) |
int8 |
Byte (-128 to 127) |
int16 |
Integer (-32768 to 32767) |
int32 |
Integer (-2147483648 to 2147483647) |
int64 |
Integer (-9223372036854775808 to 9223372036854775807) |
uint8 |
Unsigned integer (0 to 255) |
uint16 |
Unsigned integer (0 to 65535) |
uint32 |
Unsigned integer (0 to 4294967295) |
uint64 |
Unsigned integer (0 to 18446744073709551615) |
float_ |
Shorthand for float64. |
float16 |
Half precision float: sign bit, 5 bits exponent, 10 bits mantissa |
float32 |
Single precision float: sign bit, 8 bits exponent, 23 bits mantissa |
float64 |
Double precision float: sign bit, 11 bits exponent, 52 bits mantissa |
complex_ |
Shorthand for complex128. |
complex64 |
Complex number, represented by two 32-bit floats |
complex128 |
Complex number, represented by two 64-bit floats |
Data manipulation in Python is nearly synonymous with NumPy array manipulation: even newer tools like Pandas are built around the NumPy array.
In [19]:
np.random.seed(0) # seed for reproducibility
x1 = np.random.randint(10, size=6) # One-dimensional array
x2 = np.random.randint(10, size=(3, 4)) # Two-dimensional array
x3 = np.random.randint(10, size=(3, 4, 5)) # Three-dimensional array
Each array has attributes ndim (the number of dimensions), shape (the size of each dimension), size (the total size of the array) and dtype (the data type of the array):
In [21]:
print("x3 ndim: ", x3.ndim)
print("x3 shape:", x3.shape)
print("x3 size: ", x3.size)
print("dtype:", x3.dtype)
In [22]:
x1
Out[22]:
In [23]:
x1[0]
Out[23]:
In [25]:
x1[-1] # To index from the end of the array, you can use negative indices.
Out[25]:
In a multi-dimensional array, items can be accessed using a comma-separated tuple of indices:
In [26]:
x2
Out[26]:
In [27]:
x2[0, 0]
Out[27]:
In [28]:
x2[2, -1]
Out[28]:
Values can also be modified using any of the above index notation:
In [29]:
x2[0, 0] = 12
x2
Out[29]:
Keep in mind that, unlike Python lists, NumPy arrays have a fixed type.
In [30]:
x1[0] = 3.14159 # this will be truncated!
x1
Out[30]:
Just as we can use square brackets to access individual array elements, we can also use them to access subarrays with the slice notation, marked by the colon (:) character.
The NumPy slicing syntax follows that of the standard Python list; to access a slice of an array x, use this:
x[start:stop:step]
If any of these are unspecified, they default to the values start=0, stop=size of dimension, step=1.
In [31]:
x = np.arange(10)
x
Out[31]:
In [32]:
x[:5] # first five elements
Out[32]:
In [33]:
x[5:] # elements after index 5
Out[33]:
In [34]:
x[4:7] # middle sub-array
Out[34]:
In [35]:
x[::2] # every other element
Out[35]:
In [36]:
x[1::2] # every other element, starting at index 1
Out[36]:
A potentially confusing case is when the step value is negative.
In this case, the defaults for start and stop are swapped.
This becomes a convenient way to reverse an array:
In [37]:
x[::-1] # all elements, reversed
Out[37]:
In [38]:
x[5::-2] # reversed every other from index 5
Out[38]:
In [39]:
x2
Out[39]:
In [40]:
x2[:2, :3] # two rows, three columns
Out[40]:
In [41]:
x2[:3, ::2] # all rows, every other column
Out[41]:
In [42]:
x2[::-1, ::-1]
Out[42]:
In [43]:
print(x2[:, 0]) # first column of x2
In [44]:
print(x2[0, :]) # first row of x2
In [45]:
print(x2[0]) # equivalent to x2[0, :]
In [46]:
x2
Out[46]:
In [47]:
x2_sub = x2[:2, :2]
x2_sub
Out[47]:
In [49]:
x2_sub[0, 0] = 99 # if we modify this subarray, the original array is changed too
x2
Out[49]:
It is sometimes useful to instead explicitly copy the data within an array or a subarray. This can be most easily done with the copy() method.
In [110]:
np.arange(1, 10).reshape((3, 3))
Out[110]:
In [111]:
x = np.array([1, 2, 3])
x.reshape((1, 3)) # row vector via reshape
Out[111]:
In [112]:
x[np.newaxis, :] # row vector via newaxis
Out[112]:
In [113]:
x.reshape((3, 1)) # column vector via reshape
Out[113]:
In [114]:
x[:, np.newaxis] # column vector via newaxis
Out[114]:
In [105]:
x = np.array([1, 2, 3])
y = np.array([3, 2, 1])
np.concatenate([x, y])
Out[105]:
In [106]:
z = [99, 99, 99]
np.concatenate([x, y, z])
Out[106]:
In [107]:
grid = np.array([[1, 2, 3],
[4, 5, 6]])
In [108]:
np.concatenate([grid, grid]) # concatenate along the first axis
Out[108]:
In [109]:
np.concatenate([grid, grid], axis=1) # concatenate along the second axis (zero-indexed)
Out[109]:
For working with arrays of mixed dimensions, it can be clearer to use the np.vstack (vertical stack) and np.hstack (horizontal stack) functions:
In [63]:
x = np.array([1, 2, 3])
grid = np.array([[9, 8, 7],
[6, 5, 4]])
np.vstack([x, grid]) # vertically stack the arrays
Out[63]:
In [64]:
y = np.array([[99],
[99]])
np.hstack([grid, y]) # horizontally stack the arrays
Out[64]:
In [65]:
x = [1, 2, 3, 99, 99, 3, 2, 1]
x1, x2, x3 = np.split(x, [3, 5])
print(x1, x2, x3)
In [66]:
grid = np.arange(16).reshape((4, 4))
grid
Out[66]:
In [67]:
np.vsplit(grid, [2])
Out[67]:
In [68]:
np.hsplit(grid, [2])
Out[68]:
Python's default implementation (known as CPython) does some operations very slowly, this is in part due to the dynamic, interpreted nature of the language.
The relative sluggishness of Python generally manifests itself in situations where many small operations are being repeated – for instance looping over arrays to operate on each element.
For example, pretend to compute the reciprocal of values contained in a array:
In [69]:
np.random.seed(0)
def compute_reciprocals(values):
output = np.empty(len(values))
for i in range(len(values)):
output[i] = 1.0 / values[i]
return output
values = np.random.randint(1, 10, size=5)
compute_reciprocals(values)
Out[69]:
If we measure the execution time of this code for a large input, we see that this operation is very slow, perhaps surprisingly so!
In [70]:
big_array = np.random.randint(1, 100, size=1000000)
%timeit compute_reciprocals(big_array)
It takes $2.63$ seconds to compute these million operations and to store the result.
It turns out that the bottleneck here is not the operations themselves, but the type-checking and function dispatches that CPython must do at each cycle of the loop.
If we were working in compiled code instead, this type specification would be known before the code executes and the result could be computed much more efficiently.
In [72]:
%timeit (1.0 / big_array)
Vectorized operations in NumPy are implemented via ufuncs, whose main purpose is to quickly execute repeated operations on values in NumPy arrays.
Ufuncs are extremely flexible – before we saw an operation between a scalar and an array, but we can also operate between two arrays:
In [75]:
np.arange(5) / np.arange(1, 6)
Out[75]:
And ufunc operations are not limited to one-dimensional arrays–they can also act on multi-dimensional arrays as well:
In [76]:
x = np.arange(9).reshape((3, 3))
2 ** x
Out[76]:
Any time you see such a loop in a Python script, you should consider whether it can be replaced with a vectorized expression.
In [78]:
x = np.arange(4)
print("x =", x)
print("x + 5 =", x + 5)
print("x - 5 =", x - 5)
print("x * 2 =", x * 2)
print("x / 2 =", x / 2)
print("x // 2 =", x // 2) # floor division
print("-x = ", -x)
print("x ** 2 = ", x ** 2)
print("x % 2 = ", x % 2)
In [80]:
-(0.5*x + 1) ** 2 # can be strung together also
Out[80]:
In [81]:
theta = np.linspace(0, np.pi, 3)
In [83]:
print("theta = ", theta)
print("sin(theta) = ", np.sin(theta))
print("cos(theta) = ", np.cos(theta))
print("tan(theta) = ", np.tan(theta))
In [84]:
x = [1, 2, 3]
print("x =", x)
print("e^x =", np.exp(x))
print("2^x =", np.exp2(x))
print("3^x =", np.power(3, x))
In [86]:
x = [1, 2, 4, 10]
print("x =", x)
print("ln(x) =", np.log(x))
print("log2(x) =", np.log2(x))
print("log10(x) =", np.log10(x))
In [87]:
x = np.arange(5)
y = np.empty(5)
np.multiply(x, 10, out=y)
print(y)
In [88]:
y = np.zeros(10)
np.power(2, x, out=y[::2])
print(y)
In [89]:
x = np.arange(1, 6)
np.multiply.outer(x, x)
Out[89]:
In [90]:
L = np.random.random(100)
sum(L)
Out[90]:
In [91]:
np.sum(L)
Out[91]:
In [92]:
big_array = np.random.rand(1000000)
%timeit sum(big_array)
%timeit np.sum(big_array)
In [93]:
min(big_array), max(big_array)
Out[93]:
In [94]:
np.min(big_array), np.max(big_array)
Out[94]:
In [95]:
%timeit min(big_array)
%timeit np.min(big_array)
In [97]:
big_array.min(), big_array.max(), big_array.sum()
Out[97]:
In [99]:
M = np.random.random((3, 4))
M
Out[99]:
In [100]:
M.sum() # By default, each NumPy aggregation function works on the whole array
Out[100]:
In [101]:
M.min(axis=0) # specifying the axis along which the aggregate is computed
Out[101]:
In [102]:
M.max(axis=1) # find the maximum value within each row
Out[102]:
Additionally, most aggregates have a NaN-safe counterpart that computes the result while ignoring missing values, which are marked by the special IEEE floating-point NaN value
| Function Name | NaN-safe Version | Description |
|---|---|---|
np.sum |
np.nansum |
Compute sum of elements |
np.prod |
np.nanprod |
Compute product of elements |
np.mean |
np.nanmean |
Compute mean of elements |
np.std |
np.nanstd |
Compute standard deviation |
np.var |
np.nanvar |
Compute variance |
np.min |
np.nanmin |
Find minimum value |
np.max |
np.nanmax |
Find maximum value |
np.argmin |
np.nanargmin |
Find index of minimum value |
np.argmax |
np.nanargmax |
Find index of maximum value |
np.median |
np.nanmedian |
Compute median of elements |
np.percentile |
np.nanpercentile |
Compute rank-based statistics of elements |
np.any |
N/A | Evaluate whether any elements are true |
np.all |
N/A | Evaluate whether all elements are true |
In [2]:
a = np.array([0, 1, 2])
b = np.array([5, 5, 5])
a + b
Out[2]:
Broadcasting allows these types of binary operations to be performed on arrays of different sizes:
In [3]:
a + 5
Out[3]:
We can think of this as an operation that stretches or duplicates the value 5 into the array [5, 5, 5], and adds the results; the advantage of NumPy's broadcasting is that this duplication of values does not actually take place.
We can similarly extend this to arrays of higher dimensions:
In [4]:
M = np.ones((3, 3))
M
Out[4]:
In [5]:
M + a
Out[5]:
Here the one-dimensional array a is stretched, or broadcast across the second dimension in order to match the shape of M.
More complicated cases can involve broadcasting of both arrays:
In [6]:
a = np.arange(3)
b = np.arange(3)[:, np.newaxis]
a, b
Out[6]:
In [7]:
a + b
Out[7]:
Broadcasting in NumPy follows a strict set of rules to determine the interaction between the two arrays:
In [8]:
X = np.random.random((10, 3))
In [9]:
Xmean = X.mean(0)
Xmean
Out[9]:
In [10]:
X_centered = X - Xmean
In [12]:
X_centered.mean(0) # To double-check, we can check that the centered array has near 0 means.
Out[12]:
In [117]:
steps = 500
x = np.linspace(0, 5, steps) # # x and y have 500 steps from 0 to 5
y = np.linspace(0, 5, steps)[:, np.newaxis]
z = np.sin(x) ** 10 + np.cos(10 + y * x) * np.cos(x)
In [116]:
%matplotlib inline
import matplotlib.pyplot as plt
plt.imshow(z, origin='lower', extent=[0, 5, 0, 5], cmap='viridis')
plt.colorbar();
Masking comes up when you want to extract, modify, count, or otherwise manipulate values in an array based on some criterion: for example, you might wish to count all values greater than a certain value, or perhaps remove all outliers that are above some threshold. In NumPy, Boolean masking is often the most efficient way to accomplish these types of tasks.
In [15]:
x = np.array([1, 2, 3, 4, 5])
In [16]:
x < 3 # less than
Out[16]:
In [17]:
x > 3 # greater than
Out[17]:
In [18]:
x != 3 # not equal
Out[18]:
In [19]:
(2 * x) == (x ** 2)
Out[19]:
Just as in the case of arithmetic ufuncs, these will work on arrays of any size and shape:
In [20]:
rng = np.random.RandomState(0)
x = rng.randint(10, size=(3, 4))
x
Out[20]:
In [21]:
x < 6
Out[21]:
In [22]:
np.count_nonzero(x < 6) # how many values less than 6?
Out[22]:
In [23]:
np.sum(x < 6)
Out[23]:
In [24]:
np.sum(x < 6, axis=1) # how many values less than 6 in each row?
Out[24]:
In [26]:
np.any(x > 8) # are there any values greater than 8?
Out[26]:
In [25]:
np.any(x < 0) # are there any values less than zero?
Out[25]:
In [27]:
np.all(x < 10) # are all values less than 10?
Out[27]:
In [28]:
np.all(x < 8, axis=1) # are all values in each row less than 8?
Out[28]:
In [29]:
x
Out[29]:
In [30]:
x < 5
Out[30]:
In [31]:
x[x < 5]
Out[31]:
What is returned is a one-dimensional array filled with all the values that meet this condition; in other words, all the values in positions at which the mask array is True.
We saw how to access and modify portions of arrays using simple indices (e.g., arr[0]), slices (e.g., arr[:5]), and Boolean masks (e.g., arr[arr > 0]).
We'll look at another style of array indexing, known as fancy indexing, that is like the simple indexing we've already seen, but we pass arrays of indices in place of single scalars.
Fancy indexing is conceptually simple: it means passing an array of indices to access multiple array elements at once:
In [32]:
rand = np.random.RandomState(42)
x = rand.randint(100, size=10)
x
Out[32]:
In [33]:
[x[3], x[7], x[2]] # Suppose we want to access three different elements.
Out[33]:
In [34]:
ind = [3, 7, 4]
x[ind] # Alternatively, we can pass a single list or array of indices
Out[34]:
When using fancy indexing, the shape of the result reflects the shape of the index arrays rather than the shape of the array being indexed:
In [35]:
ind = np.array([[3, 7],
[4, 5]])
x[ind]
Out[35]:
Fancy indexing also works in multiple dimensions:
In [36]:
X = np.arange(12).reshape((3, 4))
X
Out[36]:
Like with standard indexing, the first index refers to the row, and the second to the column:
In [37]:
row = np.array([0, 1, 2])
col = np.array([2, 1, 3])
X[row, col]
Out[37]:
The pairing of indices in fancy indexing follows all the broadcasting rules that we've already seen:
In [38]:
X[row[:, np.newaxis], col]
Out[38]:
each row value is matched with each column vector, exactly as we saw in broadcasting of arithmetic operations
In [39]:
row[:, np.newaxis] * col
Out[39]:
Remember: with fancy indexing that the return value reflects the broadcasted shape of the indices, rather than the shape of the array being indexed.
In [41]:
X
Out[41]:
In [42]:
X[2, [2, 0, 1]] # combine fancy and simple indices
Out[42]:
In [43]:
X[1:, [2, 0, 1]] # combine fancy indexing with slicing
Out[43]:
In [44]:
mask = np.array([1, 0, 1, 0], dtype=bool)
X[row[:, np.newaxis], mask] # combine fancy indexing with masking
Out[44]:
In [119]:
mean = [0, 0]
cov = [[1, 2],
[2, 5]]
X = rand.multivariate_normal(mean, cov, 100)
X.shape
Out[119]:
In [120]:
plt.scatter(X[:, 0], X[:, 1]);
Let's use fancy indexing to select 20 random points. We'll do this by first choosing 20 random indices with no repeats, and use these indices to select a portion of the original array:
In [48]:
indices = np.random.choice(X.shape[0], 20, replace=False)
indices
Out[48]:
In [49]:
selection = X[indices] # fancy indexing here
selection.shape
Out[49]:
Now to see which points were selected, let's over-plot large circles at the locations of the selected points:
In [52]:
plt.scatter(X[:, 0], X[:, 1], alpha=0.3);
In [54]:
x = np.arange(10)
i = np.array([2, 1, 8, 4])
x[i] = 99
x
Out[54]:
In [55]:
x[i] -= 10 # use any assignment-type operator for this
x
Out[55]:
Notice, though, that repeated indices with these operations can cause some potentially unexpected results:
In [56]:
x = np.zeros(10)
x[[0, 0]] = [4, 6]
x
Out[56]:
Where did the 4 go? The result of this operation is to first assign x[0] = 4, followed by x[0] = 6.
The result, of course, is that x[0] contains the value 6.
In [57]:
i = [2, 3, 3, 4, 4, 4]
x[i] += 1
x
Out[57]:
You might expect that x[3] would contain the value 2, and x[4] would contain the value 3, as this is how many times each index is repeated. Why is this not the case?
Conceptually, this is because x[i] += 1 is meant as a shorthand of x[i] = x[i] + 1. x[i] + 1 is evaluated, and then the result is assigned to the indices in x.
With this in mind, it is not the augmentation that happens multiple times, but the assignment, which leads to the rather nonintuitive results.
In [59]:
x = np.zeros(10)
np.add.at(x, i, 1)
x
Out[59]:
The at() method does an in-place application of the given operator at the specified indices (here, i) with the specified value (here, 1).
Another method that is similar in spirit is the reduceat() method of ufuncs, which you can read about in the NumPy documentation.
In [121]:
np.random.seed(42)
x = np.random.randn(100)
# compute a histogram by hand
bins = np.linspace(-5, 5, 20)
counts = np.zeros_like(bins)
# find the appropriate bin for each x
i = np.searchsorted(bins, x)
# add 1 to each of these bins
np.add.at(counts, i, 1)
In [131]:
# The counts now reflect the number of points
# within each bin–in other words, a histogram:
line, = plt.plot(bins, counts);
line.set_drawstyle("steps")
In [65]:
print("NumPy routine:")
%timeit counts, edges = np.histogram(x, bins)
print("Custom routine:")
%timeit np.add.at(counts, np.searchsorted(bins, x), 1)
Our own one-line algorithm is several times faster than the optimized algorithm in NumPy! How can this be?
If you dig into the np.histogram source code (you can do this in IPython by typing np.histogram??), you'll see that it's quite a bit more involved than the simple search-and-count that we've done; this is because NumPy's algorithm is more flexible, and particularly is designed for better performance when the number of data points becomes large...
In [66]:
x = np.random.randn(1000000)
print("NumPy routine:")
%timeit counts, edges = np.histogram(x, bins)
print("Custom routine:")
%timeit np.add.at(counts, np.searchsorted(bins, x), 1)
What this comparison shows is that algorithmic efficiency is almost never a simple question. An algorithm efficient for large datasets will not always be the best choice for small datasets, and vice versa.
The key to efficiently using Python in data-intensive applications is knowing about general convenience routines like np.histogram and when they're appropriate, but also knowing how to make use of lower-level functionality when you need more pointed behavior.
In [70]:
x = np.array([2, 1, 4, 3, 5])
np.sort(x)
Out[70]:
In [71]:
x
Out[71]:
A related function is argsort, which instead returns the indices of the sorted elements:
In [72]:
i = np.argsort(x)
i
Out[72]:
The first element of this result gives the index of the smallest element, the second value gives the index of the second smallest, and so on.
These indices can then be used (via fancy indexing) to construct the sorted array if desired:
In [73]:
x[i]
Out[73]:
In [75]:
rand = np.random.RandomState(42)
X = rand.randint(0, 10, (4, 6))
X
Out[75]:
In [76]:
np.sort(X, axis=0) # sort each column of X
Out[76]:
In [77]:
np.sort(X, axis=1) # sort each row of X
Out[77]:
Keep in mind that this treats each row or column as an independent array, and any relationships between the row or column values will be lost!
Sometimes we're not interested in sorting the entire array, but simply want to find the k smallest values in the array. np.partition takes an array and a number K; the result is a new array with the smallest K values to the left of the partition, and the remaining values to the right, in arbitrary order:
In [78]:
x = np.array([7, 2, 3, 1, 6, 5, 4])
np.partition(x, 3)
Out[78]:
Note that the first three values in the resulting array are the three smallest in the array, and the remaining array positions contain the remaining values.
Within the two partitions, the elements have arbitrary order.
Similarly to sorting, we can partition along an arbitrary axis of a multidimensional array:
In [79]:
np.partition(X, 2, axis=1)
Out[79]:
The result is an array where the first two slots in each row contain the smallest values from that row, with the remaining values filling the remaining slots.
Finally, just as there is a np.argsort that computes indices of the sort, there is a np.argpartition that computes indices of the partition.
In [85]:
X = rand.rand(50, 2)
In [86]:
plt.scatter(X[:, 0], X[:, 1], s=100);
In [95]:
# compute the distance between each pair of points
dist_sq = np.sum((X[:, np.newaxis, :] - X[np.newaxis, :, :]) ** 2, axis=-1)
dist_sq.shape, np.all(dist_sq.diagonal() == 0)
Out[95]:
With the pairwise square-distances converted, we can now use np.argsort to sort along each row.
The leftmost columns will then give the indices of the nearest neighbors:
In [98]:
nearest = np.argsort(dist_sq, axis=1)
nearest[:,0]
Out[98]:
Notice that the first column is order because each point's closest neighbor is itself.
If we're simply interested in the nearest $k$ neighbors, all we need is to partition each row so that the smallest $k + 1$ squared distances come first, with larger distances filling the remaining positions of the array:
In [99]:
K = 2
nearest_partition = np.argpartition(dist_sq, K + 1, axis=1)
In [104]:
plt.scatter(X[:, 0], X[:, 1], s=100)
K = 2 # draw lines from each point to its two nearest neighbors
for i in range(X.shape[0]):
for j in nearest_partition[i, :K+1]:
plt.plot(*zip(X[j], X[i]), color='black')
At first glance, it might seem strange that some of the points have more than two lines coming out of them: this is due to the fact that if point A is one of the two nearest neighbors of point B, this does not necessarily imply that point B is one of the two nearest neighbors of point A.
You might be tempted to do the same type of operation by manually looping through the data and sorting each set of neighbors individually. The beauty of our approach is that it's written in a way that's agnostic to the size of the input data: we could just as easily compute the neighbors among 100 or 1,000,000 points in any number of dimensions, and the code would look the same.