Move randomized algorithms

[haampie]

Maybe I'm totally biased on this, but I feel like randomized algorithms do not belong in this package and can better be moved to say RandomizedAlgorithms.jl.

These randomized methods seem to be mainly based on this review paper, which is in my humble opinion unjustifiably optimistic about its results and extraordinary pessimistic about 'classical methods'. Let me just rant about this article for a while; take this quote:

Any comparison between randomized sampling schemes and Krylov variants becomes complicated because of the fact that “basic” Krylov schemes such as Lanczos or Arnoldi are inherently unstable. To obtain numerical robustness, we must incorporate sophisticated modifications such as restarts, reorthogonalization procedures, etc. Constructing a high-quality implementation is sufficiently hard that the authors of a popular book on “numerical recipes” qualify their treatment of spectral computations as follows [...]

Obviously "restarts, reorthogonalization procedures, etc" are part of the method; the Arnoldi method is just the way people refer to IRAM. So they attack classical methods on being unstable, yet also admit there is a solution to that problem, and then just complain that the solution is too hard to implement correctly.

Next, they continue by admitting that their methods have the same problems:

Randomized sampling does not eliminate the difficulties referred to in this quotation; however it reduces the task of computing a partial spectral decomposition of a very large matrix to the task of computing a full decomposition of a small dense matrix.

The 'however' is not making any sense at all. Really?! That is what the randomized methods bring? This is literally what Henk van der Vorsts book on Krylov Subspace Methods mentions in the first chapter. The Hessenberg matrix in GMRES is exactly this type of small matrix.

Finally they agree they bring nothing new:

The comparison is further complicated by the fact that there is significant overlap between the two sets of ideas. Algorithm 4.3 is conceptually similar to a “block Lanczos method” with a random starting matrix. Indeed, we believe that there are significant opportunities for cross-fertilization in this area. Hybrid schemes that combine the best ideas from both fields may perform very well."

What's more, some of the algorithms they present are inherently unstable. Take the 'Prototype for Randomized SVD' algorithm. They basically say: "just do orthogonal subspace iteration, but skip the orthogonalization part. And by the way, if you want stability, actually don't forgot to do the orthogonalization step." I don't get why these things got published.

Maybe their objectives are completely different. They probably have extremely large and dense matrices and want to exploit BLAS3 all over the place. They probably just want 1 correct digit to a solution and can therefore afford unstable methods.

I think the objective of IterativeSolvers.jl are somewhat different: providing accurate solutions for problems mainly involving large sparse matrices. This might not be the case for least-squares problems, but still I think these randomized methods are not a great alternative.

[denis-bz]

Would you know of any real problems for *SVD on the web ? I'm skeptical about randomized SVD too (except to approximate random matrices :) but comparing some methods on 2 real problems might show us what works / what doesn't.

Thanks, cheers

[haampie]

Some people at Facebook seems to be fan of it: https://research.fb.com/fast-randomized-svd/

Other than that you can always generate a random matrix and judge for yourself.