EigenSolver.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 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_EIGENSOLVER_H
27 #define EIGEN_EIGENSOLVER_H
28 
29 #include "./RealSchur.h"
30 
31 namespace Eigen {
32 
79 template<typename _MatrixType> class EigenSolver
80 {
81  public:
82 
84  typedef _MatrixType MatrixType;
85 
86  enum {
91  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
92  };
93 
95  typedef typename MatrixType::Scalar Scalar;
97  typedef typename MatrixType::Index Index;
98 
105  typedef std::complex<RealScalar> ComplexScalar;
106 
113 
120 
129 
137  : m_eivec(size, size),
138  m_eivalues(size),
139  m_isInitialized(false),
140  m_eigenvectorsOk(false),
141  m_realSchur(size),
142  m_matT(size, size),
143  m_tmp(size)
144  {}
145 
161  EigenSolver(const MatrixType& matrix, bool computeEigenvectors = true)
162  : m_eivec(matrix.rows(), matrix.cols()),
163  m_eivalues(matrix.cols()),
164  m_isInitialized(false),
165  m_eigenvectorsOk(false),
166  m_realSchur(matrix.cols()),
167  m_matT(matrix.rows(), matrix.cols()),
168  m_tmp(matrix.cols())
169  {
170  compute(matrix, computeEigenvectors);
171  }
172 
194 
214  {
215  eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
216  eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
217  return m_eivec;
218  }
219 
239 
259  {
260  eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
261  return m_eivalues;
262  }
263 
291  EigenSolver& compute(const MatrixType& matrix, bool computeEigenvectors = true);
292 
294  {
295  eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
296  return m_realSchur.info();
297  }
298 
299  private:
300  void doComputeEigenvectors();
301 
302  protected:
309 
312 };
313 
314 template<typename MatrixType>
316 {
317  eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
318  Index n = m_eivalues.rows();
319  MatrixType matD = MatrixType::Zero(n,n);
320  for (Index i=0; i<n; ++i)
321  {
322  if (internal::isMuchSmallerThan(internal::imag(m_eivalues.coeff(i)), internal::real(m_eivalues.coeff(i))))
323  matD.coeffRef(i,i) = internal::real(m_eivalues.coeff(i));
324  else
325  {
326  matD.template block<2,2>(i,i) << internal::real(m_eivalues.coeff(i)), internal::imag(m_eivalues.coeff(i)),
327  -internal::imag(m_eivalues.coeff(i)), internal::real(m_eivalues.coeff(i));
328  ++i;
329  }
330  }
331  return matD;
332 }
333 
334 template<typename MatrixType>
336 {
337  eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
338  eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
339  Index n = m_eivec.cols();
340  EigenvectorsType matV(n,n);
341  for (Index j=0; j<n; ++j)
342  {
343  if (internal::isMuchSmallerThan(internal::imag(m_eivalues.coeff(j)), internal::real(m_eivalues.coeff(j))) || j+1==n)
344  {
345  // we have a real eigen value
346  matV.col(j) = m_eivec.col(j).template cast<ComplexScalar>();
347  matV.col(j).normalize();
348  }
349  else
350  {
351  // we have a pair of complex eigen values
352  for (Index i=0; i<n; ++i)
353  {
354  matV.coeffRef(i,j) = ComplexScalar(m_eivec.coeff(i,j), m_eivec.coeff(i,j+1));
355  matV.coeffRef(i,j+1) = ComplexScalar(m_eivec.coeff(i,j), -m_eivec.coeff(i,j+1));
356  }
357  matV.col(j).normalize();
358  matV.col(j+1).normalize();
359  ++j;
360  }
361  }
362  return matV;
363 }
364 
365 template<typename MatrixType>
366 EigenSolver<MatrixType>& EigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvectors)
367 {
368  assert(matrix.cols() == matrix.rows());
369 
370  // Reduce to real Schur form.
371  m_realSchur.compute(matrix, computeEigenvectors);
372  if (m_realSchur.info() == Success)
373  {
374  m_matT = m_realSchur.matrixT();
375  if (computeEigenvectors)
376  m_eivec = m_realSchur.matrixU();
377 
378  // Compute eigenvalues from matT
379  m_eivalues.resize(matrix.cols());
380  Index i = 0;
381  while (i < matrix.cols())
382  {
383  if (i == matrix.cols() - 1 || m_matT.coeff(i+1, i) == Scalar(0))
384  {
385  m_eivalues.coeffRef(i) = m_matT.coeff(i, i);
386  ++i;
387  }
388  else
389  {
390  Scalar p = Scalar(0.5) * (m_matT.coeff(i, i) - m_matT.coeff(i+1, i+1));
391  Scalar z = internal::sqrt(internal::abs(p * p + m_matT.coeff(i+1, i) * m_matT.coeff(i, i+1)));
392  m_eivalues.coeffRef(i) = ComplexScalar(m_matT.coeff(i+1, i+1) + p, z);
393  m_eivalues.coeffRef(i+1) = ComplexScalar(m_matT.coeff(i+1, i+1) + p, -z);
394  i += 2;
395  }
396  }
397 
398  // Compute eigenvectors.
399  if (computeEigenvectors)
400  doComputeEigenvectors();
401  }
402 
403  m_isInitialized = true;
404  m_eigenvectorsOk = computeEigenvectors;
405 
406  return *this;
407 }
408 
409 // Complex scalar division.
410 template<typename Scalar>
411 std::complex<Scalar> cdiv(Scalar xr, Scalar xi, Scalar yr, Scalar yi)
412 {
413  Scalar r,d;
414  if (internal::abs(yr) > internal::abs(yi))
415  {
416  r = yi/yr;
417  d = yr + r*yi;
418  return std::complex<Scalar>((xr + r*xi)/d, (xi - r*xr)/d);
419  }
420  else
421  {
422  r = yr/yi;
423  d = yi + r*yr;
424  return std::complex<Scalar>((r*xr + xi)/d, (r*xi - xr)/d);
425  }
426 }
427 
428 
429 template<typename MatrixType>
430 void EigenSolver<MatrixType>::doComputeEigenvectors()
431 {
432  const Index size = m_eivec.cols();
433  const Scalar eps = NumTraits<Scalar>::epsilon();
434 
435  // inefficient! this is already computed in RealSchur
436  Scalar norm(0);
437  for (Index j = 0; j < size; ++j)
438  {
439  norm += m_matT.row(j).segment((std::max)(j-1,Index(0)), size-(std::max)(j-1,Index(0))).cwiseAbs().sum();
440  }
441 
442  // Backsubstitute to find vectors of upper triangular form
443  if (norm == 0.0)
444  {
445  return;
446  }
447 
448  for (Index n = size-1; n >= 0; n--)
449  {
450  Scalar p = m_eivalues.coeff(n).real();
451  Scalar q = m_eivalues.coeff(n).imag();
452 
453  // Scalar vector
454  if (q == Scalar(0))
455  {
456  Scalar lastr(0), lastw(0);
457  Index l = n;
458 
459  m_matT.coeffRef(n,n) = 1.0;
460  for (Index i = n-1; i >= 0; i--)
461  {
462  Scalar w = m_matT.coeff(i,i) - p;
463  Scalar r = m_matT.row(i).segment(l,n-l+1).dot(m_matT.col(n).segment(l, n-l+1));
464 
465  if (m_eivalues.coeff(i).imag() < 0.0)
466  {
467  lastw = w;
468  lastr = r;
469  }
470  else
471  {
472  l = i;
473  if (m_eivalues.coeff(i).imag() == 0.0)
474  {
475  if (w != 0.0)
476  m_matT.coeffRef(i,n) = -r / w;
477  else
478  m_matT.coeffRef(i,n) = -r / (eps * norm);
479  }
480  else // Solve real equations
481  {
482  Scalar x = m_matT.coeff(i,i+1);
483  Scalar y = m_matT.coeff(i+1,i);
484  Scalar denom = (m_eivalues.coeff(i).real() - p) * (m_eivalues.coeff(i).real() - p) + m_eivalues.coeff(i).imag() * m_eivalues.coeff(i).imag();
485  Scalar t = (x * lastr - lastw * r) / denom;
486  m_matT.coeffRef(i,n) = t;
487  if (internal::abs(x) > internal::abs(lastw))
488  m_matT.coeffRef(i+1,n) = (-r - w * t) / x;
489  else
490  m_matT.coeffRef(i+1,n) = (-lastr - y * t) / lastw;
491  }
492 
493  // Overflow control
494  Scalar t = internal::abs(m_matT.coeff(i,n));
495  if ((eps * t) * t > Scalar(1))
496  m_matT.col(n).tail(size-i) /= t;
497  }
498  }
499  }
500  else if (q < Scalar(0) && n > 0) // Complex vector
501  {
502  Scalar lastra(0), lastsa(0), lastw(0);
503  Index l = n-1;
504 
505  // Last vector component imaginary so matrix is triangular
506  if (internal::abs(m_matT.coeff(n,n-1)) > internal::abs(m_matT.coeff(n-1,n)))
507  {
508  m_matT.coeffRef(n-1,n-1) = q / m_matT.coeff(n,n-1);
509  m_matT.coeffRef(n-1,n) = -(m_matT.coeff(n,n) - p) / m_matT.coeff(n,n-1);
510  }
511  else
512  {
513  std::complex<Scalar> cc = cdiv<Scalar>(0.0,-m_matT.coeff(n-1,n),m_matT.coeff(n-1,n-1)-p,q);
514  m_matT.coeffRef(n-1,n-1) = internal::real(cc);
515  m_matT.coeffRef(n-1,n) = internal::imag(cc);
516  }
517  m_matT.coeffRef(n,n-1) = 0.0;
518  m_matT.coeffRef(n,n) = 1.0;
519  for (Index i = n-2; i >= 0; i--)
520  {
521  Scalar ra = m_matT.row(i).segment(l, n-l+1).dot(m_matT.col(n-1).segment(l, n-l+1));
522  Scalar sa = m_matT.row(i).segment(l, n-l+1).dot(m_matT.col(n).segment(l, n-l+1));
523  Scalar w = m_matT.coeff(i,i) - p;
524 
525  if (m_eivalues.coeff(i).imag() < 0.0)
526  {
527  lastw = w;
528  lastra = ra;
529  lastsa = sa;
530  }
531  else
532  {
533  l = i;
534  if (m_eivalues.coeff(i).imag() == RealScalar(0))
535  {
536  std::complex<Scalar> cc = cdiv(-ra,-sa,w,q);
537  m_matT.coeffRef(i,n-1) = internal::real(cc);
538  m_matT.coeffRef(i,n) = internal::imag(cc);
539  }
540  else
541  {
542  // Solve complex equations
543  Scalar x = m_matT.coeff(i,i+1);
544  Scalar y = m_matT.coeff(i+1,i);
545  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;
546  Scalar vi = (m_eivalues.coeff(i).real() - p) * Scalar(2) * q;
547  if ((vr == 0.0) && (vi == 0.0))
548  vr = eps * norm * (internal::abs(w) + internal::abs(q) + internal::abs(x) + internal::abs(y) + internal::abs(lastw));
549 
550  std::complex<Scalar> cc = cdiv(x*lastra-lastw*ra+q*sa,x*lastsa-lastw*sa-q*ra,vr,vi);
551  m_matT.coeffRef(i,n-1) = internal::real(cc);
552  m_matT.coeffRef(i,n) = internal::imag(cc);
553  if (internal::abs(x) > (internal::abs(lastw) + internal::abs(q)))
554  {
555  m_matT.coeffRef(i+1,n-1) = (-ra - w * m_matT.coeff(i,n-1) + q * m_matT.coeff(i,n)) / x;
556  m_matT.coeffRef(i+1,n) = (-sa - w * m_matT.coeff(i,n) - q * m_matT.coeff(i,n-1)) / x;
557  }
558  else
559  {
560  cc = cdiv(-lastra-y*m_matT.coeff(i,n-1),-lastsa-y*m_matT.coeff(i,n),lastw,q);
561  m_matT.coeffRef(i+1,n-1) = internal::real(cc);
562  m_matT.coeffRef(i+1,n) = internal::imag(cc);
563  }
564  }
565 
566  // Overflow control
567  using std::max;
568  Scalar t = (max)(internal::abs(m_matT.coeff(i,n-1)),internal::abs(m_matT.coeff(i,n)));
569  if ((eps * t) * t > Scalar(1))
570  m_matT.block(i, n-1, size-i, 2) /= t;
571 
572  }
573  }
574 
575  // We handled a pair of complex conjugate eigenvalues, so need to skip them both
576  n--;
577  }
578  else
579  {
580  eigen_assert(0 && "Internal bug in EigenSolver"); // this should not happen
581  }
582  }
583 
584  // Back transformation to get eigenvectors of original matrix
585  for (Index j = size-1; j >= 0; j--)
586  {
587  m_tmp.noalias() = m_eivec.leftCols(j+1) * m_matT.col(j).segment(0, j+1);
588  m_eivec.col(j) = m_tmp;
589  }
590 }
591 
592 } // end namespace Eigen
593 
594 #endif // EIGEN_EIGENSOLVER_H