[SciPy-user] sparse SVD

Rob rob.patro@gmail....
Wed May 6 22:46:49 CDT 2009


Kenneth Arnold wrote:
> On Wed, May 6, 2009 at 7:00 PM, Rob Patro <rob.patro@gmail.com> wrote:
>   
>> This looks very interesting indeed.  In particular, the matlab
>> implementation of sparse SVD has one feature of which I make extensive
>> use; this is the ability to request singular values around a spectral
>> shift, sigma, of my choosing.  Does SVDLIBC offer such an option?  How
>> difficult for me would it be to expose such an option in divisi?
>>     
>
> I have not heard of such an operation; what does that accomplish for
> you? (i.e., should we try it?)
>
> I'm not familiar with SVDLIBC at all beyond the part that we wrap. I
> looked at the docs a bit and didn't see an obvious option for doing
> what you describe, but that could be as much because I don't know the
> thing you're describing! SVDLIBC is simply a C port (not done by us)
> of the SVDLIB Fortran module.
>
> But if you did find it in SVDLIBC, the Cython (or Pyrex; I don't think
> we use anything cython-specific) wrappers make it quite easy to expose
> those things to Python.
>
> -Ken
> _______________________________________________
> SciPy-user mailing list
> SciPy-user@scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>   
Ken,

  Generally, such an operation is known as a spectral transform.  The
reason one might wish to
do this is as follows.  You wish to obtain enough singular value
triplets from the decomposition
of a large sparse matrix, M, to represent it to some particularly
desired precision.  In this case, computing
the entire decomposition is very  wasteful.  However,  one does not
know, a priori, how many singular triplets
will be required.  So you obtain a small number (say, 50) of the highest
energy singular triplets
.  If you can't represent M to the desired precision, you call the svd
routine again, but this time passing the
a spectral shift closest to the smallest singular value from the
previous iteration.  This returns 50 more singular
triplets which,  when  composed with the first 50, constitute the first
100 singular triplets of M; as if you had called
the SVD routine initially and requested 100 triplets.  You can iterate
this procedure, requesting a small number of triplets with a shifted
spectrum each time until you achieve your desired precision.

 In addition to not having to compute singular triplets you don't need,
this method usually has another benefit.  Often, the sparse eigenvalue
decomposition which power sparse SVD libraries are super-linear in the
number of requested eigenpairs.  Thus, it is often faster to request 50
eigenvectors (singular triplets) twice than it is to request 100
eigenvectors (singular triplets) once.  The situation only gets worse
when you are requesting many hundreds or thousands of triplets.  For
example, in our particular application, using the shift-invert spectral
transform to iteratively obtain singular triplets, we obtained more than
an order of magnitude speed up (we were doing SVDs on hundreds of
moderately sized sparse matrices).

--Rob


More information about the SciPy-user mailing list