Sil - 1 year ago 88
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!

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

Recommended from our users: Dynamic Network Monitoring from WhatsUp Gold from IPSwitch. Free Download