带状矩阵的 LU 分解的 R 代码

Question

使用 R 实现计算带状矩阵的 LU 分解的算法（可以忽略旋转）。函数的输入应该是矩阵 A 及其带宽 w（您的算法不需要确定带宽，您可以假设它将由用户提供）。输出应该是两个矩阵 L 和 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))

输出错误：A[(i + 1):end_col, (i + 1):end_col] - L[(i + 1):end_col, i] %*% 中的错误： 非一致性数组

Answer 1

如果将 drop=FALSE 参数添加到 L 和 U 的索引中，您将获得更进一步的结果。 R 的

[

函数将删除单列或行矩阵的维度。

  ... - 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

但我不知道这是否正确。 （它确实具有带状结构，所以可能是。）