Dat 605 HW2 - LU Decomposition Exercise

Problem set 1

library(matlib)
library(MASS)

• The special case where A\(^{T}\)A \(=\) AA\(^{T}\) thru an example (Identity Matrix) in R

m <- matrix(c(1, 0, 0, 0,1,0,0,0,1), nrow = 3, ncol = 3, byrow=F)
m

##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0    1    0
## [3,]    0    0    1

m%*%t(m)

##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0    1    0
## [3,]    0    0    1

t(m)%*%m

##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0    1    0
## [3,]    0    0    1

• But in general, The regular case of A\(^{T}\)A \(!=\) AA\(^{T}\)

A = matrix(c(1,2,3,1,1,1,2,0,1),nrow=3)
A

##      [,1] [,2] [,3]
## [1,]    1    1    2
## [2,]    2    1    0
## [3,]    3    1    1

A%*%t(A)

##      [,1] [,2] [,3]
## [1,]    6    3    6
## [2,]    3    5    7
## [3,]    6    7   11

t(A)%*%A

##      [,1] [,2] [,3]
## [1,]   14    6    5
## [2,]    6    3    3
## [3,]    5    3    5

Problem Set 2

• LU Decomposition Function

LUDecomposition <- function(A, n){

  U = A
  L = diag(n)
  p <- n
 
# Check to see if 1st pivot is 1, if not make it so 
if(A[1,1]> 0) {
U[1,] <- A[1,]/A[1,1]

} else {
if(A[1,1] == 0) {
U[1,] <- A[1,]+1
} else {
U[1,] <- A[1,]/A[1,1]
}
}
 
# populating off centers of Upper & lower matrix   
 i <- 2
 j <- i
 while (j < p+1) {
   
  L[j,i-1 ] <- U[j,i-1]  

  U[j, ] <- U[i-1, ] * U[j, i-1] - U[j, ] 
  j <- j+1
 }
 
 j <-p
 i <-j
 
 # Check to see if 2nd pivot is 1, if not make it so
 
 if(U[i-1,i-1] == 1 ){
   U[i-1,] <- U[i-1,]/U[i-1,i-1]
   
 } else if (U[i-1,i-1] >0){
   U[i-1,] <- U[i-1,]/U[i-1,i-1]
   
 } else if (U[i-1,i-1] <0  ){
   U[i-1,] <- 1*(U[i-1,]/U[i-1,i-1])
   
 } else 
   U[i-1,] <- U[i-1,]/U[i-1,i-1]
   

 while (j < p+1) 
  {
  L[j,i-1] <- U[j,i-1] 

  U[j, ] <- U[i-1, ] *U[j, i-1] - U[j, ] 
  j<-j+1
  }

print("L Matrix: " )
print(L)
print("U Matrix: ")
print(U)  
print("To verify that A = LU:")
Acheck <- L%*%U
return (Acheck )
  
}

n<-3
#A = matrix(c(2,-1,-2,-4,6,3,-4,-2,8),nrow=3,byrow=T)

A = matrix(c(1,2,3,1,1,1,2,0,1),nrow=n)

LUDecomposition(A, n)

## [1] "L Matrix: "
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    2    1    0
## [3,]    3    2    1
## [1] "U Matrix: "
##      [,1] [,2] [,3]
## [1,]    1    1    2
## [2,]    0    1    4
## [3,]    0    0    3
## [1] "To verify that A = LU:"

##      [,1] [,2] [,3]
## [1,]    1    1    2
## [2,]    2    3    8
## [3,]    3    5   17

print("Original Matrix A")

## [1] "Original Matrix A"

print(A)

##      [,1] [,2] [,3]
## [1,]    1    1    2
## [2,]    2    1    0
## [3,]    3    1    1

LU decomposition is Not Unique

Observation:

• I found the decomposed Lower & Upper Matrix as shown above

• However, It seems that LU decomposition is Not Unique unless some kind of constraint/s are placed on it