11
$\begingroup$

I'm working on a header-only matrix library to provide some reasonable degree of linear algebra capability in as simple a package as possible, and I'm trying to survey what the current state of the art is re: computing the SVD of a complex matrix.

I'm doing a two-phase decomposition, bidiagonalization followed by singular value computation. Right now I'm using the householder method for the bidiagonalization (I believe LAPACK uses this as well), and I think that's about as good as it gets currently (unless someone knows of an $\mathcal{O}(N^2)$ algorithm for it..).

The singular value computation is next on my list, and I'm somewhat out of the loop on what the common algorithms are for doing this. I read here that research was heading towards a inverse-iteration method that guarantees orthogonality with $\mathcal{O}(N)$ complexity. I'd be interested in hearing about that or other advances.

Improve this question
$\endgroup$
1
  • $\begingroup$ is there doc for your header-only matrix lib (apart from the .h) ? Also please add tag "svd". $\endgroup$
    – denis
    Jul 31, 2013 at 11:27

3 Answers 3

Reset to default
6
$\begingroup$

"Randomized algorithms" have recently become quite popular for partial svds. A header only implementation can be downloaded here: http://code.google.com/p/redsvd/

A review of the current methods can be found here: http://arxiv.org/abs/0909.4061

For full svds I am not sure if you can do better than Householder.

Improve this answer
$\endgroup$
6
  • $\begingroup$ That sounds very interesting, I'll have to take a look at that survey paper, thanks! $\endgroup$
    – gct
    Mar 19, 2012 at 14:13
  • $\begingroup$ OP is interested in algorithms for dense matrices. I don't think randomized algorithms are competitive in that setting, if they work at all. $\endgroup$
    – Federico Poloni
    May 15, 2015 at 19:06
  • $\begingroup$ This post indicated that randomized methods work just fine on dense matrices research.facebook.com/blog/294071574113354/fast-randomized-svd $\endgroup$
    – dranxo
    May 15, 2015 at 19:15
  • $\begingroup$ @dranxo There are no accuracy comparisons at all on that post, and the timing results do not look very meticulous. Also, randomized algorithms are based on projection + exact solution of a small-scale problem. This means OP would need an implementation of a "standard method" for the resulting small-scale problem anyway. $\endgroup$
    – Federico Poloni
    May 15, 2015 at 21:40
  • $\begingroup$ Fair enough. Though I am a bit confused why we should think that these methods only work on sparse matrices. Right out the abstract of Joel Tropp's paper: "For a dense input matrix, randomized algorithms require O(mn log(k)) floating-point operations (flops) in contrast with O(mnk) for classical algorithms." arxiv.org/pdf/0909.4061v2.pdf $\endgroup$
    – dranxo
    May 15, 2015 at 21:47
4
$\begingroup$

I read here that research was heading towards a inverse-iteration method that guarantees orthogonality with $\mathcal{O}(N)$ complexity. I'd be interested in hearing about that or other advances.

(I wanted to just make a few comments since I don't have the time to write out details, but it got rather big for the comment box.)

That I believe would be the MRRR (multiple relatively robust representations) algorithm of Dhillon and Parlett. This is rooted in previous work by Fernando, which in turn was inspired by a problem posed by Jim Wilkinson in his monumental book on eigenvalue problems. The "inverse iteration" portion for obtaining singular vectors is rooted in the concept of "twisted factorizations" by Fernando, which make use of factoring matrices into $\mathbf L\mathbf D\mathbf L^\top$ and $\mathbf U\mathbf D\mathbf U^\top$ decompositions.

The "singular value" portion of the algorithm, on the other hand, comes from the (shifted) differential quotient difference (dqd(s)) algorithm, which is a culmination of previous work by Fernando, Parlett, Demmel and Kahan (with inspiration from Heinz Rutishauser).

As you might know SVD methods usually proceed with a bidiagonal decomposition first before the singular values are obtained from the bidiagonal matrix. Unfortunately I'm not too updated on the current best method for the front-end bidiagonal decomposition; last I checked, the usual method is to start with QR decomposition with column pivoting and then apply orthogonal transformations appropriately to the triangular factor to finally obtain the bidiagonal decomposition.

I understand that I've been skimpy with the details; I'll try to flesh out this answer further once I have access to my library...

Improve this answer
$\endgroup$
2
  • $\begingroup$ Matrix to bi-diagonal form, do a column then a row, repeat down the diagonal: use givens or householder to zero the column up to the diagonal, then do the same for the row to the super-diagonal. $\endgroup$
    – adam W
    Oct 21, 2012 at 15:09
  • $\begingroup$ Ignore my previous comment, I was thinking tri-diagonal form. Sorry. Bi-diagonalization is not trivial in a similarity (it would actually reveal the eigenvalues), but your reference is not doing a similarity, it is doing something else; left and right multiply with different orthogonal matrices. $UAV$ is bi-diagonal with $UU^\top=I$ and $VV^\top=I$, which can be done as you say with QR first, but not so easily described in a comment. I might be interested if you do flesh out the answer further (but I could also figure it out since my studies are headed this direction at the moment). $\endgroup$
    – adam W
    Feb 25, 2013 at 16:28
0
$\begingroup$

There are the PROPACK and nu-TRLan libraries.

http://soi.stanford.edu/~rmunk/PROPACK/

http://crd-legacy.lbl.gov/~kewu/trlan/

Improve this answer
$\endgroup$
1
  • 2
    $\begingroup$ Here, the poster is asking for algorithms, rather than libraries; could you instead talk about the algorithms being used in these libraries, their computational complexity, and why these libraries are state-of-the-art? $\endgroup$
    – Geoff Oxberry
    Mar 19, 2012 at 10:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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