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?

Answer 1

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

Answer 2

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

Answer 3

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

Answer 4

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