Introduction
This blog is meant to be playground for beginners to learn about numpy by trying real codes. I have kept texts content as little as possible and code examples as much as possible.
This is also meant to be quick future reference guide of numpy features you already learned. The output of each cell has details describing results.
Prerequisite
 Basic programming knowledge
 Some familiarity with python (loops, arrays etc.)
What Will We Cover
Basics
 Creating Array
 Understanding Structure of Numpy Array (dimension, shape and strides)
 Data Types and Casting
 Indexing Methods
 Array Operations
Advance
What is Numpy?
Numpy is fundamental computing library for python. It supports NDimensional array and provides simple and efficient array operations.
NumPy is library of algorithms written in the C language which stores data in a contiguous block of memory, independent of other builtin Python objects and can operate on this memory without any type checking or other python overheads. NumPy arrays also use much less memory than builtin Python sequences.
Why Use Numpy?
Python was not initially designed for numerical computation. As python is interpreted language it is inherently slower than compiled languages like C. So numpy fills that gap here, following are few advantages of using numpy
 It provides efficient multidimensional array operations, both memory and computation wise
 It provides fast mathematical operations on entire array without need to use loop
 It also provides scientific operations related to linear algebra, statistics, Fourier transform and more
 It provides tool for interoperability with c and c++
How to Play With Numpy?
I would recommend two ways to play around with Numpy
 kaggle or Google Colab: you can jump right into coding without needing any setup
 Jupyter Notebook: you need to install jupyter notebook and then install numpy library using pip (numpy may be already installed if you have anaconda or miniconda)
Following This Tutorial
You can try this tutorial on kaggle or google colab by forking or making copy of following notebooks respectively
Fork this kaggle notebook
or
Make copy of this colab notebook
Some Tips On Jupyter Notebook
 To see auto completion suggestions press
tab
 To see parameters of a function press
shift + tab
after typing function name and '('. eg. typenp.asarray(
then pressshift + tab
 To view doc string use '?' like
np.asarray?
then pressshift + enter
Enough of reading!!
Lets get our hand dirty
Creating Array
From python list
import numpy as np
print(np.array([1,2,3,4]))
print('\n', 'array of 16 bit integers')
print(np.array([1,2,3,4], dtype=np.int16))
print('\n', '2 dimensional array')
print(np.array([[1,2,3], [4,5,6]]))
[1 2 3 4]
array of 16 bit integers
[1 2 3 4]
2 dimensional array
[[1 2 3]
[4 5 6]]
Numpy's Methods
print('Numpy array from range')
print(np.arange(3,8))
print('\n', '2D 3X3 array of zeros')
print(np.zeros((3,3)))
print('\n', '2D 2X3 array of ones')
print(np.ones((2,3)))
print('\n', 'Triangular array with ones at and below diagonal')
print(np.tri(3, 4))
print('\n', 'Index matrix with ones at diagonal')
print(np.eye(3))
print('\n', '20 equally spaced values between 1 and 5')
print(np.linspace(1, 5, 20))
Numpy array from range
[3 4 5 6 7]
2D 3X3 array of zeros
[[0. 0. 0.]
[0. 0. 0.]
[0. 0. 0.]]
2D 2X3 array of ones
[[1. 1. 1.]
[1. 1. 1.]]
Triangular array with ones at and below diagonal
[[1. 0. 0. 0.]
[1. 1. 0. 0.]
[1. 1. 1. 0.]]
Index matrix with ones at diagonal
[[1. 0. 0.]
[0. 1. 0.]
[0. 0. 1.]]
20 equally spaced values between 1 and 5
[1. 1.21052632 1.42105263 1.63157895 1.84210526 2.05263158
2.26315789 2.47368421 2.68421053 2.89473684 3.10526316 3.31578947
3.52631579 3.73684211 3.94736842 4.15789474 4.36842105 4.57894737
4.78947368 5. ]
Using np.random
print('3X2 array of uniformly distributed number between 0 and 1')
print(np.random.rand(3,2))
print('\n', 'Normally distributed random numbers with mean=0 and std=1')
print(np.random.randn(3,3))
print('\n', 'Randomly choose integers from a range (>=5, <11)')
print(np.random.randint(5, 11, size=(2,2)))
print('\n', "Randomly selects a permutation from array")
print(np.random.permutation([2,3,4,5,6]))
print('\n', "This is equivalent to rolling dice 10 times and counting \
occurance of getting each side")
print(np.random.multinomial(10, [1/6]*6))
3X2 array of uniformly distributed number between 0 and 1
[[0.32564739 0.97679242]
[0.26588925 0.89020385]
[0.18024366 0.90879681]]
Normally distributed random numbers with mean=0 and std=1
[[1.18661884 0.43561077 1.21316858]
[1.47847545 0.69296328 1.01348937]
[0.03562709 1.90675623 0.44639003]]
Randomly choose integers from a range (>=5, <11)
[[8 7]
[8 9]]
Randomly selects a permutation from array
[6 4 3 2 5]
This is equivalent to rolling dice 10 times and counting occurance of getting each side
[3 2 2 2 0 1]
Understanding Structure of Numpy Array (dimension, shape and strides)
import numpy as np
arr = np.array([[1,2,3], [2,3,1], [3,3,3]])
print('Number of array dimensions')
print(arr.ndim)
print('\nShape of array is tuple giving size of each dimension')
print(arr.shape)
print('\nstrides gives byte steps to be moved in memory to get to next \
index in each dimension')
print(arr.strides)
print('\nByte size of each item')
print(arr.itemsize)
Number of array dimensions
2
Shape of array is tuple giving size of each dimension
(3, 3)
strides gives byte steps to be moved in memory to get to next index in each dimension
(24, 8)
Byte size of each item
8
More on Strides
print('Slice indexing is done by changing strides, as in examples below')
print('Strides of original array')
print(arr.strides)
print('\n', 'Slice with step of 2 is done by multiplying stride(byte step size) by 2 in that dimension')
print(arr[::2].strides)
print('\n', 'Reverse index will negate the stride')
print(arr[::1].strides)
print('\n', 'Transpose will swap the stride of the dimensions')
print(arr.T.strides)
Slice indexing is done by changing strides, as in examples below
Strides of original array
(24, 8)
Slice with step of 2 is done by multiplying stride(byte step size) by 2 in that dimension
(48, 8)
Reverse index will negate the stride
(24, 8)
Transpose will swap the stride of the dimensions
(8, 24)
Some Stride Tricks: Inner product by changing strides
It is very rare that you may want to use these tricks but it helps us understand how indexing in numpy works
as_strided
function returns a view to an array with different strides and shape
from numpy.lib.stride_tricks import as_strided
arr1 = np.arange(5)
print('arr1: ', arr1)
arr2 = np.arange(3)
print('arr2: ', arr2)
print('\n', 'Adding a dimension with stride 0 allows us to repeat array in that dimension without making copy')
print('\n', 'Making stride 0 for rows repeats rows.')
print('As step size is zero to move to next row it will give same row repeatedly')
r_arr1 = as_strided(arr1, strides=(0,arr1.itemsize), shape=(len(arr2),len(arr1)))
print(r_arr1)
print('\n', 'Making stride 0 for columns repeats columns.')
r_arr2 = as_strided(arr1, strides=(arr2.itemsize, 0), shape=(len(arr2),len(arr1)))
print(r_arr2, '\n')
print('Inner product: product of every value of arr1 to every value of arr2')
print(r_arr1 * r_arr2)
arr1: [0 1 2 3 4]
arr2: [0 1 2]
Adding a dimension with stride 0 allows us to repeat array in that dimension without making copy
Making stride 0 for rows repeats rows.
As step size is zero to move to next row it will give same row repeatedly
[[0 1 2 3 4]
[0 1 2 3 4]
[0 1 2 3 4]]
Making stride 0 for columns repeats columns.
[[0 0 0 0 0]
[1 1 1 1 1]
[2 2 2 2 2]]
Inner product: product of every value of arr1 to every value of arr2
[[0 0 0 0 0]
[0 1 2 3 4]
[0 2 4 6 8]]
Using Broadcast
print('Above example is equivalent to using broadcast to do inner product')
print(arr1[np.newaxis, :] * arr2[:, np.newaxis])
print('arr1[np.newaxis, :].strides => ', arr1[np.newaxis, :].strides)
print('arr2[:, np.newaxis].strides => ', arr2[:, np.newaxis].strides)
Above example is equivalent to using broadcast to do inner product
[[0 0 0 0 0]
[0 1 2 3 4]
[0 2 4 6 8]]
arr1[np.newaxis, :].strides => (0, 8)
arr2[:, np.newaxis].strides => (8, 0)
Data Types and Casting
Notes
 Numpy array can store items of only one data type

np_array.dtype
attribute will give dtype of the array  following table shows some common datatypes with their string names
Numpy Attribute  String Name  Description

np.int8, np.int16, np.int32, np.int64  '<i1', '<i2', '<i4', '<i8'  signed int
np.uint8, np.uint16, np.uint32, np.uint64  '<u1', '<u2', '<u4', '<u8'  unsigned int
np.float16, np.float32, np.float64, np.float128  '<f2', '<f4', '<f8', '<f16' floats
np.string_ 'S1', 'S10', 'S255'  string of bytes (ascii)
np.str 'U1', 'U255'  string of unicode characters
np.datetime64 'M8'  date time
np.Object 'O'  python object
np.bool '?'  boolean
 Break down of string name '<u8': here '<' means littleendian byte order, 'u' means unsigned int and '8' means 8 bytes. Other options for byte order are '>' big endian and '=' system default
 All of the array initialization functions discussed above takes 'dtype' parameter to set datatype of the array eg:
np.random.randint(5, 11, size=(2,2), dtype=np.int8)
Casting
import numpy as np
arr = np.arange(5, dtype='<f4')
print('arr: ', arr)
print('\n', 'Cast to integer using astype function which will make copy of the array')
display(arr.astype(np.int8))
print('\n', 'By default casting is unsafe which will ignore the overflow. e.g. `2e10` is converted to 0')
arr[3] = 2e10
print(arr.astype('<i1'))
print('\n', 'Casting from string to float')
sarr = np.array("1 2 3 4 5.0".split())
print(sarr)
print(sarr.astype('<f4'))
print('\n', 'Use casting="safe" for doing safe casting, which will raise error if overflow')
# print(arr.astype('<i1', casting='safe'))
arr: [0. 1. 2. 3. 4.]
Cast to integer using astype function which will make copy of the array
array([0, 1, 2, 3, 4], dtype=int8)
By default casting is unsafe which will ignore the overflow. e.g. 2e10 is converted to 0
[0 1 2 0 4]
Casting from string to float
['1' '2' '3' '4' '5.0']
[1. 2. 3. 4. 5.]
Use casting="safe" for doing safe casting, which will raise error if overflow
Reshaping
 ndarray can be reshaped to any shape as long as total number of element are same
arr = np.arange(20)
print('arr: ', arr)
print('\n', 'reshape 1D arr of length 20 to shape (4,5)')
print(arr.reshape(4,5))
print('\n', 'One item of shape tuple can be 1 in which case the item will be calculated by numpy')
print('For total size to be 20 missing value must be 5')
print(arr.reshape(2,2,1))
arr: [ 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19]
reshape 1D arr of length 20 to shape (4,5)
[[ 0 1 2 3 4]
[ 5 6 7 8 9]
[10 11 12 13 14]
[15 16 17 18 19]]
One item of shape tuple can be 1 in which case the item will be calculated by numpy
For total size to be 20 missing value must be 5
[[[ 0 1 2 3 4]
[ 5 6 7 8 9]]
[[10 11 12 13 14]
[15 16 17 18 19]]]
Array View With Different dtype

arr.view()
method gives new view for same data with new dtype. Creating view with different dtype is not same as casting. eg. if we have ndarray of np.float32 ('<f4') creating view with dtype byte ('<i8') will read 4 bytes float data as individual bytes
arr = np.arange(5, dtype='<i2')
print('arr: ', arr)
print('\n', 'View with dtype "<i1" for array of dtype "<i2" will breakdown items to bytes')
print(arr.view('<i1'))
print('\n', 'Changing littleendian to bigendian will change value as they use different byte order')
print(arr.view('>i2'))
print('\n', 'Following will give individual bytes in memory of each items')
arr = np.arange(5, dtype='<f2')
print(arr)
print(arr.view('<i1'))
arr: [0 1 2 3 4]
View with dtype "<i1" for array of dtype "<i2" will breakdown items to bytes
[0 0 1 0 2 0 3 0 4 0]
Changing littleendian to bigendian will change value as they use different byte order
[ 0 256 512 768 1024]
Following will give individual bytes in memory of each items
[0. 1. 2. 3. 4.]
[ 0 0 0 60 0 64 0 66 0 68]
Indexing Methods
Integer and Slice Indexing
 This method of indexing is similar to indexing used in python list
 Slicing always create view to the array ie. does not copy the array
import numpy as np
arr = np.arange(20)
print("arr: ", arr)
print('\n', 'Get item at index 4(5th item) of the array')
print(arr[4])
print('\n', 'Assign 0 to index 4 of array')
arr[4] = 0
print(arr)
print('\n', 'Get items in the index range 4 to 10 not including 10')
print(arr[4:10])
print('\n', 'Set 1 to alternate items starting at index 4 to 10 ')
arr[4:10:2] = 1
print(arr)
arr: [ 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19]
Get item at index 4(5th item) of the array
4
Assign 0 to index 4 of array
[ 0 1 2 3 0 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19]
Get items in the index range 4 to 10 not including 10
[0 5 6 7 8 9]
Set 1 to alternate items starting at index 4 to 10
[ 0 1 2 3 1 5 1 7 1 9 10 11 12 13 14 15 16 17 18 19]
Slice Indexing in 2D Array
 For multidimensional array slice index can be separated using comma
arr = np.arange(20).reshape(4,5)
print('arr:\n', arr)
print('\n', 'Set 0 to first 3 rows and and last 2 columns')
arr[:3, 2:] = 1
print(arr)
arr:
[[ 0 1 2 3 4]
[ 5 6 7 8 9]
[10 11 12 13 14]
[15 16 17 18 19]]
Set 0 to first 3 rows and and last 2 columns
[[ 0 1 2 1 1]
[ 5 6 7 1 1]
[10 11 12 1 1]
[15 16 17 18 19]]
Boolean Indexing
 Boolean array of same shape as original array (or broadcastable to the shape) can be used as index. Which will select the items where index value is true
 Boolean array can also be used to filter array with certain conditions
 Boolean indexing will return copy instead of view to the array
arr = np.arange(6).reshape(2,3)
print('arr:\n', arr)
print('\n', 'Following index will gives last two items of 1st row and 1st element of 2nd row')
indx = np.array([[False, True, True], [True, False,False]])
arr[indx]
print('\n', 'Boolean index to filter values greater than 3 from arr')
filter_indx = arr>3
print('Filter Index:\n', filter_indx)
print('\n', 'Set 3 to values greater than 3 in arr')
arr[filter_indx] = 3
print(arr)
arr:
[[0 1 2]
[3 4 5]]
Following index will gives last two items of 1st row and 1st element of 2nd row
Boolean index to filter values greater than 3 from arr
Filter Index:
[[False False False]
[False True True]]
Set 3 to values greater than 3 in arr
[[0 1 2]
[3 3 3]]
Fancy Indexing
 Fancy Indexing means using array of index(integer) as index to get all items at once
 Fancy Indexing will also return copy instead of view to the array
import numpy as np
arr = np.arange(10)
print('arr:\n', arr)
print('\n', 'Get items at indexes 3,5 and 7 at once')
print(arr[[3,5,7]])
print('\n', 'Sorting arr based on another array "values"')
np.random.seed(5)
values = np.random.rand(10)
print('values:\n', values)
print('\n', 'np.argsort instead of returning sorted values will return array of indexes which will sort the array')
indexes = np.argsort(values)
print('indexes:\n', indexes)
print('Sorted array:\n', arr[indxes])
print('\n', 'You can also use fancy indexing to get same item multiple times')
print(arr[[0,1,1,2,2,2,3,3,3,3]])
arr:
[0 1 2 3 4 5 6 7 8 9]
Get items at indexes 3,5 and 7 at once
[3 5 7]
Sorting arr based on another array "values"
values:
[0.22199317 0.87073231 0.20671916 0.91861091 0.48841119 0.61174386
0.76590786 0.51841799 0.2968005 0.18772123]
np.argsort instead of returning sorted values will return array of indexes which will sort the array
indexes:
[9 2 0 8 4 7 5 6 1 3]

NameError Traceback (most recent call last)
<ipythoninput1430ac3b9b19e2> in <module>()
14 indexes = np.argsort(values)
15 print('indexes:\n', indexes)
> 16 print('Sorted array:\n', arr[indxes])
17
18 print('\n', 'You can also use fancy indexing to get same item multiple times')
NameError: name 'indxes' is not defined
Tuple Indexing
 Multidimensional array can be indexed using tuple of integer array of equal length, where each array in tuple will index corresponding dimension
 If number of indexarray in tuple is less than the dimension of array being indexed they will be used to index lower dimension (i.e. dimension starting from 0 to length of tuple)
arr2 = np.arange(15).reshape(5,3)
print('arr2:\n', arr2)
print('\n', 'Will give items at index (4,0) and (1,2)')
indx = ([4,1],[0,2])
print(arr2[indx])
print('\n', 'Tuple of length one will return rows')
indx = ([4,1],)
print(arr2[indx])
arr2:
[[ 0 1 2]
[ 3 4 5]
[ 6 7 8]
[ 9 10 11]
[12 13 14]]
Will give items at index (4,0) and (1,2)
[12 5]
Tuple of length one will return rows
[[12 13 14]
[ 3 4 5]]
Assignment With Advance Indexing
 Advance Indexing(i.e. Boolean, Fancy and Tuple) returns copy instead of view to the indexed array. But direct assignment using those index will change the original array, this feature is for convenience. But if we chain the indexing it may behaves in a way which seems to be somewhat unexpected
import numpy as np
arr = np.arange(10)
print('arr: ', arr)
print('\n', 'Direct assignment will change the original array')
arr[[3,5,7]] = 1
print(arr)
print('\n', 'When we chain the indexing it will not work')
arr[[3,5,7]][0] = 2
print(arr)
print('\n', 'But chaining index will work with slicing indexing')
arr[3:8:2][0] = 2
print(arr)
arr: [0 1 2 3 4 5 6 7 8 9]
Direct assignment will change the original array
[ 0 1 2 1 4 1 6 1 8 9]
When we chain the indexing it will not work
[ 0 1 2 1 4 1 6 1 8 9]
But chaining index will work with slicing indexing
[ 0 1 2 2 4 1 6 1 8 9]
Mixed Indexing
 In multi dimensional array we can use different indexing method (slicing, boolean and fancy) for each dimension at same time
 For mixture of boolean and fancy index to work, number of True's in boolean index must be equal to length of fancy index
arr = np.arange(64).reshape(4,4,4)
print('arr: ', arr)
print('\n', 'Following mixed indexing will select 1st and 3rd item in 0th dimension')
print('and item at index 0 and 2 at 1st dimension and item at index >=2')
print(arr[[True, False, True, False], [0,2], 2:])
arr: [[[ 0 1 2 3]
[ 4 5 6 7]
[ 8 9 10 11]
[12 13 14 15]]
[[16 17 18 19]
[20 21 22 23]
[24 25 26 27]
[28 29 30 31]]
[[32 33 34 35]
[36 37 38 39]
[40 41 42 43]
[44 45 46 47]]
[[48 49 50 51]
[52 53 54 55]
[56 57 58 59]
[60 61 62 63]]]
Following mixed indexing will select 1st and 3rd item in 0th dimension
and item at index 0 and 2 at 1st dimension and item at index >=2
[[ 2 3]
[42 43]]
Array Operations
Simple Array Operations
 Numpy provides simple syntax to perform mathematical and logical operations between array of compatible shapes. Here compatible shape means, shape of one of the array can be expanded to match shape of other using broadcast rule, which we'll discuss below
 For this section we'll only see two cases
 Arrays has same shape in which case operation will be element wise
 One of operand is scalar in which case operation will be done between scalar and each element of array
 These operations between arrays are called vectorization and are way faster than same operation using loop.
 Vectorization are faster because it is implemented in C and don't have overheads like type checking etc.
import numpy as np
print('Evaluate expression (x1*x2  3*x1 + 30) for x1 and x2 in range 0 to 10')
x1 = np.linspace(0,10,20)
x2 = np.linspace(0, 10, 20)
z = x1*x2  3*x1 + 30
print(z)
print('\n', 'Spatial distance between corresponding points in two array')
p1 = np.random.rand(20,2)
p2 = np.random.rand(20,2)
'''np.sum will add values along given axis (dimension). If shape of array is (3,4,5)')
then axis 0,1 and 2 corresponds to dimension with length 3, 4 and 5 respectively'''
d = np.sum((p1p2)**2, axis=1)
print(d)
print('\n', 'Element wise comparison, ">=" will give boolean array with True where element')
print('of p2 is greater than or equal to p1')
r = p2>=p1
print(r)
print('\n', 'Element wise logical operation, "&" will give True where point of p2 is ahead')
print('in both x and y direction from corresponding point in p1')
print(r[:,0] & r[:,1])
Evaluate expression (x1*x2  3*x1 + 30) for x1 and x2 in range 0 to 10
[ 30. 28.69806094 27.9501385 27.75623269 28.11634349
29.03047091 30.49861496 32.52077562 35.09695291 38.22714681
41.91135734 46.14958449 50.94182825 56.28808864 62.18836565
68.64265928 75.65096953 83.2132964 91.32963989 100. ]
Spatial distance between corresponding points in two array
[0.54052263 0.17505988 0.59108818 0.41593393 0.03548522 0.29946201
0.84649163 0.24975051 0.90016153 0.54062043 0.00097261 0.39826495
0.64710327 0.40655563 0.00531519 0.94567232 0.33333277 0.01713418
0.53797027 0.48080742]
Element wise comparison, ">=" will give boolean array with True where element
of p2 is greater than or equal to p1
[[ True False]
[False False]
[False False]
[ True True]
[ True False]
[ True True]
[ True True]
[ True False]
[ True False]
[False True]
[ True False]
[False False]
[ True False]
[ True False]
[False False]
[ True False]
[False False]
[ True False]
[False True]
[False False]]
Element wise logical operation, "&" will give True where point of p2 is ahead
in both x and y direction from corresponding point in p1
[False False False True False True True False False False False False
False False False False False False False False]
Functions for Array Operations
 Numpy also has function version of above operations like
np.add, np.substract, np.divide, np.greater_equal, np.logical_and
and more  Array operations we see in section above using operators like +, * are operator overloaded version of function operation
 Function version of the operation will give us extra parameters to customize, one of commonly used parameter is
out
. It isNone
by default, which will create a new array for result.  If We pass an array with shape and dtype matching to expected result to
out
parameter result will be filled to the passed array. It will be efficient memory wise if we are doing multiple operations
import numpy as np
print('Evaluate expression (x1*x2  3*x1 + 30) using functions')
x1 = np.linspace(0,10,20)
x2 = np.linspace(0, 10, 20)
'''Create empty output array with expected shape'''
z = np.empty_like(x1)
'''Code is not very clean as using operator but it will perform very well memory wise'''
np.multiply(x1, x2, out=z)
np.subtract(z, 3*x1, out=z)
np.add(z, 30, out=z)
print(z)
Evaluate expression (x1*x2  3*x1 + 30) using functions
[ 30. 28.69806094 27.9501385 27.75623269 28.11634349
29.03047091 30.49861496 32.52077562 35.09695291 38.22714681
41.91135734 46.14958449 50.94182825 56.28808864 62.18836565
68.64265928 75.65096953 83.2132964 91.32963989 100. ]
Broadcasting
Rules for broadcasting
When doing array operations between two array's whose shape doesn't match exactly then few simple steps are taken which changes shapes to match each other if they are compatible.
 Check Array Dimensions: If dimensions doesn't match added to left of array with smaller dimension
 Match Shape On Each Dimension: If shape in any dimension doesn't match and shape is 1 for one of the array then repeat it to match shape of other array in that dimension
 Raise Error if Dimension and Shape Not Matched: If dimension and shape don't match till this step then raise error
Let's Visualize the Broadcasting Rule With Custom Implementation
lets do our custom implementation to visualize in code how the broadcast rule works
import numpy as np
arr1 = np.arange(20).reshape(10,2)
arr2 = np.random.rand(2)
arr3 = arr2.copy()
print('arr1.shape: ', arr1.shape)
print('arr2.shape: ', arr2.shape)
# Step 1: Check Array Dimensions
print('\n', 'arr1 has dimension 2 and arr2 has dimension 1, so add 1 dimension to\
left side of arr2')
# np.newaxis is convenient way of adding new dimension
arr2 = arr2[np.newaxis, :]
print('arr1.shape: ', arr1.shape)
print('arr2.shape: ', arr2.shape)
# Step 2: Match Shape On Each Dimension
print('\n', 'Now in axis=0 arr1 has 10 items and arr2 has one item, so repeat it 10\
times to match arr2')
arr2 = np.repeat(arr2, 10, axis=0)
print('arr1.shape: ', arr1.shape)
print('arr2.shape: ', arr2.shape)
print('\n', 'Now both array has same dimension and shape, we can multiply them')
print('arr1*arr2:\n', arr1*arr2)
print('\n', 'Lets see if broadcasting also produce same result')
print('arr1*arr3:\n', arr1*arr3)
arr1.shape: (10, 2)
arr2.shape: (2,)
arr1 has dimension 2 and arr2 has dimension 1, so add 1 dimension toleft side of arr2
arr1.shape: (10, 2)
arr2.shape: (1, 2)
Now in axis=0 arr1 has 10 items and arr2 has one item, so repeat it 10times to match arr2
arr1.shape: (10, 2)
arr2.shape: (10, 2)
Now both array has same dimension and shape, we can multiply them
arr1*arr2:
[[ 0. 0.11111075]
[ 1.71941377 0.33333225]
[ 3.43882755 0.55555375]
[ 5.15824132 0.77777525]
[ 6.8776551 0.99999675]
[ 8.59706887 1.22221824]
[10.31648264 1.44443974]
[12.03589642 1.66666124]
[13.75531019 1.88888274]
[15.47472397 2.11110424]]
Lets see if broadcasting also produce same result
arr1*arr3:
[[ 0. 0.11111075]
[ 1.71941377 0.33333225]
[ 3.43882755 0.55555375]
[ 5.15824132 0.77777525]
[ 6.8776551 0.99999675]
[ 8.59706887 1.22221824]
[10.31648264 1.44443974]
[12.03589642 1.66666124]
[13.75531019 1.88888274]
[15.47472397 2.11110424]]
Let's Try Few Examples With Shapes
you can try it by creating arrays of given shape and doing some operation between them
Before Broadcast Step 1  Step 2 and 3
Shapes of arr1 and arr2   Shapes of result

(3, 1, 5); (4, 1)  (3, 1, 5); (1, 4, 1)  (3, 4, 5)
(10,); (1, 10)  (10, 1); (1, 10)  (10, 10)
(2, 2, 2); (2, 3)  (2, 2, 2); (1, 2, 3)  Not Broadcastable
(2, 2, 2, 1); (2, 3)  (2, 2, 2, 1); (1, 1, 2, 3)  (2, 2, 2, 2, 3)
Some Usage of Broadcast
Evaluate Linear Equation Using Broadcast
print("Let's evaluate equation c1*x1 + c2*x2 + c3*x3 for 100 points at once")
points = np.random.rand(100,3)
coefficients = np.array([5, 2, 11])
results = np.sum(points*coefficients, axis=1)
print('results first 10:\n', results[:10])
print('results.shape: ', results.shape)
Let's evaluate equation c1*x1 + c2*x2 + c3*x3 for 100 points at once
results first 10:
[ 6.35385279 0.85639146 12.87683079 5.99433896 4.50873972 10.44691041
3.87407211 6.62954602 11.00386582 10.09247866]
results.shape: (100,)
Find Common Elements Between Arrays
np.random.seed(5)
## Get 20 random value from 0 to 99
arr1 = np.random.choice(50, 20, replace=False)
arr2 = np.random.choice(50, 15, replace=False)
print("arr1: ", arr1)
print("arr2: ", arr2)
print('\n', 'arr1 and arr2 are 1D arrays of length 20, 15 respectively.')
print('To make them broadcastable Change shape of arr2 to (15, 1)')
arr2 = arr2.reshape(15, 1)
print('\n', 'Then both arrays will be broadcasted to (15, 20) matrix with all possible pairs')
comparison = (arr1 == arr2)
print('\n', 'comparison.shape: ', comparison.shape)
print('\n', 'Elements of arr1 also in arr2: ', arr1[comparison.any(axis=0)])
arr1: [42 29 6 19 28 17 2 43 3 21 31 4 32 0 23 5 48 34 37 26]
arr2: [40 37 41 48 4 20 10 18 34 28 19 32 17 22 23]
arr1 and arr2 are 1D arrays of length 20, 15 respectively.
To make them broadcastable Change shape of arr2 to (15, 1)
Then both arrays will be broadcasted to (15, 20) matrix with all possible pairs
comparison.shape: (15, 20)
Elements of arr1 also in arr2: [19 28 17 4 32 23 48 34 37]
Find knearest Neighbors
import numpy as np
np.random.seed(5)
points = np.random.rand(20, 2)
print('To calculate distance between every pair of points make copy of points ')
print('with shape (20, 1, 2) which will broadcast both array to shape (20, 20, 2)', '\n')
cp_points = points.reshape(20, 1, 2)
## calculate x2x1, y2y1
diff = (cp_points  points)
print('diff.shape: ', diff.shape)
## calculate (x2x1)**2 + (y2y1)**
distance_matrix = np.sum(diff**2, axis=2)
print('distance_matrix.shape: ', distance_matrix.shape, '\n')
## sort by distance along axis 1 and take top 4, one of which is the point itself
top_3 = np.argsort(distance_matrix, axis=1)[:,:4]
print("Get the points with it's 3 nearest neighbors")
points[top_3]
To calculate distance between every pair of points make copy of points
with shape (20, 1, 2) which will broadcast both array to shape (20, 20, 2)
diff.shape: (20, 20, 2)
distance_matrix.shape: (20, 20)
Get the points with it's 3 nearest neighbors
array([[[0.22199317, 0.87073231],
[0.20671916, 0.91861091],
[0.16561286, 0.96393053],
[0.08074127, 0.7384403 ]],
[[0.20671916, 0.91861091],
[0.22199317, 0.87073231],
[0.16561286, 0.96393053],
[0.08074127, 0.7384403 ]],
[[0.48841119, 0.61174386],
[0.62878791, 0.57983781],
[0.69984361, 0.77951459],
[0.76590786, 0.51841799]],
[[0.76590786, 0.51841799],
[0.62878791, 0.57983781],
[0.69984361, 0.77951459],
[0.87993703, 0.27408646]],
[[0.2968005 , 0.18772123],
[0.32756395, 0.1441643 ],
[0.28468588, 0.25358821],
[0.44130922, 0.15830987]],
[[0.08074127, 0.7384403 ],
[0.02293309, 0.57766286],
[0.22199317, 0.87073231],
[0.20671916, 0.91861091]],
[[0.44130922, 0.15830987],
[0.32756395, 0.1441643 ],
[0.41423502, 0.29607993],
[0.2968005 , 0.18772123]],
[[0.87993703, 0.27408646],
[0.96022672, 0.18841466],
[0.76590786, 0.51841799],
[0.5999292 , 0.26581912]],
[[0.41423502, 0.29607993],
[0.28468588, 0.25358821],
[0.44130922, 0.15830987],
[0.2968005 , 0.18772123]],
[[0.62878791, 0.57983781],
[0.48841119, 0.61174386],
[0.76590786, 0.51841799],
[0.69984361, 0.77951459]],
[[0.5999292 , 0.26581912],
[0.41423502, 0.29607993],
[0.44130922, 0.15830987],
[0.87993703, 0.27408646]],
[[0.28468588, 0.25358821],
[0.2968005 , 0.18772123],
[0.32756395, 0.1441643 ],
[0.41423502, 0.29607993]],
[[0.32756395, 0.1441643 ],
[0.2968005 , 0.18772123],
[0.44130922, 0.15830987],
[0.28468588, 0.25358821]],
[[0.16561286, 0.96393053],
[0.20671916, 0.91861091],
[0.22199317, 0.87073231],
[0.08074127, 0.7384403 ]],
[[0.96022672, 0.18841466],
[0.87993703, 0.27408646],
[0.5999292 , 0.26581912],
[0.76590786, 0.51841799]],
[[0.02430656, 0.20455555],
[0.28468588, 0.25358821],
[0.2968005 , 0.18772123],
[0.32756395, 0.1441643 ]],
[[0.69984361, 0.77951459],
[0.62878791, 0.57983781],
[0.63979518, 0.9856244 ],
[0.76590786, 0.51841799]],
[[0.02293309, 0.57766286],
[0.00164217, 0.51547261],
[0.08074127, 0.7384403 ],
[0.22199317, 0.87073231]],
[[0.00164217, 0.51547261],
[0.02293309, 0.57766286],
[0.08074127, 0.7384403 ],
[0.02430656, 0.20455555]],
[[0.63979518, 0.9856244 ],
[0.69984361, 0.77951459],
[0.48841119, 0.61174386],
[0.62878791, 0.57983781]]])
Vectorization
 In numpy vectorization means performing optimized operations on sequence of same type of data.
 In addition to having clean structure of code vectorized operations are also very performant because codes are compiled an avoids overhead of python, like type checking, memory management etc.
 The examples we see on Broadcast section above are also good example of vectorization
Vectorization Vs Loop
Let's say we have a polynomial equation of degree 10 of single variable like a1*x + a2*x^2 + a3*x^3 ... + a10*x^10
. Let's try evaluating the equation for large number of x using python only and numpy vectorization see how they compare
import numpy as np
np.random.seed(32)
X = np.random.rand(10000)
coefficients = np.random.randn(10)*20 + 50
def evaluate_polynomial_loop():
result_loop = np.empty_like(X)
for i in range(X.shape[0]):
exp_part = 1
total = 0
for j in range(coefficients.shape[0]):
exp_part *= X[i]
total+=coefficients[j]*exp_part
result_loop[i] = total
return result_loop
def evaluate_polynomial_vect():
## repeates x's in 10 columns
exponents = X[:, np.newaxis] + np.zeros((1, coefficients.shape[0]))
exponents.cumprod(axis=1, out=exponents)
result_vect = np.sum(exponents * coefficients, axis=1)
return result_vect
print('Verify that both gives same result')
print('Loop:\n', evaluate_polynomial_loop()[:10])
print('Vectorization:\n', evaluate_polynomial_vect()[:10])
Verify that both gives same result
Loop:
[222.57782534 30.62439847 59.69953776 373.52687079 123.89007218
179.70369976 6.49315699 321.685257 73.14575517 69.71437596]
Vectorization:
[222.57782534 30.62439847 59.69953776 373.52687079 123.89007218
179.70369976 6.49315699 321.685257 73.14575517 69.71437596]
Compare
For fair comparison I used numpy array in both whose indexing is much faster than python list. By comparison we see the vectorization is about 80 times faster
%timeit evaluate_polynomial_loop()
114 ms Β± 5.02 ms per loop (mean Β± std. dev. of 7 runs, 10 loops each)
%timeit evaluate_polynomial_vect()
1.19 ms Β± 79 Β΅s per loop (mean Β± std. dev. of 7 runs, 1000 loops each)
Ufunc and Numba
Ufunc: Also called Universal Functions are vectorized wrapper for a function. Ufunc can operate on ndarray and support broadcasting and type casting. np.add, np.multiply
etc are examples of Ufunc which are implemented in C. We can create custom Ufunc using np.frompyfunc
or using numba.
Numba: Numba is just in time compiler which generate optimized machine code from pure python array and numerical functions. You can use numba.jit
decorator on a function which makes the function to be compiled on its 1st run. You can use numba.vectorize
decorator to convert python function to Ufunc.
Let's compare different implementations of adding two big arrays as follow
Create Big Arrays
arr1 = np.arange(1000000, dtype='int64')
arr2 = np.arange(1000000, dtype='int64')
Implementation Using Python Loop
def add_arr(arr1, arr2):
assert len(arr1)==len(arr2), "array must have same length"
result = np.empty_like(arr1)
for i in range(len(arr1)):
result[i] = arr1[i] + arr2[i]
return result
%timeit _ = add_arr(arr1, arr2)
596 ms Β± 24.8 ms per loop (mean Β± std. dev. of 7 runs, 1 loop each)
Creating Ufunc Using np.frompyfunc
def add(a, b):
return a+b
vect_add = np.frompyfunc(add,2,1)
%timeit _ = vect_add(arr1, arr2)
195 ms Β± 12.9 ms per loop (mean Β± std. dev. of 7 runs, 1 loop each)
Using Numba JIT
 'nopython=True' means convert all to machine level code, if can't convert raise error
import numba
@numba.jit(nopython=True)
def add_arr_jit(arr1, arr2):
assert len(arr1)==len(arr2), "array must have same length"
result = np.empty_like(arr1)
for i in range(len(arr1)):
result[i] = arr1[i] + arr2[i]
return result
_ = add_arr_jit(arr1, arr2)
%timeit _ = add_arr_jit(arr1, arr2)
3.66 ms Β± 467 Β΅s per loop (mean Β± std. dev. of 7 runs, 100 loops each)
Creating Ufunc Using numba.vectorize
 'numba.vectorize' takes signature of function being converted as parameter. 'int64(int64,int64)' means takes 2 'int64' parameters and returns 'int64'
import numba
@numba.vectorize(['int64(int64,int64)'], nopython=True)
def vect_add(a, b):
return a+b
%timeit _ = vect_add(arr1, arr2)
2.93 ms Β± 569 Β΅s per loop (mean Β± std. dev. of 7 runs, 100 loops each)
Conclusion solution using numba.jit and numba.vectorize
are performing much better. You can also check how numpy vectorization compares with these
More for Exploration
Some Useful Functions

np.where
: element wiseif .. then .. Else

np.select
: select values from multiple arrays based on multiple conditions 
np.concatenate, np.vstack, np.r_, np.hstack, np.c_
: join multiple ndarray row wise, column wise or in a given axis 
np.ravel, np.flatten
: converts multidimensional array to 1D array 
np.roll
: do circular shift of the array along given axis
Set Operations
 *
np.unique(x)
: * gives array of unique elements on the array  *
Intersect1d(x, y)
: * gives 1D array of elements common to both arrays  *
Union1d(x, y)
: * gives 1D array of unique elements from both arrays  *
In1d(x, y)
: * check if each element of x is also present on y and returns array of length equal to x with boolean values  *
Setdiff1d(x, y)
: * gives elements of x not in y  *
Setxor1d(x, y)
: * give elements that is either in x or y but not in both
Save and Load ndarray From Disk

np.save("filename.npy", x)
: save single numpy array to disk 
np.load("filename.npy")
: load single numpy array from disk  *
np.savez("filename.npz", key1=arr1, key2=arr2)
: * saves multiple arrays with given key 
np.savetxt("filename.npy", x)
: save single numpy array to disk as delimited text file 
np.loadtxt("filename.npy")
: load single numpy array from text file
Memory Mapping
To work with large numpy array that doesn't fit in RAM you can use numpy.memmap function to map the array to file in disk. It will transparently loads only segments of array needed for current operations.

np.memmap(filename, dtype, mode, shape)
: create a memory mapped array to a given file 
mmap.flush()
: flush all in memory changes to disk
Thank you
Heartily thanks for reading through the blog. Hope it was helpful for getting you up and running in numpy. Any comment, suggestion and constructive criticism are welcome. If you like the content please donβt forget to like π
Top comments (0)