mirror of
https://gitlab.com/libeigen/eigen.git
synced 2025-08-11 11:19:02 +08:00
remove SVD class (was bad code taked from elsewhere)
Use JacobiSVD for now. We do plan to reintroduce a bidiagonalizing SVD asap.
This commit is contained in:
Benoit Jacob
parent
dbedc70012
commit
8eb0fc1e72
@ -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"
|
||||
|
||||
|
||||
@ -331,8 +331,6 @@ template<typename Derived> class MatrixBase
|
||||
|
||||
/////////// SVD module ///////////
|
||||
|
||||
SVD<PlainObject> svd() const;
|
||||
|
||||
/////////// Geometry module ///////////
|
||||
|
||||
template<typename OtherDerived>
|
||||
|
||||
@ -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;
|
||||
|
||||
@ -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
|
||||
|
||||
@ -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);
|
||||
|
||||
@ -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
|
||||
@ -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)
|
||||
|
||||
@ -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()
|
||||
|
||||
123
test/svd.cpp
123
test/svd.cpp
@ -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>() );
|
||||
}
|
||||
Loading…
x
Reference in New Issue
Block a user