This question was moved from Cross-Validated: https://stats.stackexchange.com/questions/274042/pseudoinverse-of-large-sparse-matrix-in-r

I am trying to calculate the pseudoinverse of a large sparse matrix in R using the singular value decomposition. The matrix is roughly 240,000 x 240,000, and I have it stored as type dgCMatrix. I have tried using pinv, ginv, and other standard pseudoinverse functions, but they error out due to memory constraints. I then tried to opt for the sparse matrix svd provided by the package irlba, which I was then going to use to compute the pseudoinverse using the standard formula after converting all outputs to sparse matrices. My code is here:

lim = 40
digits = 4
SVD =irlba(L,lim)
tU = round(SVD$u,digits)
nonZeroU = which(abs(U)>0,arr.ind = T)
sparseU = sparseMatrix(i=nonZeroU[,2],j=nonZeroU[,1],x = U[nonZeroU])
V = round(SVD$v,digits)
nonZeroV = which(abs(V)>0,arr.ind = T)
sparseV = sparseMatrix(i=nonZeroV[,1],j=nonZeroV[,2],x = U[nonZeroV])
D = as(Diagonal(x=1/SVD$d),"sparseMatrx")
pL =D%*%sparseU
pL = sparseV%*%pL

I am able to get to the last line without an issue, but then I get an error due to memory constraints that says

Error in sparseV %*% pL : 
  Cholmod error 'problem too large' at file ../Core/cholmod_dense.c, line 105

Of course I could piece together the pseudoinverse entry by entry using a for loop and vector multiplications, but I would like to be able to calculate it using a simple function that takes advantage of the sparsity of the resultant pseudoinverse matrix. Due to the nature of the original matrix (it is a graph laplacian), I know that the pseudoinverse should also be a sparse matrix. Any help would be greatly appreciated!

  • $\begingroup$ How sparse is your matrix, how sparse do you expect the pseudoinverse to be? A dense 240000^2 matrix of 4 byte numbers requires more than 200GB for storage, let alone for handling. Also, in many cases the inverse of a sparse matrix ends up being a dense matrix. $\endgroup$ – DrHansGruber Apr 17 '17 at 22:26
  • $\begingroup$ The matrix is very sparse-- <1%. The matrix L is a graph laplacian, so I expect the pseudoinverse to be equally sparse (I am saying this from empirical observation but haven't proven it yet) $\endgroup$ – Paul Apr 17 '17 at 23:33
  • 6
    $\begingroup$ Why do you expect the inverse of a Laplacian to be sparse? That's normally not true. $\endgroup$ – cfh Apr 18 '17 at 6:42
  • 1
    $\begingroup$ Do you really need the inverse/pseudoinverse matrix? For your use, you have considered the possibility of solving the linear system for applying the inverse transformation? $\endgroup$ – Mauro Vanzetto Apr 18 '17 at 7:44
  • 1
    $\begingroup$ If the earlier steps work well, can't you keep the matrix in factored form, evaluating matrix-vector products one at a time, like $V(D(Ux))$? Also, $L^{+}_{ij}$ has an interpretation in terms of the resistivity distance between $i$ and $j$, so I expect that for any connected component of the graph $L^{+}$ would be dense. Do you have a specific reason for thinking it would be sparse here? $\endgroup$ – Kirill Apr 18 '17 at 20:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.