diff --git a/gtsam/base/Matrix.cpp b/gtsam/base/Matrix.cpp index d576771e1..e0cb055b7 100644 --- a/gtsam/base/Matrix.cpp +++ b/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 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(); } /* ************************************************************************* */ diff --git a/gtsam/base/Matrix.h b/gtsam/base/Matrix.h index d01d6598e..1b4b58414 100644 --- a/gtsam/base/Matrix.h +++ b/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 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 diff --git a/gtsam/base/tests/testMatrix.cpp b/gtsam/base/tests/testMatrix.cpp index 0e8dff261..be46b3ddf 100644 --- a/gtsam/base/tests/testMatrix.cpp +++ b/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);