So in the first part of this series, I covered the very basics of using Numpy arrays. This post will be about more advanced and elaborate constructs possible with ndarrays.
Broadcasting
Arrays which almost have the same shape can be molded into different shapes for the ease of using ufuncs, so that all arrays passed to the ufunc have the same size. Here are the rules that numpy uses to decided how an array is broadcasted:
- If all arrays don't have the same number of dimensions, dimensions of size 1 are appended to the smaller-dimensional array.
- Any dimensions having a length of 1, the rest of the dimensions are propagated along that 1-length dimension. So the 1-length dimension has entirely identical elements (this is the operation that's called broadcasting, and rule 1 is actually a special case of this).
After using these rules, the size of both arrays need to match. In other words, broadcasting can only help you expand 1-length dimensions. It can't reshape a 2x6 array into some other shape for example.
There isn't a specific function for broadcasting, broadcasting is done implicitly by ordinary operations. It usually involves combining an array with smaller dimensions with a greater-dimensional array. The vectorization that happens as a result is performed by fast compiled C code.
This is an example of an operation with no broadcasting:
>>> a = np.array([1.0, 2.0, 3.0])
>>> b = np.array([2.0, 2.0, 2.0])
>>> a * b
array([ 2., 4., 6.])
This is what a simple operation with broadcasting would look like:
>>> a = np.array([1.0, 2.0, 3.0])
>>> b = 2.0
>>> a * b
array([ 2., 4., 6.])
So as you can see, it didn't matter whether b was [2, 2, 2] or just 2, because the latter was expanded (broadcasted) into the former.
An invalid broadcasting attempt will throw ValueError: operands could not be broadcast together.
As I mentioned, broadcasting works with arrays of any size as long as one array has fewer dimensions than the other. Here's a practical example from the numpy documentation showing the need for broadcasting for one-dimensional array to three-dimensional arrays.
Image (3d array): 256 x 256 x 3
Scale (1d array): 3
Result (3d array): 256 x 256 x 3
In this example each of the image colors (Red, Green and Blue, the length-3 dimension) are scaled by the corresponding value in a three-element one-dimensional array, which would look something like [2.0, 1.3, 2.2] for example.
Perhaps this example will help visualize the broadcasting operation better:
A (4d array): 8 x 1 x 6 x 1
B (3d array): 7 x 1 x 5
Result (4d array): 8 x 7 x 6 x 5
You can see that the length-1 dimensions are cancelled out here. Extra dimensions are added to the smaller array if necessary. This is the application of both of the broadcasting rules in practice.
Another practical example of broadcasting is taking the outer product of two arrays:
>>> a = np.array([0.0, 10.0, 20.0, 30.0])
>>> b = np.array([1.0, 2.0, 3.0])
>>> a[:, np.newaxis] + b
array([[ 1., 2., 3.],
[ 11., 12., 13.],
[ 21., 22., 23.],
[ 31., 32., 33.]])
np.newaxis added another dimension into a, creating a two-dimensional array, when it used to be one-dimensional. Because we didn't assign this to a, with =, it is still one-dimensional.
More indexing - Integer arrays
You can index numpy arrays with arrays of integers and arrays of boolean values. Other scientific languages such as Matlab have this capability and the reason why someone would want to do that is if they want to put the same index multiple times in the same slice.
To talk about indexing with arrays of integers first:
>>> a = np.arange(12)**2 # the first 12 square numbers
>>> i = np.array( [ 1,1,3,8,5 ] ) # an array of indices
>>> a[i] # the elements of a at the positions i
array([ 1, 1, 9, 64, 25])
>>>
>>> j = np.array( [ [ 3, 4], [ 9, 7 ] ] ) # a bidimensional array of indices
>>> a[j] # the same shape as j
array([[ 9, 16],
[81, 49]])
So you see, indexing [1] several times will retrieve the position in that index multiple times, that's where the first two 1's in the output come from. From there you can shape the array any way you like.
You can influence the dimensionality of the resulting array by putting more dimensions in the index array:
>>> palette = np.array( [ [0,0,0], # black
... [255,0,0], # red
... [0,255,0], # green
... [0,0,255], # blue
... [255,255,255] ] ) # white
>>> image = np.array( [ [ 0, 1, 2, 0 ], # each value corresponds to a color in the palette
... [ 0, 3, 4, 0 ] ] )
>>> palette[image] # the (2,4,3) color image
array([[[ 0, 0, 0],
[255, 0, 0],
[ 0, 255, 0],
[ 0, 0, 0]],
[[ 0, 0, 0],
[ 0, 0, 255],
[255, 255, 255],
[ 0, 0, 0]]])
*Each index refers to elements in the array's first dimension. In the above example the first dimension of palette is rows as you can clearly see whole arrays like [0, 0, 0] are being indexed. The array palette is a two-dimensional array. Here's an example that indexes a 3D array using this method:
In [1]: palette = np.array( [ [ [0,0,0], # black
...: [255,0,0], # red
...: [0,255,0], # green
...: [0,0,255], # blue
...: [255,255,255] ],
...: [ [1,1,1], # black
...: [255,1,1], # red
...: [1,255,1], # green
...: [1,1,255], # blue
...: [255,255,255] ] ]) # white
In [2]: image = np.array( [ [ 0, 1 ], # each value corresponds to a 2D array across depth dimension
...: [ 1, 1 ] ] )
In [3]: palette[image]
Out[3]:
array([[[[ 0, 0, 0],
[255, 0, 0],
[ 0, 255, 0],
[ 0, 0, 255],
[255, 255, 255]],
[[ 1, 1, 1],
[255, 1, 1],
[ 1, 255, 1],
[ 1, 1, 255],
[255, 255, 255]]],
[[[ 1, 1, 1],
[255, 1, 1],
[ 1, 255, 1],
[ 1, 1, 255],
[255, 255, 255]],
[[ 1, 1, 1],
[255, 1, 1],
[ 1, 255, 1],
[ 1, 1, 255],
[255, 255, 255]]]])
As you can see the first dimension is now the depth dimension so that's where the indices refer to, and is also why large 2D arrays are being indexed this time around.
Good ol' matrix indexing
Did you think I would forget this one? 😉 Obviously an array for a numerical library must have matrix-like indexing support. And there is such a feature, it is a[i, j] where i and j are numbers, or integer arrays (both should have the same shape). By no means is this limited to two-dimension indexing. On the contrary, you use this indexing format with arrays of any dimension, like a[i1, i2, i3, i4, ..., in], with integer arrays support. Again, if you use integer arrays to index these, they all must be the same size.
Slice/range notation (:) works with any of the indices.
More indexing examples
Hopefully by showing you more examples that are in numpy's documentation it will make it clear how integer array indexing works.
>>> a = np.arange(12).reshape(3,4)
>>> a
array([[ 0, 1, 2, 3],
[ 4, 5, 6, 7],
[ 8, 9, 10, 11]])
>>> i = np.array( [ [0,1], # indices for the first dim of a
... [1,2] ] )
>>> j = np.array( [ [2,1], # indices for the second dim
... [3,3] ] )
>>>
>>> a[i,j] # i and j must have equal shape
array([[ 2, 5],
[ 7, 11]])
>>>
>>> a[i,2]
array([[ 2, 6],
[ 6, 10]])
>>>
>>> a[:,j] # i.e., a[ : , j]
array([[[ 2, 1],
[ 3, 3]],
[[ 6, 5],
[ 7, 7]],
[[10, 9],
[11, 11]]])
Here, I can take an array of indices for the first dimension of a (i), and an array of indices for its second dimension (j), and get the array elements in any shape I want. This is a very flexible way of indexing.
We can also pass [i, j] as a list (but not an ndarray) and get the same results.
>>> l = [i,j]
>>> a[l] # equivalent to a[i,j]
array([[ 2, 5],
[ 7, 11]])
A few sections above, we noted that indexing using an ndarray will take elements the first dimension of the array only. Therefore the result of the following indexing operation is expected.
>>> s = np.array( [i,j] )
>>> a[s] # not what we want
Traceback (most recent call last):
File "<stdin>", line 1, in ?
IndexError: index (3) out of range (0<=index<=2) in dimension 0
>>>
>>> a[tuple(s)] # same as a[i,j]
array([[ 2, 5],
[ 7, 11]])
Notice how indexing worked properly when we created a tuple out of the ndarray. That's because the tuple constructor converted the ndarray into a tuple, which indexes equivalently to a list.
You can selectively assign a few elements in the array using integer arrays.
>>> a = np.arange(5)
>>> a
array([0, 1, 2, 3, 4])
# Don't put the same index multiple times or it will assign to that value
# repeatedly, old values that were assigned will be lost.
>>> a[[1,3,4]] = 0
>>> a
array([0, 0, 2, 0, 0])
Avoid using integer array assignments with operation-with-assignment operators such as +=, -=, *= and others, because they won't operate on duplicate indices multiple times. For example:
>>> a = np.arange(5)
>>> a[[0,0,2]]+=1 # Don't do this, += will ignore the second 0
>>> a
array([1, 1, 3, 3, 4])
More indexing - Boolean arrays
Instead of indexing my arrays of integer indices you can also create an array of the same shape as the one you are indexing, only with True and False values. The elements in the same positions as True will be selected, the ones in the same positions as False won't. The result is a one-dimensional array with the elements corresponding to True:
>>> a = np.arange(12).reshape(3,4)
>>> b = a > 4 # See what just happened there?
>>> b # b is a boolean with a's shape
array([[False, False, False, False],
[False, True, True, True],
[ True, True, True, True]])
>>> a[b] # 1d array with the selected elements
array([ 5, 6, 7, 8, 9, 10, 11])
>>> a[b] = 0 # All elements of 'a' higher than 4 become 0
>>> a
array([[0, 1, 2, 3],
[4, 0, 0, 0],
[0, 0, 0, 0]])
A single relational expression can create a boolean array. Therefore they are easier to construct than integer arrays.
To index with boolean arrays using the matrix indexing, you make a set of boolean 1D arrays, each the length of the axis it's going to index, and just index the array using these boolean arrays instead of numbers:
>>> a = np.arange(12).reshape(3,4)
>>> b1 = np.array([False,True,True]) # first dim selection
>>> b2 = np.array([True,False,True,False]) # second dim selection
>>>
>>> a[b1,:] # selecting rows
array([[ 4, 5, 6, 7],
[ 8, 9, 10, 11]])
>>>
>>> a[b1] # same thing
array([[ 4, 5, 6, 7],
[ 8, 9, 10, 11]])
>>>
>>> a[:,b2] # selecting columns
array([[ 0, 2],
[ 4, 6],
[ 8, 10]])
>>>
>>> a[b1,b2] # a weird thing to do
array([ 4, 10])
Don't be fooled by the last statement. Intuitively you would expect b1 and b2 with the shape they have, to select:
| True | False | True | False |
|---|---|---|---|
| False | ~0~ | ~1~ | ~2~ |
| True | 4 | ~5~ | 6 |
| True | 8 | ~9~ | 10 |
But actually it doesn't work like that. It will take the positions of the first True values in all the arrays, and create a matrix index out of those positions to select the first element, and then do the same for the second, third, etc. elements. So all of the 1D boolean arrays must have the same number of True values. This is what 1D boolean indexing a three-dimensional array looks like:
In [38]: a = np.arange(24).reshape(3,4,2)
In [40]: a
Out[40]:
array([[[ 0, 1],
[ 2, 3],
[ 4, 5],
[ 6, 7]],
[[ 8, 9],
[10, 11],
[12, 13],
[14, 15]],
[[16, 17],
[18, 19],
[20, 21],
[22, 23]]])
In [42]: a[np.array([True, True, False]), np.array([True, True, False, False]),
...: np.array([True, True])]
Out[42]: array([ 0, 11])
Again, this does not select specific dimension elements based on the corresponding 1D array. It uses the position of the nth True in all 1D arrays to fetch the nth value of the resulting array.
Creating "meshes" from 1D arrays using _ix()
If you have a set of vectors (specifically, 1D arrays) you need to do algebra on, and you will be working in as many dimensions as you have vectors, the _ix() function can help you stretch the shape of each vector into an n-dimensional array.
_ix() takes all of your vectors as its arguments (don't pass it as a list of vectors), and then places elements of each vector along the dimension number which is the same number as the position you put the vector in the argument list of _ix(). The rest of the dimensions are size 1.
What I mean is, suppose you have the following three vectors which you will buff into 3D arrays:
>>> a = np.array([2,3,4,5])
>>> b = np.array([8,5,4])
>>> c = np.array([5,4,6,8,3])
>>> ax,bx,cx = np.ix_(a,b,c) # _ix() returns a tuple of n-dimensional arrays
Then ax.shape will be (4, 1, 1), bx.shape will be (1, 3, 1) and cx.shape will be (1, 1, 5). You can then combine the arrays in 3D:
>>> result = ax+bx*cx
>>> result
array([[[42, 34, 50, 66, 26],
[27, 22, 32, 42, 17],
[22, 18, 26, 34, 14]],
[[43, 35, 51, 67, 27],
[28, 23, 33, 43, 18],
[23, 19, 27, 35, 15]],
[[44, 36, 52, 68, 28],
[29, 24, 34, 44, 19],
[24, 20, 28, 36, 16]],
[[45, 37, 53, 69, 29],
[30, 25, 35, 45, 20],
[25, 21, 29, 37, 17]]])
I hesitate to claim that _ix() broadcasted the vectors into resulting arrays because broadcasting is a process done implicitly by numpy, not by calling an explicit function for that, but you'll see numpy docs call it that anyway, so take their word for it.
Reduce functions
This post originally contained references to reduce functions, but nevertheless, I'd like to take a detour to explain what a reduce function is. A reduce function applies an operation or some other function to each value of the array (and the result of the previous reduce) until all of the elements have eventually been passed.
(A little bit of functional programming ahead...)
So for any operation, you could reduce it starting from the end of the array or from the beginning of the array. If the operation is associative then starting at either ends will produce the same result. There is also an initial value that is combined with the first folded element using this operation. That element is usually the identity element for which applying the operation on any other element returns that other element. For example, 0 is the identity element of addition (because anything plus 0 is itself), 1 is the identity element of multiplication and exponentiation, and there is also an identity element for matrix multiplication. I should make it clear at this point that starting from the beginning of the array is known as folding left and starting at the end is called folding right.
In Python (3), functions can be reduced using functools.reduce(func, list, initval) folding left, or functools.reduce(lambda x,y: func(y,x), reversed(list), initval) folding right.
Back to _ix()
ufuncs have their own way of defining the reduce process:
# No functools here
# Fold left
>>> def ufunc_reduce(ufct, *vectors):
... vs = np.ix_(*vectors)
... r = ufct.identity
... for v in vs:
... r = ufct(r,v)
... return r
>>> ufunc_reduce(np.add,a,b,c)
array([[[15, 14, 16, 18, 13],
[12, 11, 13, 15, 10],
[11, 10, 12, 14, 9]],
[[16, 15, 17, 19, 14],
[13, 12, 14, 16, 11],
[12, 11, 13, 15, 10]],
[[17, 16, 18, 20, 15],
[14, 13, 15, 17, 12],
[13, 12, 14, 16, 11]],
[[18, 17, 19, 21, 16],
[15, 14, 16, 18, 13],
[14, 13, 15, 17, 12]]])
Only numpy operations have the identity member, normal python operations don't have it or numpy's other special properties.
Some linear algebra properties and methods of ndarray
This section is not an exhaustive coverage of all linear algebra members in ndarray, simply because linear algebra is a big subject. I hope I can elaborate on them more in future posts. Assuming a is an ndarray:
-
a.Tanda.transpose()is the transpose of the array -
np.linalg.inv(a)is the inverse ofa -
np.eye(n)is the identity matrix of shape(n, n) -
@operator is matrix multiplication (Python 3.5+) -
np.trace(a)is the trace ofa -
np.linalg.solve(a, y)solves a system of equations, i.e. forms the equationa @ x = yand attempts to find the solution vectorx. -
np.linalg.eig(a)finds the eigenvalues ofa.
I have not checked if these functions only work on 2D arrays. If that's true it wouldn't be surprising as many linear algebra operations aren't defined for N-D arrays.
Jim Hefferon, Linear Algebra is a good freely available tutorial on linear algebra if you don't know about it yet.
And we're done
There are more numpy posts incoming, so stay tuned for more in the coming weeks.
Much material was sourced from the numpy manual.
If you see anything incorrect here please let me know so I can fix it.
Image by Arek Socha from Pixabay
Top comments (0)