STK::SymEigen< SquareArray > Class Template Reference

The class SymEigen compute the eigenvalue Decomposition of a symmetric ArrayXX. More...

#include <STK_SymEigen.h>

Inheritance diagram for STK::SymEigen< SquareArray >:

Public Types
typedef ISymEigen< SymEigen< SquareArray > >	Base

Public Member Functions
	SymEigen ()
	Default Constructor.
	SymEigen (SquareArray const &data, bool ref=false)
	Constructor.
template<class Derived >
	SymEigen (ExprBase< Derived > const &data)
	constructor.
	SymEigen (SymEigen const &S)
	Copy constructor.
virtual	~SymEigen ()
	virtual destructor
bool	runImpl ()
	Diagonalization of eigenVectors_.
SymEigen &	operator= (const SymEigen &S)
	Operator = : overwrite the SymEigen with S.
	SymEigen (CSquareX const &data, bool ref)
Public Member Functions inherited from STK::ISymEigen< SymEigen< SquareArray > >
	~ISymEigen ()
	virtual destructor
ISymEigen &	operator= (ISymEigen const &eigen)
	Operator = : overwrite the ISymEigen with eigen.
Range const &	range () const
Type const &	norm () const
int const &	rank () const
Type const &	det () const
Type const &	trace () const
CArraySquare< Type, size_ > const &	rotation () const
CArraySquare< Type, size_ > const &	eigenVectors () const
CArrayVector< Type, size_ > const &	eigenValues () const
void	setData (ExprBase< OtherDerived > const &data)
	overloading of setData.
ArraySquare &	ginv (ArraySquare &res) const
	Compute the generalized inverse of the symmetric matrix and put the result in res.
ArraySquare &	ginvsqrt (ArraySquare &res) const
	Compute the generalized square root inverse of the symmetric matrix and put the result in res.
ArraySquare &	gsqrt (ArraySquare &res) const
	Compute the square root of the symmetric matrix and put the result in res.
virtual bool	run ()
	Find the eigenvalues and eigenvectors of the matrix.
Public Member Functions inherited from STK::IRunnerBase
String const &	error () const
	get the last error message.
Public Member Functions inherited from STK::IRecursiveTemplate< Derived >
Derived &	asDerived ()
	static cast : return a reference of this with a cast to the derived class.
Derived const &	asDerived () const
	static cast : return a const reference of this with a cast to the derived class.
Derived *	asPtrDerived ()
	static cast : return a ptr on a Derived of this with a cast to the derived class.
Derived const *	asPtrDerived () const
	static cast : return a ptr on a constant Derived of this with a cast to the derived class.
Derived *	clone () const
	create a leaf using the copy constructor of the Derived class.
Derived *	clone (bool isRef) const
	create a leaf using the copy constructor of the Derived class and a flag determining if the clone is a reference or not.

Private Member Functions
void	tridiagonalize ()
	compute the tri-diagonalization of eigenVectors_
void	compHouse ()
	compute the Householder matrix and P
void	diagonalize ()
	computing the diagonalization of eigenValues_ and F_

Private Attributes
VectorX	F_
	Temporary vector.

Additional Inherited Members
Protected Types inherited from STK::ISymEigen< SymEigen< SquareArray > >
enum
typedef IRunnerBase	Base
typedef hidden::AlgebraTraits< SymEigen< SquareArray > >::SquareArray	SquareArray
typedef hidden::Traits< SquareArray >::Type	Type
Protected Member Functions inherited from STK::ISymEigen< SymEigen< SquareArray > >
	ISymEigen ()
	Default constructor.
	ISymEigen (SquareArray const &data, bool ref=false)
	Constructor The original data set can be overwritten by the eigenvectors if it is stored in a CSquareXd.
	ISymEigen (ExprBase< OtherDerived > const &data)
	template constructor
	ISymEigen (ISymEigen const &eigen)
	Copy constructor.
void	finalizeStep ()
	finalize the computation by computing the trace, rank, trace norm and determinant of the matrix.
Protected Member Functions inherited from STK::IRunnerBase
	IRunnerBase ()
	default constructor
	IRunnerBase (IRunnerBase const &runner)
	copy constructor
virtual	~IRunnerBase ()
	destructor
virtual void	update ()
	update the runner.
Protected Member Functions inherited from STK::IRecursiveTemplate< Derived >
	IRecursiveTemplate ()
	constructor.
	~IRecursiveTemplate ()
	destructor.
Protected Attributes inherited from STK::ISymEigen< SymEigen< SquareArray > >
Range	range_
	range of the original data set.
Type	trace_
	trace norm
Type	norm_
	trace norm
int	rank_
	rank
Type	det_
	determinant
CArraySquare< Type, size_ >	eigenVectors_
	Square matrix or the eigenvectors.
CArrayVector< Type, size_ >	eigenValues_
	Array of the eigenvalues.
CVectorXi	SupportEigenVectors_
	Array for the support of the eigenvectors.
Protected Attributes inherited from STK::IRunnerBase
String	msg_error_
	String with the last error message.
bool	hasRun_
	true if run has been used, false otherwise

Detailed Description

template<class SquareArray>
class STK::SymEigen< SquareArray >

The class SymEigen compute the eigenvalue Decomposition of a symmetric ArrayXX.

The decomposition of a symmetric matrix require

• Input: A symmetric matrix A of size (n,n)
• Output:
  1. P Array of size (n,n).
  2. D Vector of dimension n
  3. The matrix A can be copied or overwritten by the class.

The 2-norm (operator norm) of the matrix is given. if the 2-norm is less than the arithmetic precision of the type Real, the rank is not full. Thus the user can be faced with a deficient rank matrix and with a norm and a determinant very small (but not exactly 0.0).

See also

STK::ISymEigen, STK::lapack::SymEigen

Definition at line 85 of file STK_SymEigen.h.

Member Typedef Documentation

Base

template<class SquareArray >

typedef ISymEigen<SymEigen<SquareArray> > STK::SymEigen< SquareArray >::Base

Definition at line 88 of file STK_SymEigen.h.

Constructor & Destructor Documentation

SymEigen() [1/5]

template<class SquareArray >

STK::SymEigen< SquareArray >::SymEigen

(

)

Default Constructor.

Definition at line 144 of file STK_SymEigen.h.

: Base() {}

SymEigen() [2/5]

template<class SquareArray >

STK::SymEigen< SquareArray >::SymEigen	(	SquareArray const &	data,
		bool	ref = false
	)

Constructor.

Parameters

data	reference on a symmetric square matrix
ref	true if we overwrite the data set, false otherwise

Note

data can be a reference if and only if it is a CSquareX

Definition at line 150 of file STK_SymEigen.h.

                               : Base(data)
{}

SymEigen() [3/5]

template<class SquareArray >

template<class Derived >

STK::SymEigen< SquareArray >::SymEigen

(

ExprBase< Derived > const &

data

)

constructor.

Parameters

data	A reference on a symmetric expression matrix to decompose.

Definition at line 167 of file STK_SymEigen.h.

: Base(data) {}

SymEigen() [4/5]

template<class SquareArray >

STK::SymEigen< SquareArray >::SymEigen

(

SymEigen< SquareArray > const &

S

)

Copy constructor.

Parameters

S	the EigenValue to copy

Definition at line 172 of file STK_SymEigen.h.

: Base(S) {}

~SymEigen()

template<class SquareArray >

virtual STK::SymEigen< SquareArray >::~SymEigen ( )

inlinevirtual

virtual destructor

Definition at line 114 of file STK_SymEigen.h.

{}

SymEigen() [5/5]

STK::SymEigen< CSquareX >::SymEigen ( CSquareX const &  data, bool  ref  )

inline

Definition at line 158 of file STK_SymEigen.h.

                                   : Base(data, ref)
{}

Member Function Documentation

compHouse()

template<class SquareArray >

void STK::SymEigen< SquareArray >::compHouse ( )

private

compute the Householder matrix and P

Definition at line 292 of file STK_SymEigen.h.

{
  // compute eigenVectors_
  // iter0 is the column of the Householder vector
  // iter is the current column to compute
  for ( int iter0=range_.lastIdx()-1, iter=range_.lastIdx(), iter1=range_.lastIdx()+1
      ; iter0>=range_.begin()
      ; iter0--, iter--, iter1--)
  {
    // reference on the Householder vector
    typename hidden::Traits<CSquareX>::Col v(eigenVectors_.col(Range(iter, range_.lastIdx(), 0), iter0));
    // Get Beta
    Real beta = v[iter];
    if (beta)
    {
      // Compute Col iter -> eigenVectors_ e_{iter}= e_{iter}+ beta v
      eigenVectors_(iter, iter) = 1.0 + beta /* *1.0 */;
      for (int i=iter1; i<range_.end(); i++)
      { eigenVectors_(i, iter) = beta * v[i];}
 
      // Update the other Cols
      for (int i=iter1; i<range_.end(); i++)
      {
        // compute beta*<v, eigenVectors_^i>
        Real aux = 0.0;
        for (int j=iter1; j<range_.end(); j++)
        { aux += eigenVectors_(j, i) * v[j]; }
        aux *= beta;
        // range_.begin() row (iter)
        eigenVectors_(iter, i) = aux;
        // other rows
        for (int j=iter1; j<range_.end(); j++)
        { eigenVectors_(j, i) += aux * v[j];}
      }
    }
    else // beta = 0, nothing to do
    { eigenVectors_(iter,iter)=1.0;
      for (int j=iter1; j<range_.end(); j++)
      { eigenVectors_(iter, j) =0.0; eigenVectors_(j, iter) = 0.0;}
    }
  }
  // range_.begin() row and range_.begin() col
  eigenVectors_(range_.begin(), range_.begin()) = 1.0;
  for (int j=range_.begin()+1; j<range_.end(); j++)
  { eigenVectors_(range_.begin(), j) = 0.0; eigenVectors_(j, range_.begin()) = 0.0;}
}

diagonalize()

template<class SquareArray >

void STK::SymEigen< SquareArray >::diagonalize ( )

private

computing the diagonalization of eigenValues_ and F_

Definition at line 341 of file STK_SymEigen.h.

{
  // Diagonalisation of eigenVectors_
  for (int iend=range_.lastIdx(); iend>=range_.begin()+1; iend--)
  {
    int iter;
    for (iter=0; iter<MAXITER; iter++) // fix the number of iterations max
    {
      // check cv for the last element
      Real sum = std::abs(eigenValues_[iend]) + std::abs(eigenValues_[iend-1]);
      // if the first element of the small subdiagonal
      // is 0. we stop the QR iterations and increment iend
      if (std::abs(F_[iend-1])+sum == sum)
      { F_[iend-1] = 0.0 ; break;}
      // look for a single small subdiagonal element to split the matrix
      int ibeg = iend-1;
      while (ibeg>range_.begin())
      {
        ibeg--;
        // if a subdiagonal is zero, we get a sub matrix unreduced
        //sum = std::abs(eigenValues_[ibeg])+std::abs(eigenValues_[ibeg+1]);
        if (std::abs(F_[ibeg])+std::abs(eigenValues_[ibeg]) == std::abs(eigenValues_[ibeg]))
        { F_[ibeg] = 0.; ibeg++; break;}
      };
      // QR rotations between the rows/cols ibeg et iend
      // Computation of the Wilkinson's shift
      Real aux = (eigenValues_[iend-1] - eigenValues_[iend])/(2.0 * F_[iend-1]);
      // initialisation of the matrix of rotation
      // y is the current F_[k-1],
      Real y = eigenValues_[ibeg]-eigenValues_[iend] + F_[iend-1]/sign(aux, STK::norm(aux,1.0));
      // z is the element to annulate
      Real z       = F_[ibeg];
      // Fk is the temporary F_[k]
      Real Fk      = z;
      // temporary DeltaD(k)
      Real DeltaDk = 0.0;
      // Index of the columns
      int k,k1;
      // Givens rotation to restore tridiaonal form
      for (k=ibeg, k1=ibeg+1; k<iend; k++, k1++)
      {
        // underflow ? we have a tridiagonal form exit now
        if (z == 0.0) { F_[k]=Fk;  break;}
        // Rotation columns (k,k+1)
        F_[k-1] = (aux = STK::norm(y,z));    // compute final F_[k-1]
        // compute cosinus and sinus
        Real cosinus = y/aux, sinus = -z/aux;
        Real Dk   = eigenValues_[k] + DeltaDk;      // compute current D[k]
        Real aux1 = 2.0 * cosinus * Fk + sinus * (Dk - eigenValues_[k1]);
        // compute DeltaD(k+1)
        eigenValues_[k]     = Dk - (DeltaDk =  sinus * aux1);  // compute eigenValues_[k]
        y         = cosinus*aux1 - Fk;    // temporary F_[k]
        Fk        = cosinus * F_[k1];     // temporary F_[k+1]
        z         = -sinus * F_[k1];      // temporary z
        // update eigenVectors_
        rightGivens(eigenVectors_, k1, k, cosinus, sinus);
      }
      // k=iend if z!=0 and k<iend if z==0
      eigenValues_[k] += DeltaDk ; F_[k-1] = y;
      // restore F[ibeg-1]
      F_[ibeg-1] = 0.;
    } // iter
    if (iter == MAXITER)
    { this->msg_error_ = _T("Warning, max iteration reached in SymEigen::diagonalize()\n");}
    // We have to sort the eigenvalues : we use a basic strategy
    Real z = eigenValues_[iend];        // current value
    for (int i=iend+1; i<eigenValues_.end(); i++)
    { if (eigenValues_[i] > z)           // if the ith eigenvalue is greater
      { eigenValues_.swap(i-1, i);       // swap the cols
        eigenVectors_.swapCols(i-1, i);
      }
      else break;
    } // end sort
  } // iend
  // sort first value
  Real z = eigenValues_[range_.begin()];        // current value
  for (int i=range_.begin()+1; i<eigenValues_.end(); i++)
  { if (eigenValues_[i] > z)           // if the ith eigenvalue is greater
    {
      eigenValues_.swap(i-1, i);       // swap the cols
      eigenVectors_.swapCols(i-1, i);
    }
    else break;
  } // end sort
  // clear auxiliary array
  F_.clear();
}

References _T, MAXITER, STK::IRunnerBase::msg_error_, STK::norm(), STK::rightGivens(), STK::sign(), and STK::sum().

operator=()

template<class SquareArray >

SymEigen & STK::SymEigen< SquareArray >::operator= ( const SymEigen< SquareArray > &  S )

inline

Operator = : overwrite the SymEigen with S.

Parameters

S	SymEigen to copy

Returns

a reference on this

Definition at line 124 of file STK_SymEigen.h.

    {
      Base::operator=(S);
      return *this;
    }

References STK::ISymEigen< SymEigen< SquareArray > >::operator=().

runImpl()

template<class SquareArray >

bool STK::SymEigen< SquareArray >::runImpl

(

)

Diagonalization of eigenVectors_.

Returns

true if no error occur, false otherwise

Definition at line 177 of file STK_SymEigen.h.

{
#ifdef STK_ALGEBRA_VERBOSE
    stk_cout << _T("Entering in SymEigen::run()\n");
#endif
  try
  {
    // copy data
    eigenValues_.resize(eigenVectors_.range());
    norm_ = 0.;
    rank_ = 0;
    det_ = 0.;
#ifdef STK_ALGEBRA_VERBOSE
    stk_cout << _T("calling SymEigen::tridiagonalize()\n");
#endif
    // tridiagonalize eigenVectors_
    tridiagonalize();
#ifdef STK_ALGEBRA_VERBOSE
    stk_cout << _T("calling SymEigen::compHouse()\n");
#endif
    // compute eigenVectors_
    compHouse();
#ifdef STK_ALGEBRA_VERBOSE
    stk_cout << _T("calling SymEigen::diagonalize()\n");
#endif
    // Diagonalize
    diagonalize();
#ifdef STK_ALGEBRA_VERBOSE
    stk_cout << _T("calling SymEigen::finalize()\n");
#endif
    // compute rank, norm and determinant
    this->finalizeStep();
  }
  catch (Exception const& e)
  {
    this->msg_error_ = e.error();
    return false;
  }
  return true;
}

References _T, STK::Exception::error(), STK::IRegression< Array, Array, Weight >::finalizeStep(), STK::IRunnerBase::msg_error_, and stk_cout.

tridiagonalize()

template<class SquareArray >

void STK::SymEigen< SquareArray >::tridiagonalize ( )

private

compute the tri-diagonalization of eigenVectors_

Definition at line 223 of file STK_SymEigen.h.

{
  // Upper diagonal values
  F_.resize(Range(range_.begin()-1, range_.lastIdx(), 0));
  F_.front() = 0.0; F_.back() =0.0;
  // initial range of the Householder vectors
  Range range1(Range(range_.begin()+1, range_.lastIdx(), 0));
  Range range2(Range(range_.begin()+2, range_.lastIdx(), 0));
  // Bidiagonalization of eigenVectors_
  // loop on the cols and rows
  for ( int i=range_.begin(), i1=range_.begin()+1, i2=range_.begin()+2; i<range_.lastIdx()
      ; i++, i1++, i2++, range1.incFirst(1), range2.incFirst(1)
      )
  {
    // ref on the current column iter in the range iter1:last
    typename hidden::Traits<CSquareX>::Col v(eigenVectors_.col(range1, i));
    // Compute Householder vector and get subdiagonal element
    F_[i] = house(v);
    // Save diagonal element
    eigenValues_[i] = eigenVectors_(i,i);
    // get beta
    Real beta = v.front();
    if (beta)
    {
      // ref on the current column iter1 in the range iter1:last
      typename hidden::Traits<CSquareX>::Col M1(eigenVectors_.col(range1, i1));
      // aux1 will contain <v,p>
      Real aux1 = 0.0;
      // apply left and right Householder to eigenVectors_
      // compute D(iter1:range_.last_) = p and aux1 = <p,v>
      // Computation of p_iter1 = beta * eigenVectors_ v using the lower part of eigenVectors_
      // save p_iter1 in the unused part of eigenValues_ and compute p'v
      aux1 += (eigenValues_[i1] = beta*(M1[i1] + M1.sub(range2).dot(v.sub(range2)))) /* *1.0 */;
      // other cols
      for (int k=i2; k<range_.end(); k++)
      {
        // Computation of p_i = beta * eigenVectors_ v using the lower part of eigenVectors_
        // save p_i in the unusued part of eigenValues_ and compute p'v
        Real aux = M1[k] /* *1.0 */;
        for (int j=i2;  j<=k;     j++) { aux += eigenVectors_(k,j)*v[j];}
        for (int j=k+1; j<range_.end(); j++) { aux += eigenVectors_(j,k)*v[j];}
        // save p_i in the unusued part of eigenValues_ and compute p'v
        aux1 += (eigenValues_[k] = beta*aux) * v[k];
      }
      // update diagonal element M_ii+= 2 v_i * q_i = 2* q_i (i=iter1)
      // aux = q_iter1 and aux1 = beta*<p,v>/2 (we don't save aux in eigenValues_)
      Real aux = (eigenValues_[i1] + (aux1 *= (beta/2.0)));
      M1[i1] += 2.0*aux;
      // update lower part: all rows
      // compute q_i and update the lower part of eigenVectors_
      for (int k=i2; k<range_.end(); k++)
      {
        // get q_i and save it in eigenValues_i=q_i = p_i + <p,v> * beta * v_i/2
        eigenValues_[k] += aux1 * v[k];
        // Compute eigenVectors_ + u q' + q u',
        // update the row i, range_.begin() element
        M1[k] += eigenValues_[k] /* *1.0 */+ v[k] * aux;
        // update the row i: all cols under the diagonal
        for (int j=i2; j<=k; j++)
          eigenVectors_(k,j) += v[j] * eigenValues_[k] + v[k] * eigenValues_[j];
      }
    }
  }
  // Last col
  eigenValues_[range_.lastIdx()] = eigenVectors_(range_.lastIdx(),range_.lastIdx());
}

References STK::house().

Member Data Documentation

F_

template<class SquareArray >

VectorX STK::SymEigen< SquareArray >::F_

private

Temporary vector.

Values of the subdiagonal.

Definition at line 132 of file STK_SymEigen.h.

---

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

• STK_SymEigen.h