cleaned up svd

gtsam/base/Matrix.cpp

@@ -621,22 +621,10 @@ Matrix inverse_square_root(const Matrix& A) {
 
 /* ************************************************************************* */
 void svd(const Matrix& A, Matrix& U, Vector& S, Matrix& V) {
-	const size_t m = A.rows(), n = A.cols();
-	if (m < n) {
-		V = trans(A);
-		svd(V, S, U); // A'=V*diag(s)*U'
-	} else {
-		U = A; // copy
-		svd(U, S, V); // call in-place version
-	}
-}
-
-/* ************************************************************************* */
-void svd(Matrix& A, Vector& S, Matrix& V) {
 	Eigen::JacobiSVD<Matrix> svd(A, Eigen::ComputeThinU | Eigen::ComputeThinV);
+	U = svd.matrixU();
 	S = svd.singularValues();
-	A = svd.matrixU() * -1.0; // sign issues in tests - still valid, though
-	V = svd.matrixV() * -1.0;
+	V = svd.matrixV();
 }
 
 /* ************************************************************************* */

gtsam/base/Matrix.h

@@ -411,26 +411,16 @@ Matrix cholesky_inverse(const Matrix &A);
  * SVD computes economy SVD A=U*S*V'
  * @param A an m*n matrix
  * @param U output argument: rotation matrix
- * @param S output argument: vector of singular values, sorted by default, pass false as last argument to avoid sorting!!!
+ * @param S output argument: sorted vector of singular values
  * @param V output argument: rotation matrix
- * if m > n then U*S*V' = (m*n)*(n*n)*(n*n) (the m-n columns of U are of no use)
- * if m < n then U*S*V' = (m*m)*(m*m)*(m*n) (the n-m columns of V are of no use)
+ * if m > n then U*S*V' = (m*n)*(n*n)*(n*n)
+ * if m < n then U*S*V' = (m*m)*(m*m)*(m*n)
  * Careful! The dimensions above reflect V', not V, which is n*m if m<n.
  * U is a basis in R^m, V is a basis in R^n
  * You can just pass empty matrices U,V, and vector S, they will be re-allocated.
  */
 void svd(const Matrix& A, Matrix& U, Vector& S, Matrix& V);
 
-/*
- * In place SVD, will throw an exception when m < n
- * @param A an m*n matrix is modified to contain U
- * @param S output argument: vector of singular values, sorted by default, pass false as last argument to avoid sorting!!!
- * @param V output argument: rotation matrix
- * if m > n then U*S*V' = (m*n)*(n*n)*(n*n) (the m-n columns of U are of no use)
- * You can just pass empty matrix V and vector S, they will be re-allocated.
- */
-void svd(Matrix& A, Vector& S, Matrix& V);
-
 /**
  * Direct linear transform algorithm that calls svd
  * to find a vector v that minimizes the algebraic error A*v

gtsam/base/tests/testMatrix.cpp

@@ -979,9 +979,9 @@ TEST( matrix, svd2 )
 	Matrix U, V;
 	Vector s;
 
-	Matrix expectedU = Matrix_(3, 2, 0.,1.,0.,0.,-1.,0.);
+	Matrix expectedU = Matrix_(3, 2, 0.,-1.,0.,0.,1.,0.);
 	Vector expected_s = Vector_(2, 3.,2.);
-	Matrix expectedV = Matrix_(2, 2, -1.,0.,0.,-1.);
+	Matrix expectedV = Matrix_(2, 2, 1.,0.,0.,1.);
 
 	svd(sampleA, U, s, V);
 
@@ -1023,17 +1023,15 @@ TEST( matrix, svd4 )
 	    0.1270,    0.0975);
 
 	Matrix expectedU = Matrix_(3,2,
-								 /// right column - originally had opposite sign
-	   -0.7397,   -0.6724,
-	   -0.6659,    0.7370,
-	   -0.0970,    0.0689);
+	   0.7397,   0.6724,
+	   0.6659,   -0.7370,
+	   0.0970,   -0.0689);
 
 	Vector expected_s = Vector_(2, 1.6455, 0.1910);
 
 	Matrix expectedV = Matrix_(2,2,
-								/// right column - originally had opposite sign
-	   -0.7403,    0.6723,
-	   -0.6723,   -0.7403);
+	   0.7403,   -0.6723,
+	   0.6723,   0.7403);
 
 	svd(A, U, s, V);
 	Matrix reconstructed = U * diag(s) * trans(V);