SVD bugs in gsl version .9.1?

SVD bugs in gsl version .9.1?

[Devin Balkcom]

Hi,

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

compile with:

gcc svd_test.c -lgsl -lgslcblas -lm

I'm using Redhat 6.2, but I've also had similar problems under Redhat
7.2 and Debian 2.2.
Hope you can help, thanks for a great piece of software!

-Devin

/************** svd_test.c *****************/

#include <gsl/gsl_linalg.h>
#include <stdio.h>

void vector_print(char *s, gsl_vector *v);
void matrix_print(char *s, gsl_matrix *M);

int main(char argc, char **argv) {

  gsl_matrix *A, *V, *X;
  gsl_vector *s, *work;


  A = gsl_matrix_alloc(3, 3);

  gsl_matrix_set_zero(A);

  gsl_matrix_set(A, 0, 0, -.212);
  gsl_matrix_set(A, 0, 1, .977);
  gsl_matrix_set(A, 0, 2, .155);

  matrix_print("A before SVD", A);

  V = gsl_matrix_alloc(3, 3);
  s = gsl_vector_alloc(3);

  X = gsl_matrix_alloc(3, 3);
  work = gsl_vector_alloc(3);

  /* does not return:
  gsl_linalg_SV_decomp(A, V, s, work);
  */

  /* also does not return:
  gsl_linalg_SV_decomp_mod(A, X, V, s, work);
  */

  /* this seems to work: */
  gsl_linalg_SV_decomp_jacobi(A, V, s);

  vector_print("s", s);
  matrix_print("A after SVD (U)", A);
  matrix_print("V", V);



  return 0;

}

void matrix_print(char *s, gsl_matrix *M) {
  int i, j;

  printf("%s (%d x %d):\n", s, M->size1, M->size2);
  for(i = 0; i < M->size1; i++) {
    for(j = 0; j < M->size2; j++) {
      printf("%7.3f ", gsl_matrix_get(M, i, j));
    }
    printf("\n");
  }
  printf("\n");

}

void vector_print(char *s, gsl_vector *v) {
  int i;
  printf("%s (%d-vector):  ", s, v->size);
  for(i = 0; i < v->size; i++) printf("%f ", gsl_vector_get(v, i));
  printf("\n\n");
}

Re: GSL for MSVC++ (Warning message report)

[Brian Gough]

Sun-Gi Hong writes:
 >  I tested the demo project included in the self-installing version
 > of gsl 0.9-1 for MSVC++. Of course, it works well, but there's a
 > problem when I changed the filename of "main.c" to "main.cpp".
 >  In the project setting, I checked "Disable language extensions"
 > but the following warning messages appear.
 >
 > _____________________________________________________________
 > main.cpp ..\include\gsl/gsl_vector_long_double.h(88) : warning
 > C4190: '<Unknown>' has C-linkage specified, but returns UDT
 > 'gsl_vector_long_double_const_view' which is incompatible with C
 >         ..\include\gsl/gsl_vector_long_double.h(55) : see
 >         declaration of 'gsl_vector_long_double_const_view'
 > ...  Linking...
 >  demo.exe - 0 error(s), 202 warning(s)
 > ______________________________________________________________
 >  Thank you for your consideration.
 >

Thanks for the bug report.. I've uploaded a new setup.exe which
includes the latest header files.  These should fix the problem.  If
not, let me know.

regards 

Brian Gough

Re: SVD bugs in gsl version .9.1?

[Brian Gough]

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;
>   }

Re: SVD bugs in gsl version .9.1?

[Brian Gough]

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

Thanks for the bug report.  

It is a bug in the SVD.... which is not handling the 2x2 case [0,X;0,0]
correctly, so it never converges for that matrix. I'll fix that.

regards
Brian Gough

Re: GSL for MSVC++ (Warning message report)

[Brian Gough]

For reference the size of the new file is

-rwxr--r--   1 bjg      bjg       3887274 Aug 30 10:38 /home/bjg/ftp/setup.exe

to make sure you get the correct version.