bugGNU Scientific Library - Bugs: bug #28767, gsl_linalg_SV_decomp gives NaN...

 
 

You are not allowed to post comments on this tracker with your current authentication level.

bug #28767: gsl_linalg_SV_decomp gives NaN when tb, tab, dt = 0

Submitter:  -Deleted Account- <bjg>
Submitted:  Fri 29 Jan 2010 10:04:47 PM UTC
   
 
Category:  Runtime error Severity:  3 - Normal
Operating System:  Status:  Fixed
Assigned to:  None Open/Closed:  Closed
Release:  1.13

Wed 24 Feb 2010 05:43:59 PM UTC, comment #2: 

fixed by commit 7946d551c30943f38bfcd1c31e109919e12fbe5a

    fix for bug #28767: gsl_linalg_SV_decomp gives NaN when tb, tab, dt = 0
   
    nans were generated in gsl_linalg_householder_transform from overflow
    in 1/(alpha-beta) when alpha=0 and beta is subnormal.
   
    nans were also generated in trailing_eigenvalue when dt=0 and tab=0.
    It was also possible to generate negative mu from cancellation error.
    Since we know the roots are non-negative and have the underlying
    expressions for them we can calculate them without cancellation error
    using only positive terms.
   
    in gsl_linalg_SV_decomp we can avoid a lot of numerical problems by
    rescaling the bidiagonal matrix before carrying out the implicit shift
    procedure. Following LAPACK's QR eigenvalue routine, we rescale if the
    max abs value A is in the range A > sqrt(DBL_MAX) or 0 < A < sqrt(DBL_MIN).
   
    added a test case for bug #28767, added check for nans when none
    present in input, added subnormal matrix that triggers potential
    overflow in gsl_linalg_householder_transform.


-Deleted Account- <bjg>
Fri 29 Jan 2010 10:05:28 PM UTC, comment #1: 

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

-Deleted Account- <bjg>
Fri 29 Jan 2010 10:04:47 PM UTC, original submission:  

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



-Deleted Account- <bjg>

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attached Files
file #19564:  svdbug.c added by bjg (1KiB - text/plain)

 

Depends on the following items: None found

Items that depend on this one: None found

 

CC list is empty

 

Follow 3 latest changes.

Date Changed by Updated Field Previous Value => Replaced by
2010-02-24 bjg StatusConfirmed Fixed
    Open/ClosedOpen Closed
2010-01-29 bjg Attached File- Added svdbug.c, #19564

Back to the top

Powered by Savane 3.13-02a9.
Corresponding source code