According to the man page of svds, provided by MATLAB, svds is currently based on "Augmented Lanczos Bidiagonalization Algorithm" when it comes to the top-k SVD computation of a large-scale sparse matrix, instead of using the sampling technique answered by loup blanc (you can easily check it by typing edit svds in your MATLAB prompt). Please refer to the following paper:

 - Baglama, James, and Lothar Reichel. "Augmented implicitly restarted Lanczos bidiagonalization methods." SIAM Journal on Scientific Computing 27.1 (2005): 19-42.

This is considered as an anytime iterative algorithm, i.e., it iteratively computes and updates the target top-k singular triplets until convergence. See Algorithm 3.1 of the above paper. 

BTW, it is painful to strictly analyze the time complexity of Algorithm 3.1, since the algorithm is not that intuitive to capture the whole procedure and which part is the main bottleneck of the algorithm. 

At a glance, it is considered as O(T(|A|k + k^3 + c)) where T is # of iterations, |A| is # of non-zeros in the input sparse matrix A, k is the target number of the largest singular values, and c is the other computational cost for each step. Note that this might be wrong since it is roughly estimated based on sparse matrix multiplication and QR decomposition.