bug with dgesdd

Open discussion regarding features, bugs, issues, vendors, etc.

bug with dgesdd

Postby Bruno » Sat Jun 19, 2010 10:39 am

Hi,

in the octave bug list there is a matrix example where dgesdd performs badly
(see http://savannah.gnu.org/bugs/?func=deta ... 87#options).
I don't know if lapack team is aware of this problem (at least I didn't see it
on this forum).

I have transformed the matrix in a ascii format and written a small f90 program
which reads it and use dgesdd or dgesvd but I am unable to attach
the 2 files (what is possible to attach to a message for this forum ?) So I
have put the 2 files in my homepage :

http://www.iecn.u-nancy.fr/~pincon/

On the given matrix :

1/ dgesvd works quite well, I got :

|| U*diag(sigma)*VT - A ||_1 / || A ||_1 near 10^(-14)

2/ but with dgesdd, first the last singular value is negative !!! and
the previous relative error is about 48.4 !

I have tested using 2 differents blas lib and got the same results.

hth
Bruno
Bruno
 
Posts: 9
Joined: Fri Mar 07, 2008 3:41 am

Re: bug with dgesdd

Postby admin » Mon Jun 21, 2010 12:10 pm

Thanks Bruno!
We actually knew that there was a serious bug in dgesdd but we couldn't find a matrix to reproduce it.
For the current problems found in LAPACK, you can always refer to http://www.netlib.org/lapack/Errata/
The dgesdd problem is bug0025 as you can see.
The good news is that I was able to reproduce it and so I confirm the problem.
Your code and your matrix are really great and straightforward, so we will start debugging the dgesdd routine.
We keep you posted.
Merci beaucoup
Julie et Julien Langou
admin
Site Admin
 
Posts: 468
Joined: Wed Dec 08, 2004 7:07 pm

Re: bug with dgesdd

Postby Bruno » Mon Jun 21, 2010 2:46 pm

The matrix has been posted on the octave bug tracker
by Jaroslav Hajek, I have only translated it in ascii and
wrote the small f90 code.

Thanks for the ref bug,

Bruno
Bruno
 
Posts: 9
Joined: Fri Mar 07, 2008 3:41 am

Re: bug with dgesdd

Postby Julien Langou » Mon Jun 21, 2010 4:56 pm

Hello,

Thanks for the bug report!

So I have taken a look at all this.

Good news: I think I have found where the bug is.
Bad news: I think the routine should have failed with INFO = 1.
So does not really help. But that helps somewhat!

So line 353 from 374 of DBDSC reads:
Code: Select all
            IF( ICOMPQ.EQ.2 ) THEN
               CALL DLASD0( NSIZE, SQRE, D( START ), E( START ),
     $                      U( START, START ), LDU, VT( START, START ),
     $                      LDVT, SMLSIZ, IWORK, WORK( WSTART ), INFO )
            ELSE
               CALL DLASDA( ICOMPQ, SMLSIZ, NSIZE, SQRE, D( START ),
     $                      E( START ), Q( START+( IU+QSTART-2 )*N ), N,
     $                      Q( START+( IVT+QSTART-2 )*N ),
     $                      IQ( START+K*N ), Q( START+( DIFL+QSTART-2 )*
     $                      N ), Q( START+( DIFR+QSTART-2 )*N ),
     $                      Q( START+( Z+QSTART-2 )*N ),
     $                      Q( START+( POLES+QSTART-2 )*N ),
     $                      IQ( START+GIVPTR*N ), IQ( START+GIVCOL*N ),
     $                      N, IQ( START+PERM*N ),
     $                      Q( START+( GIVNUM+QSTART-2 )*N ),
     $                      Q( START+( IC+QSTART-2 )*N ),
     $                      Q( START+( IS+QSTART-2 )*N ),
     $                      WORK( WSTART ), IWORK, INFO )
               IF( INFO.NE.0 ) THEN
                  RETURN
               END IF
            END IF

and I strongly believe the check on INFO should be outside the ELSE statement.

In your case, DBDSDC is calling DLASD0 twice (well with vecLib not with reference BLAS !!!), the first call returns an INFO = 1 so everything should stop. The second call returns INFO = 0 and overrides the first INFO. This is wrong.

So I'll send this to the lapackers and hopefully we will patch this soon.

Julien
Julien Langou
 
Posts: 727
Joined: Thu Dec 09, 2004 12:32 pm
Location: Denver, CO, USA

Re: bug with dgesdd

Postby Bruno » Fri Jun 25, 2010 8:17 am

Thanks Julien. I have another question : when dgesdd "normally" outputs
with info=0, can we considered that the singular values and vectors have
near about the same accuray than with dgesvd ? (is there some paper
about the stability of dgessd ?)

Cordiales salutations
Bruno
Bruno
 
Posts: 9
Joined: Fri Mar 07, 2008 3:41 am

Re: bug with dgesdd

Postby Julien Langou » Thu Jul 01, 2010 2:57 am

Hello, update: Last week, I forwarded your question around because I do not know the answer for sure. I did not receive any answer yet. Bottom line there should be not much difference in term of the quality of the answer. The first step for both is a reduction to bidiagonal form. The second step whether you use GESDD or GESVD has about the same accuracy or more accurate than the first. So at the end, you'll get the same accuracy more or less. Cheers, Julien.
Julien Langou
 
Posts: 727
Joined: Thu Dec 09, 2004 12:32 pm
Location: Denver, CO, USA

Re: bug with dgesdd

Postby Bruno » Mon Jul 12, 2010 3:56 am

Thanks for the information Juien (and sorry for the delay I suppose
I forgot the dgesdd pb...).

Finally this means that for a singular elements computations one can
try first dgesdd and then dgesvd in case the first fails.

Regards
Bruno
Bruno
 
Posts: 9
Joined: Fri Mar 07, 2008 3:41 am

Confused with square root of matrix and dgesdd results

Postby vukicevic » Tue Aug 03, 2010 7:19 pm

Hi. I'm trying to find square root of the folowing matrix:

A =
(9,-16)
(-16,3)

I use dgesdd function, and I've got:

U =
(-0.769509, 0.638636)
(0.638636, 0.769509)

S =
(22.2788, 0)
(0, 10.2788)

VT =
(-0.769509, 0.638636)
(-0.638636, -0.769509)

Until now everything seems ok, because

U * S * VT = A =
(9,-16)
(-16,3)

Now it should be possible to get square root of A simply by taking square roots of eigenvalues in matrix S:

S**(1/2) =
(4.72004, 0)
(0, 3.20606)

and square root of A should be:

A**(1/2) = U * S**(1/2) * VT

But when I had tried to check the result, I've got:

A**(1/2) * A**(1/2) =
(17.3846, -5.89723)
(-5.89723, 15.1731)

which is obviously wrong.

I'm a litle confused with that result because the same procedure works fine with e.g.

B =
(2.2,0.4)
(0.4,2.8)

I suppose that I can try to solve the problem with dsyevr function, but it is interesting to know what is wrong.

Best regards, Vlatko
vukicevic
 
Posts: 1
Joined: Tue Aug 03, 2010 6:38 pm

Re: bug with dgesdd

Postby Julien Langou » Wed Aug 04, 2010 2:35 pm

You will not be able to compute the square root of a matrix using an SVD decomposition!
Julien Langou
 
Posts: 727
Joined: Thu Dec 09, 2004 12:32 pm
Location: Denver, CO, USA


Return to User Discussion

Who is online

Users browsing this forum: Google [Bot] and 1 guest

cron