Section 18 Matrix Multiplication
18.1 Multiplying matrices
First, let’s define a few matrices. We use a trick here. By putting the assignment in parentheses, it both assigns the matrix and displays it.
(A = cbind(c(1,2,3),c(4,5,6),c(1,1,-1)))## [,1] [,2] [,3]
## [1,] 1 4 1
## [2,] 2 5 1
## [3,] 3 6 -1(B = cbind(c(1,-1,1),c(1,1,1),c(0,2,1)))## [,1] [,2] [,3]
## [1,] 1 1 0
## [2,] -1 1 2
## [3,] 1 1 1(C = cbind(c(2,1,1),c(1,0,1),c(1,-3,1),c(3,2,1)))## [,1] [,2] [,3] [,4]
## [1,] 2 1 1 3
## [2,] 1 0 -3 2
## [3,] 1 1 1 1Multiply these matrices to get \(AB\) using the %*% command. As seen here:
A %*% B## [,1] [,2] [,3]
## [1,] -2 6 9
## [2,] -2 8 11
## [3,] -4 8 11Note that \(BA\) is not the same as \(AB\).
B %*% A## [,1] [,2] [,3]
## [1,] 3 9 2
## [2,] 7 13 -2
## [3,] 6 15 1That is, matrix multiplication is not commutative. It can be the case that \(AB\) equals \(BA\) but most if the time these are not equal.
We can multiply \(BC\) to get
B %*% C## [,1] [,2] [,3] [,4]
## [1,] 3 1 -2 5
## [2,] 1 1 -2 1
## [3,] 4 2 -1 6But \(CB\) does not make sense, since \(C\) is 3 x 4 and \(B\) is 3 x 3. The inner dimensions to not match. R tells us that the matrices are non-conformable.
C %*% B## Error in C %*% B: non-conformable arguments18.2 Transpose
The transpose of a matrix is computed by t(A). For example
t(A)## [,1] [,2] [,3]
## [1,] 1 2 3
## [2,] 4 5 6
## [3,] 1 1 -1t(B)## [,1] [,2] [,3]
## [1,] 1 -1 1
## [2,] 1 1 1
## [3,] 0 2 1Note that the transpose is order reversing. That is \[ (AB)^T = B^T A^T. \] We can see that here
## [,1] [,2] [,3]
## [1,] -2 -2 -4
## [2,] 6 8 8
## [3,] 9 11 11## [,1] [,2] [,3]
## [1,] -2 -2 -4
## [2,] 6 8 8
## [3,] 9 11 1118.3 Identity Matrix
The command diag(n) gives the n x n identity matrix. This is denoted \(I_n\). For example, below are \(I_3\), \(I_4\), and \(I_5.\)
diag(3)## [,1] [,2] [,3]
## [1,] 1 0 0
## [2,] 0 1 0
## [3,] 0 0 1diag(4)## [,1] [,2] [,3] [,4]
## [1,] 1 0 0 0
## [2,] 0 1 0 0
## [3,] 0 0 1 0
## [4,] 0 0 0 1diag(5)## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 0 0 0 0
## [2,] 0 1 0 0 0
## [3,] 0 0 1 0 0
## [4,] 0 0 0 1 0
## [5,] 0 0 0 0 1The identity matrix has the same effect as multiplyng by 1. For example,
diag(3) %*% A## [,1] [,2] [,3]
## [1,] 1 4 1
## [2,] 2 5 1
## [3,] 3 6 -1A %*% diag(3)## [,1] [,2] [,3]
## [1,] 1 4 1
## [2,] 2 5 1
## [3,] 3 6 -118.4 Writing Loops
Work on problem 3.8 in the homework. In this problem, you will write out a loop in R to apply a matrix over and over again. Here is an example for you to use to compare with.
First, define a vector \(v\) and a matrix \(A\) and multiply \(Av\).
v = c(10,90)
A = cbind(c(.9,.1),c(.5,.5))
A %*% v## [,1]
## [1,] 54
## [2,] 46If we want to appy \(A\) to \(v\) again and again, we can keep multiplying by \(A\):
A %*% A %*% v## [,1]
## [1,] 71.6
## [2,] 28.4A %*% A %*% A %*% v## [,1]
## [1,] 78.64
## [2,] 21.36A %*% A %*% A %*% A %*% v## [,1]
## [1,] 81.456
## [2,] 18.544It is better to write a loop to do something like this. The following loop will multiply by A 10 times. Try changing the 10 to 4 (and compare with above). Try changing the 10 to 100.
N = 10 # the number of times that we will multiply by A
w = c(10,90) # the starting value
for (i in 1:N) { # loop as the index i goes from 1 to N
w = A %*% w # multiply Aw and replace w with this new vector
}
w # show the final value of w## [,1]
## [1,] 83.32564
## [2,] 16.6743618.5 Matrix Inverses
Give a square matrix
## [,1] [,2] [,3]
## [1,] 10 -6 1
## [2,] -2 1 0
## [3,] -7 5 -1You find its inverse using solve:
## [,1] [,2] [,3]
## [1,] 1 1 1
## [2,] 2 3 2
## [3,] 3 8 2## [,1] [,2] [,3]
## [1,] 1.000000e+00 0.000000e+00 0
## [2,] 1.776357e-15 1.000000e+00 0
## [3,] 3.552714e-15 1.776357e-15 1Note that if you look closely at the last matrix, it is within round off error of the identity. For example, the number 1.776357 e-15 is 0.000000000000001776357. To make it look more like the identity, you can use the command zapsmall, which converts very small numbers like this to 0.
## [,1] [,2] [,3]
## [1,] 1 0 0
## [2,] 0 1 0
## [3,] 0 0 1What happens if you try to invert a matrix that is not invertible. For example in the following matrix, the third column is the sum of the first two, so the columns are linearly dependent. If we try to inverti it
A = cbind(c(10,-2,-7),c(-6,1,5),c(4,-1,-2))
A## [,1] [,2] [,3]
## [1,] 10 -6 4
## [2,] -2 1 -1
## [3,] -7 5 -2solve(A)## Error in solve.default(A): system is computationally singular: reciprocal condition number = 5.84328e-18We get the error that the matrix is singular. This is another term for non-invertible.