This is the mail archive of the
gsl-discuss@sources.redhat.com
mailing list for the GSL project.
Problem with Singular Value Decomposition Algorithm
- To: <gsl-discuss at sources dot redhat dot com>
- Subject: Problem with Singular Value Decomposition Algorithm
- From: "Jim Love" <Jim dot Love at asml dot com>
- Date: Wed, 12 Sep 2001 10:30:29 -0400
I have downloaded the latest beta release and the SVD algorithm produces the wrong answers. It appears that columns are swapped in the output and possibly a sign problem.
Here was my test array:
1 1 0.975
1 -1 0.975
-1 -1 -0.925
-1 1 -1.025
The correct answer for the S vector is: 2.7940 2.0000 0.0358
The gls output was: 2.0000 2.7940 0.0358
The Correct Q matrix is:
-0.7155 0.0256 -0.6981
0.0183 0.9997 0.0179
-0.6983 -0.0000 0.7158
The gls output was:
-0.025633 -0.715538 -0.698103
-0.999671 0.018347 0.017900
-0.000000 -0.698332 0.715774
A similar problem was seen in the U matrix.
Any ideas? Is this caused by my implementation or is it a real bug?
Here is the source code to the program I made to test the LIB:
#include "stdafx.h"
#include <math.h>
#include <stdio.h>
#include <C:\Program Files\gsl-0-9-2\gsl\gsl_matrix.h>
#include <C:\Program Files\gsl-0-9-2\gsl\gsl_blas.h>
#include <C:\Program Files\gsl-0-9-2\gsl\gsl_linalg.h>
int main(void)
{
gsl_matrix *a = gsl_matrix_alloc (4, 3);
gsl_matrix *u = gsl_matrix_alloc (4, 4);
gsl_matrix *v = gsl_matrix_alloc (3, 3);
gsl_vector *s = gsl_vector_alloc (3);
gsl_vector *work = gsl_vector_alloc (3);
gsl_matrix_set_zero(a);
gsl_matrix_set (a, 0, 0, 1);
gsl_matrix_set (a, 0, 1, 1);
gsl_matrix_set (a, 0, 2, 0.975);
gsl_matrix_set (a, 1, 0, 1);
gsl_matrix_set (a, 1, 1, -1);
gsl_matrix_set (a, 1, 2, 0.975);
gsl_matrix_set (a, 2, 0, -1);
gsl_matrix_set (a, 2, 1, -1);
gsl_matrix_set (a, 2, 2, -0.925);
gsl_matrix_set (a, 3, 0, -1);
gsl_matrix_set (a, 3, 1, 1);
gsl_matrix_set (a, 3, 2, -1.025);
printf ("%f %f %f \n",gsl_matrix_get (a, 0, 0), gsl_matrix_get (a, 0, 1), gsl_matrix_get (a, 0, 2));
printf ("%f %f %f \n",gsl_matrix_get (a, 1, 0), gsl_matrix_get (a, 1, 1), gsl_matrix_get (a, 1, 2));
printf ("%f %f %f \n",gsl_matrix_get (a, 2, 0), gsl_matrix_get (a, 2, 1), gsl_matrix_get (a, 2, 2));
printf ("%f %f %f \n",gsl_matrix_get (a, 3, 0), gsl_matrix_get (a, 3, 1), gsl_matrix_get (a, 3, 2));
gsl_linalg_SV_decomp (a, v, s, work);
printf ("\n");
printf ("%f %f %f \n",gsl_matrix_get (a, 0, 0), gsl_matrix_get (a, 0, 1), gsl_matrix_get (a, 0, 2));
printf ("%f %f %f \n",gsl_matrix_get (a, 1, 0), gsl_matrix_get (a, 1, 1), gsl_matrix_get (a, 1, 2));
printf ("%f %f %f \n",gsl_matrix_get (a, 2, 0), gsl_matrix_get (a, 2, 1), gsl_matrix_get (a, 2, 2));
printf ("%f %f %f \n",gsl_matrix_get (a, 3, 0), gsl_matrix_get (a, 3, 1), gsl_matrix_get (a, 3, 2));
printf ("\n");
printf ("%f %f %f \n",gsl_vector_get (s, 0), gsl_vector_get (s, 1), gsl_vector_get (s, 2));
printf ("\n");
printf ("%f %f %f \n",gsl_matrix_get (v, 0, 0), gsl_matrix_get (v, 0, 1), gsl_matrix_get (v, 0, 2));
printf ("%f %f %f \n",gsl_matrix_get (v, 1, 0), gsl_matrix_get (v, 1, 1), gsl_matrix_get (v, 1, 2));
printf ("%f %f %f \n",gsl_matrix_get (v, 2, 0), gsl_matrix_get (v, 2, 1), gsl_matrix_get (v, 2, 2));
gsl_matrix_free (a);
gsl_matrix_free (u);
gsl_matrix_free (v);
gsl_vector_free (work);
gsl_vector_free (s);
return 0;
}
James A. Love
X4477
Pager 1-800-286-1188 Pin# 400659