Why is Matlab's SVD faster on skinny matrices than on fat matrices?

Question

I noticed something odd today. I have a matrix X that is very skinny (20800 x 200), double precision real numbers, not sparse, and I want the SVD of it quickly. Matlab does it rather fast:

> tic; [U,S,V] = svd(X,'econ'); toc
Elapsed time is 0.280848 seconds.

But if I ask for the SVD of its transpose, which is an extremely fat matrix, it is considerably slower.

> Xt = X';
> tic; [UU,SS,VV] = svd(Xt,'econ'); toc
Elapsed time is 0.722308 seconds.

Any ideas why this is? It seems quite odd, since if I wanted the SVD of a fat matrix, this means I can do it faster by taking the transpose, finding the SVD of this skinny matrix, and then swapping the "U"s and "V"s.

My guess is that it's because Matlab uses column major order, so in the skinny case, the "U" matrix is the large one and whatever routine that operates on it gets to use a stride-length of 1, whereas in the opposite case, the Matlab implementation calls something with a non-unit stride-length which is less efficient in terms of memory calls.

But even if there is a good reason, it begs the question, why doesn't Matlab check for fat matrices and just take the SVD of the transpose? The transpose operator is blazing fast. e.g.,

> tic; [VV,SS,UU] = svd(Xt','econ'); toc
Elapsed time is 0.293725 seconds.

gives me the same VV,UU,SS as above, but much faster.

Answer 1

As I already mentioned in the comments here is a possible answer which is backed by some experiments from AlexE in a further comment.

The SVD in MATLAB uses the DGESVD from LAPACK, which is based on the ideas of Gene Golub. The main thing is that is is implemented on matrices in Fortran, i.e. columwise storage. In this way processing values in the same colunm is cheap in the sense of memory transfers because the prefetch system of the CPU detects a continuous data stream and elements from the same column are most likely available to the CPU before they are requested by the algorithm. Regarding the Fortran storage scheme for matrices and how the algorithm is implemented, accessing elements in a column is more efficient as accessing elements in a row. Furthermore having a huge number of rows, these operations take advantage out of the threading of the BLAS backend (which is the MKL in the case of MATLAB).

Applied to the problem of computing the SVD of tall and skinny matrix many operations take advantages of the above mentioned techniques to accelerate the computation. If one have the transpose of such a matrix the data access scheme is not longer that regular and the number of expensive row accesses increase. Furthermore, the number of rows is too small to accelerate computations by additional parallelization in a column.

Answer 2

In my opinion, and with my limited knowledge of the Matlab internals, your understanding is correct. Matlab uses a suboptimal algorithm, and I don't see other plausible reasons than "they didn't think of that" or "they didn't care optimizing it".

(Another related question to which I don't have a good answer is "why doesn't LAPACK's SVD routine DGESVD accept a TRANS parameter, unlike several others?".)

Answer 3

A common trick to accelerate SVD of skinny matrices is to firstly apply a QR decomposition to it and than you do a SVD in the smaller matrix R. This gives you A = QR = QSVD, where QS is orthogonal, then QSVD is an SVD decomposition. Maybe, MATLAB does check the shape of the matrix before applying the SVD and use some similar trick.