1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2008-2010 Gael Guennebaud <gael.guennebaud@inria.fr>
5 // Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
6 //
7 // Eigen is free software; you can redistribute it and/or
8 // modify it under the terms of the GNU Lesser General Public
9 // License as published by the Free Software Foundation; either
10 // version 3 of the License, or (at your option) any later version.
11 //
12 // Alternatively, you can redistribute it and/or
13 // modify it under the terms of the GNU General Public License as
14 // published by the Free Software Foundation; either version 2 of
15 // the License, or (at your option) any later version.
16 //
17 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
18 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
19 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
20 // GNU General Public License for more details.
21 //
22 // You should have received a copy of the GNU Lesser General Public
23 // License and a copy of the GNU General Public License along with
24 // Eigen. If not, see <http://www.gnu.org/licenses/>.
29 #include "./EigenvaluesCommon.h"
30 #include "./Tridiagonalization.h"
32 template<typename _MatrixType>
35 /** \eigenvalues_module \ingroup Eigenvalues_Module
36   *
37   *
39   *
40   * \brief Computes eigenvalues and eigenvectors of selfadjoint matrices
41   *
42   * \tparam _MatrixType the type of the matrix of which we are computing the
43   * eigendecomposition; this is expected to be an instantiation of the Matrix
44   * class template.
45   *
46   * A matrix \f$A \f$ is selfadjoint if it equals its adjoint. For real
47   * matrices, this means that the matrix is symmetric: it equals its
48   * transpose. This class computes the eigenvalues and eigenvectors of a
49   * selfadjoint matrix. These are the scalars \f$\lambda \f$ and vectors
50   * \f$v \f$ such that \f$Av = \lambda v \f$.  The eigenvalues of a
51   * selfadjoint matrix are always real. If \f$D \f$ is a diagonal matrix with
52   * the eigenvalues on the diagonal, and \f$V \f$ is a matrix with the
53   * eigenvectors as its columns, then \f$A = V D V^{-1} \f$ (for selfadjoint
54   * matrices, the matrix \f$V \f$ is always invertible). This is called the
55   * eigendecomposition.
56   *
57   * The algorithm exploits the fact that the matrix is selfadjoint, making it
58   * faster and more accurate than the general purpose eigenvalue algorithms
59   * implemented in EigenSolver and ComplexEigenSolver.
60   *
61   * Only the \b lower \b triangular \b part of the input matrix is referenced.
62   *
63   * Call the function compute() to compute the eigenvalues and eigenvectors of
64   * a given matrix. Alternatively, you can use the
65   * SelfAdjointEigenSolver(const MatrixType&, int) constructor which computes
66   * the eigenvalues and eigenvectors at construction time. Once the eigenvalue
67   * and eigenvectors are computed, they can be retrieved with the eigenvalues()
68   * and eigenvectors() functions.
69   *
70   * The documentation for SelfAdjointEigenSolver(const MatrixType&, int)
71   * contains an example of the typical use of this class.
72   *
73   * To solve the \em generalized eigenvalue problem \f$Av = \lambda Bv \f$ and
74   * the likes, see the class GeneralizedSelfAdjointEigenSolver.
75   *
76   * \sa MatrixBase::eigenvalues(), class EigenSolver, class ComplexEigenSolver
77   */
79 {
80   public:
82     typedef _MatrixType MatrixType;
83     enum {
84       Size = MatrixType::RowsAtCompileTime,
85       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
86       Options = MatrixType::Options,
87       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
88     };
90     /** \brief Scalar type for matrices of type \p _MatrixType. */
91     typedef typename MatrixType::Scalar Scalar;
92     typedef typename MatrixType::Index Index;
94     /** \brief Real scalar type for \p _MatrixType.
95       *
96       * This is just \c Scalar if #Scalar is real (e.g., \c float or
97       * \c double), and the type of the real part of \c Scalar if #Scalar is
98       * complex.
99       */
100     typedef typename NumTraits<Scalar>::Real RealScalar;
102     /** \brief Type for vector of eigenvalues as returned by eigenvalues().
103       *
104       * This is a column vector with entries of type #RealScalar.
105       * The length of the vector is the size of \p _MatrixType.
106       */
107     typedef typename internal::plain_col_type<MatrixType, RealScalar>::type RealVectorType;
108     typedef Tridiagonalization<MatrixType> TridiagonalizationType;
110     /** \brief Default constructor for fixed-size matrices.
111       *
112       * The default constructor is useful in cases in which the user intends to
113       * perform decompositions via compute(). This constructor
114       * can only be used if \p _MatrixType is a fixed-size matrix; use
115       * SelfAdjointEigenSolver(Index) for dynamic-size matrices.
116       *
119       */
121         : m_eivec(),
122           m_eivalues(),
123           m_subdiag(),
124           m_isInitialized(false)
125     { }
127     /** \brief Constructor, pre-allocates memory for dynamic-size matrices.
128       *
129       * \param [in]  size  Positive integer, size of the matrix whose
130       * eigenvalues and eigenvectors will be computed.
131       *
132       * This constructor is useful for dynamic-size matrices, when the user
133       * intends to perform decompositions via compute(). The \p size
134       * parameter is only used as a hint. It is not an error to give a wrong
135       * \p size, but it may impair performance.
136       *
137       * \sa compute() for an example
138       */
140         : m_eivec(size, size),
141           m_eivalues(size),
142           m_subdiag(size > 1 ? size - 1 : 1),
143           m_isInitialized(false)
144     {}
146     /** \brief Constructor; computes eigendecomposition of given matrix.
147       *
148       * \param[in]  matrix  Selfadjoint matrix whose eigendecomposition is to
149       *    be computed. Only the lower triangular part of the matrix is referenced.
150       * \param[in]  options Can be #ComputeEigenvectors (default) or #EigenvaluesOnly.
151       *
152       * This constructor calls compute(const MatrixType&, int) to compute the
153       * eigenvalues of the matrix \p matrix. The eigenvectors are computed if
154       * \p options equals #ComputeEigenvectors.
155       *
158       *
159       * \sa compute(const MatrixType&, int)
160       */
161     SelfAdjointEigenSolver(const MatrixType& matrix, int options = ComputeEigenvectors)
162       : m_eivec(matrix.rows(), matrix.cols()),
163         m_eivalues(matrix.cols()),
164         m_subdiag(matrix.rows() > 1 ? matrix.rows() - 1 : 1),
165         m_isInitialized(false)
166     {
167       compute(matrix, options);
168     }
170     /** \brief Computes eigendecomposition of given matrix.
171       *
172       * \param[in]  matrix  Selfadjoint matrix whose eigendecomposition is to
173       *    be computed. Only the lower triangular part of the matrix is referenced.
174       * \param[in]  options Can be #ComputeEigenvectors (default) or #EigenvaluesOnly.
175       * \returns    Reference to \c *this
176       *
177       * This function computes the eigenvalues of \p matrix.  The eigenvalues()
178       * function can be used to retrieve them.  If \p options equals #ComputeEigenvectors,
179       * then the eigenvectors are also computed and can be retrieved by
180       * calling eigenvectors().
181       *
182       * This implementation uses a symmetric QR algorithm. The matrix is first
183       * reduced to tridiagonal form using the Tridiagonalization class. The
184       * tridiagonal matrix is then brought to diagonal form with implicit
185       * symmetric QR steps with Wilkinson shift. Details can be found in
186       * Section 8.3 of Golub \& Van Loan, <i>%Matrix Computations</i>.
187       *
188       * The cost of the computation is about \f$9n^3 \f$ if the eigenvectors
189       * are required and \f$4n^3/3 \f$ if they are not required.
190       *
191       * This method reuses the memory in the SelfAdjointEigenSolver object that
192       * was allocated when the object was constructed, if the size of the
193       * matrix does not change.
194       *
197       *
198       * \sa SelfAdjointEigenSolver(const MatrixType&, int)
199       */
200     SelfAdjointEigenSolver& compute(const MatrixType& matrix, int options = ComputeEigenvectors);
202     /** \brief Returns the eigenvectors of given matrix.
203       *
204       * \returns  A const reference to the matrix whose columns are the eigenvectors.
205       *
206       * \pre The eigenvectors have been computed before.
207       *
208       * Column \f$k \f$ of the returned matrix is an eigenvector corresponding
209       * to eigenvalue number \f$k \f$ as returned by eigenvalues().  The
210       * eigenvectors are normalized to have (Euclidean) norm equal to one. If
211       * this object was used to solve the eigenproblem for the selfadjoint
212       * matrix \f$A \f$, then the matrix returned by this function is the
213       * matrix \f$V \f$ in the eigendecomposition \f$A = V D V^{-1} \f$.
214       *
217       *
218       * \sa eigenvalues()
219       */
220     const MatrixType& eigenvectors() const
221     {
222       eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
223       eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
224       return m_eivec;
225     }
227     /** \brief Returns the eigenvalues of given matrix.
228       *
229       * \returns A const reference to the column vector containing the eigenvalues.
230       *
231       * \pre The eigenvalues have been computed before.
232       *
233       * The eigenvalues are repeated according to their algebraic multiplicity,
234       * so there are as many eigenvalues as rows in the matrix. The eigenvalues
235       * are sorted in increasing order.
236       *
239       *
240       * \sa eigenvectors(), MatrixBase::eigenvalues()
241       */
242     const RealVectorType& eigenvalues() const
243     {
244       eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
245       return m_eivalues;
246     }
248     /** \brief Computes the positive-definite square root of the matrix.
249       *
250       * \returns the positive-definite square root of the matrix
251       *
252       * \pre The eigenvalues and eigenvectors of a positive-definite matrix
253       * have been computed before.
254       *
255       * The square root of a positive-definite matrix \f$A \f$ is the
256       * positive-definite matrix whose square equals \f$A \f$. This function
257       * uses the eigendecomposition \f$A = V D V^{-1} \f$ to compute the
258       * square root as \f$A^{1/2} = V D^{1/2} V^{-1} \f$.
259       *
262       *
263       * \sa operatorInverseSqrt(),
264       *     \ref MatrixFunctions_Module "MatrixFunctions Module"
265       */
266     MatrixType operatorSqrt() const
267     {
268       eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
269       eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
270       return m_eivec * m_eivalues.cwiseSqrt().asDiagonal() * m_eivec.adjoint();
271     }
273     /** \brief Computes the inverse square root of the matrix.
274       *
275       * \returns the inverse positive-definite square root of the matrix
276       *
277       * \pre The eigenvalues and eigenvectors of a positive-definite matrix
278       * have been computed before.
279       *
280       * This function uses the eigendecomposition \f$A = V D V^{-1} \f$ to
281       * compute the inverse square root as \f$V D^{-1/2} V^{-1} \f$. This is
282       * cheaper than first computing the square root with operatorSqrt() and
283       * then its inverse with MatrixBase::inverse().
284       *
287       *
288       * \sa operatorSqrt(), MatrixBase::inverse(),
289       *     \ref MatrixFunctions_Module "MatrixFunctions Module"
290       */
291     MatrixType operatorInverseSqrt() const
292     {
293       eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
294       eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
295       return m_eivec * m_eivalues.cwiseInverse().cwiseSqrt().asDiagonal() * m_eivec.adjoint();
296     }
298     /** \brief Reports whether previous computation was successful.
299       *
300       * \returns \c Success if computation was succesful, \c NoConvergence otherwise.
301       */
302     ComputationInfo info() const
303     {
304       eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
305       return m_info;
306     }
308     /** \brief Maximum number of iterations.
309       *
310       * Maximum number of iterations allowed for an eigenvalue to converge.
311       */
312     static const int m_maxIterations = 30;
314     #ifdef EIGEN2_SUPPORT
315     SelfAdjointEigenSolver(const MatrixType& matrix, bool computeEigenvectors)
316       : m_eivec(matrix.rows(), matrix.cols()),
317         m_eivalues(matrix.cols()),
318         m_subdiag(matrix.rows() > 1 ? matrix.rows() - 1 : 1),
319         m_isInitialized(false)
320     {
321       compute(matrix, computeEigenvectors);
322     }
324     SelfAdjointEigenSolver(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors = true)
325         : m_eivec(matA.cols(), matA.cols()),
326           m_eivalues(matA.cols()),
327           m_subdiag(matA.cols() > 1 ? matA.cols() - 1 : 1),
328           m_isInitialized(false)
329     {
330       static_cast<GeneralizedSelfAdjointEigenSolver<MatrixType>*>(this)->compute(matA, matB, computeEigenvectors ? ComputeEigenvectors : EigenvaluesOnly);
331     }
333     void compute(const MatrixType& matrix, bool computeEigenvectors)
334     {
335       compute(matrix, computeEigenvectors ? ComputeEigenvectors : EigenvaluesOnly);
336     }
338     void compute(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors = true)
339     {
340       compute(matA, matB, computeEigenvectors ? ComputeEigenvectors : EigenvaluesOnly);
341     }
342     #endif // EIGEN2_SUPPORT
344   protected:
345     MatrixType m_eivec;
346     RealVectorType m_eivalues;
347     typename TridiagonalizationType::SubDiagonalType m_subdiag;
348     ComputationInfo m_info;
349     bool m_isInitialized;
350     bool m_eigenvectorsOk;
351 };
353 /** \internal
354   *
355   * \eigenvalues_module \ingroup Eigenvalues_Module
356   *
357   * Performs a QR step on a tridiagonal symmetric matrix represented as a
358   * pair of two vectors \a diag and \a subdiag.
359   *
360   * \param matA the input selfadjoint matrix
361   * \param hCoeffs returned Householder coefficients
362   *
363   * For compilation efficiency reasons, this procedure does not use eigen expression
364   * for its arguments.
365   *
366   * Implemented from Golub's "Matrix Computations", algorithm 8.3.2:
367   * "implicit symmetric QR step with Wilkinson shift"
368   */
369 namespace internal {
370 template<int StorageOrder,typename RealScalar, typename Scalar, typename Index>
371 static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index start, Index end, Scalar* matrixQ, Index n);
372 }
374 template<typename MatrixType>
376 ::compute(const MatrixType& matrix, int options)
377 {
378   eigen_assert(matrix.cols() == matrix.rows());
381           && "invalid option parameter");
382   bool computeEigenvectors = (options&ComputeEigenvectors)==ComputeEigenvectors;
383   Index n = matrix.cols();
384   m_eivalues.resize(n,1);
386   if(n==1)
387   {
388     m_eivalues.coeffRef(0,0) = internal::real(matrix.coeff(0,0));
389     if(computeEigenvectors)
390       m_eivec.setOnes(n,n);
391     m_info = Success;
392     m_isInitialized = true;
393     m_eigenvectorsOk = computeEigenvectors;
394     return *this;
395   }
397   // declare some aliases
398   RealVectorType& diag = m_eivalues;
399   MatrixType& mat = m_eivec;
401   // map the matrix coefficients to [-1:1] to avoid over- and underflow.
402   RealScalar scale = matrix.cwiseAbs().maxCoeff();
403   if(scale==Scalar(0)) scale = 1;
404   mat = matrix / scale;
405   m_subdiag.resize(n-1);
406   internal::tridiagonalization_inplace(mat, diag, m_subdiag, computeEigenvectors);
408   Index end = n-1;
409   Index start = 0;
410   Index iter = 0; // number of iterations we are working on one element
412   while (end>0)
413   {
414     for (Index i = start; i<end; ++i)
415       if (internal::isMuchSmallerThan(internal::abs(m_subdiag[i]),(internal::abs(diag[i])+internal::abs(diag[i+1]))))
416         m_subdiag[i] = 0;
418     // find the largest unreduced block
419     while (end>0 && m_subdiag[end-1]==0)
420     {
421       iter = 0;
422       end--;
423     }
424     if (end<=0)
425       break;
427     // if we spent too many iterations on the current element, we give up
428     iter++;
429     if(iter > m_maxIterations) break;
431     start = end - 1;
432     while (start>0 && m_subdiag[start-1]!=0)
433       start--;
435     internal::tridiagonal_qr_step<MatrixType::Flags&RowMajorBit ? RowMajor : ColMajor>(diag.data(), m_subdiag.data(), start, end, computeEigenvectors ? m_eivec.data() : (Scalar*)0, n);
436   }
438   if (iter <= m_maxIterations)
439     m_info = Success;
440   else
441     m_info = NoConvergence;
443   // Sort eigenvalues and corresponding vectors.
444   // TODO make the sort optional ?
445   // TODO use a better sort algorithm !!
446   if (m_info == Success)
447   {
448     for (Index i = 0; i < n-1; ++i)
449     {
450       Index k;
451       m_eivalues.segment(i,n-i).minCoeff(&k);
452       if (k > 0)
453       {
454         std::swap(m_eivalues[i], m_eivalues[k+i]);
455         if(computeEigenvectors)
456           m_eivec.col(i).swap(m_eivec.col(k+i));
457       }
458     }
459   }
461   // scale back the eigen values
462   m_eivalues *= scale;
464   m_isInitialized = true;
465   m_eigenvectorsOk = computeEigenvectors;
466   return *this;
467 }
469 namespace internal {
470 template<int StorageOrder,typename RealScalar, typename Scalar, typename Index>
471 static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index start, Index end, Scalar* matrixQ, Index n)
472 {
473   // NOTE this version avoids over & underflow, however since the matrix is prescaled, overflow cannot occur,
474   // and underflows should be meaningless anyway. So I don't any reason to enable this version, but I keep
475   // it here for reference:
476 //   RealScalar td = (diag[end-1] - diag[end])*RealScalar(0.5);
477 //   RealScalar e = subdiag[end-1];
478 //   RealScalar mu = diag[end] - (e / (td + (td>0 ? 1 : -1))) * (e / hypot(td,e));
479   RealScalar td = (diag[end-1] - diag[end])*RealScalar(0.5);
480   RealScalar e2 = abs2(subdiag[end-1]);
481   RealScalar mu = diag[end] - e2 / (td + (td>0 ? 1 : -1) * sqrt(td*td + e2));
482   RealScalar x = diag[start] - mu;
483   RealScalar z = subdiag[start];
484   for (Index k = start; k < end; ++k)
485   {
486     JacobiRotation<RealScalar> rot;
487     rot.makeGivens(x, z);
489     // do T = G' T G
490     RealScalar sdk = rot.s() * diag[k] + rot.c() * subdiag[k];
491     RealScalar dkp1 = rot.s() * subdiag[k] + rot.c() * diag[k+1];
493     diag[k] = rot.c() * (rot.c() * diag[k] - rot.s() * subdiag[k]) - rot.s() * (rot.c() * subdiag[k] - rot.s() * diag[k+1]);
494     diag[k+1] = rot.s() * sdk + rot.c() * dkp1;
495     subdiag[k] = rot.c() * sdk - rot.s() * dkp1;
498     if (k > start)
499       subdiag[k - 1] = rot.c() * subdiag[k-1] - rot.s() * z;
501     x = subdiag[k];
503     if (k < end - 1)
504     {
505       z = -rot.s() * subdiag[k+1];
506       subdiag[k + 1] = rot.c() * subdiag[k+1];
507     }
509     // apply the givens rotation to the unit matrix Q = Q * G
510     if (matrixQ)
511     {
512       // FIXME if StorageOrder == RowMajor this operation is not very efficient
513       Map<Matrix<Scalar,Dynamic,Dynamic,StorageOrder> > q(matrixQ,n,n);
514       q.applyOnTheRight(k,k+1,rot);
515     }
516   }
517 }
518 } // end namespace internal