1.10. Linear Algebra with Numpy and Scipy#

Reference: Chapter 1 of Computational Nuclear Engineering and Radiological Science Using Python, R. McClarren (2018)

1.10.1. Learning Objectives#

After studying this notebook, completing the activities, and asking questions in class, you should be able to:

Create, manipulate, and use NumPy arrays
Explain scoping rules for arrays in Python
Perform element-wise and matrix operations with arrays
Access elements of an array with indices and slices
Iterate over elements of an array with for loops

1.10.2. Getting Started with NumPy Arrays#

So far, we have discussed lists, which are containers that hold any type of Python object such as strings, floats, integers, other lists, etc.

As engineers, we often want to store floating point numbers in vectors and matrices to perform linear algebra calculations. In this class, we will use NumPy. For a less engineering specific tutorial on NumPy, see https://docs.scipy.org/doc/numpy-1.15.0/user/quickstart.html

The basic unit in numpy is a multi-dimensional array:

A one-dimensional (1-D) array is a vector
A 2-D array is a matrix
A 3-D array is a vector of matrices
etc

We’ll start by loading the numpy module (a.k.a. library).

# Import the numpy module (a.k.a library) and give it the "nickname" np
# Warning: Do not modify this line. The autograde assumes numpy is loades as 'np'
import numpy as np

We can now write np. to call functions in the numpy library (a.k.a. module). You’ll see we place a few standard import statements at the top of most notebooks in this class.

# Create a vector
a_vector = np.array([1,2,3,4])

# Create a matrix
a_matrix = np.array([(1,2,3),(4,5,6),(7,8,9)])

# Print to the screen
print("The vector",a_vector)
print("The matrix\n",a_matrix)

The vector [1 2 3 4]
The matrix
 [[1 2 3]
 [4 5 6]
 [7 8 9]]

Arrays have several attributes that you can use to find out information regarding the array you’re working with.

#shape tells you the shape
print("The shape of a_vector is ", a_vector.shape)
print("The shape of a_matrix is ", a_matrix.shape)

The shape of a_vector is  (4,)
The shape of a_matrix is  (3, 3)

#ndim tells you the dimensionality of an array
print("The dimension of a_vector is ", a_vector.ndim)
print("The dimension of a_matrix is ", a_matrix.ndim)

The dimension of a_vector is  1
The dimension of a_matrix is  2

#size is the total number of elements = 
#the product of the number of elements in each dimension
print("The size of a_vector is ", a_vector.size,"= ",
      a_vector.shape[0])
print("The size of a_matrix is ", a_matrix.size,"=",
      a_matrix.shape[0],"*",a_matrix.shape[1])

The size of a_vector is  4 =  4
The size of a_matrix is  9 = 3 * 3

Notice that a_matrix.shape is a list. We can access the number of rows with a_matrix.shape[0] and the number of columns with a_matrix.shape[1].

Home Activity

Create the array shown below and store it in the variable my_array.

\[\begin{split} A = \begin{bmatrix} 1 & 2 \\ 3 & 1 \\ 9 & 2 \end{bmatrix} \end{split}\]

# Add your solution here

# Removed autograder test. You may delete this cell.

You can also re-size an array after you create it:

A = np.array([2,4,6,8])
print("A is now a vector",A,"\n")
A = A.reshape(2,2)
print("A is now a matrix\n",A)
print("\nNow let's shape A back to a vector")
print(A.reshape(4))

A is now a vector [2 4 6 8] 

A is now a matrix
 [[2 4]
 [6 8]]

Now let's shape A back to a vector
[2 4 6 8]

Notice how I needed to assign A with the reshaped array. We’ll talk more about this later.

We can resize my_array by stacking the 1st row followed by the 2nd row followed by the 3rd row:

print("Before:")
print(my_array)

print("\nAfter stacking:")
print(my_array.reshape(6))

Before:
[[1 2]
 [3 1]
 [9 2]]

After stacking:
[1 2 3 1 9 2]

We can resize my_array by stacking the 2nd column under the first column.

my_array.reshape(6,order='F')

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

I want to highlight two things in the above example:

NumPy does not distinguish between a row vector and a column vector.
Specifying order='F' stacks by columns, which seems strange at first glance.

Regarding 2., a quick check of the documentation, https://docs.scipy.org/doc/numpy/reference/generated/numpy.reshape.html, reveals that order='C' is the default. Moreover, 'C' if shorthand for C-style stacking and 'F' if short for Fortan-style stacking. C and Fortran are two mature (read “old”) and popular programming languages for scientific computing. They use different conventions for storing arrays in memory, which means different stacking styles. The main take away is, when in doubt, peek at the documentation.

1.10.3. Creating Arrays in Neat Ways#

The arange is a Numpy variant of the range we saw earlier.

#let's make a vector from 0 to 2*pi in intervals of 0.1
dx = 0.1
X = np.arange(0,2*np.pi,dx)
print(X)

[0.  0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.  1.1 1.2 1.3 1.4 1.5 1.6 1.7
8 1.9 2.  2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3.  3.1 3.2 3.3 3.4 3.5
6 3.7 3.8 3.9 4.  4.1 4.2 4.3 4.4 4.5 4.6 4.7 4.8 4.9 5.  5.1 5.2 5.3
4 5.5 5.6 5.7 5.8 5.9 6.  6.1 6.2]

Typically, however, the linspace function is actually what you want.

X = np.linspace(start = 0, stop = 2*np.pi, num = 62)
print(X)

[0.         0.10300304 0.20600608 0.30900911 0.41201215 0.51501519
61801823 0.72102126 0.8240243  0.92702734 1.03003038 1.13303342
23603645 1.33903949 1.44204253 1.54504557 1.64804861 1.75105164
85405468 1.95705772 2.06006076 2.16306379 2.26606683 2.36906987
47207291 2.57507595 2.67807898 2.78108202 2.88408506 2.9870881
09009113 3.19309417 3.29609721 3.39910025 3.50210329 3.60510632
70810936 3.8111124  3.91411544 4.01711848 4.12012151 4.22312455
32612759 4.42913063 4.53213366 4.6351367  4.73813974 4.84114278
94414582 5.04714885 5.15015189 5.25315493 5.35615797 5.459161
56216404 5.66516708 5.76817012 5.87117316 5.97417619 6.07717923
18018227 6.28318531]

Notice how it starts and ends exactly where I told it to.

If you want an array of zeros or ones, NumPy has a function for you.

zero_vector = np.zeros(10) #vector of length 10
zero_matrix = np.zeros((4,4)) #4 by 4 matrix
print("The zero vector:",zero_vector)
print("The zero matrix\n",zero_matrix)

The zero vector: [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
The zero matrix
 [[0. 0. 0. 0.]
 [0. 0. 0. 0.]
 [0. 0. 0. 0.]
 [0. 0. 0. 0.]]

ones_vector = np.ones(10) #vector of length 10
ones_matrix = np.ones((4,4)) #4 by 4 matrix
print("The ones vector:",ones_vector)
print("The ones matrix\n",ones_matrix)

The ones vector: [1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
The ones matrix
 [[1. 1. 1. 1.]
 [1. 1. 1. 1.]
 [1. 1. 1. 1.]
 [1. 1. 1. 1.]]

Home Activity

Create a 6 by 5 matrix of zeros and store in variable zeros65

# Add your solution here

# Removed autograder test. You may delete this cell.

You can also create an array with random entries between 0 and 1:

random_matrix = np.random.rand(2,3) #random 2 x 3 matrix
print("Here's a random 2 x 3 matrix\n",random_matrix)

print("\nAnother example")

#make a random array between two numbers
print(np.random.uniform(low=-5,high=6,size=(3,3)))

#make random integers

print("\nA third example")
print(np.random.randint(low=1,high=6,size=10))

Here's a random 2 x 3 matrix
 [[0.65069491 0.98923999 0.10681656]
 [0.33684833 0.28231416 0.19195413]]

Another example
[[-4.05548427  4.47772229 -4.60616325]
 [ 5.63124228 -4.07368184 -2.71655742]
 [ 0.18000039  2.4395264   1.68427497]]

A third example
[4 1 4 3 2 2 4 3 2 2]

Finally, you may want to create an identity matrix:

#3 x 3 identity matrix
identity33 = np.identity(3)
print(identity33)

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

Home Activity

Read the documentation for vstack and hstack. Then create the following array in Python and store it in variable M. Hint: How can you use vstack, identity, and zeros together?

https://docs.scipy.org/doc/numpy/reference/generated/numpy.vstack.html

https://docs.scipy.org/doc/numpy/reference/generated/numpy.hstack.html

\[\begin{split} M = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \\ 0 & 0 & 0 \end{bmatrix}\end{split}\]

# Add your solution here

# Removed autograder test. You may delete this cell.

1.10.4. Operations on Arrays#

NumPy also defines arithmetic operations mostly in the way that you would expect.

#vector addition
x = np.ones(3) #3-vector of ones
y = 3*np.ones(3)-1 #3-vector of 2's
print(x,"+",y,"=",x+y)
print(x,"-",y,"=",x-y)

[1. 1. 1.] + [2. 2. 2.] = [3. 3. 3.]
[1. 1. 1.] - [2. 2. 2.] = [-1. -1. -1.]

Multiplication and division are element-wise, which may not be what you expect, but it is helpful.

array([1., 1., 1.])

y = np.array([1.0,2.0,3.0])
print(x,"*",y,"=",x*y)
print(x,"/",y,"=",x/y)

[1. 1. 1.] * [1. 2. 3.] = [1. 2. 3.]
[1. 1. 1.] / [1. 2. 3.] = [1.         0.5        0.33333333]

If you want the dot product, you have to use the dot function:

print(x,".",y,"=",np.dot(x,y))

[1. 1. 1.] . [1. 2. 3.] = 6.0

Home Activity

Predict what will happen when you run the code below.

x = np.array([1,2.5,5])
y = np.array([4,6,7,9.0])
print(x+y)

Home Activity

Write a sentence with at least 50 characters describing why the behavior you observed makes sense. If it does not make sense, write out the question you will ask in class. Store your answer in the string my_answer.

my_answer = 'ehferbgrweavurebfvbidbv oisrevfbsrwiuevbgirwosgbvisbv . wefigwiuegfiuw'

# Removed autograder test. You may delete this cell.

Matrices work about the same way

silly_matrix = np.array([(1,2,3),(1,2,3),(1,2,3)])
print("The sum of\n",identity33,"\nand\n",
      silly_matrix,"\nis\n",identity33+silly_matrix)

The sum of
 [[1. 0. 0.]
 [0. 1. 0.]
 [0. 0. 1.]] 
and
 [[1 2 3]
 [1 2 3]
 [1 2 3]] 
is
 [[2. 2. 3.]
 [1. 3. 3.]
 [1. 2. 4.]]

Multiplication and division are element-wise:

identity33 * silly_matrix

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

identity33 / silly_matrix

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

The dot function will give you the matrix product:

print("The matrix product of\n",identity33,
      "\nand\n",silly_matrix,"\nis\n",np.dot(identity33,silly_matrix))

The matrix product of
 [[1. 0. 0.]
 [0. 1. 0.]
 [0. 0. 1.]] 
and
 [[1 2 3]
 [1 2 3]
 [1 2 3]] 
is
 [[1. 2. 3.]
 [1. 2. 3.]
 [1. 2. 3.]]

#matrix times a vector
y = np.array([1.0,2.0,3.0])
print(silly_matrix,"times", y, "is")
print(np.dot(silly_matrix,y))

[[1 2 3]
 [1 2 3]
 [1 2 3]] times [1. 2. 3.] is
[14. 14. 14.]

Writing np.dot can get tedious. You can also use the @ symbol for matrix multiplication:

print(silly_matrix @ y)

[14. 14. 14.]

1.10.5. Universal Functions#

NumPy also provides universal functions that operate on each element of an array. Common mathematical functions are defined in this way.

#recall we defined X as a linspace from 0 to 2pi
print(X)

# print an empty line
print("\n")

#taking the sin(X) should be one whole sine wave
print(np.sin(X))

[0.         0.10300304 0.20600608 0.30900911 0.41201215 0.51501519
 0.61801823 0.72102126 0.8240243  0.92702734 1.03003038 1.13303342
 1.23603645 1.33903949 1.44204253 1.54504557 1.64804861 1.75105164
 1.85405468 1.95705772 2.06006076 2.16306379 2.26606683 2.36906987
 2.47207291 2.57507595 2.67807898 2.78108202 2.88408506 2.9870881
 3.09009113 3.19309417 3.29609721 3.39910025 3.50210329 3.60510632
 3.70810936 3.8111124  3.91411544 4.01711848 4.12012151 4.22312455
 4.32612759 4.42913063 4.53213366 4.6351367  4.73813974 4.84114278
 4.94414582 5.04714885 5.15015189 5.25315493 5.35615797 5.459161
 5.56216404 5.66516708 5.76817012 5.87117316 5.97417619 6.07717923
 6.18018227 6.28318531]


[ 0.00000000e+00  1.02820997e-01  2.04552066e-01  3.04114832e-01
  4.00453906e-01  4.92548068e-01  5.79421098e-01  6.60152121e-01
  7.33885366e-01  7.99839245e-01  8.57314628e-01  9.05702263e-01
  9.44489229e-01  9.73264374e-01  9.91722674e-01  9.99668468e-01
  9.97017526e-01  9.83797952e-01  9.60149874e-01  9.26323968e-01
  8.82678798e-01  8.29677014e-01  7.67880446e-01  6.97944155e-01
  6.20609482e-01  5.36696194e-01  4.47093793e-01  3.52752087e-01
  2.54671120e-01  1.53890577e-01  5.14787548e-02 -5.14787548e-02
 -1.53890577e-01 -2.54671120e-01 -3.52752087e-01 -4.47093793e-01
 -5.36696194e-01 -6.20609482e-01 -6.97944155e-01 -7.67880446e-01
 -8.29677014e-01 -8.82678798e-01 -9.26323968e-01 -9.60149874e-01
 -9.83797952e-01 -9.97017526e-01 -9.99668468e-01 -9.91722674e-01
 -9.73264374e-01 -9.44489229e-01 -9.05702263e-01 -8.57314628e-01
 -7.99839245e-01 -7.33885366e-01 -6.60152121e-01 -5.79421098e-01
 -4.92548068e-01 -4.00453906e-01 -3.04114832e-01 -2.04552066e-01
 -1.02820997e-01 -2.44929360e-16]

import matplotlib.pyplot as plt

#the next line is only needed in iPython notebooks
%matplotlib inline

# Create a plot. We'll learn more about this below
plt.plot(X,np.sin(X));
plt.axis([0,2*np.pi,-1.05,1.05])
plt.show()

Home Activity

Copy the code from above to below. Then modify to plot cosine instead of sine.

# Add your solution here

1.10.6. Copying Arrays#

When you assign a new variable name to an existing array, it is the same as giving two names for the object.

a = np.array([1.0,2,3,4,5,6])
print(a)
b = a # this will make a and b different names for the same array
b[2] = 2.56 #change b at position 2, also changes a
print("The value of array a is",a)
print("The value of array b is",b)

[1. 2. 3. 4. 5. 6.]
The value of array a is [1.   2.   2.56 4.   5.   6.  ]
The value of array b is [1.   2.   2.56 4.   5.   6.  ]

Warning to MATLAB users: MATLAB would not change a with the assignment b[2] = 2.56. Python and MATLAB behave differently in this respect.

When you pass an array to a function, the function does not copy the array, it just assigns that array another name as in the previous example.

a1 = ["apples","oranges"]
b1 = a1
b1[1] = "grapes"
print(a1)
print(b1)

['apples', 'grapes']
['apples', 'grapes']

def devious_function(func_array):
    func_array[0] = -1.0e6 #changes the value of array passed in

a = np.array([1.0,2,3,4,5,6])
print("Before the function a =",a)
devious_function(a)
print("After the function a =",a)

Before the function a = [1. 2. 3. 4. 5. 6.]
After the function a = [-1.e+06  2.e+00  3.e+00  4.e+00  5.e+00  6.e+00]

def less_devious_function(func_float):
    func_float = 9.9 #changes local copy of func_float
a = 1
print("Before the function a =",a)
less_devious_function(a)
print("After the function a =",a)

Before the function a = 1
After the function a = 1

This is different than what we saw previously for passing floats, integers, booleans, and strings to functions.

Why would Python do this, besides to just confuse me!?! you ask.

Imagine a NumPy array has millions of elements. It is best for efficient memory usage if functions do not make multiple copies of millions of elements.

Summary:

Floats, integers, booleans, strings: You do NOT need to worry about .copy(). These are called primative types in Python.
Lists, arrays, dictionaries, and anything else. You do need to worry about .copy().

To make a real copy of the entire elements of an array you need to use the copy function on the array.

a = np.array([1.0,2,3,4,5,6])
print(a)
#this will make a and b different copies for the same array
b = a.copy()
b[2] = 2.56 #change b at position 2, will not change a
print("\nThe value of array a is",a)
print("\nThe value of array b is",b)

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

The value of array a is [1. 2. 3. 4. 5. 6.]

The value of array b is [1.   2.   2.56 4.   5.   6.  ]

1.10.7. Indexing, Slicing, and Iterating#

Oftentimes we’ll want to access more than one element of an array at a time. We can do this by slicing.

#bring these guys back
import numpy as np
a_vector = np.array([1,2,3,4])
a_matrix = np.array([(1,2,3),(4,5,6),(7,8,9)])
print("The vector",a_vector)
print("\nThe matrix\n",a_matrix)
print("\na_vector[:] =",a_vector[:]) #single colon gives everything
print("\na_vector[1:3] =",a_vector[1:3]) #print out position 1 to position 2 (same as for lists)
print("\nFor a matrix, we can slice in each dimension")
print("a_matrix[0,:] =",a_matrix[0,:]) #every column in row 0
print("a_matrix[0,1:3] =",a_matrix[0,1:3]) #columns 1 and 2 in row 0
print("a_matrix[:,2] =",a_matrix[:,2]) #every row in column 2

The vector [1 2 3 4]

The matrix
 [[1 2 3]
 [4 5 6]
 [7 8 9]]

a_vector[:] = [1 2 3 4]

a_vector[1:3] = [2 3]

For a matrix, we can slice in each dimension
a_matrix[0,:] = [1 2 3]
a_matrix[0,1:3] = [2 3]
a_matrix[:,2] = [3 6 9]

Iteration over a matrix will give you everything in a row:

a_matrix = np.array([(1,2,3),(4,5,6),(7,8,9)])
count = 0
for i in a_matrix:
    print("Row",count,"of a_matrix is",i)
    count += 1
    
    
for column in a_matrix.transpose():
    print("Column",count,"of a_matrix is",column)
    count += 1

Row 0 of a_matrix is [1 2 3]
Row 1 of a_matrix is [4 5 6]
Row 2 of a_matrix is [7 8 9]
Column 3 of a_matrix is [1 4 7]
Column 4 of a_matrix is [2 5 8]
Column 5 of a_matrix is [3 6 9]

To iterate over every element in the matrix you’ll need two for loops:

a_matrix = np.array([(1,2,3),(4,5,6),(7,8,9)])
row_count = 0
col_count = 0
for row in a_matrix:
    col_count = 0
    for col in row:
        print("Row",row_count,"Column",col_count,"of a_matrix is",col)
        col_count += 1
    row_count += 1

Row 0 Column 0 of a_matrix is 1
Row 0 Column 1 of a_matrix is 2
Row 0 Column 2 of a_matrix is 3
Row 1 Column 0 of a_matrix is 4
Row 1 Column 1 of a_matrix is 5
Row 1 Column 2 of a_matrix is 6
Row 2 Column 0 of a_matrix is 7
Row 2 Column 1 of a_matrix is 8
Row 2 Column 2 of a_matrix is 9

Home Activity

Using a single for loop, print the diagonal elements of matrix A. Store the diagonal elements in the list diagonal.

Tip: Let my_list be a Python list. Remember that my_list.append(1.0) will add 1.0 to the list my_list.

N = 3
A = np.random.rand(N,N)
print("A =\n",A)

# create any empty list
diagonal = []

#print the diagonal of the matrix using a single for loop
print("\nThe diagonal elements of A")

# Add your solution here

Warning: You need to programmically extract the diagonals. You cannot manually copy them. When the autograder runs, it will generate a new A matrix.

# Removed autograder test. You may delete this cell.

NumPy also has a built-in function to extract the diagonal of a matrix.

print("\nThe diagonal of A as a vector")
print(np.diag(A))

The diagonal of A as a vector
[0.80080604 0.50448948 0.00712939]

1.10.8. Complex Numbers#

Numpy can handle complex numbers, but you need to tell it that you have complex numbers by adding a parameter to the function that creates your array (or matrix, or vector, etc.).

cArray = np.zeros(5, dtype="complex")
print(cArray)

[0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j]

dtype tells NumPy that the array will hold complex numbers. It is not needed if it is obvious you are using complex numbers.

cArray2 = np.array([1+1j,-1])
print(cArray2)

[ 1.+1.j -1.+0.j]

You must be careful when you manipulate an array of floats into complex numbers, say by taking the square root

fArray = np.array([-1,1])
print("Wrong way gives", np.sqrt(fArray))

Wrong way gives [nan  1.]

C:\Users\ebrea\AppData\Local\Temp\ipykernel_22088\1555985613.py:2: RuntimeWarning: invalid value encountered in sqrt
  print("Wrong way gives", np.sqrt(fArray))

print("Right way gives", np.sqrt(fArray, dtype="complex"))

Right way gives [0.+1.j 1.+0.j]

1.10.9. Integrating Everything Together#

Class Activity

With a partner, predict the output of Python code.

# Add your solution here

Class Activity

With a partner, create a Python function that does the following:

Takes the size of a matrix, \(N\), as the input. Rounds to an integer.
If \(N < 1\), prints a warning message and returns nothing. Otherwise, continues.
Generates a matrix of size \(N \times N\) filled with random floating point numbers from -5 to 5
Print the matrix to the screen
Returns a vector that is the sum of the columns of the matrix

def np_activity(N):
    '''Generates a random N x N matrix, prints to screen,
    and returns a vector that is the sum of the columns.
    
    Arguments:
        YOU NEED TO FILL THIS IN
    
    Returns:
        YOU NEED TO FILL THIS IN
    '''

Discussion Question: what would we need to change to access the matrix we created outside our function?

Data and Computing for Chemical Engineers

Linear Algebra with Numpy and Scipy

Contents

1.10. Linear Algebra with Numpy and Scipy#

1.10.1. Learning Objectives#

1.10.2. Getting Started with NumPy Arrays#

1.10.3. Creating Arrays in Neat Ways#

1.10.4. Operations on Arrays#

1.10.5. Universal Functions#

1.10.6. Copying Arrays#

1.10.7. Indexing, Slicing, and Iterating#

1.10.8. Complex Numbers#

1.10.9. Integrating Everything Together#