0

I need to solve more than thousand matrices in a list. However, I get the error Lapack routine dgesv: system is exactly singular. My problem is that my input data are non singular matrices, but following the calculations I need to do on the matrices some of them get singular. A reproducible example with a subset of my data is however not possible as it would be to long (I tried already). Here a basic example of my problem (A would be a matrix after some calculations and R the next calculation I need to do):

A=matrix(c(1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1), nrow=4)
R = solve(diag(4)-A)

Do you have ideas how to "solve" this problem, maybe other function? Or how to transform the singular matrices very very slightly in order to not create totally different results? Thanks

EDIT: according to @Roman Susi I include the function (with all calculations) I have to do:

function(Z, p) {
  imp <- as.vector(cbind(imp=rowSums(Z)))
  exp <- as.vector(t(cbind(exp=colSums(Z))))
  x = p + imp
  ac = p + imp - exp  
  einsdurchx = 1/as.vector(x)    
  einsdurchx[is.infinite(einsdurchx)] <- 0
  A = Z %*% diag(einsdurchx)
  R = solve(diag(length(p))-A) %*% diag(p)    
  C = ac * einsdurchx
  R_bar = diag(as.vector(C)) %*% R
  rR_bar = round(R_bar)
  return(rR_bar)
}

The problem is in line 8 of the function calculating solve(diag(length(p))-A). Here I can provide new example data for Z and p, however in this example it works fine as I was not able to recreate an example which brings the error:

p = c(200, 1000, 100, 10)
Z = matrix(
  c(0,0,100,200,0,0,0,0,50,350,0,50,50,200,200,0),
  nrow = 4, 
  ncol = 4,
  byrow = T)

So, the question according to @Roman Susi is: Is there a way to change the calculations before so that det(diag(length(p))-A) never gets 0 in order to solve the equation? I hope you can understand what I want:) Ideas, thanks. Edit2: Maybe asked easier: How to avoid singularity within this function (at least before line 8)?

Josh O'Brien
  • 159,210
  • 26
  • 366
  • 455
N.Varela
  • 910
  • 1
  • 11
  • 25
  • `?try` or `?tryCatch` – rawr Nov 06 '15 at 19:22
  • From the question it is not clear if this is a math problem rather than programming one. If a singular matrix to invert is an intermediate result, then maybe a different formula will allow to avoid "division by zero"? Because if singularity can't be avoided, it is inherent in the problem, and the problem should be reformulated broader or accuracy increased, etc, it depends. – Roman Susi Nov 06 '15 at 19:27
  • When you pass a matrix to `solve`, your telling it to find the inverse matrix, but singular matrices don't have an inverse. I think you're asking it to do the impossible.... – Jthorpe Nov 06 '15 at 19:36

2 Answers2

3

The generalized inverse ginv in the MASS package can handle singular matrices but whether it makes sense or not for your problem would have to be determined.

 library(MASS)
 ginv(diag(4) - A)

giving:

     [,1] [,2] [,3] [,4]
[1,]    0    0    0    0
[2,]    0    0    0    0
[3,]    0    0    0    0
[4,]    0    0    0    0

There is also the Ginv function in the ibdreg package.

G. Grothendieck
  • 254,981
  • 17
  • 203
  • 341
  • This isn't an inverse in the general sense -- just a programming convenience / inconvenience (I'd always rather have an error than an answer that may or may not make sense!) – Jthorpe Nov 06 '15 at 19:39
  • ...or perhaps "whether it makes sense or not would have to be overdetermined". Sorry, couldn't resist. – Brian B Nov 06 '15 at 19:39
  • @Jthorpe, It's an inverse on the subset of the domain consisting of the orthogonal complement to the nullspace. The example given in the question is not particularly illustrative since the nullspace is the entire space in that case but that will not be the situation in general. – G. Grothendieck Nov 06 '15 at 22:41
0

R's QR decomposition functions may have your answer. They provide a way to solve linear equations robustly. QR decomposition does not provide an inverse, but rather a matrix decomposition that can often be used where an inverse would be used.

For rectangular matrices QR decomposition can be used to find the least-squares fit. For square, (near) singular matrices, qr() detects this (near) singularity, and qr.coef() can then be used to obtain a solution without any errors but possibly some NAs (which can be converted to zeros).

Paul Harrison
  • 455
  • 3
  • 6