12. Array operations in numpy#
Saleh Rezaeiravesh, saleh.rezaeiravesh@manchester.ac.uk
Department of Mechanical and Aerospace Engineering, The University of Manchester, Manchester, UK
12.1. Intended Learning Outcomes#
After reading and running this notebook, the students should be able to:
Apply basic operations (e.g. addition, subtraction, multiplication) to
numpyarrays.Use
numpyto compute the transpose, summation of elements, etc of arrays.Use
numpyto evaluate arithmetic expressions involving arrays.
12.2. Import numpy#
At the beginning of your Jupyter notebook or Python script script (can be saved as a .py file), we import numpy to have access to its classes, functions, and attributes. The imported library can be named by anything, here we choose np.
import numpy as np
12.3. Addition and subtraction#
As we know to add or subtract two or more matrices, the shape (dimension and size) of the matrices must be the same. Keeping this in mind, we can simply use + and - signs for addition and subtraction of numpy arrays, respectively.
12.3.1. 1D Arrays#
Consider two vectors (1d arrays) \(a_1\) and \(a_2\):
a1=np.array([-1,3,4,5])
a2=np.array([5,-1,2,8])
For addition and subtraction of these two vectors in numpy, we have:
s = a1 + a2
d = a1 - a2
print('a1+a2=',s)
print('a1-a2=',d)
a1+a2= [ 4 2 6 13]
a1-a2= [-6 4 2 -3]
The results are vectors with the same shape as a1 and a2. Note that the operations have been applied to corresponding elements of \(a_1\) and \(a_2\).
12.3.2. 2D Arrays#
The same procedure is applicable to multidimensional arrays. For instance, for 2D matrices \(\mathbf{A}\) and \(\mathbf{B}\) of size \(m\times n\), we have \(\mathbf{C} = \mathbf{A}+\mathbf{B}\) where
Let’s consider \(\mathbf{B}_1\) and \(\mathbf{B}_2\) that are both \(2\times 3\) matrices:
The addition and subtraction of these two arrays are written as B1+B2 and B1-B2, respectively.
B1 = np.array([[-1,3,4],[2,3,6]])
B2 = np.array([[2,-5,1],[12,2,4]])
S1 = B1+B2
D1 = B1-B2
print('B1+B2=',S1)
print('B1-B2=',D1)
B1+B2= [[ 1 -2 5]
[14 5 10]]
B1-B2= [[ -3 8 3]
[-10 1 2]]
Example: Instead of using the built-in operator +, write a Python function that returns the summation of two 2D matrices.
For this, we write a nested loop to add corresponding elements of arrays \(\mathbf{A}\) and \(\mathbf{B}\) and assign the result to a new matrix \(\mathbf{C}\):
def sumArray(A,B):
"""
Returns sum of mxn matrices A and B
Args:
`A`: 2D numpy array of size m x n
`B`: 2D numpy array of size m x n
Returns
`C`: 2D numpy array of size m x n
"""
#Get the size of matrix A
m=A.shape[0]
n=A.shape[1]
if m != B.shape[0] or n != B.shape[1]: #check if A and B have the same size
raise ValueError('A and B should have the same size!')
# initialise matrix C where C=A+B
C = np.zeros((m,n))
# do the addition elementwise
for i in range(m): #loop over the rows of A and B
for j in range(n): #loop over the columns of A and B
C[i,j] = A[i,j] + B[i,j]
return C
Let’s test this function and compare its output to the built-in numpy operation +:
print('B1+B2 (direct implementation)=',sumArray(B1,B2))
print('B1+B2 (numpy) =',B1+B2)
B1+B2 (direct implementation)= [[ 1. -2. 5.]
[14. 5. 10.]]
B1+B2 (numpy) = [[ 1 -2 5]
[14 5 10]]
Discussion: From the above implementation, can you say how many elementary operations are required for computing \(\mathbf{A}+\mathbf{B}\), where both matrices are of size \(m\times n\)?
12.4. Matrix transpose#
A transpose of a matrix is obtained by switching its rows and columns. To find the transpose of a numpy array A, we can use either of the following commands:
A.Tnumpy.transpose(A)
For instance, for the \(2\times 3\) array \(\mathbf{B}_1\) defined as
we should get
Let’s check this:
print(B1)
[[-1 3 4]
[ 2 3 6]]
print(B1.T)
print(np.transpose(B1))
[[-1 2]
[ 3 3]
[ 4 6]]
[[-1 2]
[ 3 3]
[ 4 6]]
Clearly, the transposed matrix has the size \(3\times 2\):
print('B1 shape:',B1.shape)
print('B1.T shape:',B1.T.shape)
B1 shape: (2, 3)
B1.T shape: (3, 2)
Note: The transpose of a 1D numpy array (vector) is the same as the array. If we want to convert a row vector to a column one, or vice versa, we should first make the vector two-dimensional and then transpose it.
Let’s consider a1 that is a vector:
print('a1 as a vector:', a1)
print(a1.shape) #shape of a1
print(a1.T.shape) #shape of a1.T
a1 as a vector: [-1 3 4 5]
(4,)
(4,)
Clearly, the shape of the array has not changed. Therefore, we, first, make it two dimensional (the result is a column matrix):
print('a1 as a column matrix:')
aa1 = a1[:,None]
print(aa1)
print(aa1.shape)
a1 as a column matrix:
[[-1]
[ 3]
[ 4]
[ 5]]
(4, 1)
Now, the transpose of this is:
print('a1 as a row matrix:')
aa1T = aa1.T
print(aa1T)
print(aa1T.shape)
a1 as a row matrix:
[[-1 3 4 5]]
(1, 4)
Example: Above, we learned how to use a built-in numpy operator to find the transpose of a matrix. Write a Python function that switches the rows and columns of a 2D matrix and return its transpose.
Assume A is a m by n array. We write a nested loop over rows and columns of A, and assign A[i,j] to B[j,i], where \(B=A^T\):
def trans(A):
m = A.shape[0] #number of rows of A
n = A.shape[1] #number of columns of A
B = np.zeros((n,m)) #initialize the transpose of A
for i in range(m): #loop over rows of A
for j in range(n): #loop over rows of B
B[j,i] = A[i,j]
return B
The output of this function is the same as numpy transpose:
print(B1.T) #numpy transpose
print(trans(B1)) #our implementation
[[-1 2]
[ 3 3]
[ 4 6]]
[[-1. 2.]
[ 3. 3.]
[ 4. 6.]]
12.4.1. Conjugate transpose#
Consider an \(m\times n\) matrix \(\mathbf{A}\) with complex elements. Its conjugate or Hermitian transpose is \(\mathbf{A}^H\) of size \(n\times m\), where \(a^H_{ij} = \bar{a}_{ji}\) where \(\bar{}\) is the complex conjugate operator.
Recall: The complex conjugate pair of \(z=a+bj\) is \(\bar{z}=a-bj\).
In numpy, we can compute the complex conjugate of an array by first applying numpy.conjugate() and then transpose. See the example below.
C = np.array([[1., 2.-3.j, 4.],[5.+2.j,1.+6.j, 3-2.j]])
print("C = ")
print(C, C.shape)
print("C.H = ")
CH= np.conjugate(C).T
print(CH, C.shape)
C =
[[1.+0.j 2.-3.j 4.+0.j]
[5.+2.j 1.+6.j 3.-2.j]] (2, 3)
C.H =
[[1.-0.j 5.-2.j]
[2.+3.j 1.-6.j]
[4.-0.j 3.+2.j]] (2, 3)
12.5. Multiplication of arrays#
For 2D arrays \(\mathbf{A}\) and \(\mathbf{B}\) of size \(m\times k\) and \(k\times n\), respectively. The multiplication \(A\) by \(B\) is a 2D matrix of size \(m\times n\), where
12.5.1. numpy operator#
For multiplication of arrays in numpy, we use the operator @. But, we must be careful about checking the shape of arrays which are going to be multiplied to make sure the operation is mathematically possible (i.e. there is no mismatch between the relevant sizes).
This means, for instance, if \(\mathbf{A}\) is \(2\times 3\) and \(\mathbf{B}\) is \(3\times 4\), then only \(\mathbf{A\,B}\) can be done, and not \(\mathbf{B\,A}\):
A = np.array([[2.,-5.,1.],[12.,2.,4.]])
B = np.array([[4.,2.,-7.,11.],[-2.,2.,6.,0.],[19.,-3.,-8.,4.]])
print('Shape of A:',A.shape)
print('Shape of B:',B.shape)
Shape of A: (2, 3)
Shape of B: (3, 4)
\(\mathbf{C = A\,B}\) should be a \(2\times 4\) matrix:
C = A @ B
print(C)
print('Shape of AB:',C.shape)
[[ 37. -9. -52. 26.]
[ 120. 16. -104. 148.]]
Shape of AB: (2, 4)
But if we try \(\mathbf{BA}\), an error message will appear saying there is a “mismatch in the dimensions”:
print(B@A)
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[19], line 1
----> 1 print(B@A)
ValueError: matmul: Input operand 1 has a mismatch in its core dimension 0, with gufunc signature (n?,k),(k,m?)->(n?,m?) (size 2 is different from 4)
Note: The operator @ is a shorthand for numpy.matmul(...). Their results are the same:
C2 = np.matmul(A,B)
print(C2)
[[ 37. -9. -52. 26.]
[ 120. 16. -104. 148.]]
12.5.2. Direct implementation#
Above, we use the Python built-in operator @ for matrix-matrix multiplication.
Write a Python function that multiplies two input 2D matrix \(\mathbf{A}\) and \(\mathbf{B}\) without using numpy built-in operators. Your function should check if the size of the matrices match and allow their multiplication.
We follow the definition provided above for multiplication of 2D arrays \(\mathbf{A}\) and \(\mathbf{B}\) of size \(m\times k\) and \(k\times n\), respectively
def matrixMultip(A,B):
"""
Returns multplication of A by B, where A and B are 2D arrays
"""
m = A.shape[0] #number of rows of A
n = B.shape[1] #number of columns of B
k = A.shape[1] #number of columns of A OR rows of B
#check if the size of A and B match for computing AB
if A.shape[1] != B.shape[0]:
raise ValueError("Dimensions of A and B do not match for multiplication of A by B!")
#initialise matrix C
C = np.zeros((m,n)) # C is m x n
#do the matrix-matrix multiplication elementwise
for i in range(m): #loop over rows of C and A
for j in range(n): #loop over rows of columns of C and B
for l in range(k): #loop over columns of A = rows of B
C[i,j] = C[i,j] + A[i,l] * B[l,j]
return C
First, validate the result of our function compared to the numpy operator @. We can investigate this for two random matrices:
A = np.random.rand(4,2)
B = np.random.rand(2,3)
Cn = A @ B
Cm = matrixMultip(A,B)
print(Cn)
print(Cn.shape)
[[0.26593672 0.74907718 0.49796522]
[0.05690651 0.17315821 0.11126285]
[0.37289859 0.58039038 0.52637157]
[0.55859206 0.73024624 0.73759543]]
(4, 3)
print(Cm)
print(Cm.shape)
[[0.26593672 0.74907718 0.49796522]
[0.05690651 0.17315821 0.11126285]
[0.37289859 0.58039038 0.52637157]
[0.55859206 0.73024624 0.73759543]]
(4, 3)
Next, let’s compare the cost of computing multiplication of \(\mathbf{A}\) by \(\mathbf{B}\) using our function and the numpy. To this end, we create arrays of arbitrary sizes, and measure the system time before and after running the multiplication. For this, we import library time.
import time as tm
A = np.random.rand(100,150)
B = np.random.rand(150,200)
numpy operation:
t1 = tm.time() #system time before running A@B
Cn = A @ B
t2 = tm.time() #system time after running A@B
Calling our function (Note: depending on the size of the arrays, it may take a while):
t3 = tm.time() #system time before calling matrixMultip
Cm = matrixMultip(A,B)
t4 = tm.time() #system time before calling matrixMultip
print('Time for multuiplication of A by B using numpy:', t2-t1)
print('Time for multuiplication of A by B using our function:', t4-t3)
Time for multuiplication of A by B using numpy: 0.0011703968048095703
Time for multuiplication of A by B using our function: 0.9372291564941406
Clearly, our implementation takes much longer than the numpy operation. This is becasue of the nested loop in our function.
Exercise: Repeat the above experiment for different sizes of the arrays. How does the elapsed times scale with the size of the arrays?
Matrix vector multiplication
If \(\mathbf{A}\) is a \(m\times n\) matrix and \(\textbf{b}\) is \(n\times 1\), then \(\mathbf{C} = \mathbf{Ab}\) is a \(m\times 1\) matrix.
Consider, for instance, \(\mathbf{A}\) is a \(2\times 3\) matrix, and \(b\) is a vector of size \(3\). For implementing the multiplication \(\mathbf{A \,b}\) in numpy, A@b works if b is either a \(3\times 1\) column vector of just a 1D array of size \(3\). For these two cases, the resulting A @ b has the shape \(2\times 1\) and \(2\), respectively.
i) if \(b\) is a 1D array:
A = np.array([[2,-5,1],[12,2,4]])
b = np.array([-4,7,1]) #1D array
print('Shape of A:',A.shape)
print('Shape of b:',b.shape)
Shape of A: (2, 3)
Shape of b: (3,)
Let’s compute A @ b:
c = A @ b
print(c)
print('Ab shape:',c.shape)
[-42 -30]
Ab shape: (2,)
ii) if \(b\) is a 2D array:
b = np.array([-4,7,1])[:,None] #2D array 3 x1
print('Shape of A:',A.shape)
print('Shape of b:',b.shape)
Shape of A: (2, 3)
Shape of b: (3, 1)
c2 = A @ b
print(c2)
print('Ab shape:',c2.shape)
[[-42]
[-30]]
Ab shape: (2, 1)
Exercise: Implement the matrix-vector multiplication in a function. Compare your results and its cost with the case of using @ operator.
12.6. Inner (dot) product#
In many cases, we would like to have an inner or dot product of arrays. But, we should note that:
For 1D arrays (vectors) \(\mathbf{a}\) and \(\mathbf{b}\) of equal lengths,
numpy.dot(a,b)is a scalar (exactly like the mathematical dot product). A similar result can be obtained froma @ b(See Example 1 below).For higher-dimensional matrices \(\mathbf{A}\) and \(\mathbf{B}\),
numpy.dot(A,B)is identical toA @ B(see Example 2 below).
Example 1: dot product of two vectors (1D arrays).
a=np.array([1,2,4])
b=np.array([4,2,7])
c=np.dot(a,b)
print('shape of a and b:',a.shape,b.shape)
print('dot(a,b) =',c)
print('a@b =', a @ b)
shape of a and b: (3,) (3,)
dot(a,b) = 36
a@b = 36
Example 2: dot product of two 2D arrays
A = np.array([[2,-5,1],[12,2,4]])
B = np.array([[4,2,-7,11],[-2,2,6,0],[19,-3,-8,4]])
print('Shape of A:',A.shape)
print('Shape of B:',B.shape)
C1=np.dot(A,B) #for multi-dim arrays, dot products acts the same as @
print(C1)
print('Shape of dot(A,B)',C1.shape)
print('A @ B =',A @ B)
Shape of A: (2, 3)
Shape of B: (3, 4)
[[ 37 -9 -52 26]
[ 120 16 -104 148]]
Shape of dot(A,B) (2, 4)
A @ B = [[ 37 -9 -52 26]
[ 120 16 -104 148]]
12.7. Minimum, maximum, and summation#
Let A be a numpy array of any dimension;
numpy.max(A)returns the element ofAwith the maximum valuenumpy.min(A): returns the element ofAwith the minimum valuenumpy.sum(A): returns the sum of the elements ofA
As explained in Arrays in numpy notebook, we can specify the axis along which the minimum element, maximum element, and sum of the elements of a numpy array are computed. Recall that for a 2D array:
axis=0: Operates across the rows (vertically).axis=1: Operates across the columns (horizontally).
Example: apply the above operators to 1D arrays
a = np.random.rand(10)
print('a =',a)
print('min(a)=',np.min(a))
print('max(a)=',np.max(a))
print('sum(a)=',np.sum(a))
a = [0.70645876 0.19254541 0.83272147 0.77627379 0.08946537 0.80463202
0.73221394 0.5177645 0.03357301 0.11236578]
min(a)= 0.033573011901882
max(a)= 0.8327214703665241
sum(a)= 4.798014056797415
Example: Apply the above operators to 2D arrays:
A = np.array([[1., -5, 6],[-1, 2, 8]])
print(A)
print('min(A)=',np.min(A))
print('max(A)=',np.max(A))
print('sum(A)=',np.sum(A))
[[ 1. -5. 6.]
[-1. 2. 8.]]
min(A)= -5.0
max(A)= 8.0
sum(A)= 11.0
Let’s repeat the above commands along axis=0 of A, to find its minimum/maximum element and sum of its elements across the rows:
print('min(A)=',np.min(A,axis=0))
print('max(A)=',np.max(A,axis=0))
print('sum(A)=',np.sum(A,axis=0))
min(A)= [-1. -5. 6.]
max(A)= [1. 2. 8.]
sum(A)= [ 0. -3. 14.]
The above applied across the columns of A:
print('min(A)=',np.min(A,axis=1))
print('max(A)=',np.max(A,axis=1))
print('sum(A)=',np.sum(A,axis=1))
min(A)= [-5. -1.]
max(A)= [6. 8.]
sum(A)= [2. 9.]
12.8. Scalar-array operations#
We can simply multiply or divide a numpy array of any shape by a scalar using * and /, respectively. By these, all elements of the array are multiplied/divided by the scalar. The shape of the resulting array is the same as the original one.
Example: 1D array
print('a:',a)
print('5*a:', 5*a)
print('a/2.3:', a/2.3)
a: [0.70645876 0.19254541 0.83272147 0.77627379 0.08946537 0.80463202
0.73221394 0.5177645 0.03357301 0.11236578]
5*a: [3.5322938 0.96272705 4.16360735 3.88136896 0.44732687 4.0231601
3.6610697 2.5888225 0.16786506 0.5618289 ]
a/2.3: [0.30715598 0.0837154 0.36205281 0.33751034 0.03889799 0.34984001
0.31835389 0.225115 0.01459696 0.04885469]
Example: 2D array
print('A:',A)
print('5*A:', 5*A)
print('A/2.3:', A/2.3)
A: [[ 1. -5. 6.]
[-1. 2. 8.]]
5*A: [[ 5. -25. 30.]
[ -5. 10. 40.]]
A/2.3: [[ 0.43478261 -2.17391304 2.60869565]
[-0.43478261 0.86956522 3.47826087]]
Note: In contrast to Matlab, we do not need to use .* or ./ for the above operations.
We can combine different operations applied to numpy arrays.
For instance, Consider \(\mathbf{A}\) and \(\mathbf{B}\) are two \(3\times 2\) arrays. We can implement the mathematical expression
in numpy simply by writing C=(4.3*A+B/3)@A.T.
This shows how high-level, prctical and nice numpy is :)
Example
A = np.array([[2,-5,1],[12,2,4]])
B = np.array([[-12,4,-1],[6,8,-7]])
C = (4.3*A + B/3) @ A.T
print(C)
print(C.shape)
[[114. 30.73333333]
[ 65.73333333 725.2 ]]
(2, 2)