Scaling/Performance of Matlab's svds function (Lanczos bidiagonalization)

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 $5000times5000$, 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))

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

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

Parameters:
  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 =

   55.5512

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

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

Parameters:
  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 =

   45.7203