diff --git a/Eigen/QR b/Eigen/QR index 71a72b65c..b893e104e 100644 --- a/Eigen/QR +++ b/Eigen/QR @@ -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 +// +// 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 . + +#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 class HessenbergDecomposition +{ + public: + + typedef _MatrixType MatrixType; + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits::Real RealScalar; + + enum { + Size = MatrixType::RowsAtCompileTime, + SizeMinusOne = MatrixType::RowsAtCompileTime==Dynamic + ? Dynamic + : MatrixType::RowsAtCompileTime-1}; + + typedef Matrix CoeffVectorType; + typedef Matrix DiagonalType; + typedef Matrix SubDiagonalType; + + typedef typename NestByValue >::RealReturnType DiagonalReturnType; + + typedef typename NestByValue > > >::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 +void HessenbergDecomposition::_compute(MatrixType& matA, CoeffVectorType& hCoeffs) +{ + assert(matA.rows()==matA.cols()); + int n = matA.rows(); + for (int i = 0; i(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::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 HessenbergDecomposition::MatrixType +HessenbergDecomposition::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 HessenbergDecomposition::MatrixType +HessenbergDecomposition::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().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 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 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::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 Tridiagonalization::MatrixType +Tridiagonalization::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().setZero(); + matT.corner(BottomLeft,n-2, n-2).template part().setZero(); + return matT; +} + /** Performs a full decomposition in place */ template void Tridiagonalization::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 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 . +// 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 void qr(const MatrixType& m) MatrixType a = MatrixType::random(rows,cols); QR 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 tridiag(b); + VERIFY_IS_APPROX(b, tridiag.matrixQ() * tridiag.matrixT() * tridiag.matrixQ().adjoint()); + + // check hessenberg decomposition + HessenbergDecomposition 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)) ); } }