This is the mail archive of the gsl-discuss@sources.redhat.com mailing list for the GSL project.


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]

Re: SVD bugs in gsl version .9.1?


Brian Gough writes:
 > Devin Balkcom writes:
 >  > I'm using gsl to take singular value decompositions.  I have had
 >  > problems with all three of the
 >  > implementations.   I have written a test program.
 >  > 
 >  > - gsl_linalg_SV_decomp() does not return for the test matrix.
 >  > - gsl_linalg_SV_decomp_mod() also does not return.
 >  > - gsl_linalg_SV_decomp_jacobi() does return, although I have had it fail
 >  > in other circumstances,
 >  >     with an incorrect result.  (I will send test code for this case
 >  > later.)

Ok now fixed, your program should work with the patch below or the
latest cvs version.  Thanks for the bug report... keep them coming.

regards
Brian Gough


Index: ChangeLog
===================================================================
RCS file: /cvs/gsl/gsl/linalg/ChangeLog,v
retrieving revision 1.30
diff -r1.30 ChangeLog
0a1,8
> Wed Aug 29 16:34:50 2001  Brian Gough  <bjg@network-theory.co.uk>
> 
> 	* svd.c (gsl_linalg_SV_decomp_jacobi): make sure all singular
>  	vectors are zero, not just first.
> 
> 	* svdstep.c (svd2): added explicit calculation of 2x2 svd, fixes
>  	bug that prevents convergence.
> 
Index: svd.c
===================================================================
RCS file: /cvs/gsl/gsl/linalg/svd.c,v
retrieving revision 1.16
diff -r1.16 svd.c
556c556,557
< 	    if (norm == 0 || (j > 0 && norm <= tolerance * prev_norm))
---
> 	    if (norm == 0.0 || prev_norm == 0.0 
>                 || (j > 0 && norm <= tolerance * prev_norm))
561c562
< 		prev_norm = 0;
---
> 		prev_norm = 0.0;
Index: svdstep.c
===================================================================
RCS file: /cvs/gsl/gsl/linalg/svdstep.c,v
retrieving revision 1.2
diff -r1.2 svdstep.c
85c85,114
< qrstep (gsl_vector * d, gsl_vector * f, gsl_matrix * U, gsl_matrix * V)
---
> create_schur (double d0, double f0, double d1, double * c, double * s)
> {
>   double apq = 2.0 * d0 * f0;
>   
>   if (apq != 0.0)
>     {
>       double t;
>       double tau = (f0*f0 + (d1 + d0)*(d1 - d0)) / apq;
>       
>       if (tau >= 0.0)
>         {
>           t = 1.0/(tau + hypot(1.0, tau));
>         }
>       else
>         {
>           t = -1.0/(-tau + hypot(1.0, tau));
>         }
> 
>       *c = 1.0 / hypot(1.0, t);
>       *s = t * (*c);
>     }
>   else
>     {
>       *c = 1.0;
>       *s = 0.0;
>     }
> }
> 
> static void
> svd2 (gsl_vector * d, gsl_vector * f, gsl_matrix * U, gsl_matrix * V)
86a116,118
>   size_t i;
>   double c, s, a11, a12, a21, a22;
> 
89,92d120
<   const size_t n = d->size;
<   double y, z;
<   double ak, bk, zk, ap, bp, aq, bq;
<   size_t i, k;
96c124
< 
---
>   
98d125
<   double f1 = (n > 2) ? gsl_vector_get (f, 1) : 0.0;
100c127,227
<   double mu = trailing_eigenvalue (d, f);
---
>   if (d0 == 0.0)
>     {
>       /* Eliminate off-diagonal element in [0,f1;0,d1] */
> 
>       create_givens (f0, d1, &c, &s);
> 
>       /* compute B <= G^T B */
> 
>       gsl_vector_set (f, 0, c * f0 - s * d1);
>       gsl_vector_set (d, 1, s * f0 + c * d1);
>       
>       /* Compute U <= U G */
> 
>       for (i = 0; i < M; i++)
> 	{
> 	  double Uip = gsl_matrix_get (U, i, 0);
> 	  double Uiq = gsl_matrix_get (U, i, 1);
> 	  gsl_matrix_set (U, i, 0, c * Uip - s * Uiq);
> 	  gsl_matrix_set (U, i, 1, s * Uip + c * Uiq);
> 	}
> 
>       return;
>     }
>   else if (d1 == 0.0)
>     {
>       /* Eliminate off-diagonal element in [d0,f1;0,0] */
> 
>       create_givens (d0, f0, &c, &s);
> 
>       /* compute B <= B G */
> 
>       gsl_vector_set (d, 0, d0 * c - f0 * s);
>       gsl_vector_set (f, 0, d0 * s + f0 * c);
> 
>       /* Compute V <= V G */
> 
>       for (i = 0; i < N; i++)
>         {
>           double Vip = gsl_matrix_get (V, i, 0);
>           double Viq = gsl_matrix_get (V, i, 1);
>           gsl_matrix_set (V, i, 0, c * Vip - s * Viq);
>           gsl_matrix_set (V, i, 1, s * Vip + c * Viq);
>         }
> 
>       return;
>     }
>   else
>     {
>       /* Make columns orthogonal, A = [d0, f1; 0, d1] * G */
>       
>       create_schur (d0, f0, d1, &c, &s);
>       
>       /* compute B <= B G */
>       
>       a11 = c * d0 - s * f0;
>       a21 = - s * d1;
>       
>       a12 = s * d0 + c * f0;
>       a22 = c * d1;
>       
>       /* Compute V <= V G */
>       
>       for (i = 0; i < N; i++)
>         {
>           double Vip = gsl_matrix_get (V, i, 0);
>           double Viq = gsl_matrix_get (V, i, 1);
>           gsl_matrix_set (V, i, 0, c * Vip - s * Viq);
>           gsl_matrix_set (V, i, 1, s * Vip + c * Viq);
>         }
>       
>       /* Eliminate off-diagonal elements, work with the column that
>          has the largest norm */
>       
>       if (fabs(a11) + fabs(a21) > fabs(a12) + fabs(a22))
>         {
>           create_givens (a11, a21, &c, &s);
>         } 
>       else
>         {
>           create_givens (a22, -a12, &c, &s);
>         }
>       
>       /* compute B <= G^T B */
>       
>       gsl_vector_set (d, 0, c * a11 - s * a21);
>       gsl_vector_set (f, 0, c * a12 - s * a22);
>       gsl_vector_set (d, 1, s * a12 + c * a22);
>       
>       /* Compute U <= U G */
>       
>       for (i = 0; i < M; i++)
>         {
>           double Uip = gsl_matrix_get (U, i, 0);
>           double Uiq = gsl_matrix_get (U, i, 1);
>           gsl_matrix_set (U, i, 0, c * Uip - s * Uiq);
>           gsl_matrix_set (U, i, 1, s * Uip + c * Uiq);
>         }
> 
>       return;
>     }
> }
102,103d228
<   y = d0 * d0 - mu;
<   z = d0 * f0;
105d229
<   /* Set up the recurrence for Givens rotations on a bidiagonal matrix */
107,111c231,239
<   ak = 0;
<   bk = 0;
< 
<   ap = d0;
<   bp = f0;
---
> static void
> qrstep (gsl_vector * d, gsl_vector * f, gsl_matrix * U, gsl_matrix * V)
> {
>   const size_t M = U->size1;
>   const size_t N = V->size1;
>   const size_t n = d->size;
>   double y, z;
>   double ak, bk, zk, ap, bp, aq, bq;
>   size_t i, k;
113,114c241,271
<   aq = d1;
<   bq = f1;
---
>   if (n == 2)
>     {
>       svd2 (d, f, U, V);
>       return;
>     }
> 
>   {
>     double d0 = gsl_vector_get (d, 0);
>     double f0 = gsl_vector_get (f, 0);
>     
>     double d1 = gsl_vector_get (d, 1);
>     double f1 = (n > 2) ? gsl_vector_get (f, 1) : 0.0;
>     
>     {
>       double mu = trailing_eigenvalue (d, f);
>     
>       y = d0 * d0 - mu;
>       z = d0 * f0;
>     }
>     
>     /* Set up the recurrence for Givens rotations on a bidiagonal matrix */
>     
>     ak = 0;
>     bk = 0;
>     
>     ap = d0;
>     bp = f0;
>     
>     aq = d1;
>     bq = f1;
>   }


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]