GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-04-24 02:29:33 +08:00
Code Issues Packages Projects Releases Wiki Activity

Allow user to compute only the eigenvalues and not the eigenvectors.

Browse Source
This commit is contained in:
Jitse Niesen 2010-05-31 18:17:47 +01:00
parent 609941380a
commit 8dc947821b
11 changed files with 235 additions and 155 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

127
Eigen/src/Eigenvalues/ComplexEigenSolver.h
View File

@ -56,7 +56,10 @@
template<typename _MatrixType> class ComplexEigenSolver
{
public:
/** \brief Synonym for the template parameter \p _MatrixType. */
typedef _MatrixType MatrixType;
enum {
RowsAtCompileTime = MatrixType::RowsAtCompileTime,
ColsAtCompileTime = MatrixType::ColsAtCompileTime,
@ -65,12 +68,12 @@ template<typename _MatrixType> class ComplexEigenSolver
MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
};
/** \brief Scalar type for matrices of type \p _MatrixType. */
/** \brief Scalar type for matrices of type #MatrixType. */
typedef typename MatrixType::Scalar Scalar;
typedef typename NumTraits<Scalar>::Real RealScalar;
typedef typename MatrixType::Index Index;
/** \brief Complex scalar type for \p _MatrixType.
/** \brief Complex scalar type for #MatrixType.
*
* This is \c std::complex<Scalar> if #Scalar is real (e.g.,
* \c float or \c double) and just \c Scalar if #Scalar is
@ -81,14 +84,14 @@ template<typename _MatrixType> class ComplexEigenSolver
/** \brief Type for vector of eigenvalues as returned by eigenvalues().
*
* This is a column vector with entries of type #ComplexScalar.
* The length of the vector is the size of \p _MatrixType.
* The length of the vector is the size of #MatrixType.
*/
typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options, MaxColsAtCompileTime, 1> EigenvalueType;
/** \brief Type for matrix of eigenvectors as returned by eigenvectors().
*
* This is a square matrix with entries of type #ComplexScalar.
* The size is the same as the size of \p _MatrixType.
* The size is the same as the size of #MatrixType.
*/
typedef Matrix<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, ColsAtCompileTime> EigenvectorType;
@ -102,6 +105,7 @@ template<typename _MatrixType> class ComplexEigenSolver
m_eivalues(),
m_schur(),
m_isInitialized(false),
m_eigenvectorsOk(false),
m_matX()
{}
@ -116,40 +120,46 @@ template<typename _MatrixType> class ComplexEigenSolver
m_eivalues(size),
m_schur(size),
m_isInitialized(false),
m_eigenvectorsOk(false),
m_matX(size, size)
{}
/** \brief Constructor; computes eigendecomposition of given matrix.
*
* \param[in] matrix Square matrix whose eigendecomposition is to be computed.
* \param[in] computeEigenvectors If true, both the eigenvectors and the
* eigenvalues are computed; if false, only the eigenvalues are
* computed.
*
* This constructor calls compute() to compute the eigendecomposition.
*/
ComplexEigenSolver(const MatrixType& matrix)
ComplexEigenSolver(const MatrixType& matrix, bool computeEigenvectors = true)
: m_eivec(matrix.rows(),matrix.cols()),
m_eivalues(matrix.cols()),
m_schur(matrix.rows()),
m_isInitialized(false),
m_eigenvectorsOk(false),
m_matX(matrix.rows(),matrix.cols())
{
compute(matrix);
compute(matrix, computeEigenvectors);
}
/** \brief Returns the eigenvectors of given matrix.
*
* \returns A const reference to the matrix whose columns are the eigenvectors.
*
* It is assumed that either the constructor
* ComplexEigenSolver(const MatrixType& matrix) or the member
* function compute(const MatrixType& matrix) has been called
* before to compute the eigendecomposition of a matrix. This
* function returns a matrix whose columns are the
* eigenvectors. Column \f$ k \f$ is an eigenvector
* corresponding to eigenvalue number \f$ k \f$ as returned by
* eigenvalues(). The eigenvectors are normalized to have
* (Euclidean) norm equal to one. The matrix returned by this
* function is the matrix \f$ V \f$ in the eigendecomposition \f$
* A = V D V^{-1} \f$, if it exists.
* \pre Either the constructor
* ComplexEigenSolver(const MatrixType& matrix, bool) or the member
* function compute(const MatrixType& matrix, bool) has been called before
* to compute the eigendecomposition of a matrix, and
* \p computeEigenvectors was set to true (the default).
*
* This function returns a matrix whose columns are the eigenvectors. Column
* \f$ k \f$ is an eigenvector corresponding to eigenvalue number \f$ k
* \f$ as returned by eigenvalues(). The eigenvectors are normalized to
* have (Euclidean) norm equal to one. The matrix returned by this
* function is the matrix \f$ V \f$ in the eigendecomposition \f$ A = V D
* V^{-1} \f$, if it exists.
*
* Example: \include ComplexEigenSolver_eigenvectors.cpp
* Output: \verbinclude ComplexEigenSolver_eigenvectors.out
@ -157,6 +167,7 @@ template<typename _MatrixType> class ComplexEigenSolver
const EigenvectorType& eigenvectors() const
{
ei_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
ei_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
return m_eivec;
}
@ -164,11 +175,12 @@ template<typename _MatrixType> class ComplexEigenSolver
*
* \returns A const reference to the column vector containing the eigenvalues.
*
* It is assumed that either the constructor
* ComplexEigenSolver(const MatrixType& matrix) or the member
* function compute(const MatrixType& matrix) has been called
* before to compute the eigendecomposition of a matrix. This
* function returns a column vector containing the
* \pre Either the constructor
* ComplexEigenSolver(const MatrixType& matrix, bool) or the member
* function compute(const MatrixType& matrix, bool) has been called before
* to compute the eigendecomposition of a matrix.
*
* This function returns a column vector containing the
* eigenvalues. Eigenvalues are repeated according to their
* algebraic multiplicity, so there are as many eigenvalues as
* rows in the matrix.
@ -185,10 +197,14 @@ template<typename _MatrixType> class ComplexEigenSolver
/** \brief Computes eigendecomposition of given matrix.
*
* \param[in] matrix Square matrix whose eigendecomposition is to be computed.
* \param[in] computeEigenvectors If true, both the eigenvectors and the
* eigenvalues are computed; if false, only the eigenvalues are
* computed.
*
* This function computes the eigenvalues and eigenvectors of \p
* matrix. The eigenvalues() and eigenvectors() functions can be
* used to retrieve the computed eigendecomposition.
* This function computes the eigenvalues of the complex matrix \p matrix.
* The eigenvalues() function can be used to retrieve them. If
* \p computeEigenvectors is true, then the eigenvectors are also computed
* and can be retrieved by calling eigenvectors().
*
* The matrix is first reduced to Schur form using the
* ComplexSchur class. The Schur decomposition is then used to
@ -201,19 +217,20 @@ template<typename _MatrixType> class ComplexEigenSolver
* Example: \include ComplexEigenSolver_compute.cpp
* Output: \verbinclude ComplexEigenSolver_compute.out
*/
void compute(const MatrixType& matrix);
void compute(const MatrixType& matrix, bool computeEigenvectors = true);
protected:
EigenvectorType m_eivec;
EigenvalueType m_eivalues;
ComplexSchur<MatrixType> m_schur;
bool m_isInitialized;
bool m_eigenvectorsOk;
EigenvectorType m_matX;
};
template<typename MatrixType>
void ComplexEigenSolver<MatrixType>::compute(const MatrixType& matrix)
void ComplexEigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvectors)
{
// this code is inspired from Jampack
assert(matrix.cols() == matrix.rows());
@ -222,40 +239,45 @@ void ComplexEigenSolver<MatrixType>::compute(const MatrixType& matrix)
// Step 1: Do a complex Schur decomposition, A = U T U^*
// The eigenvalues are on the diagonal of T.
m_schur.compute(matrix);
m_schur.compute(matrix, computeEigenvectors);
m_eivalues = m_schur.matrixT().diagonal();
// Step 2: Compute X such that T = X D X^(-1), where D is the diagonal of T.
// The matrix X is unit triangular.
m_matX = EigenvectorType::Zero(n, n);
for(Index k=n-1 ; k>=0 ; k--)
if(computeEigenvectors)
{
m_matX.coeffRef(k,k) = ComplexScalar(1.0,0.0);
// Compute X(i,k) using the (i,k) entry of the equation X T = D X
for(Index i=k-1 ; i>=0 ; i--)
// Step 2: Compute X such that T = X D X^(-1), where D is the diagonal of T.
// The matrix X is unit triangular.
m_matX = EigenvectorType::Zero(n, n);
for(Index k=n-1 ; k>=0 ; k--)
{
m_matX.coeffRef(i,k) = -m_schur.matrixT().coeff(i,k);
if(k-i-1>0)
m_matX.coeffRef(i,k) -= (m_schur.matrixT().row(i).segment(i+1,k-i-1) * m_matX.col(k).segment(i+1,k-i-1)).value();
ComplexScalar z = m_schur.matrixT().coeff(i,i) - m_schur.matrixT().coeff(k,k);
if(z==ComplexScalar(0))
m_matX.coeffRef(k,k) = ComplexScalar(1.0,0.0);
// Compute X(i,k) using the (i,k) entry of the equation X T = D X
for(Index i=k-1 ; i>=0 ; i--)
{
// If the i-th and k-th eigenvalue are equal, then z equals 0.
// Use a small value instead, to prevent division by zero.
ei_real_ref(z) = NumTraits<RealScalar>::epsilon() * matrixnorm;
m_matX.coeffRef(i,k) = -m_schur.matrixT().coeff(i,k);
if(k-i-1>0)
m_matX.coeffRef(i,k) -= (m_schur.matrixT().row(i).segment(i+1,k-i-1) * m_matX.col(k).segment(i+1,k-i-1)).value();
ComplexScalar z = m_schur.matrixT().coeff(i,i) - m_schur.matrixT().coeff(k,k);
if(z==ComplexScalar(0))
{
// If the i-th and k-th eigenvalue are equal, then z equals 0.
// Use a small value instead, to prevent division by zero.
ei_real_ref(z) = NumTraits<RealScalar>::epsilon() * matrixnorm;
}
m_matX.coeffRef(i,k) = m_matX.coeff(i,k) / z;
}
m_matX.coeffRef(i,k) = m_matX.coeff(i,k) / z;
}
// Step 3: Compute V as V = U X; now A = U T U^* = U X D X^(-1) U^* = V D V^(-1)
m_eivec.noalias() = m_schur.matrixU() * m_matX;
// .. and normalize the eigenvectors
for(Index k=0 ; k<n ; k++)
{
m_eivec.col(k).normalize();
}
}
// Step 3: Compute V as V = U X; now A = U T U^* = U X D X^(-1) U^* = V D V^(-1)
m_eivec.noalias() = m_schur.matrixU() * m_matX;
// .. and normalize the eigenvectors
for(Index k=0 ; k<n ; k++)
{
m_eivec.col(k).normalize();
}
m_isInitialized = true;
m_eigenvectorsOk = computeEigenvectors;
// Step 4: Sort the eigenvalues
for (Index i=0; i<n; i++)
@ -266,7 +288,8 @@ void ComplexEigenSolver<MatrixType>::compute(const MatrixType& matrix)
{
k += i;
std::swap(m_eivalues[k],m_eivalues[i]);
m_eivec.col(i).swap(m_eivec.col(k));
if(computeEigenvectors)
m_eivec.col(i).swap(m_eivec.col(k));
}
}
}

100
Eigen/src/Eigenvalues/EigenSolver.h
View File

@ -57,16 +57,16 @@
* this variant of the eigendecomposition the pseudo-eigendecomposition.
*
* Call the function compute() to compute the eigenvalues and eigenvectors of
* a given matrix. Alternatively, you can use the
* EigenSolver(const MatrixType&) constructor which computes the eigenvalues
* and eigenvectors at construction time. Once the eigenvalue and eigenvectors
* are computed, they can be retrieved with the eigenvalues() and
* a given matrix. Alternatively, you can use the
* EigenSolver(const MatrixType&, bool) constructor which computes the
* eigenvalues and eigenvectors at construction time. Once the eigenvalue and
* eigenvectors are computed, they can be retrieved with the eigenvalues() and
* eigenvectors() functions. The pseudoEigenvalueMatrix() and
* pseudoEigenvectors() methods allow the construction of the
* pseudo-eigendecomposition.
*
* The documentation for EigenSolver(const MatrixType&) contains an example of
* the typical use of this class.
* The documentation for EigenSolver(const MatrixType&, bool) contains an
* example of the typical use of this class.
*
* \note The implementation is adapted from
* <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a> (public domain).
@ -78,7 +78,9 @@ template<typename _MatrixType> class EigenSolver
{
public:
/** \brief Synonym for the template parameter \p _MatrixType. */
typedef _MatrixType MatrixType;
enum {
RowsAtCompileTime = MatrixType::RowsAtCompileTime,
ColsAtCompileTime = MatrixType::ColsAtCompileTime,
@ -87,12 +89,12 @@ template<typename _MatrixType> class EigenSolver
MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
};
/** \brief Scalar type for matrices of type \p _MatrixType. */
/** \brief Scalar type for matrices of type #MatrixType. */
typedef typename MatrixType::Scalar Scalar;
typedef typename NumTraits<Scalar>::Real RealScalar;
typedef typename MatrixType::Index Index;
/** \brief Complex scalar type for \p _MatrixType.
/** \brief Complex scalar type for #MatrixType.
*
* This is \c std::complex<Scalar> if #Scalar is real (e.g.,
* \c float or \c double) and just \c Scalar if #Scalar is
@ -103,27 +105,27 @@ template<typename _MatrixType> class EigenSolver
/** \brief Type for vector of eigenvalues as returned by eigenvalues().
*
* This is a column vector with entries of type #ComplexScalar.
* The length of the vector is the size of \p _MatrixType.
* The length of the vector is the size of #MatrixType.
*/
typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> EigenvalueType;
/** \brief Type for matrix of eigenvectors as returned by eigenvectors().
*
* This is a square matrix with entries of type #ComplexScalar.
* The size is the same as the size of \p _MatrixType.
* The size is the same as the size of #MatrixType.
*/
typedef Matrix<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime> EigenvectorsType;
/** \brief Default constructor.
*
* The default constructor is useful in cases in which the user intends to
* perform decompositions via EigenSolver::compute(const MatrixType&).
* perform decompositions via EigenSolver::compute(const MatrixType&, bool).
*
* \sa compute() for an example.
*/
EigenSolver() : m_eivec(), m_eivalues(), m_isInitialized(false), m_realSchur(), m_matT(), m_tmp() {}
/** \brief Default Constructor with memory preallocation
/** \brief Default constructor with memory preallocation
*
* Like the default constructor but with preallocation of the internal data
* according to the specified problem \a size.
@ -133,6 +135,7 @@ template<typename _MatrixType> class EigenSolver
: m_eivec(size, size),
m_eivalues(size),
m_isInitialized(false),
m_eigenvectorsOk(false),
m_realSchur(size),
m_matT(size, size),
m_tmp(size)
@ -141,6 +144,9 @@ template<typename _MatrixType> class EigenSolver
/** \brief Constructor; computes eigendecomposition of given matrix.
*
* \param[in] matrix Square matrix whose eigendecomposition is to be computed.
* \param[in] computeEigenvectors If true, both the eigenvectors and the
* eigenvalues are computed; if false, only the eigenvalues are
* computed.
*
* This constructor calls compute() to compute the eigenvalues
* and eigenvectors.
@ -150,23 +156,26 @@ template<typename _MatrixType> class EigenSolver
*
* \sa compute()
*/
EigenSolver(const MatrixType& matrix)
EigenSolver(const MatrixType& matrix, bool computeEigenvectors = true)
: m_eivec(matrix.rows(), matrix.cols()),
m_eivalues(matrix.cols()),
m_isInitialized(false),
m_eigenvectorsOk(false),
m_realSchur(matrix.cols()),
m_matT(matrix.rows(), matrix.cols()),
m_tmp(matrix.cols())
{
compute(matrix);
compute(matrix, computeEigenvectors);
}
/** \brief Returns the eigenvectors of given matrix.
*
* \returns %Matrix whose columns are the (possibly complex) eigenvectors.
*
* \pre Either the constructor EigenSolver(const MatrixType&) or the
* member function compute(const MatrixType&) has been called before.
* \pre Either the constructor
* EigenSolver(const MatrixType&,bool) or the member function
* compute(const MatrixType&, bool) has been called before, and
* \p computeEigenvectors was set to true (the default).
*
* Column \f$ k \f$ of the returned matrix is an eigenvector corresponding
* to eigenvalue number \f$ k \f$ as returned by eigenvalues(). The
@ -185,9 +194,10 @@ template<typename _MatrixType> class EigenSolver
*
* \returns Const reference to matrix whose columns are the pseudo-eigenvectors.
*
* \pre Either the constructor EigenSolver(const MatrixType&) or
* the member function compute(const MatrixType&) has been called
* before.
* \pre Either the constructor
* EigenSolver(const MatrixType&,bool) or the member function
* compute(const MatrixType&, bool) has been called before, and
* \p computeEigenvectors was set to true (the default).
*
* The real matrix \f$ V \f$ returned by this function and the
* block-diagonal matrix \f$ D \f$ returned by pseudoEigenvalueMatrix()
@ -201,6 +211,7 @@ template<typename _MatrixType> class EigenSolver
const MatrixType& pseudoEigenvectors() const
{
ei_assert(m_isInitialized && "EigenSolver is not initialized.");
ei_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
return m_eivec;
}
@ -208,8 +219,9 @@ template<typename _MatrixType> class EigenSolver
*
* \returns A block-diagonal matrix.
*
* \pre Either the constructor EigenSolver(const MatrixType&) or the
* member function compute(const MatrixType&) has been called before.
* \pre Either the constructor
* EigenSolver(const MatrixType&,bool) or the member function
* compute(const MatrixType&, bool) has been called before.
*
* The matrix \f$ D \f$ returned by this function is real and
* block-diagonal. The blocks on the diagonal are either 1-by-1 or 2-by-2
@ -226,8 +238,9 @@ template<typename _MatrixType> class EigenSolver
*
* \returns A const reference to the column vector containing the eigenvalues.
*
* \pre Either the constructor EigenSolver(const MatrixType&) or the
* member function compute(const MatrixType&) has been called before.
* \pre Either the constructor
* EigenSolver(const MatrixType&,bool) or the member function
* compute(const MatrixType&, bool) has been called before.
*
* The eigenvalues are repeated according to their algebraic multiplicity,
* so there are as many eigenvalues as rows in the matrix.
@ -247,34 +260,40 @@ template<typename _MatrixType> class EigenSolver
/** \brief Computes eigendecomposition of given matrix.
*
* \param[in] matrix Square matrix whose eigendecomposition is to be computed.
* \param[in] computeEigenvectors If true, both the eigenvectors and the
* eigenvalues are computed; if false, only the eigenvalues are
* computed.
* \returns Reference to \c *this
*
* This function computes the eigenvalues and eigenvectors of \p matrix.
* The eigenvalues() and eigenvectors() functions can be used to retrieve
* the computed eigendecomposition.
* This function computes the eigenvalues of the real matrix \p matrix.
* The eigenvalues() function can be used to retrieve them. If
* \p computeEigenvectors is true, then the eigenvectors are also computed
* and can be retrieved by calling eigenvectors().
*
* The matrix is first reduced to real Schur form using the RealSchur
* class. The Schur decomposition is then used to compute the eigenvalues
* and eigenvectors.
*
* The cost of the computation is dominated by the cost of the Schur
* decomposition, which is very approximately \f$ 25n^3 \f$ where
* \f$ n \f$ is the size of the matrix.
* The cost of the computation is dominated by the cost of the
* Schur decomposition, which is very approximately \f$ 25n^3 \f$
* (where \f$ n \f$ is the size of the matrix) if \p computeEigenvectors
* is true, and \f$ 10n^3 \f$ if \p computeEigenvectors is false.
*
* This method reuses of the allocated data in the EigenSolver object.
*
* Example: \include EigenSolver_compute.cpp
* Output: \verbinclude EigenSolver_compute.out
*/
EigenSolver& compute(const MatrixType& matrix);
EigenSolver& compute(const MatrixType& matrix, bool computeEigenvectors = true);
private:
void computeEigenvectors();
void doComputeEigenvectors();
protected:
MatrixType m_eivec;
EigenvalueType m_eivalues;
bool m_isInitialized;
bool m_eigenvectorsOk;
RealSchur<MatrixType> m_realSchur;
MatrixType m_matT;
@ -286,7 +305,7 @@ template<typename MatrixType>
MatrixType EigenSolver<MatrixType>::pseudoEigenvalueMatrix() const
{
ei_assert(m_isInitialized && "EigenSolver is not initialized.");
Index n = m_eivec.cols();
Index n = m_eivalues.rows();
MatrixType matD = MatrixType::Zero(n,n);
for (Index i=0; i<n; ++i)
{
@ -306,6 +325,7 @@ template<typename MatrixType>
typename EigenSolver<MatrixType>::EigenvectorsType EigenSolver<MatrixType>::eigenvectors() const
{
ei_assert(m_isInitialized && "EigenSolver is not initialized.");
ei_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
Index n = m_eivec.cols();
EigenvectorsType matV(n,n);
for (Index j=0; j<n; ++j)
@ -332,14 +352,15 @@ typename EigenSolver<MatrixType>::EigenvectorsType EigenSolver<MatrixType>::eige
}
template<typename MatrixType>
EigenSolver<MatrixType>& EigenSolver<MatrixType>::compute(const MatrixType& matrix)
EigenSolver<MatrixType>& EigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvectors)
{
assert(matrix.cols() == matrix.rows());
// Reduce to real Schur form.
m_realSchur.compute(matrix);
m_realSchur.compute(matrix, computeEigenvectors);
m_matT = m_realSchur.matrixT();
m_eivec = m_realSchur.matrixU();
if (computeEigenvectors)
m_eivec = m_realSchur.matrixU();
// Compute eigenvalues from matT
m_eivalues.resize(matrix.cols());
@ -362,9 +383,12 @@ EigenSolver<MatrixType>& EigenSolver<MatrixType>::compute(const MatrixType& matr
}
// Compute eigenvectors.
computeEigenvectors();
if (computeEigenvectors)
doComputeEigenvectors();
m_isInitialized = true;
m_eigenvectorsOk = computeEigenvectors;
return *this;
}
@ -389,7 +413,7 @@ std::complex<Scalar> cdiv(Scalar xr, Scalar xi, Scalar yr, Scalar yi)
template<typename MatrixType>
void EigenSolver<MatrixType>::computeEigenvectors()
void EigenSolver<MatrixType>::doComputeEigenvectors()
{
const Index size = m_eivec.cols();
const Scalar eps = NumTraits<Scalar>::epsilon();
@ -404,7 +428,7 @@ void EigenSolver<MatrixType>::computeEigenvectors()
// Backsubstitute to find vectors of upper triangular form
if (norm == 0.0)
{
return;
return;
}
for (Index n = size-1; n >= 0; n--)

6
Eigen/src/Eigenvalues/MatrixBaseEigenvalues.h
View File

@ -37,7 +37,7 @@ struct ei_eigenvalues_selector
{
typedef typename Derived::PlainObject PlainObject;
PlainObject m_eval(m);
return ComplexEigenSolver<PlainObject>(m_eval).eigenvalues();
return ComplexEigenSolver<PlainObject>(m_eval, false).eigenvalues();
}
};
@ -49,7 +49,7 @@ struct ei_eigenvalues_selector<Derived, false>
{
typedef typename Derived::PlainObject PlainObject;
PlainObject m_eval(m);
return EigenSolver<PlainObject>(m_eval).eigenvalues();
return EigenSolver<PlainObject>(m_eval, false).eigenvalues();
}
};
@ -101,7 +101,7 @@ SelfAdjointView<MatrixType, UpLo>::eigenvalues() const
{
typedef typename SelfAdjointView<MatrixType, UpLo>::PlainObject PlainObject;
PlainObject thisAsMatrix(*this);
return SelfAdjointEigenSolver<PlainObject>(thisAsMatrix).eigenvalues();
return SelfAdjointEigenSolver<PlainObject>(thisAsMatrix, false).eigenvalues();
}

84
Eigen/src/Eigenvalues/RealSchur.h
View File

@ -50,13 +50,13 @@
* the eigendecomposition of a matrix.
*
* Call the function compute() to compute the real Schur decomposition of a
* given matrix. Alternatively, you can use the RealSchur(const MatrixType&)
* given matrix. Alternatively, you can use the RealSchur(const MatrixType&, bool)
* constructor which computes the real Schur decomposition at construction
* time. Once the decomposition is computed, you can use the matrixU() and
* matrixT() functions to retrieve the matrices U and T in the decomposition.
*
* The documentation of RealSchur(const MatrixType&) contains an example of
* the typical use of this class.
* The documentation of RealSchur(const MatrixType&, bool) contains an example
* of the typical use of this class.
*
* \note The implementation is adapted from
* <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a> (public domain).
@ -98,41 +98,46 @@ template<typename _MatrixType> class RealSchur
m_matU(size, size),
m_workspaceVector(size),
m_hess(size),
m_isInitialized(false)
m_isInitialized(false),
m_matUisUptodate(false)
{ }
/** \brief Constructor; computes real Schur decomposition of given matrix.
*
* \param[in] matrix Square matrix whose Schur decomposition is to be computed.
* \param[in] matrix Square matrix whose Schur decomposition is to be computed.
* \param[in] computeU If true, both T and U are computed; if false, only T is computed.
*
* This constructor calls compute() to compute the Schur decomposition.
*
* Example: \include RealSchur_RealSchur_MatrixType.cpp
* Output: \verbinclude RealSchur_RealSchur_MatrixType.out
*/
RealSchur(const MatrixType& matrix)
RealSchur(const MatrixType& matrix, bool computeU = true)
: m_matT(matrix.rows(),matrix.cols()),
m_matU(matrix.rows(),matrix.cols()),
m_workspaceVector(matrix.rows()),
m_hess(matrix.rows()),
m_isInitialized(false)
m_isInitialized(false),
m_matUisUptodate(false)
{
compute(matrix);
compute(matrix, computeU);
}
/** \brief Returns the orthogonal matrix in the Schur decomposition.
*
* \returns A const reference to the matrix U.
*
* \pre Either the constructor RealSchur(const MatrixType&) or the member
* function compute(const MatrixType&) has been called before to compute
* the Schur decomposition of a matrix.
* \pre Either the constructor RealSchur(const MatrixType&, bool) or the
* member function compute(const MatrixType&, bool) has been called before
* to compute the Schur decomposition of a matrix, and \p computeU was set
* to true (the default value).
*
* \sa RealSchur(const MatrixType&) for an example
* \sa RealSchur(const MatrixType&, bool) for an example
*/
const MatrixType& matrixU() const
{
ei_assert(m_isInitialized && "RealSchur is not initialized.");
ei_assert(m_matUisUptodate && "The matrix U has not been computed during the RealSchur decomposition.");
return m_matU;
}
@ -140,11 +145,11 @@ template<typename _MatrixType> class RealSchur
*
* \returns A const reference to the matrix T.
*
* \pre Either the constructor RealSchur(const MatrixType&) or the member
* function compute(const MatrixType&) has been called before to compute
* the Schur decomposition of a matrix.
* \pre Either the constructor RealSchur(const MatrixType&, bool) or the
* member function compute(const MatrixType&, bool) has been called before
* to compute the Schur decomposition of a matrix.
*
* \sa RealSchur(const MatrixType&) for an example
* \sa RealSchur(const MatrixType&, bool) for an example
*/
const MatrixType& matrixT() const
{
@ -154,19 +159,21 @@ template<typename _MatrixType> class RealSchur
/** \brief Computes Schur decomposition of given matrix.
*
* \param[in] matrix Square matrix whose Schur decomposition is to be computed.
* \param[in] matrix Square matrix whose Schur decomposition is to be computed.
* \param[in] computeU If true, both T and U are computed; if false, only T is computed.
*
* The Schur decomposition is computed by first reducing the matrix to
* Hessenberg form using the class HessenbergDecomposition. The Hessenberg
* matrix is then reduced to triangular form by performing Francis QR
* iterations with implicit double shift. The cost of computing the Schur
* decomposition depends on the number of iterations; as a rough guide, it
* may be taken to be \f$25n^3\f$ flops.
* may be taken to be \f$25n^3\f$ flops if \a computeU is true and
* \f$10n^3\f$ flops if \a computeU is false.
*
* Example: \include RealSchur_compute.cpp
* Output: \verbinclude RealSchur_compute.out
*/
void compute(const MatrixType& matrix);
void compute(const MatrixType& matrix, bool computeU = true);
private:
@ -175,38 +182,39 @@ template<typename _MatrixType> class RealSchur
ColumnVectorType m_workspaceVector;
HessenbergDecomposition<MatrixType> m_hess;
bool m_isInitialized;
bool m_matUisUptodate;
typedef Matrix<Scalar,3,1> Vector3s;
Scalar computeNormOfT();
Index findSmallSubdiagEntry(Index iu, Scalar norm);
void splitOffTwoRows(Index iu, Scalar exshift);
void splitOffTwoRows(Index iu, bool computeU, Scalar exshift);
void computeShift(Index iu, Index iter, Scalar& exshift, Vector3s& shiftInfo);
void initFrancisQRStep(Index il, Index iu, const Vector3s& shiftInfo, Index& im, Vector3s& firstHouseholderVector);
void performFrancisQRStep(Index il, Index im, Index iu, const Vector3s& firstHouseholderVector, Scalar* workspace);
void performFrancisQRStep(Index il, Index im, Index iu, bool computeU, const Vector3s& firstHouseholderVector, Scalar* workspace);
};
template<typename MatrixType>
void RealSchur<MatrixType>::compute(const MatrixType& matrix)
void RealSchur<MatrixType>::compute(const MatrixType& matrix, bool computeU)
{
assert(matrix.cols() == matrix.rows());
// Step 1. Reduce to Hessenberg form
// TODO skip Q if skipU = true
m_hess.compute(matrix);
m_matT = m_hess.matrixH();
m_matU = m_hess.matrixQ();
if (computeU)
m_matU = m_hess.matrixQ();
// Step 2. Reduce to real Schur form
m_workspaceVector.resize(m_matU.cols());
m_workspaceVector.resize(m_matT.cols());
Scalar* workspace = &m_workspaceVector.coeffRef(0);
// The matrix m_matT is divided in three parts.
// Rows 0,...,il-1 are decoupled from the rest because m_matT(il,il-1) is zero.
// Rows il,...,iu is the part we are working on (the active window).
// Rows iu+1,...,end are already brought in triangular form.
Index iu = m_matU.cols() - 1;
Index iu = m_matT.cols() - 1;
Index iter = 0; // iteration count
Scalar exshift = 0.0; // sum of exceptional shifts
Scalar norm = computeNormOfT();
@ -226,7 +234,7 @@ void RealSchur<MatrixType>::compute(const MatrixType& matrix)
}
else if (il == iu-1) // Two roots found
{
splitOffTwoRows(iu, exshift);
splitOffTwoRows(iu, computeU, exshift);
iu -= 2;
iter = 0;
}
@ -237,18 +245,19 @@ void RealSchur<MatrixType>::compute(const MatrixType& matrix)
iter = iter + 1; // (Could check iteration count here.)
Index im;
initFrancisQRStep(il, iu, shiftInfo, im, firstHouseholderVector);
performFrancisQRStep(il, im, iu, firstHouseholderVector, workspace);
performFrancisQRStep(il, im, iu, computeU, firstHouseholderVector, workspace);
}
}
m_isInitialized = true;
m_matUisUptodate = computeU;
}
/** \internal Computes and returns vector L1 norm of T */
template<typename MatrixType>
inline typename MatrixType::Scalar RealSchur<MatrixType>::computeNormOfT()
{
const Index size = m_matU.cols();
const Index size = m_matT.cols();
// FIXME to be efficient the following would requires a triangular reduxion code
// Scalar norm = m_matT.upper().cwiseAbs().sum()
// + m_matT.bottomLeftCorner(size-1,size-1).diagonal().cwiseAbs().sum();
@ -277,9 +286,9 @@ inline typename MatrixType::Index RealSchur<MatrixType>::findSmallSubdiagEntry(I
/** \internal Update T given that rows iu-1 and iu decouple from the rest. */
template<typename MatrixType>
inline void RealSchur<MatrixType>::splitOffTwoRows(Index iu, Scalar exshift)
inline void RealSchur<MatrixType>::splitOffTwoRows(Index iu, bool computeU, Scalar exshift)
{
const Index size = m_matU.cols();
const Index size = m_matT.cols();
// The eigenvalues of the 2x2 matrix [a b; c d] are
// trace +/- sqrt(discr/4) where discr = tr^2 - 4*det, tr = a + d, det = ad - bc
@ -300,7 +309,8 @@ inline void RealSchur<MatrixType>::splitOffTwoRows(Index iu, Scalar exshift)
m_matT.rightCols(size-iu+1).applyOnTheLeft(iu-1, iu, rot.adjoint());
m_matT.topRows(iu+1).applyOnTheRight(iu-1, iu, rot);
m_matT.coeffRef(iu, iu-1) = Scalar(0);
m_matU.applyOnTheRight(iu-1, iu, rot);
if (computeU)
m_matU.applyOnTheRight(iu-1, iu, rot);
}
if (iu > 1)
@ -375,12 +385,12 @@ inline void RealSchur<MatrixType>::initFrancisQRStep(Index il, Index iu, const V
/** \internal Perform a Francis QR step involving rows il:iu and columns im:iu. */
template<typename MatrixType>
inline void RealSchur<MatrixType>::performFrancisQRStep(Index il, Index im, Index iu, const Vector3s& firstHouseholderVector, Scalar* workspace)
inline void RealSchur<MatrixType>::performFrancisQRStep(Index il, Index im, Index iu, bool computeU, const Vector3s& firstHouseholderVector, Scalar* workspace)
{
assert(im >= il);
assert(im <= iu-2);
const Index size = m_matU.cols();
const Index size = m_matT.cols();
for (Index k = im; k <= iu-2; ++k)
{
@ -406,7 +416,8 @@ inline void RealSchur<MatrixType>::performFrancisQRStep(Index il, Index im, Inde
// These Householder transformations form the O(n^3) part of the algorithm
m_matT.block(k, k, 3, size-k).applyHouseholderOnTheLeft(ess, tau, workspace);
m_matT.block(0, k, std::min(iu,k+3) + 1, 3).applyHouseholderOnTheRight(ess, tau, workspace);
m_matU.block(0, k, size, 3).applyHouseholderOnTheRight(ess, tau, workspace);
if (computeU)
m_matU.block(0, k, size, 3).applyHouseholderOnTheRight(ess, tau, workspace);
}
}
@ -420,7 +431,8 @@ inline void RealSchur<MatrixType>::performFrancisQRStep(Index il, Index im, Inde
m_matT.coeffRef(iu-1, iu-2) = beta;
m_matT.block(iu-1, iu-1, 2, size-iu+1).applyHouseholderOnTheLeft(ess, tau, workspace);
m_matT.block(0, iu-1, iu+1, 2).applyHouseholderOnTheRight(ess, tau, workspace);
m_matU.block(0, iu-1, size, 2).applyHouseholderOnTheRight(ess, tau, workspace);
if (computeU)
m_matU.block(0, iu-1, size, 2).applyHouseholderOnTheRight(ess, tau, workspace);
}
// clean up pollution due to round-off errors

2
doc/snippets/ComplexEigenSolver_eigenvalues.cpp
View File

@ -1,4 +1,4 @@
MatrixXcf ones = MatrixXcf::Ones(3,3);
ComplexEigenSolver<MatrixXcf> ces(ones);
ComplexEigenSolver<MatrixXcf> ces(ones, /* computeEigenvectors = */ false);
cout << "The eigenvalues of the 3x3 matrix of ones are:"
<< endl << ces.eigenvalues() << endl;

4
doc/snippets/EigenSolver_compute.cpp
View File

@ -1,6 +1,6 @@
EigenSolver<MatrixXf> es;
MatrixXf A = MatrixXf::Random(4,4);
es.compute(A);
es.compute(A, /* computeEigenvectors = */ false);
cout << "The eigenvalues of A are: " << es.eigenvalues().transpose() << endl;
es.compute(A + MatrixXf::Identity(4,4)); // re-use es to compute eigenvalues of A+I
es.compute(A + MatrixXf::Identity(4,4), false); // re-use es to compute eigenvalues of A+I
cout << "The eigenvalues of A+I are: " << es.eigenvalues().transpose() << endl;

2
doc/snippets/EigenSolver_eigenvalues.cpp
View File

@ -1,4 +1,4 @@
MatrixXd ones = MatrixXd::Ones(3,3);
EigenSolver<MatrixXd> es(ones);
EigenSolver<MatrixXd> es(ones, false);
cout << "The eigenvalues of the 3x3 matrix of ones are:"
<< endl << es.eigenvalues() << endl;

4
doc/snippets/RealSchur_compute.cpp
View File

@ -1,6 +1,6 @@
MatrixXf A = MatrixXf::Random(4,4);
RealSchur<MatrixXf> schur(4);
schur.compute(A);
schur.compute(A, /* computeU = */ false);
cout << "The matrix T in the decomposition of A is:" << endl << schur.matrixT() << endl;
schur.compute(A.inverse());
schur.compute(A.inverse(), /* computeU = */ false);
cout << "The matrix T in the decomposition of A^(-1) is:" << endl << schur.matrixT() << endl;

24
test/eigensolver_complex.cpp
View File

@ -2,6 +2,7 @@
// for linear algebra. Eigen itself is part of the KDE project.
//
// Copyright (C) 2008-2009 Gael Guennebaud <g.gael@free.fr>
// Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
@ -66,7 +67,10 @@ template<typename MatrixType> void eigensolver(const MatrixType& m)
// Note: If MatrixType is real then a.eigenvalues() uses EigenSolver and thus
// another algorithm so results may differ slightly
verify_is_approx_upto_permutation(a.eigenvalues(), ei1.eigenvalues());
ComplexEigenSolver<MatrixType> eiNoEivecs(a, false);
VERIFY_IS_APPROX(ei1.eigenvalues(), eiNoEivecs.eigenvalues());
// Regression test for issue #66
MatrixType z = MatrixType::Zero(rows,cols);
ComplexEigenSolver<MatrixType> eiz(z);
@ -76,11 +80,15 @@ template<typename MatrixType> void eigensolver(const MatrixType& m)
VERIFY_IS_APPROX(id.operatorNorm(), RealScalar(1));
}
template<typename MatrixType> void eigensolver_verify_assert()
template<typename MatrixType> void eigensolver_verify_assert(const MatrixType& m)
{
ComplexEigenSolver<MatrixType> eig;
VERIFY_RAISES_ASSERT(eig.eigenvectors())
VERIFY_RAISES_ASSERT(eig.eigenvalues())
VERIFY_RAISES_ASSERT(eig.eigenvectors());
VERIFY_RAISES_ASSERT(eig.eigenvalues());
MatrixType a = MatrixType::Random(m.rows(),m.cols());
eig.compute(a, false);
VERIFY_RAISES_ASSERT(eig.eigenvectors());
}
void test_eigensolver_complex()
@ -92,10 +100,10 @@ void test_eigensolver_complex()
CALL_SUBTEST_4( eigensolver(Matrix3f()) );
}
CALL_SUBTEST_1( eigensolver_verify_assert<Matrix4cf>() );
CALL_SUBTEST_2( eigensolver_verify_assert<MatrixXcd>() );
CALL_SUBTEST_3(( eigensolver_verify_assert<Matrix<std::complex<float>, 1, 1> >() ));
CALL_SUBTEST_4( eigensolver_verify_assert<Matrix3f>() );
CALL_SUBTEST_1( eigensolver_verify_assert(Matrix4cf()) );
CALL_SUBTEST_2( eigensolver_verify_assert(MatrixXcd(14,14)) );
CALL_SUBTEST_3( eigensolver_verify_assert(Matrix<std::complex<float>, 1, 1>()) );
CALL_SUBTEST_4( eigensolver_verify_assert(Matrix3f()) );
// Test problem size constructors
CALL_SUBTEST_5(ComplexEigenSolver<MatrixXf>(10));

32
test/eigensolver_generic.cpp
View File

@ -2,6 +2,7 @@
// for linear algebra.
//
// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
// Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
@ -60,19 +61,26 @@ template<typename MatrixType> void eigensolver(const MatrixType& m)
ei1.eigenvectors() * ei1.eigenvalues().asDiagonal());
VERIFY_IS_APPROX(a.eigenvalues(), ei1.eigenvalues());
EigenSolver<MatrixType> eiNoEivecs(a, false);
VERIFY_IS_APPROX(ei1.eigenvalues(), eiNoEivecs.eigenvalues());
VERIFY_IS_APPROX(ei1.pseudoEigenvalueMatrix(), eiNoEivecs.pseudoEigenvalueMatrix());
MatrixType id = MatrixType::Identity(rows, cols);
VERIFY_IS_APPROX(id.operatorNorm(), RealScalar(1));
}
template<typename MatrixType> void eigensolver_verify_assert()
template<typename MatrixType> void eigensolver_verify_assert(const MatrixType& m)
{
MatrixType tmp;
EigenSolver<MatrixType> eig;
VERIFY_RAISES_ASSERT(eig.eigenvectors())
VERIFY_RAISES_ASSERT(eig.pseudoEigenvectors())
VERIFY_RAISES_ASSERT(eig.pseudoEigenvalueMatrix())
VERIFY_RAISES_ASSERT(eig.eigenvalues())
VERIFY_RAISES_ASSERT(eig.eigenvectors());
VERIFY_RAISES_ASSERT(eig.pseudoEigenvectors());
VERIFY_RAISES_ASSERT(eig.pseudoEigenvalueMatrix());
VERIFY_RAISES_ASSERT(eig.eigenvalues());
MatrixType a = MatrixType::Random(m.rows(),m.cols());
eig.compute(a, false);
VERIFY_RAISES_ASSERT(eig.eigenvectors());
VERIFY_RAISES_ASSERT(eig.pseudoEigenvectors());
}
void test_eigensolver_generic()
@ -88,11 +96,11 @@ void test_eigensolver_generic()
CALL_SUBTEST_4( eigensolver(Matrix2d()) );
}
CALL_SUBTEST_1( eigensolver_verify_assert<Matrix4f>() );
CALL_SUBTEST_2( eigensolver_verify_assert<MatrixXd>() );
CALL_SUBTEST_4( eigensolver_verify_assert<Matrix2d>() );
CALL_SUBTEST_5( eigensolver_verify_assert<MatrixXf>() );
CALL_SUBTEST_1( eigensolver_verify_assert(Matrix4f()) );
CALL_SUBTEST_2( eigensolver_verify_assert(MatrixXd(17,17)) );
CALL_SUBTEST_3( eigensolver_verify_assert(Matrix<double,1,1>()) );
CALL_SUBTEST_4( eigensolver_verify_assert(Matrix2d()) );
// Test problem size constructors
CALL_SUBTEST_6(EigenSolver<MatrixXf>(10));
CALL_SUBTEST_5(EigenSolver<MatrixXf>(10));
}

5
test/schur_real.cpp
View File

@ -73,6 +73,11 @@ template<typename MatrixType> void schur(int size = MatrixType::ColsAtCompileTim
RealSchur<MatrixType> rs2(A);
VERIFY_IS_EQUAL(rs1.matrixT(), rs2.matrixT());
VERIFY_IS_EQUAL(rs1.matrixU(), rs2.matrixU());
// Test computation of only T, not U
RealSchur<MatrixType> rsOnlyT(A, false);
VERIFY_IS_EQUAL(rs1.matrixT(), rsOnlyT.matrixT());
VERIFY_RAISES_ASSERT(rsOnlyT.matrixU());
}
void test_schur_real()