Assignment 2 DATA605

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 upper 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,] #Extract all columns but only for the jth row
        multiplier <- U[i, j] / row[j] #Calculate the multiplier
        U[i,] <- U[i,] - (multiplier * row) #Use the multiplier to swap out entries in the U matrix
        L[i,j] <- multiplier #Use the multiplier calculated to swap entries in the L matrix
      }
    }
  }
  
  return(list(L=L, U=U))
}

#Below matrix was used as an example in the LU decomposition video of Week 2 study materials and is used here to illustrate results.

A = matrix(c(1,-2,3,4,8,4,-3,5,7),nrow=3,ncol = 3)
LU(A)

## $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 obtained using the man-made function with those obtained using R’s built-in package.

We can double-check using the lu-decomposition package of R and see that results are the same in both cases.

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