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, 2013 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, 2013 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, 2013 at 15:05

4 Answers 4

Reset to default
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, 2013 at 10:05
  • 1
    $\begingroup$ Looks like scipy has wrappers. I'll add them to answer body. $\endgroup$
    – Max Hutchinson
    Jun 10, 2013 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, 2013 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.