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:
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
tic
[Uk, Sk, Vk] = svds(A,k);
Ahat = Vk*diag(1./diag(Sk))*Uk';
test = Ahat * b;
time_k = toc
results = [results; k time_k];
end
plot(results(:,1), results(:,2))
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$