diff options
| -rw-r--r-- | Eigen/SVD | 1 | ||||
| -rw-r--r-- | Eigen/src/Core/MatrixBase.h | 2 | ||||
| -rw-r--r-- | Eigen/src/Core/util/ForwardDeclarations.h | 1 | ||||
| -rw-r--r-- | Eigen/src/Geometry/Transform.h | 26 | ||||
| -rw-r--r-- | Eigen/src/Geometry/Umeyama.h | 2 | ||||
| -rw-r--r-- | Eigen/src/SVD/SVD.h | 604 | ||||
| -rw-r--r-- | test/CMakeLists.txt | 1 | ||||
| -rw-r--r-- | test/nomalloc.cpp | 3 | ||||
| -rw-r--r-- | test/svd.cpp | 123 |
9 files changed, 26 insertions, 737 deletions
@@ -23,7 +23,6 @@ namespace Eigen { */ #include "src/misc/Solve.h" -#include "src/SVD/SVD.h" #include "src/SVD/JacobiSVD.h" #include "src/SVD/UpperBidiagonalization.h" diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index 2765db5b3..f324a0870 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -331,8 +331,6 @@ template<typename Derived> class MatrixBase /////////// SVD module /////////// - SVD<PlainObject> svd() const; - /////////// Geometry module /////////// template<typename OtherDerived> diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h index 1729432ca..f6cf3fe41 100644 --- a/Eigen/src/Core/util/ForwardDeclarations.h +++ b/Eigen/src/Core/util/ForwardDeclarations.h @@ -169,7 +169,6 @@ template<typename MatrixType> struct ei_inverse_impl; template<typename MatrixType> class HouseholderQR; template<typename MatrixType> class ColPivHouseholderQR; template<typename MatrixType> class FullPivHouseholderQR; -template<typename MatrixType> class SVD; template<typename MatrixType, int QRPreconditioner = ColPivHouseholderQRPreconditioner> class JacobiSVD; template<typename MatrixType, int UpLo = Lower> class LLT; template<typename MatrixType, int UpLo = Lower> class LDLT; diff --git a/Eigen/src/Geometry/Transform.h b/Eigen/src/Geometry/Transform.h index 58bb2d0c0..e8099495d 100644 --- a/Eigen/src/Geometry/Transform.h +++ b/Eigen/src/Geometry/Transform.h @@ -928,7 +928,18 @@ template<typename Scalar, int Dim, int Mode> template<typename RotationMatrixType, typename ScalingMatrixType> void Transform<Scalar,Dim,Mode>::computeRotationScaling(RotationMatrixType *rotation, ScalingMatrixType *scaling) const { - linear().svd().computeRotationScaling(rotation, scaling); + JacobiSVD<LinearMatrixType> svd(linear(), ComputeFullU | ComputeFullV); + + Scalar x = (svd.matrixU() * svd.matrixV().adjoint()).determinant(); // so x has absolute value 1 + VectorType sv(svd.singularValues()); + sv.coeffRef(0) *= x; + if(scaling) scaling->lazyAssign(svd.matrixV() * sv.asDiagonal() * svd.matrixV().adjoint()); + if(rotation) + { + LinearMatrixType m(svd.matrixU()); + m.col(0) /= x; + rotation->lazyAssign(m * svd.matrixV().adjoint()); + } } /** decomposes the linear part of the transformation as a product rotation x scaling, the scaling being @@ -946,7 +957,18 @@ template<typename Scalar, int Dim, int Mode> template<typename ScalingMatrixType, typename RotationMatrixType> void Transform<Scalar,Dim,Mode>::computeScalingRotation(ScalingMatrixType *scaling, RotationMatrixType *rotation) const { - linear().svd().computeScalingRotation(scaling, rotation); + JacobiSVD<LinearMatrixType> svd(linear(), ComputeFullU | ComputeFullV); + + Scalar x = (svd.matrixU() * svd.matrixV().adjoint()).determinant(); // so x has absolute value 1 + VectorType sv(svd.singularValues()); + sv.coeffRef(0) *= x; + if(scaling) scaling->lazyAssign(svd.matrixU() * sv.asDiagonal() * svd.matrixU().adjoint()); + if(rotation) + { + LinearMatrixType m(svd.matrixU()); + m.col(0) /= x; + rotation->lazyAssign(m * svd.matrixV().adjoint()); + } } /** Convenient method to set \c *this from a position, orientation and scale diff --git a/Eigen/src/Geometry/Umeyama.h b/Eigen/src/Geometry/Umeyama.h index fa6bfcbee..04449f623 100644 --- a/Eigen/src/Geometry/Umeyama.h +++ b/Eigen/src/Geometry/Umeyama.h @@ -141,7 +141,7 @@ umeyama(const MatrixBase<Derived>& src, const MatrixBase<OtherDerived>& dst, boo // Eq. (38) const MatrixType sigma = one_over_n * dst_demean * src_demean.transpose(); - SVD<MatrixType> svd(sigma); + JacobiSVD<MatrixType> svd(sigma, ComputeFullU | ComputeFullV); // Initialize the resulting transformation with an identity matrix... TransformationMatrixType Rt = TransformationMatrixType::Identity(m+1,m+1); diff --git a/Eigen/src/SVD/SVD.h b/Eigen/src/SVD/SVD.h deleted file mode 100644 index e8c970d4f..000000000 --- a/Eigen/src/SVD/SVD.h +++ /dev/null @@ -1,604 +0,0 @@ -// This file is part of Eigen, a lightweight C++ template library -// for linear algebra. -// -// Copyright (C) 2008-2009 Gael Guennebaud <gael.guennebaud@inria.fr> -// -// Eigen is free software; you can redistribute it and/or -// modify it under the terms of the GNU Lesser General Public -// License as published by the Free Software Foundation; either -// version 3 of the License, or (at your option) any later version. -// -// Alternatively, you can redistribute it and/or -// modify it under the terms of the GNU General Public License as -// published by the Free Software Foundation; either version 2 of -// the License, or (at your option) any later version. -// -// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY -// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the -// GNU General Public License for more details. -// -// You should have received a copy of the GNU Lesser General Public -// License and a copy of the GNU General Public License along with -// Eigen. If not, see <http://www.gnu.org/licenses/>. - -#ifndef EIGEN_SVD_H -#define EIGEN_SVD_H - -template<typename MatrixType, typename Rhs> struct ei_svd_solve_impl; - -/** \ingroup SVD_Module - * - * - * \class SVD - * - * \brief Standard SVD decomposition of a matrix and associated features - * - * \param MatrixType the type of the matrix of which we are computing the SVD decomposition - * - * This class performs a standard SVD decomposition of a real matrix A of size \c M x \c N. - * - * Requires M >= N, in other words, at least as many rows as columns. - * - * \sa MatrixBase::SVD() - */ -template<typename _MatrixType> class SVD -{ - public: - typedef _MatrixType MatrixType; - typedef typename MatrixType::Scalar Scalar; - typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar; - typedef typename MatrixType::Index Index; - - enum { - RowsAtCompileTime = MatrixType::RowsAtCompileTime, - ColsAtCompileTime = MatrixType::ColsAtCompileTime, - PacketSize = ei_packet_traits<Scalar>::size, - AlignmentMask = int(PacketSize)-1, - MinSize = EIGEN_SIZE_MIN_PREFER_DYNAMIC(RowsAtCompileTime, ColsAtCompileTime), - MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, - MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime, - MatrixOptions = MatrixType::Options - }; - - typedef typename ei_plain_col_type<MatrixType>::type ColVector; - typedef typename ei_plain_row_type<MatrixType>::type RowVector; - - typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, MatrixOptions, MaxRowsAtCompileTime, MaxRowsAtCompileTime> MatrixUType; - typedef Matrix<Scalar, ColsAtCompileTime, ColsAtCompileTime, MatrixOptions, MaxColsAtCompileTime, MaxColsAtCompileTime> MatrixVType; - typedef ColVector SingularValuesType; - - /** - * \brief Default Constructor. - * - * The default constructor is useful in cases in which the user intends to - * perform decompositions via SVD::compute(const MatrixType&). - */ - SVD() : m_matU(), m_matV(), m_sigma(), m_isInitialized(false) {} - - /** \brief Default Constructor with memory preallocation - * - * Like the default constructor but with preallocation of the internal data - * according to the specified problem \a size. - * \sa JacobiSVD() - */ - SVD(Index rows, Index cols) : m_matU(rows, rows), - m_matV(cols,cols), - m_sigma(std::min(rows, cols)), - m_workMatrix(rows, cols), - m_rv1(cols), - m_isInitialized(false) - { - ei_assert(rows >= cols && "SVD is only defined if rows>=cols."); - } - - SVD(const MatrixType& matrix) : m_matU(matrix.rows(), matrix.rows()), - m_matV(matrix.cols(),matrix.cols()), - m_sigma(std::min(matrix.rows(), matrix.cols())), - m_workMatrix(matrix.rows(), matrix.cols()), - m_rv1(matrix.cols()), - m_isInitialized(false) - { - compute(matrix); - } - - /** \returns a solution of \f$ A x = b \f$ using the current SVD decomposition of A. - * - * \param b the right-hand-side of the equation to solve. - * - * \note_about_checking_solutions - * - * \note_about_arbitrary_choice_of_solution - * - * \sa MatrixBase::svd(), - */ - template<typename Rhs> - inline const ei_solve_retval<SVD, Rhs> - solve(const MatrixBase<Rhs>& b) const - { - ei_assert(m_isInitialized && "SVD is not initialized."); - return ei_solve_retval<SVD, Rhs>(*this, b.derived()); - } - - const MatrixUType& matrixU() const - { - ei_assert(m_isInitialized && "SVD is not initialized."); - return m_matU; - } - - const SingularValuesType& singularValues() const - { - ei_assert(m_isInitialized && "SVD is not initialized."); - return m_sigma; - } - - const MatrixVType& matrixV() const - { - ei_assert(m_isInitialized && "SVD is not initialized."); - return m_matV; - } - - SVD& compute(const MatrixType& matrix); - - template<typename UnitaryType, typename PositiveType> - void computeUnitaryPositive(UnitaryType *unitary, PositiveType *positive) const; - template<typename PositiveType, typename UnitaryType> - void computePositiveUnitary(PositiveType *positive, UnitaryType *unitary) const; - template<typename RotationType, typename ScalingType> - void computeRotationScaling(RotationType *unitary, ScalingType *positive) const; - template<typename ScalingType, typename RotationType> - void computeScalingRotation(ScalingType *positive, RotationType *unitary) const; - - inline Index rows() const - { - ei_assert(m_isInitialized && "SVD is not initialized."); - return m_rows; - } - - inline Index cols() const - { - ei_assert(m_isInitialized && "SVD is not initialized."); - return m_cols; - } - - protected: - // Computes (a^2 + b^2)^(1/2) without destructive underflow or overflow. - inline static Scalar pythag(Scalar a, Scalar b) - { - Scalar abs_a = ei_abs(a); - Scalar abs_b = ei_abs(b); - if (abs_a > abs_b) - return abs_a*ei_sqrt(Scalar(1.0)+ei_abs2(abs_b/abs_a)); - else - return (abs_b == Scalar(0.0) ? Scalar(0.0) : abs_b*ei_sqrt(Scalar(1.0)+ei_abs2(abs_a/abs_b))); - } - - inline static Scalar sign(Scalar a, Scalar b) - { - return (b >= Scalar(0.0) ? ei_abs(a) : -ei_abs(a)); - } - - protected: - /** \internal */ - MatrixUType m_matU; - /** \internal */ - MatrixVType m_matV; - /** \internal */ - SingularValuesType m_sigma; - MatrixType m_workMatrix; - RowVector m_rv1; - bool m_isInitialized; - Index m_rows, m_cols; -}; - -/** Computes / recomputes the SVD decomposition A = U S V^* of \a matrix - * - * \note this code has been adapted from Numerical Recipes, third edition. - * - * \returns a reference to *this - */ -template<typename MatrixType> -SVD<MatrixType>& SVD<MatrixType>::compute(const MatrixType& matrix) -{ - const Index m = m_rows = matrix.rows(); - const Index n = m_cols = matrix.cols(); - - m_matU.resize(m, m); - m_matU.setZero(); - m_sigma.resize(n); - m_matV.resize(n,n); - m_workMatrix = matrix; - - Index max_iters = 30; - - MatrixVType& V = m_matV; - MatrixType& A = m_workMatrix; - SingularValuesType& W = m_sigma; - - bool flag; - Index i=0,its=0,j=0,k=0,l=0,nm=0; - Scalar anorm, c, f, g, h, s, scale, x, y, z; - bool convergence = true; - Scalar eps = NumTraits<Scalar>::dummy_precision(); - - m_rv1.resize(n); - - g = scale = anorm = 0; - // Householder reduction to bidiagonal form. - for (i=0; i<n; i++) - { - l = i+2; - m_rv1[i] = scale*g; - g = s = scale = 0.0; - if (i < m) - { - scale = A.col(i).tail(m-i).cwiseAbs().sum(); - if (scale != Scalar(0)) - { - for (k=i; k<m; k++) - { - A(k, i) /= scale; - s += A(k, i)*A(k, i); - } - f = A(i, i); - g = -sign( ei_sqrt(s), f ); - h = f*g - s; - A(i, i)=f-g; - for (j=l-1; j<n; j++) - { - s = A.col(j).tail(m-i).dot(A.col(i).tail(m-i)); - f = s/h; - A.col(j).tail(m-i) += f*A.col(i).tail(m-i); - } - A.col(i).tail(m-i) *= scale; - } - } - W[i] = scale * g; - g = s = scale = 0.0; - if (i+1 <= m && i+1 != n) - { - scale = A.row(i).tail(n-l+1).cwiseAbs().sum(); - if (scale != Scalar(0)) - { - for (k=l-1; k<n; k++) - { - A(i, k) /= scale; - s += A(i, k)*A(i, k); - } - f = A(i,l-1); - g = -sign(ei_sqrt(s),f); - h = f*g - s; - A(i,l-1) = f-g; - m_rv1.tail(n-l+1) = A.row(i).tail(n-l+1)/h; - for (j=l-1; j<m; j++) - { - s = A.row(i).tail(n-l+1).dot(A.row(j).tail(n-l+1)); - A.row(j).tail(n-l+1) += s*m_rv1.tail(n-l+1).transpose(); - } - A.row(i).tail(n-l+1) *= scale; - } - } - anorm = std::max( anorm, (ei_abs(W[i])+ei_abs(m_rv1[i])) ); - } - // Accumulation of right-hand transformations. - for (i=n-1; i>=0; i--) - { - //Accumulation of right-hand transformations. - if (i < n-1) - { - if (g != Scalar(0.0)) - { - for (j=l; j<n; j++) //Double division to avoid possible underflow. - V(j, i) = (A(i, j)/A(i, l))/g; - for (j=l; j<n; j++) - { - s = V.col(j).tail(n-l).dot(A.row(i).tail(n-l)); - V.col(j).tail(n-l) += s * V.col(i).tail(n-l); - } - } - V.row(i).tail(n-l).setZero(); - V.col(i).tail(n-l).setZero(); - } - V(i, i) = 1.0; - g = m_rv1[i]; - l = i; - } - // Accumulation of left-hand transformations. - for (i=std::min(m,n)-1; i>=0; i--) - { - l = i+1; - g = W[i]; - if (n-l>0) - A.row(i).tail(n-l).setZero(); - if (g != Scalar(0.0)) - { - g = Scalar(1.0)/g; - if (m-l) - { - for (j=l; j<n; j++) - { - s = A.col(j).tail(m-l).dot(A.col(i).tail(m-l)); - f = (s/A(i,i))*g; - A.col(j).tail(m-i) += f * A.col(i).tail(m-i); - } - } - A.col(i).tail(m-i) *= g; - } - else - A.col(i).tail(m-i).setZero(); - ++A(i,i); - } - // Diagonalization of the bidiagonal form: Loop over - // singular values, and over allowed iterations. - for (k=n-1; k>=0; k--) - { - for (its=0; its<max_iters; its++) - { - flag = true; - for (l=k; l>=0; l--) - { - // Test for splitting. - nm = l-1; - // Note that rv1[1] is always zero. - //if ((double)(ei_abs(rv1[l])+anorm) == anorm) - if (l==0 || ei_abs(m_rv1[l]) <= eps*anorm) - { - flag = false; - break; - } - //if ((double)(ei_abs(W[nm])+anorm) == anorm) - if (ei_abs(W[nm]) <= eps*anorm) - break; - } - if (flag) - { - c = 0.0; //Cancellation of rv1[l], if l > 0. - s = 1.0; - for (i=l ;i<k+1; i++) - { - f = s*m_rv1[i]; - m_rv1[i] = c*m_rv1[i]; - //if ((double)(ei_abs(f)+anorm) == anorm) - if (ei_abs(f) <= eps*anorm) - break; - g = W[i]; - h = pythag(f,g); - W[i] = h; - h = Scalar(1.0)/h; - c = g*h; - s = -f*h; - V.applyOnTheRight(i,nm,PlanarRotation<Scalar>(c,s)); - } - } - z = W[k]; - if (l == k) //Convergence. - { - if (z < 0.0) { // Singular value is made nonnegative. - W[k] = -z; - V.col(k) = -V.col(k); - } - break; - } - if (its+1 == max_iters) - { - convergence = false; - } - x = W[l]; // Shift from bottom 2-by-2 minor. - nm = k-1; - y = W[nm]; - g = m_rv1[nm]; - h = m_rv1[k]; - f = ((y-z)*(y+z) + (g-h)*(g+h))/(Scalar(2.0)*h*y); - g = pythag(f,1.0); - f = ((x-z)*(x+z) + h*((y/(f+sign(g,f)))-h))/x; - c = s = 1.0; - //Next QR transformation: - for (j=l; j<=nm; j++) - { - i = j+1; - g = m_rv1[i]; - y = W[i]; - h = s*g; - g = c*g; - - z = pythag(f,h); - m_rv1[j] = z; - c = f/z; - s = h/z; - f = x*c + g*s; - g = g*c - x*s; - h = y*s; - y *= c; - V.applyOnTheRight(i,j,PlanarRotation<Scalar>(c,s)); - - z = pythag(f,h); - W[j] = z; - // Rotation can be arbitrary if z = 0. - if (z!=Scalar(0)) - { - z = Scalar(1.0)/z; - c = f*z; - s = h*z; - } - f = c*g + s*y; - x = c*y - s*g; - A.applyOnTheRight(i,j,PlanarRotation<Scalar>(c,s)); - } - m_rv1[l] = 0.0; - m_rv1[k] = f; - W[k] = x; - } - } - - // sort the singular values: - { - for (Index i=0; i<n; i++) - { - Index k; - W.tail(n-i).maxCoeff(&k); - if (k != 0) - { - k += i; - std::swap(W[k],W[i]); - A.col(i).swap(A.col(k)); - V.col(i).swap(V.col(k)); - } - } - } - m_matU.leftCols(n) = A; - m_matU.rightCols(m-n).setZero(); - - // Gram Schmidt orthogonalization to fill up U - for (int col = A.cols(); col < A.rows(); ++col) - { - typename MatrixUType::ColXpr colVec = m_matU.col(col); - colVec(col) = 1; - for (int prevCol = 0; prevCol < col; ++prevCol) - { - typename MatrixUType::ColXpr prevColVec = m_matU.col(prevCol); - colVec -= colVec.dot(prevColVec)*prevColVec; - } - // Here we can run into troubles when colVec.norm() = 0. - m_matU.col(col) = colVec.normalized(); - } - - m_isInitialized = true; - return *this; -} - -template<typename _MatrixType, typename Rhs> -struct ei_solve_retval<SVD<_MatrixType>, Rhs> - : ei_solve_retval_base<SVD<_MatrixType>, Rhs> -{ - EIGEN_MAKE_SOLVE_HELPERS(SVD<_MatrixType>,Rhs) - - template<typename Dest> void evalTo(Dest& dst) const - { - ei_assert(rhs().rows() == dec().rows()); - - for (Index j=0; j<cols(); ++j) - { - Matrix<Scalar,MatrixType::RowsAtCompileTime,1> aux = dec().matrixU().adjoint() * rhs().col(j); - - for (Index i = 0; i < dec().singularValues().size(); ++i) - { - Scalar si = dec().singularValues().coeff(i); - if(si == RealScalar(0)) - aux.coeffRef(i) = Scalar(0); - else - aux.coeffRef(i) /= si; - } - aux.tail(aux.size() - dec().singularValues().size()).setZero(); - - const Index minsize = std::min(dec().rows(),dec().cols()); - dst.col(j).head(minsize) = aux.head(minsize); - if(dec().cols()>dec().rows()) dst.col(j).tail(cols()-minsize).setZero(); - dst.col(j) = dec().matrixV() * dst.col(j); - } - } -}; - -/** Computes the polar decomposition of the matrix, as a product unitary x positive. - * - * If either pointer is zero, the corresponding computation is skipped. - * - * Only for square matrices. - * - * \sa computePositiveUnitary(), computeRotationScaling() - */ -template<typename MatrixType> -template<typename UnitaryType, typename PositiveType> -void SVD<MatrixType>::computeUnitaryPositive(UnitaryType *unitary, - PositiveType *positive) const -{ - ei_assert(m_isInitialized && "SVD is not initialized."); - ei_assert(m_matU.cols() == m_matV.cols() && "Polar decomposition is only for square matrices"); - if(unitary) *unitary = m_matU * m_matV.adjoint(); - if(positive) *positive = m_matV * m_sigma.asDiagonal() * m_matV.adjoint(); -} - -/** Computes the polar decomposition of the matrix, as a product positive x unitary. - * - * If either pointer is zero, the corresponding computation is skipped. - * - * Only for square matrices. - * - * \sa computeUnitaryPositive(), computeRotationScaling() - */ -template<typename MatrixType> -template<typename UnitaryType, typename PositiveType> -void SVD<MatrixType>::computePositiveUnitary(UnitaryType *positive, - PositiveType *unitary) const -{ - ei_assert(m_isInitialized && "SVD is not initialized."); - ei_assert(m_matU.rows() == m_matV.rows() && "Polar decomposition is only for square matrices"); - if(unitary) *unitary = m_matU * m_matV.adjoint(); - if(positive) *positive = m_matU * m_sigma.asDiagonal() * m_matU.adjoint(); -} - -/** decomposes the matrix as a product rotation x scaling, the scaling being - * not necessarily positive. - * - * If either pointer is zero, the corresponding computation is skipped. - * - * This method requires the Geometry module. - * - * \sa computeScalingRotation(), computeUnitaryPositive() - */ -template<typename MatrixType> -template<typename RotationType, typename ScalingType> -void SVD<MatrixType>::computeRotationScaling(RotationType *rotation, ScalingType *scaling) const -{ - ei_assert(m_isInitialized && "SVD is not initialized."); - ei_assert(m_matU.rows() == m_matV.rows() && "Polar decomposition is only for square matrices"); - Scalar x = (m_matU * m_matV.adjoint()).determinant(); // so x has absolute value 1 - Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> sv(m_sigma); - sv.coeffRef(0) *= x; - if(scaling) scaling->lazyAssign(m_matV * sv.asDiagonal() * m_matV.adjoint()); - if(rotation) - { - MatrixType m(m_matU); - m.col(0) /= x; - rotation->lazyAssign(m * m_matV.adjoint()); - } -} - -/** decomposes the matrix as a product scaling x rotation, the scaling being - * not necessarily positive. - * - * If either pointer is zero, the corresponding computation is skipped. - * - * This method requires the Geometry module. - * - * \sa computeRotationScaling(), computeUnitaryPositive() - */ -template<typename MatrixType> -template<typename ScalingType, typename RotationType> -void SVD<MatrixType>::computeScalingRotation(ScalingType *scaling, RotationType *rotation) const -{ - ei_assert(m_isInitialized && "SVD is not initialized."); - ei_assert(m_matU.rows() == m_matV.rows() && "Polar decomposition is only for square matrices"); - Scalar x = (m_matU * m_matV.adjoint()).determinant(); // so x has absolute value 1 - Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> sv(m_sigma); - sv.coeffRef(0) *= x; - if(scaling) scaling->lazyAssign(m_matU * sv.asDiagonal() * m_matU.adjoint()); - if(rotation) - { - MatrixType m(m_matU); - m.col(0) /= x; - rotation->lazyAssign(m * m_matV.adjoint()); - } -} - - -/** \svd_module - * \returns the SVD decomposition of \c *this - */ -template<typename Derived> -inline SVD<typename MatrixBase<Derived>::PlainObject> -MatrixBase<Derived>::svd() const -{ - return SVD<PlainObject>(derived()); -} - -#endif // EIGEN_SVD_H diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 052d4ec13..e774c874d 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -90,7 +90,6 @@ ei_add_test(schur_complex) ei_add_test(eigensolver_selfadjoint " " "${GSL_LIBRARIES}") ei_add_test(eigensolver_generic " " "${GSL_LIBRARIES}") ei_add_test(eigensolver_complex) -ei_add_test(svd) ei_add_test(jacobi) ei_add_test(jacobisvd) ei_add_test(geo_orthomethods) diff --git a/test/nomalloc.cpp b/test/nomalloc.cpp index c1965a8ac..b4c0ff741 100644 --- a/test/nomalloc.cpp +++ b/test/nomalloc.cpp @@ -125,8 +125,7 @@ void ctms_decompositions() Eigen::FullPivHouseholderQR<Matrix> fpQR; fpQR.compute(A); // SVD module - Eigen::JacobiSVD<Matrix> jSVD; jSVD.compute(A); - Eigen::SVD<Matrix> svd; svd.compute(A); + Eigen::JacobiSVD<Matrix> jSVD; jSVD.compute(A, ComputeFullU | ComputeFullV); } void test_nomalloc() diff --git a/test/svd.cpp b/test/svd.cpp deleted file mode 100644 index 1452e39bf..000000000 --- a/test/svd.cpp +++ /dev/null @@ -1,123 +0,0 @@ -// This file is part of Eigen, a lightweight C++ template library -// for linear algebra. -// -// Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr> -// -// Eigen is free software; you can redistribute it and/or -// modify it under the terms of the GNU Lesser General Public -// License as published by the Free Software Foundation; either -// version 3 of the License, or (at your option) any later version. -// -// Alternatively, you can redistribute it and/or -// modify it under the terms of the GNU General Public License as -// published by the Free Software Foundation; either version 2 of -// the License, or (at your option) any later version. -// -// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY -// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the -// GNU General Public License for more details. -// -// You should have received a copy of the GNU Lesser General Public -// License and a copy of the GNU General Public License along with -// Eigen. If not, see <http://www.gnu.org/licenses/>. - -#include "main.h" -#include <Eigen/SVD> -#include <Eigen/LU> - -template<typename MatrixType> void svd(const MatrixType& m) -{ - /* this test covers the following files: - SVD.h - */ - typename MatrixType::Index rows = m.rows(); - typename MatrixType::Index cols = m.cols(); - - typedef typename MatrixType::Scalar Scalar; - typedef typename NumTraits<Scalar>::Real RealScalar; - MatrixType a = MatrixType::Random(rows,cols); - Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> x(cols,1), x2(cols,1); - - { - SVD<MatrixType> svd(a); - MatrixType sigma = MatrixType::Zero(rows,cols); - MatrixType matU = MatrixType::Zero(rows,rows); - MatrixType matV = MatrixType::Zero(cols,cols); - - sigma.diagonal() = svd.singularValues(); - matU = svd.matrixU(); - VERIFY_IS_UNITARY(matU); - matV = svd.matrixV(); - VERIFY_IS_UNITARY(matV); - VERIFY_IS_APPROX(a, matU * sigma * matV.transpose()); - } - - - if (rows>=cols) - { - SVD<MatrixType> svd(a); - Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> b = Matrix<Scalar, MatrixType::ColsAtCompileTime, 1>::Random(rows,1); - Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> x = svd.solve(b); - // evaluate normal equation which works also for least-squares solutions - VERIFY_IS_APPROX(a.adjoint()*a*x,a.adjoint()*b); - } - - - if(rows==cols) - { - SVD<MatrixType> svd(a); - MatrixType unitary, positive; - svd.computeUnitaryPositive(&unitary, &positive); - VERIFY_IS_APPROX(unitary * unitary.adjoint(), MatrixType::Identity(unitary.rows(),unitary.rows())); - VERIFY_IS_APPROX(positive, positive.adjoint()); - for(int i = 0; i < rows; i++) VERIFY(positive.diagonal()[i] >= 0); // cheap necessary (not sufficient) condition for positivity - VERIFY_IS_APPROX(unitary*positive, a); - - svd.computePositiveUnitary(&positive, &unitary); - VERIFY_IS_APPROX(unitary * unitary.adjoint(), MatrixType::Identity(unitary.rows(),unitary.rows())); - VERIFY_IS_APPROX(positive, positive.adjoint()); - for(int i = 0; i < rows; i++) VERIFY(positive.diagonal()[i] >= 0); // cheap necessary (not sufficient) condition for positivity - VERIFY_IS_APPROX(positive*unitary, a); - } -} - -template<typename MatrixType> void svd_verify_assert() -{ - MatrixType tmp; - - SVD<MatrixType> svd; - VERIFY_RAISES_ASSERT(svd.solve(tmp)) - VERIFY_RAISES_ASSERT(svd.matrixU()) - VERIFY_RAISES_ASSERT(svd.singularValues()) - VERIFY_RAISES_ASSERT(svd.matrixV()) - VERIFY_RAISES_ASSERT(svd.computeUnitaryPositive(&tmp,&tmp)) - VERIFY_RAISES_ASSERT(svd.computePositiveUnitary(&tmp,&tmp)) - VERIFY_RAISES_ASSERT(svd.computeRotationScaling(&tmp,&tmp)) - VERIFY_RAISES_ASSERT(svd.computeScalingRotation(&tmp,&tmp)) - - VERIFY_RAISES_ASSERT(SVD<MatrixXf>(10, 20)) -} - -void test_svd() -{ - for(int i = 0; i < g_repeat; i++) - { - CALL_SUBTEST_1( svd(Matrix3f()) ); - CALL_SUBTEST_2( svd(Matrix4d()) ); - - int cols = ei_random<int>(2,50); - int rows = cols + ei_random<int>(0,50); - CALL_SUBTEST_3( svd(MatrixXf(rows,cols)) ); - CALL_SUBTEST_4( svd(MatrixXd(rows,cols)) ); - - //complex are not implemented yet - //CALL_SUBTEST(svd(MatrixXcd(6,6)) ); - //CALL_SUBTEST(svd(MatrixXcf(3,3)) ); - } - - CALL_SUBTEST_1( svd_verify_assert<Matrix3f>() ); - CALL_SUBTEST_2( svd_verify_assert<Matrix4d>() ); - CALL_SUBTEST_3( svd_verify_assert<MatrixXf>() ); - CALL_SUBTEST_4( svd_verify_assert<MatrixXd>() ); -} |