Inverse of large matrices (ILP64?)

[pmgurman]

Hi All,

I have an odd problem and I don't know if this is the correct place to ask, so apologies if not but here goes...

Some of the matrices that I have to invert are >50k per axis (i.e. 50k x 50k) and symmetric. The naive approach is np.linalg.inv(), but this means storing the square matrix. I know that LAPACK has inversion algorithms that work on the one triangle only, but once I exceed the limit of 32 bit integers to address, they crash (n (n+1) / 2 > 2^32-?). Here is some sample code and a crash message:

      1 def invp(B):
----> 2     Bp = sp.linalg.lapack.dtrttf(B)[0]
      3     Bpc = sp.linalg.lapack.dpftrf(B.shape[0], Bp)[0]
      4     Bpi = sp.linalg.lapack.dpftri(B.shape[0], Bpc)[0]
      5     Bi = sp.linalg.lapack.dtfttr(B.shape[0], Bpi)[0]

ValueError: failed to create intent(cache|hide)|optional array-- must have defined 
dimensions but got (-897458648,) 

I believe this is due to my scipy not using ILP64 but I don't know. I am using Anaconda and Scipy version 1.6.2.

This might help, if i try to force ILP64, I get the following.

sp.linalg.get_lapack_funcs('trttf', (B,), ilp64=True)

RuntimeError: LAPACK ILP64 routine requested, but Scipy compiled only with 32-bit BLAS

Note, I understand that it is best to avoid matrix inversion, however, the need still arises. Given the size of the matrices I am working with here, almost having the memory usage would help :). Any help or guidance as to where to seek answers is greatly appreciated.

[mdhaber]

Are you able to say what you want to do with the inverse? I've wanted to compile a list of reasonable uses of numerical matrix inverses.

[pmgurman]

In my case, the calculation can be written as (AA' + B)^-1 + C, with this result then added into a sub-block of another matrix. In my case, B and C are symmetric so I wanted to try to utilise the rectangular full packed format functions in LAPACK to reduce memory usage, but it seems to crash when the matrix requires 64 integers to address, which is why I was mentioning ILP64. Maybe the title is misleading?

@mdhaber, I agree with you that there are reasonable reasons to invert matrices explicitly, which mostly revolve around the need to do more complex calculations with the results of an inverse. While you could use, say, scipy.linalg.solve with columns of the identity matrix to calculate through clockwise columns, this would get messy for something like A^-1 A^-1.

[mdhaber]

Thanks!

i agree with you that there are reasonable reasons to invert matrices

Oh, i don't know firsthand or have my own opinion. I just always remember reading or hearing that it is almost never a good idea to invert, but those sources would never give examples of the exceptions to that rule.

Not sure that i can help here. @ilayn maybe?

[ilayn]

First the technical part, the standard build of SciPy if I am not mistaken is built with ILP64 disabled. Hence you need to use a BLAS implementation that uses ILP64. OpenBLAS can be built to use 64-bit integers with the INTERFACE64=1 flag. And in conda I think MKL has a similar flag. Then you can build SciPy with it and that error will disappear. You can then use the 64bit interface.

However even then, matrix inversion would be, I am guessing, useless since at these sizes the entries wildly swing due to numerical error build-up. You can let LAPACK do the hardwork with dsytri computing the inverse directly, it won't help with the precision but at least you don't need to shift LU representation around.

Instead as Matt said, you can avoid almost always the explicit inverses. But if you really have to, one obvious way is of course Schur complement which the typical matrix inversion lemma applies https://en.wikipedia.org/wiki/Schur_complement The idea is to reduce the size of the blocks that are inverted and carry out explicit expressions. But you have to somehow make sure that some inverse of the subblock exists

Another case is when B + AA' is pos def you can (and should) dpotri instead.

[dschmitz89]

Interesting that I stumble upon this now. I know that in practice very often large matrices are indeed inverted, especially in industry. Why? On one hand most people do not have a really deep understanding of numerical linear algebra (including me ;) ) and often enough it still works.

When it comes to really massive matrices, algorithms that process the matrix in parts such as A graph theoretic approach to matrix inversion by partitioning can be used successfully.

[matospiso]

50k per axis (i.e. 50k x 50k) and symmetric

I'm very late to reply but we have an industry ML use case for large-scale (10^6+) dimensional linear regression, and we compute approximate closed-form solution (approximate = we approximate the inverse)

In particular, we compute a sparse factorized inverse of the Gramian of a large, sparse matrix X in two steps:

1. compute incomplete Cholesky factorization of its regularized Gram matrix X^TX + aI
2. compute approximate inverse of the obtained lower triangular factor The computed inverse factor and its transpose then approximate a much denser inverse of the matrix X^TX + aI.

Check out https://github.com/glami/sansa/blob/main/src/sansa/model.py#L89-L114, you can extract this part to use in your computations.