11
$\begingroup$

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.

$\endgroup$
11
  • 1
    $\begingroup$ In (1), it looks like Python isn't checking the return values of the Lapack functions it calls. $\endgroup$ – Federico Poloni Dec 15 '20 at 20:14
  • 6
    $\begingroup$ 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? $\endgroup$ – Brian Borchers Dec 15 '20 at 21:26
  • 1
    $\begingroup$ @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. $\endgroup$ – SAM Dec 16 '20 at 14:57
  • 2
    $\begingroup$ 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.) $\endgroup$ – Brian Borchers Dec 16 '20 at 19:54
  • 2
    $\begingroup$ 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. $\endgroup$ – Thijs Steel Dec 22 '20 at 19:25
0
$\begingroup$

Did you try with Dask ?

https://examples.dask.org/machine-learning/svd.html

You can manage very large matrices. There is also a nice blog post about it https://blog.dask.org/2020/05/13/large-svds

$\endgroup$
1
  • 2
    $\begingroup$ Unless I am misunderstanding the links suggest dask is good for square matrices in the case where you only need a few singular values rather than when you need the full decomposition - which might be the case here, the question doesn't say, but I think it worth flagging up $\endgroup$ – Ian Bush Dec 23 '20 at 9:47

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.