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 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 namespace Eigen {
17 
64 template<typename _MatrixType> class EigenSolver
65 {
66  public:
67 
69  typedef _MatrixType MatrixType;
70 
71  enum {
72  RowsAtCompileTime = MatrixType::RowsAtCompileTime,
73  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
74  Options = MatrixType::Options,
75  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
76  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
77  };
78 
80  typedef typename MatrixType::Scalar Scalar;
81  typedef typename NumTraits<Scalar>::Real RealScalar;
82  typedef typename MatrixType::Index Index;
83 
90  typedef std::complex<RealScalar> ComplexScalar;
91 
97  typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> EigenvalueType;
98 
104  typedef Matrix<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime> EigenvectorsType;
105 
113  EigenSolver() : m_eivec(), m_eivalues(), m_isInitialized(false), m_realSchur(), m_matT(), m_tmp() {}
114 
121  EigenSolver(Index size)
122  : m_eivec(size, size),
123  m_eivalues(size),
124  m_isInitialized(false),
125  m_eigenvectorsOk(false),
126  m_realSchur(size),
127  m_matT(size, size),
128  m_tmp(size)
129  {}
130 
146  EigenSolver(const MatrixType& matrix, bool computeEigenvectors = true)
147  : m_eivec(matrix.rows(), matrix.cols()),
148  m_eivalues(matrix.cols()),
149  m_isInitialized(false),
150  m_eigenvectorsOk(false),
151  m_realSchur(matrix.cols()),
152  m_matT(matrix.rows(), matrix.cols()),
153  m_tmp(matrix.cols())
154  {
155  compute(matrix, computeEigenvectors);
156  }
157 
178  EigenvectorsType eigenvectors() const;
179 
198  const MatrixType& pseudoEigenvectors() const
199  {
200  eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
201  eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
202  return m_eivec;
203  }
204 
223  MatrixType pseudoEigenvalueMatrix() const;
224 
243  const EigenvalueType& eigenvalues() const
244  {
245  eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
246  return m_eivalues;
247  }
248 
276  EigenSolver& compute(const MatrixType& matrix, bool computeEigenvectors = true);
277 
278  ComputationInfo info() const
279  {
280  eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
281  return m_realSchur.info();
282  }
283 
284  private:
285  void doComputeEigenvectors();
286 
287  protected:
288  MatrixType m_eivec;
289  EigenvalueType m_eivalues;
290  bool m_isInitialized;
291  bool m_eigenvectorsOk;
292  RealSchur<MatrixType> m_realSchur;
293  MatrixType m_matT;
294 
295  typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ColumnVectorType;
296  ColumnVectorType m_tmp;
297 };
298 
299 template<typename MatrixType>
300 MatrixType EigenSolver<MatrixType>::pseudoEigenvalueMatrix() const
301 {
302  eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
303  Index n = m_eivalues.rows();
304  MatrixType matD = MatrixType::Zero(n,n);
305  for (Index i=0; i<n; ++i)
306  {
307  if (internal::isMuchSmallerThan(internal::imag(m_eivalues.coeff(i)), internal::real(m_eivalues.coeff(i))))
308  matD.coeffRef(i,i) = internal::real(m_eivalues.coeff(i));
309  else
310  {
311  matD.template block<2,2>(i,i) << internal::real(m_eivalues.coeff(i)), internal::imag(m_eivalues.coeff(i)),
312  -internal::imag(m_eivalues.coeff(i)), internal::real(m_eivalues.coeff(i));
313  ++i;
314  }
315  }
316  return matD;
317 }
318 
319 template<typename MatrixType>
320 typename EigenSolver<MatrixType>::EigenvectorsType EigenSolver<MatrixType>::eigenvectors() const
321 {
322  eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
323  eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
324  Index n = m_eivec.cols();
325  EigenvectorsType matV(n,n);
326  for (Index j=0; j<n; ++j)
327  {
328  if (internal::isMuchSmallerThan(internal::imag(m_eivalues.coeff(j)), internal::real(m_eivalues.coeff(j))) || j+1==n)
329  {
330  // we have a real eigen value
331  matV.col(j) = m_eivec.col(j).template cast<ComplexScalar>();
332  matV.col(j).normalize();
333  }
334  else
335  {
336  // we have a pair of complex eigen values
337  for (Index i=0; i<n; ++i)
338  {
339  matV.coeffRef(i,j) = ComplexScalar(m_eivec.coeff(i,j), m_eivec.coeff(i,j+1));
340  matV.coeffRef(i,j+1) = ComplexScalar(m_eivec.coeff(i,j), -m_eivec.coeff(i,j+1));
341  }
342  matV.col(j).normalize();
343  matV.col(j+1).normalize();
344  ++j;
345  }
346  }
347  return matV;
348 }
349 
350 template<typename MatrixType>
351 EigenSolver<MatrixType>& EigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvectors)
352 {
353  assert(matrix.cols() == matrix.rows());
354 
355  // Reduce to real Schur form.
356  m_realSchur.compute(matrix, computeEigenvectors);
357  if (m_realSchur.info() == Success)
358  {
359  m_matT = m_realSchur.matrixT();
360  if (computeEigenvectors)
361  m_eivec = m_realSchur.matrixU();
362 
363  // Compute eigenvalues from matT
364  m_eivalues.resize(matrix.cols());
365  Index i = 0;
366  while (i < matrix.cols())
367  {
368  if (i == matrix.cols() - 1 || m_matT.coeff(i+1, i) == Scalar(0))
369  {
370  m_eivalues.coeffRef(i) = m_matT.coeff(i, i);
371  ++i;
372  }
373  else
374  {
375  Scalar p = Scalar(0.5) * (m_matT.coeff(i, i) - m_matT.coeff(i+1, i+1));
376  Scalar z = internal::sqrt(internal::abs(p * p + m_matT.coeff(i+1, i) * m_matT.coeff(i, i+1)));
377  m_eivalues.coeffRef(i) = ComplexScalar(m_matT.coeff(i+1, i+1) + p, z);
378  m_eivalues.coeffRef(i+1) = ComplexScalar(m_matT.coeff(i+1, i+1) + p, -z);
379  i += 2;
380  }
381  }
382 
383  // Compute eigenvectors.
384  if (computeEigenvectors)
385  doComputeEigenvectors();
386  }
387 
388  m_isInitialized = true;
389  m_eigenvectorsOk = computeEigenvectors;
390 
391  return *this;
392 }
393 
394 // Complex scalar division.
395 template<typename Scalar>
396 std::complex<Scalar> cdiv(Scalar xr, Scalar xi, Scalar yr, Scalar yi)
397 {
398  Scalar r,d;
399  if (internal::abs(yr) > internal::abs(yi))
400  {
401  r = yi/yr;
402  d = yr + r*yi;
403  return std::complex<Scalar>((xr + r*xi)/d, (xi - r*xr)/d);
404  }
405  else
406  {
407  r = yr/yi;
408  d = yi + r*yr;
409  return std::complex<Scalar>((r*xr + xi)/d, (r*xi - xr)/d);
410  }
411 }
412 
413 
414 template<typename MatrixType>
415 void EigenSolver<MatrixType>::doComputeEigenvectors()
416 {
417  const Index size = m_eivec.cols();
418  const Scalar eps = NumTraits<Scalar>::epsilon();
419 
420  // inefficient! this is already computed in RealSchur
421  Scalar norm(0);
422  for (Index j = 0; j < size; ++j)
423  {
424  norm += m_matT.row(j).segment((std::max)(j-1,Index(0)), size-(std::max)(j-1,Index(0))).cwiseAbs().sum();
425  }
426 
427  // Backsubstitute to find vectors of upper triangular form
428  if (norm == 0.0)
429  {
430  return;
431  }
432 
433  for (Index n = size-1; n >= 0; n--)
434  {
435  Scalar p = m_eivalues.coeff(n).real();
436  Scalar q = m_eivalues.coeff(n).imag();
437 
438  // Scalar vector
439  if (q == Scalar(0))
440  {
441  Scalar lastr(0), lastw(0);
442  Index l = n;
443 
444  m_matT.coeffRef(n,n) = 1.0;
445  for (Index i = n-1; i >= 0; i--)
446  {
447  Scalar w = m_matT.coeff(i,i) - p;
448  Scalar r = m_matT.row(i).segment(l,n-l+1).dot(m_matT.col(n).segment(l, n-l+1));
449 
450  if (m_eivalues.coeff(i).imag() < 0.0)
451  {
452  lastw = w;
453  lastr = r;
454  }
455  else
456  {
457  l = i;
458  if (m_eivalues.coeff(i).imag() == 0.0)
459  {
460  if (w != 0.0)
461  m_matT.coeffRef(i,n) = -r / w;
462  else
463  m_matT.coeffRef(i,n) = -r / (eps * norm);
464  }
465  else // Solve real equations
466  {
467  Scalar x = m_matT.coeff(i,i+1);
468  Scalar y = m_matT.coeff(i+1,i);
469  Scalar denom = (m_eivalues.coeff(i).real() - p) * (m_eivalues.coeff(i).real() - p) + m_eivalues.coeff(i).imag() * m_eivalues.coeff(i).imag();
470  Scalar t = (x * lastr - lastw * r) / denom;
471  m_matT.coeffRef(i,n) = t;
472  if (internal::abs(x) > internal::abs(lastw))
473  m_matT.coeffRef(i+1,n) = (-r - w * t) / x;
474  else
475  m_matT.coeffRef(i+1,n) = (-lastr - y * t) / lastw;
476  }
477 
478  // Overflow control
479  Scalar t = internal::abs(m_matT.coeff(i,n));
480  if ((eps * t) * t > Scalar(1))
481  m_matT.col(n).tail(size-i) /= t;
482  }
483  }
484  }
485  else if (q < Scalar(0) && n > 0) // Complex vector
486  {
487  Scalar lastra(0), lastsa(0), lastw(0);
488  Index l = n-1;
489 
490  // Last vector component imaginary so matrix is triangular
491  if (internal::abs(m_matT.coeff(n,n-1)) > internal::abs(m_matT.coeff(n-1,n)))
492  {
493  m_matT.coeffRef(n-1,n-1) = q / m_matT.coeff(n,n-1);
494  m_matT.coeffRef(n-1,n) = -(m_matT.coeff(n,n) - p) / m_matT.coeff(n,n-1);
495  }
496  else
497  {
498  std::complex<Scalar> cc = cdiv<Scalar>(0.0,-m_matT.coeff(n-1,n),m_matT.coeff(n-1,n-1)-p,q);
499  m_matT.coeffRef(n-1,n-1) = internal::real(cc);
500  m_matT.coeffRef(n-1,n) = internal::imag(cc);
501  }
502  m_matT.coeffRef(n,n-1) = 0.0;
503  m_matT.coeffRef(n,n) = 1.0;
504  for (Index i = n-2; i >= 0; i--)
505  {
506  Scalar ra = m_matT.row(i).segment(l, n-l+1).dot(m_matT.col(n-1).segment(l, n-l+1));
507  Scalar sa = m_matT.row(i).segment(l, n-l+1).dot(m_matT.col(n).segment(l, n-l+1));
508  Scalar w = m_matT.coeff(i,i) - p;
509 
510  if (m_eivalues.coeff(i).imag() < 0.0)
511  {
512  lastw = w;
513  lastra = ra;
514  lastsa = sa;
515  }
516  else
517  {
518  l = i;
519  if (m_eivalues.coeff(i).imag() == RealScalar(0))
520  {
521  std::complex<Scalar> cc = cdiv(-ra,-sa,w,q);
522  m_matT.coeffRef(i,n-1) = internal::real(cc);
523  m_matT.coeffRef(i,n) = internal::imag(cc);
524  }
525  else
526  {
527  // Solve complex equations
528  Scalar x = m_matT.coeff(i,i+1);
529  Scalar y = m_matT.coeff(i+1,i);
530  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;
531  Scalar vi = (m_eivalues.coeff(i).real() - p) * Scalar(2) * q;
532  if ((vr == 0.0) && (vi == 0.0))
533  vr = eps * norm * (internal::abs(w) + internal::abs(q) + internal::abs(x) + internal::abs(y) + internal::abs(lastw));
534 
535  std::complex<Scalar> cc = cdiv(x*lastra-lastw*ra+q*sa,x*lastsa-lastw*sa-q*ra,vr,vi);
536  m_matT.coeffRef(i,n-1) = internal::real(cc);
537  m_matT.coeffRef(i,n) = internal::imag(cc);
538  if (internal::abs(x) > (internal::abs(lastw) + internal::abs(q)))
539  {
540  m_matT.coeffRef(i+1,n-1) = (-ra - w * m_matT.coeff(i,n-1) + q * m_matT.coeff(i,n)) / x;
541  m_matT.coeffRef(i+1,n) = (-sa - w * m_matT.coeff(i,n) - q * m_matT.coeff(i,n-1)) / x;
542  }
543  else
544  {
545  cc = cdiv(-lastra-y*m_matT.coeff(i,n-1),-lastsa-y*m_matT.coeff(i,n),lastw,q);
546  m_matT.coeffRef(i+1,n-1) = internal::real(cc);
547  m_matT.coeffRef(i+1,n) = internal::imag(cc);
548  }
549  }
550 
551  // Overflow control
552  using std::max;
553  Scalar t = (max)(internal::abs(m_matT.coeff(i,n-1)),internal::abs(m_matT.coeff(i,n)));
554  if ((eps * t) * t > Scalar(1))
555  m_matT.block(i, n-1, size-i, 2) /= t;
556 
557  }
558  }
559 
560  // We handled a pair of complex conjugate eigenvalues, so need to skip them both
561  n--;
562  }
563  else
564  {
565  eigen_assert(0 && "Internal bug in EigenSolver"); // this should not happen
566  }
567  }
568 
569  // Back transformation to get eigenvectors of original matrix
570  for (Index j = size-1; j >= 0; j--)
571  {
572  m_tmp.noalias() = m_eivec.leftCols(j+1) * m_matT.col(j).segment(0, j+1);
573  m_eivec.col(j) = m_tmp;
574  }
575 }
576 
577 } // end namespace Eigen
578 
579 #endif // EIGEN_EIGENSOLVER_H