PSY 613, Spring 2015: Lab #1 - Matrix Algebra

Note: For Problem Set #1, you can do your work by hand on paper, and then just scan it and turn in the pdf. For all future assignments, we prefer editable document files (like .doc or .docx), so it’s easier for us to provide feedback on your work.

Notation and Terminology:

\(a_{rc}\) refers to the entry at row \(r\) and column \(c\) in matrix \(a\).
“Order” refers to the dimensions of a matrix: the number of rows and the number of columns \((r,c)\).

Types of Matrices

  1. Rectangular: \(r \neq c\)
    \[ a = \begin{bmatrix} 1 & 5 & 7 \\ 2 & 1 & 2 \end{bmatrix} \]
  2. Square: \(r = c\)
    \[ a = \begin{bmatrix} 2 & 9\\ 4 & 6 \end{bmatrix} \]
  3. Vector: A matrix where the row or column (not both) is 1.
    \[ a = \begin{bmatrix} 5 & 6 & 9 \end{bmatrix} \]
  4. Diagonal: A square matrix where all of the elements equal zero except for those making up the principal diagonal.
    \[ a = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 4 & 0 \\ 0 & 0 & 3 \end{bmatrix} \]
  5. Identity: The diagonal matrix with 1s along the principal diagonal.
    \[ a = \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix} \]
  6. Null: A matrix that consists entirely of 0s.
    \[ 0 = \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix} \]

Arithmetic operations on the elements of matrices:

Equality

Two matrices are equal if (1) they are of the same order and (2) all elements in corresponding positions (e.g., \(a_{11}\) & \(b_{11}\)) are the same.

Which matrices are equal?

\[ a = \begin{bmatrix} 1 & 2 & 3 \\ 1 & 2 & 3 \end{bmatrix}\ b = \begin{bmatrix} 1 & 1 \\ 2 & 2 \\ 3 & 3 \end{bmatrix}\ c = \begin{bmatrix} 1 & 1 \\ 2 & 4 \\ 3 & 5 \end{bmatrix}\ d = \begin{bmatrix} 1 & 1 \\ 2 & 2 \\ 3 & 3 \end{bmatrix} \]

We can use R for matrices too!

#Create these matrices using the matrix function.

#byrow = T tells the function to read the first 3 numbers as row 1,
#  and the second three numbers as row 2. Set byrow = F to see what
#  happens.  


a <- matrix(data = c(1,2,3, 1,2,3), nrow = 2, ncol = 3, byrow = T) 
print(a)
##      [,1] [,2] [,3]
## [1,]    1    2    3
## [2,]    1    2    3
b <- matrix(data = c(1,1, 2,2, 3,3), nrow = 3, ncol = 2, byrow = T) 
print(b)
##      [,1] [,2]
## [1,]    1    1
## [2,]    2    2
## [3,]    3    3
c <- matrix(data = c(1,1, 2,4, 3,5), nrow = 3, ncol = 2, byrow = T) 
print(c)
##      [,1] [,2]
## [1,]    1    1
## [2,]    2    4
## [3,]    3    5
d <- matrix(data = c(1,1, 2,2, 3,3), nrow = 3, ncol = 2, byrow = T) 
print(d)
##      [,1] [,2]
## [1,]    1    1
## [2,]    2    2
## [3,]    3    3
#First, check if the order is the same using dim(), which returns the dimensions of the matrix (rows, columns).
# we get two results for each test, 1) same number of rows? 2) same number of columns?
dim(a)==dim(b)
## [1] FALSE FALSE
dim(a)==dim(c)
## [1] FALSE FALSE
dim(a)==dim(d)
## [1] FALSE FALSE
# Nope. Now let's check b against c and d.
dim(b)==dim(c)
## [1] TRUE TRUE
dim(b)==dim(d)
## [1] TRUE TRUE
#Dimensions of b and d are the same, so let's test equality of elements.
b==c
##      [,1]  [,2]
## [1,] TRUE  TRUE
## [2,] TRUE FALSE
## [3,] TRUE FALSE
b==d
##      [,1] [,2]
## [1,] TRUE TRUE
## [2,] TRUE TRUE
## [3,] TRUE TRUE
#Bingo!

#You might also want to play around with identical() and all.equal().

Transpose

When transposing a matrix, simply exchange the row and column for each element in the matrix.

Example:

  1. \[ a = \begin{bmatrix} 1 & 2 \\ 5 & 1 \\ 7 & 2 \end{bmatrix} = \begin{bmatrix} 1_{11} & 2_{12} \\ 5_{21} & 1_{22} \\ 7_{31} & 2_{32} \end{bmatrix} \]
  2. So,
  • \(1_{11} \rightarrow 1_{11}\)
  • \(5_{21} \rightarrow 5_{12}\)
  • \(7_{31} \rightarrow 7_{13}\)
  • \(2_{12} \rightarrow 2_{21}\)
  • \(1_{22} \rightarrow 1_{22}\)
  • \(2_{32} \rightarrow 2_{23}\)
  1. Since the order of the above matrix is 3x2, what will the order of its transpose be?
    \[ a' = \begin{bmatrix} 1 & 5 & 7 \\ 2 & 1 & 2 \end{bmatrix} \]

To transpose in r, just use the function t().

#We'll use our matrix from before, c
print(c)
##      [,1] [,2]
## [1,]    1    1
## [2,]    2    4
## [3,]    3    5
t(c)
##      [,1] [,2] [,3]
## [1,]    1    2    3
## [2,]    1    4    5

Trace

The sum of the diagonal elements in a square matrix

Example:

\[ a= \begin{bmatrix} 2 & 9 \\ 4 & 6 \end{bmatrix} \] \[TR(a)=2+6=8\]

a <- matrix(data = c(2,9, 4,6), nrow = 2, ncol = 2, byrow = T) 
print(a)
##      [,1] [,2]
## [1,]    2    9
## [2,]    4    6
sum(diag(a))
## [1] 8

Addition/subtraction

You can only add/subtract matrices if they have the same _______.

\[ a=\begin{bmatrix} 6 & 1 \\ 2 & 10 \end{bmatrix}, b=\begin{bmatrix} 2 & 1 \\ 1 & 6 \end{bmatrix}, c=\begin{bmatrix} 4 & 1 \\ 3 & 2 \end{bmatrix} \]

  1. \(a+b=?\)
  2. Is the following statement true? \(a + b = b + a\)
  3. Is the following statement true? \((a + b) + c = a + (b + c)\)
  4. \(a - b = ?\)
  5. Is the following statement true? \(a - b = b - a\)
  6. Is the following statement true? \((a - b) - c = a - (b - c)\)

Try the above problems in r using +, -, and ==.

a <- matrix(data = c(6,1, 2,10), nrow = 2, ncol = 2, byrow = T) 
b <- matrix(data = c(2,1, 1,6), nrow = 2, ncol = 2, byrow = T) 
c <- matrix(data = c(4,1, 3,2), nrow = 2, ncol = 2, byrow = T) 

Multiplication

1:

\[ a=\begin{bmatrix} 1 & 3 \\ 2 & 1 \\ 4 & 6 \end{bmatrix}, b=\begin{bmatrix} 1 & 5 & 8 & 3 \\ 4 & 2 & 6 & 4 \end{bmatrix}\]

  1. Is \(ab\) conformable?
  2. Is \(ba\) conformable?
  3. For \(ab\), what will the order of the resulting matrix be?
[ab=

\

ab=

]

2:

\[ a=\begin{bmatrix} 1 & 2 \\ 3 & 2 \\ 4 & 1 \end{bmatrix}, b=\begin{bmatrix} 2 & 3 & 4 \\ 1 & 2 & 4 \end{bmatrix}\]

  1. Is \(ab\) conformable?
  2. is \(ba\) conformable?
  3. Does \(ab = ba\)? In other words, does the matrix product possess the commutative property?

Note: the matrix product does possess the associative property: \(a(bc)=(ab)c\).

Matrix multiplication in r: use %*%.

a <- matrix(data = c(1,2, 3,2, 4,1), nrow = 3, ncol = 2, byrow = T) 
b <- matrix(data = c(2,3,4, 1,2,4), nrow = 2, ncol = 3, byrow = T)
a %*% b
##      [,1] [,2] [,3]
## [1,]    4    7   12
## [2,]    8   13   20
## [3,]    9   14   20

3:

\[ c=\begin{bmatrix} 2 & 5 \\ 4 & 3 \end{bmatrix} \]

  1. If c is multiplied by its identity matrix (\(I_{2}\)), what will the resulting matrix be?
  2. Does \(cI_{2} = I_{2}c\)?

Note: for your homework, you can’t just say \(cI_{2} = c\), you need to go through the steps of the multiplication to prove it.

Use diag(1,r,c) to create the identity matrix.

c <- matrix(data = c(2,5, 4,3), nrow = 2, ncol = 2, byrow = T)
mI <- diag(1,dim(c)[1],dim(c)[2])
print(c)
##      [,1] [,2]
## [1,]    2    5
## [2,]    4    3
print(mI)
##      [,1] [,2]
## [1,]    1    0
## [2,]    0    1
c %*% mI
##      [,1] [,2]
## [1,]    2    5
## [2,]    4    3

Row echelon form

Tells you the number of linearly independent rows in a matrix (which represents the row rank of a matrix). For example, a row rank of 2 means that 2 rows are independent.

\[a=\begin{bmatrix} 1 & 2 & 5 & 4 \\ 3 & 2 & 4 & 3 \\ 6 & 4 & 8 & 6 \\ 2 & 2 & 8 & 1 \end{bmatrix} \]

  1. What is the row echelon form of \(a\)?
    Note: When you’re trying to make a 1, multiply things to the row you’re working on. When you’re trying to make 0’s, multiply things to the other rows, then add them to the row you’re working on.
    Step 1: Make the entry in R1 and C1 (1,1) into 1.
    Already done!

Step 2: Make the remaining elements in C1 into 0.
- For R2:
- Add -3 times R1 to R2. In other words, multiply R1 by -3, and then add those elements to the corresponding elements in R2. - R1 x -3: -3, -6, -15, -12 (Note: this is only a temporary change and does NOT go into your new matrix.) - Above + R2: 0, -4, -11, -9 (This DOES go into your new matrix.) - For R3: - Add -6 times R1 to R3. - R1 x -6: -6, -12, -30, -24 - Above + R3: 0, -4, -11, -9 - For R4: - Add -2 times R1 to R4. - R1 x -2: -2, -4, -10, -8 - Above + R4: 0, -2, -2, -7

\[a=\begin{bmatrix} 1 & 2 & 5 & 4 \\ 0 & -4 & -11 & -9 \\ 0 & -8 & -22 & -18 \\ 0 & -2 & -2 & -7 \end{bmatrix} \]

Step 3: Make the entry in (2,2) into 1. - Multiply R2 by -1/4.

\[a=\begin{bmatrix} 1 & 2 & 5 & 4 \\ 0 & 1 & 11/4 & 9/4 \\ 0 & -8 & -22 & -18 \\ 0 & -2 & -2 & -7 \end{bmatrix} \]

Step 4: Make the remaining elements in C2 into 0. - For R1: - Add -2 times R2 to R1 - R2 x -2: 0, -2, -11/2, -9/2 - Above + R1: 1, 0, -1/2, -1/2 - For R3: - Add 8 times R2 to R3 - R2 x 8: 0, 8, 22, 18 - Above + R3: 0, 0, 0, 0
- For R4: - Add 2 times R2 to R4 - R2 x 2: 0, 2, 11/2, 9/2 - Above + R4: 0, 0, 7/2, -5/2

\[a=\begin{bmatrix} 1 & 0 & -1/2 & -1/2 \\ 0 & 1 & 11/4 & 9/4 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 7/2 & -5/2 \end{bmatrix} \]

Uh-oh. There’s a 0 at (3,3), which means there’s nothing we could multiply to that row to get a 1 there!! This indicates that there’s at least one row in this matrix that is not indpendent. Since we haven’t really finished testing R4, we don’t know whether there are 3 indpendent rows here, or only 2. To be able to test R4 fully, swap the zeros row down to the end. Weird, huh? But remember, all we’re testing is how many rows are independent; the order of the rows actually doesn’t matter at all, so it’s fine to swap them around as much as you like.

Step 5: Swap R3 and R4 to move the zeros to the bottom.

\[a=\begin{bmatrix} 1 & 0 & -1/2 & -1/2 \\ 0 & 1 & 11/4 & 9/4 \\ 0 & 0 & 7/2 & -5/2 \\ 0 & 0 & 0 & 0 \end{bmatrix} \]

Step 6: Make the entry in (3,3) into 1. - Multiply R3 by 2/7.

\[a=\begin{bmatrix} 1 & 0 & -1/2 & -1/2 \\ 0 & 1 & 11/4 & 9/4 \\ 0 & 0 & 1 & -10/14 \\ 0 & 0 & 0 & 0 \end{bmatrix} \]

Step 7: Make the remaining elements in C3 into 0. - For R1: - Add 1/2 times R3 to R1 - R3 x 1/2: 0, 0, 1/2, -5/14 - Above + R1: 1, 0, 0, -6/7 - For R2: - Add -11/4 times R3 to R2 - R3 x -11/4: 0, 0, -11/4, 55/28 - Above + R2: 0, 1, 0, 59/14 - For R4: - Already done!

\[a=\begin{bmatrix} 1 & 0 & 0 & -6/7 \\ 0 & 1 & 0 & 59/14 \\ 0 & 0 & 1 & -10/14 \\ 0 & 0 & 0 & 0 \end{bmatrix} \]

Note that this process is designed to tell you how many independent rows there are, not how individiual rows are related to each other. For example, you don’t know from the above calculations that row three from the original matrix is the one that’s not independent (and in fact that’s not really the case - you can see from carefully examining the original matrix that R3 = R2*2, so there’s no unique information in R3, but you could just as easily say that R2 = R3/2, so there’s no unique information in R2). By definition, independence or non-independence of a given row is relative to the other rows.

Not surprisingly, there is an R package that provides a function that produces the row echelon form of matrix. If you don’t already have the pracma package installed, install it first using install.packages(“pracma”)

a <- matrix(data = c(1,2,5,4, 3,2,4,3, 6,4,8,6,  2,2,8,1), nrow = 4, ncol = 4, byrow = T)

library(pracma)
rref(a)
##      [,1] [,2] [,3]       [,4]
## [1,]    1    0    0 -0.8571429
## [2,]    0    1    0  4.2142857
## [3,]    0    0    1 -0.7142857
## [4,]    0    0    0  0.0000000
# if you want to display the numbers as fractions instead of decimals, use the fractions() function from the MASS package
library(MASS)
fractions(rref(a))
##      [,1]  [,2]  [,3]  [,4] 
## [1,]     1     0     0  -6/7
## [2,]     0     1     0 59/14
## [3,]     0     0     1  -5/7
## [4,]     0     0     0     0

Inverse

A matrix (\(a\)) multiplied by its inverse (\(a^-1\)) equals the identity matrix.
For certain matrices, you can calculate the inverse by hand.
1. A 2x2 matrix

\[a=\begin{bmatrix} a & b \\ c & d \end{bmatrix} \] \[a^-1=\begin{bmatrix} d/(a*d-b*c) & -b/(a*d-b*c) \\ -c/(a*d-b*c) & a/(a*d-b*c) \end{bmatrix} \]

Example: What does \(a^-1\) equal? \[a=\begin{bmatrix} 3 & 2 \\ 2 & 4 \end{bmatrix} \]

In r, you can get the inverse of a square matrix with the function solve()

?solve
a <- matrix(data = c(3,2,2,4), nrow = 2, byrow = T)
solve(a)
##       [,1]   [,2]
## [1,]  0.50 -0.250
## [2,] -0.25  0.375

Verify that \(aa^-1 = I_{2}\)

a <- matrix(data = c(3,2,2,4), nrow = 2, byrow = T)
a_inv <- solve(a)
a %*% a_inv
##              [,1] [,2]
## [1,] 1.000000e+00    0
## [2,] 1.110223e-16    1
# Note that sometimes when R does matrix algebra there's a little weirdness with rounding. If you expected to see a 0 and instead you get a super small number (like 1e-16), assume it really is just zero. 
  1. Diagonal matrix: A square matrix where all of the elements are equal to 0 except the elements that make up the principal diagonal.

\[b=\begin{bmatrix} 2 & 0 & 0 & 0 \\ 0 & 6 & 0 & 0 \\ 0 & 0 & 8 & 0 \\ 0 & 0 & 0 & 4 \end{bmatrix} \] \[b^-1=\begin{bmatrix} 1/2 & 0 & 0 & 0 \\ 0 & 1/6 & 0 & 0 \\ 0 & 0 & 1/8 & 0 \\ 0 & 0 & 0 & 1/4 \end{bmatrix} \]

Verify that \(bb^-1 = I_{4}\)

b <- matrix(data = c(2,0,0,0,0,6,0,0,0,0,8,0,0,0,0,4), nrow = 4, byrow = T)
b_inv <- solve(b)
b %*% b_inv
##      [,1] [,2] [,3] [,4]
## [1,]    1    0    0    0
## [2,]    0    1    0    0
## [3,]    0    0    1    0
## [4,]    0    0    0    1
# Yay!

Check out this guide to matrix algebra in r