How to solve singular matrices?

Question

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

Answer 1

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

Answer 2

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.