6
$\begingroup$

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.

$\endgroup$
  • 2
    $\begingroup$ This is a quite reasonable question to ask of the MATLAB developers, but not an appropriate one for this group- you're asking questions about the internal design of MATLAB and the rational for a particular choice (not transposing the matrix) that can really only be explained by someone involved in writing MATLAB. $\endgroup$ – Brian Borchers Mar 30 '15 at 4:12
  • 2
    $\begingroup$ One possible reason for not taking the transpose is that there might not be enough available memory for $X$, $X^{T}$, and the SVD. Also note that MATLAB uses LAPACK routines internally. You could test whether this same performance hit occurs in a C or Fortran program that calls the LAPACK library. $\endgroup$ – Brian Borchers Mar 30 '15 at 4:16
  • 6
    $\begingroup$ 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 colum is cheap and thats done in the implementation. Furthermore having a huge number of row, these operations take advantage out of the threading of the BLAS backend ( which is the MKL in the case of MATLAB). $\endgroup$ – M.K. aka Grisu Mar 30 '15 at 19:40
  • $\begingroup$ @BrianBorchers, good point about the memory. I suppose they can't assume you have the memory available. Grisu, yes, that's why I mentioned the stride-length being 1. $\endgroup$ – Stephen Mar 30 '15 at 20:35
  • 1
    $\begingroup$ To back Grisu's suggestion a bit more, here are timings obtained from NumPy's SVD (linked against OpenBLAS): a SVD for a C-contiguous matrix needs 0.71s, whereas the (still C-contiguous) transpose needs 1.23s. Changing storage order to Fortran storage, the picture reverses. $\endgroup$ – AlexE Mar 31 '15 at 13:30
4
$\begingroup$

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.

$\endgroup$
  • 3
    $\begingroup$ This may explain why the SVD implementation used by MATLAB is slower on short fat matrices, but it does not answer the OP's question "why isn't SVD of fat matrices in MATLAB implemented as the faster version [V,S,U]=svd(X')". $\endgroup$ – Federico Poloni Apr 1 '15 at 18:49
2
$\begingroup$

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?".)

$\endgroup$
0
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Yes, but @Federico Poloni's comment to Grisu's post still applies $\endgroup$ – Stephen Apr 1 '16 at 6:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.