Python is a high-level open-source language. But the Python world is inhabited by many packages or libraries that provide useful things like array operations, plotting functions, and much more. We can import libraries of functions to expand the capabilities of Python in our programs.
OK! We’ll start by importing a few libraries to help us out. First: our favorite library is
NumPy, providing a bunch of useful array operations (similar to MATLAB). We will use it a lot! The second library we need is
Matplotlib, a 2D plotting library which we will use to plot our results.
import numpy as np from matplotlib import pyplot as plt
We are importing one library named
numpy (that we will call
np)and we are importing a module called
pyplot (we refer to this library as
plt) of a big library called
matplotlib. To use a function belonging to one of these libraries, we have to tell Python where to look for it. For that, each function name is written following the library name, with a dot in between. So if we want to use the NumPy function
linspace(), which creates an array with equally spaced numbers between a start and end, we call it by writing:
import numpy as np from matplotlib import pyplot as plt myarray = np.linspace(0,5,10) print(myarray)
Python doesn’t require explicitly declared variable types like Fortran and other languages.
a =5 #a is an integer b= 'Five' #b is a string c = 5.0 #c is a floating point number print(type(a)) print(type(b)) print(type(c))
Note: if you divide an integer by an integer that yields a remainder, the result will be converted to a float.
Whitespace in Python
Python uses indents and whitespace to group statements together. To write a short loop in Fortran, you might use:
b= "Five" do i = 1, 5 print(b) end do
Python does not use
end do blocks like Fortran, so the same program as above is written in Python as follows:
b= 'Five' #b is a string for i in range(5): print(b)
If you have nested for-loops, there is a further indent for the inner loop.
for i in range(3): for j in range(3): print(i, j) print("This statement is within the i-loop, but not the j-loop")
NumPy, you can look at portions of arrays in the same way as in Matlab, with a few extra tricks thrown in. Let’s take an array of values from 1 to 5.
vec = np.array([1, 2, 3, 4, 5]) print(vec)
Python uses a zero-based index, so let’s look at the first and last element in the array
print(vec) print(vec, vec)
it has this output:
[1 2 3 4 5] 1 5
Arrays can also be ‘sliced‘, grabbing a range of values. Let’s look at the first three elements:
vec = np.array([1, 2, 3, 4, 5]) print(vec) print(vec, vec) print(vec[0:3])
Note here, the slice is inclusive on the front end and exclusive on the back, so the above command gives us the values of
vec, but not
Assigning Array Variables
One of the strange little quirks/features in Python that often confuses people comes up when assigning and comparing arrays of values. Here is a quick example. Let’s start by defining a 1-D array called
vecx = np.linspace(1,5,5) print(vecx)
linspace command return evenly spaced numbers over a specified interval. In this case, for
vecx we have:
[1. 2. 3. 4. 5.]
OK, so we have an array
vecx, with the values 1 through 5. I want to make a copy of that array, called
vecy, so I’ll try the following:
vecy = vecx print(vecy)
and we have:
[1. 2. 3. 4. 5.]
vecx has the values 1 through 5 and now so does
vecy. Now that I have a backup of
vecx, I can change its values without worrying about losing data (or so I may think!).
vecx=10 print(vecx) print(vecy)
And that’s how things go wrong! When you use a statement like
vecx=vecy, rather than copying all the values of vecx
into a new array called
vecy, Python just creates an alias (or a pointer) called
vecy and tells it to route us to
vecx. So if we change a value in
vecy will reflect that change (technically, this is called assignment by reference).
If you want to make a true copy of the array, you have to tell Python to
copy every element of
vecx into a new array. Let’s call it
vecx = np.linspace(1,5,5) print(vecx) vecy = vecx vecz=vecx.copy() vecx=10 print(vecx) print(vecy) print(vecz)
A matrix is a 2D array, where each element in the array has 2 indices. For example,
[[1, 2], [3, 4]] is a matrix, and the index of
1 is (0,0).
We can prove this using Python and
import numpy as np A = [[1, 2], [3, 4]] print(np.array(A)[0,0])
The simple form of matrix multiplication is called scalar multiplication, multiplying a scalar by a matrix.
Scalar multiplication is generally easy. Each value in the input matrix is multiplied by the scalar, and the output has the same shape as the input matrix.
import numpy as np b=7 A = [[1, 2], [3, 4]] C=np.dot(b,A) print(np.array(C))
There are a few things to keep in mind.
- Order matters now. \(A \times B \neq B \times A\);
- Matrices can be multiplied if the number of columns in the 1st equals the number of rows in the 2nd;
- Multiplication is the dot product of rows and columns. Rows of the 1st matrix with columns of the 2nd.
Let’s replicate the result in Python.
import numpy as np A = [[1, 2], [3, 4]] B = [[5, 6], [7, 8]] C = np.dot(A,B) print(np.array(C))
The inner product takes two vectors of equal size and returns a single number (scalar). This is calculated by multiplying the corresponding elements in each vector and adding up all of those products. In numpy, vectors are defined as one-dimensional numpy arrays.
To get the inner product, we can use either
np.dot(). Both give the same results. The inputs for these functions are two vectors and they should be the same size.
import numpy as np # Vectors as 1D numpy arrays a = np.array([1, 2, 3]) b = np.array([4, 5, 6]) print("a= ", a) print("b= ", b) print("\ninner:", np.inner(a, b)) print("dot:", np.dot(a, b))
The transpose of a matrix is found by switching its rows with its columns. We can use np.transpose() function or NumPy ndarray.transpose() method or ndarray.T (a special method which does not require parentheses) to get the transpose. All give the same output.
import numpy as np a = np.array([[1, 2], [3, 4], [5, 6]]) print("a = ") print(a) print("\nWith np.transpose(a) function") print(np.transpose(a)) print("\nWith ndarray.transpose() method") print(a.transpose()) print("\nWith ndarray.T short form") print(a.T)
The trace is the sum of diagonal elements in a square matrix. There are two methods to calculate the trace. We can simply use the
trace() method of an ndarray object or get the diagonal elements first and then get the sum.
import numpy as np a = np.array([[2, 2, 1], [1, 3, 1], [1, 2, 2]]) print("a = ") print(a) print("\nTrace:", a.trace()) print("Trace:", sum(a.diagonal()))
The rank of a matrix is the dimensions of the vector space spanned (generated) by its columns or rows. In other words, it can be defined as the maximum number of linearly independent column vectors or row vectors.
The rank of a matrix can be found using the
matrix_rank() function which comes from the
numpy linalg package.
import numpy as np a = np.arange(1, 10) a.shape = (3, 3) print("a = ") print(a) rank = np.linalg.matrix_rank(a) print("\nRank:", rank)
The determinant of a square matrix can be calculated
det() function which also comes from the numpy linalg package. If the determinant is 0, that matrix is not invertible. It is known as a singular matrix in algebra terms.
import numpy as np a = np.array([[2, 2, 1], [1, 3, 1], [1, 2, 2]]) print("a = ") print(a) det = np.linalg.det(a) print("\nDeterminant:", np.round(det))
The true inverse of a square matrix can be found using the
inv() function of the numpy linalg package. If the determinant of a square matrix is not 0, it has a true inverse.
import numpy as np a = np.array([[2, 2, 1], [1, 3, 1], [1, 2, 2]]) print("a = ") print(a) det = np.linalg.det(a) print("\nDeterminant:", np.round(det)) inv = np.linalg.inv(a) print("\nInverse of a = ") print(inv)
Flatten is a simple method to transform a matrix into a one-dimensional numpy array. For this, we can use the flatten() method of an ndarray object.
import numpy as np a = np.arange(1, 10) a.shape = (3, 3) print("a = ") print(a) print("\nAfter flattening") print("------------------") print(a.flatten())
Eigenvalues and eigenvectors
Let A be an n x n matrix. A scalar \(\lambda\) is called an eigenvalue of A if there is a non-zero vector x satisfying the following equation.
\(A x=\lambda x\)
The vector x is called the eigenvector of A corresponding to \(\lambda\).
In numpy, both eigenvalues and eigenvectors can be calculated simultaneously using the
import numpy as np a = np.array([[2, 2, 1], [1, 3, 1], [1, 2, 2]]) print("a = ") print(a) w, v = np.linalg.eig(a) print("\nEigenvalues:") print(w) print("\nEigenvectors:") print(v)
The sum of eigenvalues (1+5+1=7) is equal to the trace (2+3+2=7) of the same matrix! The product of the eigenvalues (1x5x1=5) is equal to the determinant (5) of the same matrix!