##Problem Set 1 ##(1) Show that ATA≠AAT in general. (Proof and demonstration) ##On page 180 of the book, the author states that matrix multiplication is not commutative. Therefore ATA≠AAT.
##If matrix multiplication were commutative then it means that switching the order of multiplcation should ALWAYS give the same answer.
##The demonstration below shows that a 2 x 2 matrix A gives 2 different answers for ATA and AAT.
A = matrix(c(3,3,5,6), nrow=2, ncol=2)
A_T <- t(A)
A
## [,1] [,2]
## [1,] 3 5
## [2,] 3 6
A_T
## [,1] [,2]
## [1,] 3 3
## [2,] 5 6
#ATA results in:
A_T %*% A
## [,1] [,2]
## [1,] 18 33
## [2,] 33 61
AA^T
A %*% A_T
## [,1] [,2]
## [1,] 34 39
## [2,] 39 45
ATA=AAT AA=AA
On page 167, it states that Suppose that A is a symmetric matrix. Then A is square.
Below is a demonstratation that shows ATA=AAT for a symmetric matrix A.
A = matrix(c(2,3,3,1), nrow=2, ncol=2)
A_T <- t(A)
A
## [,1] [,2]
## [1,] 2 3
## [2,] 3 1
A_T
## [,1] [,2]
## [1,] 2 3
## [2,] 3 1
A_T %*% A
## [,1] [,2]
## [1,] 13 9
## [2,] 9 10
A %*% A_T
## [,1] [,2]
## [1,] 13 9
## [2,] 9 10
Problem Set 2 Write an R function to factorize a square matrix A into LU or LDU, whichever you prefer.
Define some functions to help with the LU decomposition.
#Selects which row to multiple to zero out a specific position
select_row <- function(A, ref_row, ref_col){
for(i in 1:nrow(A)){
#print(paste("looking at row: ", i))
#Case when ref_row and current row are the same or the value at A[i, ref_col] is zero.
#We need to look for another row to perform the row operations with.
if(i==ref_row | A[i,ref_col]==0) next #skip to next iteration
if(is_prevColZero(A, i, ref_col)==TRUE){
#we found a row we can use for the row operations
return(i)
}
}
return(-1)#we did not find a row we can use
}
#Returns true if column is 1 or if A[i,1] through A[i,j-1] are all zeros.
is_prevColZero <- function(A, ref_row, ref_col){
#Case: when col is 1, no need to check for zero
#by default we are turning TRUE in this case since it is not relevant
if(ref_col==1)
return(TRUE)
for(j in 1:(ref_col-1)){
#print(paste("checking ", A[ref_row, j]))
if(A[ref_row, j] != 0)
return(FALSE)
}
return(TRUE)
}
#Returns the multiplier to be used to generate a zero
get_multiplier <- function(A, selected_row, ref_row, ref_col){
c <- A[ref_row, ref_col]
mult <- (-1)* A[ref_row, ref_col] * (1/A[selected_row, ref_col])
return(mult)
}
#Returns the LU decomposition as a list.
#item 1 in the list is "L"
#item 2 in the list is "U"
factor_LU <- function(A,n){
U <- A
L <- diag(n)
for(i in 1:nrow(U)){
#print(paste("row: ", i))
for(j in 1:i){
if(j < i){
ref_number = U[i,j]
ref_row = i
ref_col = j
#print(paste("ref_row: ",ref_row, ", ref_col: ", ref_col, ", with ref_number: ", ref_number))
selected_row <- select_row(U,ref_row,ref_col)
if(selected_row < 0){
print("This function is unable to process this matrix.")
break
}
mult <- get_multiplier(U,selected_row, ref_row, ref_col)
#print(mult)
U[ref_row, ] <- mult * U[selected_row,] + U[ref_row, ]
L[ref_row, ref_col] <- (-1) * mult
}
}
}
return(list(L,U))
}
LU decomposition for this given 3x3 matrix A.
n = 3
A <- matrix(c(1,4,-3,-2,8,5,3,4,7),nrow=n,ncol=n,byrow = TRUE)
#A <- matrix(c(1,1,-1,1,-2,3,2,3,1),nrow=n,ncol=n,byrow = TRUE)
LU_decomp <- factor_LU(A,n)
L <- LU_decomp[[1]]
U <- LU_decomp[[2]]
A
## [,1] [,2] [,3]
## [1,] 1 4 -3
## [2,] -2 8 5
## [3,] 3 4 7
L
## [,1] [,2] [,3]
## [1,] 1 0.0 0
## [2,] -2 1.0 0
## [3,] 3 -0.5 1
U
## [,1] [,2] [,3]
## [1,] 1 4 -3.0
## [2,] 0 16 -1.0
## [3,] 0 0 15.5
L %*% U
## [,1] [,2] [,3]
## [1,] 1 4 -3
## [2,] -2 8 5
## [3,] 3 4 7