* SVD bugs in gsl version .9.1?
@ 2001-12-19 13:20 Devin Balkcom
2001-12-19 13:20 ` Brian Gough
0 siblings, 1 reply; 5+ messages in thread
From: Devin Balkcom @ 2001-12-19 13:20 UTC (permalink / raw)
To: gsl-discuss
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");
}
^ permalink raw reply [flat|nested] 5+ messages in thread
* Re: GSL for MSVC++ (Warning message report)
[not found] ` <000b01c13115$f384b780$ad96f88f@sghong>
@ 2001-12-19 13:20 ` Brian Gough
2001-12-19 13:20 ` Brian Gough
0 siblings, 1 reply; 5+ messages in thread
From: Brian Gough @ 2001-12-19 13:20 UTC (permalink / raw)
To: Sun-Gi Hong; +Cc: gsl-discuss
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
^ permalink raw reply [flat|nested] 5+ messages in thread
* Re: SVD bugs in gsl version .9.1?
2001-12-19 13:20 ` Brian Gough
@ 2001-12-19 13:20 ` Brian Gough
[not found] ` <000b01c13115$f384b780$ad96f88f@sghong>
0 siblings, 1 reply; 5+ messages in thread
From: Brian Gough @ 2001-12-19 13:20 UTC (permalink / raw)
To: Devin Balkcom, gsl-discuss
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;
> }
^ permalink raw reply [flat|nested] 5+ messages in thread
* Re: SVD bugs in gsl version .9.1?
2001-12-19 13:20 SVD bugs in gsl version .9.1? Devin Balkcom
@ 2001-12-19 13:20 ` Brian Gough
2001-12-19 13:20 ` Brian Gough
0 siblings, 1 reply; 5+ messages in thread
From: Brian Gough @ 2001-12-19 13:20 UTC (permalink / raw)
To: Devin Balkcom; +Cc: gsl-discuss
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
^ permalink raw reply [flat|nested] 5+ messages in thread
* Re: GSL for MSVC++ (Warning message report)
2001-12-19 13:20 ` GSL for MSVC++ (Warning message report) Brian Gough
@ 2001-12-19 13:20 ` Brian Gough
0 siblings, 0 replies; 5+ messages in thread
From: Brian Gough @ 2001-12-19 13:20 UTC (permalink / raw)
To: Sun-Gi Hong, gsl-discuss
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.
^ permalink raw reply [flat|nested] 5+ messages in thread
end of thread, other threads:[~2001-12-19 13:20 UTC | newest]
Thread overview: 5+ messages (download: mbox.gz / follow: Atom feed)
-- links below jump to the message on this page --
2001-12-19 13:20 SVD bugs in gsl version .9.1? Devin Balkcom
2001-12-19 13:20 ` Brian Gough
2001-12-19 13:20 ` Brian Gough
[not found] ` <000b01c13115$f384b780$ad96f88f@sghong>
2001-12-19 13:20 ` GSL for MSVC++ (Warning message report) Brian Gough
2001-12-19 13:20 ` Brian Gough
This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox;
as well as URLs for read-only IMAP folder(s) and NNTP newsgroup(s).