198812 Session 2: Matrix algebra
Math to code mapping.
The purpose of the document it to show how mathematical conpects learned in the Computer Programming Prerequisites course could be implemented in a programming language. This supplementary material should help you to transfer theoretical concepts to programming. If you have programming skills you may use the material during the course. If not, just come back to it, once you are advanced enough in programming.
- This document is basesd on the lecture notes by Sébastien Court. It structured in the same way as the lecture notes. It consists only of sections, where numeric examples were given.
- It provides code in Python and R for examples used there. The code serves only demonstration purposes to show a possible implementation. We assumed the correctness of input data.
- For multiple examples of the same calculation, only one is taken. You can easily adapt code for missing calculations.
For this session, code in basic Python, Python with NumPy package and in R is given. Vector and matrix algebra is relatively easy to implement using numpy (some courses from module 3) or R (198803). In most of the cases just by calling dedicated functions. In basic Python (198801), the required functionality in many cases must be implemented from scratch.
Detailed description on how to work with matrices in R is provided in our textbook.
1 Definition and basic operations
Defining an example matrix.
\[ A = \begin{pmatrix} 7.9 & 3.2 & 6 \\ 1.3 & 0 & -4 \\ \end{pmatrix} \]
We define it as a list of lists:
a = [[7.9, 3.2, 6],[1.3, 0, -4]]; a[[7.9, 3.2, 6], [1.3, 0, -4]]We define it a an array:
import numpy
A = numpy.array(a)A =array([[ 7.9, 3.2, 6. ],
[ 1.3, 0. , -4. ]])shape = (2, 3)We define it as matrix:
a <- matrix(data = c(7.9, 3.2, 6,1.3, 0, -4), nrow = 2, byrow = TRUE) [,1] [,2] [,3]
[1,] 7.9 3.2 6
[2,] 1.3 0.0 -4Basic Operations
We take the addition as an example operation.
\[ \begin{pmatrix} 2 & 3 \\ 0 & -6 \\ \end{pmatrix} + \begin{pmatrix} -2.5 & 0 \\ 12.3 & -7.4 \\ \end{pmatrix} = \begin{pmatrix} -0.5 & 3 \\ 12.3 & -13.4 \\ \end{pmatrix} \]
We need a nested for loop to add two matrices.
a[[2, 3], [0, -6]]b[[-2.5, 0], [12.3, -7.4]]a + b =
nrow = len(a)
ncol = len(a[0])
result = []
for i in range(nrow):
row = []
for j in range(ncol):
row.append(a[i][j] + b[i][j])
result.append(row)
result [[-0.5, 3], [12.3, -13.4]]We just add two arrays with + or the numpy.add() function:
Aarray([[ 2, 3],
[ 0, -6]])Barray([[-2.5, 0. ],
[12.3, -7.4]])A + Barray([[ -0.5, 3. ],
[ 12.3, -13.4]])numpy.add(A, B)array([[ -0.5, 3. ],
[ 12.3, -13.4]])We define it as matrix:
a [,1] [,2]
[1,] 2 3
[2,] 0 -6b [,1] [,2]
[1,] -2.5 0.0
[2,] 12.3 -7.4a + b [,1] [,2]
[1,] -0.5 3.0
[2,] 12.3 -13.42 Subsetting and notation
2.1-3 Subsetting notation
In math and programming languages we may access single elements of composite data. The notation on how to do it depends on the language.
Matrix:
\[ A = \begin{pmatrix} 3 & 6 & 0 \\ 8 & 7 & 1 \\ \end{pmatrix} \]
Elements:
\[ \begin{matrix} a_{11} = 3, & a_{12} = 5, & a_{13} = 0 \\ a_{21} = 8, & a_{22} = 7, & a_{23} = 1 \\ \end{matrix} \]
Matrix (list of lists):
a[[3, 6, 0], [8, 7, 1]]Elements:
a[0][0] = 3 a[0][1] = 6 a[0][2] = 0
a[1][0] = 8 a[1][1] = 7 a[1][2] = 1 Matrix (array):
Aarray([[3, 6, 0],
[8, 7, 1]])Elements:
A[0, 0] = 3 A[0, 1] = 6 A[0, 2] = 0
A[1, 0] = 8 A[1, 1] = 7 A[1, 2] = 1 An array may be accessed the same way as a list, too.
A[0][0] = 3 A[0][1] = 6 A[0][2] = 0
A[1][0] = 8 A[1][1] = 7 A[1][2] = 1 Matrix:
a [,1] [,2] [,3]
[1,] 3 6 0
[2,] 8 7 1Elements:
a[1, 1] = 3 a[1, 2] = 6 a[1, 3] = 0
a[2, 1] = 8 a[2, 2] = 7 a[2, 3] = 1 2.4 Subsetting sub-matrices
In math and programming languages we may access parts of composite data.
\[ \begin{eqnarray*} A[2, 1] & = & 8 &&& \text{- a single element} \\ A[:,1] & = & \displaystyle \left(\begin{matrix}3 \\ 8\end{matrix}\right) &&& \text{- a submatrix, the first column} \\ A[2,:] & = & (8, 7, 1) &&& \text{- a submatrix, the second row}\\ A[:,(1, 3)] & = & \displaystyle \left(\begin{matrix}3 & 0 \\ 8 & 1 \end{matrix}\right) &&& \text{- a submatrix}\\ \end{eqnarray*} \]
Matrix (list of lists):
a[[3, 6, 0], [8, 7, 1]]Selecting a column must be done programmatically, the shortest is via the list comprehension.
[row[0] for row in a][3, 8]Selecting a row with indexing.
a[1][8, 7, 1]Selecting a submatrix must be done programmatically, too.
[[row[0], row[2]] for row in a][[3, 0], [8, 1]]Matrix (array):
Aarray([[3, 6, 0],
[8, 7, 1]])Selecting a column
A[:, 0]array([3, 8])Selecting a row
A[1, :]array([8, 7, 1])Selecting a submatrix
A[:, [0, 2]]array([[3, 0],
[8, 1]])Matrix:
a [,1] [,2] [,3]
[1,] 3 6 0
[2,] 8 7 1Selecting a column
a[, 1][1] 3 8Selecting a row
a[2, ][1] 8 7 1Selecting a submatrix
a[, c(1, 3)] [,1] [,2]
[1,] 3 0
[2,] 8 13 The product Matrix/Vector
3.1-4 … In any dimension
\[ \begin{eqnarray*} & & A = \left( \begin{matrix} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 9 \end{matrix} \right), \quad X = \left( \begin{matrix} -2 \\ 1 \\ 0 \end{matrix} \right), \\[5pt] & & AX = \left( \begin{matrix} 1\times (-2) + 2\times 1 + 3 \times 0\\ 4\times (-2)+ 5\times 1 + 6 \times 0 \\ 7\times (-2) + 8 \times 1 + 9 \times 0 \end{matrix} \right) = \left( \begin{matrix} 0 \\ -3 \\ -6 \end{matrix} \right). \end{eqnarray*} \]
a[[1, 2, 3], [4, 5, 6], [7, 8, 9]]x[-2, 1, 0]a * x =
ncol = len(x)
nrow = len(a)
result = [0] * ncol
for i in range(nrow):
for j in range(ncol):
result[i] += a[i][j] * x[j]
result [0, -3, -6]Aarray([[1, 2, 3],
[4, 5, 6],
[7, 8, 9]])Xarray([-2, 1, 0])numpy.dot(A, X)array([ 0, -3, -6])a [,1] [,2] [,3]
[1,] 1 2 3
[2,] 4 5 6
[3,] 7 8 9x[1] -2 1 0a %*% x [,1]
[1,] 0
[2,] -3
[3,] -64 The product Matrix/Matrix
4.1 Definition
\[ \begin{eqnarray*} & & A = \left( \begin{matrix} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 9 \end{matrix} \right), \quad B = \left( \begin{matrix} -2 & 5 & -1 \\ 1 & 0 & -3 \\ 7 & 0 & -9 \end{matrix} \right), \\[5pt] & & AB = \left( \begin{matrix} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 9 \end{matrix} \right) \left( \begin{matrix} -2 & 5 & -1 \\ 1 & 0 & -3 \\ 7 & 0 & -9 \end{matrix} \right) = \left( \begin{matrix} 21 & 5 & -34 \\ 39 & 20 & -73 \\ 57 & 35 & -112 \end{matrix} \right). \end{eqnarray*} \]
a[[1, 2, 3], [4, 5, 6], [7, 8, 9]]b[[-2, 5, -1], [1, 0, -3], [7, 0, -9]]As we will need this code for later examples, we will define our multiply() function.
def multiply(a, b):
nrow_a = len(a)
nrow_b = len(b)
ncol_b = len(b[0])
# preparing the result matrix filled with 0
result = []
zeros = [0] * ncol_b
for _ in range(nrow_a):
result.append(zeros[:])
# a * b
for i in range(nrow_a):
for j in range(ncol_b):
for k in range(nrow_b):
result[i][j] += a[i][k] * b[k][j]
return result
multiply(a, b)[[21, 5, -34], [39, 20, -73], [57, 35, -112]]Aarray([[1, 2, 3],
[4, 5, 6],
[7, 8, 9]])numpy.dot(A, B)array([[ 21, 5, -34],
[ 39, 20, -73],
[ 57, 35, -112]])a [,1] [,2] [,3]
[1,] 1 2 3
[2,] 4 5 6
[3,] 7 8 9b [,1] [,2] [,3]
[1,] -2 5 -1
[2,] 1 0 -3
[3,] 7 0 -9a %*% b [,1] [,2] [,3]
[1,] 21 5 -34
[2,] 39 20 -73
[3,] 57 35 -1125 Inversion of a matrix
5.1 The identity matrix
\[ \begin{eqnarray*} I_3 & = & \left( \begin{matrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{matrix} \right) \end{eqnarray*} \]
We may either define the required identity matrix manually or write a for loop to create it.
[[1, 0, 0], [0, 1, 0], [0, 0, 1]][[1, 0, 0], [0, 1, 0], [0, 0, 1]]n = 3
result = []
for i in range(n):
row = []
for j in range(n):
row.append(int(i == j))
result.append(row)
result [[1, 0, 0], [0, 1, 0], [0, 0, 1]]numpy.identity(3, dtype=int)array([[1, 0, 0],
[0, 1, 0],
[0, 0, 1]])diag(3) [,1] [,2] [,3]
[1,] 1 0 0
[2,] 0 1 0
[3,] 0 0 15.2 Invertible matrices
For
\[ \begin{eqnarray*} A = \left( \begin{matrix} 4 & 1 \\ -1 & -0.2 \end{matrix} \right), & & B = \left( \begin{matrix} -1 & -5 \\ 5 & 20 \end{matrix} \right) \end{eqnarray*} \]
one can verify that
\[ \begin{equation*} AB \quad = \quad BA \quad = \quad I_2. \end{equation*} \]
a[[4, 1], [-1, -0.2]]b[[-1, -5], [5, 20]]multiply(a, b)[[1, 0], [0.0, 1.0]]multiply(b, a)[[1, 0.0], [0, 1.0]]Aarray([[ 4. , 1. ],
[-1. , -0.2]])Barray([[-1, -5],
[ 5, 20]])numpy.dot(A, B)array([[ 1.00000000e+00, 0.00000000e+00],
[-5.55111512e-17, 1.00000000e+00]])numpy.dot(B, A)array([[1.00000000e+00, 5.55111512e-17],
[0.00000000e+00, 1.00000000e+00]])Due to impreciseness we have insignificant difference of the obrained and expected results. Let us just round the obtained results.
numpy.around(numpy.dot(A, B), decimals=15)array([[ 1., 0.],
[-0., 1.]])numpy.around(numpy.dot(B, A), decimals=15)array([[1., 0.],
[0., 1.]])a [,1] [,2]
[1,] 4 1.0
[2,] -1 -0.2b [,1] [,2]
[1,] -1 -5
[2,] 5 20a %*% b [,1] [,2]
[1,] 1 0
[2,] 0 1b %*% a [,1] [,2]
[1,] 1 0
[2,] 0 16 Other products
6.1 Componentwise
\[ \begin{eqnarray*} & & A = \left(\begin{matrix} 3 & -1 \\ -4 & 5 \end{matrix} \right), \quad B = \left(\begin{matrix} 0 & 5 \\ -1 & 7 \end{matrix} \right), \quad A \ast B = \left(\begin{matrix} 0 & -5 \\ 4 & 35 \end{matrix} \right). \end{eqnarray*} \]
a[[3, -1], [-4, -5]]b[[0, 5], [-1, 7]]a * b =
nrow = len(a)
ncol = len(a[0])
result = []
for i in range(nrow):
row = []
for j in range(ncol):
row.append(a[i][j]*b[i][j])
result.append(row)
result[[0, -5], [4, -35]]Aarray([[ 3, -1],
[-4, -5]])Barray([[ 0, 5],
[-1, 7]])A * Barray([[ 0, -5],
[ 4, -35]])a [,1] [,2]
[1,] 3 -1
[2,] -4 -5b [,1] [,2]
[1,] 0 5
[2,] -1 7a * b [,1] [,2]
[1,] 0 -5
[2,] 4 -356.2 The tensor product
\[ \begin{eqnarray*} & & v = \left(\begin{matrix} -2 \\ 3 \\ 0 \\ 1 \end{matrix} \right), \quad w = \left(\begin{matrix} 1 & 4 & -7 & 9 \end{matrix} \right), \\[5pt] & & v \otimes w = \left(\begin{matrix} -2 & -8 & 14 & -18 \\ 3 & 12 & -21 & 27 \\ 0 & 0 & 0 & 0 \\ 1 & 4 & -7 & 9 \end{matrix} \right). \end{eqnarray*} \]
v[-2, 3, 0, 1]w[1, 4, -7, 9]v \(\otimes\) w =
nrow = len(v)
ncol = len(w)
result = []
for i in range(nrow):
row = []
for j in range(ncol):
row.append(v[i]*w[j])
result.append(row)
result[[-2, -8, 14, -18], [3, 12, -21, 27], [0, 0, 0, 0], [1, 4, -7, 9]]We use the numpy.tensordot() function. The axis determine the type of tensors to be are calculated. We use axis=0 to get the tensor product.
Varray([-2, 3, 0, 1])Warray([ 1, 4, -7, 9])numpy.tensordot(V, W, 0) array([[ -2, -8, 14, -18],
[ 3, 12, -21, 27],
[ 0, 0, 0, 0],
[ 1, 4, -7, 9]])It is possible to use outer() function or the tensor operator %o% from the tensor package.
v[1] -2 3 0 1w[1] 1 4 -7 9outer(v, w) [,1] [,2] [,3] [,4]
[1,] -2 -8 14 -18
[2,] 3 12 -21 27
[3,] 0 0 0 0
[4,] 1 4 -7 9library("tensor")
v %o% w [,1] [,2] [,3] [,4]
[1,] -2 -8 14 -18
[2,] 3 12 -21 27
[3,] 0 0 0 0
[4,] 1 4 -7 9