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@inria.fr>
5 // Copyright (C) 2010 Jitse Niesen <jitse@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)
m_matT(size,size)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 */
matrixU()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 */
matrixT()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 */
info()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>
compute(const MatrixType & matrix,bool computeU)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>
computeNormOfT()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>
findSmallSubdiagEntry(Index iu,Scalar norm)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>
splitOffTwoRows(Index iu,bool computeU,Scalar exshift)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>
computeShift(Index iu,Index iter,Scalar & exshift,Vector3s & shiftInfo)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>
initFrancisQRStep(Index il,Index iu,const Vector3s & shiftInfo,Index & im,Vector3s & firstHouseholderVector)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>
performFrancisQRStep(Index il,Index im,Index iu,bool computeU,const Vector3s & firstHouseholderVector,Scalar * workspace)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