Sil - 1 year ago 66

R Question

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 Source

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