diff options
| -rw-r--r-- | Eigen/QR | 1 | ||||
| -rwxr-xr-x | Eigen/src/QR/HessenbergDecomposition.h | 243 | ||||
| -rwxr-xr-x | Eigen/src/QR/Tridiagonalization.h | 24 | ||||
| -rw-r--r-- | test/eigensolver.cpp | 2 | ||||
| -rw-r--r-- | test/nomalloc.cpp | 2 | ||||
| -rw-r--r-- | test/qr.cpp | 20 |
6 files changed, 287 insertions, 5 deletions
@@ -9,6 +9,7 @@ namespace Eigen { #include "src/QR/Tridiagonalization.h" #include "src/QR/EigenSolver.h" #include "src/QR/SelfAdjointEigenSolver.h" +#include "src/QR/HessenbergDecomposition.h" } // namespace Eigen diff --git a/Eigen/src/QR/HessenbergDecomposition.h b/Eigen/src/QR/HessenbergDecomposition.h new file mode 100755 index 000000000..0cfd61832 --- /dev/null +++ b/Eigen/src/QR/HessenbergDecomposition.h @@ -0,0 +1,243 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. Eigen itself is part of the KDE project. +// +// Copyright (C) 2008 Gael Guennebaud <g.gael@free.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_HESSENBERGDECOMPOSITION_H +#define EIGEN_HESSENBERGDECOMPOSITION_H + +/** \class HessenbergDecomposition + * + * \brief Reduces a squared matrix to an Hessemberg form + * + * \param MatrixType the type of the matrix of which we are computing the Hessenberg decomposition + * + * This class performs an Hessenberg decomposition of a matrix \f$ A \f$ such that: + * \f$ A = Q H Q^* \f$ where \f$ Q \f$ is unitary and \f$ H \f$ a Hessenberg matrix. + * + * \sa class Tridiagonalization, class Qr + */ +template<typename _MatrixType> class HessenbergDecomposition +{ + public: + + typedef _MatrixType MatrixType; + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits<Scalar>::Real RealScalar; + + enum { + Size = MatrixType::RowsAtCompileTime, + SizeMinusOne = MatrixType::RowsAtCompileTime==Dynamic + ? Dynamic + : MatrixType::RowsAtCompileTime-1}; + + typedef Matrix<Scalar, SizeMinusOne, 1> CoeffVectorType; + typedef Matrix<RealScalar, Size, 1> DiagonalType; + typedef Matrix<RealScalar, SizeMinusOne, 1> SubDiagonalType; + + typedef typename NestByValue<DiagonalCoeffs<MatrixType> >::RealReturnType DiagonalReturnType; + + typedef typename NestByValue<DiagonalCoeffs< + NestByValue<Block< + MatrixType,SizeMinusOne,SizeMinusOne> > > >::RealReturnType SubDiagonalReturnType; + + HessenbergDecomposition() + {} + + HessenbergDecomposition(int rows, int cols) + : m_matrix(rows,cols), m_hCoeffs(rows-1) + {} + + HessenbergDecomposition(const MatrixType& matrix) + : m_matrix(matrix), + m_hCoeffs(matrix.cols()-1) + { + _compute(m_matrix, m_hCoeffs); + } + + /** Computes or re-compute the Hessenberg decomposition for the matrix \a matrix. + * + * This method allows to re-use the allocated data. + */ + void compute(const MatrixType& matrix) + { + m_matrix = matrix; + m_hCoeffs.resize(matrix.rows()-1,1); + _compute(m_matrix, m_hCoeffs); + } + + /** \returns the householder coefficients allowing to + * reconstruct the matrix Q from the packed data. + * + * \sa packedMatrix() + */ + CoeffVectorType householderCoefficients(void) const { return m_hCoeffs; } + + /** \returns the internal result of the decomposition. + * + * The returned matrix contains the following information: + * - the upper part and lower sub-diagonal represent the Hessenberg matrix H + * - the rest of the lower part contains the Householder vectors that, combined with + * Householder coefficients returned by householderCoefficients(), + * allows to reconstruct the matrix Q as follow: + * Q = H_{N-1} ... H_1 H_0 + * where the matrices H are the Householder transformation: + * H_i = (I - h_i * v_i * v_i') + * where h_i == householderCoefficients()[i] and v_i is a Householder vector: + * v_i = [ 0, ..., 0, 1, M(i+2,i), ..., M(N-1,i) ] + * + * See LAPACK for further details on this packed storage. + */ + const MatrixType& packedMatrix(void) const { return m_matrix; } + + MatrixType matrixQ(void) const; + MatrixType matrixH(void) const; + + private: + + static void _compute(MatrixType& matA, CoeffVectorType& hCoeffs); + + protected: + MatrixType m_matrix; + CoeffVectorType m_hCoeffs; +}; + + +/** \internal + * Performs a tridiagonal decomposition of \a matA in place. + * + * \param matA the input selfadjoint matrix + * \param hCoeffs returned Householder coefficients + * + * The result is written in the lower triangular part of \a matA. + * + * Implemented from Golub's "Matrix Computations", algorithm 8.3.1. + * + * \sa packedMatrix() + */ +template<typename MatrixType> +void HessenbergDecomposition<MatrixType>::_compute(MatrixType& matA, CoeffVectorType& hCoeffs) +{ + assert(matA.rows()==matA.cols()); + int n = matA.rows(); + for (int i = 0; i<n-2; ++i) + { + // let's consider the vector v = i-th column starting at position i+1 + + // start of the householder transformation + // squared norm of the vector v skipping the first element + RealScalar v1norm2 = matA.col(i).end(n-(i+2)).norm2(); + + if (ei_isMuchSmallerThan(v1norm2,static_cast<Scalar>(1))) + { + hCoeffs.coeffRef(i) = 0.; + } + else + { + Scalar v0 = matA.col(i).coeff(i+1); + RealScalar beta = ei_sqrt(ei_abs2(v0)+v1norm2); + if (ei_real(v0)>=0.) + beta = -beta; + matA.col(i).end(n-(i+2)) *= (Scalar(1)/(v0-beta)); + matA.col(i).coeffRef(i+1) = beta; + Scalar h = (beta - v0) / beta; + // end of the householder transformation + + // Apply similarity transformation to remaining columns, + // i.e., A = H' A H where H = I - h v v' and v = matA.col(i).end(n-i-1) + matA.col(i).coeffRef(i+1) = 1; + + // first let's do A = H A + matA.corner(BottomRight,n-i-1,n-i-1) -= ((ei_conj(h) * matA.col(i).end(n-i-1)) * + (matA.col(i).end(n-i-1).adjoint() * matA.corner(BottomRight,n-i-1,n-i-1)).lazy()).lazy(); + + // now let's do A = A H + matA.corner(BottomRight,n,n-i-1) -= ((matA.corner(BottomRight,n,n-i-1) * matA.col(i).end(n-i-1)).lazy() * + (h * matA.col(i).end(n-i-1).adjoint())).lazy(); + + matA.col(i).coeffRef(i+1) = beta; + hCoeffs.coeffRef(i) = h; + } + } + if (NumTraits<Scalar>::IsComplex) + { + // Householder transformation on the remaining single scalar + int i = n-2; + Scalar v0 = matA.coeff(i+1,i); + + RealScalar beta = ei_sqrt(ei_abs2(v0)); + if (ei_real(v0)>=0.) + beta = -beta; + Scalar h = (beta - v0) / beta; + hCoeffs.coeffRef(i) = h; + + // A = H* A + matA.corner(BottomRight,n-i-1,n-i) -= ei_conj(h) * matA.corner(BottomRight,n-i-1,n-i); + + // A = A H + matA.col(n-1) -= h * matA.col(n-1); + } + else + { + hCoeffs.coeffRef(n-2) = 0; + } +} + +/** reconstructs and returns the matrix Q */ +template<typename MatrixType> +typename HessenbergDecomposition<MatrixType>::MatrixType +HessenbergDecomposition<MatrixType>::matrixQ(void) const +{ + int n = m_matrix.rows(); + MatrixType matQ = MatrixType::identity(n,n); + for (int i = n-2; i>=0; i--) + { + Scalar tmp = m_matrix.coeff(i+1,i); + m_matrix.const_cast_derived().coeffRef(i+1,i) = 1; + + matQ.corner(BottomRight,n-i-1,n-i-1) -= + ((m_hCoeffs.coeff(i) * m_matrix.col(i).end(n-i-1)) * + (m_matrix.col(i).end(n-i-1).adjoint() * matQ.corner(BottomRight,n-i-1,n-i-1)).lazy()).lazy(); + + m_matrix.const_cast_derived().coeffRef(i+1,i) = tmp; + } + return matQ; +} + +/** constructs and returns the matrix H. + * Note that the matrix H is equivalent to the upper part of the packed matrix + * (including the lower sub-diagonal). Therefore, it might be often sufficient + * to directly use the packed matrix instead of creating a new one. + */ +template<typename MatrixType> +typename HessenbergDecomposition<MatrixType>::MatrixType +HessenbergDecomposition<MatrixType>::matrixH(void) const +{ + // FIXME should this function (and other similar) rather take a matrix as argument + // and fill it (avoids temporaries) + int n = m_matrix.rows(); + MatrixType matH = m_matrix; + matH.corner(BottomLeft,n-2, n-2).template part<Lower>().setZero(); + return matH; +} + +#endif // EIGEN_HESSENBERGDECOMPOSITION_H diff --git a/Eigen/src/QR/Tridiagonalization.h b/Eigen/src/QR/Tridiagonalization.h index e5b412c46..e76fbad96 100755 --- a/Eigen/src/QR/Tridiagonalization.h +++ b/Eigen/src/QR/Tridiagonalization.h @@ -29,7 +29,7 @@ * * \brief Trigiagonal decomposition of a selfadjoint matrix * - * \param MatrixType the type of the matrix of which we are computing the eigen decomposition + * \param MatrixType the type of the matrix of which we are performing the tridiagonalization * * This class performs a tridiagonal decomposition of a selfadjoint matrix \f$ A \f$ such that: * \f$ A = Q T Q^* \f$ where \f$ Q \f$ is unitatry and \f$ T \f$ a real symmetric tridiagonal matrix @@ -81,7 +81,7 @@ template<typename _MatrixType> class Tridiagonalization void compute(const MatrixType& matrix) { m_matrix = matrix; - m_hCoeffs.resize(matrix.rows()-1); + m_hCoeffs.resize(matrix.rows()-1, 1); _compute(m_matrix, m_hCoeffs); } @@ -111,6 +111,7 @@ template<typename _MatrixType> class Tridiagonalization const MatrixType& packedMatrix(void) const { return m_matrix; } MatrixType matrixQ(void) const; + MatrixType matrixT(void) const; const DiagonalReturnType diagonal(void) const; const SubDiagonalReturnType subDiagonal(void) const; @@ -252,6 +253,25 @@ Tridiagonalization<MatrixType>::subDiagonal(void) const .nestByValue().diagonal().nestByValue().real(); } +/** constructs and returns the tridiagonal matrix T. + * Note that the matrix T is equivalent to the diagonal and sub-diagonal of the packed matrix. + * Therefore, it might be often sufficient to directly use the packed matrix, or the vector + * expressions returned by diagonal() and subDiagonal() instead of creating a new matrix. + */ +template<typename MatrixType> +typename Tridiagonalization<MatrixType>::MatrixType +Tridiagonalization<MatrixType>::matrixT(void) const +{ + // FIXME should this function (and other similar) rather take a matrix as argument + // and fill it (avoids temporaries) + int n = m_matrix.rows(); + MatrixType matT = m_matrix; + matT.corner(TopRight,n-1, n-1).diagonal() = subDiagonal().conjugate(); + matT.corner(TopRight,n-2, n-2).template part<Upper>().setZero(); + matT.corner(BottomLeft,n-2, n-2).template part<Lower>().setZero(); + return matT; +} + /** Performs a full decomposition in place */ template<typename MatrixType> void Tridiagonalization<MatrixType>::decomposeInPlace(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ) diff --git a/test/eigensolver.cpp b/test/eigensolver.cpp index 820446ac9..cb519657c 100644 --- a/test/eigensolver.cpp +++ b/test/eigensolver.cpp @@ -54,9 +54,11 @@ template<typename MatrixType> void eigensolver(const MatrixType& m) void test_eigensolver() { for(int i = 0; i < 1; i++) { + // very important to test a 3x3 matrix since we provide a special path for it CALL_SUBTEST( eigensolver(Matrix3f()) ); CALL_SUBTEST( eigensolver(Matrix4d()) ); CALL_SUBTEST( eigensolver(MatrixXd(7,7)) ); CALL_SUBTEST( eigensolver(MatrixXcd(6,6)) ); + CALL_SUBTEST( eigensolver(MatrixXcd(3,3)) ); } } diff --git a/test/nomalloc.cpp b/test/nomalloc.cpp index 6d40222a4..d6c566bef 100644 --- a/test/nomalloc.cpp +++ b/test/nomalloc.cpp @@ -22,7 +22,9 @@ // License and a copy of the GNU General Public License along with // Eigen. If not, see <http://www.gnu.org/licenses/>. +// this hack is needed to make this file compiles with -pedantic (gcc) #define throw(X) +// discard vectorization since operator new is not called in that case #define EIGEN_DONT_VECTORIZE 1 #include "main.h" diff --git a/test/qr.cpp b/test/qr.cpp index 8facbac96..716ec94c9 100644 --- a/test/qr.cpp +++ b/test/qr.cpp @@ -39,17 +39,31 @@ template<typename MatrixType> void qr(const MatrixType& m) MatrixType a = MatrixType::random(rows,cols); QR<MatrixType> qrOfA(a); - VERIFY_IS_APPROX(a, qrOfA.matrixQ() * qrOfA.matrixR()); VERIFY_IS_NOT_APPROX(a+MatrixType::identity(rows, cols), qrOfA.matrixQ() * qrOfA.matrixR()); + + SquareMatrixType b = a.adjoint() * a; + + // check tridiagonalization + Tridiagonalization<SquareMatrixType> tridiag(b); + VERIFY_IS_APPROX(b, tridiag.matrixQ() * tridiag.matrixT() * tridiag.matrixQ().adjoint()); + + // check hessenberg decomposition + HessenbergDecomposition<SquareMatrixType> hess(b); + VERIFY_IS_APPROX(b, hess.matrixQ() * hess.matrixH() * hess.matrixQ().adjoint()); + VERIFY_IS_APPROX(tridiag.matrixT(), hess.matrixH()); + b = SquareMatrixType::random(cols,cols); + hess.compute(b); + VERIFY_IS_APPROX(b, hess.matrixQ() * hess.matrixH() * hess.matrixQ().adjoint()); } void test_qr() { for(int i = 0; i < 1; i++) { CALL_SUBTEST( qr(Matrix2f()) ); - CALL_SUBTEST( qr(Matrix3d()) ); + CALL_SUBTEST( qr(Matrix4d()) ); CALL_SUBTEST( qr(MatrixXf(12,8)) ); -// CALL_SUBTEST( qr(MatrixXcd(17,7)) ); // complex numbers are not supported yet + CALL_SUBTEST( qr(MatrixXcd(5,5)) ); + CALL_SUBTEST( qr(MatrixXcd(7,3)) ); } } |