198812 Session 5: Matrix algebra (cont.)

Math to code mapping.

Author

Joanna Chimiak-Opoka

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.

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.

Matrices in R

Detailed description on how to work with matrices in R is provided in our textbook.

1 Symmetric matrices

1.1 Transposition

\[ \begin{eqnarray*} A = \left(\begin{matrix} 1 & 2 & 3 \\ 4 & 5 & 6 \end{matrix} \right) \quad \Rightarrow \quad A^T = \left(\begin{matrix} 1 & 4 \\ 2 & 5 \\ 3 & 6 \end{matrix} \right) \\ \end{eqnarray*} \]

This time we will define a function, as we will need this code later to check symmetry of a matrix.

def transpose(a):
    nrow = len(a)
    ncol = len(a[0])
 
    result = []
    for j in range(ncol):
        row = [0] * nrow
        result.append(row[:])
 
    for i in range(nrow):
        for j in range(ncol):
            result[j][i] = a[i][j]
    return result

Let us check how it works for our example.

a
[[1, 2, 3], [4, 5, 6]]
transpose(a)
[[1, 4], [2, 5], [3, 6]]
A
array([[1, 2, 3],
       [4, 5, 6]])
A.T
array([[1, 4],
       [2, 5],
       [3, 6]])
a
     [,1] [,2] [,3]
[1,]    1    2    3
[2,]    4    5    6
t(a)
     [,1] [,2]
[1,]    1    4
[2,]    2    5
[3,]    3    6

1.2 Symmetric matrices

\[ \begin{eqnarray*} A = A^T = \left(\begin{matrix} 1 & 8 & 7 \\ 8 & 5 & 3 \\ 7 & 3 & 5 \end{matrix} \right) \\ \end{eqnarray*} \]

We will use our transpose() function.

a
[[1, 8, 7], [8, 5, 3], [7, 3, 5]]
transpose(a)
[[1, 8, 7], [8, 5, 3], [7, 3, 5]]
a == transpose(a)
True
A 
array([[1, 8, 7],
       [8, 5, 3],
       [7, 3, 5]])
A.T
array([[1, 8, 7],
       [8, 5, 3],
       [7, 3, 5]])
numpy.all(A == A.T) 
True
a
     [,1] [,2] [,3]
[1,]    1    8    7
[2,]    8    5    3
[3,]    7    3    5
t(a)
     [,1] [,2] [,3]
[1,]    1    8    7
[2,]    8    5    3
[3,]    7    3    5
all(a==t(a))
[1] TRUE

1.3 Skew-symmetric matrices

\[ \begin{eqnarray*} A = -A^T = \left(\begin{matrix} 0 & -8 & 7 \\ 8 & 0 & -3 \\ -7 & 3 & 0 \end{matrix} \right) \\ \end{eqnarray*} \]

At first, we will define our negate() function.

def negate(a):
    nrow = len(a)
    ncol = len(a[0])
    result = []
    for i in range(nrow):
        row = []
        for j in range(ncol):
            row.append(-a[i][j])
        result.append(row)
    return result

Let us check how our transpose() and negate() functions work for our example.

a
[[0, -8, 7], [8, 0, -3], [-7, 3, 0]]
at = transpose(a); at
[[0, 8, -7], [-8, 0, 3], [7, -3, 0]]
atn = negate(at); atn 
[[0, -8, 7], [8, 0, -3], [-7, 3, 0]]
a == atn 
True
A 
array([[ 0, -8,  7],
       [ 8,  0, -3],
       [-7,  3,  0]])
-A.T
array([[ 0, -8,  7],
       [ 8,  0, -3],
       [-7,  3,  0]])
numpy.all(A == -A.T) 
True
a
     [,1] [,2] [,3]
[1,]    0   -8    7
[2,]    8    0   -3
[3,]   -7    3    0
-t(a)
     [,1] [,2] [,3]
[1,]    0   -8    7
[2,]    8    0   -3
[3,]   -7    3    0
all(a==-t(a))
[1] TRUE

2 Diagonalization

2.1 Diagonal matrices

\[ \begin{equation*} D_1 = \left(\begin{matrix} -1 & 0 & 0 \\ 0 & 7 & 0 \\ 0 & 0 & 5\end{matrix}\right) \\ \end{equation*} \]

x = [-1, 7, 5]

nrow = ncol = len(x)
result = []
row = [0] * ncol
for i in range(nrow):
    result.append(row[:]) 
    result[i][i] = x[i]

result
[[-1, 0, 0], [0, 7, 0], [0, 0, 5]]
numpy.diag([-1, 7, 5])
array([[-1,  0,  0],
       [ 0,  7,  0],
       [ 0,  0,  5]])
diag(c(-1, 7, 5))
     [,1] [,2] [,3]
[1,]   -1    0    0
[2,]    0    7    0
[3,]    0    0    5

2.2 Eigenvalues and eigenvectors

For

\[ \begin{equation*} A = \left( \begin{matrix} -3 & 2 & 4 \\ -2 & 1 & 4 \\ -2 & 2 & 3 \end{matrix} \right), \quad X = \left( \begin{matrix} 1\\1\\1 \end{matrix} \right). \end{equation*} \]

we can verify \(AX = 3X\).

a
[[-3, 2, 4], [-2, 1, 4], [-2, 2, 3]]
x
[1, 1, 1]

a * x =

ncol = len(x)
nrow = len(a)
result_1= [0] * ncol

for i in range(nrow):
    for j in range(ncol):
        result_1[i] += a[i][j] * x[j]
result_1
[3, 3, 3]

3 * x =

result_2 = [3*n for n in x]; result_2
[3, 3, 3]
result_1 == result_2
True
A
array([[-3,  2,  4],
       [-2,  1,  4],
       [-2,  2,  3]])
X
array([1, 1, 1])
numpy.dot(A, X) 
array([3, 3, 3])
3 * X
array([3, 3, 3])
numpy.all(numpy.dot(A, X) == 3 * X)
True
a
     [,1] [,2] [,3]
[1,]   -3    2    4
[2,]   -2    1    4
[3,]   -2    2    3
x
     [,1]
[1,]    1
[2,]    1
[3,]    1
a %*% x 
     [,1]
[1,]    3
[2,]    3
[3,]    3
3 * x
     [,1]
[1,]    3
[2,]    3
[3,]    3
all(a %*% x  == 3 * x)
[1] TRUE

2.3 Matrices that are diagonalizable

\[ \begin{equation*} A = \left(\begin{matrix}1 & 2 \\ 3 & 2 \end{matrix}\right) = PDP^{-1} = \left(\begin{matrix}2 & 1 \\ 3 & -1 \end{matrix}\right) \left(\begin{matrix}4 & 0 \\ 0 & -1 \end{matrix}\right) \left(\begin{matrix}1/5 & 1/5 \\ 3/5 & -2/5 \end{matrix}\right) . \end{equation*} \]

a
[[1, 2], [3, 2]]
p
[[2, 1], [3, -1]]
d
[[4, 0], [0, -1]]
pi
[[0.2, 0.2], [0.6, -0.4]]

We use our multiply() function from Session 2 with rounding added.

result = multiply(multiply(p, d), pi); result
[[1.0, 2.0], [3.0, 2.0]]
a == result 
True
A
array([[1, 2],
       [3, 2]])
P
array([[ 2,  1],
       [ 3, -1]])
D
array([[ 4,  0],
       [ 0, -1]])
Pi = numpy.linalg.inv(P); Pi
array([[ 0.2,  0.2],
       [ 0.6, -0.4]])
result = numpy.dot(numpy.dot(P, D), Pi) 
result 
array([[1., 2.],
       [3., 2.]])

As elements of result are real numbers and elements of A are whole numbers, we will use the allclose() function to check if two arrays are element-wise equal within a tolerance.

numpy.allclose(A, result)
True
library(matlib) 
a
     [,1] [,2]
[1,]    1    2
[2,]    3    2
p
     [,1] [,2]
[1,]    2    1
[2,]    3   -1
d
     [,1] [,2]
[1,]    4    0
[2,]    0   -1
pi = inv(p)
p %*% d %*% pi
     [,1] [,2]
[1,]    1    2
[2,]    3    2
all(a == p %*% d %*% pi)
[1] FALSE

2.4 Orthogonal matrices

For

\[ \begin{eqnarray*} R_1 = \left(\begin{matrix}0.6 & -0.8 \\ 0.8 & 0.6 \end{matrix} \right) \\ \end{eqnarray*} \]

we can verify \(RR^T=I\).

We use our multiply() and transpose() functions.

r1
[[0.6, -0.8], [0.8, 0.6]]
multiply(r1, transpose(r1))
[[1.0, 0.0], [0.0, 1.0]]
R1
array([[ 0.6, -0.8],
       [ 0.8,  0.6]])
numpy.dot(R1, R1.T) 
array([[ 1.00000000e+00, -2.66453526e-17],
       [-2.66453526e-17,  1.00000000e+00]])
numpy.around(numpy.dot(R1, R1.T), decimals=15) 
array([[ 1., -0.],
       [-0.,  1.]])
r1
     [,1] [,2]
[1,]  0.6 -0.8
[2,]  0.8  0.6
r1 %*% t(r1)
     [,1] [,2]
[1,]    1    0
[2,]    0    1