Sil Sil - 7 days ago 10
R Question

Recursive function for implementing block-wise inversion technique in R

I am very new to linear algebra and I am trying to implement a recursive function which inverts any matrix using
the block-wise inversion technique from scratch without using the R library "solve".

This question was already answer in the following post: function for matrix

However, it did not work for me and I tried to implement my own version:

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))
}
}


This version works successfully for matrix with even number of columns, but not with matrix with odd number of columns. I cannot understand how to implement it correctly. Any suggestion?

Thanks!

Answer

In addition to points made in the comments, you incorrectly negate newD. Your function with changes now works with both odd-row and even-row matrices:

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

For the record, I found this by incrementally checking larger matrices. Something like:

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

At 2, it was equal, and it started breaking above that. I used debugonce(matrixInversion) to check inside the first two conditional blocks, but they were fine. I then compared your equations with the originally-cited wiki-page, and noticed your incorrect calculation. Of course, it was the last one I checked :-)

Comments