Instruction on Mathematical Programming in Python for COMP24112

Disclaimer: This is not a tutorial for complete beginners in Python. It assumes you already know the language at a very basic level. You should go through the official introductory tutorial first (https://docs.python.org/3/tutorial/introduction.html), if you haven’t already. You can use the following cell as an interactive console to learn hands-on, although we recommend you to use Spyder’s one.

Write your code here

Table of Contents

• Control structures
• Matrix programming
• Creating vectors and matrices
• Accessing elements of arrays
• Creating random matrices
• Matrix attributes and functions
• Operators on matrices
• Plotting

1. Control structures

We will show the basics constructs for control flow in Python. If you want to know more we redirect you to the official guide: https://docs.python.org/3/tutorial/controlflow.html.

If

x = 5

if x < 10:
print(‘Less than 10.’)

if x < 10:
print(‘Less than 10.’)
else:
print(‘At least 10.’)

if x < 10:
print(‘Less than 10.’)
elif x < 20:
print(‘Less than 20.’)
elif x < 30:
print(‘Less than 30.’)
else:
print(‘At least 30.’)

For

for i in range(5): # Note how 5 is excluded
print(i, end=’ ')
print()

for i in range(3, 7): # Again, 7 is excluded
print(i, end=’ ')
print()

for i in range(3, 12, 2): # Step 2
print(i, end=’ ')
print()

for i in range(5):
if i > 2:
break
print(i, end=’ ')
print()

for i in range(5):
if i > 7: # This will never occur.
break
print(i, end=’ ')
else: # Python allows an else clause to the for loop that is executed when no break occurred.
print(‘Went through the whole loop!’)
print()

User-defined functions

Defining functions for code reuse and structure is of primary importance. Generally, copy-pasting code signals that you should probably write a function instead. Here are the most common, yet most useful, ways of defining a function in Python.
def my_function():
# Write your code here
print(‘This function is not very useful.’)

my_function()
def my_sum(a, b):
s = a + b
return s

print(my_sum(5, 9))
def division_and_reminder(a, b):
d = a // b # // is the integer division operator
r = a % b # % is the modulo operator
return d, r

d, r = division_and_reminder(13, 3)
print(d, r)
def exp(x, base=10): # You can assign values as defaults
return base ** x

print(exp(2))
print(exp(2, base=3))

2. Matrix programming

A matrix is a rectangular array of numbers, like a two-dimensional array in C or JAVA. An ‘’ n n n by m m m" (or n × m n\times m n×m) matrix has n n n rows and m m m columns. Special meaning is sometimes attached to 1 by 1 matrices, which are called ‘‘scalars’’ (ordinary numbers, basically), and to matrices with only one row or only one column, which are called ‘‘vectors’’.

In Python, the most common way of dealing with matrices is using NumPy’s arrays. They represent multi-dimensional arrays of a fixed size. This means that they can be used to represent scalars (0D), vectors (1D), matrices (2D) and even higher-order arrays (but you will not need to use them in this course).

First of all, we must include NumPy.
import numpy as np

2.1 Creating vectors and matrices

To manually instantiate a vector in NumPy use the following syntax:
a = np.array([1, 10, 100])
print(a)
If you have gone through the Python introductory tutorial (as you should), you have probably realised that we have called the np.array() function by passing a list as an argument. This is why the sequence of number is delimited by brackets [ ] and elements are separated by a comma.

To define a matrix you can run:
A = np.array([[0, 1, 2],
[3, 4, 5]
])

Note: new lines are not important here. We could write: np.array([[0, 1, 2], [3, 4, 5]]).

print(A)
We are instantiating a list of lists first, and then passing it as an argument to np.array(). Each of the inner lists defines a row. Note that each row must have the same number of elements, otherwise the resulting object will not be a matrix:

THIS IS WRONG!

print(np.array([[0, 1, 2],
[3, 4]
]))
However, you will rarely need to define a vector or a matrix manually. Most of the time you will use one of the built-in functions to create specific types of arrays:
range_v = np.arange(10) # Array of integers from 0 to 9
print(range_v)
X = np.zeros([2, 3]) # 2x3 matrix full of zeros
print(X)
y = np.ones(4) # array of length 4 full of ones
print(y)
I = np.eye(3) # 3x3 identity matrix
print(I)
There are a few things to note here. First, the shape must be passed as a list in order to create a matrix: np.zeros(2, 3) would return an error. On the other hand, for a vector you can pass a single scalar to represent the length, as we did in y = np.ones(4). Second, the default data type is float. In previous examples the resulting arrays where integers because NumPy automatically inferred it from the constants we used. Compare these:
int_vector = np.array([0, 1, 2]) # All constants are integers.
print(int_vector)
float_vector = np.array([0., 1., 2.5]) # All constants are floats. Note: “1.” is the same as “1.0” in Python.
print(float_vector)
another_float_vector = np.array([0, 1.7, 2]) # Integers are converted to floats when mixed.
print(another_float_vector)
To force NumPy to use integers you can use the keyword argument dtype:
y_float = np.ones(4)
print(y_float)
y_int = np.ones(4, dtype=int)
print(y_int)
Finally, np.zeros() and np.ones() have their counterparts that create an array copying dimensions and data type from an existing one:
z = np.zeros_like(y_int) # y_int was defined previously as an array of integers of length 4, so z is too.
print(z)
Y = np.ones_like(X) # X was a 2x3 matrix of zeros.
print(Y)

2.2 Accessing elements of arrays

You can access elements of a NumPy array exactly like you would with a list or string:
a = np.array([5, 50, 500])
print(a[0], a[1], a[2], a[-1])
For matrices you can do the same, but it’s much clearer and more convenient to use NumPy indexing:
A = np.array([[0, 1, 2],
[3, 4, 5]
])
print(A[0][0], A[0][2], A[1][-2]) # List indexing. You should not use this.
print(A[0, 0], A[0, 2], A[1, -2]) # NumPy indexing. Use this.
You can use the syntax start:stop to select a part of a matrix. Note that start is included, while stop is excluded. This is the same as Python’s slicing that you have seen in the tutorial.
B = np.array([[ 0, 1, 2, 3],
[ 4, 5, 6, 7],
[ 8, 9, 10, 11]
])
print(B[0:2, 1:3])
If you omit start, Python will start from the beginning. If you omit stop, it will continue until the end. If you omit both, it will do both, thus selecting everything.
print(B[:2, 1:3]) # Note: this is the same as before, because 0 is the starting value.
print(B[:2, 1:])
print(B[:2, :])
You can also use negative indices, as usual:

These are the same because 3 is the last column index.

print(B[:2, :3])
print(B[:2, :-1])
You can also specify a step by using the start:stop:step syntax, if you want non-unit spacing (note that step is the last element!):
C = np.array([[ 0, 1, 2, 3, 4, 5, 6, 7],
[ 8, 9, 10, 11, 12, 13, 14, 15]
])
print(C[:, ::2])
print(C[:, 1:6:2])
If you use a negative step you go backwards:
print(C[:, ::-1]) # From last element to first
print(C[:, ::-2]) # Same as above, but with step two
print(C[:, 6:2:-2]) # Also specify start and end
You can also use lists or other arrays of integers to index.
print(B)
print(B[0::2, :]) # Even numbered rows.
print(B[[0, 2], :]) # Rows 0 and 2. Since matrix B has 3 rows, this is effectively the same as above.
row_inds = np.array([0, 2])
print(B[row_inds, :]) # NumPy arrays work as well.
Note that NumPy supports a more advanced way of indexing, that allows you to use matrices as indices as well. The example below shows one example where advanced indexing is useful. You can read more about it on the NumPy documentation.
a = np.array([5, 50, 500]) # a is the array that contains the values we are interested in.
I = np.array([[0, 1, 2], [2, 0, 1]]) # I is a 2x3 matrix containing indices over array a.

R is a 2x3 matrix (same shape as I) containing elements taken from array a.

For example, R[0, 1] is defined as a[I[0, 1]].

R = a[I]

print(a)
print(I)
print()
print®

2.3 Creating random matrices

In Python there are two common ways of dealing with randomness: use the standard random library or NumPy’s one.
import random # Standard library
print(random.random())

Now we use NumPy’s random library.

print(np.random.rand())
It is generally more convenient to use NumPy because it makes it easier to create and handle random vectors and matrices.

We just used np.random.rand() without arguments and saw that it returns a single scalar. We can use it to create vectors or matrices as well:
print(np.random.rand(5)) # vector of length 5
print(np.random.rand(2, 3)) # 2x3 matrix
np.random.rand() samples from a uniform distribution over [ 0 , 1 ) [0, 1) [0,1). You can use np.random.randn() to sample from a standard normal/Gaussian distribution instead:
print(np.random.randn(2, 3))
For integers you can use np.random.randint(). There are a few ways of calling it:
print(np.random.randint(5)) # Single integer in [0, 5).
print(np.random.randint(3, 10)) # Single integer in [3, 10).
print(np.random.randint(6, size=10)) # Vector of length 10 with integers in [0, 6).
print(np.random.randint(3, 10, size=[3, 4])) # 3x4 matrix of integers in [3, 10).
To perform permutations you can use np.random.permutation(). You can use it in two ways:
print(np.random.permutation(4)) # Random permutation of [0, 1, 2, 3]
print(np.random.permutation([93, 2, 1919, -38])) # Random permutation of the input list.
You can select just some elements of the result if you want a sample:

Sample 4 different numbers in [0, 10).

rnd = np.random.permutation(10)
print(rnd[:4])

2.4 Matrix attributes and functions

NumPy arrays have quite a few useful attributes.

Data type

y_float = np.ones(4)
print(y_float.dtype)
y_int = np.ones(4, dtype=int)
print(y_int.dtype)

Shape and size

A = np.array([[0, 1, 2],
[3, 4, 5]
])
print(A.shape) # Get all dimensions.
print(A.shape[0]) # Get only the first dimension.
print(A.size) # This is the number of elements in the matrix, i.e., the product of its dimensions.
Please note that, in NumPy, vectors are NOT exactly equivalent to matrices with a dimension equal to one.
a = np.array([0, 1, 2])
print(a.shape)
b = np.array([[0, 1, 2]]) # Note the double brackets.
print(b.shape)

Transposition

print(A.T)
For vectors, transposition only makes sense with 2-dimensional ones:
print(a)
print(a.T) # This is exactly the same.
print(b)
print(b.T) # This is not!

Diagonals

I = np.eye(3)
print(I.diagonal())
print(np.diag(I)) # Alternative way
np.diag() can also be used to construct a diagonal matrix, like in Matlab:
D = np.diag([1, 2, 4, 8])
print(D)

Computation functions

B = np.array([[ 0, 1, 2, 3],
[ 4, 5, 6, 7],
[ 8, 9, 10, 11]
])
print(np.min(B))
print(np.max(B))
print(np.sum(B))
print(np.mean(B))
Note that all of these can be called with the keyword argument axis to perform the computation along the given axis instead of over the whole matrix. Using np.sum() as an example:
print(np.sum(B, axis=0)) # Over rows, i.e., compute the sum of each column.
print(np.sum(B, axis=1)) # Same, but over columns.

Other useful functions

Inverse

print(np.linalg.inv(D))

Elementwise square root

print(np.sqrt(B))

Sort

R = np.random.randint(10, size=[3, 5])
print®
print(np.sort(R, axis=0)) # Over rows, so it sorts each column.

2.5 Operators on matrices

NumPy array support most Python arithmetic and logic operators. We will show the most useful using these two matrices:
M = np.array([[ 1, 2, 3, 4],
[ 5, 6, 7, 8],
[ 9, 10, 11, 12],
[13, 14, 15, 16],
])
D = np.diag([1, 2, 4, 8])

Arithmetic operators

Standard arithmetic operators still work on NumPy arrays in a component-by-component or element-by-element fashion:
print(D + M)
print(D - M)
print(D * M)
print(D / M)
Please note the automatic casting to float when necessary, in the division just performed.
print(D ** 2) # Square. ** is Python’s power operator.
Note that ** can be used as an alternative method to perform the square root:
print(D ** 0.5)

Logical operators

Logical operators also work component-by-component:
print(D > M) # All other comparison operator such as <, <=, ==, >= and != work as well.
print(D > 2)
print(M > 10)

Logical AND

print((D > 2) & (M > 10)) # Parentheses are important here!

Logical OR

print((D > 2) | (M > 10))

Logical NOT

print(~(D > 2)) # Obviously, in this case we could have written D <= 2.
Note that you should never use the == operator to compare floats or array of floats (Python or in any other language). Use NumPy’s allclose() instead:
print((.1 + .2) == .3)
print(np.allclose(.1 + .2, .3))

Matrix multiplication

print(D)
print(M)
print(D @ M)
print(M @ D) # Remember that matrix multiplication is not commutative!
When multiplying a matrix with a vector (remember, vectors only have one dimension), NumPy automatically treats it as a column or row vector as appropriate:
v = np.array([1, 2, 3, 4])
print(M @ v)
print(v @ M)
Compare it with the following:
v_row = np.array([[1, 2, 3, 4]])
print(v_row @ M) # This is the same…
#print(M @ v_row) # … but this raises an error!
Matrix multiplication can be also used to compute the dot product between two vectors:
u = np.array([10, 20, 30, 40])
print(u @ v)

Concatenating

You can concatenate multiple vectors using the function np.concatenate():
print(v)
vv = np.concatenate([v, v])
print(vv)
print(np.concatenate([v, v, v])) # You can concatenate as many arrays as you want.
This also works on matrices:
print(np.concatenate([D, D]))
You may have notices that the default behaviour is concatenating over the first dimension, which for matrix means along rows. You can change it through the keyword axis:
#print(np.concatenate([D, D], axis=0)) # This is the default behaviour, i.e., the same as above.
print(np.concatenate([D, D], axis=1)) # This concatenates along columns instead.
You can also concatenate matrix and vectors, as long as it makes sense:
print(v_row.shape, D.shape)
print(np.concatenate([v_row, D], axis=0)) # This works.
#print(np.concatenate([v_row, D], axis=1)) # This doesn’t.
However, you cannot concatenate a 1-dimensional vector with a matrix:
print(v.shape, D.shape)
#print(np.concatenate([v, D], axis=0)) # This raises an error.

Broadcasting

Broadcasting is a handy feature of NumPy that relaxes shape constraints. The simplest case of broadcasting is using a constant with a binary matrix operator:
print(D)
print(D + 5)
As you can see, the constant is added to every element of the matrix and you don’t need to create a 4$\times$4 matrix full of 5 to get the same result.

Broadcasting also applies to dimensions of size 1. For example:
print(v_row)
print(D.shape, v_row.shape)
print(D + v_row)
Since vector v_row was added to the 4$ \times 4 m a t r i x ‘ D ‘ h a s b e e n i m p l i c i t l y e x p a n d e d t o a 4 4 matrix `D` has been implicitly expanded to a 4 4matrix‘D‘hasbeenimplicitlyexpandedtoa4 \times $4 matrix by duplicating it across dimension 0, which has size 1. If you want to check the result explicitly you can do:
print(np.broadcast_to(v_row, D.shape))
When broadcasting to a higher dimension, empty dimensions will be prepended to the array. For instance, adding a vector with shape (n,) to a two-dimensional matrix is the same as adding a (row) vector with shape (1, n) - note that this will be broadcast differently from a (column) vector with shape (n, 1)! This also means that a vector with shape (n,) can be broadcast to (m, n) but not to (n, m).

Consider for example:
print(v)
print(v.shape)
print(D + v)
Thus, you need to be careful when doing this. You might mess up when dealing with column vectors:
print(D + v_row) # Adding a row vector
print(D + v_row.T) # Adding a column vector. The result is different.
print()
print(D + v) # Adding a 1-dimensional vector. It is treated as a row vector.
print(D + v.T) # Still adding a 1-dimensional vector!

Adding and removing dimensions

NumPy allows you to easily expand an array by adding dimensions of size 1:
print(v.shape)

v_exp_row = v[None, :] # Add dimension 0.
print(v_exp_row.shape)
print(v_exp_row) # Content is the same.

v_exp_column = v[:, None] # Add dimension 1.
print(v_exp_column.shape)
This is mostly useful to 1) avoid mistakes when broadcasting and 2) concatenating:
#print(np.concatenate([v, D], axis=0)) # This raises an error.
print(np.concatenate([v[None, :], D], axis=0)) # This works.
You can remove dimensions of size 1 by using the squeeze function:
print(v_row.shape)
v_row_reduced = np.squeeze(v_row)
print(v_row_reduced.shape)
print(v_row_reduced) # Content is the same.

3. Plotting

The matplotlib library contains a range of useful functions for data visualisation.
from matplotlib import pyplot as plt

Standard plot

x = np.arange(10)
y = x ** 2
plt.plot(x, y)
plt.xlabel(‘x’)
plt.ylabel(‘square’)
plt.title(‘A plot’)
plt.show()

Change plot style

plt.plot(x, y, ‘r–o’) # r stands for ‘red’, -- means dashed line, o enables round markers.
plt.xlabel(‘x’)
plt.ylabel(‘square’)
plt.show()
You can find a comprehensive list of possible options in the official documentation (scroll down a bit).

Multiple plots

In same figure. Note that Matplotlib automatically assigns different colours.

from math import pi
x = np.arange(0, 2*pi, 0.01) # From 0 to 2pi with step 0.01.
y_cos = np.cos(x)
plt.plot(x, y_cos, label=‘Cosine’)

y_sin = np.sin(x)
plt.plot(x, y_sin, label=‘Sine’)

plt.xlabel(‘x’)
plt.legend()

plt.show()

In different figures

plt.plot(x, y_cos, ‘–’, label=‘Cosine’)
plt.xlabel(‘x’)
plt.legend()

plt.figure()

plt.plot(x, y_sin, ‘r’, label=‘Sine’)
plt.xlabel(‘x’)
plt.legend()

plt.show()

Error bars

x = np.arange(0, 5, 0.5)
y = x ** 2
y_err = np.random.rand(y.size) * 5 # A vector of random numbers in [0, 5) with the same length as y
plt.errorbar(x, y, yerr=y_err, ecolor=‘r’)
plt.show()