Problem 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 flights use a technique called Kalman filtering. At the heart of Kalman Filtering is a Matrix Factorization operation.Kalman Filters are solving linear systems of equations when they track your flight using radars.

Write an R function to factorize a square matrix A into LU or LDU, whichever you prefer. Please submit your response in an R Markdown document using our class naming convention, E.g. LFulton_Assignment2_PS2.png

You don’t have to worry about permuting rows of A and you can assume that A is less than 5x5, if you need to hard-code any variables in your code. If you are doing the entire assignment in R, then please submit only one markdown document for both the problems.

Solution:

We can decompose a given square matrix A into an L (a lower triangular matrix) and U (an upper triangular matrix) such that A = LU. In an upper triangular matrix, all elements above the main diagonal are 1 and in a lower diagonal matrix, all elements below the main diagonal are 1.

LU = function(A){
  
#Test if the input matrix A has equal number of rows and columns, i.e. the matrix is square. Otherwise, return an error message.
  
  if(dim(A)[1]!=dim(A)[2]){
    return('Not a square matrix')
  }
  
  # Initialize variables
  
  n <- nrow(A)   #nrow returns the number of rows in matrix A
  U <- A         #Assign original matrix to the uppper triangular matrix
  D <- diag(n)   #Construct a diagonal matrix of dimensions entered (mxm of input square matrix)
  L <- diag(n)

  # Obtaining the upper and the lower triangular matrix
  
  for (j in c(1:n)){
    for(i in c(2:n)){
      if(i > j){
        row <- U[j,]; print(row) #Extract all columns but only for the jth row
        multiplier <- U[i, j] / row[j] #Calculate the multiplier
        U[i,] <- U[i,] - (multiplier * row);print(U) #Use the multiplier to swap out entries in U matrix
        L[i,j] <- multiplier; print(L[i,j]); print(multiplier); print(L) #Use the multiplier calculated to swap entries in L matrix
      }
    }
  }
  
  return(list(L=L, U=U))
}

#Below matrix was used in the LU decomposition video of Week 2 study materials

A = matrix(c(1,-2,3,4,8,4,-3,5,7),nrow=3,ncol = 3)
LU(A)
## [1]  1  4 -3
##      [,1] [,2] [,3]
## [1,]    1    4   -3
## [2,]    0   16   -1
## [3,]    3    4    7
## [1] -2
## [1] -2
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]   -2    1    0
## [3,]    0    0    1
## [1]  1  4 -3
##      [,1] [,2] [,3]
## [1,]    1    4   -3
## [2,]    0   16   -1
## [3,]    0   -8   16
## [1] 3
## [1] 3
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]   -2    1    0
## [3,]    3    0    1
## [1]  0 16 -1
##      [,1] [,2] [,3]
## [1,]    1    4 -3.0
## [2,]    0   16 -1.0
## [3,]    0    0 15.5
## [1] -0.5
## [1] -0.5
##      [,1] [,2] [,3]
## [1,]    1  0.0    0
## [2,]   -2  1.0    0
## [3,]    3 -0.5    1
## $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

Double-check results using built-in R package.

We can double-check using the lu-decomposition package of R.

library(matrixcalc)
luA <- lu.decomposition(A)
L <- luA$L
U <- luA$U
print( L )
##      [,1] [,2] [,3]
## [1,]    1  0.0    0
## [2,]   -2  1.0    0
## [3,]    3 -0.5    1
print( U )
##      [,1] [,2] [,3]
## [1,]    1    4 -3.0
## [2,]    0   16 -1.0
## [3,]    0    0 15.5
print( L %*% U )
##      [,1] [,2] [,3]
## [1,]    1    4   -3
## [2,]   -2    8    5
## [3,]    3    4    7
print(A)
##      [,1] [,2] [,3]
## [1,]    1    4   -3
## [2,]   -2    8    5
## [3,]    3    4    7
#Results are same as obtained with function above.