#Randomly generated 5 by 5 matrix
ma <- matrix(sample.int(500,size=25,replace=FALSE),nrow=5,ncol=5)
#LU decomposition function
lu_decomp_square <- function(x){
#make sure that x is a square matrix, and larger than 1
if(nrow(x) != ncol(x) | (nrow(x) == 1 & ncol(x) == 1)){
message("This function is only for square matrices greater than 1 by 1")
}
else{
#lower diagonal
l<-diag(nrow(x))
#cycle through each row of the matrix
for(i in 1:nrow(x)){
#cycle through each column of the matrix
for(j in 1:ncol(x)){
#only want the indexes where column index is less than row index
if (j < i){
#formula for row reduction
zero_maker <-x[i,j]/x[j,j]
#add the multiple to the proper index of the lower triangle matrix
l[i,j]<-zero_maker
#row reduce
x[i,]<- x[i,]- zero_maker*x[j,]
}
}
}
}
#initial matrix
print(ma)
#lower matrix
print('L')
print(l)
#upper matrix
print('U')
print(x)
reconstruction<-l%*%x
#The reconstructed matrix, should be same as initial matrix
print("L*U")
reconstruction
}
lu_decomp_square(ma)
## [,1] [,2] [,3] [,4] [,5]
## [1,] 312 420 44 302 45
## [2,] 353 55 213 109 343
## [3,] 113 95 425 325 476
## [4,] 93 145 352 451 342
## [5,] 180 147 456 173 210
## [1] "L"
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1.0000000 0.00000000 0.0000000 0.000000 0
## [2,] 1.1314103 1.00000000 0.0000000 0.000000 0
## [3,] 0.3621795 0.13592677 1.0000000 0.000000 0
## [4,] 0.2980769 -0.04713959 0.8958335 1.000000 0
## [5,] 0.5769231 0.22681922 1.0173595 -1.556176 1
## [1] "U"
## [,1] [,2] [,3] [,4] [,5]
## [1,] 3.120000e+02 420.0000 44.0000 302.0000 45.00000
## [2,] 0.000000e+00 -420.1923 163.2179 -232.6859 292.08654
## [3,] 0.000000e+00 0.0000 386.8784 247.2500 419.99954
## [4,] 0.000000e+00 0.0000 0.0000 128.5172 -33.89428
## [5,] 2.842171e-14 0.0000 0.0000 0.0000 -362.24834
## [1] "L*U"
## [,1] [,2] [,3] [,4] [,5]
## [1,] 312 420 44 302 45
## [2,] 353 55 213 109 343
## [3,] 113 95 425 325 476
## [4,] 93 145 352 451 342
## [5,] 180 147 456 173 210