Problem Set 1

(1) Show that A^TA != AA^T in general.

The commutative law is usually broken for matrix multiplication, so that AB != BA.

Here is a specific example to illustrate that general principal for matrix A and transposed matrix A:

# Declare matrix A
A <- matrix(c(1,0,1,2,1,0,6,4,0), nrow = 3, byrow = T)
A
##      [,1] [,2] [,3]
## [1,]    1    0    1
## [2,]    2    1    0
## [3,]    6    4    0
# Transpose matrix A, for A^T
t(A)
##      [,1] [,2] [,3]
## [1,]    1    2    6
## [2,]    0    1    4
## [3,]    1    0    0
# Transposed matrix A times matrix A (A^TA)
ATA <- t(A) %*% A
ATA
##      [,1] [,2] [,3]
## [1,]   41   26    1
## [2,]   26   17    0
## [3,]    1    0    1
# Matrix A times transposed matrix A (AA^T)
AAT <- A %*% t(A)
AAT
##      [,1] [,2] [,3]
## [1,]    2    2    6
## [2,]    2    5   16
## [3,]    6   16   52

So, A^TA != AA^T.

ATA == AAT
##       [,1]  [,2]  [,3]
## [1,] FALSE FALSE FALSE
## [2,] FALSE FALSE FALSE
## [3,] FALSE FALSE FALSE

(2) For a special type of square matrix A, we get A^TA = AA^T . Under what conditions could this be true? (Hint: The Identity matrix I is an example of such a matrix).

This could be true with a symmetric matrix, as it is, by definition, equal to its transpose. As such, you are essentially multiplying matrix A by itself both times, so the resulting matrices will be equal.

Here is a specific example using a symmetric matrix:

# Declare matrix A
A <- matrix(c(3,2,4,2,0,2,4,2,3), nrow = 3, byrow = T)
A
##      [,1] [,2] [,3]
## [1,]    3    2    4
## [2,]    2    0    2
## [3,]    4    2    3
# Transpose matrix A, for A^T
t(A)
##      [,1] [,2] [,3]
## [1,]    3    2    4
## [2,]    2    0    2
## [3,]    4    2    3
# Transposed matrix A times matrix A (A^TA)
ATA <- t(A) %*% A
ATA
##      [,1] [,2] [,3]
## [1,]   29   14   28
## [2,]   14    8   14
## [3,]   28   14   29
# Matrix A times transposed matrix A (AA^T)
AAT <- A %*% t(A)
AAT
##      [,1] [,2] [,3]
## [1,]   29   14   28
## [2,]   14    8   14
## [3,]   28   14   29

Here, A^TA = AA^T.

ATA == AAT
##      [,1] [,2] [,3]
## [1,] TRUE TRUE TRUE
## [2,] TRUE TRUE TRUE
## [3,] TRUE TRUE TRUE

Problem Set 2

Write an R function to factorize a square matrix A into LU or LDU, whichever you prefer.

The function factorize checks for the size of a square matrix, and then factorizes it into LU.
#Declare the function
factorize <- function(A) {

#Check size of matrix - first 2x2, then 3x3, then 4x4
if (nrow(A) == 2) {

  #2x2 matrices

  #Subract multipler E21
  E21 <- -A[2,1] / A[1,1]
  L <- matrix(c(1,0,E21,1),nrow=2, byrow = T)
  U <- L %*% A
  #U

  print(A)
  print(L)
  print(U)

  L <- solve(L)

  (L %*% U == A)

} else if (nrow(A) == 3) {

  #3x3 matrices

  #Subract multipler E21
  E21 <- -A[2,1] / A[1,1]
  mE21 <- matrix(c(1,0,0,E21,1,0,0,0,1), nrow=3, byrow = T)
  A2 <- mE21 %*% A
  #A2

  #Subract multipler E31
  E31 <- -A2[3,1] / A2[1,1]
  mE31 <- matrix(c(1,0,0,0,1,0,E31,0,1), nrow=3, byrow = T)
  A3 <- mE31 %*% A2
  #A3

  #Subtract multipler E32
  E32 <- -A3[3,2] / A3[2,2]
  mE32 <- matrix(c(1,0,0,0,1,0,0,E32,1), nrow=3, byrow = T)
  U <- mE32 %*% A3
  #U

  L <- solve(mE21) %*% solve(mE31) %*% solve(mE32)

  A <- L %*% U

  print(A)
  print(L)
  print(U)

  (L %*% U == A)

} else if (nrow(A) == 4) {

  #4x4 matrices
  
  #Subract multipler E21
  E21 <- -A[2,1] / A[1,1]
  mE21 <- matrix(c(1,0,0,0,E21,1,0,0,0,0,1,0,0,0,0,1), nrow=4, byrow = T)
  A2 <- mE21 %*% A
  #A2
  
  #Subract multipler E31
  E31 <- -A2[3,1] / A2[1,1]
  mE31 <- matrix(c(1,0,0,0,0,1,0,0,E31,0,1,0,0,0,0,1), nrow=4, byrow = T)
  A3 <- mE31 %*% A2
  #A3
  
  #Subtract multipler E32
  E32 <- -A3[3,2] / A3[2,2]
  mE32 <- matrix(c(1,0,0,0,0,1,0,0,0,E32,1,0,0,0,0,1), nrow=4, byrow = T)
  A4 <- mE32 %*% A3
  #A4
  
  #Subract multipler E41
  E41 <- -A4[4,1] / A4[1,1]
  mE41 <- matrix(c(1,0,0,0,0,1,0,0,0,0,1,0,E41,0,0,1), nrow=4, byrow = T)
  A5 <- mE41 %*% A4
  #A5
  
  #Subtract multipler E42
  E42 <- -A5[4,2] / A5[2,2]
  mE42 <- matrix(c(1,0,0,0,0,1,0,0,0,0,1,0,0,E42,0,1), nrow=4, byrow = T)
  A6 <- mE42 %*% A5
  #A6
  
  #Subtract multipler E43
  E43 <- -A6[4,3] / A6[3,3]
  mE43 <- matrix(c(1,0,0,0,0,1,0,0,0,0,1,0,0,0,E43,1), nrow=4, byrow = T)
  U <- mE43 %*% A6
  #U
  
  L <- solve(mE21) %*% solve(mE31) %*% solve(mE32) %*% solve(mE41) %*% solve(mE42) %*% solve(mE43)
  
  A <- L %*% U
  
  print(A)
  print(L)
  print(U)
  
  (L %*% U == A)
}
  
} #final function right bracket
Test with 2x2 matrix
A <- matrix(c(2,1,6,8),nrow=2, byrow = T)
factorize(A)
##      [,1] [,2]
## [1,]    2    1
## [2,]    6    8
##      [,1] [,2]
## [1,]    1    0
## [2,]   -3    1
##      [,1] [,2]
## [1,]    2    1
## [2,]    0    5
##      [,1] [,2]
## [1,] TRUE TRUE
## [2,] TRUE TRUE
Test with 3x3 matrix
A <- matrix(c(1,1,2,2,1,0,3,1,1),nrow=3, byrow = T)
factorize(A)
##      [,1] [,2] [,3]
## [1,]    1    1    2
## [2,]    2    1    0
## [3,]    3    1    1
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    2    1    0
## [3,]    3    2    1
##      [,1] [,2] [,3]
## [1,]    1    1    2
## [2,]    0   -1   -4
## [3,]    0    0    3
##      [,1] [,2] [,3]
## [1,] TRUE TRUE TRUE
## [2,] TRUE TRUE TRUE
## [3,] TRUE TRUE TRUE
Test with 4x4 matrix
A <- matrix(c(1,1,1,1,1,2,3,4,1,3,6,10,1,4,10,20),nrow=4, byrow = T)
factorize(A)
##      [,1] [,2] [,3] [,4]
## [1,]    1    1    1    1
## [2,]    1    2    3    4
## [3,]    1    3    6   10
## [4,]    1    4   10   20
##      [,1] [,2] [,3] [,4]
## [1,]    1    0    0    0
## [2,]    1    1    0    0
## [3,]    1    2    1    0
## [4,]    1    3    3    1
##      [,1] [,2] [,3] [,4]
## [1,]    1    1    1    1
## [2,]    0    1    2    3
## [3,]    0    0    1    3
## [4,]    0    0    0    1
##      [,1] [,2] [,3] [,4]
## [1,] TRUE TRUE TRUE TRUE
## [2,] TRUE TRUE TRUE TRUE
## [3,] TRUE TRUE TRUE TRUE
## [4,] TRUE TRUE TRUE TRUE