This is the mail archive of the
gsl-discuss@sources.redhat.com
mailing list for the GSL project.
Re: Problem with Singular Value Decomposition Algorithm
- To: <bjg at network-theory dot co dot uk>, <gsl-discuss at sources dot redhat dot com>
- Subject: Re: Problem with Singular Value Decomposition Algorithm
- From: "Jim Love" <Jim dot Love at asml dot com>
- Date: Thu, 13 Sep 2001 11:42:24 -0400
This code fixes the order problem of the S vector and the other matrix, but their is still a sign problem. Using this matrix for A:
1.000000 1.000000 0.975000
1.000000 -1.000000 0.975000
-1.000000 -1.000000 -0.925000
-1.000000 1.000000 -1.025000
The modified code produces:
s:
2.793961 2.000000 0.035791
This is correct!
For V:
-0.715538 -0.025633 0.698103
0.018347 -0.999671 -0.017900
-0.698332 -0.000000 -0.715774
This is NOT correct!
The correct answer for V is:
-0.7155 0.0256 -0.6981
0.0183 0.9997 0.0179
-0.6983 -0.0000 0.7158
U is also wrong: the program outputs:
-0.493230 -0.512652 -0.493875
-0.506363 0.487019 0.506368
0.480733 0.512652 -0.506047
0.518861 -0.487019 0.493554
The correct U is:
-0.4932 0.5127 0.4939
-0.5064 -0.4870 -0.5064
0.4807 -0.5127 0.5060
0.5189 0.4870 -0.4936
Note last column missing for both solutions for U.
James A. Love
X4477
Pager 1-800-286-1188 Pin# 400659
>>> Brian Gough <bjg@network-theory.co.uk> 09/12/01 06:35PM >>>
Patch below fixes the problem, I'll do more testing on that routine
tomorrow.
Index: svdstep.c
===================================================================
RCS file: /cvs/gsl/gsl/linalg/svdstep.c,v
retrieving revision 1.3
diff -u -r1.3 svdstep.c
--- svdstep.c 2001/08/29 15:37:21 1.3
+++ svdstep.c 2001/09/12 22:33:01
@@ -126,14 +126,14 @@
if (d0 == 0.0)
{
- /* Eliminate off-diagonal element in [0,f1;0,d1] */
+ /* Eliminate off-diagonal element in [0,f1;0,d1] to make [d,0;0,0] */
create_givens (f0, d1, &c, &s);
- /* compute B <= G^T B */
+ /* compute B <= G^T B X, where X = [0,1;1,0] */
- gsl_vector_set (f, 0, c * f0 - s * d1);
- gsl_vector_set (d, 1, s * f0 + c * d1);
+ gsl_vector_set (d, 0, c * f0 - s * d1);
+ gsl_vector_set (f, 1, s * f0 + c * d1);
/* Compute U <= U G */
@@ -145,6 +145,10 @@
gsl_matrix_set (U, i, 1, s * Uip + c * Uiq);
}
+ /* Compute V <= V X */
+
+ gsl_matrix_swap_columns (V, 0, 1);
+
return;
}
else if (d1 == 0.0)
@@ -194,17 +198,24 @@
gsl_matrix_set (V, i, 1, s * Vip + c * Viq);
}
- /* Eliminate off-diagonal elements, work with the column that
- has the largest norm */
+ /* Eliminate off-diagonal elements, bring column with largest
+ norm to first column */
- if (fabs(a11) + fabs(a21) > fabs(a12) + fabs(a22))
+ if (hypot(a11, a21) < hypot(a12,a22))
{
- create_givens (a11, a21, &c, &s);
+ double t1, t2;
+
+ /* B <= B X */
+
+ t1 = a11; a11 = a12; a12 = t1;
+ t2 = a21; a21 = a22; a22 = t2;
+
+ /* V <= V X */
+
+ gsl_matrix_swap_columns(V, 0, 1);
}
- else
- {
- create_givens (a22, -a12, &c, &s);
- }
+
+ create_givens (a11, a21, &c, &s);
/* compute B <= G^T B */