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

bug #645: patch from Tobias Wood implementing the extraction of eigenvectors in GeneralizedEigenSolver

Browse Source
This commit is contained in:
Gael Guennebaud 2016-08-23 17:37:38 +02:00
parent 504a4404f1
commit 00b2666853
2 changed files with 125 additions and 66 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

182
Eigen/src/Eigenvalues/GeneralizedEigenSolver.h
View File

@ -1,8 +1,9 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2012 Gael Guennebaud <gael.guennebaud@inria.fr>
// Copyright (C) 2012-2016 Gael Guennebaud <gael.guennebaud@inria.fr>
// Copyright (C) 2010,2012 Jitse Niesen <jitse@maths.leeds.ac.uk>
// Copyright (C) 2016 Tobias Wood <tobias@spinicist.org.uk>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
@ -89,7 +90,7 @@ template<typename _MatrixType> class GeneralizedEigenSolver
*/
typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> VectorType;
/** \brief Type for vector of complex scalar values eigenvalues as returned by betas().
/** \brief Type for vector of complex scalar values eigenvalues as returned by alphas().
*
* This is a column vector with entries of type #ComplexScalar.
* The length of the vector is the size of #MatrixType.
@ -114,7 +115,14 @@ template<typename _MatrixType> class GeneralizedEigenSolver
*
* \sa compute() for an example.
*/
GeneralizedEigenSolver() : m_eivec(), m_alphas(), m_betas(), m_isInitialized(false), m_realQZ(), m_matS(), m_tmp() {}
GeneralizedEigenSolver()
: m_eivec(),
m_alphas(),
m_betas(),
m_valuesOkay(false),
m_vectorsOkay(false),
m_realQZ()
{}
/** \brief Default constructor with memory preallocation
*
@ -126,11 +134,9 @@ template<typename _MatrixType> class GeneralizedEigenSolver
: m_eivec(size, size),
m_alphas(size),
m_betas(size),
m_isInitialized(false),
m_eigenvectorsOk(false),
m_realQZ(size),
m_matS(size, size),
m_tmp(size)
m_valuesOkay(false),
m_vectorsOkay(false),
m_realQZ(size)
{}
/** \brief Constructor; computes the generalized eigendecomposition of given matrix pair.
@ -149,33 +155,29 @@ template<typename _MatrixType> class GeneralizedEigenSolver
: m_eivec(A.rows(), A.cols()),
m_alphas(A.cols()),
m_betas(A.cols()),
m_isInitialized(false),
m_eigenvectorsOk(false),
m_realQZ(A.cols()),
m_matS(A.rows(), A.cols()),
m_tmp(A.cols())
m_valuesOkay(false),
m_vectorsOkay(false),
m_realQZ(A.cols())
{
compute(A, B, computeEigenvectors);
}
/* \brief Returns the computed generalized eigenvectors.
*
* \returns %Matrix whose columns are the (possibly complex) eigenvectors.
* \returns %Matrix whose columns are the (possibly complex) right eigenvectors.
* i.e. the eigenvectors that solve (A - l*B)x = 0. The ordering matches the eigenvalues.
*
* \pre Either the constructor
* GeneralizedEigenSolver(const MatrixType&,const MatrixType&, bool) or the member function
* compute(const MatrixType&, 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
* eigenvectors are normalized to have (Euclidean) norm equal to one. The
* matrix returned by this function is the matrix \f$ V \f$ in the
* generalized eigendecomposition \f$ A = B V D V^{-1} \f$, if it exists.
*
* \sa eigenvalues()
*/
// EigenvectorsType eigenvectors() const;
EigenvectorsType eigenvectors() const {
eigen_assert(m_vectorsOkay && "Eigenvectors for GeneralizedEigenSolver were not calculated.")
return m_eivec;
}
/** \brief Returns an expression of the computed generalized eigenvalues.
*
@ -197,7 +199,7 @@ template<typename _MatrixType> class GeneralizedEigenSolver
*/
EigenvalueType eigenvalues() const
{
eigen_assert(m_isInitialized && "GeneralizedEigenSolver is not initialized.");
eigen_assert(m_valuesOkay && "GeneralizedEigenSolver is not initialized.");
return EigenvalueType(m_alphas,m_betas);
}
@ -208,7 +210,7 @@ template<typename _MatrixType> class GeneralizedEigenSolver
* \sa betas(), eigenvalues() */
ComplexVectorType alphas() const
{
eigen_assert(m_isInitialized && "GeneralizedEigenSolver is not initialized.");
eigen_assert(m_valuesOkay && "GeneralizedEigenSolver is not initialized.");
return m_alphas;
}
@ -219,7 +221,7 @@ template<typename _MatrixType> class GeneralizedEigenSolver
* \sa alphas(), eigenvalues() */
VectorType betas() const
{
eigen_assert(m_isInitialized && "GeneralizedEigenSolver is not initialized.");
eigen_assert(m_valuesOkay && "GeneralizedEigenSolver is not initialized.");
return m_betas;
}
@ -250,7 +252,7 @@ template<typename _MatrixType> class GeneralizedEigenSolver
ComputationInfo info() const
{
eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
eigen_assert(m_valuesOkay && "EigenSolver is not initialized.");
return m_realQZ.info();
}
@ -270,29 +272,13 @@ template<typename _MatrixType> class GeneralizedEigenSolver
EIGEN_STATIC_ASSERT(!NumTraits<Scalar>::IsComplex, NUMERIC_TYPE_MUST_BE_REAL);
}
MatrixType m_eivec;
EigenvectorsType m_eivec;
ComplexVectorType m_alphas;
VectorType m_betas;
bool m_isInitialized;
bool m_eigenvectorsOk;
bool m_valuesOkay, m_vectorsOkay;
RealQZ<MatrixType> m_realQZ;
MatrixType m_matS;
typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ColumnVectorType;
ColumnVectorType m_tmp;
};
//template<typename MatrixType>
//typename GeneralizedEigenSolver<MatrixType>::EigenvectorsType GeneralizedEigenSolver<MatrixType>::eigenvectors() const
//{
// eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
// eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
// Index n = m_eivec.cols();
// EigenvectorsType matV(n,n);
// // TODO
// return matV;
//}
template<typename MatrixType>
GeneralizedEigenSolver<MatrixType>&
GeneralizedEigenSolver<MatrixType>::compute(const MatrixType& A, const MatrixType& B, bool computeEigenvectors)
@ -302,28 +288,70 @@ GeneralizedEigenSolver<MatrixType>::compute(const MatrixType& A, const MatrixTyp
using std::sqrt;
using std::abs;
eigen_assert(A.cols() == A.rows() && B.cols() == A.rows() && B.cols() == B.rows());
m_valuesOkay = false;
m_vectorsOkay = false;
// Reduce to generalized real Schur form:
// A = Q S Z and B = Q T Z
m_realQZ.compute(A, B, computeEigenvectors);
if (m_realQZ.info() == Success)
{
m_matS = m_realQZ.matrixS();
const MatrixType &matT = m_realQZ.matrixT();
if (computeEigenvectors)
m_eivec = m_realQZ.matrixZ().transpose();
// Compute eigenvalues from matS
// Temp space for the untransformed eigenvectors
VectorType v;
ComplexVectorType cv;
// Resize storage
m_alphas.resize(A.cols());
m_betas.resize(A.cols());
if (computeEigenvectors) {
m_eivec.resize(A.cols(), A.cols());
v.resize(A.cols());
cv.resize(A.cols());
}
// Grab some references
const MatrixType &mZT = m_realQZ.matrixZ().transpose();
const MatrixType &mS = m_realQZ.matrixS();
const MatrixType &mT = m_realQZ.matrixT();
Index i = 0;
while (i < A.cols())
{
if (i == A.cols() - 1 || m_matS.coeff(i+1, i) == Scalar(0))
if (i == A.cols() - 1 || mS.coeff(i+1, i) == Scalar(0))
{
m_alphas.coeffRef(i) = m_matS.coeff(i, i);
m_betas.coeffRef(i) = matT.coeff(i,i);
// Real eigenvalue
m_alphas.coeffRef(i) = mS.diagonal().coeff(i);
m_betas.coeffRef(i) = mT.diagonal().coeff(i);
if (computeEigenvectors)
{
v.setConstant(Scalar(0.0));
v.coeffRef(i) = Scalar(1.0);
// For singular eigenvalues do nothing more
if(abs(m_betas.coeffRef(i)) >= (std::numeric_limits<RealScalar>::min)())
{
// Non-singular eigenvalue
const Scalar alpha = real(m_alphas.coeffRef(i));
const Scalar beta = m_betas.coeffRef(i);
for (Index j = i-1; j >= 0; j--)
{
const Index st = j+1;
const Index sz = i-j;
if (j > 0 && mS.coeff(j, j-1) != Scalar(0))
{ // 2x2 block
Matrix<Scalar, 2, 1> RHS;
RHS[0] = -v.segment(st,sz).transpose().cwiseProduct(beta*mS.block(j-1,st,1,sz) - alpha*mT.block(j-1,st,1,sz)).sum();
RHS[1] = -v.segment(st,sz).transpose().cwiseProduct(beta*mS.block(j,st,1,sz) - alpha*mT.block(j,st,1,sz)).sum();
Matrix<Scalar, 2, 2> LHS;
LHS << beta*mS.coeffRef(j-1,j-1) - alpha*mT.coeffRef(j-1,j-1), beta*mS.coeffRef(j-1,j) - alpha*mT.coeffRef(j-1,j),
beta*mS.coeffRef(j,j-1) - alpha*mT.coeffRef(j,j-1), beta*mS.coeffRef(j,j) - alpha*mT.coeffRef(j,j);
v.segment(j-1,2) = LHS.partialPivLu().solve(RHS);
j--;
}
else
{
v.coeffRef(j) = -v.segment(st,sz).transpose().cwiseProduct(beta*mS.block(j,st,1,sz) - alpha*mT.block(j,st,1,sz)).sum() / (beta*mS.coeffRef(j,j) - alpha*mT.coeffRef(j,j));
}
}
}
m_eivec.col(i).real() = (mZT * v).normalized();
m_eivec.col(i).imag().setConstant(0);
}
++i;
}
else
@ -333,27 +361,49 @@ GeneralizedEigenSolver<MatrixType>::compute(const MatrixType& A, const MatrixTyp
// T = [a 0]
// [0 b]
RealScalar a = matT.diagonal().coeff(i),
b = matT.diagonal().coeff(i+1);
RealScalar a = mT.diagonal().coeff(i),
b = mT.diagonal().coeff(i+1);
const RealScalar beta = m_betas.coeffRef(i) = m_betas.coeffRef(i+1) = a*b;
// ^^ NOTE: using diagonal()(i) instead of coeff(i,i) workarounds a MSVC bug.
Matrix<RealScalar,2,2> S2 = m_matS.template block<2,2>(i,i) * Matrix<Scalar,2,1>(b,a).asDiagonal();
Matrix<RealScalar,2,2> S2 = mS.template block<2,2>(i,i) * Matrix<Scalar,2,1>(b,a).asDiagonal();
Scalar p = Scalar(0.5) * (S2.coeff(0,0) - S2.coeff(1,1));
Scalar z = sqrt(abs(p * p + S2.coeff(1,0) * S2.coeff(0,1)));
m_alphas.coeffRef(i) = ComplexScalar(S2.coeff(1,1) + p, z);
m_alphas.coeffRef(i+1) = ComplexScalar(S2.coeff(1,1) + p, -z);
const ComplexScalar alpha = ComplexScalar(S2.coeff(1,1) + p, (beta > 0) ? z : -z);
m_alphas.coeffRef(i) = conj(alpha);
m_alphas.coeffRef(i+1) = alpha;
m_betas.coeffRef(i) =
m_betas.coeffRef(i+1) = a*b;
if (computeEigenvectors) {
// Compute eigenvector in position (i+1) and then position (i) is just the conjugate
cv.setConstant(ComplexScalar(0.0, 0.0));
cv.coeffRef(i+1) = ComplexScalar(1.0, 0.0);
cv.coeffRef(i) = -(beta*mS.coeffRef(i,i+1) - alpha*mT.coeffRef(i,i+1)) / (beta*mS.coeffRef(i,i) - alpha*mT.coeffRef(i,i));
for (Index j = i-1; j >= 0; j--) {
const Index st = j+1;
const Index sz = i+1-j;
if (j > 0 && mS.coeff(j, j-1) != Scalar(0)) { // 2x2 block
Matrix<ComplexScalar, 2, 1> RHS;
RHS[0] = -cv.segment(st,sz).transpose().cwiseProduct(beta*mS.block(j-1,st,1,sz) - alpha*mT.block(j-1,st,1,sz)).sum();
RHS[1] = -cv.segment(st,sz).transpose().cwiseProduct(beta*mS.block(j,st,1,sz) - alpha*mT.block(j,st,1,sz)).sum();
Matrix<ComplexScalar, 2, 2> LHS;
LHS << beta*mS.coeffRef(j-1,j-1) - alpha*mT.coeffRef(j-1,j-1), beta*mS.coeffRef(j-1,j) - alpha*mT.coeffRef(j-1,j),
beta*mS.coeffRef(j,j-1) - alpha*mT.coeffRef(j,j-1), beta*mS.coeffRef(j,j) - alpha*mT.coeffRef(j,j);
cv.segment(j-1,2) = LHS.partialPivLu().solve(RHS);
j--;
} else {
cv.coeffRef(j) = -cv.segment(st,sz).transpose().cwiseProduct(beta*mS.block(j,st,1,sz) - alpha*mT.block(j,st,1,sz)).sum() / (beta*mS.coeffRef(j,j) - alpha*mT.coeffRef(j,j));
}
}
m_eivec.col(i+1) = (mZT * cv).normalized();
m_eivec.col(i) = m_eivec.col(i+1).conjugate();
}
i += 2;
}
}
m_valuesOkay = true;
m_vectorsOkay = computeEigenvectors;
}
m_isInitialized = true;
m_eigenvectorsOk = false;//computeEigenvectors;
return *this;
}

9
test/eigensolver_generalized_real.cpp
View File

@ -42,6 +42,11 @@ template<typename MatrixType> void generalized_eigensolver_real(const MatrixType
VectorType realEigenvalues = eig.eigenvalues().real();
std::sort(realEigenvalues.data(), realEigenvalues.data()+realEigenvalues.size());
VERIFY_IS_APPROX(realEigenvalues, symmEig.eigenvalues());
// check eigenvectors
typename GeneralizedEigenSolver<MatrixType>::EigenvectorsType D = eig.eigenvalues().asDiagonal();
typename GeneralizedEigenSolver<MatrixType>::EigenvectorsType V = eig.eigenvectors();
VERIFY_IS_APPROX(spdA*V, spdB*V*D);
}
// non symmetric case:
@ -54,6 +59,10 @@ template<typename MatrixType> void generalized_eigensolver_real(const MatrixType
tmp /= tmp.norm();
VERIFY_IS_MUCH_SMALLER_THAN( std::abs(tmp.determinant()), Scalar(1) );
}
// check eigenvectors
typename GeneralizedEigenSolver<MatrixType>::EigenvectorsType D = eig.eigenvalues().asDiagonal();
typename GeneralizedEigenSolver<MatrixType>::EigenvectorsType V = eig.eigenvectors();
VERIFY_IS_APPROX(a*V, b*V*D);
}
// regression test for bug 1098