Eigen  3.4.90 (git rev 67eeba6e720c5745abc77ae6c92ce0a44aa7b7ae)
EigenSolver.h
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,2012 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_EIGENSOLVER_H
12 #define EIGEN_EIGENSOLVER_H
13 
14 #include "./RealSchur.h"
15 
16 #include "./InternalHeaderCheck.h"
17 
18 namespace Eigen {
19 
66 template<typename MatrixType_> class EigenSolver
67 {
68  public:
69 
71  typedef MatrixType_ MatrixType;
72 
73  enum {
74  RowsAtCompileTime = MatrixType::RowsAtCompileTime,
75  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
76  Options = MatrixType::Options,
77  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
78  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
79  };
80 
82  typedef typename MatrixType::Scalar Scalar;
83  typedef typename NumTraits<Scalar>::Real RealScalar;
84  typedef Eigen::Index Index;
85 
92  typedef std::complex<RealScalar> ComplexScalar;
93 
100 
107 
115  EigenSolver() : m_eivec(), m_eivalues(), m_isInitialized(false), m_eigenvectorsOk(false), m_realSchur(), m_matT(), m_tmp() {}
116 
123  explicit EigenSolver(Index size)
124  : m_eivec(size, size),
125  m_eivalues(size),
126  m_isInitialized(false),
127  m_eigenvectorsOk(false),
128  m_realSchur(size),
129  m_matT(size, size),
130  m_tmp(size)
131  {}
132 
148  template<typename InputType>
149  explicit EigenSolver(const EigenBase<InputType>& matrix, bool computeEigenvectors = true)
150  : m_eivec(matrix.rows(), matrix.cols()),
151  m_eivalues(matrix.cols()),
152  m_isInitialized(false),
153  m_eigenvectorsOk(false),
154  m_realSchur(matrix.cols()),
155  m_matT(matrix.rows(), matrix.cols()),
156  m_tmp(matrix.cols())
157  {
158  compute(matrix.derived(), computeEigenvectors);
159  }
160 
182 
202  {
203  eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
204  eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
205  return m_eivec;
206  }
207 
227 
247  {
248  eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
249  return m_eivalues;
250  }
251 
279  template<typename InputType>
280  EigenSolver& compute(const EigenBase<InputType>& matrix, bool computeEigenvectors = true);
281 
284  {
285  eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
286  return m_info;
287  }
288 
291  {
292  m_realSchur.setMaxIterations(maxIters);
293  return *this;
294  }
295 
298  {
299  return m_realSchur.getMaxIterations();
300  }
301 
302  private:
303  void doComputeEigenvectors();
304 
305  protected:
306 
307  static void check_template_parameters()
308  {
309  EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
310  EIGEN_STATIC_ASSERT(!NumTraits<Scalar>::IsComplex, NUMERIC_TYPE_MUST_BE_REAL);
311  }
312 
313  MatrixType m_eivec;
314  EigenvalueType m_eivalues;
315  bool m_isInitialized;
316  bool m_eigenvectorsOk;
317  ComputationInfo m_info;
318  RealSchur<MatrixType> m_realSchur;
319  MatrixType m_matT;
320 
321  typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ColumnVectorType;
322  ColumnVectorType m_tmp;
323 };
324 
325 template<typename MatrixType>
327 {
328  eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
329  const RealScalar precision = RealScalar(2)*NumTraits<RealScalar>::epsilon();
330  Index n = m_eivalues.rows();
331  MatrixType matD = MatrixType::Zero(n,n);
332  for (Index i=0; i<n; ++i)
333  {
334  if (internal::isMuchSmallerThan(numext::imag(m_eivalues.coeff(i)), numext::real(m_eivalues.coeff(i)), precision))
335  matD.coeffRef(i,i) = numext::real(m_eivalues.coeff(i));
336  else
337  {
338  matD.template block<2,2>(i,i) << numext::real(m_eivalues.coeff(i)), numext::imag(m_eivalues.coeff(i)),
339  -numext::imag(m_eivalues.coeff(i)), numext::real(m_eivalues.coeff(i));
340  ++i;
341  }
342  }
343  return matD;
344 }
345 
346 template<typename MatrixType>
348 {
349  eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
350  eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
351  const RealScalar precision = RealScalar(2)*NumTraits<RealScalar>::epsilon();
352  Index n = m_eivec.cols();
353  EigenvectorsType matV(n,n);
354  for (Index j=0; j<n; ++j)
355  {
356  if (internal::isMuchSmallerThan(numext::imag(m_eivalues.coeff(j)), numext::real(m_eivalues.coeff(j)), precision) || j+1==n)
357  {
358  // we have a real eigen value
359  matV.col(j) = m_eivec.col(j).template cast<ComplexScalar>();
360  matV.col(j).normalize();
361  }
362  else
363  {
364  // we have a pair of complex eigen values
365  for (Index i=0; i<n; ++i)
366  {
367  matV.coeffRef(i,j) = ComplexScalar(m_eivec.coeff(i,j), m_eivec.coeff(i,j+1));
368  matV.coeffRef(i,j+1) = ComplexScalar(m_eivec.coeff(i,j), -m_eivec.coeff(i,j+1));
369  }
370  matV.col(j).normalize();
371  matV.col(j+1).normalize();
372  ++j;
373  }
374  }
375  return matV;
376 }
377 
378 template<typename MatrixType>
379 template<typename InputType>
381 EigenSolver<MatrixType>::compute(const EigenBase<InputType>& matrix, bool computeEigenvectors)
382 {
383  check_template_parameters();
384 
385  using std::sqrt;
386  using std::abs;
387  using numext::isfinite;
388  eigen_assert(matrix.cols() == matrix.rows());
389 
390  // Reduce to real Schur form.
391  m_realSchur.compute(matrix.derived(), computeEigenvectors);
392 
393  m_info = m_realSchur.info();
394 
395  if (m_info == Success)
396  {
397  m_matT = m_realSchur.matrixT();
398  if (computeEigenvectors)
399  m_eivec = m_realSchur.matrixU();
400 
401  // Compute eigenvalues from matT
402  m_eivalues.resize(matrix.cols());
403  Index i = 0;
404  while (i < matrix.cols())
405  {
406  if (i == matrix.cols() - 1 || m_matT.coeff(i+1, i) == Scalar(0))
407  {
408  m_eivalues.coeffRef(i) = m_matT.coeff(i, i);
409  if(!(isfinite)(m_eivalues.coeffRef(i)))
410  {
411  m_isInitialized = true;
412  m_eigenvectorsOk = false;
413  m_info = NumericalIssue;
414  return *this;
415  }
416  ++i;
417  }
418  else
419  {
420  Scalar p = Scalar(0.5) * (m_matT.coeff(i, i) - m_matT.coeff(i+1, i+1));
421  Scalar z;
422  // Compute z = sqrt(abs(p * p + m_matT.coeff(i+1, i) * m_matT.coeff(i, i+1)));
423  // without overflow
424  {
425  Scalar t0 = m_matT.coeff(i+1, i);
426  Scalar t1 = m_matT.coeff(i, i+1);
427  Scalar maxval = numext::maxi<Scalar>(abs(p),numext::maxi<Scalar>(abs(t0),abs(t1)));
428  t0 /= maxval;
429  t1 /= maxval;
430  Scalar p0 = p/maxval;
431  z = maxval * sqrt(abs(p0 * p0 + t0 * t1));
432  }
433 
434  m_eivalues.coeffRef(i) = ComplexScalar(m_matT.coeff(i+1, i+1) + p, z);
435  m_eivalues.coeffRef(i+1) = ComplexScalar(m_matT.coeff(i+1, i+1) + p, -z);
436  if(!((isfinite)(m_eivalues.coeffRef(i)) && (isfinite)(m_eivalues.coeffRef(i+1))))
437  {
438  m_isInitialized = true;
439  m_eigenvectorsOk = false;
440  m_info = NumericalIssue;
441  return *this;
442  }
443  i += 2;
444  }
445  }
446 
447  // Compute eigenvectors.
448  if (computeEigenvectors)
449  doComputeEigenvectors();
450  }
451 
452  m_isInitialized = true;
453  m_eigenvectorsOk = computeEigenvectors;
454 
455  return *this;
456 }
457 
458 
459 template<typename MatrixType>
460 void EigenSolver<MatrixType>::doComputeEigenvectors()
461 {
462  using std::abs;
463  const Index size = m_eivec.cols();
464  const Scalar eps = NumTraits<Scalar>::epsilon();
465 
466  // inefficient! this is already computed in RealSchur
467  Scalar norm(0);
468  for (Index j = 0; j < size; ++j)
469  {
470  norm += m_matT.row(j).segment((std::max)(j-1,Index(0)), size-(std::max)(j-1,Index(0))).cwiseAbs().sum();
471  }
472 
473  // Backsubstitute to find vectors of upper triangular form
474  if (norm == Scalar(0))
475  {
476  return;
477  }
478 
479  for (Index n = size-1; n >= 0; n--)
480  {
481  Scalar p = m_eivalues.coeff(n).real();
482  Scalar q = m_eivalues.coeff(n).imag();
483 
484  // Scalar vector
485  if (q == Scalar(0))
486  {
487  Scalar lastr(0), lastw(0);
488  Index l = n;
489 
490  m_matT.coeffRef(n,n) = Scalar(1);
491  for (Index i = n-1; i >= 0; i--)
492  {
493  Scalar w = m_matT.coeff(i,i) - p;
494  Scalar r = m_matT.row(i).segment(l,n-l+1).dot(m_matT.col(n).segment(l, n-l+1));
495 
496  if (m_eivalues.coeff(i).imag() < Scalar(0))
497  {
498  lastw = w;
499  lastr = r;
500  }
501  else
502  {
503  l = i;
504  if (m_eivalues.coeff(i).imag() == Scalar(0))
505  {
506  if (w != Scalar(0))
507  m_matT.coeffRef(i,n) = -r / w;
508  else
509  m_matT.coeffRef(i,n) = -r / (eps * norm);
510  }
511  else // Solve real equations
512  {
513  Scalar x = m_matT.coeff(i,i+1);
514  Scalar y = m_matT.coeff(i+1,i);
515  Scalar denom = (m_eivalues.coeff(i).real() - p) * (m_eivalues.coeff(i).real() - p) + m_eivalues.coeff(i).imag() * m_eivalues.coeff(i).imag();
516  Scalar t = (x * lastr - lastw * r) / denom;
517  m_matT.coeffRef(i,n) = t;
518  if (abs(x) > abs(lastw))
519  m_matT.coeffRef(i+1,n) = (-r - w * t) / x;
520  else
521  m_matT.coeffRef(i+1,n) = (-lastr - y * t) / lastw;
522  }
523 
524  // Overflow control
525  Scalar t = abs(m_matT.coeff(i,n));
526  if ((eps * t) * t > Scalar(1))
527  m_matT.col(n).tail(size-i) /= t;
528  }
529  }
530  }
531  else if (q < Scalar(0) && n > 0) // Complex vector
532  {
533  Scalar lastra(0), lastsa(0), lastw(0);
534  Index l = n-1;
535 
536  // Last vector component imaginary so matrix is triangular
537  if (abs(m_matT.coeff(n,n-1)) > abs(m_matT.coeff(n-1,n)))
538  {
539  m_matT.coeffRef(n-1,n-1) = q / m_matT.coeff(n,n-1);
540  m_matT.coeffRef(n-1,n) = -(m_matT.coeff(n,n) - p) / m_matT.coeff(n,n-1);
541  }
542  else
543  {
544  ComplexScalar cc = ComplexScalar(Scalar(0),-m_matT.coeff(n-1,n)) / ComplexScalar(m_matT.coeff(n-1,n-1)-p,q);
545  m_matT.coeffRef(n-1,n-1) = numext::real(cc);
546  m_matT.coeffRef(n-1,n) = numext::imag(cc);
547  }
548  m_matT.coeffRef(n,n-1) = Scalar(0);
549  m_matT.coeffRef(n,n) = Scalar(1);
550  for (Index i = n-2; i >= 0; i--)
551  {
552  Scalar ra = m_matT.row(i).segment(l, n-l+1).dot(m_matT.col(n-1).segment(l, n-l+1));
553  Scalar sa = m_matT.row(i).segment(l, n-l+1).dot(m_matT.col(n).segment(l, n-l+1));
554  Scalar w = m_matT.coeff(i,i) - p;
555 
556  if (m_eivalues.coeff(i).imag() < Scalar(0))
557  {
558  lastw = w;
559  lastra = ra;
560  lastsa = sa;
561  }
562  else
563  {
564  l = i;
565  if (m_eivalues.coeff(i).imag() == RealScalar(0))
566  {
567  ComplexScalar cc = ComplexScalar(-ra,-sa) / ComplexScalar(w,q);
568  m_matT.coeffRef(i,n-1) = numext::real(cc);
569  m_matT.coeffRef(i,n) = numext::imag(cc);
570  }
571  else
572  {
573  // Solve complex equations
574  Scalar x = m_matT.coeff(i,i+1);
575  Scalar y = m_matT.coeff(i+1,i);
576  Scalar vr = (m_eivalues.coeff(i).real() - p) * (m_eivalues.coeff(i).real() - p) + m_eivalues.coeff(i).imag() * m_eivalues.coeff(i).imag() - q * q;
577  Scalar vi = (m_eivalues.coeff(i).real() - p) * Scalar(2) * q;
578  if ((vr == Scalar(0)) && (vi == Scalar(0)))
579  vr = eps * norm * (abs(w) + abs(q) + abs(x) + abs(y) + abs(lastw));
580 
581  ComplexScalar cc = ComplexScalar(x*lastra-lastw*ra+q*sa,x*lastsa-lastw*sa-q*ra) / ComplexScalar(vr,vi);
582  m_matT.coeffRef(i,n-1) = numext::real(cc);
583  m_matT.coeffRef(i,n) = numext::imag(cc);
584  if (abs(x) > (abs(lastw) + abs(q)))
585  {
586  m_matT.coeffRef(i+1,n-1) = (-ra - w * m_matT.coeff(i,n-1) + q * m_matT.coeff(i,n)) / x;
587  m_matT.coeffRef(i+1,n) = (-sa - w * m_matT.coeff(i,n) - q * m_matT.coeff(i,n-1)) / x;
588  }
589  else
590  {
591  cc = ComplexScalar(-lastra-y*m_matT.coeff(i,n-1),-lastsa-y*m_matT.coeff(i,n)) / ComplexScalar(lastw,q);
592  m_matT.coeffRef(i+1,n-1) = numext::real(cc);
593  m_matT.coeffRef(i+1,n) = numext::imag(cc);
594  }
595  }
596 
597  // Overflow control
598  Scalar t = numext::maxi<Scalar>(abs(m_matT.coeff(i,n-1)),abs(m_matT.coeff(i,n)));
599  if ((eps * t) * t > Scalar(1))
600  m_matT.block(i, n-1, size-i, 2) /= t;
601 
602  }
603  }
604 
605  // We handled a pair of complex conjugate eigenvalues, so need to skip them both
606  n--;
607  }
608  else
609  {
610  eigen_assert(0 && "Internal bug in EigenSolver (INF or NaN has not been detected)"); // this should not happen
611  }
612  }
613 
614  // Back transformation to get eigenvectors of original matrix
615  for (Index j = size-1; j >= 0; j--)
616  {
617  m_tmp.noalias() = m_eivec.leftCols(j+1) * m_matT.col(j).segment(0, j+1);
618  m_eivec.col(j) = m_tmp;
619  }
620 }
621 
622 } // end namespace Eigen
623 
624 #endif // EIGEN_EIGENSOLVER_H
Computes eigenvalues and eigenvectors of general matrices.
Definition: EigenSolver.h:67
MatrixType::Scalar Scalar
Scalar type for matrices of type MatrixType.
Definition: EigenSolver.h:82
EigenSolver(const EigenBase< InputType > &matrix, bool computeEigenvectors=true)
Constructor; computes eigendecomposition of given matrix.
Definition: EigenSolver.h:149
const EigenvalueType & eigenvalues() const
Returns the eigenvalues of given matrix.
Definition: EigenSolver.h:246
EigenSolver(Index size)
Default constructor with memory preallocation.
Definition: EigenSolver.h:123
MatrixType pseudoEigenvalueMatrix() const
Returns the block-diagonal matrix in the pseudo-eigendecomposition.
Definition: EigenSolver.h:326
ComputationInfo info() const
Definition: EigenSolver.h:283
EigenSolver & setMaxIterations(Index maxIters)
Sets the maximum number of iterations allowed.
Definition: EigenSolver.h:290
EigenvectorsType eigenvectors() const
Returns the eigenvectors of given matrix.
Definition: EigenSolver.h:347
std::complex< RealScalar > ComplexScalar
Complex scalar type for MatrixType.
Definition: EigenSolver.h:92
EigenSolver()
Default constructor.
Definition: EigenSolver.h:115
Matrix< ComplexScalar, ColsAtCompileTime, 1, Options &~RowMajor, MaxColsAtCompileTime, 1 > EigenvalueType
Type for vector of eigenvalues as returned by eigenvalues().
Definition: EigenSolver.h:99
Eigen::Index Index
Definition: EigenSolver.h:84
EigenSolver & compute(const EigenBase< InputType > &matrix, bool computeEigenvectors=true)
Computes eigendecomposition of given matrix.
const MatrixType & pseudoEigenvectors() const
Returns the pseudo-eigenvectors of given matrix.
Definition: EigenSolver.h:201
Index getMaxIterations()
Returns the maximum number of iterations.
Definition: EigenSolver.h:297
Matrix< ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime > EigenvectorsType
Type for matrix of eigenvectors as returned by eigenvectors().
Definition: EigenSolver.h:106
MatrixType_ MatrixType
Synonym for the template parameter MatrixType_.
Definition: EigenSolver.h:71
Scalar & coeffRef(Index rowId, Index colId)
Definition: PlainObjectBase.h:187
RealSchur & setMaxIterations(Index maxIters)
Sets the maximum number of iterations allowed.
Definition: RealSchur.h:208
Index getMaxIterations()
Returns the maximum number of iterations.
Definition: RealSchur.h:215
ComputationInfo
Definition: Constants.h:442
@ NumericalIssue
Definition: Constants.h:446
@ Success
Definition: Constants.h:444
Namespace containing all symbols from the Eigen library.
Definition: Core:139
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:59
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_isfinite_op< typename Derived::Scalar >, const Derived > isfinite(const Eigen::ArrayBase< Derived > &x)
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_abs_op< typename Derived::Scalar >, const Derived > abs(const Eigen::ArrayBase< Derived > &x)
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_sqrt_op< typename Derived::Scalar >, const Derived > sqrt(const Eigen::ArrayBase< Derived > &x)
Definition: EigenBase.h:32
Derived & derived()
Definition: EigenBase.h:48
EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT
Definition: EigenBase.h:65
EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT
Definition: EigenBase.h:62
Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
Definition: NumTraits.h:231