Help needed with bdsqr

[asd__asdqwe]

Hello,

I'm trying to compute a partial SVD of a rectangular matrix A. I tried to adapt an example of the MKL that uses gesvd. While I get the same singular values, I'm not able to compute the correct left singular vectors. Any held would be greatly appreciated.

Thank you.

[mecej4]

In your C program you are using the CLAPACK interface, which was designed for use in C programs to make calls to Fortran Lapack routines compiled with the Fortran77-C translator f2c. In this interface, multidimensional arrays have to follow the Fortran column-major convention.

The results given by gesvd are correct and agree with those given by Matlab for the matrix

    8.79   6.11  -9.15   9.57  -3.49   9.84 
    9.93   6.91  -7.93   1.64   4.02   0.15 
    9.83   5.04   4.86   8.83   9.80  -8.99 

In the first part of your code, you have computed the SVD using a different set of routines, namely, dgebrd_ and dbdsqr_. These routines return P^T * VT and U * Q rather than VT and U. It is not clear to me whether you cared about the singular vectors returned by gebrd/bdsqr, or whether you asked only about the singular vectors returned by gesvd.

[asd__asdqwe]

Thanks for your answer.

I completely agree that the results with the second set of routines is right (gesvd), it is the MKL example. However, I would like to compute the singular vectors with bdsqr, and there, I can't make it work. Unfortunatly, there is no example with bdsqr inside the mkl/examples folder.

So in short, how to fix the first part of the code to return the same two first singular vectors ?

Thank you.

[mecej4]

Please read the documentation for the routines that you wish to use. The array u must be initialized to a unit matrix before calling the subroutine bdsqr, if nru is not zero. In your code you set equal to zero the elements of array What do you think?, but failed to set the diagonal elements to 1. Similarly, the array vt must be initialized to a unit matrix if ncvt is not zero.

[asd__asdqwe]

Right, sorry about that. I now compute the array u using dorgbr to get the singular vectors of the original system. I'm still not quite there, but at least it is not all zeros. If you have any clue on what I might be doing wrong, please let me know.

Thank you for your help.

[mecej4]

You can always set up the computation to obtain the min(m,n) left singular vectors and discard all but the first column of the resulting matrix U (whose columns are the left singular vectors), if you call the ?gesvd set of routines. I do not know offhand what to do when you want only one singular vector and you want to obtain it at minimum expense. I have some ideas, but they are untested.

Situations where you do not want to discard results obtained with significant CPU time spent are probably going to involve large, perhaps sparse, matrices. For such matrices, you would probably not use the dense matrix routines of Lapack. At any rate, time spent searching the Web and libraries for materials on solving for individual eigenvectors and singular vectors would be a good investment.

You might, on the other hand, have a need for finding singular vectors of a sequence of thousands or millions of small or medium dense matrices. The selection of routines for this purpose would probably be different.

Is it out of question for you to compute min(m,n) left singular vectors, keep the first and discard the rest? Perhaps you could expand on how the singular vector calculation plays a role in the overall computation that you have in mind. You gave the impression that you wanted to use the computational routines ?gebrd, ?orgbr, ?bdsqr rather than the driver routine ?gesvd. If so, why?

[mecej4]

Here is a corrected and shortened version of your program. See also http://www.nag.com/numeric/CL/nagdoc_cl24/examples/source/f08kfce.c and http://www.nag.com/numeric/fl/nagdoc_fl23/xhtml/F08/f08mef.xml. The MKL version of the latter is provided as  .../MKL/examples/lapack/source/dorgbrx.f with an Intel compiler installation, after the examples zip file examples_core_f.zip has been extracted. For some reason, no corresponding example is included in examples_core_c.zip.

[Added on 1-11-2015] If you intend to call Lapack from C, you may find it advantageous to use the LAPACKE interface, which is a bit more C-like and relieves you from having to pass workspace arrays as arguments. The second attachment, dgebrdx.cpp, is the same program as the previous one, modified to use this interface.

[asd__asdqwe]

It could be possible to keep the first one(s), but the idea is to do SVD on 100.000x10.000 matrices, and only keep ~100 singular vectors, so I would like to be able to only keep compute 100 singular vectors right away. I thought it was possible with ?bdsqr.

[mecej4]

I do not think that with direct methods you can compute selected singular values/vectors. In fact, the manual page for ?bdsqr says, "u(ldu,*) contains an nru by n unit matrix U. The second dimension of u must be at least max(1, n)." Thus, the choice as far as singular vectors are concerned is "all or none". 

Is your 10000 X 10000 matrix sparse or dense?

Please take a look at ARPACK and ARPACK-ng; the example driver program at https://github.com/opencollab/arpack-ng/blob/master/EXAMPLES/SVD/dsvd.f appears to do exactly what you seek.

[asd__asdqwe]

My matrix is dense, so I guess there is a better alternative than ARPACK. Anyway, thank you for your help.