Eigen-Contrib  5.0.1
 
Loading...
Searching...
No Matches
Eigen::DPR1EigenSolver< RealScalar_ > Class Template Reference

#include <contrib/Eigen/src/StructuredMatrices/DPR1EigenSolver.h>

Detailed Description

template<typename RealScalar_>
class Eigen::DPR1EigenSolver< RealScalar_ >

Direct O(n^2) eigensolver for real symmetric diagonal-plus-rank-one matrices \( A = D + \rho\, z z^T \), via the secular equation.

This is the standalone version of the kernel at the heart of the divide-and-conquer symmetric eigensolvers (LAPACK's xLAED2/3/4): after deflation – entries with negligible \( |z_i| \) are eigenpairs of the diagonal already, and (nearly) equal diagonal entries are combined by Givens rotations whose dropped coupling is below a backward-stability threshold – the surviving eigenvalues are the roots of the secular equation

\[ f(\lambda) = 1 + \rho \sum_i \frac{z_i^2}{d_i - \lambda} = 0, \]

one in each interval between consecutive poles. Each root is bracketed and bisected in coordinates shifted to its nearest pole, so every distance \( \lambda - d_i \) is retained as an exact data difference plus a small offset instead of a cancellation-prone subtraction of close numbers. Eigenvectors are then built not from the original z but from the Gu-Eisenstat vector \( \hat z \) – the one for which the computed roots are exact secular eigenvalues – which is what makes the computed eigenvector matrix numerically orthogonal without any reorthogonalization.

The total cost is O(n^2): O(n log(1/eps)) per bisected root in the common case (the iteration cap is sized to the scalar's full exponent range, so even roots subnormally close to their pole resolve), O(n) per Gu-Eisenstat weight, and O(n) per eigenvector (the deflation rotations are replayed instead of accumulated into a dense matrix).

Both signs of \( \rho \) are supported (negative \( \rho \) is handled by negating the matrix), as are \( \rho = 0 \), zero z, repeated diagonal entries and any ordering of d. The problem is rescaled internally by the exact power of two that brings \( \max(\|D\|_\infty, |\rho| \|z\|^2) \) into [1/2, 1), which changes no bits of the result when \( \rho \|z\|^2 \) and the scaled poles are normal. InvalidInput is reported for non-finite input, for an update \( \rho \|z\|^2 \) that overflows, and for a computed eigenvalue that is not finite.

VectorXd lambda = es.eigenvalues(); // ascending
MatrixXd V = es.eigenvectors(); // orthogonal
Matrix< double, Dynamic, Dynamic > MatrixXd
Matrix< double, Dynamic, 1 > VectorXd
Template Parameters
RealScalar_one of float, double, or long double.

References:

  • M. Gu and S. C. Eisenstat, "A stable and efficient algorithm for the rank-one modification of the symmetric eigenproblem," SIAM J. Matrix Anal. Appl., 15(4):1266-1276, 1994.
  • P. H. Sterbenz, "Floating-Point Computation", Prentice-Hall, 1974. Scaling by a power of two is exact, the property the problem scaling relies on.
See also
class SelfAdjointEigenSolver

Public Member Functions

DPR1EigenSolver & compute (const VectorType &d, RealScalar rho, const VectorType &z, int options=ComputeEigenvectors)
 
 DPR1EigenSolver ()=default
 
 DPR1EigenSolver (const VectorType &d, RealScalar rho, const VectorType &z, int options=ComputeEigenvectors)
 
const VectorType & eigenvalues () const
 
const MatrixType & eigenvectors () const
 
ComputationInfo info () const
 

Constructor & Destructor Documentation

◆ DPR1EigenSolver() [1/2]

template<typename RealScalar_>
Eigen::DPR1EigenSolver< RealScalar_ >::DPR1EigenSolver ( )
default

Default constructor; call compute before querying results.

◆ DPR1EigenSolver() [2/2]

template<typename RealScalar_>
Eigen::DPR1EigenSolver< RealScalar_ >::DPR1EigenSolver ( const VectorType & d,
RealScalar rho,
const VectorType & z,
int options = ComputeEigenvectors )
inline

Computes the eigendecomposition of diag(d) + rho*z*z^T. options is ComputeEigenvectors (the default) or EigenvaluesOnly.

Member Function Documentation

◆ compute()

template<typename RealScalar_>
DPR1EigenSolver< RealScalar_ > & Eigen::DPR1EigenSolver< RealScalar_ >::compute ( const VectorType & d,
RealScalar rho,
const VectorType & z,
int options = ComputeEigenvectors )

Computes the eigendecomposition of diag(d) + rho*z*z^T.

See also
DPR1EigenSolver()

◆ eigenvalues()

template<typename RealScalar_>
const VectorType & Eigen::DPR1EigenSolver< RealScalar_ >::eigenvalues ( ) const
inline
Returns
the eigenvalues, sorted in increasing order.

◆ eigenvectors()

template<typename RealScalar_>
const MatrixType & Eigen::DPR1EigenSolver< RealScalar_ >::eigenvectors ( ) const
inline
Returns
the orthogonal matrix of eigenvectors; column k matches eigenvalues()[k].
Precondition
compute was called with ComputeEigenvectors.

◆ info()

template<typename RealScalar_>
ComputationInfo Eigen::DPR1EigenSolver< RealScalar_ >::info ( ) const
inline
Returns
Success if the decomposition succeeded, NoConvergence if a secular root could not be fully resolved, InvalidInput if the input was non-finite, \( \rho \|z\|^2 \) overflows, or a computed eigenvalue is not finite (the eigenvalues are then NaN).

The documentation for this class was generated from the following file: