DGESVD and machine precision

[gorcee]

Hello,

I have encountered an issue with the DGESVD routine in LAPACK that I don't quite know how to resolve.

I am attempting to compute the singular value decomposition of a matrix whose entries are small. I can compute the singular values in MATLAB with the svd() function, and I get correct answers. However, when using DGESVD in a MATLAB mex function, only the first singular value is correct.

I suspect that it has to do with the fact that the singular values are small and far apart. MATLAB's results for my test matrix give the following as the singular values:

6.725622826239752e-023
4.355711003542524e-040
8.444188487308763e-042
4.893839131447302e-042

Obviously, relatively speaking, the second value is under the traditional double-precision machine epsilon relative to the first value.

Is there a way to balance or scale the matrix to avoid this issue?

[sven]

As you yourself note, you cannot expect any correct figures in the small singular values since they satisfy sigma(i) < eps*sigma(1). By the way, the MATLAB svd function calls the LAPACK routine DGESVD (see 'doc svd').

As to whether or not you can scale depends upon your application. For example, if you are solving a full rank least squares problem, then you can column scale the matrix A since, if we put

r = b - A*x,

then

r = b - A*D*D^{-1}*x.

Thus we can minimize r = b - B*y, where

B = A*D and x = D*y.

So, we can choose D to be diagonal such that it makes the column norms of B roughly equal.

Sven Hammarling.

[gorcee]

Sven,

Thank you for your reply.

I was aware that MATLAB called DGESVD. I'm still curious as to why MATLAB generates the results it does whereas calling DGESVD gives different results. MATLAB must be pre-conditioning the matrix somehow, but I am uncertain as to how.

To provide some further detail, I am computing the SVD of a covariance matrix whose entries are small. The covariance matrix is formed by a Kalman filter innovation rho, such that C = cov(rho * rho^T). I am decorellating the residual by taking the svd of C = U * S * V^T, and saying that rho' = U^T * rho.

I am uncertain as to the overall effect these numerical discrepancies will have on the implementation of the algorithm I am developing. Nonetheless, it would be convenient to get close MATLAB's value, at least for aiding debugging.

Thanks again.

[Julien Langou]

Not sure what Matlab is doing behind the hood. (wood / hood ?)
There might be some pre-processing. No idea. Two easy possible differences:
* they might use a different BLAS library (should not have much impact for small matrices, will have for large, say more than 100)
* they might use DGESDD (as opposed to DGESVD)
* they might compute singular value only (as opposed to singular value and singular vectors)
--julien.

[gorcee]

Julien,

Hood is correct. Technically the American saying is "under the hood", but that's alright ;)

The MATLAB documentation indicates that they use DGESVD. But I did try using DGESDD instead. However, it doesn't give the same results, and is not really feasible for small matrices anyway (not to mention the limited memory available in an embedded application).

MATLAB is returning the singular vectors, but I am uncertain if they are doing this with more than one call to DGESVD. This is a possibility that I can try. It should be that with a symmetric matrix that the decomposition C = U S V^T that U = V, although the values coming out of MATLAB for U and V are just quite not the same. This is almost certainly an artifact of numerical precision.

Thank you again for your responses, it has given me a few avenues to look for the bugs in my code.

-Ed Gorcenski

[sven]

You don't say what your matrix C is, or what results DGESVD gives, but nevertheless if

sigma(i) < eps*sigma(1), i > 1,

then you cannot expect any correct figures in the sigma(i), i > 1. Indeed, numerically your matrix is singular. If MATLAB gives the 'correct' result then that is fortuitous. Results can differ due to different optimizations between compilers and compiler switches.

Do I understand that your matrix C is symmetric? If so, why not do a spectral factorization of C rather than an SVD? But, better still, if your matrix C is the form X * X^T, and you have X, it would be better to take the SVD of X,

X = U * S * V^T so that C = U * S^2 * U^T.

Since the singular values of X are the square roots of those of C, you would still have some accuracy in the small singular values.

Best wishes,

Sven.

[BirchRain]

Dear Sven

I'm a newcomer to use CLAPACK. I have a project need to use SVD filter. So I tried calculating SVD by using CLAPACK. However I found the results are great different from those in MATLAB. I called the routine dgesdd and sgesdd, but they gave similar results. I tried the example data in MATLAB svd help doc. A[4][2]= { 1, 2, 3, 4, 5, 6, 7, 8}. dgesdd gives the results of sv are: 14.2274, 1.2573, but MATLAB gives：14.2691, 0.6268. so I doubt why they are so different, Are the results believable? How could I improve the results precision?

It would be appreciated very much if you may give some help on this.

Bernard