src/Eigenvalues/ComplexSchur.h

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2009 Claire Maurice
0005 // Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr>
0006 // Copyright (C) 2010,2012 Jitse Niesen <jitse@maths.leeds.ac.uk>
0007 //
0008 // This Source Code Form is subject to the terms of the Mozilla
0009 // Public License v. 2.0. If a copy of the MPL was not distributed
0010 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0011
0012 #ifndef EIGEN_COMPLEX_SCHUR_H
0013 #define EIGEN_COMPLEX_SCHUR_H
0014
0015 #include "./HessenbergDecomposition.h"
0016
0017 namespace Eigen {
0018
0019 namespace internal {
0020 template<typename MatrixType, bool IsComplex> struct complex_schur_reduce_to_hessenberg;
0021 }
0022
0023 /** \eigenvalues_module \ingroup Eigenvalues_Module
0024   *
0025   *
0026   * \class ComplexSchur
0027   *
0028   * \brief Performs a complex Schur decomposition of a real or complex square matrix
0029   *
0030   * \tparam _MatrixType the type of the matrix of which we are
0031   * computing the Schur decomposition; this is expected to be an
0032   * instantiation of the Matrix class template.
0033   *
0034   * Given a real or complex square matrix A, this class computes the
0035   * Schur decomposition: \f$ A = U T U^*\f$ where U is a unitary
0036   * complex matrix, and T is a complex upper triangular matrix.  The
0037   * diagonal of the matrix T corresponds to the eigenvalues of the
0038   * matrix A.
0039   *
0040   * Call the function compute() to compute the Schur decomposition of
0041   * a given matrix. Alternatively, you can use the
0042   * ComplexSchur(const MatrixType&, bool) constructor which computes
0043   * the Schur decomposition at construction time. Once the
0044   * decomposition is computed, you can use the matrixU() and matrixT()
0045   * functions to retrieve the matrices U and V in the decomposition.
0046   *
0047   * \note This code is inspired from Jampack
0048   *
0049   * \sa class RealSchur, class EigenSolver, class ComplexEigenSolver
0050   */
0051 template<typename _MatrixType> class ComplexSchur
0052 {
0053   public:
0054     typedef _MatrixType MatrixType;
0055     enum {
0056       RowsAtCompileTime = MatrixType::RowsAtCompileTime,
0057       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
0058       Options = MatrixType::Options,
0059       MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
0060       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
0061     };
0062
0063     /** \brief Scalar type for matrices of type \p _MatrixType. */
0064     typedef typename MatrixType::Scalar Scalar;
0065     typedef typename NumTraits<Scalar>::Real RealScalar;
0066     typedef Eigen::Index Index; ///< \deprecated since Eigen 3.3
0067
0068     /** \brief Complex scalar type for \p _MatrixType.
0069       *
0070       * This is \c std::complex<Scalar> if #Scalar is real (e.g.,
0071       * \c float or \c double) and just \c Scalar if #Scalar is
0072       * complex.
0073       */
0074     typedef std::complex<RealScalar> ComplexScalar;
0075
0076     /** \brief Type for the matrices in the Schur decomposition.
0077       *
0078       * This is a square matrix with entries of type #ComplexScalar.
0079       * The size is the same as the size of \p _MatrixType.
0080       */
0081     typedef Matrix<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime> ComplexMatrixType;
0082
0083     /** \brief Default constructor.
0084       *
0085       * \param [in] size  Positive integer, size of the matrix whose Schur decomposition will be computed.
0086       *
0087       * The default constructor is useful in cases in which the user
0088       * intends to perform decompositions via compute().  The \p size
0089       * parameter is only used as a hint. It is not an error to give a
0090       * wrong \p size, but it may impair performance.
0091       *
0092       * \sa compute() for an example.
0093       */
0094     explicit ComplexSchur(Index size = RowsAtCompileTime==Dynamic ? 1 : RowsAtCompileTime)
0095       : m_matT(size,size),
0096         m_matU(size,size),
0097         m_hess(size),
0098         m_isInitialized(false),
0099         m_matUisUptodate(false),
0100         m_maxIters(-1)
0101     {}
0102
0103     /** \brief Constructor; computes Schur decomposition of given matrix.
0104       *
0105       * \param[in]  matrix    Square matrix whose Schur decomposition is to be computed.
0106       * \param[in]  computeU  If true, both T and U are computed; if false, only T is computed.
0107       *
0108       * This constructor calls compute() to compute the Schur decomposition.
0109       *
0110       * \sa matrixT() and matrixU() for examples.
0111       */
0112     template<typename InputType>
0113     explicit ComplexSchur(const EigenBase<InputType>& matrix, bool computeU = true)
0114       : m_matT(matrix.rows(),matrix.cols()),
0115         m_matU(matrix.rows(),matrix.cols()),
0116         m_hess(matrix.rows()),
0117         m_isInitialized(false),
0118         m_matUisUptodate(false),
0119         m_maxIters(-1)
0120     {
0121       compute(matrix.derived(), computeU);
0122     }
0123
0124     /** \brief Returns the unitary matrix in the Schur decomposition.
0125       *
0126       * \returns A const reference to the matrix U.
0127       *
0128       * It is assumed that either the constructor
0129       * ComplexSchur(const MatrixType& matrix, bool computeU) or the
0130       * member function compute(const MatrixType& matrix, bool computeU)
0131       * has been called before to compute the Schur decomposition of a
0132       * matrix, and that \p computeU was set to true (the default
0133       * value).
0134       *
0135       * Example: \include ComplexSchur_matrixU.cpp
0136       * Output: \verbinclude ComplexSchur_matrixU.out
0137       */
0138     const ComplexMatrixType& matrixU() const
0139     {
0140       eigen_assert(m_isInitialized && "ComplexSchur is not initialized.");
0141       eigen_assert(m_matUisUptodate && "The matrix U has not been computed during the ComplexSchur decomposition.");
0142       return m_matU;
0143     }
0144
0145     /** \brief Returns the triangular matrix in the Schur decomposition.
0146       *
0147       * \returns A const reference to the matrix T.
0148       *
0149       * It is assumed that either the constructor
0150       * ComplexSchur(const MatrixType& matrix, bool computeU) or the
0151       * member function compute(const MatrixType& matrix, bool computeU)
0152       * has been called before to compute the Schur decomposition of a
0153       * matrix.
0154       *
0155       * Note that this function returns a plain square matrix. If you want to reference
0156       * only the upper triangular part, use:
0157       * \code schur.matrixT().triangularView<Upper>() \endcode
0158       *
0159       * Example: \include ComplexSchur_matrixT.cpp
0160       * Output: \verbinclude ComplexSchur_matrixT.out
0161       */
0162     const ComplexMatrixType& matrixT() const
0163     {
0164       eigen_assert(m_isInitialized && "ComplexSchur is not initialized.");
0165       return m_matT;
0166     }
0167
0168     /** \brief Computes Schur decomposition of given matrix.
0169       *
0170       * \param[in]  matrix  Square matrix whose Schur decomposition is to be computed.
0171       * \param[in]  computeU  If true, both T and U are computed; if false, only T is computed.
0172
0173       * \returns    Reference to \c *this
0174       *
0175       * The Schur decomposition is computed by first reducing the
0176       * matrix to Hessenberg form using the class
0177       * HessenbergDecomposition. The Hessenberg matrix is then reduced
0178       * to triangular form by performing QR iterations with a single
0179       * shift. The cost of computing the Schur decomposition depends
0180       * on the number of iterations; as a rough guide, it may be taken
0181       * on the number of iterations; as a rough guide, it may be taken
0182       * to be \f$25n^3\f$ complex flops, or \f$10n^3\f$ complex flops
0183       * if \a computeU is false.
0184       *
0185       * Example: \include ComplexSchur_compute.cpp
0186       * Output: \verbinclude ComplexSchur_compute.out
0187       *
0188       * \sa compute(const MatrixType&, bool, Index)
0189       */
0190     template<typename InputType>
0191     ComplexSchur& compute(const EigenBase<InputType>& matrix, bool computeU = true);
0192
0193     /** \brief Compute Schur decomposition from a given Hessenberg matrix
0194      *  \param[in] matrixH Matrix in Hessenberg form H
0195      *  \param[in] matrixQ orthogonal matrix Q that transform a matrix A to H : A = Q H Q^T
0196      *  \param computeU Computes the matriX U of the Schur vectors
0197      * \return Reference to \c *this
0198      *
0199      *  This routine assumes that the matrix is already reduced in Hessenberg form matrixH
0200      *  using either the class HessenbergDecomposition or another mean.
0201      *  It computes the upper quasi-triangular matrix T of the Schur decomposition of H
0202      *  When computeU is true, this routine computes the matrix U such that
0203      *  A = U T U^T =  (QZ) T (QZ)^T = Q H Q^T where A is the initial matrix
0204      *
0205      * NOTE Q is referenced if computeU is true; so, if the initial orthogonal matrix
0206      * is not available, the user should give an identity matrix (Q.setIdentity())
0207      *
0208      * \sa compute(const MatrixType&, bool)
0209      */
0210     template<typename HessMatrixType, typename OrthMatrixType>
0211     ComplexSchur& computeFromHessenberg(const HessMatrixType& matrixH, const OrthMatrixType& matrixQ,  bool computeU=true);
0212
0213     /** \brief Reports whether previous computation was successful.
0214       *
0215       * \returns \c Success if computation was successful, \c NoConvergence otherwise.
0216       */
0217     ComputationInfo info() const
0218     {
0219       eigen_assert(m_isInitialized && "ComplexSchur is not initialized.");
0220       return m_info;
0221     }
0222
0223     /** \brief Sets the maximum number of iterations allowed.
0224       *
0225       * If not specified by the user, the maximum number of iterations is m_maxIterationsPerRow times the size
0226       * of the matrix.
0227       */
0228     ComplexSchur& setMaxIterations(Index maxIters)
0229     {
0230       m_maxIters = maxIters;
0231       return *this;
0232     }
0233
0234     /** \brief Returns the maximum number of iterations. */
0235     Index getMaxIterations()
0236     {
0237       return m_maxIters;
0238     }
0239
0240     /** \brief Maximum number of iterations per row.
0241       *
0242       * If not otherwise specified, the maximum number of iterations is this number times the size of the
0243       * matrix. It is currently set to 30.
0244       */
0245     static const int m_maxIterationsPerRow = 30;
0246
0247   protected:
0248     ComplexMatrixType m_matT, m_matU;
0249     HessenbergDecomposition<MatrixType> m_hess;
0250     ComputationInfo m_info;
0251     bool m_isInitialized;
0252     bool m_matUisUptodate;
0253     Index m_maxIters;
0254
0255   private:
0256     bool subdiagonalEntryIsNeglegible(Index i);
0257     ComplexScalar computeShift(Index iu, Index iter);
0258     void reduceToTriangularForm(bool computeU);
0259     friend struct internal::complex_schur_reduce_to_hessenberg<MatrixType, NumTraits<Scalar>::IsComplex>;
0260 };
0261
0262 /** If m_matT(i+1,i) is neglegible in floating point arithmetic
0263   * compared to m_matT(i,i) and m_matT(j,j), then set it to zero and
0264   * return true, else return false. */
0265 template<typename MatrixType>
0266 inline bool ComplexSchur<MatrixType>::subdiagonalEntryIsNeglegible(Index i)
0267 {
0268   RealScalar d = numext::norm1(m_matT.coeff(i,i)) + numext::norm1(m_matT.coeff(i+1,i+1));
0269   RealScalar sd = numext::norm1(m_matT.coeff(i+1,i));
0270   if (internal::isMuchSmallerThan(sd, d, NumTraits<RealScalar>::epsilon()))
0271   {
0272     m_matT.coeffRef(i+1,i) = ComplexScalar(0);
0273     return true;
0274   }
0275   return false;
0276 }
0277
0278
0279 /** Compute the shift in the current QR iteration. */
0280 template<typename MatrixType>
0281 typename ComplexSchur<MatrixType>::ComplexScalar ComplexSchur<MatrixType>::computeShift(Index iu, Index iter)
0282 {
0283   using std::abs;
0284   if (iter == 10 || iter == 20)
0285   {
0286     // exceptional shift, taken from http://www.netlib.org/eispack/comqr.f
0287     return abs(numext::real(m_matT.coeff(iu,iu-1))) + abs(numext::real(m_matT.coeff(iu-1,iu-2)));
0288   }
0289
0290   // compute the shift as one of the eigenvalues of t, the 2x2
0291   // diagonal block on the bottom of the active submatrix
0292   Matrix<ComplexScalar,2,2> t = m_matT.template block<2,2>(iu-1,iu-1);
0293   RealScalar normt = t.cwiseAbs().sum();
0294   t /= normt;     // the normalization by sf is to avoid under/overflow
0295
0296   ComplexScalar b = t.coeff(0,1) * t.coeff(1,0);
0297   ComplexScalar c = t.coeff(0,0) - t.coeff(1,1);
0298   ComplexScalar disc = sqrt(c*c + RealScalar(4)*b);
0299   ComplexScalar det = t.coeff(0,0) * t.coeff(1,1) - b;
0300   ComplexScalar trace = t.coeff(0,0) + t.coeff(1,1);
0301   ComplexScalar eival1 = (trace + disc) / RealScalar(2);
0302   ComplexScalar eival2 = (trace - disc) / RealScalar(2);
0303   RealScalar eival1_norm = numext::norm1(eival1);
0304   RealScalar eival2_norm = numext::norm1(eival2);
0305   // A division by zero can only occur if eival1==eival2==0.
0306   // In this case, det==0, and all we have to do is checking that eival2_norm!=0
0307   if(eival1_norm > eival2_norm)
0308     eival2 = det / eival1;
0309   else if(eival2_norm!=RealScalar(0))
0310     eival1 = det / eival2;
0311
0312   // choose the eigenvalue closest to the bottom entry of the diagonal
0313   if(numext::norm1(eival1-t.coeff(1,1)) < numext::norm1(eival2-t.coeff(1,1)))
0314     return normt * eival1;
0315   else
0316     return normt * eival2;
0317 }
0318
0319
0320 template<typename MatrixType>
0321 template<typename InputType>
0322 ComplexSchur<MatrixType>& ComplexSchur<MatrixType>::compute(const EigenBase<InputType>& matrix, bool computeU)
0323 {
0324   m_matUisUptodate = false;
0325   eigen_assert(matrix.cols() == matrix.rows());
0326
0327   if(matrix.cols() == 1)
0328   {
0329     m_matT = matrix.derived().template cast<ComplexScalar>();
0330     if(computeU)  m_matU = ComplexMatrixType::Identity(1,1);
0331     m_info = Success;
0332     m_isInitialized = true;
0333     m_matUisUptodate = computeU;
0334     return *this;
0335   }
0336
0337   internal::complex_schur_reduce_to_hessenberg<MatrixType, NumTraits<Scalar>::IsComplex>::run(*this, matrix.derived(), computeU);
0338   computeFromHessenberg(m_matT, m_matU, computeU);
0339   return *this;
0340 }
0341
0342 template<typename MatrixType>
0343 template<typename HessMatrixType, typename OrthMatrixType>
0344 ComplexSchur<MatrixType>& ComplexSchur<MatrixType>::computeFromHessenberg(const HessMatrixType& matrixH, const OrthMatrixType& matrixQ, bool computeU)
0345 {
0346   m_matT = matrixH;
0347   if(computeU)
0348     m_matU = matrixQ;
0349   reduceToTriangularForm(computeU);
0350   return *this;
0351 }
0352 namespace internal {
0353
0354 /* Reduce given matrix to Hessenberg form */
0355 template<typename MatrixType, bool IsComplex>
0356 struct complex_schur_reduce_to_hessenberg
0357 {
0358   // this is the implementation for the case IsComplex = true
0359   static void run(ComplexSchur<MatrixType>& _this, const MatrixType& matrix, bool computeU)
0360   {
0361     _this.m_hess.compute(matrix);
0362     _this.m_matT = _this.m_hess.matrixH();
0363     if(computeU)  _this.m_matU = _this.m_hess.matrixQ();
0364   }
0365 };
0366
0367 template<typename MatrixType>
0368 struct complex_schur_reduce_to_hessenberg<MatrixType, false>
0369 {
0370   static void run(ComplexSchur<MatrixType>& _this, const MatrixType& matrix, bool computeU)
0371   {
0372     typedef typename ComplexSchur<MatrixType>::ComplexScalar ComplexScalar;
0373
0374     // Note: m_hess is over RealScalar; m_matT and m_matU is over ComplexScalar
0375     _this.m_hess.compute(matrix);
0376     _this.m_matT = _this.m_hess.matrixH().template cast<ComplexScalar>();
0377     if(computeU)
0378     {
0379       // This may cause an allocation which seems to be avoidable
0380       MatrixType Q = _this.m_hess.matrixQ();
0381       _this.m_matU = Q.template cast<ComplexScalar>();
0382     }
0383   }
0384 };
0385
0386 } // end namespace internal
0387
0388 // Reduce the Hessenberg matrix m_matT to triangular form by QR iteration.
0389 template<typename MatrixType>
0390 void ComplexSchur<MatrixType>::reduceToTriangularForm(bool computeU)
0391 {
0392   Index maxIters = m_maxIters;
0393   if (maxIters == -1)
0394     maxIters = m_maxIterationsPerRow * m_matT.rows();
0395
0396   // The matrix m_matT is divided in three parts.
0397   // Rows 0,...,il-1 are decoupled from the rest because m_matT(il,il-1) is zero.
0398   // Rows il,...,iu is the part we are working on (the active submatrix).
0399   // Rows iu+1,...,end are already brought in triangular form.
0400   Index iu = m_matT.cols() - 1;
0401   Index il;
0402   Index iter = 0; // number of iterations we are working on the (iu,iu) element
0403   Index totalIter = 0; // number of iterations for whole matrix
0404
0405   while(true)
0406   {
0407     // find iu, the bottom row of the active submatrix
0408     while(iu > 0)
0409     {
0410       if(!subdiagonalEntryIsNeglegible(iu-1)) break;
0411       iter = 0;
0412       --iu;
0413     }
0414
0415     // if iu is zero then we are done; the whole matrix is triangularized
0416     if(iu==0) break;
0417
0418     // if we spent too many iterations, we give up
0419     iter++;
0420     totalIter++;
0421     if(totalIter > maxIters) break;
0422
0423     // find il, the top row of the active submatrix
0424     il = iu-1;
0425     while(il > 0 && !subdiagonalEntryIsNeglegible(il-1))
0426     {
0427       --il;
0428     }
0429
0430     /* perform the QR step using Givens rotations. The first rotation
0431        creates a bulge; the (il+2,il) element becomes nonzero. This
0432        bulge is chased down to the bottom of the active submatrix. */
0433
0434     ComplexScalar shift = computeShift(iu, iter);
0435     JacobiRotation<ComplexScalar> rot;
0436     rot.makeGivens(m_matT.coeff(il,il) - shift, m_matT.coeff(il+1,il));
0437     m_matT.rightCols(m_matT.cols()-il).applyOnTheLeft(il, il+1, rot.adjoint());
0438     m_matT.topRows((std::min)(il+2,iu)+1).applyOnTheRight(il, il+1, rot);
0439     if(computeU) m_matU.applyOnTheRight(il, il+1, rot);
0440
0441     for(Index i=il+1 ; i<iu ; i++)
0442     {
0443       rot.makeGivens(m_matT.coeffRef(i,i-1), m_matT.coeffRef(i+1,i-1), &m_matT.coeffRef(i,i-1));
0444       m_matT.coeffRef(i+1,i-1) = ComplexScalar(0);
0445       m_matT.rightCols(m_matT.cols()-i).applyOnTheLeft(i, i+1, rot.adjoint());
0446       m_matT.topRows((std::min)(i+2,iu)+1).applyOnTheRight(i, i+1, rot);
0447       if(computeU) m_matU.applyOnTheRight(i, i+1, rot);
0448     }
0449   }
0450
0451   if(totalIter <= maxIters)
0452     m_info = Success;
0453   else
0454     m_info = NoConvergence;
0455
0456   m_isInitialized = true;
0457   m_matUisUptodate = computeU;
0458 }
0459
0460 } // end namespace Eigen
0461
0462 #endif // EIGEN_COMPLEX_SCHUR_H