R Markdown

Problem 1

##1. Show that A^T A = AA^T in general. (Proof and demonstration.)

A <- matrix(1:9, 3, 3)
A

##      [,1] [,2] [,3]
## [1,]    1    4    7
## [2,]    2    5    8
## [3,]    3    6    9

AT <- t(A)
AT

##      [,1] [,2] [,3]
## [1,]    1    2    3
## [2,]    4    5    6
## [3,]    7    8    9

x <- A %*% AT
x

##      [,1] [,2] [,3]
## [1,]   66   78   90
## [2,]   78   93  108
## [3,]   90  108  126

y <- AT %*% A
y

##      [,1] [,2] [,3]
## [1,]   14   32   50
## [2,]   32   77  122
## [3,]   50  122  194

x != y

##      [,1] [,2] [,3]
## [1,] TRUE TRUE TRUE
## [2,] TRUE TRUE TRUE
## [3,] TRUE TRUE TRUE

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

A = diag(3)

AT = t(A)

x <- A %*% AT

y <- AT %*% A

x==y

##      [,1] [,2] [,3]
## [1,] TRUE TRUE TRUE
## [2,] TRUE TRUE TRUE
## [3,] TRUE TRUE TRUE

Problem 2

###2. Matrix factorization is a very important problem. There are supercomputers built just to do matrix factorizations. Every second you are on an airplane, matrices are being factorized. Radars that track fights use a technique called Kalman filltering. At the heart of Kalman Filtering is a Matrix Factorization operation. Kalman Filters are solving linear systems of equations when they track your light using radars. Write an R function to factorize a square matrix A into LU or LDU, whichever you prefer.

#Calculate Upper triangular matrix U and Lower matrix L

lu_function <- function(matrix) {

  #Return if the input matrix is not a square one
  if(nrow(matrix) != ncol(matrix)) {
    print('Please provide a square matrix input')
  } else {
    # Initialize a diagonal matrix for the purpose of lower triangular matrix
    L = diag(1, nrow(matrix))
    
    # Create an upper triangular matrix
    U = matrix
    
    #Loop over the matrix
    for(c in 1:(ncol(U)-1)) {
      v1 <- U[c,c]
      for (r in 1:(nrow(U)-c)) {
        v2 <- U[c+r, c]
        L[c+r,c] = v2/v1
        U[c+r,] = (-v2/v1)*U[c,] + U[c+r,]
      }
    }
    print("Lower Matrix(L) = ")
    print(L)
    
    print("Upper Matrix(U) = ")
    print(U)
    
    A2 = L%*%U
    
    #Test whether the product of L and U is equal to original matrix
    matrix==A2
  }
}

lu_function(matrix(c(2,1,-6,4,-4,-9,-4,3,5), 3, 3))

## [1] "Lower Matrix(L) = "
##      [,1] [,2] [,3]
## [1,]  1.0  0.0    0
## [2,]  0.5  1.0    0
## [3,] -3.0 -0.5    1
## [1] "Upper Matrix(U) = "
##      [,1] [,2] [,3]
## [1,]    2    4 -4.0
## [2,]    0   -6  5.0
## [3,]    0    0 -4.5

##      [,1] [,2] [,3]
## [1,] TRUE TRUE TRUE
## [2,] TRUE TRUE TRUE
## [3,] TRUE TRUE TRUE