Re: Problem with Singular Value Decomposition Algorithm

Re: Problem with Singular Value Decomposition Algorithm

[Jim Love]

In every text that I have read dealing with the process of SVD, they have always placed the eigenvalues in a decreasing order.  The stated reason is to provide a standard solution and methodology.  Why would anybody build an API that does something different??  If you are using the SVD to find the least squares fit of a plane or a cylinder for example, order matters (if you want the right answer hehehe).  Yes I can add to my code to re-order the output, but so will most everybody else that uses this function.  So it should be done in the API.  Just a suggestion.

James A. Love
X4477
Pager 1-800-286-1188 Pin# 400659

>>> James Theiler <jt@lanl.gov> 09/12/01 10:51AM >>>
On Wed, 12 Sep 2001, Jim Love wrote:

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

I don't believe this is a bug at all.
Both answers are correct.  For SVD, the algorithm is
asked to find matrices U,S,Q such that U.S.Q' equals
the original matrix.  If you swap the columns of Q
and swap the equivalent rows of U, and also swap the
corresponding elements of S, you have another solution.

Also if you multiply one of the columns of Q by -1,
and multiply the corresponding row of U by -1, you
will get another equivalent solution.

Sometimes you want the solution with the eigenvalues (the
diagonal elements of the S matrix -- represented as a vector)
sorted numerically, as the one you cite as "correct" but
that is a convenience that can be performed after the fact.
And in fact, I notice that the documentation says that
eigenvalues are "generally chosen" to form a non-decreasing
sequence.  Maybe worth a sentence to say that these algorithms
do not always do that, or else add a utility that does this
postprocessing step.

jt

---------------------------------------------
James Theiler                     jt@lanl.gov 
MS-D436, NIS-2, LANL        tel: 505/665-5682
Los Alamos, NM 87545        fax: 505/665-4414
----- Space and Remote Sensing Sciences -----

Re: Problem with Singular Value Decomposition Algorithm

[Brian Gough]

Jim Love writes:
 > In every text that I have read dealing with the process of SVD, they have always placed the eigenvalues in a decreasing order.  The stated reason is to provide a standard solution and methodology.  Why would anybody build an API that does something different??  If you are using the SVD to find the least squares fit of a plane or a cylinder for example, order matters (if you want the right answer hehehe).  Yes I can add to my code to re-order the output, but so will most everybody else that uses this function.  So it should be done in the API.  Just a suggestion.
 > 

Thanks for the bug report.  

I'm pretty sure this is caused by me adding some new code to
explicitly handle the 2x2 case in the past few weeks.  I overlooked
the ordering of the singular values for that routine and I can see
that it will be unsorted. That would explain why the first two
elements are wrong --- they are handled by the 2x2 routine. I'll add
the appropriate sorting to that case.

regards
Brian Gough

Re: Problem with Singular Value Decomposition Algorithm

[Brian Gough]

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

Re: Problem with Singular Value Decomposition Algorithm

[Jim Love]

I guess the best way to explain this is via a physical example.  If I have a plane that is tilted some angle about the Y-axis and I make 4 measurements of the distance from the xy plane to the surface at 4 points.  

1,1 I get 1
1,-1 I get 1
-1, -1 I get -0.9
-1, 1 I get -1

You can then process the data to set up the least square fit to find the "best fit" plane by making a matrix:

x1-mean(x) y1-mean(y) z1-mean(z)
x2-mean(x) y2-mean(y) z2-mean(z)
x3-mean(x) y3-mean(y) z3-mean(z)
x4-mean(x) y4-mean(y) z4-mean(z)

so in this example you would have

1.000000 1.000000 0.975000
1.000000 -1.000000 0.975000
-1.000000 -1.000000 -0.925000
-1.000000 1.000000 -1.025000

find the svd of this matrix.

The solution of the problem is the vector from v the corresponds to the smallest eigenvalue

in our case:

min eigenvalue =  0.035791

the vector from v is -0.698332 -0.000000 0.715774

the equation of a plane is 

a(x-xo) + b(y-yo) + c(z-zo) = 0

in this case:

a = -0.698332
b = -0.000000
c = 0.715774

therefore

-0.698332*x + -0.000000*y + 0.715774*z = 0

substituting back our initial points we get:

1,1 = 1.000633
1,-1 = 1.000633
-1, -1 = -0.950633
1,-1 = -0.950633

You can see by inspection that this is correct.

Now the solution that you API provides is:

a = -0.698332
b = -0.000000
c = -0.715774

-0.698332*x + -0.000000*y + -0.715774*z = 0

1,1 = -0.950633
1,-1 = -0.950633
-1,-1 = 1.000633
1,-1 = 1.000633

This is clearly not the "best fit" plane to the data and is clearly wrong.  

The sign is not arbitrary for this case or any use of SVD.


James A. Love
X4477
Pager 1-800-286-1188 Pin# 400659

>>> Brian Gough <bjg@network-theory.co.uk> 09/13/01 12:21PM >>>
Columns 2 and 3 have the opposite sign, but this is the arbitrary
minus-sign factor referred to earlier.  The results satisfy m =
u*diag(s)*v' and u'*u = I, v'*v = I numerically --- so the
decomposition looks ok.  Let me know if there's something I missing
here.

regards
Brian Gough

Jim Love writes:
 > 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.
 > 

Re: Problem with Singular Value Decomposition Algorithm

[James Theiler]

i think there may be some confusion about transposes; see below.

jt

On Thu, 13 Sep 2001, Jim Love wrote:

] I guess the best way to explain this is via a physical example.
] If I have a plane that is tilted some angle about the Y-axis and I
] make 4 measurements of the distance from the xy plane to the
] surface at 4 points.
]
] 1,1 I get 1
] 1,-1 I get 1
] -1, -1 I get -0.9
] -1, 1 I get -1
]
] You can then process the data to set up the least square fit to
] find the "best fit" plane by making a matrix:
]
] x1-mean(x) y1-mean(y) z1-mean(z)
] x2-mean(x) y2-mean(y) z2-mean(z)
] x3-mean(x) y3-mean(y) z3-mean(z)
] x4-mean(x) y4-mean(y) z4-mean(z)
]
] so in this example you would have
]
] 1.000000 1.000000 0.975000
] 1.000000 -1.000000 0.975000
] -1.000000 -1.000000 -0.925000
] -1.000000 1.000000 -1.025000
]
] find the svd of this matrix.
]
] The solution of the problem is the vector from v the corresponds
] to the smallest eigenvalue
]
] in our case:
]
] min eigenvalue =  0.035791
]
] the vector from v is -0.698332 -0.000000 0.715774

isn't it 0.698103, -0.017900, -0.715774
ie, the third column and not the third row?

note also, if you use the columns of V
and not the rows, then the sign ambiguity
doesn't matter.

]
] the equation of a plane is
]
] a(x-xo) + b(y-yo) + c(z-zo) = 0
]
] in this case:
]
] a = -0.698332
] b = -0.000000
] c = 0.715774
]
] therefore
]
] -0.698332*x + -0.000000*y + 0.715774*z = 0
]
] substituting back our initial points we get:
]
] 1,1 = 1.000633
] 1,-1 = 1.000633
] -1, -1 = -0.950633
] 1,-1 = -0.950633   <== you mean -1,1 ??

let's see z = zo - [a(x-xo)+b(y-yo)]/c
if i use

a = 0.698103
b = -0.017900
c = -0.715774

then i get

 1  1    0.97530415
 1 -1    1.0253199
-1 -1   -0.92530415
-1  1   -0.97531993

which also seems reasonable to me.



]
] You can see by inspection that this is correct.
]
] Now the solution that you API provides is:
]
] a = -0.698332
] b = -0.000000
] c = -0.715774
]
] -0.698332*x + -0.000000*y + -0.715774*z = 0
]
] 1,1 = -0.950633
] 1,-1 = -0.950633
] -1,-1 = 1.000633
] 1,-1 = 1.000633
]
] This is clearly not the "best fit" plane to the data and is clearly wrong.
]
] The sign is not arbitrary for this case or any use of SVD.
]
]
] James A. Love
] X4477
] Pager 1-800-286-1188 Pin# 400659
]
] >>> Brian Gough <bjg@network-theory.co.uk> 09/13/01 12:21PM >>>
] Columns 2 and 3 have the opposite sign, but this is the arbitrary
] minus-sign factor referred to earlier.  The results satisfy m =
] u*diag(s)*v' and u'*u = I, v'*v = I numerically --- so the
] decomposition looks ok.  Let me know if there's something I missing
] here.
]
] regards
] Brian Gough
]
] Jim Love writes:
]  > 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 Theiler                     jt@lanl.gov
MS-D436, NIS-2, LANL        tel: 505/665-5682
Los Alamos, NM 87545        fax: 505/665-4414
----- Space and Remote Sensing Sciences -----

Re: Problem with Singular Value Decomposition Algorithm

[Brian Gough]

James Theiler writes:
 > i think there may be some confusion about transposes; see below.

 > isn't it 0.698103, -0.017900, -0.715774
 > ie, the third column and not the third row?

Yes, the right singular vectors are given by the columns of V. The
third column is the appropriate one.  

In this case V happens to be almost symmetric so the difference
between the two is less obvious than it would be for a random matrix.

Re: Problem with Singular Value Decomposition Algorithm

[Peter Schmitteckert]

Dear All,

please excuse me, I have used gsl only once and havn't read the documents,
but I got curious by the follwing comment:

On Thursday 13 September 2001 22:33, Brian Gough wrote:

> Yes, the right singular vectors are given by the columns of V. The
> third column is the appropriate one.
>

Isn't that counter-intuitive fo a C programmer ?
In the original post I found no hints on a desired
Fortran semantic for the memory layout.  But I would 
expect that the default semantics of a C matrix library is
C style, so eigen vectors / singular vectors should be in
the rows, otherwise they are not in a contiguous memory
block.

Best wishes
Peter

P.S: 
I just looked up the documentation:

> Matrices are stored in row-major order, meaning that each row of elements 
> forms a contiguous block in memory. This is the standard "C-language 
> ordering" of two-dimensional arrays. Note that FORTRAN stores arrays in 
> column-major order. The number of rows is size1. The range of valid 

therefore, returnning vectors in the coulmn should be considered
as a bug. Have I missed something ?
Note, that the documentation says: "The matrix Q contains the elements of Q 
in untransposed form".

Re: Problem with Singular Value Decomposition Algorithm

[Brian Gough]

Peter Schmitteckert writes:
 > > Matrices are stored in row-major order, meaning that each row of elements 
 > > forms a contiguous block in memory. This is the standard "C-language 
 > > ordering" of two-dimensional arrays. Note that FORTRAN stores arrays in 
 > > column-major order. The number of rows is size1. The range of valid 
 > 
 > therefore, returnning vectors in the coulmn should be considered
 > as a bug. Have I missed something ?
 > Note, that the documentation says: "The matrix Q contains the elements of Q 
 > in untransposed form".

The routine does what the documentation says... so this is not a bug ;-)

GSL vectors use strides, so both rows and columns can be accessed as
vectors without problems.

As for efficiency -- that will depend on the application's access
pattern, which could be either row or column oriented, so there is
always a winner and a loser.

regards
Brian Gough

Re: Problem with Singular Value Decomposition Algorithm

[Jim Love]

You have made a mistake,  the equation for a z at point x,y of the "best fit" plane is:

z = (a(x-mean(x))+b(y-mean(y))-c(mean(z)))/-c

so using your values of a,b,c i get

1,1 = -0.975
1,-1 = -1.025

no need to go and calc the others - this is clearly wrong.

James A. Love
X4477
Pager 1-800-286-1188 Pin# 400659

>>> James Theiler <jt@lanl.gov> 09/13/01 03:28PM >>>
i think there may be some confusion about transposes; see below.

jt

On Thu, 13 Sep 2001, Jim Love wrote:

] I guess the best way to explain this is via a physical example.
] If I have a plane that is tilted some angle about the Y-axis and I
] make 4 measurements of the distance from the xy plane to the
] surface at 4 points.
]
] 1,1 I get 1
] 1,-1 I get 1
] -1, -1 I get -0.9
] -1, 1 I get -1
]
] You can then process the data to set up the least square fit to
] find the "best fit" plane by making a matrix:
]
] x1-mean(x) y1-mean(y) z1-mean(z)
] x2-mean(x) y2-mean(y) z2-mean(z)
] x3-mean(x) y3-mean(y) z3-mean(z)
] x4-mean(x) y4-mean(y) z4-mean(z)
]
] so in this example you would have
]
] 1.000000 1.000000 0.975000
] 1.000000 -1.000000 0.975000
] -1.000000 -1.000000 -0.925000
] -1.000000 1.000000 -1.025000
]
] find the svd of this matrix.
]
] The solution of the problem is the vector from v the corresponds
] to the smallest eigenvalue
]
] in our case:
]
] min eigenvalue =  0.035791
]
] the vector from v is -0.698332 -0.000000 0.715774

isn't it 0.698103, -0.017900, -0.715774
ie, the third column and not the third row?

note also, if you use the columns of V
and not the rows, then the sign ambiguity
doesn't matter.

]
] the equation of a plane is
]
] a(x-xo) + b(y-yo) + c(z-zo) = 0
]
] in this case:
]
] a = -0.698332
] b = -0.000000
] c = 0.715774
]
] therefore
]
] -0.698332*x + -0.000000*y + 0.715774*z = 0
]
] substituting back our initial points we get:
]
] 1,1 = 1.000633
] 1,-1 = 1.000633
] -1, -1 = -0.950633
] 1,-1 = -0.950633   <== you mean -1,1 ??

let's see z = zo - [a(x-xo)+b(y-yo)]/c
if i use

a = 0.698103
b = -0.017900
c = -0.715774

then i get

 1  1    0.97530415
 1 -1    1.0253199
-1 -1   -0.92530415
-1  1   -0.97531993

which also seems reasonable to me.



]
] You can see by inspection that this is correct.
]
] Now the solution that you API provides is:
]
] a = -0.698332
] b = -0.000000
] c = -0.715774
]
] -0.698332*x + -0.000000*y + -0.715774*z = 0
]
] 1,1 = -0.950633
] 1,-1 = -0.950633
] -1,-1 = 1.000633
] 1,-1 = 1.000633
]
] This is clearly not the "best fit" plane to the data and is clearly wrong.
]
] The sign is not arbitrary for this case or any use of SVD.
]
]
] James A. Love
] X4477
] Pager 1-800-286-1188 Pin# 400659
]
] >>> Brian Gough <bjg@network-theory.co.uk> 09/13/01 12:21PM >>>
] Columns 2 and 3 have the opposite sign, but this is the arbitrary
] minus-sign factor referred to earlier.  The results satisfy m =
] u*diag(s)*v' and u'*u = I, v'*v = I numerically --- so the
] decomposition looks ok.  Let me know if there's something I missing
] here.
]
] regards
] Brian Gough
]
] Jim Love writes:
]  > 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 Theiler                     jt@lanl.gov 
MS-D436, NIS-2, LANL        tel: 505/665-5682
Los Alamos, NM 87545        fax: 505/665-4414
----- Space and Remote Sensing Sciences -----

Problem with Singular Value Decomposition Algorithm

[Jim Love]

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

Re: Problem with Singular Value Decomposition Algorithm

[James Theiler]

On Wed, 12 Sep 2001, Jim Love wrote:

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

I don't believe this is a bug at all.
Both answers are correct.  For SVD, the algorithm is
asked to find matrices U,S,Q such that U.S.Q' equals
the original matrix.  If you swap the columns of Q
and swap the equivalent rows of U, and also swap the
corresponding elements of S, you have another solution.

Also if you multiply one of the columns of Q by -1,
and multiply the corresponding row of U by -1, you
will get another equivalent solution.

Sometimes you want the solution with the eigenvalues (the
diagonal elements of the S matrix -- represented as a vector)
sorted numerically, as the one you cite as "correct" but
that is a convenience that can be performed after the fact.
And in fact, I notice that the documentation says that
eigenvalues are "generally chosen" to form a non-decreasing
sequence.  Maybe worth a sentence to say that these algorithms
do not always do that, or else add a utility that does this
postprocessing step.

jt

---------------------------------------------
James Theiler                     jt@lanl.gov
MS-D436, NIS-2, LANL        tel: 505/665-5682
Los Alamos, NM 87545        fax: 505/665-4414
----- Space and Remote Sensing Sciences -----

Re: Problem with Singular Value Decomposition Algorithm

[Jim Love]

I still think their is a problem with your API.  I have done this problem using a pencil and paper, using Matlab, using Mathematica, and using the LAPACK DGESVD routine.  These 4 implementations all produce the same exact solution for the U S and Q matrix - only yours is different.  How do you explain this?

James A. Love
X4477
Pager 1-800-286-1188 Pin# 400659

Re: Problem with Singular Value Decomposition Algorithm

[Brian Gough]

Jim Love writes:
 > I still think their is a problem with your API.  I have done this problem using a pencil and paper, using Matlab, using Mathematica, and using the LAPACK DGESVD routine.  These 4 implementations all produce the same exact solution for the U S and Q matrix - only yours is different.  How do you explain this?
 > 

Matlab uses LAPACK.  I believe Mathematica might use that too now.  So
that's probably why they are giving the same results.

Provided the results returned by GSL satisfy u'*u = I, v'*v = I and
u*diag(s)*v' = A then it's a valid singular value decomposition.  From
my own testing this checks out ok, and also with the numbers you sent
in your emails.  If you have an example that doesn't then that would
be a bug.

Re: Problem with Singular Value Decomposition Algorithm

[Brian Gough]

Columns 2 and 3 have the opposite sign, but this is the arbitrary
minus-sign factor referred to earlier.  The results satisfy m =
u*diag(s)*v' and u'*u = I, v'*v = I numerically --- so the
decomposition looks ok.  Let me know if there's something I missing
here.

regards
Brian Gough

Jim Love writes:
 > 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.
 > 

Re: Problem with Singular Value Decomposition Algorithm

[Jim Love]

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