How to solve singular matrices?

r matrix singular
N.Varela picture N.Varela · Nov 6, 2015 · Viewed 9.2k times · Source

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

G. Grothendieck picture G. Grothendieck · Nov 6, 2015

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.

Related questions

R memory management / cannot allocate vector of size n Mb

I am running into issues trying to use large objects in R. For example: > memory.limit(4000) > a = matrix(NA, 1500000, 60) > a = matrix(NA, 2500000, 60) > a = matrix(NA, 3500000, 60) Error: cannot allocate vector of size 801.1 Mb > a = matrix(NA, 2500000, 60) …

Read More
Subscript out of bounds - general definition and solution?

When working with R I frequently get the error message "subscript out of bounds". For example: # Load necessary libraries and data library(igraph) library(NetData) data(kracknets, package = "NetData") # Reduce dataset to nonzero edges krack_full_nonzero_edges <- …

Read More
Select rows of a matrix that meet a condition

In R with a matrix: one two three four [1,] 1 6 11 16 [2,] 2 7 12 17 [3,] 3 8 11 18 [4,] 4 9 11 19 [5,] 5 10 15 20 I want to extract the submatrix whose rows have column three = 11. That is: one two three four [1,] 1 6 11 16 [3,] 3 8 11 18 [4,] 4 9 11 19 I want to do this without looping. I am new to R so this …

Read More