NumPy

NumPy is the core Python package for numerical computing. The main features of NumPy are:

$N$-dimensional array object ndarray
Vectorized operations and functions which broadcast across arrays for fast computation

To get started with NumPy, let's adopt the standard convention and import it using the name np:

import numpy as np

NumPy Arrays

The fundamental object provided by the NumPy package is the ndarray. We can think of a 1D (1-dimensional) ndarray as a list, a 2D (2-dimensional) ndarray as a matrix, a 3D (3-dimensional) ndarray as a 3-tensor (or a "cube" of numbers), and so on. See the NumPy tutorial for more about NumPy arrays.

Creating Arrays

The function numpy.array creates a NumPy array from a Python sequence such as a list, a tuple or a list of lists. For example, create a 1D NumPy array from a Python list:

a = np.array([1,2,3,4,5])
print(a)

[1 2 3 4 5]

Notice that when we print a NumPy array it looks a lot like a Python list except that the entries are separated by spaces whereas entries in a Python list are separated by commas:

print([1,2,3,4,5])

[1, 2, 3, 4, 5]

Notice also that a NumPy array is displayed slightly differently when output by a cell (as opposed to being explicitly printed to output by the print function):

array([1, 2, 3, 4, 5])

Use the built-in function type to verify the type:

type(a)

numpy.ndarray

Create a 2D NumPy array from a Python list of lists:

M = np.array([[1,2,3],[4,5,6]])
print(M)

[[1 2 3]
 [4 5 6]]

type(M)

numpy.ndarray

Create an $n$-dimensional NumPy array from nested Python lists. For example, the following is a 3D NumPy array:

N = np.array([ [[1,2],[3,4]] , [[5,6],[7,8]] , [[9,10],[11,12]] ])
print(N)

[[[ 1  2]
  [ 3  4]]

 [[ 5  6]
  [ 7  8]]

 [[ 9 10]
  [11 12]]]

There are several NumPy functions for creating arrays:

Function	Description
`numpy.array(a)`	Create $n$-dimensional NumPy array from sequence `a`
`numpy.linspace(a,b,N)`	Create 1D NumPy array with `N` equally spaced values from `a` to `b` (inclusively)
`numpy.arange(a,b,step)`	Create 1D NumPy array with values from `a` to `b` (exclusively) incremented by `step`
`numpy.zeros(N)`	Create 1D NumPy array of zeros of length $N$
`numpy.zeros((n,m))`	Create 2D NumPy array of zeros with $n$ rows and $m$ columns
`numpy.ones(N)`	Create 1D NumPy array of ones of length $N$
`numpy.ones((n,m))`	Create 2D NumPy array of ones with $n$ rows and $m$ columns
`numpy.eye(N)`	Create 2D NumPy array with $N$ rows and $N$ columns with ones on the diagonal (ie. the identity matrix of size $N$)

Create a 1D NumPy array with 11 equally spaced values from 0 to 1:

x = np.linspace(0,1,11)
print(x)

[0.  0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1. ]

Create a 1D NumPy array with values from 0 to 20 (exclusively) incremented by 2.5:

y = np.arange(0,20,2.5)
print(y)

[ 0.   2.5  5.   7.5 10.  12.5 15.  17.5]

These are the functions that we'll use most often when creating NumPy arrays. The function numpy.linspace works best when we know the number of points we want in the array, and numpy.arange works best when we know step size between values in the array.

Create a 1D NumPy array of zeros of length 5:

z = np.zeros(5)
print(z)

[0. 0. 0. 0. 0.]

Create a 2D NumPy array of zeros with 2 rows and 5 columns:

M = np.zeros((2,5))
print(M)

[[0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0.]]

Create a 1D NumPy array of ones of length 7:

w = np.ones(7)
print(w)

[1. 1. 1. 1. 1. 1. 1.]

Create a 2D NumPy array of ones with 3 rows and 2 columns:

N = np.ones((3,2))
print(N)

[[1. 1.]
 [1. 1.]
 [1. 1.]]

Create the identity matrix of size 10:

I = np.eye(10)
print(I)

[[1. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 1. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 1. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 1. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 1. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 1. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 1. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 1. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 1. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 1.]]

Array Datatype

NumPy arrays are homogeneous: all entries in the array are the same datatype. We will only work with numeric arrays and our arrays will contain either integers, floats, complex numbers or booleans. There are different kinds of datatypes provided by NumPy for different applications but we'll mostly be working with the default integer type numpy.int64 and the default float type numpy.float64. These are very similar to the built-in Python datatypes int and float but with some differences that we won't go into. Check out the NumPy documentation on numeric datatypes for more information.

The most important point for now is to know how to determine if a NumPy array contains integers elements or float elements. We can access the datatype of a NumPy array by its .dtype attribute. For example, create a 2D NumPy array from a list of lists of integers:

A = np.array([[1,2,3],[4,5,6]])
print(A)

[[1 2 3]
 [4 5 6]]

We expect the datatype of A to be integers and we verify:

A.dtype

dtype('int64')

Most of the other NumPy functions which create arrays use the numpy.float64 datatype by default. For example, using numpy.linspace:

u = np.linspace(0,1,5)
print(u)

[0.   0.25 0.5  0.75 1.  ]

u.dtype

dtype('float64')

Notice that numbers are printed with a decimal point when the datatype of the NumPy array is any kind of float.

Dimension, Shape and Size

We can think of a 1D NumPy array as a list of numbers, a 2D NumPy array as a matrix, a 3D NumPy array as a cube of numbers, and so on. Given a NumPy array, we can find out how many dimensions it has by accessing its .ndim attribute. The result is a number telling us how many dimensions it has.

For example, create a 2D NumPy array:

A = np.array([[1,2],[3,4],[5,6]])
print(A)

[[1 2]
 [3 4]
 [5 6]]

A.ndim

The result tells us that A has 2 dimensions. The first dimension corresponds to the vertical direction counting the rows and the second dimension corresponds to the horizontal direction counting the columns.

We can find out how many rows and columns A has by accessing its .shape attribute:

A.shape

(3, 2)

The result is a tuple (3,2) of length 2 which means that A is a 2D array with 3 rows and 2 columns.

We can also find out how many entries A has in total by accessing its .size attribute:

A.size

This is the expected result since we know that A has 3 rows and 2 columns and therefore 2(3) = 6 total entries.

Create a 1D NumPy array and inspect its dimension, shape and size:

r = np.array([9,3,1,7])
print(r)

[9 3 1 7]

r.ndim

r.shape

(4,)

r.size

The variable r is assigned to a 1D NumPy array of length 4. Notice that r.shape is a tuple with a single entry (4,).

Slicing and Indexing

Accessing the entries in an array is called indexing and accessing rows and columns (or subarrays) is called slicing. See the NumPy documentation for more information about indexing and slicing.

Create a 1D NumPy array:

v = np.linspace(0,5,11)
print(v)

[0.  0.5 1.  1.5 2.  2.5 3.  3.5 4.  4.5 5. ]

Access the entries in a 1D array using the square brackets notation just like a Python list. For example, access the entry at index 3:

v[3]

1.5

Notice that NumPy array indices start at 0 just like Python sequences.

Create a 2D array of integers:

B = np.array([[6, 5, 3, 1, 1],[1, 0, 4, 0, 1],[5, 9, 2, 2, 9]])
print(B)

[[6 5 3 1 1]
 [1 0 4 0 1]
 [5 9 2 2 9]]

Access the entries in a 2D array using the square brackets with 2 indices. In particular, access the entry at row index 1 and column index 2:

B[1,2]

Access the top left entry in the array:

B[0,0]

Negative indices work for NumPy arrays as they do for Python sequences. Access the bottom right entry in the array:

B[-1,-1]

Access the row at index 2 using the colon : syntax:

B[2,:]

array([5, 9, 2, 2, 9])

Access the column at index 3:

B[:,3]

array([1, 0, 2])

Select the subarray of rows at index 1 and 2, and columns at index 2, 3 and 4:

subB = B[1:3,2:5]
print(subB)

[[4 0 1]
 [2 2 9]]

Slices of NumPy arrays are again NumPy arrays but possibly of a different dimension:

subB.ndim

subB.shape

(2, 3)

type(subB)

numpy.ndarray

The variable subB is assigned to a 2D NumPy array of shape 2 by 2.

Let's do the same for the column at index 2:

colB = B[:,2]
print(colB)

[3 4 2]

colB.ndim

colB.shape

(3,)

type(colB)

numpy.ndarray

The variable colB is assigned to a 1D NumPy array of length 3.

Stacking

We can build bigger arrays out of smaller arrays by stacking along different dimensions using the functions numpy.hstack and numpy.vstack.

Stack 3 different 1D NumPy arrays of length 3 vertically forming a 3 by 3 matrix:

x = np.array([1,1,1])
y = np.array([2,2,2])
z = np.array([3,3,3])
vstacked = np.vstack((x,y,z))
print(vstacked)

[[1 1 1]
 [2 2 2]
 [3 3 3]]

Stack 1D NumPy arrays horizontally to create another 1D array:

hstacked = np.hstack((x,y,z))
print(hstacked)

[1 1 1 2 2 2 3 3 3]

Use numpy.hstack and numpy.vstack to build the matrix $T$ where

$$ T = \begin{bmatrix} 1 & 1 & 2 & 2 \ 1 & 1 & 2 & 2 \ 3 & 3 & 4 & 4 \ 3 & 3 & 4 & 4 \end{bmatrix} $$

A = np.ones((2,2))
B = 2*np.ones((2,2))
C = 3*np.ones((2,2))
D = 4*np.ones((2,2))
A_B = np.hstack((A,B))
print(A_B)

[[1. 1. 2. 2.]
 [1. 1. 2. 2.]]

C_D = np.hstack((C,D))
print(C_D)

[[3. 3. 4. 4.]
 [3. 3. 4. 4.]]

T = np.vstack((A_B,C_D))
print(T)

[[1. 1. 2. 2.]
 [1. 1. 2. 2.]
 [3. 3. 4. 4.]
 [3. 3. 4. 4.]]

Copies versus Views

Under construction

Operations and Functions

Array Operations

Arithmetic operators including addition +, subtraction -, multiplication *, division / and exponentiation ** are applied to arrays elementwise. For addition and substraction, these are the familiar vector operations we see in linear algebra:

v = np.array([1,2,3])
w = np.array([1,0,-1])

v + w

array([2, 2, 2])

v - w

array([0, 2, 4])

In the same way, array multiplication and division are performed element by element:

v * w

array([ 1,  0, -3])

w / v

array([ 1.        ,  0.        , -0.33333333])

Notice that the datatype of both v and w is numpy.int64 however division w / v returns an array with datatype numpy.float64.

The exponent operator ** also acts element by element in the array:

v ** 2

array([1, 4, 9])

Let's see these operations for 2D arrays:

A = np.array([[3,1],[2,-1]])
B = np.array([[2,-2],[5,1]])

A + B

array([[ 5, -1],
       [ 7,  0]])

A - B

array([[ 1,  3],
       [-3, -2]])

A / B

array([[ 1.5, -0.5],
       [ 0.4, -1. ]])

A * B

array([[ 6, -2],
       [10, -1]])

A ** 2

array([[9, 1],
       [4, 1]])

Notice that array multiplication and exponentiation are performed elementwise.

In Python 3.5+, the symbol @ computes matrix multiplication for NumPy arrays:

A @ B

array([[11, -5],
       [-1, -5]])

Matrix powers are performed by the function numpy.linalg.matrix_power. It's a long function name and so it's convenient to import it with a shorter name:

from numpy.linalg import matrix_power as mpow

Compute $A^3$:

mpow(A,3)

array([[37,  9],
       [18,  1]])

Equivalently, use the @ operator to compute $A^3$:

A @ A @ A

array([[37,  9],
       [18,  1]])

Broadcasting

We know from linear algebra that we can only add matrices of the same size. Braodcasting is a set of NumPy rules which relaxes this constraint and allows us to combine a smaller array with a bigger when it makes sense.

For example, suppose we want to create a 1D NumPy array of $y$ values for $x=0.0,0.25,0.5,0.75,1.0$ for the function $y = x^2 + 1$. From what we've seen so far, it makes sense to create x, then x**2 and then add an array of ones [1. 1. 1. 1. 1.]:

x = np.array([0,0.25,0.5,0.75,1.0])
y = x**2 + np.array([1,1,1,1,1])
print(y)

[1.     1.0625 1.25   1.5625 2.    ]

An example of broadcasting in NumPy is the following equivalent operation:

x = np.array([0,0.25,0.5,0.75,1.0])
y = x**2 + 1
print(y)

[1.     1.0625 1.25   1.5625 2.    ]

The number 1 is a scalar and we are adding it to a 1D NumPy array of length 5. The broadcasting rule in this case is to broadcast the scalar value 1 across the larger array. The result is a simpler syntax for a very comman operation.

Let's try another example. What happens when we try to add a 1D NumPy array of length 4 to a 2D NumPy array of size 3 by 4?

u = np.array([1,2,3,4])
print(u)

[1 2 3 4]

A = np.array([[1,1,1,1],[2,2,2,2],[3,3,3,3]])
print(A)

[[1 1 1 1]
 [2 2 2 2]
 [3 3 3 3]]

result = A + u
print(result)

[[2 3 4 5]
 [3 4 5 6]
 [4 5 6 7]]

The 1D NumPy array is broadcast across the 2D array because the length of the first dimension in each array are equal!

Array Functions

There are many array functions we can use to compute with NumPy arrays. The following is a partial list and we'll look closer at mathematical functions in the next section.


`numpy.sum`	`numpy.prod`	`numpy.mean`
`numpy.max`	`numpy.min`	`numpy.std`
`numpy.argmax`	`numpy.argmin`	`numpy.var`

Create a 1D NumPy array with random values and compute:

arr = np.array([8,-2,4,7,-3])
print(arr)

[ 8 -2  4  7 -3]

Compute the mean of the values in the array:

np.mean(arr)

2.8

Verify the mean once more:

m = np.sum(arr) / arr.size
print(m)

2.8

Find the index of the maximum element in the array:

max_i = np.argmax(arr)
print(max_i)

Verify the maximum value in the array:

np.max(arr)

arr[max_i]

Array functions apply to 2D arrays as well (and $N$-dimensional arrays in general) with the added feature that we can choose to apply array functions to the entire array, down the columns or across the rows (or any axis).

Create a 2D NumPy array with random values and compute the sum of all the entries:

M = np.array([[2,4,2],[2,1,1],[3,2,0],[0,6,2]])
print(M)

[[2 4 2]
 [2 1 1]
 [3 2 0]
 [0 6 2]]

np.sum(M)

The function numpy.sum also takes a keyword argument axis which determines along which dimension to compute the sum:

np.sum(M,axis=0) # Sum of the columns

array([ 7, 13,  5])

np.sum(M,axis=1) # Sum of the rows

array([8, 4, 5, 8])

Mathematical Functions

Mathematical functions in NumPy are called universal functions and are vectorized. Vectorized functions operate elementwise on arrays producing arrays as output and are built to compute values across arrays very quickly. The following is a partial list of mathematical functions:


`numpy.sin`	`numpy.cos`	`numpy.tan`
`numpy.exp`	`numpy.log`	`numpy.log10`
`numpy.arcsin`	`numpy.arccos`	`numpy.arctan`

Compute the values $\sin(2 \pi x)$ for $x = 0,0.25,0.5\dots,1.75$:

x = np.arange(0,1.25,0.25)
print(x)

[0.   0.25 0.5  0.75 1.  ]

np.sin(2*np.pi*x)

array([ 0.0000000e+00,  1.0000000e+00,  1.2246468e-16, -1.0000000e+00,
       -2.4492936e-16])

We expect the array [0. 1. 0. -1. 0.] however there is (as always with floating point numbers) some rounding errors in the result. In numerical computing, we can interpret a number such as $10^{-16}$ as $0$.

Compute the values $\log_{10}(x)$ for $x = 1,10,100,1000,10000$:

x = np.array([1,10,100,1000,10000])
print(x)

[    1    10   100  1000 10000]

np.log10(x)

array([0., 1., 2., 3., 4.])

Note that we can also evaluate mathematical functions with scalar values:

np.sin(0)

0.0

NumPy also provides familiar mathematical constants such as $\pi$ and $e$:

np.pi

3.141592653589793

np.e

2.718281828459045

For example, verify the limit

$$ \lim_{x \to \infty} \arctan(x) = \frac{\pi}{2} $$

by evaluating $\arctan(x)$ for some (arbitrary) large value $x$:

np.arctan(10000)

1.5706963267952299

np.pi/2

1.5707963267948966

Random Number Generators

The subpackage numpy.random contains functions to generate NumPy arrays of random numbers sampled from different distributions. The following is a partial list of distributions:

Function	Description
`numpy.random.rand(d1,...,dn)`	Create a NumPy array (with shape `(d1,...,dn)`) with entries sampled uniformly from `[0,1)`
`numpy.random.randn(d1,...,dn)`	Create a NumPy array (with shape `(d1,...,dn)`) with entries sampled from the standard normal distribution
`numpy.random.randint(a,b,size)`	Create a NumPy array (with shape `size`) with integer entries from `low` (inclusive) to `high` (exclusive)

Sample a random number from the uniform distribution:

np.random.rand()

0.6695906195141056

Sample 3 random numbers:

np.random.rand(3)

array([0.40770395, 0.12158461, 0.72083088])

Create 2D NumPy array of random samples:

np.random.rand(2,4)

array([[0.61244625, 0.32645792, 0.55859886, 0.97613741],
       [0.57227614, 0.14315638, 0.49034299, 0.00099473]])

Random samples from the standard normal distribution:

np.random.randn()

1.4440026351051256

np.random.randn(3)

array([ 0.6172098 , -1.67631666, -2.20365265])

np.random.randn(3,1)

array([[ 0.29643549],
       [-0.44039303],
       [-1.52246126]])

Random integers sampled uniformly from various intervals:

np.random.randint(-10,10)

np.random.randint(0,2,(4,8))

array([[1, 0, 0, 1, 0, 0, 0, 0],
       [1, 1, 0, 0, 0, 0, 1, 0],
       [1, 0, 0, 0, 0, 0, 1, 0],
       [1, 0, 0, 0, 1, 1, 1, 1]])

np.random.randint(-9,10,(5,2))

array([[ 3,  9],
       [ 5, -6],
       [ 7, -7],
       [ 0, -5],
       [ 0,  0]])

Examples

Brute Force Optimization

Find the absolute maximum and minimum values of the function

$$ f(x) = x\sin(x)+\cos(4x) $$

on the interval $[0,2\pi]$. We know that the maximum and minimum values must occur at either the endpoints $x=0,2\pi$ or at critical points where $f'(x)=0$. However, the derivative is given by

$$ f'(x) = \sin(x) + x \cos(x) - 4\sin(4x) $$

and the equation $f'(x) = 0$ is impossible to solve explicitly.

Instead, create a 1D NumPy array of $x$ values from $0$ to $2\pi$ of length $N$ (for some arbitrarily large value $N$) and use the functions numpy.min and numpy.max to find maximum and minimum $y$ values, and the functions numpy.argmin and numpy.argmax to find the indices of the corresponding $x$ values.

N = 10000
x = np.linspace(0,2*np.pi,N)
y = x * np.sin(x) + np.cos(4*x)
y_max = np.max(y)
y_min = np.min(y)
x_max = x[np.argmax(y)]
x_min = x[np.argmin(y)]
print('Absolute maximum value is y =',y_max,'at x =',x_max)
print('Absolute minimum value is y =',y_min,'at x =',x_min)

Absolute maximum value is y = 2.5992726072887007 at x = 1.628136126702901
Absolute minimum value is y = -5.129752039182 at x = 5.34187001663503

Riemann Sums

Write a function called exp_int which takes input parameters $b$ and $N$ and returns the (left) Riemann sum

$$ \int_0^b e^{-x^2} dx \approx \sum_{k=0}^{N-1} e^{-x_k^2} \Delta x $$

for $\Delta x = b/N$ and the partition $x_k=k \, \Delta x$, $k=0,\dots,N$.

def exp_int(b,N):
    "Compute left Riemann sum of exp(-x^2) from 0 to b with N subintervals."
    x = np.linspace(0,b,N+1)
    x_left_endpoints = x[:-1]
    Delta_x = b/N
    I = Delta_x * np.sum(np.exp(-x_left_endpoints**2))
    return I

The infinite integral satisfies the beautiful identity

$$ \int_0^{\infty} e^{-x^2} dx = \frac{\sqrt{\pi}}{2} $$

Compute the integral with large values of $b$ and $N$:

exp_int(100,100000)

0.886726925452758

Compare to the true value:

np.pi**0.5/2

0.8862269254527579

Infinite Products

The cosine function has the following infinite product representation

$$ \cos x = \prod_{k = 1}^{\infty} \left(1 - \frac{4 x^2}{\pi^2 (2k - 1)^2} \right) $$

Write a function called cos_product which takes input parameters $x$ and $N$ and returns the $N$th partial product

$$ \prod_{k = 1}^{N} \left(1 - \frac{4 x^2}{\pi^2 (2k - 1)^2} \right) $$

def cos_product(x,N):
    "Compute the product \prod_{k=1}^N (1 - 4x^2/(pi^2 (2k - 1)^2)."
    k = np.arange(1,N+1)
    terms = 1 - 4*x**2 / (np.pi**2 * (2*k - 1)**2)
    return np.prod(terms)

Verify our function using values for which we know the result. For example, $\cos(0)=1$, $\cos(\pi)=-1$ and $\cos(\pi/4) = \frac{1}{\sqrt{2}}$.

cos_product(0,10)

1.0

cos_product(np.pi,10000)

-1.0001000050002433

cos_product(np.pi/4,10000000)

0.7071067856245614

1/2**0.5

0.7071067811865475

Matrix Multiplication

Under construction

Exercises

Exercise 1. The natural log satisfies the following definite integral

$$ \int_1^e \frac{\ln x \ dx}{(1 + \ln x)^2} = \frac{e}{2} - 1 $$

Write a function called log_integral which takes input parameters $c$ and $N$ and returns the value of the (right) Riemann sum

$$ \int_1^c \frac{\ln x \ dx}{(1 + \ln x)^2} \approx \sum_{k=1}^N \frac{\ln x_k \ \Delta x}{(1 + \ln x_k)^2} \ , \ \ \Delta x = \frac{c - 1}{N} $$

for the partition $x_k = 1 + k \Delta x$, for $k = 0, \dots , N$.

Exercise 2. Write a function called k_sum which takes input parameters k and N and returns the partial sum

$$ \sum_{n=1}^{N} \frac{1}{n^k} $$

Verify your function by comparing to the infinite series identity

$$ \sum_{n=1}^{\infty} \frac{(-1)^{n+1}}{n^2} = \frac{\pi^2}{12} $$

Exercise 3. Write a function called dot which takes 3 inputs M, i and j where M is a square NumPy array and the function returns the dot product of the $i$th row and the $j$th column of $M$.