SGESDD or DGESDD with huge matrices

Open discussion regarding features, bugs, issues, vendors, etc.

SGESDD or DGESDD with huge matrices

Postby flappinski » Thu May 21, 2009 9:15 am

Hi,
I am new to lapack and especially I am not a mathematician. So sorry for stupid questions.
I want to extract proncipal components out of a m-by-m matrix. Without lapack it took for weeks on relatively big datasets (1000x1000). So I introduced lapack, what improved the speed significantly (only hours).
But now I get a memory problem. with DGESDD the program crashes even before getting to this function (while initializing the variables) on a 16GB machine. So I switched to SGESDD. Since it uses float values, the initialization is fine. But I get an error, since the variable lwork is 2.805.682.880 for a matrix of 20020x20020. This figure is to high for an integer and gets a negative value.
I am using a 64bit machine.
Is anybody out there who can give me solution to this issue?
thanks a bunch,
Stephan
flappinski
 
Posts: 3
Joined: Thu May 21, 2009 9:00 am

Re: SGESDD or DGESDD with huge matrices

Postby Julien Langou » Thu May 21, 2009 10:12 am

That's really huge matrices you are speaking about there ... You are right there is a problem in this case. You have hit one of the few LAPACK routine that requires huge workspace. For the Divide and Conquer SVD routine, the required workspace is larger than the input matrix itself (in the square case). Our LWORK mechanism was not meant to handle that huge workspace with that huge matrices, in case of huge workspace (Divide and Conquer) and huge matrices, then you can potentially overflow the 32-bit INTEGER in which the computation are done for the workspace mechanism query (LWORK=-1). This is a good point. Moreover, assuming that you allocate the relevant workspace size, I have no idea how you can say it to LAPACK xGESDD with LWORK .... LWORK will overflow as well. Oups problem. If you can recompile LAPACK quickly by forcing all the INTEGER to be 64-bit by default (and not 32) that will fix the problem (I think, I do not know how to do this though, some compilers have a flag for this I believe).

So I am adding this to the bug list. Thanks for reporting this. See:
http://www.netlib.org/lapack/Errata/
(send me an email for last name and affiliation if relevant)

For information, in your case (JOBZ='A', M=N=20020) using the default ILAENV, the value of the required workspace size (to use the defualt block size) is: 1202541340. This is smaller (slightly) than 2^31, so I am not sure why you see some overflow at this level. Anyway, this is a good point.)

You wrote that LAPACK for an SVD of an 1,000-by-1,000 matrix was taking hours. You meant 10,000-by-10,000? Otherwise that's kind of slow, are you sure you are using an optimized BLAS library?

Best wishes,
Julien
Julien Langou
 
Posts: 727
Joined: Thu Dec 09, 2004 12:32 pm
Location: Denver, CO, USA

Re: SGESDD or DGESDD with huge matrices

Postby flappinski » Thu May 21, 2009 10:35 am

1) You are right, I meant 10.000 x 10.000.
2) I calculated lwork for 20020 after the old formula (3*m*m + 4*m*m + 4*m). Meanwhile, I found, that you corrected this, though, the old formula is still found in the net. but I will encounter very soon matrices with 40.000 entries.....
3) thanks for your fast and clarifying answer.
4) can you inform me, if you find a solution about this issue? (actually I don't know how to compile with 64bit Integers, I just use the shared library). Is there maybe another function, that finds principal components using less memory (I could spend more time, as long I don't get back to weeks)

Stephan Ripke
Broad Institute, Boston
flappinski
 
Posts: 3
Joined: Thu May 21, 2009 9:00 am

Re: SGESDD or DGESDD with huge matrices

Postby Julien Langou » Thu May 21, 2009 10:44 am

You can definetly give a try to DGESVD. It uses the Golub Kahan bidiagonal QR algorithm (while DGESDD uses Divide and Conquer on the bidiagonal matrix). The workspace will be of the order of M*NB where NB is the "optimal" block size. Way smaller than Divide and Conquer. The routine is in general slower though. Both xGESVD and xGESDD are O(n^3) method (due to bidiagonalization) so they have the same "order" in term of time.
Julien.
Julien Langou
 
Posts: 727
Joined: Thu Dec 09, 2004 12:32 pm
Location: Denver, CO, USA

Re: SGESDD or DGESDD with huge matrices

Postby Julien Langou » Thu May 21, 2009 10:52 am

By the way, just to be clear, to obtain "good" performance from LAPACK, you need to
(1) use the workspace size LWORK returned by a workspace query. So you do LWORK=-1, get optimal LWORK, allocate, call again. You do not use the formula given in the header of the LAPACK routine (e.g. for DGESVD that would be: LWORK >= MAX(1,3*MIN(M,N)+MAX(M,N),5*MIN(M,N))) If you use this that'll be slow.
(2) use an optimized BLAS library.
(1) and (2) gives "good" performance. If you do not have (1) or (2), you get bad performance.
--julien.
Julien Langou
 
Posts: 727
Joined: Thu Dec 09, 2004 12:32 pm
Location: Denver, CO, USA

Re: SGESDD or DGESDD with huge matrices

Postby flappinski » Thu May 21, 2009 3:44 pm

thanks.,
I first tried with the minimum workspace (since it is almost passing the integer threshhold). And even it is not optimal, it is still fast enough. i am absolutely happy with it.
I am only afraid of the even bigger datasets. The next one I already see is 35000 x 35000. So I will have to try DGESVD, I guess.
thanks a lot,
Stephan
flappinski
 
Posts: 3
Joined: Thu May 21, 2009 9:00 am


Return to User Discussion

Who is online

Users browsing this forum: No registered users and 1 guest

cron