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

1


## 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}$$

or

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