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
16namespace Eigen {
17
64template<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
299template<typename MatrixType>
300MatrixType 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
319template<typename MatrixType>
320typename 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
350template<typename MatrixType>
351EigenSolver<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.
395template<typename Scalar>
396std::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
414template<typename MatrixType>
415void 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
Eigen::EigenSolver
Computes eigenvalues and eigenvectors of general matrices.
Definition EigenSolver.h:65
Eigen::EigenSolver::compute
EigenSolver & compute(const MatrixType &matrix, bool computeEigenvectors=true)
Computes eigendecomposition of given matrix.
Definition EigenSolver.h:351
Eigen::EigenSolver::eigenvectors
EigenvectorsType eigenvectors() const
Returns the eigenvectors of given matrix.
Definition EigenSolver.h:320
Eigen::EigenSolver::ComplexScalar
std::complex< RealScalar > ComplexScalar
Complex scalar type for MatrixType.
Definition EigenSolver.h:90
Eigen::EigenSolver::Scalar
MatrixType::Scalar Scalar
Scalar type for matrices of type MatrixType.
Definition EigenSolver.h:80
Eigen::EigenSolver::EigenvectorsType
Matrix< ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime > EigenvectorsType
Type for matrix of eigenvectors as returned by eigenvectors().
Definition EigenSolver.h:104
Eigen::EigenSolver::EigenvalueType
Matrix< ComplexScalar, ColsAtCompileTime, 1, Options &~RowMajor, MaxColsAtCompileTime, 1 > EigenvalueType
Type for vector of eigenvalues as returned by eigenvalues().
Definition EigenSolver.h:97
Eigen::EigenSolver::EigenSolver
EigenSolver()
Default constructor.
Definition EigenSolver.h:113
Eigen::EigenSolver::EigenSolver
EigenSolver(const MatrixType &matrix, bool computeEigenvectors=true)
Constructor; computes eigendecomposition of given matrix.
Definition EigenSolver.h:146
Eigen::EigenSolver::MatrixType
_MatrixType MatrixType
Synonym for the template parameter _MatrixType.
Definition EigenSolver.h:69
Eigen::EigenSolver::pseudoEigenvectors
const MatrixType & pseudoEigenvectors() const
Returns the pseudo-eigenvectors of given matrix.
Definition EigenSolver.h:198
Eigen::EigenSolver::pseudoEigenvalueMatrix
MatrixType pseudoEigenvalueMatrix() const
Returns the block-diagonal matrix in the pseudo-eigendecomposition.
Definition EigenSolver.h:300
Eigen::EigenSolver::EigenSolver
EigenSolver(Index size)
Default constructor with memory preallocation.
Definition EigenSolver.h:121
Eigen::EigenSolver::eigenvalues
const EigenvalueType & eigenvalues() const
Returns the eigenvalues of given matrix.
Definition EigenSolver.h:243
Eigen::Matrix
The matrix class, also used for vectors and row-vectors.
Definition Matrix.h:129
Eigen::RealSchur::info
ComputationInfo info() const
Reports whether previous computation was successful.
Definition RealSchur.h:170
Eigen::ComputationInfo
ComputationInfo
Definition Constants.h:367
Eigen::Success
@ Success
Definition Constants.h:369