Also, since you computing the SVD of ATQ , which is presumably a small matrix even if A is large and sparse, wouldn’t you just use the dense SVD on this?
The whole idea came from the RSVDPACK paper, chapter 3.4. The claim is that for very large matrices,
eigen(A' * A) is much faster than
svd(A). And this is easy to check:
A = randn(100000, 1000);
@time eigen(A' * A);
20.607551 seconds (12 allocations: 1.527 GiB, 0.16% gc time)
2.657308 seconds (17 allocations: 30.877 MiB)
And it gets worse for larger matrices.
As for the conditioning problem that’s what they write:
A disadvantage of using this method is that the matrix BB^∗ has essentially the square of the condition number of B. As a result, very small singular values of A near machine precision may not be properly resolved. This is an issue only if A is expected to have very small singular values amongst σ_1, . . . , σ_k.
That’s a bit concerning, but at the end,
svd has really bad asymptotic from the matrix size, and we need to run it on
A' * Q, which has low number of columns, but as many rows as
A. Given that I personally need the method to work with matrices of >10M rows, using normal
svd is not an option…
Even for sparse matrices, ATA squares the condition number. I’m not sure why you consider KrylovKit a “heavy dependency” and not IterativeSolvers.jl?
IterativeSolvers just doesn’t have any real dependencies, and
KrylovKit has more of them, including
GPUArraysCore. Though it’s a minor thing.
How does your function compare to calling
KrylovKit.svdsolve directly on A ?
On the synthetic example
A = randn(100000, 1000) for 50 top values,
eigSVD returns the same eigenvalues as
eigsolve(A' * A, ...), but
eigSVD takes 2.5s,
eigsolve 1.7s, while
svdsolve takes 74s. Surprisingly,
svdsolve is much slower than full
However, on these tests I observed that
IterativeSolvers.lobpcg sometimes fails, while
eigsolve behaves much stabler. So I guess I’ll switch to
KrylovKit.eigsolve at the end.
Note also that your use of
randn means that you will automatically promote to
Float64. Better to use
randn(eltype(A), n, k+s) so that it uses the type of
While I see how it helps for consistency, I’m not sure that would be correct in this case. When we do multiplication of large matrices, it’s a lot of sums inside, which can lead to overflowing. As an example, below values in
A are in the range of 0-1000,
Q are standard normal, but
A' * Q consists fully of
Inf. And if we use low-precision floats to save memory, type promotion makes total sense to me:
A is large, so it’s
A' * Q are much smaller, so we can afford to make them
Float64. Same logic as with sparsity.
A = rand(Float16, 100000, 200) .* 1000;
Q = randn(eltype(A), size(A, 1), 10);
(A' * Q)[1:2, 1:5]
Inf Inf Inf -Inf Inf
Inf Inf Inf -Inf -Inf