This lecture¶
- NumPy arrays ‐ vectors
- NumPy arrays ‐ matrices
- Linear transformations ‐ rotation matrix in the plane
Before we start: Imports¶
We will work with NumPy (for arrays and computations) and Matplotlib (for plotting):
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline
The history saving thread hit an unexpected error (OperationalError('attempt to write a readonly database')).History will not be written to the database.
Every function/object from these packages starts with np. or plt., respectively
Motivation¶
Many technical applications ultimately lead to problems of the following form. Given $A \in \mathbb{R}^{n \times n}$ and $b \in \mathbb{R}^n$, $n>0$, find $x \in \mathbb{R}^n$ such that
$$ A \, x = b. $$- Weather forecast and renewable energies
- Numerical simulation in general (droplet simulation)
- Banking and financial applications
- Data Science
- and many more
NumPy arrays ‐ vectors¶
Lists are almost like vectors but operations on lists are not linear algebra operations.
l1 = [1, 2]
l2 = [3, 4]
print(f"{l1 + l2 = }")
print(f"{3 * l1 = }")
l1 + l2 = [1, 2, 3, 4] 3 * l1 = [1, 2, 1, 2, 1, 2]
v = np.array([3, -4]) # a vector in the plane
print(f"{v = }")
print("Elementwise operations:")
print(f"{2*v = }")
print(f"{v*v = }")
print(f"{v/2 = }")
print(f"{v/v = }\n")
print("Euclidian norm ||v||:")
print(f"using {np.linalg.norm(v) = }")
print(f"using {np.sqrt(np.dot(v,v)) = }\n")
print("Scalar product <v, v>:")
print(f"using {np.dot(v, v) = }")
print(f"using {v @ v = }\n")
v = array([ 3, -4]) Elementwise operations: 2*v = array([ 6, -8]) v*v = array([ 9, 16]) v/2 = array([ 1.5, -2. ]) v/v = array([1., 1.]) Euclidian norm ||v||: using np.linalg.norm(v) = 5.0 using np.sqrt(np.dot(v,v)) = 5.0 Scalar product <v, v>: using np.dot(v, v) = 25 using v @ v = 25
Vectors are similar to lists¶
l1 = [1, 2, 3]
v1 = np.array(l1)
print(f"{l1 = }")
print(f"{v1 = }")
print(type(l1))
print(type(v1))
l1 = [1, 2, 3] v1 = array([1, 2, 3]) <class 'list'> <class 'numpy.ndarray'>
Vectors are similar to lists (cont)¶
- Access vectors via their indices
v = np.array([1, 2, 3])
print(v[0])
1
- The length of a vector is still obtained by the function
len
print(len(v))
3
Vectors are similar to lists (cont)¶
- Parts of vectors using slices
v = np.array([1, 2, 3])
print(v[:2])
[1 2]
- Replace parts of vectors using slices
v[:2] = [10, 20]
print(v)
[10 20 3]
Vectors are not lists!¶
Operations are not the same:
- Operations
+and*are different - More operations are defined:
-,/ - Numpy functions act elementwise on vectors:
sin,exp,sqrt, etc. - Scalar product with
dotor@ - Vectors have a fixed size: no
appendmethod!
Vectors are not lists! (cont)¶
Two names are bound to the same list:
l1 = [1, 2]
l2 = l1
l1[0] = 10
print(l2)
[10, 2]
To make an actual copy, you can use the copy method, or make a slice.
l1 = [1, 2]
l2 = l1
l3 = l1.copy()
l4 = l1[:] # a slice of the entire list
l1[0] = 10
print(f"{l2 = }")
print(f"{l3 = }")
print(f"{l4 = }")
l2 = [10, 2] l3 = [1, 2] l4 = [1, 2]
Vectors are not lists! (cont)¶
- Vector slices, on the other hand, are views:
l1 = [1, 2, 3]
v1 = np.array([1, 2, 3])
l2 = l1[:2]
l2[0] = 10
print(f"now {l1 = }")
v3 = v1.copy()[:2]
v3[0] = 10
print(f"now {v1 = }")
v2 = v1[:2]
v2[0] = 10
print(f"now {v1 = }")
v3 = v1.copy()[:2]
v3[0] = 10
now l1 = [1, 2, 3] now v1 = array([1, 2, 3]) now v1 = array([10, 2, 3])
Vectors are not lists! (cont)¶
- Only one type throughout the entire vector
v1 = np.array([1., 2., 3.]) # don't forget the dots!
v2 = np.array([2, 0, 1])
print(v2)
print(v1*v2) # what about v1/v2
print(3*v1 + 2*v2)
print(v2)
print(type(v1[0]))
print(type(v2[0]))
[2 0 1] [2. 0. 3.] [ 7. 6. 11.] [2 0 1] <class 'numpy.float64'> <class 'numpy.int64'>
More examples¶
v1 = np.array([1., 2., 3.])
v2 = np.array([2, 0, 1.])
print(type(v2[0]))
# s = v1[0]*v2[0] + v1[1]*v2[1] + v1[2] * v2[2]
# print(f"scalar product <v1, v2> = {v1 @ v2 = }")
# access
v1[0] = 10
# print(f"{v1 = }")
# slices
v1[:2] = [0, 1]
# print(f"{v1 = }")
# v1[0:2] = [1, 2, 3] # error
l1 = [0, 1, 2]
l1[:2] = [1, 2, 3]
# print(f"{l1 = }")
<class 'numpy.float64'>
Creating vectors ‐ linspace¶
The linspace function is convenient for creating equally spaced arrays.
xs = np.linspace(0, 10, 200) # 200 points between 0 and 10
print(type(xs[0]))
print(xs[0]) # the first point is 0
print(xs[-1]) # the last point is 10
<class 'numpy.float64'> 0.0 10.0
plt.plot(xs, np.sin(xs))
plt.xlabel("x")
plt.ylabel("f(x)")
Text(0, 0.5, 'f(x)')
Creating vectors ‐ zeros, ones, full¶
Some handy functions to quickly create vectors:
v1 = np.zeros(5)
print("v1 =", v1)
v2 = np.ones(3)
print("v2 =", v2)
v3 = np.full(4, 3.14) # fill it with the value 3.14
print("v3 =", v3)
v4 = np.full(5, 0.0)
print("v4 =", v4)
v1 = [0. 0. 0. 0. 0.] v2 = [1. 1. 1.] v3 = [3.14 3.14 3.14 3.14] v4 = [0. 0. 0. 0. 0.]
Creating vectors using dtype¶
v1 = np.zeros(5, dtype = complex)
print("v1 =", v1)
v2 = np.ones(3, dtype = int)
print("v2 =", v2)
v3 = np.full(4, "hello")
print("v3 =", v3)
print(type(v3[0]))
v1 = [0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j] v2 = [1 1 1] v3 = ['hello' 'hello' 'hello' 'hello'] <class 'numpy.str_'>
Creating random vectors¶
Fill a vector with random numbers in the interval [0, 1).
v1 = np.random.random(5)
print(v1)
[0.26190383 0.85589526 0.83860531 0.97295582 0.01607221]
We can also "shuffle" the values in an array.
v2 = np.array([10, 20, 30, 40, 50, 60, 80])
np.random.shuffle(v2)
print(v2)
[10 60 20 40 30 80 50]
Concatenating vectors¶
Since the + operation is redefined, we need other means to concatenate vectors.
The command hstack([v1, v2, ..., vn]) concatenates the vectors
$v_1, v_2, \ldots, v_n$.
v1 = [1, 2, 3]
v2 = [4, 5]
v3 = [6, 7, 8, 9]
L = v1 + v2 + v3
print(L)
v = np.hstack([v1, v2, v3])
# w = np.vstack([v1, v2, v3])
print(v)
[1, 2, 3, 4, 5, 6, 7, 8, 9] [1 2 3 4 5 6 7 8 9]
Example¶
Compute the distance between the unit vectors in $\mathbb{R}^2$, i.e. $a = (1,0)^T$ and $b = (0,1)^T$. Compute the distance between $a$ and $b$. 5 min, let's go!
a = np.array([1, 0])
b = np.array([0, 1])
np.linalg.norm(a - b, 2)
1.4142135623730951
Functions and arrays¶
As long as +, -, *, / and standard mathematical numpy functions are used, a function value can be calculated elementwise.
def my_func(x, A, omega):
return A * np.sin(omega * x) / x
x = np.linspace(-100, 100, 200) # using an even number of values to avoid zero
A = [2 for i in range(200)] # full(100, 2)
y = my_func(x, A, 0.5)
plt.plot(x, y)
plt.xlabel("x")
plt.ylabel("f(x)");
Vectorized functions¶
Note that not all functions may be applied to vectors. For instance this one:
def const(x):
return 1
x = np.linspace(-100, 100, 50)
y = const(x)
#plt.plot(x, y)
# what is the value of y?
print(y)
1
We will see later how to automatically vectorize a function so that it works on vectors.
# An identity matrix
id_matrix = np.array([[1., 0.], [0., 1.]])
# Python allows this:
id_matrix = np.array([[1., 0],
[0., 1.]])
# Which is more readable?
print(id_matrix)
[[1. 0.] [0. 1.]]
Accessing matrix entries¶
Matrix elements are accessed with two indices.
M = np.array([[1., 2.], [3., 4.]])
print(M[0, 0]) # first row, first column
print(M[-1, 0]) # last row, first column
# Don't do this: (performance)
print(M[0][0])
1.0 3.0 1.0
Compare this to accessing list elements.
L = [[1., 2.], [3., 4.]]
print(L[0][0])
print(L[-1][0])
1.0 3.0
Creating matrices ‐ eye, zeros, full¶
eye(n)is the identity matrix of size $n$.zeros([n, m])fills an $n\times m$ matrix with zeroesfull([n, m], value)fills an $n\times m$ matrix withvalue
id_matrix = np.eye(2)
M0 = np.zeros((2, 3)) # using a list [2, 3]
M1 = np.full((2, 3), 7.) # using a tuple (2, 3)
print(id_matrix)
print(M0)
print(M1)
[[1. 0.] [0. 1.]] [[0. 0. 0.] [0. 0. 0.]] [[7. 7. 7.] [7. 7. 7.]]
Creating random matrices¶
M1 = np.random.random((3, 4))
print(M1)
M2 = np.random.randint(1, 7, (5,))
print(M2)
[[0.23206038 0.96568178 0.82859427 0.07275405] [0.74440674 0.77756663 0.96288204 0.12766979] [0.39198518 0.2577695 0.66056364 0.41648916]] [4 2 2 4 6]
Shape¶
The shape of a matrix is the tuple of its dimensions. The shape of an $n\times m$ matrix is (n,m). It is given by the attribute shape:
M = np.eye(3)
print(M.shape)
V = np.array([1., 2., 1., 4.])
print(V.shape) # tuple with one element
(3, 3) (4,)
Tip¶
np.zeros(A.shape) returns a matrix with the same shape as A containing only zeros.
Transpose¶
The transpose of a matrix $A$ is a matrix $B$ such that
$$B_{ij} = A_{ji}.$$By transposing a matrix you switch the two shape elements.
A = np.array([[1, 2, 3],
[4, 5, 6]])
print("A.shape =", A.shape)
B = A.T # A transpose
print("B.shape =", B.shape, '\n')
B[0,0] = 11
print("A =")
print(A, '\n')
print("B =")
print(B)
A.shape = (2, 3) B.shape = (3, 2) A = [[11 2 3] [ 4 5 6]] B = [[11 4] [ 2 5] [ 3 6]]
1D vectors and 2D row/column matrices¶
(here D stands for tensor dimension: the length of shape)
v = np.array([1, 2, 3])
Mcol = np.array([[1],
[2],
[3]])
Mrow = np.array([[1, 2, 3]])
print("v.shape = ", v.shape)
print("Mcol.shape = ", Mcol.shape)
print("Mrow.shape = ", Mrow.shape)
v.shape = (3,) Mcol.shape = (3, 1) Mrow.shape = (1, 3)
Transposing a 1D vector doesn't change it
Matrix multiplication¶
We can use the dot function for matrix multiplication. We can also use the @ operator.
If $A$ is a $n\times m$ matrix and $B$ a $m \times p$ matrix, then $AB$ is a $n\times p$ matrix.
A = np.array([[1, 0, 0], # n = 2, m = 3
[2, -1, 3]])
B = np.array([[2, 1, 0, 4], # m = 3, p = 4
[1, 1, -1, 0],
[0, 1, 2, 2]])
AB = np.dot(A, B)
print(AB)
# same as
AB = A @ B
print(AB)
# BA = np.dot(B, A) # error
[[ 2 1 0 4] [ 3 4 7 14]] [[ 2 1 0 4] [ 3 4 7 14]]
Matrix vector multiplication¶
A matrix can be multipled with a vector using dot or the @ operator.
The multiplication is matrix multiplication, as if the vector was a column or row matrix. The result is a vector.
A = np.array([[1, 0],
[4, -1]])
v = np.array([1, 2])
vA = np.dot(v, A)
vA = v @ A
Av = np.dot(A, v)
print(vA)
print(Av)
print("vA.shape =", vA.shape, " Av.shape = ", Av.shape)
[ 9 -2] [1 2] vA.shape = (2,) Av.shape = (2,)
Dot product and @ operator¶
| vector vector | $s = \sum_i x_i y_i$ | s = x @ y or s = dot(x,y) | scalar product |
|:---|:---|:---|:---|
| matrix matrix | $C_{ij} = \sum_k A_{ik}B_{kj}$ | C = A @ B or C = dot(A, B) | matrix multiplication |
| matrix vector | $y_i = \sum_j A_{ij}x_j$ | y = A @ x or y = dot(A, x) | matrix multiplication, treating vector as column matrix |
| vector matrix | $y_j = \sum_i x_i A_{ij}$ | y = x @ A or y = dot(x, A) | matrix multiplication, treating vector as row matrix |
Linear transformations ‐ rotation matrix in the plane¶
Definition¶
A linear transformation in the plane is a function $f:\mathbb{R}^2 \mapsto \mathbb{R}^2$ such that
$f(a\mathbf{u} + b\mathbf{v}) = af (\mathbf{u}) + bf(\mathbf{v}) \text{ for all }a, b\in\mathbb{R} \text{ and } \mathbf{u}, \mathbf{v} \in \mathbb{R}^2$
Matrices and linear transformations¶
Let $x= (x_1, x_2) \in \mathbb{R}^2$ and $y= (y_1, y_2) \in \mathbb{R}^2$. Let $y = f(x)$ be the function defined by
$$ \begin{cases} y_1 = ax_1 + bx_2 \\ y_2 = cx_1 + dx_2 \end{cases} $$for some numbers $a, b, c, d$. Then matrix multiplication can be used to describe $f$. Let
$$A = \left( \begin{array}{rr} a & b \\ c & d \end{array} \right), x = \left( \begin{array}{rr} x_1 \\ x_2 \end{array} \right), y = \left( \begin{array}{rr} y_1 \\ y_2 \end{array} \right). \text{ Then } y = Ax, \text{ i.e. } \left( \begin{array}{rr} y_1 \\ y_2 \end{array} \right) = \left( \begin{array}{rr} a & b \\ c & d \end{array} \right) \left( \begin{array}{rr} x_1 \\ x_2 \end{array} \right) $$$A$ is a transformation matrix of $f$, $y = f(x) = Ax$
Finding the transformation matrix of a rotation in the plane¶
Let $\mathbf{e_1}, \mathbf{e_2}$ be two basis vectors and let $A$ be a transformation matrix that rotates a vector $\mathbf{u} = x\mathbf{e_1} + y\mathbf{e_2}$ by an angle $\alpha $. To find $A$ it suffices to consider the basis vectors.
Now $ A\mathbf{u} = xA\mathbf{e_1} + yA\mathbf{e_1} = x\mathbf{e_1}' + y\mathbf{e_2}'$
Rotate a point¶
alpha = np.pi / 2
A = np.array([[np.cos(alpha), -np.sin(alpha)],
[np.sin(alpha), np.cos(alpha)]])
x = 3
y = 4
P = np.array([x, y]) # a vector
AP = A @ P
plt.plot(x, y, 'bo')
plt.xlim(-5, 5)
plt.ylim(-5, 5)
plt.grid()
plt.plot(*AP, 'ro') # unpack to get x and y
[<matplotlib.lines.Line2D at 0x13491c650>]
Rotate many points¶
5min, let's go!
alpha = np.pi / 4
A = np.array([[np.cos(alpha), -np.sin(alpha)],
[np.sin(alpha), np.cos(alpha)]])
x = np.linspace(-10, 10, 10)
y = 3 * x
plt.plot(x, y)
plt.xlabel("x")
plt.ylabel("y")
M = np.array([x, y])
rot_M = A @ M
rot_x, rot_y = rot_M
plt.plot(rot_x, rot_y);
