# How to compute Singular value decomposition of a large matrix with Python

Language: Python3

Problem: I have a matrix Q of shape [51200 rows x 51200 cols] stored in a binary file, each of the element in this matrix has a data type of complex64. To load the data into memory I therefore need ~20GB of RAM, checked from getsizeof(Q). I do have access to a server with 120GB of RAM in a LINUX machine.

My aim is to decompose the matrix with SVD.

1. The easiest way in Python to do this is by using np.linalg.svd(Q). To do this, I first use np.fromfile() to load the Q, and then execute the svd function. The problem here is, I do not know, how much memory I exactly need to compute this function. And I do get a warning init_zgesdd failed init. Though this does not stop the computation, but in the end, the values for U,S,V* are all zeros. As checked this warning is due to memory allocation.

2. A second approach I tried is by using scipy.sparse.linalg.svds library. Since there are a lot of zeros (about 20%), I thought defining the matrix as sparse would have better memory usage. I found that while running this, the consumption of memory fluctuates from 50GB to 100GB, but it gets killed after running about 15-20 min.

3. I have also looked into on how can I decrease the precision of the matrix element. As of now I am using complex64 (that is 32 bit float for each real and imag part), I am not seeing any option for making it complex32.

I wanted to know the best way to compute SVD for such matrix.

• In (1), it looks like Python isn't checking the return values of the Lapack functions it calls. – Federico Poloni Dec 15 '20 at 20:14
• Another problem with integer overflow could be related to the LWORK parameter of cgesdd(). LWORK is roughly 52000^2, which is slightly larger than the largest possible 32 bit signed integer. Can you scale down to a slightly smaller problem (say 40000 by 40000) and see what happens? – Brian Borchers Dec 15 '20 at 21:26
• @BrianBorchers I have rescaled the matrix to 40K x 40K. And using np.linalg.svd(Q). Now I do not get any error/warning message, but the process keeps on running (waited for 20 min) with max RAM usage at 84GB (out of allocated 120GB). With the top command, I can see the process is consuming the resources. When I use svds(Q,k=10) so to get the first 10 singular values with corresponding eigen vectors, it gives the result after 10 min or so. – SAM Dec 16 '20 at 14:57
• It's pretty clear that your python is using a LAPACK library with 32 bit integers rather than 64 bit integers. If you get to a point where you really need all the singular values and vectors then you should find a python distribution that links to an optimized LAPACK/BLAS library with support for 64 bit integer parameters. Computing the SVD of a 50,000 by 50,000 matrix should take no more than about 5 hours (based on 150 seconds to do a 10,000 by 10,000 matrix on my 3 year old desktop and scaling up by (50000/10000)^3. Your machine is probably considerably faster.) – Brian Borchers Dec 16 '20 at 19:54
• if you're prepared to write a bit more difficult code, SCALAPACK might be worth giving a shot. That code is more suited to be run on larger compute servers. – Thijs Steel Dec 22 '20 at 19:25