Moore-Penrose pseudoinverse of singular rank degenerate matrix

Question

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?

Answer 1

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?

EDIT Given your broader explanation, here's what I would try:

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

Answer 2

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.