I have a simple Matlab script which aims to compute $k$ singular values of a matrix $A$. $A$ is a random dense square matrix of size $5000\times5000$, with 100 of its singular values constrained to be 0 (though that last detail does not seem to matter for my question).

I'm doing this in Matlab via [Uk, Sk, Vk] = svds(A, k);. According to the documentation, svds uses Lanczos bidiagonalization to compute these values. I looked at the function definition (edit svds) and do not see any relevant branching, e.g. using different algorithms under the hood based on different conditions. However, when I increase $k$, I get very curious scaling/performance:

Scaling of svds as a function of the number of singular values k

The docs mention

Increasing k can sometimes improve performance, especially when the matrix has repeated singular values.

But I interpret this to mean performance would be improved per $k$, rather than some huge reduction in total overall runtime.

Is this a known behavior of Lanczos bidiagonalization (an algorithm I'm not very familiar with)? Or does anyone have any speculations as to why the performance of svds is like this?

Edit: Here is a minimal version of my script so others can try to reproduce:

results = [];
A = rand(5000, 5000);
[U, S, V] = svd(A);
dS = diag(S);
dS(4900:5000) = 0;
A = U*diag(dS)*V;
b = rand(5000, 1);
for k = 100 : 100 : 4500
        [Uk, Sk, Vk] = svds(A,k);
        Ahat = Vk*diag(1./diag(Sk))*Uk';
        test = Ahat * b;
    time_k = toc
    results = [results; k time_k];
plot(results(:,1), results(:,2))
  • $\begingroup$ Can you post a minimal version of your script to see if others can reproduce? $\endgroup$ – Spencer Bryngelson Aug 2 '20 at 3:12
  • 1
    $\begingroup$ @SpencerBryngelson Added! And tried this snippet on a separate machine and see similar behavior... $\endgroup$ – davewy Aug 2 '20 at 8:13
  • 1
    $\begingroup$ Have you checked the results of the SVD are correct? $\endgroup$ – Ian Bush Aug 2 '20 at 8:39
  • 2
    $\begingroup$ I don 't have my work laptop with me so i can't check, but it would surprise me if this wasnt because of using a different algorithm. These timings would be explained by the matlab switching to a dense SVD (note that the timings stay constant). At the peak, you are asking for more than 20% of the singular values. At that point using a dense SVD is very reasonable. $\endgroup$ – Thijs Steel Aug 2 '20 at 14:41
  • $\begingroup$ What @ThijsSteel said. It seems you're using svds wrong; it is not designed to compute 30% of the singular values. It is very likely that the branching is inside the libraries, buried inside Fortran code. $\endgroup$ – Federico Poloni Aug 2 '20 at 18:43

I was able to reproduce your initial result via the snippet. However, by adding some more options to your svd call:

[Uk, Sk, Vk] = svds(A,k,'largest','display',true);

we see that the algorithm indeed changes to a dense one (@ThijsSteel).

For 300 singular values:

=== Singular value decomposition A*v = sigma*u, A'*u = sigma*v ===

Computing 300 largest singular values of 5000-by-5000 matrix A.

  Maximum number of iterations: 100
  Tolerance: 1e-10
  Subspace Dimension: 900

Find largest singular values for A*v = sigma*u, A'*u = sigma * v.

--- Start of Lanczos bidiagonalization method ---
Iteration   1: 144 of 300 singular values converged. Smallest non-converged residual 1.4e-09 (tolerance 1.0e-10).
Iteration   2: 300 of 300 singular values converged.
To check if singular value multiplicities were missed, restart the method, looking for k+1 singular values.
Iteration   3: 301 of 301 singular values converged.
No additional multiple singular values found. Successful return.

time_k =


For 2300 singular values:

=== Singular value decomposition A*v = sigma*u, A'*u = sigma*v ===

Computing 2300 largest singular values of 5000-by-5000 matrix A.

  Maximum number of iterations: 100
  Tolerance: 1e-10
  Subspace Dimension: 6900

Compute SVDS by calling SVD, because the subspace dimension is equal to the minimum matrix size.

time_k =

  • $\begingroup$ Thank you for showing how to get the debug output which makes the behavior clear here. I guess I'm a little surprised dense SVD is so much better than svds for even moderate numbers of singular values (e.g. 500 out of 5000). svds was advertised as the "fast" way to get a subset of singular values/vectors (scicomp.stackexchange.com/questions/35453/…). Is there no faster algorithm available (in Matlab or elsewhere) that beats or matches dense SVD when e.g. $k < .9*n$? $\endgroup$ – davewy Aug 2 '20 at 21:02
  • 1
    $\begingroup$ No, 90% is way too much to be using a krylov method. Your matrix is also dense, not sparse. If the dense SVD is too slow for you, you can look at GPU implementations, those are pretty fast for SVD (since the bidiagonalization phase is memory bound). $\endgroup$ – Thijs Steel Aug 3 '20 at 6:42
  • $\begingroup$ Right. My goal is to just use whatever is considered best / state of the art (from an algorithmic perspective) for estimating a (non-tiny, but < 100%) percentage of the singular components of a dense $A$. If the community feels the answer is dense SVD (whatever implementation) then that is good enough for me :) $\endgroup$ – davewy Aug 3 '20 at 8:38

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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