implement matrix determinant in R

Question

I was asked to implement function that calculates n-dimensional matrix determinant using Laplace expansion. This involves recursion. I developed this:

minor<-function(A,i,j) {
return(A[c(1:(i-1),(i+1):dim(A)[1]),c(1:(j-1),(j+1):dim(A)[2])])
}


determinantRec<-function(X,k) {
if (dim(X)[1] == 1 && dim(X)[2] == 1) return(X[1][1])
else {
s = 0
for (i in 1:dim(X)[2]) {
  s = s + X[k][i]*(-1)^(k+i)*determinantRec(minor(X,k,i),k)
}
return(s)
}
}

where k in determinantRec(X,k) function indicates which row I want to use Laplace expansion along of.

My problem is when I run determinantRec(matrix(c(1,2,3,4),nrow = 2,ncol = 2),1) this error appears:

C stack usage 7970628 is too close to the limit

What is wrong with my code?

Have you debugged code before in R? Take a look at `?browser` for some options. Your code seems to have boundary condition problems. When `k=1`, then `minor()` is called with `i=1` and `1:(i-1)` is `1:0` which is probably not what you want. — MrFlick, Oct 13 '17 at 18:34
Your `minor` function does not correctly reduce the matrix size when `i` or `j` are 1 (due to @MrFlick's comment above). As a result, `determinantRec` becomes stuck in infinite recursion. — Artem Sokolov, Oct 13 '17 at 18:36
The following question may be helpful: https://stackoverflow.com/questions/12919984/how-to-delete-specific-rows-and-columns-from-a-matrix-in-a-smarter-way — Artem Sokolov, Oct 13 '17 at 18:39

score 0 · Answer 1 · answered Sep 13 '20 at 17:30

@julia, there is one simple type in your code. Just remove the '*' at the end of the definition of 's'. And don't indent the recursion.

determinantRek<-function(X,k) {
if (dim(X)[1] == 1 && dim(X)[2] == 1)
return(X[1,1])
if (dim(X)[1] == 2 && dim(X)[2] == 2)
return(X[1,1]*X[2,2]-X[1,2]*X[2,1])
else
s = 0
for (i in 1:dim(X)[2]) {
  s = s + X[k,i]*(-1)^(k+i)
  determinantRek(X[-k,-i],k)
}
return(s)
}

Gabriel Matte · Answer 2 · 2021-09-21T01:06:16.907

I did this way and works just fine, although it is super slow, compared to the det function in base R

laplace_expansion <- function(mat){
  det1 <- function(mat){
    mat[1]*mat[4]-mat[2]*mat[3]
  }
  determinant <- 0
  for(j in 1:ncol(mat)){
    mat1 <- mat[-1,-j]
    if(nrow(mat1) == 2){
      determinant <- determinant+mat[1,j]*(-1)^(1+j)*det1(mat1)
    }else{
      val <- mat[1,j]*(-1)^(1+j)
      if(val != 0){
        determinant <- determinant+val*laplace_expansion(mat1)
      }
    }
  }
  return(determinant)
}

score 0 · Answer 3 · answered Mar 14 '22 at 22:31

This is my approach, I think it's cleaner.

deter <- function(X) {
  
  stopifnot(is.matrix(X))
  stopifnot(identical(ncol(X), nrow(X))) 
  
  if (all(dim(X) == c(1, 1))) return(as.numeric(X))
  
  i <- 1:nrow(X)
  
  out <- purrr::map_dbl(i, function(i){
    X[i, 1] * (-1)^(i + 1) * deter(X[-i, -1, drop = FALSE])
  })
  
  return(sum(out))
  
}

score -1 · Answer 4 · answered Oct 15 '17 at 18:44

Thank you @ArtemSokolov and @MrFlick for pointing the problem cause, it was it. I also discovered that this code does not calculate properly the determinant of 2x2 matrix. After all it looks like that:

determinantRek<-function(X,k) {
if (dim(X)[1] == 1 && dim(X)[2] == 1)
return(X[1,1])
if (dim(X)[1] == 2 && dim(X)[2] == 2)
return(X[1,1]*X[2,2]-X[1,2]*X[2,1])
else
s = 0
for (i in 1:dim(X)[2]) {
  s = s + X[k,i]*(-1)^(k+i)*
    determinantRek(X[-k,-i],k)
}
return(s)
}

Debuging with browser() was also helpful :)

implement matrix determinant in R

4 Answers4