Question 1

Show that \(A^TA \neq A A^T\) in general. (Proof and demonstration.)

For a 2x2 matrix suppose: \(\mathbf{A} = \left[\begin{array} {rrr} a & b \\ c & d \\ \end{array}\right]\) and \(\mathbf{A^T} = \left[\begin{array} {rrr} a & c \\ b & d \\ \end{array}\right]\) then

\(AA^T = \left[\begin{array}{rrr}a & b \\c & d \\\end{array}\right] \left[\begin{array}{rrr}a & c \\b & d \\\end{array}\right]\) \(= \left[\begin{array}{rrr}aa+bb & ac+bd \\ac+bd & cc+dd \\\end{array}\right]\)

\(A^TA = \left[\begin{array}{rrr}a & c \\b & d \\\end{array}\right] \left[\begin{array}{rrr}a & b \\c & d \\\end{array}\right]\) \(= \left[\begin{array}{rrr}aa+cc & ab+cd \\ab+cd & bb+dd \\\end{array}\right]\)

Thus \(AA^T = \left[\begin{array}{rrr}aa+bb & ac+bd \\ac+bd & cc+dd \\\end{array}\right] = \left[\begin{array}{rrr}aa+cc & ab+cd \\ab+cd & bb+dd \\\end{array}\right]\) if \(c=b\).

That is, if a matrix is equal to its transpose then trivially \(AA^T = A^TA\)

A <- matrix(c(1, 2, -3, 1), 2, 2, byrow=TRUE)
print(A%*%t(A))
##      [,1] [,2]
## [1,]    5   -1
## [2,]   -1   10
print(t(A)%*%A)
##      [,1] [,2]
## [1,]   10   -1
## [2,]   -1    5
 A <- matrix(c(1, 2, 3, 2, 1, 4, 3, 4, 1), 3, 3, byrow=TRUE)
print(A%*%t(A))
##      [,1] [,2] [,3]
## [1,]   14   16   14
## [2,]   16   21   14
## [3,]   14   14   26
print(t(A)%*%A)
##      [,1] [,2] [,3]
## [1,]   14   16   14
## [2,]   16   21   14
## [3,]   14   14   26

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).

As discussed above, a matrix whose transpose is equal to itself trivially makes the equation true. Further, a special matrix called orthogonal has the property \(AA^T = I = A^TA\) and will therefore always make this equation true.

 A <- matrix(c(0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0), 4, 4, byrow=TRUE)
print(A%*%t(A))
##      [,1] [,2] [,3] [,4]
## [1,]    1    0    0    0
## [2,]    0    1    0    0
## [3,]    0    0    1    0
## [4,]    0    0    0    1
print(t(A)%*%A)
##      [,1] [,2] [,3] [,4]
## [1,]    1    0    0    0
## [2,]    0    1    0    0
## [3,]    0    0    1    0
## [4,]    0    0    0    1

Question 2

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

LU.Decomp <- function(A){
  if(nrow(A) != ncol(A)){
    print('Function requires square matrix. Exiting...')
    return(list('L'=NA, 'U'=NA))
  }
  
  U <- A
  L <-diag(nrow(A))
  row <- 1
  column <- 1
  
  while(row < nrow(U)){
    
    #Check for 0 in pivot position. Cannot change rows. No LU decomp
    if(U[row, column] == 0){
      print('0 in pivot position. Exiting...')
      return(list('L'=NA, 'U'=NA))
    }
    
    #Find constants to multiply each row by and create mask to transform U matrix
    constants <- U[(row+1):nrow(U), column] / U[row, column]
    mask <- matrix(rep(constants, each=ncol(U)) * -1 * rep(U[row, ], length(constants)), nrow=nrow(U)-row, byrow=TRUE)
    
    #Modify L and U matrices
    L[(row+1):nrow(U), column] <- constants
    U[(row+1):nrow(A), ] <- mask + U[(row+1):nrow(U), ]
    
    row <- row + 1
    column <- column + 1
  }
  print(L %*% U)
  print(L %*% U == A)
  return(list('L'=L, 'U'=U))
}
A <- matrix(c(1, 2, 3, 1, 1, 1, 2, 0, 1), nrow=3)
print(A)
##      [,1] [,2] [,3]
## [1,]    1    1    2
## [2,]    2    1    0
## [3,]    3    1    1
LU.Decomp(A)
##      [,1] [,2] [,3]
## [1,]    1    1    2
## [2,]    2    1    0
## [3,]    3    1    1
##      [,1] [,2] [,3]
## [1,] TRUE TRUE TRUE
## [2,] TRUE TRUE TRUE
## [3,] TRUE TRUE TRUE
## $L
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    2    1    0
## [3,]    3    2    1
## 
## $U
##      [,1] [,2] [,3]
## [1,]    1    1    2
## [2,]    0   -1   -4
## [3,]    0    0    3
set.seed(123)
A <- matrix(floor(runif(25, 1, 10)), nrow=5)
print(A)
##      [,1] [,2] [,3] [,4] [,5]
## [1,]    3    1    9    9    9
## [2,]    8    5    5    3    7
## [3,]    4    9    7    1    6
## [4,]    8    5    6    3    9
## [5,]    9    5    1    9    6
LU.Decomp(A)
##      [,1] [,2] [,3] [,4] [,5]
## [1,]    3    1    9    9    9
## [2,]    8    5    5    3    7
## [3,]    4    9    7    1    6
## [4,]    8    5    6    3    9
## [5,]    9    5    1    9    6
##      [,1] [,2] [,3] [,4] [,5]
## [1,] TRUE TRUE TRUE TRUE TRUE
## [2,] TRUE TRUE TRUE TRUE TRUE
## [3,] TRUE TRUE TRUE TRUE TRUE
## [4,] TRUE TRUE TRUE TRUE TRUE
## [5,] TRUE TRUE TRUE TRUE TRUE
## $L
##          [,1]      [,2]        [,3]      [,4] [,5]
## [1,] 1.000000 0.0000000  0.00000000  0.000000    0
## [2,] 2.666667 1.0000000  0.00000000  0.000000    0
## [3,] 1.333333 3.2857143  1.00000000  0.000000    0
## [4,] 2.666667 1.0000000  0.01741294  1.000000    0
## [5,] 3.000000 0.8571429 -0.16915423 -9.714286    1
## 
## $U
##      [,1]     [,2]          [,3]      [,4]       [,5]
## [1,]    3 1.000000  9.000000e+00   9.00000   9.000000
## [2,]    0 2.333333 -1.900000e+01 -21.00000 -17.000000
## [3,]    0 0.000000  5.742857e+01  58.00000  49.857143
## [4,]    0 0.000000  1.110223e-16  -1.00995   1.131841
## [5,]    0 0.000000  1.078502e-15   0.00000  13.000000