Linear Algebra in NumPy

Unit 9, Lecture 2

Numerical Methods and Statistics

Prof. Andrew White, March 30, 2020



In [1]:

    
import random
import numpy as np
import matplotlib.pyplot as plt
from math import sqrt, pi, erf
import scipy.stats
import numpy.linalg

Working with Matrices in Numpy

We saw earlier in the class how to create numpy matrices. Let's review that and learn about explicit initialization

Explicit Initialization

You can explitily set the values in your matrix by first creating a list and then converting it into a numpy array



In [2]:

    
matrix = [ [4,3], [6, 2] ]
print('As Python list:')
print(matrix)

np_matrix = np.array(matrix)

print('The shape of the array:', np.shape(np_matrix))
print('The numpy matrix/array:')
print(np_matrix)









    



As Python list:
[[4, 3], [6, 2]]
The shape of the array: (2, 2)
The numpy matrix/array:
[[4 3]
 [6 2]]

You can use multiple lines in python to specify your list. This can make the formatting cleaner



In [3]:

    
np_matrix_2 = np.array([
                        [ 4, 3],   
                        [ 1, 2],  
                        [-1, 4], 
                        [ 4, 2]   
                        ])
print(np_matrix_2)









    



[[ 4  3]
 [ 1  2]
 [-1  4]
 [ 4  2]]

Create and Set

You can also create an array and then set the elements.



In [4]:

    
np_matrix_3 = np.zeros( (2, 10) )

print(np_matrix_3)









    



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



In [5]:

    
np_matrix_3[:, 1] = 2
print(np_matrix_3)









    



[[0. 2. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 2. 0. 0. 0. 0. 0. 0. 0. 0.]]



In [6]:

    
np_matrix_3[0, :] = -1
print(np_matrix_3)









    



[[-1. -1. -1. -1. -1. -1. -1. -1. -1. -1.]
 [ 0.  2.  0.  0.  0.  0.  0.  0.  0.  0.]]



In [7]:

    
np_matrix_3[1, 6] = 43
print(np_matrix_3)









    



[[-1. -1. -1. -1. -1. -1. -1. -1. -1. -1.]
 [ 0.  2.  0.  0.  0.  0. 43.  0.  0.  0.]]



In [8]:

    
rows, columns = np.shape(np_matrix_3) #get the number of rows and columns
for i in range(columns): #Do a for loop over columns
    np_matrix_3[1, i] = i ** 2 #Set the value of the 2nd row, ith column to be i^2
    
print(np_matrix_3)









    



[[-1. -1. -1. -1. -1. -1. -1. -1. -1. -1.]
 [ 0.  1.  4.  9. 16. 25. 36. 49. 64. 81.]]

Linear Algebra

The linear algebra routines for python are in the numpy.linalg library. See here

Matrix Multiplication

Matrix multiplication is done with the dot method. Let's compare that with *



In [9]:

    
np_matrix_1 = np.random.random( (2, 4) ) #create a random 2 x 4 array
np_matrix_2 = np.random.random( (4, 1) ) #create a random 4 x 1 array

print(np_matrix_1.dot(np_matrix_2))









    



[[2.02700017]
 [1.20440067]]

So, dot correctly gives us a 2x1 matrix as expected for the two shapes

Using the special @ character:



In [10]:

    
print(np_matrix_1 @ np_matrix_2)









    



[[2.02700017]
 [1.20440067]]

The element-by-element multiplication, *, doesn't work on different sized arrays.



In [11]:

    
print(np_matrix_1 * np_matrix_2)









    



---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-11-f16136fbd191> in <module>
----> 1 print(np_matrix_1 * np_matrix_2)

ValueError: operands could not be broadcast together with shapes (2,4) (4,1)

Method vs Function

Instead of using dot as a method (it comes after a .), you can use the dot function as well. Let's see an example:



In [12]:

    
print(np_matrix_1.dot(np_matrix_2))

print(np.dot(np_matrix_1, np_matrix_2))









    



[[2.02700017]
 [1.20440067]]
[[2.02700017]
 [1.20440067]]

Matrix Rank

The rank of a matrix can be found with singular value decomposition. In numpy, we can do this simply with a call to linalg.matrix_rank



In [13]:

    
import numpy.linalg as linalg

matrix =  [ [1, 0], [0, 0] ]
np_matrix = np.array(matrix)

print(linalg.matrix_rank(np_matrix))

Matrix Inverse

The inverse of a matrix can be found using the linalg.inverse command. Consider the following system of equations:

$$\begin{array}{lr} 3 x + 2 y + z & = 5\\ 2 x - y & = 4 \\ x + y - 2z & = 12 \\ \end{array}$$

We can encode it as a matrix equation:

$$\left[\begin{array}{lcr} 3 & 2 & 1\\ 2 & -1 & 0\\ 1 & 1 & -2\\ \end{array}\right] \left[\begin{array}{l} x\\ y\\ z\\ \end{array}\right] = \left[\begin{array}{l} 5\\ 4\\ 12\\ \end{array}\right]$$$$\mathbf{A}\mathbf{x} = \mathbf{b}$$$$\mathbf{A}^{-1}\mathbf{b} = \mathbf{x}$$



In [14]:

    
#Enter the data as lists
a_matrix = [[3, 2, 1],
            [2,-1,0],
            [1,1,-2]]
b_matrix = [5, 4, 12]

#convert them to numpy arrays/matrices
np_a_matrix = np.array(a_matrix)
np_b_matrix = np.array(b_matrix).transpose()

#Solve the problem
np_a_inv = linalg.inv(np_a_matrix)
np_x_matrix = np_a_inv @ np_b_matrix

#print the solution
print(np_x_matrix)

#check to make sure the answer works
print(np_a_matrix @ np_x_matrix)









    



[ 2.47058824  0.94117647 -4.29411765]
[ 5.  4. 12.]

Computation cost for inverse

Computing a matrix inverse can be VERY expensive for large matrices. Do not exceed about 500 x 500 matrices

Eigenvectors/Eigenvalues

Before trying to understand what an eigenvector is, let's try to understand their analogue, a stationary point.

A stationary point of a function $f(x)$ is an $x$ such that:

$$x = f(x)$$

Consider this function:

$$f(x) = x - \frac{x^2 - 612}{2x}$$

If we found a stationary point, that would be mean that

$$x = x - \frac{x^2 - 612}{2x} $$

$$ x^2 = 612 $$

More generally, you can find a square root of $A$ by finding a stationary point to:

$$f(x) = x - \frac{x^2 - A}{2x} $$

In this case, you can find the stationary point by just doing $x_{i+1} = f(x_i)$ until you are stationary



In [ ]:

    
x = 1
for i in range(10):
    x = x - (x**2 - 612) / (2 * x)
    print(i, x)

Eigenvectors/Eigenvalues

Matrices are analogues of functions. They take in a vector and return a vector. $$\mathbf{A}\mathbf{x} = \mathbf{y}$$

Just like stationary points, there is sometimes a special vector which has this property:

$$\mathbf{A}\mathbf{x} = \mathbf{x}$$

Such a vector is called an eigenvector. It turns out such vectors are rarely always exists. If we instead allow a scalar, we can find a whole bunch like this:

$$\mathbf{A}\mathbf{v} =\lambda\mathbf{v}$$

These are like the stationary points above, except we are getting back our input times a constant. That means it's a particular direction that is unchanged, not the value.

Finding Eigenvectors/Eigenvalues

Eigenvalues/eigenvectors can be found easily as well in python, including for complex numbers and sparse matrices. The command linalg.eigh will return only the real eigenvalues/eigenvectors. That assumes your matrix is Hermitian, meaning it is symmetric (if your matrix is real numbers). Use eig to get general possibly complex eigenvalues Here's an easy example:

Let's consider this matrix:

$$ A = \left[\begin{array}{lr} 3 & 1\\ 1 & 3\\ \end{array}\right]$$

Imagine it as a geometry operator. It takes in a 2D vector and morphs it into another 2D vector.

$$\vec{x} = [1, 0]$$$$A \,\vec{x}^T = [3, 1]^T$$

Now is there a particular direction where $\mathbf{A}$ cannot affect it?



In [ ]:

    
A = np.array([[3,1], [1,3]])

e_values, e_vectors = np.linalg.eig(A)

print(e_vectors)
print(e_values)

So that means $v_1 = [0.7, 0.7]$ and $v_2 = [-0.7, 0.7]$. Let's find out:



In [ ]:

    
v1 = e_vectors[:,0]
v2 = e_vectors[:,1]

A @ v1

Yes, that is the same direction! And notice that it's 4 times as much as the input vector, which is what the eigenvalue is telling us.

A random matrix will almost never be Hermitian, so look out for complex numbers. In engineering, your matrices commonly be Hermitian.



In [ ]:

    
A = np.random.normal(size=(3,3))
e_values, e_vectors = linalg.eig(A)
print(e_values)
print(e_vectors)

Notice that there are compelx eigenvalues, so eigh was not correct to use