R Code for LU Factorization for a Banded Matrix

Question

Use R to implement an algorithm to compute the LU factorization of banded matrices (you can ignore pivoting). The inputs to your function should be a matrix A and its bandwidth w (your algorithm does not need to determine the bandwidth, you can assume it will be provided by the user). The outputs should be a list of two matrices L and U .

lu_with_bandwidth <- function(A, bandwidth) {
  n <- nrow(A)
  L <- diag(n)
  U <- matrix(0, nrow = n, ncol = n)
  
  for (i in 1:n) {
    start_col <- max(1, i - bandwidth)
    end_col <- min(n, i + bandwidth)
    
    U[i, i:end_col] <- A[i, i:end_col]
    
    if (start_col <= end_col) {
      L[start_col:end_col, i] <- A[start_col:end_col, i] / U[i, i]
      if (i < n) {
        A[(i+1):end_col, (i+1):end_col] <- A[(i+1):end_col, (i+1):end_col] - L[(i+1):end_col, i] %*% U[i, (i+1):end_col]
      }
    }
  }
  
  return(list(L = L, U = U))
}
B <- matrix (c(3,1,4,0,0,0,1,5,9,2,0,0,6,5,3,5,8,0,0,9,7,9,3,2,0,0,3,8,4,6,0,0,0,2,6,4), nrow=6)
print(lu_with_bandwidth(B,2))

Output error: Error in A[(i + 1):end_col, (i + 1):end_col] - L[(i + 1):end_col, i] %*% : non-conformable arrays

here is the algorithm i follow. algorithm

Answer 1

If you add the drop=FALSE parameter to the indexing of L and U you will get further. R's [ function will drop dimensions for single column or row matrices.

  ... - L[(i+1):end_col, i, drop=FALSE] %*% U[i, (i+1):end_col, drop=FALSE]

#no prettying screen scraping:
> lu_with_bandwidth <- function(A, bandwidth) {
+     n <- nrow(A)
+     L <- diag(n)
+     U <- matrix(0, nrow = n, ncol = n)
+     
+     for (i in 1:n) {
+         start_col <- max(1, i - bandwidth)
+         end_col <- min(n, i + bandwidth)
+         
+         U[i, i:end_col] <- A[i, i:end_col]
+         
+         if (start_col <= end_col) {
+             L[start_col:end_col, i] <- A[start_col:end_col, i] / U[i, i]
+             if (i < n) {cat("starting")
+                 A[(i+1):end_col, (i+1):end_col] <- A[(i+1):end_col, (i+1):end_col] - L[(i+1):end_col, i, drop=FALSE] %*% U[i, (i+1):end_col, drop=FALSE]
+             }
+         }
+     }
+     
+     return(list(L = L, U = U))
+ }
> B <- matrix (c(3,1,4,0,0,0,1,5,9,2,0,0,6,5,3,5,8,0,0,9,7,9,3,2,0,0,3,8,4,6,0,0,0,2,6,4), nrow=6)
> print(lu_with_bandwidth(B,2))
startingstartingstartingstartingstarting$L
          [,1]      [,2]       [,3]       [,4]       [,5]      [,6]
[1,] 1.0000000 0.2142857 -0.6043165  0.0000000  0.0000000 0.0000000
[2,] 0.3333333 1.0000000 -0.3021583  4.0354839  0.0000000 0.0000000
[3,] 1.3333333 1.6428571  1.0000000 -3.4910138  0.1514658 0.0000000
[4,] 0.0000000 0.4285714 -0.3741007  1.0000000  0.4605723 0.6269144
[5,] 0.0000000 0.0000000 -0.8057554 -1.4677419  1.0000000 2.8008919
[6,] 0.0000000 0.0000000  0.0000000  0.8967742 -0.1100977 1.0000000

$U
     [,1]     [,2]      [,3]      [,4]      [,5]     [,6]
[1,]    3 1.000000  6.000000  0.000000  0.000000 0.000000
[2,]    0 4.666667  3.000000  9.000000  0.000000 0.000000
[3,]    0 0.000000 -9.928571 -7.785714  3.000000 0.000000
[4,]    0 0.000000  0.000000  2.230216  9.122302 2.000000
[5,]    0 0.000000  0.000000  0.000000 19.806452 8.935484
[6,]    0 0.000000  0.000000  0.000000  0.000000 3.190228

But I don't know if that is correct or not. (It does have a banded structure so it might be.)