On page 180 of the book, the author states that matrix multiplication is not commutative. Therefore \(A^TA \ne AA^T\).
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 \(A^TA\) and \(AA^T\).
A = matrix(c(1,2,4,3), nrow=2, ncol=2)
A_T <- t(A)Matrix A:
A## [,1] [,2]
## [1,] 1 4
## [2,] 2 3Matrix \(A^T\):
A_T## [,1] [,2]
## [1,] 1 2
## [2,] 4 3\(A^TA\) results in:
A_T %*% A## [,1] [,2]
## [1,] 5 10
## [2,] 10 25\(AA^T\) resutls in:
A %*% A_T## [,1] [,2]
## [1,] 17 14
## [2,] 14 13As you can see, the results are NOT the same. This shows that \(A^TA \ne AA^T\).
The definition of a symmetric matrix on page 166 says that matrix A is symmetric if \(A = A^T\). So this means that when A is a symmetric matrix:
\[A^TA = AA^T\] \[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 \(A^TA = AA^T\) for a symmetric matrix A.
A = matrix(c(1,3,3,2), nrow=2, ncol=2)
A_T <- t(A)Matrix A:
A## [,1] [,2]
## [1,] 1 3
## [2,] 3 2Matrix \(A^T\):
A_T## [,1] [,2]
## [1,] 1 3
## [2,] 3 2\(A^TA\) results in:
A_T %*% A## [,1] [,2]
## [1,] 10 9
## [2,] 9 13\(AA^T\) resutls in:
A %*% A_T## [,1] [,2]
## [1,] 10 9
## [2,] 9 13Write an R function to factorize a square matrix A into LU or LDU, whichever you prefer.
#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))
}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]]After running the code above, here are the values for L, U for Matrix A:
A Matrix is:
A## [,1] [,2] [,3]
## [1,] 1 4 -3
## [2,] -2 8 5
## [3,] 3 4 7L Matrix is:
L## [,1] [,2] [,3]
## [1,] 1 0.0 0
## [2,] -2 1.0 0
## [3,] 3 -0.5 1U Matrix is:
U## [,1] [,2] [,3]
## [1,] 1 4 -3.0
## [2,] 0 16 -1.0
## [3,] 0 0 15.5LU Matrix is:
L %*% U## [,1] [,2] [,3]
## [1,] 1 4 -3
## [2,] -2 8 5
## [3,] 3 4 7