mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-02-23 18:20:47 +08:00
Code Issues Packages Projects Releases Wiki Activity

Fixes in Eigensolver:

Browse Source 
* eigenvectors => pseudoEigenvectors
 * added pseudoEigenvalueMatrix
 * clear the documentation
 * added respective unit test
Still missing: a proper eigenvectors() function.
This commit is contained in:
Gael Guennebaud 2008-10-01 10:17:08 +00:00
parent 618de17bf7
commit d907cd4410
2 changed files with 104 additions and 15 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

64
Eigen/src/QR/EigenSolver.h
View File

@ -58,9 +58,50 @@ template<typename _MatrixType> class EigenSolver
compute(matrix);
}
MatrixType eigenvectors(void) const { return m_eivec; }
// TODO compute the complex eigen vectors
// MatrixType eigenvectors(void) const { return m_eivec; }
EigenvalueType eigenvalues(void) const { return m_eivalues; }
/** \returns a real matrix V of pseudo eigenvectors.
*
* Let D be the block diagonal matrix with the real eigenvalues in 1x1 blocks,
* and any complex values u+iv in 2x2 blocks [u v ; -v u]. Then, the matrices D
* and V satisfy A*V = V*D.
*
* More precisely, if the diagonal matrix of the eigen values is:\n
* \f$
* \left[ \begin{array}{cccccc}
* u+iv & & & & & \\
* & u-iv & & & & \\
* & & a+ib & & & \\
* & & & a-ib & & \\
* & & & & x & \\
* & & & & & y \\
* \end{array} \right]
* \f$ \n
* then, we have:\n
* \f$
* D =\left[ \begin{array}{cccccc}
* u & v & & & & \\
* -v & u & & & & \\
* & & a & b & & \\
* & & -b & a & & \\
* & & & & x & \\
* & & & & & y \\
* \end{array} \right]
* \f$
*
* \sa pseudoEigenvalueMatrix()
*/
const MatrixType& pseudoEigenvectors() const { return m_eivec; }
/** \returns the real block diagonal matrix D of the eigenvalues.
*
* See pseudoEigenvectors() for the details.
*/
MatrixType pseudoEigenvalueMatrix() const;
/** \returns the eigenvalues as a column vector */
EigenvalueType eigenvalues() const { return m_eivalues; }
void compute(const MatrixType& matrix);
@ -74,6 +115,25 @@ template<typename _MatrixType> class EigenSolver
EigenvalueType m_eivalues;
};
template<typename MatrixType>
MatrixType EigenSolver<MatrixType>::pseudoEigenvalueMatrix() const
{
int n = m_eivec.cols();
MatrixType matD = MatrixType::Zero(n,n);
for (int i=0; i<n; i++)
{
if (ei_isMuchSmallerThan(ei_imag(m_eivalues.coeff(i)), ei_real(m_eivalues.coeff(i))))
matD.coeffRef(i,i) = ei_real(m_eivalues.coeff(i));
else
{
matD.template block<2,2>(i,i) << ei_real(m_eivalues.coeff(i)), ei_imag(m_eivalues.coeff(i)),
-ei_imag(m_eivalues.coeff(i)), ei_real(m_eivalues.coeff(i));
i++;
}
}
return matD;
}
template<typename MatrixType>
void EigenSolver<MatrixType>::compute(const MatrixType& matrix)
{

55
test/eigensolver.cpp
View File

@ -29,7 +29,7 @@
#include "gsl_helper.h"
#endif
template<typename MatrixType> void eigensolver(const MatrixType& m)
template<typename MatrixType> void selfadjointeigensolver(const MatrixType& m)
{
/* this test covers the following files:
EigenSolver.h, SelfAdjointEigenSolver.h (and indirectly: Tridiagonalization.h)
@ -69,11 +69,11 @@ template<typename MatrixType> void eigensolver(const MatrixType& m)
convert<MatrixType>(symmB, gSymmB);
convert<MatrixType>(symmA, gEvec);
gEval = GslTraits<RealScalar>::createVector(rows);
Gsl::eigen_symm(gSymmA, gEval, gEvec);
convert(gEval, _eval);
convert(gEvec, _evec);
// test gsl itself !
VERIFY((symmA * _evec).isApprox(_evec * _eval.asDiagonal().eval(), largerEps));
@ -108,13 +108,40 @@ template<typename MatrixType> void eigensolver(const MatrixType& m)
VERIFY((symmA * eiSymmGen.eigenvectors()).isApprox(
symmB * (eiSymmGen.eigenvectors() * eiSymmGen.eigenvalues().asDiagonal().eval()), largerEps));
// EigenSolver<MatrixType> eiNotSymmButSymm(covMat);
// VERIFY_IS_APPROX((covMat.template cast<Complex>()) * (eiNotSymmButSymm.eigenvectors().template cast<Complex>()),
// (eiNotSymmButSymm.eigenvectors().template cast<Complex>()) * (eiNotSymmButSymm.eigenvalues().asDiagonal()));
}
// EigenSolver<MatrixType> eiNotSymm(a);
// VERIFY_IS_APPROX(a.template cast<Complex>() * eiNotSymm.eigenvectors().template cast<Complex>(),
// eiNotSymm.eigenvectors().template cast<Complex>() * eiNotSymm.eigenvalues().asDiagonal());
template<typename MatrixType> void eigensolver(const MatrixType& m)
{
/* this test covers the following files:
EigenSolver.h
*/
int rows = m.rows();
int cols = m.cols();
typedef typename MatrixType::Scalar Scalar;
typedef typename NumTraits<Scalar>::Real RealScalar;
typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType;
typedef Matrix<RealScalar, MatrixType::RowsAtCompileTime, 1> RealVectorType;
typedef typename std::complex<typename NumTraits<typename MatrixType::Scalar>::Real> Complex;
RealScalar largerEps = 10*test_precision<RealScalar>();
MatrixType a = MatrixType::Random(rows,cols);
MatrixType a1 = MatrixType::Random(rows,cols);
MatrixType symmA = a.adjoint() * a + a1.adjoint() * a1;
// MatrixType b = MatrixType::Random(rows,cols);
// MatrixType b1 = MatrixType::Random(rows,cols);
// MatrixType symmB = b.adjoint() * b + b1.adjoint() * b1;
EigenSolver<MatrixType> ei0(symmA);
VERIFY_IS_APPROX((symmA.template cast<Complex>()) * (ei0.pseudoEigenvectors().template cast<Complex>()),
(ei0.pseudoEigenvectors().template cast<Complex>()) * (ei0.eigenvalues().asDiagonal()));
a = a + symmA;
EigenSolver<MatrixType> ei1(a);
VERIFY_IS_APPROX(a * ei1.pseudoEigenvectors(), ei1.pseudoEigenvectors() * ei1.pseudoEigenvalueMatrix());
}
@ -122,10 +149,12 @@ void test_eigensolver()
{
for(int i = 0; i < g_repeat; i++) {
// very important to test a 3x3 matrix since we provide a special path for it
CALL_SUBTEST( eigensolver(Matrix3f()) );
CALL_SUBTEST( selfadjointeigensolver(Matrix3f()) );
CALL_SUBTEST( selfadjointeigensolver(Matrix4d()) );
CALL_SUBTEST( selfadjointeigensolver(MatrixXf(7,7)) );
CALL_SUBTEST( selfadjointeigensolver(MatrixXcd(5,5)) );
CALL_SUBTEST( selfadjointeigensolver(MatrixXd(19,19)) );
CALL_SUBTEST( eigensolver(Matrix4d()) );
CALL_SUBTEST( eigensolver(MatrixXf(7,7)) );
CALL_SUBTEST( eigensolver(MatrixXcd(5,5)) );
CALL_SUBTEST( eigensolver(MatrixXd(19,19)) );
}
}