[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Bug-gsl] Re: GSL: NaN results in SVD calculation (gsl_linalg_SV_decomp)
From: |
Bruno Daniel |
Subject: |
[Bug-gsl] Re: GSL: NaN results in SVD calculation (gsl_linalg_SV_decomp) |
Date: |
Tue, 26 Jan 2010 12:48:18 +0100 |
--------------------------------------------------------------
GSL version number: 1.9
hardware and operating system, OS details:
Linux phit1 2.6.25.20-0.5-default #1 SMP 2009-08-14 01:48:11 +0200 x86_64
x86_64 x86_64 GNU/Linux
gcc version:
gcc version 4.3.1 20080507 (prerelease) [gcc-4_3-branch revision 135036]
(SUSE Linux)
genre: linalg, SVD, gsl_linalg_SV_decomp
description of the bug behavior: NaN results in gsl_linalg_SV_decomp
short program which exercises the bug: included in this E-Mail
--------------------------------------------------------------
Dear GSL developers,
I managed to reproduce the data with a randomly generated matrix of
appropriate statistics: In this 462x421 matrix, about 90% of the rows and 43%
of the columns are completely zero. Of the other elements, 60% are 1 and the
rest is zero. Here's the test program and the output. Most of the time the
result contains NaNs; sometimes GSL's SVD algorithm fails to converge:
----- svdbug.c --------------------------------------------------
/*
GSL version: 1.9
Compile and run:
gcc `gsl-config --cflags` svdbug.c `gsl-config --libs` -o svdbug && ./svdbug
*/
#include <stdio.h>
#include <stdlib.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#define rand_double() ((double)(rand()) / RAND_MAX)
int main(int argc, char** argv) {
int m = 462, n = 421;
gsl_matrix* A = gsl_matrix_alloc(m, n);
srand(time(NULL));
int i, j;
int row_completely_zero[462];
int column_completely_zero[421];
for (i = 0; i < m; i++) {
row_completely_zero[i] = rand_double() < 0.9;
}
for (i = 0; i < n; i++) {
column_completely_zero[i] = rand_double() < 0.43;
}
int n_ones = 0, n_zeros = 0;
for (i = 0; i < m; i++) {
for (j = 0; j < n; j++) {
double a = row_completely_zero[i] == 1 || column_completely_zero[j] == 1
? 0 : (rand_double() < 0.6 ? 1 : 0);
gsl_matrix_set(A, i, j, a);
if (a == 0) { n_zeros++; }
else if (a == 1) { n_ones++; }
}
}
printf("m = %d, n = %d, %d elements, %d zeros, %d ones, rest: %d\n", m, n,
m * n, n_zeros, n_ones, m * n - n_zeros - n_ones);
gsl_matrix* V = gsl_matrix_alloc(n, n);
gsl_vector* S = gsl_vector_alloc(n);
gsl_vector* work = gsl_vector_alloc(n);
gsl_linalg_SV_decomp(A, V, S, work);
printf("first 10 singular values:\n");
for (i = 0; i < 10; i++) {
if (i > 0) { printf(", "); }
printf("%g", gsl_vector_get(S, i));
}
return 0;
}
---- output ----------------------------------------------------
(phit1 ~/pj) gcc `gsl-config --cflags` svdbug.c `gsl-config --libs` -
o svdbug && ./svdbug
m = 462, n = 421, 194502 elements, 189295 zeros, 5207 ones, rest: 0
first 10 singular values:
0, 0, nan, nan, nan, nan, nan, nan, nan, nan
(phit1 ~/pj) gcc `gsl-config --cflags` svdbug.c `gsl-config --libs` -
o svdbug && ./svdbug
m = 462, n = 421, 194502 elements, 188300 zeros, 6202 ones, rest: 0
first 10 singular values:
0, nan, nan, nan, nan, nan, nan, nan, nan, nan
(phit1 ~/pj) gcc `gsl-config --cflags` svdbug.c `gsl-config --libs` -
o svdbug && ./svdbug
m = 462, n = 421, 194502 elements, 189691 zeros, 4811 ones, rest: 0
first 10 singular values:
nan, nan, nan, nan, nan, nan, nan, nan, nan, nan
(phit1 ~/pj) gcc `gsl-config --cflags` svdbug.c `gsl-config --libs` -
o svdbug && ./svdbug
m = 462, n = 421, 194502 elements, 188698 zeros, 5804 ones, rest: 0
gsl: svd.c:149: ERROR: SVD decomposition failed to converge
Default GSL error handler invoked.
[1] 29225 abort ./svdbug
(phit1 ~/pj)
Kind regards
Bruno Daniel
Bruno Daniel wrote:
> --------------------------------------------------------------
> GSL version number: 1.9
> hardware and operating system, OS details:
> Linux phit1 2.6.25.20-0.5-default #1 SMP 2009-08-14 01:48:11 +0200 x86_64
> x86_64 x86_64 GNU/Linux
> gcc version:
> gcc version 4.3.1 20080507 (prerelease) [gcc-4_3-branch revision 135036]
> (SUSE Linux)
> genre: linalg, SVD, gsl_linalg_SV_decomp
> description of the bug behavior: NaN results in gsl_linalg_SV_decomp
> short program which exercises the bug: attached
> --------------------------------------------------------------
>
> Dear GSL devolopers,
>
> your library is a great help for me and I use it very much. Thank you very
> much for your great work!
>
> Recently, I encountered a problem in the SVD routine of GSL version 1.9:
> Certain matrices A with very many zeros lead to an SVD with NaN ("not a
> number") results in the vector of singular values S as well as in the left and
> right singular vectors.
>
> I tracked the error down to the function trailing_eigenvalue in svdstep.c:
> If tb, tab and dt are all zero, the calculated eigenvalue is undefined:
>
> mu = tb - (tab * tab) / (dt + hypot (dt, tab)) = 0 - 0 / 0;
>
> This NaN result is then propagated in the function qrstep via create_givens to
> all the other results.
>
> I attach a test program reproducing the problem when applied to a certain
> matrix of 462 rows and 421 columns i.e. 194502 elements of which 189694 are
> zeros and 4808 are ones, saved in a file of 389 kbyte. The distribution of
> ones in the matrix is not random. I can send this matrix file to whoever is
> interested (This mailing list refuses attachments of this size).
>
> Here's the output of this program on my machine with the following
> architecture:
>
> Linux phit1 2.6.25.20-0.5-default #1 SMP 2009-08-14 01:48:11 +0200 x86_64
> x86_64
> x86_64 GNU/Linux
>
> --------------------------------------------------------------
>
> (phit1 ~/pj) gcc `gsl-config --cflags` svdbug.c `gsl-config --libs` -o
> svdbug &&
> ./svdbug
> 000000000000000000000000000000000000000000000000000000000000000000000000000000000
> 000000000000000000000000000000000000000000000000000000000000000000000000000000000
> 000000000000000000000000000000000000000000000000000000000000000000000000000000000
> 000000000000000000000000000000000000000000000000000000000000000000000000000000000
> 000000000000000000000000000000000000000000000000000000000000000000000000000000000
> 0000000000000000
> 000000000000000000000000000000000000000000000000000000000100010001000000000000000
> 010000000000000000000000000000000000000000000000000000000000000000000000000000000
> 000000001000000111000000000000000000000100000000000000100000000000000000000000000
> 000000000000000000000000000000000000000000000000000000000000000000000000000000000
> 000000000000000000000000000000000000000000000000000000000000000000000000000000000
> 0000000000000000
> 000000000000000000000000000000000000000000000000000000000000001000000000000000000
> 000000000000000000000000000000000000000000000000100000000000000001000000000000000
> 000000000000011000000000000000000000000000000000000000000000000000000000000000000
> 000000100000000000000000000000000000000000000000000000000000000000000000000000000
> 000000000000000000000000000000000000000000000000000000000000000000000000000000000
> 0000000000000000
> 000000100000000000000000000000000000000000000000000000000000001000100000000000000
> 000000000000100000000000000000000000000000000000000000000000000001000000000000000
> 000000000000011000000000000000000000000000000000000000000000000000000000000000000
> 000000000000000000000000000000000000000000000000000000000000000000000000000000000
> 000000000000000000000000000000000000000000000000000000000000000000000000000000000
> 0000000000000000
> 000000000000000000000000000000000000000000000000000001000000001000000000000000000
> 000000000000000000000000000000000000000000000000000000000000000001000000000000000
> 000000000000011000000000000000000000000000000000000000000000000000000000000000000
> 001000010000000000000000000000000000000000000000000000000000000000000000000000000
> 000000000000000000000000000000000000000000000000000000000000000000000000000000000
> 0000000000000000
> 000000000000000100000000000000000000000000000000000000000000100000000000000000000
> 000000000000000000000000000000000000000000000000100000000000000000000000000000000
> 000000000011100000000000000000000000000000000000000010000000000000000000000000000
> 000000000000000000000000000000000000000000000000000000000000000000000000000000000
> 000000000000000000000000000000000000000000000000000000000000000000000000000000000
> 0000000000000000
> 000000000000000100000000000000000000000000000001000000000000100000000000000000000
> 000000000000000000000000000000000000000000000000000000000000000000000000000000000
> 010000000011100000100000000000000000000000000000000100000000000000000000000000000
> 000000000000010000000000000000000000000000000000000000000000000000000000000000000
> 000000000000000000000000000000000000000000000000000000000000000000000000000000000
> 0000000000000000
> 000000000000000100000000000000000000000000000001000000000000100000000000000000100
> 000000000000000000000000000000000000000000000000000000000000000000000000000000000
> 010000000011100000100000000000000000000000000000000100000000000000000000000000000
> 000000000000000000000000000000000000000000000000000000000000000000000000000000000
> 000000000000000000000000000000000000000000000000000000000000000000000000000000000
> 0000000000000000
> 000000000000000000000000000000000000000000000000000000011000010001000000000000000
> 000000000000000000000000000000000000000000000000000000000000000000000000000000000
> 001000001000000000100000000000000000110100000000000101000100010000000000000000100
> 000000000000000000000000000000000000000000000000000000000000000000000000000000000
> 000000000000000000000000000000000000000000000000000000000000000000000000000000000
> 0000000000000000
> 000000000000000000000000000000000000010000000001000000000000100000000000000000000
> 001000000000000000000000000001000000000000000000000000000000000000000000000000000
> 010000000111100000001000000001000000000000000000001000010000000000000000000000000
> 000000000000010000000000000000000000000000000000000000000000000000000000000000000
> 000000000000000000000000000000010000000000000000000000000000000000000000000000010
> 0000000000000000
> m = 462, n = 421, 189694 zeros, 4808 ones, 194502 elements
> first 10 singular values:
> 0, 0, nan, nan, nan, nan, nan, nan, nan, nan
> (phit1 ~/pj)
> --------------------------------------------------------------
>
> Kind regards,
> Bruno Daniel
>
>
> /*
> GSL version: 1.9
>
> Compile and run:
> gcc `gsl-config --cflags` svdbug.c `gsl-config --libs` -o svdbug &&
> ./svdbug
> */
> #include <stdio.h>
> #include <gsl/gsl_matrix.h>
> #include <gsl/gsl_vector.h>
>
> int main(int argc, char** argv) {
> int m = 462, n = 421;
> gsl_matrix* A = gsl_matrix_alloc(m, n);
>
> // Read the matrix A from the file svdbug.c, show the first 10 lines and
> report
> // the number of zeros and ones.
> FILE* fp = fopen("svdbug.csv", "r");
> int i, j, n_ones = 0, n_zeros = 0;
> for (i = 0; i < m; i++) {
> for (j = 0; j < n; j++) {
> double a;
> fscanf(fp, "%lf", &a);
> gsl_matrix_set(A, i, j, a);
> if (i < 10) { printf("%g", a); }
> if (a == 0) { n_zeros++; }
> else if (a == 1) { n_ones++; }
> }
> if (i < 10) { printf("\n"); }
> }
> fclose(fp);
> printf("m = %d, n = %d, %d zeros, %d ones, %d elements\n", m, n, n_zeros,
> n_ones, m * n);
>
> gsl_matrix* V = gsl_matrix_alloc(n, n);
> gsl_vector* S = gsl_vector_alloc(n);
> gsl_vector* work = gsl_vector_alloc(n);
> gsl_linalg_SV_decomp(A, V, S, work);
> printf("first 10 singular values:\n");
> for (i = 0; i < 10; i++) {
> if (i > 0) { printf(", "); }
> printf("%g", gsl_vector_get(S, i));
> }
> return 0;
> }