title: How to solve singular matrices?
body: <p>I need to <code>solve</code> more than thousand matrices in a list. However, I get the error <code>Lapack routine dgesv: system is exactly singular</code>. 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):</p>
  
  <pre><code>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)
  </code></pre>
  
  <p>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 </p>
  
  <p><strong>EDIT:</strong> according to @Roman Susi I include the function (with all calculations) I have to do:</p>
  
  <pre><code>function(Z, p) {
    imp &lt;- as.vector(cbind(imp=rowSums(Z)))
    exp &lt;- as.vector(t(cbind(exp=colSums(Z))))
    x = p + imp
    ac = p + imp - exp  
    einsdurchx = 1/as.vector(x)    
    einsdurchx[is.infinite(einsdurchx)] &lt;- 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)
  }
  </code></pre>
  
  <p>The problem is in line 8 of the function calculating <code>solve(diag(length(p))-A)</code>. Here I can provide new example data for <code>Z</code> and <code>p</code>, however in this example it works fine as I was not able to recreate an example which brings the error:</p>
  
  <pre><code>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) 
  </code></pre>
  
  <p>So, the question according to @Roman Susi is: Is there a way to change the calculations before so that <code>det(diag(length(p))-A)</code> never gets 0 in order to <code>solve</code> the equation? I hope you can understand what I want:) Ideas, thanks. <strong>Edit2:</strong> Maybe asked easier: How to avoid singularity within this function (at least before line 8)? </p>