$\newcommand{\re}{\mathbb{R}}$ Numpy provides many useful facilities to perform numerical computations including vectors, matrices and linear algebra.
import numpy
print(numpy.pi)
print(numpy.sin(numpy.pi/2))
3.141592653589793 1.0
It is common practice to import numpy like this.
import numpy as np
Note that np
acts as an alias or short-hand for numpy
.
print(np.pi)
print(np.sin(np.pi/2))
3.141592653589793 1.0
Create an array of zeros
x = np.zeros(5)
print(x)
[0. 0. 0. 0. 0.]
These one dimensional arrays are of type ndarray
type(x)
numpy.ndarray
Create an array of ones
x = np.ones(5)
print(x)
[1. 1. 1. 1. 1.]
Create and add two arrays
x = np.array([1.0, 2.0, 3.0])
y = np.array([4.0, 5.0, 6.0])
z = x + y
print(z)
[5. 7. 9.]
Get the size of array
print(len(x))
print(x.size)
3 3
Get the shape of an array
print(x.shape)
(3,)
We see that these are arrays of reals by default. We can specify the type
a = np.zeros(5, dtype=int)
print(a)
[0 0 0 0 0]
This behaves same way as Matlab's linspace function.
Generate 10 uniformly spaced numbers in [1,10]
x = np.linspace(1,10,10)
print(x)
[ 1. 2. 3. 4. 5. 6. 7. 8. 9. 10.]
Note that this includes the end points 1 and 10. The output of linspace is an ndarray
of floats.
type(x)
numpy.ndarray
x = linspace(a,b,n)
is such that x
is an array of n
elements
x[i] = a + i*h, i=0,1,2,...,n-1, h = (b-a)/(n-1)
so that
x[0] = a, x[-1] = b
x = np.arange(1,10)
print(x)
print(type(x))
[1 2 3 4 5 6 7 8 9] <class 'numpy.ndarray'>
For $m < n$, arange(m,n)
returns
$$
m, m+1, \ldots, n-1
$$
Note the last value $n$ is not included.
We can specify a step size
x = np.arange(1,10,2)
print(x)
[1 3 5 7 9]
If the arguments are float, the returned value is also float.
x = np.arange(1.0,10.0)
print(x)
[1. 2. 3. 4. 5. 6. 7. 8. 9.]
x = np.arange(0,1,0.1)
print(x)
[0. 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9]
Create an array of ones.
x = np.ones(10)
print(x)
[1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
Maybe we want set all elements to zero, so we might try this
x = 0.0
print(x)
0.0
x
has changed from an array to a scalar !!! The correct way is this.
x = np.ones(10)
print(x)
x[:] = 0.0
print(x)
[1. 1. 1. 1. 1. 1. 1. 1. 1. 1.] [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
x = np.ones(5)
y = x
x[:] = 0.0
print(x)
print(y)
[0. 0. 0. 0. 0.] [0. 0. 0. 0. 0.]
Why did y
change ? This happened because when we do
y = x
then y
is just a pointer to x
, so that changing x
changes y
also. If we want y
to be an independent copy of x
then do this
x = np.ones(5)
y = x.copy() # or y = np.copy(x)
x[:] = 0.0
print(x)
print(y)
[0. 0. 0. 0. 0.] [1. 1. 1. 1. 1.]
2-d arrays can be considered as matrices, though Numpy has a separate matrix class.
Create an array of zeros
A = np.zeros((5,5))
print(A)
[[0. 0. 0. 0. 0.] [0. 0. 0. 0. 0.] [0. 0. 0. 0. 0.] [0. 0. 0. 0. 0.] [0. 0. 0. 0. 0.]]
We can access the elements by two indices
A[1,1] = 1.0
print(A)
[[0. 0. 0. 0. 0.] [0. 1. 0. 0. 0.] [0. 0. 0. 0. 0.] [0. 0. 0. 0. 0.] [0. 0. 0. 0. 0.]]
Create an array of ones
A = np.ones((2,3))
print(A)
[[1. 1. 1.] [1. 1. 1.]]
Create identity matrix
A = np.eye(5)
print(A)
[[1. 0. 0. 0. 0.] [0. 1. 0. 0. 0.] [0. 0. 1. 0. 0.] [0. 0. 0. 1. 0.] [0. 0. 0. 0. 1.]]
Create an array by specifying its elements
A = np.array([[1.0, 2.0],
[3.0, 4.0]])
print(A)
[[1. 2.] [3. 4.]]
Create a random array and inspect its shape
m, n = 2, 3
A = np.random.rand(m,n)
print(A)
print('A.size =',A.size)
print('A.shape =',A.shape)
print('A.shape[0] =',A.shape[0],' (number of rows)')
print('A.shape[1] =',A.shape[1],' (number of columns)')
[[0.56171896 0.66235847 0.78542121] [0.40789654 0.07574282 0.7198766 ]] A.size = 6 A.shape = (2, 3) A.shape[0] = 2 (number of rows) A.shape[1] = 3 (number of columns)
To print the elements of an array, we need two for loops, one over rows and another over the columns.
for i in range(m): # loop over rows
for j in range(n): # loop over columns
print(i,j,A[i,j])
0 0 0.5617189572734113 0 1 0.6623584720511232 0 2 0.7854212087229677 1 0 0.40789654277553433 1 1 0.07574281632417834 1 2 0.7198765978001941
To transpose a 2-d array
A = np.array([[1,2],[3,4]])
print("A ="); print(A)
B = A.T
print("B ="); print(B)
A = [[1 2] [3 4]] B = [[1 3] [2 4]]
Create
$$ A = \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 2 & 0 & 0 \\ 0 & 0 & 3 & 0 \\ 0 & 0 & 0 & 4 \end{bmatrix} $$a = np.array([1,2,3,4]) # diagonal elements
A = np.diag(a)
print(A)
[[1 0 0 0] [0 2 0 0] [0 0 3 0] [0 0 0 4]]
Create tri-diagonal matrix
a = np.array([1,2,3]) # sub-diagonal : length = n-1
b = np.array([4,5,6,7]) # main diagonal : length = n
c = np.array([-1,-2,-3]) # super-diagonal: length = n-1
A = np.diag(a,-1) + np.diag(b,0) + np.diag(c,+1)
print(A)
[[ 4 -1 0 0] [ 1 5 -2 0] [ 0 2 6 -3] [ 0 0 3 7]]
Sometimes we have an existing array and we want to create a new array of the same shape and type, but without initializing its elements.
A = np.random.rand(2,3)
B = np.empty_like(A)
print(A.shape,B.shape)
(2, 3) (2, 3)
C = np.zeros_like(A)
print(C)
[[0. 0. 0.] [0. 0. 0.]]
D = np.ones_like(A)
print(D)
[[1. 1. 1.] [1. 1. 1.]]
x = np.array([1,2,3])
print(x.shape, x)
y = x.T
print(y.shape, y)
(3,) [1 2 3] (3,) [1 2 3]
Create a row vector
x = np.array([[1,2,3]]) # row vector
print('x.shape =',x.shape)
print(x)
y = x.T
print('y.shape =',y.shape)
print(y)
x.shape = (1, 3) [[1 2 3]] y.shape = (3, 1) [[1] [2] [3]]
Create a column vector
x = np.array([[1],[2],[3]]) # column vector
print('x.shape =',x.shape)
print(x)
y = x.T
print('y.shape =',y.shape)
print(y)
x.shape = (3, 1) [[1] [2] [3]] y.shape = (1, 3) [[1 2 3]]
Functions like zeros
and ones
can return row/column vector rather than array.
x = np.ones((3,1)) # column vector
print('x ='); print(x)
y = np.ones((1,3)) # row vector
print('y ='); print(y)
x = [[1.] [1.] [1.]] y = [[1. 1. 1.]]
Row/column vectors are like row/column matrix. We have to use two indices to access the elements of such row/column vectors.
Here we access elements of a column vector.
print(x[0][0])
print(x[1][0])
print(x[2][0])
1.0 1.0 1.0
x[0]
gives an array for the first row.
print('First row of x =',x[0], type(x[0]))
print('Element in first row =',x[0][0], type(x[0][0]))
First row of x = [1.] <class 'numpy.ndarray'> Element in first row = 1.0 <class 'numpy.float64'>
Array of 10 elements
x[0] | x[1] | x[2] | x[3] | x[4] | x[5] | x[6] | x[7] | x[8] | x[9] |
---|---|---|---|---|---|---|---|---|---|
0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 |
x[-10] | x[-9] | x[-8] | x[-7] | x[-6] | x[-5] | x[-4] | x[-3] | x[-2] | x[-1] |
x = np.linspace(0,9,10)
print(x)
[0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
Get elements x[2],...,x[5]
print(x[2:6])
[2. 3. 4. 5.]
Hence x[m:n]
gives the elements x[m],x[m+1],...,x[n-1]
.
Get elements from x[5]
upto the last element
print(x[5:])
[5. 6. 7. 8. 9.]
Get the last element
print(x[-1])
9.0
Get element x[5]
upto last but one element
print(x[5:-1])
[5. 6. 7. 8.]
Access every alternate element of array
print(x[0::2])
[0. 2. 4. 6. 8.]
print(x[1::2])
[1. 3. 5. 7. 9.]
These operations work on multi dimensional arrays also.
A = np.random.rand(3,4)
print(A)
[[5.67717102e-01 2.94552656e-01 1.76235501e-04 3.49796707e-01] [5.60834819e-01 7.28843246e-04 8.59005230e-01 5.02693996e-01] [6.32336995e-01 6.88768612e-01 6.63854136e-01 5.80712788e-01]]
print(A[0,:]) # 0'th row
[5.67717102e-01 2.94552656e-01 1.76235501e-04 3.49796707e-01]
print(A[:,0]) # 0'th column
[0.5677171 0.56083482 0.632337 ]
print(A[0:2,0:3]) # print submatrix
[[5.67717102e-01 2.94552656e-01 1.76235501e-04] [5.60834819e-01 7.28843246e-04 8.59005230e-01]]
A[0,:] = 0.0 # zero out zeroth row
print(A)
[[0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00] [5.60834819e-01 7.28843246e-04 8.59005230e-01 5.02693996e-01] [6.32336995e-01 6.88768612e-01 6.63854136e-01 5.80712788e-01]]
Arithmetic operations act element-wise
x = np.array([1.0, 2.0, 3.0])
y = np.array([4.0, 5.0, 6.0])
print(x*y) # multiply
[ 4. 10. 18.]
print(x/y) # divide
[0.25 0.4 0.5 ]
print(y**x) # exponentiation
[ 4. 25. 216.]
A = np.ones((3,3))
print(A*x)
[[1. 2. 3.] [1. 2. 3.] [1. 2. 3.]]
If A
and x
are arrays, then A*x
does not give matrix-vector product. For that use dot
print(A.dot(x))
[6. 6. 6.]
or equivalently
print(np.dot(A,x))
[6. 6. 6.]
In Python3 versions, we can use @
to achieve matrix operations
print(A@x)
[6. 6. 6.]
We can of course do matrix-matrix products using dot
or @
A = np.ones((3,3))
B = 2*A
print('A =\n',A)
print('B =\n',B)
print('A*B =\n',A@B)
A = [[1. 1. 1.] [1. 1. 1.] [1. 1. 1.]] B = [[2. 2. 2.] [2. 2. 2.] [2. 2. 2.]] A*B = [[6. 6. 6.] [6. 6. 6.] [6. 6. 6.]]
On vectors, @
performs dot product, i.e., x@y = dot(x,y)
x = np.array([1,1,1])
y = np.array([2,3,4])
print(x@y)
9
x = np.array([-3,-2,-1,0,1,2,3])
print('min = ',np.min(x))
print('max = ',np.max(x))
print('abs min = ',np.abs(x).min()) # np.min(np.abs(x))
print('abs max = ',np.abs(x).max()) # np.max(np.abs(x))
print('sum = ',np.sum(x))
print('dot = ',np.dot(x,x))
print('dot = ',np.sum(x*x))
min = -3 max = 3 abs min = 0 abs max = 3 sum = 0 dot = 28 dot = 28
We can compute vector norms using numpy.linalg.norm (also see scipy.linalg.norm)
from numpy.linalg import norm
print(norm(x)) # L2 norm
print(norm(x,2)) # L2 norm
print(norm(x,1)) # L1 norm
print(norm(x,np.inf)) # Linf norm
5.291502622129181 5.291502622129181 12.0 3.0
or import the linalg
module and use methods inside it.
import numpy.linalg as la
print(la.norm(x)) # L2 norm
print(la.norm(x,2)) # L2 norm
print(la.norm(x,1)) # L1 norm
print(la.norm(x,np.inf)) # Linf norm
5.291502622129181 5.291502622129181 12.0 3.0
These methods work on 2-d arrays which are square.
A = 2*np.random.rand(3,3) - 1
print(A)
print('min = ',np.min(A))
print('max = ',np.max(A))
print('abs min = ',np.abs(A).min()) # np.min(np.abs(x))
print('abs max = ',np.abs(A).max()) # np.max(np.abs(x))
print('sum = ',np.sum(A))
print('Frob norm = ',la.norm(A)) # Frobenius norm
print('L1 norm = ',la.norm(A,1)) # L1 norm
print('L2 norm = ',la.norm(A,2)) # L2 norm
print('Linf norm = ',la.norm(A,np.inf)) # Linf norm
[[ 0.40159839 -0.8591952 0.16223542] [-0.9098026 -0.51844933 -0.04381662] [-0.22818608 -0.20678772 -0.81851499]] min = -0.9098026020175445 max = 0.4015983909742702 abs min = 0.04381661776429091 abs max = 0.9098026020175445 sum = -3.020918714630895 Frob norm = 1.6700494585979953 L1 norm = 1.5844322417075245 L2 norm = 1.1504001789154477 Linf norm = 1.4720685462097143
Let $$ x, y \in \re^n, \qquad A \in \re^{n \times n} $$ To compute the matrix-vector product $y=Ax$, we can do it element-wise $$ y_i = \sum_{j=0}^{n-1} A_{ij} x_j, \qquad 0 \le i \le n-1 $$
n = 10
x = np.random.rand(n)
A = np.random.rand(n,n)
y = np.zeros(n)
for i in range(n):
for j in range(n):
y[i] += A[i,j]*x[j]
We can verify that our result is correct by computing $\| y - A x\|$
print(np.linalg.norm(y-A@x))
7.108895957933346e-16
We can also compute the product column-wise. Let $$ A_{:,j} = \textrm{j'th column of A} $$ Then the matrix-vector product can also be written as $$ y = \sum_{j=0}^{n-1} A_{:,j} x_j $$ Warning: This may have inefficient memory access since by default, numpy arrays have column-major ordering, see below.
y[:] = 0.0
for j in range(n):
y += A[:,j]*x[j]
# Now check the result
print(np.linalg.norm(y-A@x))
7.108895957933346e-16
If $A \in \re^{m\times n}$ and $B \in \re^{n \times p}$ then $C = AB \in \re^{m \times p}$ is given by $$ C_{ij} = \sum_{k=0}^{n-1} A_{ik} B_{kj}, \qquad 0 \le i \le m-1, \quad 0 \le j \le p-1 $$
m,n,p = 10,8,6
A = np.random.rand(m,n)
B = np.random.rand(n,p)
C = np.zeros((m,p))
for i in range(m):
for j in range(p):
for k in range(n):
C[i,j] += A[i,k]*B[k,j]
Let us verify the result is correct by computing the Frobenius norm
print(np.linalg.norm(C - A@B))
1.4770556661825108e-15
Another view-point is the following $$ C_{ij} = (\textrm{i'th row of A}) \cdot (\textrm{j'th column of B}) $$
for i in range(m):
for j in range(p):
C[i,j] = A[i,:] @ B[:,j]
# Now check the result
print(np.linalg.norm(C - A@B))
0.0
Given $m$-vector $a$ and $n$-vector $b$, both considered as column vectors, their outer product
$$ A = a b^\top = \begin{bmatrix} | & | & & | \\ a b_0 & a b_1 & \ldots & a b_{n-1} \\ | & | & & | \end{bmatrix} $$is an $m \times n$ matrix with elements
$$ A_{ij} = a_j b_j, \qquad 0 \le i \le m-1, \quad 0 \le n \le n-1 $$a = np.array([1,2,3])
b = np.array([1,2])
A = np.outer(a,b)
print(A.shape)
print(A)
(3, 2) [[1 2] [2 4] [3 6]]
Numpy provides standard functions like sin
, cos
, log
, etc. which can act on arrays in an element-by-element manner. This is not the case for functions in math
module, which can only take scalar arguments.
x = np.linspace(0.0, 2.0*np.pi, 5)
y = np.sin(x)
print('x =',x)
print('y =',y)
x = [0. 1.57079633 3.14159265 4.71238898 6.28318531] y = [ 0.0000000e+00 1.0000000e+00 1.2246468e-16 -1.0000000e+00 -2.4492936e-16]
A = np.random.rand(5,5) # 5x5 matrix
Y = np.sin(A)
print('A =\n',A)
print('Y =\n',Y)
A = [[0.69799961 0.44193454 0.42009594 0.59430268 0.44412881] [0.88453597 0.60609966 0.57272725 0.68715342 0.16661875] [0.5920093 0.63078135 0.78454055 0.81219262 0.73284029] [0.27186644 0.51657847 0.28403543 0.18333706 0.2893531 ] [0.59918434 0.9679843 0.18048533 0.82641006 0.26114357]] Y = [[0.64268642 0.42768894 0.40784805 0.55993113 0.42967137] [0.77362104 0.5696662 0.54192611 0.63433919 0.16584888] [0.5580295 0.58977593 0.7065001 0.72579724 0.66898346] [0.26852979 0.49390796 0.28023166 0.18231171 0.28533228] [0.56396909 0.82374457 0.17950705 0.73550386 0.25818551]]
By default, the ordering is same as in C/C++, the inner-most index is the fastest running one. For example, if we have an array of size (2,3), they are stored in memory in this order
a[0,0], a[0,1], a[0,2], a[1,0], a[1,1], a[1,2]
i.e.,
a[0,0] --> a[0,1] --> a[0,2]
|
|----------<----------|
|
a[1,0] --> a[1,1] --> a[1,2]
a = np.array([[1,2,3], [4,5,6]])
print('Row contiguous =', a[0,:].data.contiguous)
print('Col contiguous =', a[:,0].data.contiguous)
a.flags
Row contiguous = True Col contiguous = False
C_CONTIGUOUS : True F_CONTIGUOUS : False OWNDATA : True WRITEABLE : True ALIGNED : True WRITEBACKIFCOPY : False
To get fortran style ordering, where the outer-most index is the fastest running one, which corresponds to the following layout in memory
a[0,0], a[1,0], a[0,1], a[1,1], a[0,2], a[1,2]
i.e.,
a[0,0] -- a[0,1] -- a[0,2]
| | | | |
v ^ v ^ v
| | | | |
a[1,0] -- a[1,1] -- a[1,2]
create like this
b = np.array([[1,2,3], [4,5,6]], order='F')
print('Row contiguous =', b[0,:].data.contiguous)
print('Col contiguous =', b[:,0].data.contiguous)
b.flags
Row contiguous = False Col contiguous = True
C_CONTIGUOUS : False F_CONTIGUOUS : True OWNDATA : True WRITEABLE : True ALIGNED : True WRITEBACKIFCOPY : False
x = np.linspace(0,3,4)
y = np.linspace(0,2,3)
X, Y = np.meshgrid(x,y)
print('len(x) = ',len(x))
print('len(y) = ',len(y))
print('shape X= ',X.shape)
print('shape Y= ',Y.shape)
print('x = ', x)
print('y = ', y)
print('X = '); print(X)
print('Y = '); print(Y)
len(x) = 4 len(y) = 3 shape X= (3, 4) shape Y= (3, 4) x = [0. 1. 2. 3.] y = [0. 1. 2.] X = [[0. 1. 2. 3.] [0. 1. 2. 3.] [0. 1. 2. 3.]] Y = [[0. 0. 0. 0.] [1. 1. 1. 1.] [2. 2. 2. 2.]]
If $x \in \re^m$ and $y \in \re^n$, then the output is arranged like this $$ X[i,j] = x[j], \qquad Y[i,j] = y[i], \qquad 0 \le i \le n-1, \quad 0 \le j \le m-1 $$ If we want the following arrangement $$ X[i,j] = x[i], \qquad Y[i,j] = y[j], \qquad 0 \le i \le m-1, \quad 0 \le j \le n-1 $$ we have to do the following
X, Y = np.meshgrid(x,y,indexing='ij')
print('len(x) = ',len(x))
print('len(y) = ',len(y))
print('shape X= ',X.shape)
print('shape Y= ',Y.shape)
len(x) = 4 len(y) = 3 shape X= (4, 3) shape Y= (4, 3)
The second form is useful when working with finite difference schemes on Cartesian grids, where we want to use i index running along x-axis and j index running along y-axis.
A = np.array([[1,2,3],
[4,5,6]])
print(A)
[[1 2 3] [4 5 6]]
B = np.reshape(A,2*3,order='C')
print(B)
[1 2 3 4 5 6]
A1 = np.reshape(B,(2,3),order='C')
print(A1)
[[1 2 3] [4 5 6]]
C = np.reshape(A,2*3,order='F')
print(C)
[1 4 2 5 3 6]
A2 = np.reshape(C,(2,3),order='F')
print(A2)
[[1 2 3] [4 5 6]]
Write an array to file
A = np.random.rand(3,3)
print(A)
np.savetxt('A.txt',A)
[[0.23800817 0.96267486 0.42996244] [0.31093536 0.61795882 0.65229253] [0.44645543 0.3532086 0.88249008]]
Check the contents of the file in your terminal
cat A.txt
We can also print it from within the notebook
!cat A.txt
2.380081667929159206e-01 9.626748557062282385e-01 4.299624366157210886e-01 3.109353645036705416e-01 6.179588184832025544e-01 6.522925255033611425e-01 4.464554342116026087e-01 3.532086017204271178e-01 8.824900753021190924e-01
Write two 1-D arrays as columns into file
x = np.array([1.0,2.0,3.0,4.0])
y = np.array([2.0,4.0,6.0,8.0])
np.savetxt('data.txt',np.column_stack([x,y]))
!cat data.txt
1.000000000000000000e+00 2.000000000000000000e+00 2.000000000000000000e+00 4.000000000000000000e+00 3.000000000000000000e+00 6.000000000000000000e+00 4.000000000000000000e+00 8.000000000000000000e+00
We can control the number of decimals saved, and use scientific notation
np.savetxt('data.txt',np.column_stack([x,y]),fmt='%8.4e')
!cat data.txt
1.0000e+00 2.0000e+00 2.0000e+00 4.0000e+00 3.0000e+00 6.0000e+00 4.0000e+00 8.0000e+00
We can read an existing file like this
d = np.loadtxt('data.txt')
print('Shape d =', d.shape)
x1 = d[:,0]
y1 = d[:,1]
print('x =',x1)
print('y =',y1)
Shape d = (4, 2) x = [1. 2. 3. 4.] y = [2. 4. 6. 8.]
Loops are slow in Python and it is good to avoid them. The following example show how to time code execution.
Create some random matrices
m = 1000
A = np.random.random((m,m))
B = np.random.random((m,m))
First we add two matrices with an explicit for loop.
%%timeit -r10 -n10
C = np.empty_like(A)
for i in range(m):
for j in range(m):
C[i,j] = A[i,j] + B[i,j]
230 ms ± 5.1 ms per loop (mean ± std. dev. of 10 runs, 10 loops each)
Now we use the + operator.
%%timeit -r10 -n10
C = A + B
1.1 ms ± 362 µs per loop (mean ± std. dev. of 10 runs, 10 loops each)
The loop version is significantly slower.