# Moore-Penrose pseudoinverse of singular rank degenerate matrix

I am trying to attain the Moore-Penrose pseudoinverse of a very large, very sparse, rank-degenerate, singular, and square matrix. ($75000 \times 75000$, near rank). The matrix is a graph Laplacian and I need to find the resistance distance between a large number of nodes (LU determinates are too slow). I realize the inverse will be very dense, but I would like to persist it so I can query out distance metrics ad-hoc . I have been trying some brute approaches by converting it to a dense matrix (~25GB) on a large memory machine (124GB).

Both lstsq (scipy.linalg.pinv) and some SVD approaches (numpy.linalg.pinv) quickly exhaust the memory. The only approach that has worked is a slow but memory efficient SVD (scipy.linalg.pinv2).

I am exploring what other techniques are out there that can balance memory efficiency and speed, although I have yet to try scipy.sparse.inv or scipy.sparse.svds for a thorough comparison. I don't think any iterative approaches will work well because I can't provide a decent starting point. I think this leaves QR as a prime candidate. I have been exploring some of the methods listed in this paper. I have the following questions:

1. I want to sample selected entries of the MPI to calculate the distance between a large number of graph nodes in a ad-hoc fashion. Is there a better way to sample select elements of the MPI than persisting the whole matrix?
2. Are there any algorithms particularly well suited to my matrix type?
3. Is QR factorization my best bet? Is there some other approaches I should be considering?
4. If QR, what memory efficient QR implementations exist? Or will I need to hand roll one from a paper?
• Do you actually need a pseudoinverse or you want to do something with it? Sep 12 '18 at 21:02
• @AntonMenshov I actually need the inverse. I am interested in some of its characteristics. I plan on persisting it to disk, so I can do ad-hoc analysis on it without a recalculation.
– feik
Sep 12 '18 at 21:04
• "I am interested in some of its characteristics" is a vague non-answer.. which ones? I think it's good that you've experimented with the tools you have at hand, but computing the explicit inverse of a sparse matrix is a pretty misguided thing to do, and resorting to the pseudo-inverse is even worse. Sep 12 '18 at 21:32
• To be more specific, I am trying to compute the distance metrics between a large number of nodes using the graph Laplacian.
– feik
Sep 12 '18 at 21:55

Although I'm unfamiliar with "distance metrics", there's no shortage of information about graph laplacians on the web because of their practical applicability to matrix reordering and classification/clustering problems. The nullspace of the undirected graph laplacian is well known (a vector of 1's, or one such vector per connected component if the graph is not strongly connected).

Because the structure of the laplacian is pretty well known, you might not need to resort to brute force techniques like computing a dense pseudoinverse. There are cheaper/faster tools [sparse cholesky, lanczos, lobpcg] that can interrogate these matrices in other ways [solve by it, extract the Fiedler vector or other important eigenvectors, sample selected entries of the inverse]. Maybe they (or something similar) could help?

Could you write a bit more about "distance metrics", and how the pseudoinverse helps you find them?

Instead of seeking entries the pseudoinverse $\mathbf \Gamma$ of your (singular) Laplacian $\mathbf L$, consider the exact inverse $\mathbf \Gamma'$ of a nearby (nonsingular) matrix, $\mathbf L' = \mathbf L + \sigma \mathbf I$, where $\sigma$ is some small shift parameter (say, $1\times10^{-8}$).

I would expect (hope?) that $\mathbf \Gamma'_{i,j}$ is very close to $\mathbf \Gamma_{i,j}$, but the former is vastly easier to compute, because you can just use sparse cholesky decomposition to represent $\mathbf \Gamma' = \left( \mathbf L + \sigma \mathbf I\right)^{-1}$, and then $\mathbf \Gamma'_{ij} = \mathbf e_i \left( \mathbf L + \sigma \mathbf I\right)^{-1} \mathbf e_j$. That is, direct backsolve by a vector containing a single 1 at entry $j$, then extract the $i$'th entry of the solution.

This would be pretty easy to code up just in Matlab. On a small problem (N=5K?) you could compare the approximate $\mathbf \Gamma'_{i,j}$ to the exact ($\sigma=0$) $\mathbf \Gamma_{i,j}$ as computed via the pseudoinverse. If they're pretty close, this might be all you need. (I expect they'd differ by about $\sigma \cdot \kappa$, with $\kappa$ the condition number).

• You can use the edge weights of a graph and then come with a measure of distance between graph nodes. An approach that appealed to me is resistance distance, which models the networks as if it were a electrical network. I have successfully leveraged a determinate based approach but it is too computationally slow to get values for all of the node pairs I am interested in en.wikipedia.org/wiki/Resistance_distance
– feik
Sep 13 '18 at 13:42
• Perhaps I jumped into the solution before laying out the problem, but in short I want to sample selected entries of the MPI in a efficient manner.
– feik
Sep 13 '18 at 13:50
• Very helpful. Do you have a handful of (i,j) pairs that you are trying to measure the distance between? Or maybe you are looking for pairs that have very large or very small distance? (I'd think it's the large distances that are interesting .. the (i,j) pairs with very short distances are basically encoded directly in the input adjacency matrix, right?). I admit, if you need the distance of every pair, you will be forced into this dense/explicit inverse strategy. But hopefully that's not the case. Sep 13 '18 at 13:51
• Does the accept indicate that this works? I'd be interested to hear your findings if you have a moment. Sep 13 '18 at 22:09
• Yes! (I think) I found a paper that supported a flavor of your recommendation $L^+ = (L + 1/n)^{-1} - 1/n$, where n is the dimension of the Laplacian. I was able to calculate the pseudoinverse in a a matter of hours. I am currently computing the psuedoinverse via svd for comparison but it is still cranking at the moment. Kron reduction also seemed like the next step to reduce the dimension of the problem(s).
– feik
Sep 17 '18 at 17:37

75k x 75k double-precision entries is 45 gigabytes. That fits in memory, but barely; you need to be careful.

The linear algebra routines in most languages rely on Lapack as a backend, which is a highly optimized linear algebra library written in Fortran. Most Lapack routines operate in-place: for instance, its QR routine xGEQRF does not allocate a second and a third matrix to store the results, but overwrites A with R, and then stores a compressed representation of Q in its lower triangular part (which would otherwise contain zeros). Typically, the "user-friendly" linear algebra packages like Numpy will make a copy of A, call xGEQRF, and then reconstruct Q and R as full matrices.

In your case, you have space in memory to store two matrices, but not three, so you have to avoid this. So Q, R = numpy.linalg.qr(A) won't work, but you have (barely) the space to make a Lapack call for the QR factorization, copy R to another memory location, and compute its product in-place with the implicitly-represented Q that is still in memory (there's another Lapack routine for that, xORMQR).

Moreover, computing the pseudoinverse with QR only works for full-rank matrices. If your matrix is rank-degenerate, you will have to work with the SVD, to which the same arguments apply. And, also, you will probably need to implement some form of regularization / truncation (for instance, truncated SVD, or Tikhonov, also known as ridge regression).

TL;DR: (1) don't use QR, use SVD (2) ditch the user-friendly libraries and use the Lapack interface directly.

• The following paper by Katsikis and Pappas describe a QR based approach to find the MPI: arxiv.org/pdf/1102.1845.pdf
– feik
Sep 13 '18 at 13:20
• @feik What they use is QR with column pivoting, not vanilla QR; it is a different factorization. It is a classical algorithm; it has been around since the 1960s (but it is still actual). As far as I know it is more prone to numerical problems when it comes to determining ranks, so it could be more troublesome to use for a non-full-rank matrix; see for instance Section 5.4.3 in Golub-Van Loan, Matrix computations, 4th edition. Sep 13 '18 at 13:53
• My matrix is a Laplacian which I believe is always rank n-1, if that can help with any shortcuts
– feik
Sep 13 '18 at 13:55