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