Why is it that SVD routines on nearly square matrices run significantly faster than if the matrix was highly non square?

Question

In Python / Matlab, if you run a routine for SVD on a significantly non-square matrix, X, such as X.shape = (2,15000) you will get significantly longer run time than on a matrix with X.shape = (800,1000), even though 2*15000 = 30000, and 800*1000=800000.

That is the second matrix clearly has much more terms with which to iterate through than the first matrix yet it is computed faster in standard LAPACK computational routines.

I've tried to look through some big-O complexity analyses, since I believe LAPACK's dgeSVD subroutine is using a series of bi-diagonalization and Househoulder transformations:

https://www.netlib.org/lapack/explore-html/d1/d7e/group__double_g_esing_ga84fdf22a62b12ff364621e4713ce02f2.html

But it is confusing to know exactly what it is up to since in their recent paper, they point various small tweaks that have been made over the years which makes me forming a proper conclusion as to why I have my observation as above.

https://www.netlib.org/utk/people/JackDongarra/PAPERS/siam-svd-2018.pdf

I would've thought if matrices were no square, surely they would be placed under the same routine (for example). Then the matrix with fundamentally more parameters to iterate through should take longer than the "smaller" in terms of parameter count, matrix. Or it could be that they are indeed passed through the same subroutines but the difference in matrix size could be making a larger impact than I thought (I'm thinking too naive).

Any help on intuition / links / explanations would be greatly appreciated.

EDIT:

For reference on my PC:

 a = time.time()
X = np.random.randn(800,1000)
_,_,_ = np.linalg.svd(X)
b = time.time()
print(b-a)


a = time.time()
X = np.random.randn(1000,800)
_,_,_ = np.linalg.svd(X)
b = time.time()
print(b-a)


a = time.time()
X = np.random.randn(2,15000)
_,_,_ = np.linalg.svd(X)
b = time.time()
print(b-a)


a = time.time()
X = np.random.randn(15000,2)
_,_,_ = np.linalg.svd(X)
b = time.time()
print(b-a)

Has output (seconds):

0.6534378528594971

0.6089930534362793

2.464693069458008

2.5114269256591797

Answer 1

You are asking for a full (dense) SVD, which also needs to generate the unitary components of $U$ and $V$ which correspond with the null space of your input.

for the $1000 \times 800$ case, your input has $800000$ entries, while the output has $1000^2 + 800 + 800^2 = 1640800$ entries for a total of $2440800$ entries.

For the $15000 \times 2$ case, your input has $30000$ entries, while the output has $15000^2 + 2 + 2^2 = 225000006$ entries for a total of $225030006$ entries.

There is several orders of magnitude more work required to SVD decompose the highly rectangular matrix than the nearly square one.

If you can get away with a reduced SVD, you can pass full_matrices=False and now the amount of work required for the reduced SVD is significantly smaller for the $15000 \times 2$ matrix vs. the $1000 \times 800$ matrix.

Edit:

For the above analysis I just counted the number of entries in the input and output.

A full SVD has $U$ as an $M\times M$ matrix, $\Sigma$ is a diagonal matrix which has a maximum of $\min(M,N)$ non-zero entries along the diagonal, and $V$ is a $N \times N$ matrix.

In general, SVD algorithms usually follow the following mold (assume $M \ge N$):

1. Compute the eigendecomposition $A^* A = V \Sigma^* \Sigma V^*$. This step is $O(M N^2 + N^3)$.
2. Solve $U \Sigma = A V$ for $U$. This step is $O(M^3)$.

So our total complexity is $O(M^3 + M N^2 + N^3) = O(M^3)$. If we have the reverse situation where $M < N$, we can perform the process in a slightly different fashion to get the complexity $O(M^3 + M^2 N + N^3) = O(N^3)$.

See Trefethan and Bau lecture 31 for algorithms on how to compute SVDs.