用于在 R 中实现块式反演技术的递归函数

Question

我对线性代数很陌生，我正在尝试实现一个递归函数，它使用在不使用 R 库的情况下从头开始的逐块反演技术 "solve".

此问题已在以下 post 中回答：function for matrix

然而，它对我不起作用，我尝试实现我自己的版本：

matrixInversion <- function(M){
   if(nrow(M) == 2){
      a <- M[1,1]
      b <- M[1,2]
      c <- M[2,1]
      d <- M[2,2]
      deter <- ((a*d)-(b*c))
      InverseMatrix <- ((1/deter)*matrix(c(d,-c,-b,a),nrow=2,ncol=2))
   } else {
      x <- (floor(nrow(M) / 2))
      A <- M[1:x, 1:x, drop=F]
      B <- M[1:x, -1:-x, drop=F]
      C <- M[-1:-x, 1:x, drop=F]
      D <- M[-1:-x, -1:-x, drop=F]

      Ainv <- matrixInversion(A)
      common <- matrixInversion(D - C %*% Ainv %*% B)
      newA <- Ainv+Ainv%*%B%*%common%*%C%*%Ainv
      newB <- (-Ainv)%*%B%*%common
      newC <- (-common)%*%C%*%Ainv
      newD <- (-common)

      result <- cbind(rbind(newA, newC), rbind(newB, newD))
    }
}

此版本适用于偶数列的矩阵，但不适用于奇数列的矩阵。我不明白如何正确实施它。有什么建议吗？

谢谢！

Answer 1

除了评论中提出的观点外，您还错误地否定了newD。您的函数经过更改后现在可以同时处理奇数行和偶数行矩阵：

matrixInversion <- function(M){
  if (nrow(M) == 1){
    return(1/M)
  } else if (nrow(M) == 2){
    a <- M[1,1]
    b <- M[1,2]
    c <- M[2,1]
    d <- M[2,2]
    deter <- ((a*d)-(b*c))
    return((1/deter)*matrix(c(d,-c,-b,a),nrow=2,ncol=2))
  } else {
    x <- (floor(nrow(M) / 2))
    A <- M[1:x, 1:x, drop=F]
    B <- M[1:x, -1:-x, drop=F]
    C <- M[-1:-x, 1:x, drop=F]
    D <- M[-1:-x, -1:-x, drop=F]
    Ainv <- matrixInversion(A)
    common <- matrixInversion(D - C %*% Ainv %*% B)
    newA <- Ainv+Ainv%*%B%*%common%*%C%*%Ainv
    newB <- (-Ainv)%*%B%*%common
    newC <- (-common)%*%C%*%Ainv
    newD <- common
    return(cbind(rbind(newA, newC), rbind(newB, newD)))
  }
}

## random seed not relevant here ...
m1 <- matrix(runif(150^2), nr=150)
all.equal(solve(m1), matrixInversion(m1))
# [1] TRUE
any(is.nan(matrixInversion(m1))) # based on your comment
# [1] FALSE

m2 <- matrix(runif(151^2), nr=151)
all.equal(solve(m2), matrixInversion(m2))
# [1] TRUE
any(is.nan(matrixInversion(m2)))
# [1] FALSE

郑重声明，这是我通过逐步检查更大的矩阵发现的。类似于：

m1 <- matrix(runif(150^2), nr=150)
i <- 4; all.equal(solve(m1[1:i,1:i]), matrixInversion(m1[1:i,1:i]))

在 2 时，它是相等的，并且开始突破该值。我使用 debugonce(matrixInversion) 检查前两个条件块的内部，但它们没问题。然后我将你的方程与最初引用的 wiki-page 进行了比较，并注意到你的计算不正确。当然，这是我检查的最后一个:-)

用于在 R 中实现块式反演技术的递归函数

Recursive function for implementing block-wise inversion technique in R

r

recursion

matrix

matrix-inverse