10
$\begingroup$

For model reduction, I want to compute the left singular vectors associated to the - say 20 - largest singular values of a matrix $A \in \mathbb R^{N,k}$, where $N\approx 10^6$ and $k\approx 10^3$. Unfortunately, my matrix $A$ will be dense without any structure.

If I just call the svd routine from the numpy.linalg module in Python for a random matrix of this size, I run into a memory error. This is due to the allocation of $V\in \mathbb R^{N,N}$ for the decomposition $A = VSU$.

Are there algorithms around, that avoid this pitfall? E.g. by setting up only the singular vectors assiociated with nonzero singular values.

I am ready to trade in computation time and accuracy.

Improve this question
$\endgroup$
3
  • 1
    $\begingroup$ Interesting, it seems Numpy does not know how to do a thin SVD... $\endgroup$ – J. M. Jun 7 '13 at 12:27
  • $\begingroup$ Thanks for the hint. Indeed, numpy.linalg.svd has the option full_matrices that be set to False so that only the 'nonzero' parts are computed. Nevertheless, is there a way to reduce the computation even further? $\endgroup$ – Jan Jun 7 '13 at 12:42
  • 3
    $\begingroup$ The numpy backend uses fortran code, the LAPACKE_dgesvd routine for standard svd. However, typically your matrix is C_CONTIGOUS (check with matrix.flags). Therefore it copies the data for fortran alignment. Additionally while running the lapack routine dgesvd another copy of your matrix is needed (or at least the memory for it). You can get rid of one copy if you make sure that the memory alignment is fortran style right from the beginning. $\endgroup$ – Bort Jun 7 '13 at 15:05
6
$\begingroup$

If you only want a few singular values/vectors, ARPACK should do the trick. The SVD docs aren't great, and this distribution is more up to date.

EDIT: If you want to do this in python, SciPy has a wrapper. Since your matrix is dense, you could try the block sparse row (BSR) format.

Improve this answer
$\endgroup$
2
  • $\begingroup$ I will have a look, how ARPACK integrates with python... $\endgroup$ – Jan Jun 10 '13 at 10:05
  • 1
    $\begingroup$ Looks like scipy has wrappers. I'll add them to answer body. $\endgroup$ – Max Hutchinson Jun 10 '13 at 18:46
2
$\begingroup$

Take a look at sklearn.decomposition.TruncatedSVD in scikit-learn 0.14-rc.
(I believe that the scikit-learn people follow stackoverflow.com/questions/tagged/scikit-learn, so I'd ask detailed questions there.)

(How much memory do you have ? 10$^{6+3}$ doubles is already 8G.)

Improve this answer
$\endgroup$
1
  • $\begingroup$ Thanks for your answer. By now, I do well with the scipy routines. Also, I haven't gone till $10^6 \times 10^3$ yet but to about the half of it what is still feasible for my laptop. If necessary, I can use a working machine with 32 GB RAM. $\endgroup$ – Jan Aug 10 '13 at 11:17
2
$\begingroup$

Maybe you can try this.

https://github.com/jakevdp/pypropack

This is a Python wrapper for the PROPACK package, which implements efficient partial singular value decompositions of large sparse matrices and linear operators.

Improve this answer
$\endgroup$
2
$\begingroup$

Intel MKL implements the new Jacobi-SVD algorithm. Here are the implemention details: http://www.netlib.org/lapack/lawnspdf/lawn169.pdf http://www.fernuni-hagen.de/MATHPHYS/veselic/downloads/j02.pdf

And the LAPACK routine: http://software.intel.com/sites/products/documentation/hpc/mkl/mklman/GUID-732F9EE1-BCEC-4D9B-9B93-AF5499B21140.htm#DRMAC08-1

Work size is of course adjustable. You could call C functions from Python easily using Cython, SWIG or any other wrapping mechanism.

Improve this answer
$\endgroup$

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.