SelfAdjointEigenSolver.h
Go to the documentation of this file.
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/>.
25 
26 #ifndef EIGEN_SELFADJOINTEIGENSOLVER_H
27 #define EIGEN_SELFADJOINTEIGENSOLVER_H
28 
29 #include "./Tridiagonalization.h"
30 
31 namespace Eigen {
32 
33 template<typename _MatrixType>
34 class GeneralizedSelfAdjointEigenSolver;
35 
36 namespace internal {
37 template<typename SolverType,int Size,bool IsComplex> struct direct_selfadjoint_eigenvalues;
38 }
39 
83 template<typename _MatrixType> class SelfAdjointEigenSolver
84 {
85  public:
86 
87  typedef _MatrixType MatrixType;
88  enum {
89  Size = MatrixType::RowsAtCompileTime,
92  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
93  };
94 
96  typedef typename MatrixType::Scalar Scalar;
97  typedef typename MatrixType::Index Index;
98 
106 
107  friend struct internal::direct_selfadjoint_eigenvalues<SelfAdjointEigenSolver,Size,NumTraits<Scalar>::IsComplex>;
108 
114  typedef typename internal::plain_col_type<MatrixType, RealScalar>::type RealVectorType;
116 
128  : m_eivec(),
129  m_eivalues(),
130  m_subdiag(),
131  m_isInitialized(false)
132  { }
133 
147  : m_eivec(size, size),
148  m_eivalues(size),
149  m_subdiag(size > 1 ? size - 1 : 1),
150  m_isInitialized(false)
151  {}
152 
169  : m_eivec(matrix.rows(), matrix.cols()),
170  m_eivalues(matrix.cols()),
171  m_subdiag(matrix.rows() > 1 ? matrix.rows() - 1 : 1),
172  m_isInitialized(false)
173  {
174  compute(matrix, options);
175  }
176 
207  SelfAdjointEigenSolver& compute(const MatrixType& matrix, int options = ComputeEigenvectors);
208 
224 
243  const MatrixType& eigenvectors() const
244  {
245  eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
246  eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
247  return m_eivec;
248  }
249 
266  {
267  eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
268  return m_eivalues;
269  }
270 
290  {
291  eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
292  eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
293  return m_eivec * m_eivalues.cwiseSqrt().asDiagonal() * m_eivec.adjoint();
294  }
295 
315  {
316  eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
317  eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
318  return m_eivec * m_eivalues.cwiseInverse().cwiseSqrt().asDiagonal() * m_eivec.adjoint();
319  }
320 
326  {
327  eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
328  return m_info;
329  }
330 
336  static const int m_maxIterations = 30;
337 
338  #ifdef EIGEN2_SUPPORT
339  SelfAdjointEigenSolver(const MatrixType& matrix, bool computeEigenvectors)
340  : m_eivec(matrix.rows(), matrix.cols()),
341  m_eivalues(matrix.cols()),
342  m_subdiag(matrix.rows() > 1 ? matrix.rows() - 1 : 1),
343  m_isInitialized(false)
344  {
345  compute(matrix, computeEigenvectors);
346  }
347 
348  SelfAdjointEigenSolver(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors = true)
349  : m_eivec(matA.cols(), matA.cols()),
350  m_eivalues(matA.cols()),
351  m_subdiag(matA.cols() > 1 ? matA.cols() - 1 : 1),
352  m_isInitialized(false)
353  {
354  static_cast<GeneralizedSelfAdjointEigenSolver<MatrixType>*>(this)->compute(matA, matB, computeEigenvectors ? ComputeEigenvectors : EigenvaluesOnly);
355  }
356 
357  void compute(const MatrixType& matrix, bool computeEigenvectors)
358  {
359  compute(matrix, computeEigenvectors ? ComputeEigenvectors : EigenvaluesOnly);
360  }
361 
362  void compute(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors = true)
363  {
364  compute(matA, matB, computeEigenvectors ? ComputeEigenvectors : EigenvaluesOnly);
365  }
366  #endif // EIGEN2_SUPPORT
367 
368  protected:
375 };
376 
393 namespace internal {
394 template<int StorageOrder,typename RealScalar, typename Scalar, typename Index>
395 static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index start, Index end, Scalar* matrixQ, Index n);
396 }
397 
398 template<typename MatrixType>
399 SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>
400 ::compute(const MatrixType& matrix, int options)
401 {
402  eigen_assert(matrix.cols() == matrix.rows());
403  eigen_assert((options&~(EigVecMask|GenEigMask))==0
404  && (options&EigVecMask)!=EigVecMask
405  && "invalid option parameter");
406  bool computeEigenvectors = (options&ComputeEigenvectors)==ComputeEigenvectors;
407  Index n = matrix.cols();
408  m_eivalues.resize(n,1);
409 
410  if(n==1)
411  {
412  m_eivalues.coeffRef(0,0) = internal::real(matrix.coeff(0,0));
413  if(computeEigenvectors)
414  m_eivec.setOnes(n,n);
415  m_info = Success;
416  m_isInitialized = true;
417  m_eigenvectorsOk = computeEigenvectors;
418  return *this;
419  }
420 
421  // declare some aliases
422  RealVectorType& diag = m_eivalues;
423  MatrixType& mat = m_eivec;
424 
425  // map the matrix coefficients to [-1:1] to avoid over- and underflow.
426  RealScalar scale = matrix.cwiseAbs().maxCoeff();
427  if(scale==RealScalar(0)) scale = RealScalar(1);
428  mat = matrix / scale;
429  m_subdiag.resize(n-1);
430  internal::tridiagonalization_inplace(mat, diag, m_subdiag, computeEigenvectors);
431 
432  Index end = n-1;
433  Index start = 0;
434  Index iter = 0; // total number of iterations
435 
436  while (end>0)
437  {
438  for (Index i = start; i<end; ++i)
439  if (internal::isMuchSmallerThan(internal::abs(m_subdiag[i]),(internal::abs(diag[i])+internal::abs(diag[i+1]))))
440  m_subdiag[i] = 0;
441 
442  // find the largest unreduced block
443  while (end>0 && m_subdiag[end-1]==0)
444  {
445  end--;
446  }
447  if (end<=0)
448  break;
449 
450  // if we spent too many iterations, we give up
451  iter++;
452  if(iter > m_maxIterations * n) break;
453 
454  start = end - 1;
455  while (start>0 && m_subdiag[start-1]!=0)
456  start--;
457 
458  internal::tridiagonal_qr_step<MatrixType::Flags&RowMajorBit ? RowMajor : ColMajor>(diag.data(), m_subdiag.data(), start, end, computeEigenvectors ? m_eivec.data() : (Scalar*)0, n);
459  }
460 
461  if (iter <= m_maxIterations * n)
462  m_info = Success;
463  else
464  m_info = NoConvergence;
465 
466  // Sort eigenvalues and corresponding vectors.
467  // TODO make the sort optional ?
468  // TODO use a better sort algorithm !!
469  if (m_info == Success)
470  {
471  for (Index i = 0; i < n-1; ++i)
472  {
473  Index k;
474  m_eivalues.segment(i,n-i).minCoeff(&k);
475  if (k > 0)
476  {
477  std::swap(m_eivalues[i], m_eivalues[k+i]);
478  if(computeEigenvectors)
479  m_eivec.col(i).swap(m_eivec.col(k+i));
480  }
481  }
482  }
483 
484  // scale back the eigen values
485  m_eivalues *= scale;
486 
487  m_isInitialized = true;
488  m_eigenvectorsOk = computeEigenvectors;
489  return *this;
490 }
491 
492 
493 namespace internal {
494 
495 template<typename SolverType,int Size,bool IsComplex> struct direct_selfadjoint_eigenvalues
496 {
497  static inline void run(SolverType& eig, const typename SolverType::MatrixType& A, int options)
498  { eig.compute(A,options); }
499 };
500 
501 template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,3,false>
502 {
503  typedef typename SolverType::MatrixType MatrixType;
504  typedef typename SolverType::RealVectorType VectorType;
505  typedef typename SolverType::Scalar Scalar;
506 
507  static inline void computeRoots(const MatrixType& m, VectorType& roots)
508  {
509  using std::sqrt;
510  using std::atan2;
511  using std::cos;
512  using std::sin;
513  const Scalar s_inv3 = Scalar(1.0)/Scalar(3.0);
514  const Scalar s_sqrt3 = sqrt(Scalar(3.0));
515 
516  // The characteristic equation is x^3 - c2*x^2 + c1*x - c0 = 0. The
517  // eigenvalues are the roots to this equation, all guaranteed to be
518  // real-valued, because the matrix is symmetric.
519  Scalar c0 = m(0,0)*m(1,1)*m(2,2) + Scalar(2)*m(1,0)*m(2,0)*m(2,1) - m(0,0)*m(2,1)*m(2,1) - m(1,1)*m(2,0)*m(2,0) - m(2,2)*m(1,0)*m(1,0);
520  Scalar c1 = m(0,0)*m(1,1) - m(1,0)*m(1,0) + m(0,0)*m(2,2) - m(2,0)*m(2,0) + m(1,1)*m(2,2) - m(2,1)*m(2,1);
521  Scalar c2 = m(0,0) + m(1,1) + m(2,2);
522 
523  // Construct the parameters used in classifying the roots of the equation
524  // and in solving the equation for the roots in closed form.
525  Scalar c2_over_3 = c2*s_inv3;
526  Scalar a_over_3 = (c1 - c2*c2_over_3)*s_inv3;
527  if (a_over_3 > Scalar(0))
528  a_over_3 = Scalar(0);
529 
530  Scalar half_b = Scalar(0.5)*(c0 + c2_over_3*(Scalar(2)*c2_over_3*c2_over_3 - c1));
531 
532  Scalar q = half_b*half_b + a_over_3*a_over_3*a_over_3;
533  if (q > Scalar(0))
534  q = Scalar(0);
535 
536  // Compute the eigenvalues by solving for the roots of the polynomial.
537  Scalar rho = sqrt(-a_over_3);
538  Scalar theta = atan2(sqrt(-q),half_b)*s_inv3;
539  Scalar cos_theta = cos(theta);
540  Scalar sin_theta = sin(theta);
541  roots(0) = c2_over_3 + Scalar(2)*rho*cos_theta;
542  roots(1) = c2_over_3 - rho*(cos_theta + s_sqrt3*sin_theta);
543  roots(2) = c2_over_3 - rho*(cos_theta - s_sqrt3*sin_theta);
544 
545  // Sort in increasing order.
546  if (roots(0) >= roots(1))
547  std::swap(roots(0),roots(1));
548  if (roots(1) >= roots(2))
549  {
550  std::swap(roots(1),roots(2));
551  if (roots(0) >= roots(1))
552  std::swap(roots(0),roots(1));
553  }
554  }
555 
556  static inline void run(SolverType& solver, const MatrixType& mat, int options)
557  {
558  using std::sqrt;
559  eigen_assert(mat.cols() == 3 && mat.cols() == mat.rows());
560  eigen_assert((options&~(EigVecMask|GenEigMask))==0
561  && (options&EigVecMask)!=EigVecMask
562  && "invalid option parameter");
563  bool computeEigenvectors = (options&ComputeEigenvectors)==ComputeEigenvectors;
564 
565  MatrixType& eivecs = solver.m_eivec;
566  VectorType& eivals = solver.m_eivalues;
567 
568  // map the matrix coefficients to [-1:1] to avoid over- and underflow.
569  Scalar scale = mat.cwiseAbs().maxCoeff();
570  MatrixType scaledMat = mat / scale;
571 
572  // compute the eigenvalues
573  computeRoots(scaledMat,eivals);
574 
575  // compute the eigen vectors
576  if(computeEigenvectors)
577  {
578  Scalar safeNorm2 = Eigen::NumTraits<Scalar>::epsilon();
579  safeNorm2 *= safeNorm2;
580  if((eivals(2)-eivals(0))<=Eigen::NumTraits<Scalar>::epsilon())
581  {
582  eivecs.setIdentity();
583  }
584  else
585  {
586  scaledMat = scaledMat.template selfadjointView<Lower>();
587  MatrixType tmp;
588  tmp = scaledMat;
589 
590  Scalar d0 = eivals(2) - eivals(1);
591  Scalar d1 = eivals(1) - eivals(0);
592  int k = d0 > d1 ? 2 : 0;
593  d0 = d0 > d1 ? d1 : d0;
594 
595  tmp.diagonal().array () -= eivals(k);
596  VectorType cross;
597  Scalar n;
598  n = (cross = tmp.row(0).cross(tmp.row(1))).squaredNorm();
599 
600  if(n>safeNorm2)
601  eivecs.col(k) = cross / sqrt(n);
602  else
603  {
604  n = (cross = tmp.row(0).cross(tmp.row(2))).squaredNorm();
605 
606  if(n>safeNorm2)
607  eivecs.col(k) = cross / sqrt(n);
608  else
609  {
610  n = (cross = tmp.row(1).cross(tmp.row(2))).squaredNorm();
611 
612  if(n>safeNorm2)
613  eivecs.col(k) = cross / sqrt(n);
614  else
615  {
616  // the input matrix and/or the eigenvaues probably contains some inf/NaN,
617  // => exit
618  // scale back to the original size.
619  eivals *= scale;
620 
621  solver.m_info = NumericalIssue;
622  solver.m_isInitialized = true;
623  solver.m_eigenvectorsOk = computeEigenvectors;
624  return;
625  }
626  }
627  }
628 
629  tmp = scaledMat;
630  tmp.diagonal().array() -= eivals(1);
631 
633  eivecs.col(1) = eivecs.col(k).unitOrthogonal();
634  else
635  {
636  n = (cross = eivecs.col(k).cross(tmp.row(0).normalized())).squaredNorm();
637  if(n>safeNorm2)
638  eivecs.col(1) = cross / sqrt(n);
639  else
640  {
641  n = (cross = eivecs.col(k).cross(tmp.row(1))).squaredNorm();
642  if(n>safeNorm2)
643  eivecs.col(1) = cross / sqrt(n);
644  else
645  {
646  n = (cross = eivecs.col(k).cross(tmp.row(2))).squaredNorm();
647  if(n>safeNorm2)
648  eivecs.col(1) = cross / sqrt(n);
649  else
650  {
651  // we should never reach this point,
652  // if so the last two eigenvalues are likely to ve very closed to each other
653  eivecs.col(1) = eivecs.col(k).unitOrthogonal();
654  }
655  }
656  }
657 
658  // make sure that eivecs[1] is orthogonal to eivecs[2]
659  Scalar d = eivecs.col(1).dot(eivecs.col(k));
660  eivecs.col(1) = (eivecs.col(1) - d * eivecs.col(k)).normalized();
661  }
662 
663  eivecs.col(k==2 ? 0 : 2) = eivecs.col(k).cross(eivecs.col(1)).normalized();
664  }
665  }
666  // Rescale back to the original size.
667  eivals *= scale;
668 
669  solver.m_info = Success;
670  solver.m_isInitialized = true;
671  solver.m_eigenvectorsOk = computeEigenvectors;
672  }
673 };
674 
675 // 2x2 direct eigenvalues decomposition, code from Hauke Heibel
676 template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,2,false>
677 {
678  typedef typename SolverType::MatrixType MatrixType;
679  typedef typename SolverType::RealVectorType VectorType;
680  typedef typename SolverType::Scalar Scalar;
681 
682  static inline void computeRoots(const MatrixType& m, VectorType& roots)
683  {
684  using std::sqrt;
685  const Scalar t0 = Scalar(0.5) * sqrt( abs2(m(0,0)-m(1,1)) + Scalar(4)*m(1,0)*m(1,0));
686  const Scalar t1 = Scalar(0.5) * (m(0,0) + m(1,1));
687  roots(0) = t1 - t0;
688  roots(1) = t1 + t0;
689  }
690 
691  static inline void run(SolverType& solver, const MatrixType& mat, int options)
692  {
693  eigen_assert(mat.cols() == 2 && mat.cols() == mat.rows());
694  eigen_assert((options&~(EigVecMask|GenEigMask))==0
695  && (options&EigVecMask)!=EigVecMask
696  && "invalid option parameter");
697  bool computeEigenvectors = (options&ComputeEigenvectors)==ComputeEigenvectors;
698 
699  MatrixType& eivecs = solver.m_eivec;
700  VectorType& eivals = solver.m_eivalues;
701 
702  // map the matrix coefficients to [-1:1] to avoid over- and underflow.
703  Scalar scale = mat.cwiseAbs().maxCoeff();
704  scale = (std::max)(scale,Scalar(1));
705  MatrixType scaledMat = mat / scale;
706 
707  // Compute the eigenvalues
708  computeRoots(scaledMat,eivals);
709 
710  // compute the eigen vectors
711  if(computeEigenvectors)
712  {
713  scaledMat.diagonal().array () -= eivals(1);
714  Scalar a2 = abs2(scaledMat(0,0));
715  Scalar c2 = abs2(scaledMat(1,1));
716  Scalar b2 = abs2(scaledMat(1,0));
717  if(a2>c2)
718  {
719  eivecs.col(1) << -scaledMat(1,0), scaledMat(0,0);
720  eivecs.col(1) /= sqrt(a2+b2);
721  }
722  else
723  {
724  eivecs.col(1) << -scaledMat(1,1), scaledMat(1,0);
725  eivecs.col(1) /= sqrt(c2+b2);
726  }
727 
728  eivecs.col(0) << eivecs.col(1).unitOrthogonal();
729  }
730 
731  // Rescale back to the original size.
732  eivals *= scale;
733 
734  solver.m_info = Success;
735  solver.m_isInitialized = true;
736  solver.m_eigenvectorsOk = computeEigenvectors;
737  }
738 };
739 
740 }
741 
742 template<typename MatrixType>
743 SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>
744 ::computeDirect(const MatrixType& matrix, int options)
745 {
747  return *this;
748 }
749 
750 namespace internal {
751 template<int StorageOrder,typename RealScalar, typename Scalar, typename Index>
752 static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index start, Index end, Scalar* matrixQ, Index n)
753 {
754  // NOTE this version avoids over & underflow, however since the matrix is prescaled, overflow cannot occur,
755  // and underflows should be meaningless anyway. So I don't any reason to enable this version, but I keep
756  // it here for reference:
757 // RealScalar td = (diag[end-1] - diag[end])*RealScalar(0.5);
758 // RealScalar e = subdiag[end-1];
759 // RealScalar mu = diag[end] - (e / (td + (td>0 ? 1 : -1))) * (e / hypot(td,e));
760  RealScalar td = (diag[end-1] - diag[end])*RealScalar(0.5);
761  RealScalar e2 = abs2(subdiag[end-1]);
762  RealScalar mu = diag[end] - e2 / (td + (td>0 ? 1 : -1) * sqrt(td*td + e2));
763  RealScalar x = diag[start] - mu;
764  RealScalar z = subdiag[start];
765  for (Index k = start; k < end; ++k)
766  {
768  rot.makeGivens(x, z);
769 
770  // do T = G' T G
771  RealScalar sdk = rot.s() * diag[k] + rot.c() * subdiag[k];
772  RealScalar dkp1 = rot.s() * subdiag[k] + rot.c() * diag[k+1];
773 
774  diag[k] = rot.c() * (rot.c() * diag[k] - rot.s() * subdiag[k]) - rot.s() * (rot.c() * subdiag[k] - rot.s() * diag[k+1]);
775  diag[k+1] = rot.s() * sdk + rot.c() * dkp1;
776  subdiag[k] = rot.c() * sdk - rot.s() * dkp1;
777 
778 
779  if (k > start)
780  subdiag[k - 1] = rot.c() * subdiag[k-1] - rot.s() * z;
781 
782  x = subdiag[k];
783 
784  if (k < end - 1)
785  {
786  z = -rot.s() * subdiag[k+1];
787  subdiag[k + 1] = rot.c() * subdiag[k+1];
788  }
789 
790  // apply the givens rotation to the unit matrix Q = Q * G
791  if (matrixQ)
792  {
793  // FIXME if StorageOrder == RowMajor this operation is not very efficient
794  Map<Matrix<Scalar,Dynamic,Dynamic,StorageOrder> > q(matrixQ,n,n);
795  q.applyOnTheRight(k,k+1,rot);
796  }
797  }
798 }
799 
800 } // end namespace internal
801 
802 } // end namespace Eigen
803 
804 #endif // EIGEN_SELFADJOINTEIGENSOLVER_H