Array definition

• Arrays have a fixed length
• Array elements have the same type
• We will use arrays from the Numpy library
• Optimized for numerical calculations

# Load the Numpy library
import numpy as np

a = np.array([1, 2, 4, 6, 5, 9])
print(a)

[1 2 4 6 5 9]

Are arrays faster than lists?

• Accessing array elements is easy & fast:
• Example: array “a” of N integers of size 64-bit
  • i\(^{th}\) element is at memory-address: (address of a[0]) + \(i \times 64\) bits
  • Copy/change M integers: memory chunk of \(M \times 64\) bits

Arrays are faster !

If you do not grow them !

\[ [a_1, .. , a_N] \rightarrow [a_1, .. , a_N, a_{N+1}] \]

• Extending an array of \(N\) elements with an extra element:
  • Creates a new array of \(N+1\) elements
  • Copies the old array into the new
  • Puts the extra element value in the last position

Array initialization by casting

import numpy as np

a = np.array([1, 2, 4, 6, 5, 9])
print(a)

[1 2 4 6 5 9]

a = np.array([[6, 0], [7, 4], [2, 9]])
print(a)

[[6 0]
 [7 4]
 [2 9]]

Array properties

• Shape of the array
• Array dimensions
• Type of the elements
  • Default: 64-bit floats
  • Other options: “uint32”: unsigned 32-bit integers
• Remember: all elements are of same type !

Get array properties

print(a)

[[6 0]
 [7 4]
 [2 9]]

print(a.shape)

(3, 2)

print(a.dtype)

int32

Adapt array properties: element type

print(a)
print(f"Type of the elements is: {a.dtype}")

[[6 0]
 [7 4]
 [2 9]]
Type of the elements is: int32

b = a.astype(np.complex128)
print(b)
print(f"Type of the elements is: {b.dtype}")

[[6.+0.j 0.+0.j]
 [7.+0.j 4.+0.j]
 [2.+0.j 9.+0.j]]
Type of the elements is: complex128

Array initialization

a = np.zeros((4, 3))       # 2D array with zeros
print(a)

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

a = np.ones(4)             # 1D array with ones
print(a)

[1. 1. 1. 1.]

Array initialization

a = np.full((5, 3), 1.2)   # 2D array filled with a value
print(a)

[[1.2 1.2 1.2]
 [1.2 1.2 1.2]
 [1.2 1.2 1.2]
 [1.2 1.2 1.2]
 [1.2 1.2 1.2]]

a = np.full((3, 3), [0.2, 0.3, 1.7])   # 2D array with values
print(a)

[[0.2 0.3 1.7]
 [0.2 0.3 1.7]
 [0.2 0.3 1.7]]

Array initialization: similar shape

a = np.ones((2, 4))       # 2D array with zeros
print(a)

[[1. 1. 1. 1.]
 [1. 1. 1. 1.]]

a = np.zeros_like(a)       # 2D array with zeros
print(a)

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

Array initialization: similar shape

b = np.full_like(a, 4.5)       # 2D array with zeros
print()
print(f"a = {a}")
print()
print(f"b = {b}")

a = [[0. 0. 0. 0.]
 [0. 0. 0. 0.]]

b = [[4.5 4.5 4.5 4.5]
 [4.5 4.5 4.5 4.5]]

Empty array initialization

a = np.empty((3, 2))      # empty is not empty !
print(a)

[[ 3.      0.34  ]
 [ 0.546  34.9   ]
 [ 5.3462  9.    ]]

nx, ny = a.shape
for i in range(nx):
  for j in range(ny):
    a[i, j] = i + 2*j
print(a)

[[0. 2.]
 [1. 3.]
 [2. 4.]]

Array initialization: intervals

a = np.arange(15)          # similar to range()
print(a)

[ 0  1  2  3  4  5  6  7  8  9 10 11 12 13 14]

a = np.arange(5, 7, 0.3)  # decimal float step
print(a)

[5.  5.3 5.6 5.9 6.2 6.5 6.8]

a = np.linspace(0, 5, 11)  # linearly spaced interval
print(a)

[0.  0.5 1.  1.5 2.  2.5 3.  3.5 4.  4.5 5. ]

2D domains/intervals: meshgrid

x = np.linspace(0, 5, 6)
y = 2 ** np.linspace(0, 3, 4)
xx, yy = np.meshgrid(x, y)

print(xx)

[[0. 1. 2. 3. 4. 5.]
 [0. 1. 2. 3. 4. 5.]
 [0. 1. 2. 3. 4. 5.]
 [0. 1. 2. 3. 4. 5.]]

print(yy)

[[1. 1. 1. 1. 1. 1.]
 [2. 2. 2. 2. 2. 2.]
 [4. 4. 4. 4. 4. 4.]
 [8. 8. 8. 8. 8. 8.]]

For more on domains see the Numpy meshgrid howto.

Higher dimensional arrays

• Multi-dimensional N-D arrays have N axes

tmp = np.array([[ 1,  2,  3],[ 4,  5,  6],[ 7,  8,  9]])
a = np.array([tmp, 2*tmp, 3*tmp])
print(a)

[[[ 1  2  3]
  [ 4  5  6]
  [ 7  8  9]]

 [[ 2  4  6]
  [ 8 10 12]
  [14 16 18]]

 [[ 3  6  9]
  [12 15 18]
  [21 24 27]]]

Element selection & Array slicing

• list element selection with [i]
• list slicing operator: [start:end+1:step]

a = np.array([[ 1,  2,  3],[ 4,  5,  6],[ 7,  8,  9]])
print(f"a[2, 0] = {a[2, 0]}")
print(f"a[0:2:2, :] = {a[0:2:2, :]}")

a[2, 0] = 7
a[0:2:2, :] = [[1 2 3]]

print(f"Squeezed a[0:2:2, :] = {np.squeeze(a[0:2:2, :])}")

Squeezed a[0:2:2, :] = [1 2 3]

Adapting the array shape

a = np.array([[1,5,6,2],[5,3,8,9],[1,1,1,0]])
print(a)

[[1 5 6 2]
 [5 3 8 9]
 [1 1 1 0]]

print(a.flatten())

[1 5 6 2 5 3 8 9 1 1 1 0]

print(a.reshape((2,6)))

[[1 5 6 2 5 3]
 [8 9 1 1 1 0]]

Array arithmetic operations

• similar to arithmetic operations on numbers
• element-wise
• shapes need to be compatible

a = np.array([2.5, 3, 4])   
b = np.array([1, 0.1, 10])
print(f"{a} + {b} = {a + b}")
print(f"{a} * {b} = {a * b}")
print(f"{a} / {b} = {a / b}")

[2.5 3.  4. ] + [ 1.   0.1 10. ] = [ 3.5  3.1 14. ]
[2.5 3.  4. ] * [ 1.   0.1 10. ] = [ 2.5  0.3 40. ]
[2.5 3.  4. ] / [ 1.   0.1 10. ] = [ 2.5 30.   0.4]

Array logic operations

• similar to logic operations on (boolean) numbers
• element-wise
• shapes need to be compatible

a = np.array([1, 3, 4])   
b = np.array([1, 1, 10])
print(f"{a} > {b} = {a > b}")
print(f"{a} == {b} = {a == b}")

[1 3 4] > [ 1  1 10] = [False  True False]
[1 3 4] == [ 1  1 10] = [ True False False]

Array mathematical functions

• Apply mathematical functions similar to math library
• Numpy functions, for example: math.sin() \(\rightarrow\) np.sin()

print(f"sin({a}) = {np.sin(a)}")

sin([1 3 4]) = [ 0.84147098  0.14112001 -0.7568025 ]

print(f"sqrt({a}) = {np.sqrt(a)}")

sqrt([1 3 4]) = [1.         1.73205081 2.        ]

print(f"exp({a}) = {np.exp(a)}")

exp([1 3 4]) = [ 2.71828183 20.08553692 54.59815003]

Vectors

A vector is a 1D array, e.g. a column or a row vector: \[ \vec{a} = \pmatrix{a_1\\a_2\\ \vdots \\ a_M} \qquad \vec{b} = \pmatrix{b_1 & b_2 & \dots & b_N} \]

a = np.array([[1],[3],[5]])
print(a)

[[1]
 [3]
 [5]]

Vectors

A vector is a 1D array, e.g. a column or a row vector: \[ \vec{a} = \pmatrix{a_1\\a_2\\ \vdots \\ a_M} \qquad \vec{b} = \pmatrix{b_1 & b_2 & \dots & b_N} \]

A row vector:

b = np.array([2, 4, 6, 8])
print(b)

[2 4 6 8]

The scalar product (dot-product)

The scalar product \(\vec{a} \cdot \vec{b}\) is defined as scalar value

\[ \vec{a} \cdot \vec{b} = |a| |b| \cos(\theta) \]

with \(\theta\) the angle between \(a\) and \(b\)

Alternative definition of the scalar product \(\vec{a} \cdot \vec{b}\):

\[ \vec{a} \cdot \vec{b} = \sum_{i=1}^d a_i \, b_i \] where \(d\) is the dimension: 2D, 3D, …, ND

The scalar product (dot-product)

Scalar product \(\vec{a} \cdot \vec{b}\) can be computed in Numpy using np.dot:

a = np.array([1, -1])   
b = np.array([1.5, 2.5])
s = np.dot(a, b)
print(f"{a} . {b} = {s}")

[ 1 -1] . [1.5 2.5] = -1.0

a = np.array([-1 + 1j, 1 - 1j])   
b = np.array([0, 2j])
print(f"{np.conj(a)} . {b} = {np.vdot(a, b)}")

[-1.-1.j  1.+1.j] . [0.+0.j 0.+2.j] = (-2+2j)

The scalar product (dot-product)

Scalar product \(\vec{a} \cdot \vec{b}\) can be computed in Numpy using np.dot:

a = np.array([1, -1])   
b = np.array([1.5, 2.5])
s = np.dot(a, b)
print(f"{a} . {b} = {s}")

[ 1 -1] . [1.5 2.5] = -1.0

a = np.array([-1, 2, 2, 2, 1])   
b = np.array([0, 0, -3, 3, 2])   
print(f"{a} . {b} = {np.dot(a, b)}")

[-1  2  2  2  1] . [ 0  0 -3  3  2] = 2

Norm of a vector

• Norm of a vector is the length of a vector
• Defined as square root of scalar product:

\[ |\vec{a}| = \sqrt{\vec{a}\cdot\vec{a}} = | \pmatrix{a_1 & a_1 & \dots & a_d} | = \sqrt{(\sum_{i=1}^d a_i^2)} \]

import numpy.linalg as la
a = np.array([1, 3, 0, 5])
print(f"The norm of {a} = {la.norm(a)}")

The norm of [1 3 0 5] = 5.916079783099616

Vector products (cross-product)

The vector product \(\vec{a} \times \vec{b}\) is defined as vector \(\vec{c}\) with

• Norm \(|\vec{a} \times \vec{b}| = |a| |b| \sin(\theta)\)
• Perpendicular to surface spanned by \(\vec{a}\) and \(\vec{b}\)

Vector products

Alternative definition (suitable for computers):

\[ \begin{aligned} \vec{a} \times \vec{b} & = \det\pmatrix{ \vec{e_1} & \vec{e_2} & \vec{e_3}\\ a_1 & a_2 & a_3 \\ b_1 & b_2 & b_3 \\} \\ & \\ & = \begin{vmatrix} a_2 & a_3\\ b_2 & b_3\\ \end{vmatrix}\vec{e_1} - \begin{vmatrix} a_1 & a_3\\ b_1 & b_3\\ \end{vmatrix} \vec{e_2} + \begin{vmatrix} a_2 & a_3\\ b_2 & b_3\\ \end{vmatrix} \vec{e_3} \\ & \\ & = (a_2 b_3 - a_3 b_2) \vec{e_1} - (a_1 b_3 - a_3 b_1) \vec{e_2} + (a_1 b_2 - a_2 b_1) \vec{e_2} \\ \end{aligned} \]

Vector products

Calculate \(\vec{a} \times \vec{b}\) using the Numpy function np.cross(a, b):

a = np.array([1, -1/2, 1])
b = np.array([1/2, 1, 1/2])
c = np.cross(a, b)
print(f"{a} x {b} = {c}")

[ 1.  -0.5  1. ] x [0.5 1.  0.5] = [-1.25  0.    1.25]

Array initialization: matrices

a = np.eye(3)             # 2D unity matrix
print(a)

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

# Constructing a matrix with "diagonals"
a = np.diagflat([6, 5, 4], 1) + np.diagflat([1, 2, 3], -1)
print(a)

[[0 6 0 0]
 [1 0 5 0]
 [0 2 0 4]
 [0 0 3 0]]

Matrix products

The matrix product between matrices \(A\) and \(B\) is defined as

\[ \begin{aligned} A \cdot B & = \pmatrix{ a_{11} & a_{12} & a_{13}\\ a_{21} & a_{22} & a_{23}\\ } \pmatrix{ b_{11} & b_{12} & b_{13}\\ b_{21} & b_{22} & b_{23}\\ b_{31} & b_{32} & b_{33}\\ }\\ &\\ & = \sum_{j} a_{ij} b_{jk}\\ \end{aligned} \]

• Rows of \(A\) are multiplied by columns of \(B\).
• \(A_{MN} \cdot B_{NK} \leftarrow\) No. columns of \(A\) must equal No. rows of \(B\)

Matrix products

The matrix product \(A \cdot B\) using Numpy np.matmul():

A = np.array([[1,-2,3],[0,2,4]])
B = np.array([[1,-2,3],[0,2,4],[-1,1,0]])
print(A)

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

print(B)

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

print(np.matmul(A, B))

[[-2 -3 -5]
 [-4  8  8]]

Matrix products

Alternatives to np.matmul(A, B) for \(A \cdot B\) using Numpy:

• @ operator: A @ B, and
• np.dot(A, B)

print(np.matmul(A, B))

[[-2 -3 -5]
 [-4  8  8]]

print(A @ B)

[[-2 -3 -5]
 [-4  8  8]]

print(np.dot(A, B))

[[-2 -3 -5]
 [-4  8  8]]

Determinant of a matrix

The determinant of a matrix \(A\): \[ \begin{aligned} \det(A) & = \det\pmatrix{ a_{11} & a_{12} & a_{13}\\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \\} \\ & \\ & = \begin{vmatrix} a_{22} & a_{23}\\ a_{32} & a_{33}\\ \end{vmatrix} a_{11} - \begin{vmatrix} a_{21} & a_{23}\\ a_{31} & a_{33}\\ \end{vmatrix} a_{12} + \begin{vmatrix} a_{21} & a_{22}\\ a_{31} & a_{32}\\ \end{vmatrix} a_{13} \\ &\\ & = (a_{22} a_{33} - a_{23} a_{32}) a_{11} - \dots \\ \end{aligned} \]

Determinant of a matrix

Use Numpy numpy.linalg.det() function

import numpy.linalg as la
A = np.array([[1,-2,3],[0,2,4],[-1,1,0]])
print(A)

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

Determinant is given by

print(f"Det(A) = {la.det(A)}")

Det(A) = 9.999999999999998

Inverse of a matrix

The inverse \(A^{-1}\) of a matrix \(A\) is defined as:

\[ A^{-1} A = A^{-1} A = \mathbb{1} \]

import numpy.linalg as la
A = np.array([[1,-2,3],[0,2,4],[-1,1,0]])
print(A)
print(f"Det(A) = {la.det(A)}")

[[ 1 -2  3]
 [ 0  2  4]
 [-1  1  0]]
Det(A) = 9.999999999999998

Inverse of a matrix

The inverse \(A^{-1}\) of a matrix \(A\) is defined as:

\[ A^{-1} A = A^{-1} A = \mathbb{1} \] Since \(\det(A) \neq 0 \rightarrow A^{-1}\) exists

Ainv = la.inv(A)
print(Ainv)

[[-0.4  0.3 -1.4]
 [-0.4  0.3 -0.4]
 [ 0.2  0.1  0.2]]

Inverse of a matrix

The inverse \(A^{-1}\) of a matrix \(A\) is defined as:

\[ A^{-1} A = A^{-1} A = \mathbb{1} \]

Verify \(A^{-1} A = \mathbb{1}\):

np.dot(A, Ainv)

array([[ 1.00000000e+00, -2.77555756e-17, -5.55111512e-17],
       [ 0.00000000e+00,  1.00000000e+00,  0.00000000e+00],
       [ 1.11022302e-16,  5.55111512e-17,  1.00000000e+00]])

np.round(np.dot(A, Ainv), 5) # Round numerical errors

array([[ 1., -0., -0.],
       [ 0.,  1.,  0.],
       [ 0.,  0.,  1.]])