Diary

Libraries

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)

Variables

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))

we have:

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")

Slicing Arrays

In 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 vec

print(vec)
print(vec[0], vec[4])

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[0], vec[4])
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[0], vec[1] and vec[2], but not vec[3].

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:

vecx = np.linspace(1,5,5)
print(vecx)

The 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.]

Great. So 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[2]=10
print(vecx)
print(vecy)

we have:

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 vecx then 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 vecz:

vecx = np.linspace(1,5,5)
print(vecx)

vecy = vecx
vecz=vecx.copy()

vecx[2]=10

print(vecx)
print(vecy)
print(vecz)

Matrix Multiplication

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 Numpy.

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.

In Python:

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.

1. Order matters now. \(A \times B \neq B \times A\);
2. Matrices can be multiplied if the number of columns in the 1st equals the number of rows in the 2nd;
3. 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))

Inner product

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.inner() or 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))

Transpose

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)

Trace

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()))

Rank

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)

Determinant

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))

True inverse

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

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 eig() function.

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!