Randomized Block Krylov Methods for Stronger and Faster Approximate Singular Value Decomposition
Since being analyzed by Rokhlin, Szlam, and Tygert [1]  and popularized by Halko, Martinsson, and Tropp [2], randomized Simultaneous Power Iteration has become the method of choice for approximate singular value decomposition. It is more accurate than simpler sketching algorithms, yet still converges quickly for any matrix, independently of singular value gaps. AfterÕ(1/ ) iterations, it gives a low-rank approximation within (1 + ) of optimal for spectral norm error. We give the first provable runtime improvement on Simultaneous Iteration: a randomized block Krylov method, closely related to the classic Block Lanczos algorithm, gives the same guarantees in justÕ(1/ √ ) iterations and performs substantially better experimentally. Our analysis is the first of a Krylov subspace method that does not depend on singular value gaps, which are unreliable in practice. Furthermore, while it is a simple accuracy benchmark, even (1 + ) error for spectral norm low-rank approximation does not imply that an algorithm returns high quality principal components, a major issue for data applications. We address this problem for the first time by showing that both Block Krylov Iteration and Simultaneous Iteration give nearly optimal PCA for any matrix. This result further justifies their strength over non-iterative sketching methods.
The authors discuss SVD approximation using a (randomized) Block Lanczos algorithm. The idea for the algorithm is not new (has been proposed before) but the author's claim of proven convergence guarantee is better than the (state-of-the-art) guarantee for the simultaneous power iteration.The proofs are relegated to the full version of the paper, and while that would not prevent acceptance (in the minds of this reviewer), when the major contribution is in these proofs (with concerns regarding the experiment, applicability, see later) then this makes the paper less suitable for NIPS.One compelling insight is that to compute a low rank approximation (classically done via a partial SVD) one may not necessarily need to produce exact (and orthogonal etc.) singular vectors (and then form the low rank approximator) and therefore the (numerical, orthogonality) issues arising from close singular values may not prevent the finding of a low rank approximation itself efficiently (and for any matrix). However, this is not how the authors communicated their message, case in point: low rank approximation does not get mentioned in the title and only indirectly in the abstract when referring very briefly to PCA.It is well known that the Lanczos (block Lanczos) algorithm suffers from numerical issues manifesting itself in losing orthogonality (in floating point arithmetic) that is guaranteed in exact arithmetic. A comprehensive read on this topic is "Lanczos Algorithms for Large Symmetric Eigenvalue Computations" by J.K. Cullum and R.A Willoughby (SIAM classics in applied mathematics, 41). Various degrees of reorthogonalization have been considered to address this issue, which is of course expensive and limit the number of eigenvalues that can be computed to the required accuracy. Paige in 1980 have shown what exactly is behind the loss of accuracy and it is not simply the result of roundoff error but the interaction between roundoff and convergence of some eigenvalues.The authors of this paper do not discuss the problem with loss of orthogonality and simply state they reorthogonalize. This is not adequate in this reviewer's opinion. Furthermore, the complexity computation does not consider the reorthogonalization work and therefore can be labeled as misleading. Despite the above this reviewer thinks that the effort to apply block Lanczos is important and there may be a way to achieve the complexity claimed by the authors, meaning not doing full reorthogonalization.It could be valuable to see the author's experiments regarding the implementation when no reortogonalization is used. Did they run into the numerical issues prompting the reorthogonalization?Section 2.2 the authors discuss stronger guarantees for low rank approximation. This is an interesting and valuable section with the introduction of per vector guarantee and the discussion of failures of prior low-rank approximation guarantees in practice.The authors augment their theoretical work with a valuable experiment section, that shows the improved convergence. However, as mentioned above, using iteration count hides the costs of reorthogonalization and therefore more insightful discussion is warranted to discuss practical applicability.The paper is well organized and the presentation is good. There are several typos and inconsistencies, listed below, that should be addressed by the authors to improve the menuscript.38: Why mention rank r here at all? In addition k is used for rank later (see line 72). This can just cause slight confusion.39: The authors (at least consistently) use a weird symbol to denote matrix transpose. Superscript T should work fine instead.114: This reviewer expected K to have \Pi as well as shown on line 273.119: This reviewer is confused why take the top k singular vectors of the k x k matrix M?159: "instead of multiplicatively like" The reviewer understands that the per vector guarantee uses \sigma_{k+1} and not ]\sigma_i but what does multiplicatively has to do with this?343, 373: There is more discussion need regarding reorthogonalization and its impact on computational cost.424: This was published in 2010.429: The author's first name should be: Tam\'as.The line numbers below correspond to the full version:507: Instead of just telling about the existence of p_1, the authors could give some motivation so that the reader knows somewhat why are we doing this. While this becomes clear later, it is good to give motivation to help the reader not to be puzzled.524: It is quite unexpected and (not well defined in this reviewer's opinion) to use x < O(y). This reviewer understands that the authors try to convey that there is a constant when moving to the last inequality. Using big O in this manner is casual at best.1070: This reviewer finds the formulation and proof a bit cumbersome. The use of \eta is not really necessary is it? One could say instead:     ||A-B||_2^2 - ||A-A_k||_2^2 <= ||A-B||_F^2 - ||A-A_k||_F^2The proof could be made more streamlined if continuing the above thought. This reviewer also suggests to recall (as it was mentioned before) that A_k is the best rank -k approximation and that ||A-A_k||_2^2 = \sigma_{k+1}^2, and similar for the Frobenius norm.Additional comments following the author response and reviewer discussions.I revised my review based on the discussions and the author responses. I change my rating to a 5.However, I still think that with a major revision this could be a very valuable paper. My driving concern is that once the decision is made, there is no incentive to make the changes necessary. (At least that was my experience last year when reviewing.) So I am keeping my score as a slightly below acceptance. (If this were a journal publication, my call would recommend acceptance with a major revision.)I implemented the algorithms and I believe that it has merit and could be useful in actual applications.I went through the complexity computations. Please stay with me so I can make my point. Instead of using \epsilon I find it much more informative if we look at the cost in terms of q (number of iterations) and the cost of the matrix vector computation for A (and A') that I will denote by m. (m = 2*nnz(A) which I would expect to be about c*n assuming a constant number of non-zeros per column.)Alg 1, Simultaneous iteration has cost breakdown: m(2*q+1)k    - to form K 2nk^2        - compute QR to get Q A'Q          - km form M       - dk^2 SVD of M     - 2k^3 form Z       - 2nk^2Assuming k << q < n, the dominating term will be forming of K.Alg 2, Block Lanczos has cost breakdown: m(2*q+1)k    - to form K (but also note much more memory needed) 2nq^2k^2     - compute QR to get Q A'Q          - qkm form M       - dq^2k^2 top k SVD of M - 2q^2k^3 form Z       - 2nqk^2(Small note: forming M can be done by applying A' then A and then Q instead, the term is not dominating so we do not discuss which is better.)The dominating term here appears to be the QR computation assuming m is about c*n.The achieved accuracy is C_SI* log(d)/q and (C_BL*log(d)/q)^2, where C_SI, C_BL are constants for the algorithms respectively. This shows however that the accuracy achieved per computational effort is not satisfactory:     effort               accuracy (guarantee) SI cn(2q+1)k            C_SI*log(d)/q BL 2nq^2k^2             (C_BL*log(d)/q)^2Therefore it would seem that using SI with q^2 iterations is better. Of course the above ignores that in practice performance depends on use of back storage and in general memory etc. (But I am not the one submitting the paper.)This however tells me that the approach to use QR to orthogonalize K is flawed. After all, the insight of the Lanczos method is in computing an orthogonal base efficiently. Here we use QR and not a recurrence. I for one even find the name Block Lanzos for this process misleading. It is a Krylov space based method, but no Lanczos recurrence about it as proposed. (However, this bad terminology could have gotten hold in the literature earlier and it is not the authors fault.)I believe there is an efficient way of producing a Q that is "good enough" (in presence of numerical issues) that would work in practice, but more work is needed. I note that the idea has been out there before using the Krylov space, the trick is to make it work efficiently.Some more comments below that would further improve the paper (in addition to the previous ones).l. 115: Consider defining K with going up to q-1 only. As is written K is an n x (q+1)*k matrix. Otherwise, q downstream needs to be updated to (q+1) for correctness.l 378: Compare the algorithms on a computational effort versus error plot as well as empirical time versus error plot. (The latter requires a faithfully good implementation of both algorithms, but not state-of-the art, just fair to the point being made.)Emphasize that A is accessible as an operator through its action (matrix-vector product) in the targeted applications and that k is intended to be very small compared to n.***The implementation of the algorithms and the test script for reference:SI.m: function [Z] = SI(A, k, q)    [n, d] = size(A);P = randn(d, k);B = A*A'; C = A*P; K = B*C; %keyboard;[Q, R] = qr(K); Qk = Q(:, 1:k);M = Qk'*B*Qk;[U, S, V] = svd(M);Uk = U(:, 1:k);Z = Qk*Uk; return;**** BL.m: function [Z] = BL(A, k, q)    [n, d] = size(A);P = randn(d, k);B = A*A'; C = A*P; K = zeros(n, q*k); for i = 0:q-1   K(:, i*k + 1:(i+1)*k) = C;   C = B*C; end[Q, R] = qr(K); Qk = Q(:, 1:q*k);M = Qk'*B*Qk;[U, S, V] = svd(M);Uk = U(:, 1:k);Z = Qk*Uk; return;test.m: *** rand("state", 423);n = 1000; d = 300; k = 5; tol = 1e-3; q = log(d)/sqrt(tol)*0.01 q = 10;A = randn(n, d); [U, S, V] = svd(A); Uk = U(:, 1:k); Ak = Uk*Uk'*A;Z1 = SI(A, k, q); Z2 = BL(A, k, q);norm(A - Z1*Z1'*A)/norm(A - Ak) norm(A - Z1*Z1'*A, 'fro')/norm(A - Ak, 'fro')norm(A - Z2*Z2'*A)/norm(A - Ak) norm(A - Z2*Z2'*A, 'fro')/norm(A - Ak, 'fro')