1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2008 Gael Guennebaud <gael.guennebaud (at) inria.fr> 5 // Copyright (C) 2010 Jitse Niesen <jitse (at) maths.leeds.ac.uk> 6 // 7 // This Source Code Form is subject to the terms of the Mozilla 8 // Public License v. 2.0. If a copy of the MPL was not distributed 9 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 10 11 #ifndef EIGEN_REAL_SCHUR_H 12 #define EIGEN_REAL_SCHUR_H 13 14 #include "./HessenbergDecomposition.h" 15 16 namespace Eigen { 17 18 /** \eigenvalues_module \ingroup Eigenvalues_Module 19 * 20 * 21 * \class RealSchur 22 * 23 * \brief Performs a real Schur decomposition of a square matrix 24 * 25 * \tparam _MatrixType the type of the matrix of which we are computing the 26 * real Schur decomposition; this is expected to be an instantiation of the 27 * Matrix class template. 28 * 29 * Given a real square matrix A, this class computes the real Schur 30 * decomposition: \f$ A = U T U^T \f$ where U is a real orthogonal matrix and 31 * T is a real quasi-triangular matrix. An orthogonal matrix is a matrix whose 32 * inverse is equal to its transpose, \f$ U^{-1} = U^T \f$. A quasi-triangular 33 * matrix is a block-triangular matrix whose diagonal consists of 1-by-1 34 * blocks and 2-by-2 blocks with complex eigenvalues. The eigenvalues of the 35 * blocks on the diagonal of T are the same as the eigenvalues of the matrix 36 * A, and thus the real Schur decomposition is used in EigenSolver to compute 37 * the eigendecomposition of a matrix. 38 * 39 * Call the function compute() to compute the real Schur decomposition of a 40 * given matrix. Alternatively, you can use the RealSchur(const MatrixType&, bool) 41 * constructor which computes the real Schur decomposition at construction 42 * time. Once the decomposition is computed, you can use the matrixU() and 43 * matrixT() functions to retrieve the matrices U and T in the decomposition. 44 * 45 * The documentation of RealSchur(const MatrixType&, bool) contains an example 46 * of the typical use of this class. 47 * 48 * \note The implementation is adapted from 49 * <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a> (public domain). 50 * Their code is based on EISPACK. 51 * 52 * \sa class ComplexSchur, class EigenSolver, class ComplexEigenSolver 53 */ 54 template<typename _MatrixType> class RealSchur 55 { 56 public: 57 typedef _MatrixType MatrixType; 58 enum { 59 RowsAtCompileTime = MatrixType::RowsAtCompileTime, 60 ColsAtCompileTime = MatrixType::ColsAtCompileTime, 61 Options = MatrixType::Options, 62 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, 63 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime 64 }; 65 typedef typename MatrixType::Scalar Scalar; 66 typedef std::complex<typename NumTraits<Scalar>::Real> ComplexScalar; 67 typedef typename MatrixType::Index Index; 68 69 typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> EigenvalueType; 70 typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ColumnVectorType; 71 72 /** \brief Default constructor. 73 * 74 * \param [in] size Positive integer, size of the matrix whose Schur decomposition will be computed. 75 * 76 * The default constructor is useful in cases in which the user intends to 77 * perform decompositions via compute(). The \p size parameter is only 78 * used as a hint. It is not an error to give a wrong \p size, but it may 79 * impair performance. 80 * 81 * \sa compute() for an example. 82 */ 83 RealSchur(Index size = RowsAtCompileTime==Dynamic ? 1 : RowsAtCompileTime) 84 : m_matT(size, size), 85 m_matU(size, size), 86 m_workspaceVector(size), 87 m_hess(size), 88 m_isInitialized(false), 89 m_matUisUptodate(false) 90 { } 91 92 /** \brief Constructor; computes real Schur decomposition of given matrix. 93 * 94 * \param[in] matrix Square matrix whose Schur decomposition is to be computed. 95 * \param[in] computeU If true, both T and U are computed; if false, only T is computed. 96 * 97 * This constructor calls compute() to compute the Schur decomposition. 98 * 99 * Example: \include RealSchur_RealSchur_MatrixType.cpp 100 * Output: \verbinclude RealSchur_RealSchur_MatrixType.out 101 */ 102 RealSchur(const MatrixType& matrix, bool computeU = true) 103 : m_matT(matrix.rows(),matrix.cols()), 104 m_matU(matrix.rows(),matrix.cols()), 105 m_workspaceVector(matrix.rows()), 106 m_hess(matrix.rows()), 107 m_isInitialized(false), 108 m_matUisUptodate(false) 109 { 110 compute(matrix, computeU); 111 } 112 113 /** \brief Returns the orthogonal matrix in the Schur decomposition. 114 * 115 * \returns A const reference to the matrix U. 116 * 117 * \pre Either the constructor RealSchur(const MatrixType&, bool) or the 118 * member function compute(const MatrixType&, bool) has been called before 119 * to compute the Schur decomposition of a matrix, and \p computeU was set 120 * to true (the default value). 121 * 122 * \sa RealSchur(const MatrixType&, bool) for an example 123 */ 124 const MatrixType& matrixU() const 125 { 126 eigen_assert(m_isInitialized && "RealSchur is not initialized."); 127 eigen_assert(m_matUisUptodate && "The matrix U has not been computed during the RealSchur decomposition."); 128 return m_matU; 129 } 130 131 /** \brief Returns the quasi-triangular matrix in the Schur decomposition. 132 * 133 * \returns A const reference to the matrix T. 134 * 135 * \pre Either the constructor RealSchur(const MatrixType&, bool) or the 136 * member function compute(const MatrixType&, bool) has been called before 137 * to compute the Schur decomposition of a matrix. 138 * 139 * \sa RealSchur(const MatrixType&, bool) for an example 140 */ 141 const MatrixType& matrixT() const 142 { 143 eigen_assert(m_isInitialized && "RealSchur is not initialized."); 144 return m_matT; 145 } 146 147 /** \brief Computes Schur decomposition of given matrix. 148 * 149 * \param[in] matrix Square matrix whose Schur decomposition is to be computed. 150 * \param[in] computeU If true, both T and U are computed; if false, only T is computed. 151 * \returns Reference to \c *this 152 * 153 * The Schur decomposition is computed by first reducing the matrix to 154 * Hessenberg form using the class HessenbergDecomposition. The Hessenberg 155 * matrix is then reduced to triangular form by performing Francis QR 156 * iterations with implicit double shift. The cost of computing the Schur 157 * decomposition depends on the number of iterations; as a rough guide, it 158 * may be taken to be \f$25n^3\f$ flops if \a computeU is true and 159 * \f$10n^3\f$ flops if \a computeU is false. 160 * 161 * Example: \include RealSchur_compute.cpp 162 * Output: \verbinclude RealSchur_compute.out 163 */ 164 RealSchur& compute(const MatrixType& matrix, bool computeU = true); 165 166 /** \brief Reports whether previous computation was successful. 167 * 168 * \returns \c Success if computation was succesful, \c NoConvergence otherwise. 169 */ 170 ComputationInfo info() const 171 { 172 eigen_assert(m_isInitialized && "RealSchur is not initialized."); 173 return m_info; 174 } 175 176 /** \brief Maximum number of iterations. 177 * 178 * Maximum number of iterations allowed for an eigenvalue to converge. 179 */ 180 static const int m_maxIterations = 40; 181 182 private: 183 184 MatrixType m_matT; 185 MatrixType m_matU; 186 ColumnVectorType m_workspaceVector; 187 HessenbergDecomposition<MatrixType> m_hess; 188 ComputationInfo m_info; 189 bool m_isInitialized; 190 bool m_matUisUptodate; 191 192 typedef Matrix<Scalar,3,1> Vector3s; 193 194 Scalar computeNormOfT(); 195 Index findSmallSubdiagEntry(Index iu, Scalar norm); 196 void splitOffTwoRows(Index iu, bool computeU, Scalar exshift); 197 void computeShift(Index iu, Index iter, Scalar& exshift, Vector3s& shiftInfo); 198 void initFrancisQRStep(Index il, Index iu, const Vector3s& shiftInfo, Index& im, Vector3s& firstHouseholderVector); 199 void performFrancisQRStep(Index il, Index im, Index iu, bool computeU, const Vector3s& firstHouseholderVector, Scalar* workspace); 200 }; 201 202 203 template<typename MatrixType> 204 RealSchur<MatrixType>& RealSchur<MatrixType>::compute(const MatrixType& matrix, bool computeU) 205 { 206 assert(matrix.cols() == matrix.rows()); 207 208 // Step 1. Reduce to Hessenberg form 209 m_hess.compute(matrix); 210 m_matT = m_hess.matrixH(); 211 if (computeU) 212 m_matU = m_hess.matrixQ(); 213 214 // Step 2. Reduce to real Schur form 215 m_workspaceVector.resize(m_matT.cols()); 216 Scalar* workspace = &m_workspaceVector.coeffRef(0); 217 218 // The matrix m_matT is divided in three parts. 219 // Rows 0,...,il-1 are decoupled from the rest because m_matT(il,il-1) is zero. 220 // Rows il,...,iu is the part we are working on (the active window). 221 // Rows iu+1,...,end are already brought in triangular form. 222 Index iu = m_matT.cols() - 1; 223 Index iter = 0; // iteration count 224 Scalar exshift(0); // sum of exceptional shifts 225 Scalar norm = computeNormOfT(); 226 227 if(norm!=0) 228 { 229 while (iu >= 0) 230 { 231 Index il = findSmallSubdiagEntry(iu, norm); 232 233 // Check for convergence 234 if (il == iu) // One root found 235 { 236 m_matT.coeffRef(iu,iu) = m_matT.coeff(iu,iu) + exshift; 237 if (iu > 0) 238 m_matT.coeffRef(iu, iu-1) = Scalar(0); 239 iu--; 240 iter = 0; 241 } 242 else if (il == iu-1) // Two roots found 243 { 244 splitOffTwoRows(iu, computeU, exshift); 245 iu -= 2; 246 iter = 0; 247 } 248 else // No convergence yet 249 { 250 // The firstHouseholderVector vector has to be initialized to something to get rid of a silly GCC warning (-O1 -Wall -DNDEBUG ) 251 Vector3s firstHouseholderVector(0,0,0), shiftInfo; 252 computeShift(iu, iter, exshift, shiftInfo); 253 iter = iter + 1; 254 if (iter > m_maxIterations * m_matT.cols()) break; 255 Index im; 256 initFrancisQRStep(il, iu, shiftInfo, im, firstHouseholderVector); 257 performFrancisQRStep(il, im, iu, computeU, firstHouseholderVector, workspace); 258 } 259 } 260 } 261 if(iter <= m_maxIterations * m_matT.cols()) 262 m_info = Success; 263 else 264 m_info = NoConvergence; 265 266 m_isInitialized = true; 267 m_matUisUptodate = computeU; 268 return *this; 269 } 270 271 /** \internal Computes and returns vector L1 norm of T */ 272 template<typename MatrixType> 273 inline typename MatrixType::Scalar RealSchur<MatrixType>::computeNormOfT() 274 { 275 const Index size = m_matT.cols(); 276 // FIXME to be efficient the following would requires a triangular reduxion code 277 // Scalar norm = m_matT.upper().cwiseAbs().sum() 278 // + m_matT.bottomLeftCorner(size-1,size-1).diagonal().cwiseAbs().sum(); 279 Scalar norm(0); 280 for (Index j = 0; j < size; ++j) 281 norm += m_matT.row(j).segment((std::max)(j-1,Index(0)), size-(std::max)(j-1,Index(0))).cwiseAbs().sum(); 282 return norm; 283 } 284 285 /** \internal Look for single small sub-diagonal element and returns its index */ 286 template<typename MatrixType> 287 inline typename MatrixType::Index RealSchur<MatrixType>::findSmallSubdiagEntry(Index iu, Scalar norm) 288 { 289 Index res = iu; 290 while (res > 0) 291 { 292 Scalar s = internal::abs(m_matT.coeff(res-1,res-1)) + internal::abs(m_matT.coeff(res,res)); 293 if (s == 0.0) 294 s = norm; 295 if (internal::abs(m_matT.coeff(res,res-1)) < NumTraits<Scalar>::epsilon() * s) 296 break; 297 res--; 298 } 299 return res; 300 } 301 302 /** \internal Update T given that rows iu-1 and iu decouple from the rest. */ 303 template<typename MatrixType> 304 inline void RealSchur<MatrixType>::splitOffTwoRows(Index iu, bool computeU, Scalar exshift) 305 { 306 const Index size = m_matT.cols(); 307 308 // The eigenvalues of the 2x2 matrix [a b; c d] are 309 // trace +/- sqrt(discr/4) where discr = tr^2 - 4*det, tr = a + d, det = ad - bc 310 Scalar p = Scalar(0.5) * (m_matT.coeff(iu-1,iu-1) - m_matT.coeff(iu,iu)); 311 Scalar q = p * p + m_matT.coeff(iu,iu-1) * m_matT.coeff(iu-1,iu); // q = tr^2 / 4 - det = discr/4 312 m_matT.coeffRef(iu,iu) += exshift; 313 m_matT.coeffRef(iu-1,iu-1) += exshift; 314 315 if (q >= Scalar(0)) // Two real eigenvalues 316 { 317 Scalar z = internal::sqrt(internal::abs(q)); 318 JacobiRotation<Scalar> rot; 319 if (p >= Scalar(0)) 320 rot.makeGivens(p + z, m_matT.coeff(iu, iu-1)); 321 else 322 rot.makeGivens(p - z, m_matT.coeff(iu, iu-1)); 323 324 m_matT.rightCols(size-iu+1).applyOnTheLeft(iu-1, iu, rot.adjoint()); 325 m_matT.topRows(iu+1).applyOnTheRight(iu-1, iu, rot); 326 m_matT.coeffRef(iu, iu-1) = Scalar(0); 327 if (computeU) 328 m_matU.applyOnTheRight(iu-1, iu, rot); 329 } 330 331 if (iu > 1) 332 m_matT.coeffRef(iu-1, iu-2) = Scalar(0); 333 } 334 335 /** \internal Form shift in shiftInfo, and update exshift if an exceptional shift is performed. */ 336 template<typename MatrixType> 337 inline void RealSchur<MatrixType>::computeShift(Index iu, Index iter, Scalar& exshift, Vector3s& shiftInfo) 338 { 339 shiftInfo.coeffRef(0) = m_matT.coeff(iu,iu); 340 shiftInfo.coeffRef(1) = m_matT.coeff(iu-1,iu-1); 341 shiftInfo.coeffRef(2) = m_matT.coeff(iu,iu-1) * m_matT.coeff(iu-1,iu); 342 343 // Wilkinson's original ad hoc shift 344 if (iter == 10) 345 { 346 exshift += shiftInfo.coeff(0); 347 for (Index i = 0; i <= iu; ++i) 348 m_matT.coeffRef(i,i) -= shiftInfo.coeff(0); 349 Scalar s = internal::abs(m_matT.coeff(iu,iu-1)) + internal::abs(m_matT.coeff(iu-1,iu-2)); 350 shiftInfo.coeffRef(0) = Scalar(0.75) * s; 351 shiftInfo.coeffRef(1) = Scalar(0.75) * s; 352 shiftInfo.coeffRef(2) = Scalar(-0.4375) * s * s; 353 } 354 355 // MATLAB's new ad hoc shift 356 if (iter == 30) 357 { 358 Scalar s = (shiftInfo.coeff(1) - shiftInfo.coeff(0)) / Scalar(2.0); 359 s = s * s + shiftInfo.coeff(2); 360 if (s > Scalar(0)) 361 { 362 s = internal::sqrt(s); 363 if (shiftInfo.coeff(1) < shiftInfo.coeff(0)) 364 s = -s; 365 s = s + (shiftInfo.coeff(1) - shiftInfo.coeff(0)) / Scalar(2.0); 366 s = shiftInfo.coeff(0) - shiftInfo.coeff(2) / s; 367 exshift += s; 368 for (Index i = 0; i <= iu; ++i) 369 m_matT.coeffRef(i,i) -= s; 370 shiftInfo.setConstant(Scalar(0.964)); 371 } 372 } 373 } 374 375 /** \internal Compute index im at which Francis QR step starts and the first Householder vector. */ 376 template<typename MatrixType> 377 inline void RealSchur<MatrixType>::initFrancisQRStep(Index il, Index iu, const Vector3s& shiftInfo, Index& im, Vector3s& firstHouseholderVector) 378 { 379 Vector3s& v = firstHouseholderVector; // alias to save typing 380 381 for (im = iu-2; im >= il; --im) 382 { 383 const Scalar Tmm = m_matT.coeff(im,im); 384 const Scalar r = shiftInfo.coeff(0) - Tmm; 385 const Scalar s = shiftInfo.coeff(1) - Tmm; 386 v.coeffRef(0) = (r * s - shiftInfo.coeff(2)) / m_matT.coeff(im+1,im) + m_matT.coeff(im,im+1); 387 v.coeffRef(1) = m_matT.coeff(im+1,im+1) - Tmm - r - s; 388 v.coeffRef(2) = m_matT.coeff(im+2,im+1); 389 if (im == il) { 390 break; 391 } 392 const Scalar lhs = m_matT.coeff(im,im-1) * (internal::abs(v.coeff(1)) + internal::abs(v.coeff(2))); 393 const Scalar rhs = v.coeff(0) * (internal::abs(m_matT.coeff(im-1,im-1)) + internal::abs(Tmm) + internal::abs(m_matT.coeff(im+1,im+1))); 394 if (internal::abs(lhs) < NumTraits<Scalar>::epsilon() * rhs) 395 { 396 break; 397 } 398 } 399 } 400 401 /** \internal Perform a Francis QR step involving rows il:iu and columns im:iu. */ 402 template<typename MatrixType> 403 inline void RealSchur<MatrixType>::performFrancisQRStep(Index il, Index im, Index iu, bool computeU, const Vector3s& firstHouseholderVector, Scalar* workspace) 404 { 405 assert(im >= il); 406 assert(im <= iu-2); 407 408 const Index size = m_matT.cols(); 409 410 for (Index k = im; k <= iu-2; ++k) 411 { 412 bool firstIteration = (k == im); 413 414 Vector3s v; 415 if (firstIteration) 416 v = firstHouseholderVector; 417 else 418 v = m_matT.template block<3,1>(k,k-1); 419 420 Scalar tau, beta; 421 Matrix<Scalar, 2, 1> ess; 422 v.makeHouseholder(ess, tau, beta); 423 424 if (beta != Scalar(0)) // if v is not zero 425 { 426 if (firstIteration && k > il) 427 m_matT.coeffRef(k,k-1) = -m_matT.coeff(k,k-1); 428 else if (!firstIteration) 429 m_matT.coeffRef(k,k-1) = beta; 430 431 // These Householder transformations form the O(n^3) part of the algorithm 432 m_matT.block(k, k, 3, size-k).applyHouseholderOnTheLeft(ess, tau, workspace); 433 m_matT.block(0, k, (std::min)(iu,k+3) + 1, 3).applyHouseholderOnTheRight(ess, tau, workspace); 434 if (computeU) 435 m_matU.block(0, k, size, 3).applyHouseholderOnTheRight(ess, tau, workspace); 436 } 437 } 438 439 Matrix<Scalar, 2, 1> v = m_matT.template block<2,1>(iu-1, iu-2); 440 Scalar tau, beta; 441 Matrix<Scalar, 1, 1> ess; 442 v.makeHouseholder(ess, tau, beta); 443 444 if (beta != Scalar(0)) // if v is not zero 445 { 446 m_matT.coeffRef(iu-1, iu-2) = beta; 447 m_matT.block(iu-1, iu-1, 2, size-iu+1).applyHouseholderOnTheLeft(ess, tau, workspace); 448 m_matT.block(0, iu-1, iu+1, 2).applyHouseholderOnTheRight(ess, tau, workspace); 449 if (computeU) 450 m_matU.block(0, iu-1, size, 2).applyHouseholderOnTheRight(ess, tau, workspace); 451 } 452 453 // clean up pollution due to round-off errors 454 for (Index i = im+2; i <= iu; ++i) 455 { 456 m_matT.coeffRef(i,i-2) = Scalar(0); 457 if (i > im+2) 458 m_matT.coeffRef(i,i-3) = Scalar(0); 459 } 460 } 461 462 } // end namespace Eigen 463 464 #endif // EIGEN_REAL_SCHUR_H 465
