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

We saw earlier in the class how to create numpy matrices. Let's review that and learn about 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)
```

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

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

In [4]:

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

In [5]:

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

In [6]:

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

In [7]:

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

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

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

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

The element-by-element multiplication, `*`

, doesn't work on different sized arrays.

In [11]:

```
print(np_matrix_1 * np_matrix_2)
```

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

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

The inverse of a matrix can be found using the `linalg.inverse`

command. Consider the following system of equations:

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

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

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

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:

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.

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