# What's the current state of the art regarding algorithms for the singular value decomposition?

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.

• is there doc for your header-only matrix lib (apart from the .h) ? Also please add tag "svd". Jul 31 '13 at 11:27

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.

• That sounds very interesting, I'll have to take a look at that survey paper, thanks!
– gct
Mar 19 '12 at 14:13
• OP is interested in algorithms for dense matrices. I don't think randomized algorithms are competitive in that setting, if they work at all. May 15 '15 at 19:06
• This post indicated that randomized methods work just fine on dense matrices research.facebook.com/blog/294071574113354/fast-randomized-svd May 15 '15 at 19:15
• @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. May 15 '15 at 21:40
• 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 May 15 '15 at 21:47

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...

• 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. Oct 21 '12 at 15:09
• 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). Feb 25 '13 at 16:28

There are the PROPACK and nu-TRLan libraries.

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

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

• 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? Mar 19 '12 at 10:20