From mboxrd@z Thu Jan 1 00:00:00 1970 From: Brian Gough To: Devin Balkcom , gsl-discuss@sources.redhat.com Subject: Re: SVD bugs in gsl version .9.1? Date: Wed, 19 Dec 2001 13:20:00 -0000 Message-id: <15245.3443.41485.347734@debian> References: <3B8BE5C9.FA9DE2F@ri.cmu.edu> <15244.5395.985031.513235@debian> X-SW-Source: 2001/msg00439.html 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 > > * 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; > }