JBochnik_Assignment2_PS2

function to get L and U matrices from LU decomposition of a square matrix A

getLU_Decomp <- function(A) {
  n <- nrow(A)
  if (n != ncol(A)) {
    stop("Input matrix must be square.")
  }
  
  # Initialize matrices L and U
  L <- matrix(0, n, n)
  diag(L) <- 1;
  U <- matrix(0, n, n)
  
  # Perform LU decomposition
  for (i in 1:n) {
    # Set values starting at the diagonal for U
    U[i, i:n] <- A[i, i:n]

    # Compute Lower triangular matrix (L)
    if (i < n) {
      L[(i+1):n, i] <- A[(i+1):n, i] / U[i, i]
      A[(i+1):n, (i+1):n] <- A[(i+1):n, (i+1):n] - outer(L[(i+1):n, i], U[i, (i+1):n])
    }
  }
  # Set diagonal of L to 1
  diag(L) <- 1
  
  # Return the result as a list
  result <- list(L = L, U = U)
  return(result)
}

A <- matrix(c(2, 1, -1, 8, -3, -1, 2, -11, -2, 1, 2, -3, 2, 1, 1, -1), nrow = 4)
A

##      [,1] [,2] [,3] [,4]
## [1,]    2   -3   -2    2
## [2,]    1   -1    1    1
## [3,]   -1    2    2    1
## [4,]    8  -11   -3   -1

lu<- getLU_Decomp(A)
lu

## $L
##      [,1] [,2] [,3] [,4]
## [1,]  1.0    0    0    0
## [2,]  0.5    1    0    0
## [3,] -0.5    1    1    0
## [4,]  4.0    2   -1    1
## 
## $U
##      [,1] [,2] [,3] [,4]
## [1,]    2 -3.0   -2    2
## [2,]    0  0.5    2    0
## [3,]    0  0.0   -1    2
## [4,]    0  0.0    0   -7

Check the results by checking if A = LU

L <- lu$L
U <- lu$U
LxU <-L %*% U
LxU

##      [,1] [,2] [,3] [,4]
## [1,]    2   -3   -2    2
## [2,]    1   -1    1    1
## [3,]   -1    2    2    1
## [4,]    8  -11   -3   -1

A

##      [,1] [,2] [,3] [,4]
## [1,]    2   -3   -2    2
## [2,]    1   -1    1    1
## [3,]   -1    2    2    1
## [4,]    8  -11   -3   -1

Generate 100 random matrices to test on

# Set the seed for reproducibility
set.seed(26);

# Function to generate a random square matrix of size n
generate_random_matrix <- function(n) {
  matrix(rnorm(n^2), nrow = n);
}

# Number of matrices to generate
num_matrices <- 100;

# Maximum size for matrices
max_size <- 4;

# Generate a list of random square matrices
matrix_list <- vector("list", length = num_matrices)
for (i in 1:num_matrices) {
  matrix_list[[i]] <- generate_random_matrix(sample(2:max_size, 1));
}

Test code on all the random matrices

test_results <- list();
for(i in matrix_list){
  LU <- getLU_Decomp(i);
  L <- LU$L
  U <- LU$U
  # use the all.equal function to account for precision issues (adds a tolerance to to the equality)
  if(all.equal(i, L%*%U)) {
    test_results <- append(test_results, TRUE);
  }else{
    test_results <- append(test_results, FALSE);
  }
}

test_results;

## [[1]]
## [1] TRUE
## 
## [[2]]
## [1] TRUE
## 
## [[3]]
## [1] TRUE
## 
## [[4]]
## [1] TRUE
## 
## [[5]]
## [1] TRUE
## 
## [[6]]
## [1] TRUE
## 
## [[7]]
## [1] TRUE
## 
## [[8]]
## [1] TRUE
## 
## [[9]]
## [1] TRUE
## 
## [[10]]
## [1] TRUE
## 
## [[11]]
## [1] TRUE
## 
## [[12]]
## [1] TRUE
## 
## [[13]]
## [1] TRUE
## 
## [[14]]
## [1] TRUE
## 
## [[15]]
## [1] TRUE
## 
## [[16]]
## [1] TRUE
## 
## [[17]]
## [1] TRUE
## 
## [[18]]
## [1] TRUE
## 
## [[19]]
## [1] TRUE
## 
## [[20]]
## [1] TRUE
## 
## [[21]]
## [1] TRUE
## 
## [[22]]
## [1] TRUE
## 
## [[23]]
## [1] TRUE
## 
## [[24]]
## [1] TRUE
## 
## [[25]]
## [1] TRUE
## 
## [[26]]
## [1] TRUE
## 
## [[27]]
## [1] TRUE
## 
## [[28]]
## [1] TRUE
## 
## [[29]]
## [1] TRUE
## 
## [[30]]
## [1] TRUE
## 
## [[31]]
## [1] TRUE
## 
## [[32]]
## [1] TRUE
## 
## [[33]]
## [1] TRUE
## 
## [[34]]
## [1] TRUE
## 
## [[35]]
## [1] TRUE
## 
## [[36]]
## [1] TRUE
## 
## [[37]]
## [1] TRUE
## 
## [[38]]
## [1] TRUE
## 
## [[39]]
## [1] TRUE
## 
## [[40]]
## [1] TRUE
## 
## [[41]]
## [1] TRUE
## 
## [[42]]
## [1] TRUE
## 
## [[43]]
## [1] TRUE
## 
## [[44]]
## [1] TRUE
## 
## [[45]]
## [1] TRUE
## 
## [[46]]
## [1] TRUE
## 
## [[47]]
## [1] TRUE
## 
## [[48]]
## [1] TRUE
## 
## [[49]]
## [1] TRUE
## 
## [[50]]
## [1] TRUE
## 
## [[51]]
## [1] TRUE
## 
## [[52]]
## [1] TRUE
## 
## [[53]]
## [1] TRUE
## 
## [[54]]
## [1] TRUE
## 
## [[55]]
## [1] TRUE
## 
## [[56]]
## [1] TRUE
## 
## [[57]]
## [1] TRUE
## 
## [[58]]
## [1] TRUE
## 
## [[59]]
## [1] TRUE
## 
## [[60]]
## [1] TRUE
## 
## [[61]]
## [1] TRUE
## 
## [[62]]
## [1] TRUE
## 
## [[63]]
## [1] TRUE
## 
## [[64]]
## [1] TRUE
## 
## [[65]]
## [1] TRUE
## 
## [[66]]
## [1] TRUE
## 
## [[67]]
## [1] TRUE
## 
## [[68]]
## [1] TRUE
## 
## [[69]]
## [1] TRUE
## 
## [[70]]
## [1] TRUE
## 
## [[71]]
## [1] TRUE
## 
## [[72]]
## [1] TRUE
## 
## [[73]]
## [1] TRUE
## 
## [[74]]
## [1] TRUE
## 
## [[75]]
## [1] TRUE
## 
## [[76]]
## [1] TRUE
## 
## [[77]]
## [1] TRUE
## 
## [[78]]
## [1] TRUE
## 
## [[79]]
## [1] TRUE
## 
## [[80]]
## [1] TRUE
## 
## [[81]]
## [1] TRUE
## 
## [[82]]
## [1] TRUE
## 
## [[83]]
## [1] TRUE
## 
## [[84]]
## [1] TRUE
## 
## [[85]]
## [1] TRUE
## 
## [[86]]
## [1] TRUE
## 
## [[87]]
## [1] TRUE
## 
## [[88]]
## [1] TRUE
## 
## [[89]]
## [1] TRUE
## 
## [[90]]
## [1] TRUE
## 
## [[91]]
## [1] TRUE
## 
## [[92]]
## [1] TRUE
## 
## [[93]]
## [1] TRUE
## 
## [[94]]
## [1] TRUE
## 
## [[95]]
## [1] TRUE
## 
## [[96]]
## [1] TRUE
## 
## [[97]]
## [1] TRUE
## 
## [[98]]
## [1] TRUE
## 
## [[99]]
## [1] TRUE
## 
## [[100]]
## [1] TRUE