The class EigenvaluesSymmetric compute the eigenvalue Decomposition of a symmetric Matrix. More...

#include <STK_EigenvaluesSymmetric.h>

Inheritance diagram for STK::EigenvaluesSymmetric:

Collaboration diagram for STK::EigenvaluesSymmetric:

Public Member Functions
	EigenvaluesSymmetric (MatrixSquare const *A)
	Default constructor.
	EigenvaluesSymmetric (const EigenvaluesSymmetric &S)
	Copy constructor.
virtual	~EigenvaluesSymmetric ()
	virtual destructor
virtual bool	run ()
	Diagonalization of P_.
virtual bool	run (Vector const &weights)
	weighted diagonalization of P_
EigenvaluesSymmetric &	operator= (const EigenvaluesSymmetric &S)
	Operator = : overwrite the EigenvaluesSymmetric with S.
Real	norm () const
	norm of the matrix
Integer	rank () const
	rank of the matrix
Real	det () const
	determinant of the Matrix
MatrixSquare const &	rotation () const
	get rotation matrix
Vector const &	eigenvalues () const
	get eigenvalues
Integer const &	first () const
	get the first index of the rows/columns
Integer const &	last () const
	get last index of the rows/columns
MatrixSquare *	ginv ()
	Compute the generalized inverse of the symmetric matrix.
void	ginv (MatrixSquare &res)
	Compute the generalized inverse of the symmetric matrix and put the result in res.
Protected Attributes
MatrixSquare	P_
	P_ Square matrix or the eigenvectors.
Vector	D_
	Array of the eigenvalues.
Integer	first_
	first row/col of P_
Integer	last_
	last row/col of P_
Real	norm_
	norm of the matrix
Integer	rank_
	rank of the matrix
Real	det_
	determinant of the Matrix
bool	error_
	Everything OK during computation ?
Private Types
typedef IRunnerPtr2D< Real, Matrix >	IRunnerPtrMatrix
	concrete type of the runner
Private Member Functions
void	tridiagonalize ()
	compute the tri-diagonalization of P_
void	compHouse ()
	compute the Householder matrix and P
void	diagonalize ()
	computing the diagonalization of D_ and F_
void	compEstimates ()
	compute rank, norm and determinant.
Private Attributes
Vector	F_
	Temporary vector.

Detailed Description

The class EigenvaluesSymmetric compute the eigenvalue Decomposition of a symmetric Matrix.

The decomposition of a symmetric matrix require

Input: A symmetric matrix A of size (n,n)
Output:
1. P Matrix 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 set to 0. Thus the user can be faced with a rank 0 matrix and with a norm and a determinant very small.

Definition at line 62 of file STK_EigenvaluesSymmetric.h.

Member Typedef Documentation

typedef IRunnerPtr2D<Real, Matrix> STK::EigenvaluesSymmetric::IRunnerPtrMatrix [private]

concrete type of the runner

Definition at line 65 of file STK_EigenvaluesSymmetric.h.

Constructor & Destructor Documentation

STK::EigenvaluesSymmetric::EigenvaluesSymmetric ( MatrixSquare const * A )

Default constructor.

Parameters:

A	The symmetric matrix to decompose.

Definition at line 52 of file STK_EigenvaluesSymmetric.cpp.

                                          : IRunnerPtrMatrix(p_data)
                                          , P_(p_data->range())
                                          , D_(p_data->range())
                                          , first_(p_data->first())
                                          , last_(p_data->last())
                                          , norm_(0)
                                          , rank_(0)
                                          , det_(0)
{ }

STK::EigenvaluesSymmetric::EigenvaluesSymmetric ( const EigenvaluesSymmetric & S )

Copy constructor.

Parameters:

S	the EigenValue to copy

Definition at line 64 of file STK_EigenvaluesSymmetric.cpp.

                                          : IRunnerPtrMatrix(S)
                                          , P_(S.P_)
                                          , D_(S.D_)
                                          , first_(S.first_)
                                          , last_(S.last_)
                                          , norm_(S.norm_)
                                          , rank_(S.rank_)
                                          , det_(S.det_)
{ ;}

STK::EigenvaluesSymmetric::~EigenvaluesSymmetric ( ) [virtual]

virtual destructor

Definition at line 92 of file STK_EigenvaluesSymmetric.cpp.

{ ;}

Member Function Documentation

bool STK::EigenvaluesSymmetric::run ( ) [virtual]

Diagonalization of P_.

Returns:: true if no error occur, false otherwise

Implements STK::IRunnerBase.

Definition at line 154 of file STK_EigenvaluesSymmetric.cpp.

References STK::IRecursiveTemplate< Leaf >::asLeaf(), compEstimates(), compHouse(), D_, det_, diagonalize(), STK::Exception::error(), error_, first_, STK::IContainer2D::firstCol(), last_, STK::IContainer2D::lastCol(), norm_, P_, STK::IRunnerPtr2D< Real, Matrix >::p_data_, STK::IContainer2D::rangeHo(), rank_, STK::IContainer1D::resize(), and tridiagonalize().

Referenced by STK::LocalVariance::computeAxis(), STK::GinvSymmetric::operator()(), and STK::Law::MultivariateNormal< Container1D >::update().

{
  try
  {
    // copy data
    P_ = p_data_->asLeaf();
    D_.resize(p_data_->rangeHo());
    first_ = p_data_->firstCol();
    last_ = p_data_->lastCol();
    norm_ = 0.;
    rank_ = 0;
    det_ = 0.;
    // tridiagonalize P_
    tridiagonalize();
    // compute P_
    compHouse();
    // Diagonalize
    diagonalize();
    // compute rank, norm and determinant
    compEstimates();
  }
  catch (Exception e)
  {
    error_ = e.error();
    return false;
  }
  return true;
}

Here is the call graph for this function:

bool STK::EigenvaluesSymmetric::run ( Vector const & weights ) [virtual]

weighted diagonalization of P_

Parameters:

weights the weights of each rows

Returns:: true if no error occur, false otherwise

Definition at line 185 of file STK_EigenvaluesSymmetric.cpp.

References STK::IRecursiveTemplate< Leaf >::asLeaf(), compEstimates(), compHouse(), D_, det_, diagonalize(), STK::Exception::error(), error_, STK::MatrixSquare::first(), first_, STK::MatrixSquare::last(), last_, norm_, P_, STK::IRunnerPtr2D< Real, Matrix >::p_data_, STK::MatrixSquare::range(), STK::ITContainer1D< TYPE, TContainer1D >::range(), STK::IContainer2D::rangeVe(), rank_, STK::IContainer1D::resize(), and tridiagonalize().

{
  // copy data
  if (p_data_->rangeVe() != weights.range())
  { throw runtime_error("In EigenvaluesSymmetric::run(weights) "
                             "p_data_->rangeVe() != weights.range().\n");
  }

  try
  {
    // copy data
    // copy data
    P_ = p_data_->asLeaf();
    D_.resize(P_.range());
    first_ = P_.first();
    last_ = P_.last();
    norm_ = 0.;
    rank_ = 0;
    det_ = 0.;
    // weight rows and columns of the container
    for (Integer i= first_; i <= last_; ++i)
    {
      Real aux = Real(sqrt(weights[i]));
      P_(i) *= aux;
      P_[i] *= aux;
    }

    // tridiagonalize P_
    tridiagonalize();
    // compute P_
    compHouse();
    // Diagonalize
    diagonalize();
    // compute rank, norm and determinant
    compEstimates();

    // Unweighed rows and columns of the rotation
    for (Integer i= first_; i <= last_; ++i)
    {
      Real aux = Real(sqrt(weights[i]));
      if (aux)
      {
        P_(i) /= aux;
        P_[i] /= aux;
      }
    }
  }
  catch (Exception e)
  {
    error_ = e.error();
    return false;
  }
  return true;
}

Here is the call graph for this function:

EigenvaluesSymmetric & STK::EigenvaluesSymmetric::operator= ( const EigenvaluesSymmetric & S )

Operator = : overwrite the EigenvaluesSymmetric with S.

Parameters:

S	EigenvaluesSymmetric to copy

Returns:: a reference on this

Definition at line 77 of file STK_EigenvaluesSymmetric.cpp.

References D_, det_, error_, first_, last_, norm_, P_, and rank_.

{
  P_     = S.P_;        // Matrix P
  D_     = S.D_;        // Singular values
  first_ = S.first_;    // first value
  last_  = S.last_;     // last value
  norm_  = S.norm_;     // norm of the matrix
  rank_  = S.rank_;     // rank of the matrix
  det_   = S.det_;      // determinant of the matrix
  error_ = S.error_;    // Everything OK during computation ?

  return *this;
}

Real STK::EigenvaluesSymmetric::norm ( ) const [inline]

norm of the matrix

Returns:: the norm of the matrix

Definition at line 101 of file STK_EigenvaluesSymmetric.h.

References norm_.

Referenced by diagonalize().

{ return norm_;}

Integer STK::EigenvaluesSymmetric::rank ( ) const [inline]

rank of the matrix

Returns:: the rank of the matrix

Definition at line 105 of file STK_EigenvaluesSymmetric.h.

References rank_.

Referenced by STK::Law::MultivariateNormal< Container1D >::update().

{ return rank_;}

Real STK::EigenvaluesSymmetric::det ( ) const [inline]

determinant of the Matrix

Returns:: the determinant of the Matrix

Definition at line 109 of file STK_EigenvaluesSymmetric.h.

References det_.

Referenced by STK::Law::MultivariateNormal< Container1D >::logLikelihood(), STK::Law::MultivariateNormal< Container1D >::lpdf(), and STK::Law::MultivariateNormal< Container1D >::pdf().

{ return det_;}

MatrixSquare const& STK::EigenvaluesSymmetric::rotation ( ) const [inline]

get rotation matrix

Returns:: the rotation matrix

Definition at line 113 of file STK_EigenvaluesSymmetric.h.

References P_.

Referenced by STK::LocalVariance::computeAxis(), STK::Law::MultivariateNormal< Container1D >::logLikelihood(), STK::Law::MultivariateNormal< Container1D >::lpdf(), and STK::Law::MultivariateNormal< Container1D >::update().

{ return P_;}

Vector const& STK::EigenvaluesSymmetric::eigenvalues ( ) const [inline]

get eigenvalues

Returns:: the eigenvalues

Definition at line 117 of file STK_EigenvaluesSymmetric.h.

References D_.

Referenced by STK::LocalVariance::computeAxis(), and STK::Law::MultivariateNormal< Container1D >::update().

{ return D_;}

Integer const& STK::EigenvaluesSymmetric::first ( ) const [inline]

get the first index of the rows/columns

Returns:: the index of the first row/column

Definition at line 121 of file STK_EigenvaluesSymmetric.h.

References first_.

{ return first_;}

Integer const& STK::EigenvaluesSymmetric::last ( ) const [inline]

get last index of the rows/columns

Returns:: the index of the last row/column

Definition at line 125 of file STK_EigenvaluesSymmetric.h.

References last_.

{ return last_;}

MatrixSquare * STK::EigenvaluesSymmetric::ginv ( )

Compute the generalized inverse of the symmetric matrix.

The result is allocated dynamically and is not liberated by this class.

Returns:: A pointer on the generalized inverse.

Definition at line 101 of file STK_EigenvaluesSymmetric.cpp.

References STK::abs(), D_, first_, last_, norm_, P_, STK::MatrixSquare::range(), and STK::sum().

Referenced by STK::LocalVariance::computeAxis(), and STK::GinvSymmetric::operator()().

{
  // create pseudo inverse matrix
  MatrixSquare* res = new MatrixSquare(P_.range());
  // compute tolerance
  Real tol = Arithmetic<Real>::epsilon() * norm_;
  // compute PDP'
  for (Integer i = first_; i<=last_; i++)
  {
    for (Integer j = first_; j<=last_; j++)
    {
      Real sum = 0.0;
      for (Integer k = first_; k<=last_; k++)
      {
        if (abs(D_[k]) > tol)
          sum += (P_(i, k) * P_(j, k)) / D_[k];
      }
      (*res)(j,i) = sum;
    }
  }
  return res;
}

Here is the call graph for this function:

void STK::EigenvaluesSymmetric::ginv ( MatrixSquare & res )

Compute the generalized inverse of the symmetric matrix and put the result in res.

Parameters:

res	the generalized inverse of the Matrix.

Definition at line 129 of file STK_EigenvaluesSymmetric.cpp.

References STK::abs(), D_, first_, last_, norm_, P_, STK::MatrixSquare::range(), STK::MatrixSquare::resize(), and STK::sum().

{
  // create pseudo inverse matrix
  res.resize(P_.range());
  // compute tolerance
  Real tol = Arithmetic<Real>::epsilon() * norm_;
  // compute PDP'
  for (Integer i = first_; i<=last_; i++)
  {
    for (Integer j = first_; j<=last_; j++)
    {
      Real sum = 0.0;
      for (Integer k = first_; k<=last_; k++)
      {
        if (abs(D_[k]) > tol)
          sum += (P_(i, k) * P_(j, k)) / D_[k];
      }
      res(j,i) = sum;
    }
  }
}

Here is the call graph for this function:

void STK::EigenvaluesSymmetric::tridiagonalize ( ) [private]

compute the tri-diagonalization of P_

Definition at line 244 of file STK_EigenvaluesSymmetric.cpp.

References STK::ITContainer1D< TYPE, TContainer1D >::back(), beta(), D_, F_, first_, STK::ITContainer1D< TYPE, TContainer1D >::front(), STK::house(), STK::Range::incFirst(), last_, P_, and STK::IContainer1D::resize().

Referenced by run().

{
  // Upper diagonal values
  F_.resize(Range(first_-1, last_));
  F_.front() = 0.0; F_.back() =0.0;
  // initial range of the Householder vectors
  Range range1(Range(first_+1, last_));
  // Bidiagonalisation of P_
  // loop on the cols and rows
  for ( Integer iter=first_, iter1=first_+1, iter2=first_+2
      ; iter<last_
      ; iter++, iter1++, iter2++, range1.incFirst()
      )
  {
    // ref on the current column iter in the range iter1:last_
    Vector v(P_, range1, iter);
    // Compute Householder Vector and get subdiagonal element
    F_[iter] = house(v);
    // Save diagonal element
    D_[iter] = P_(iter,iter);
    // get beta
    Real beta = v.front();
    if (beta)
    {
      // ref on the current column iter1 in the range iter1:last_
      Vector M1(P_, range1, iter1);
      // aux1 will contain <v,p>
      Real aux1 = 0.0;
      // apply left and right Householder to P_
      // compute D(iter1:last_) = p and aux1 = <p,v>
      // Computation of p_iter1 = beta * P_ v using the lower part of P_
      Real aux = M1[iter1] /* *1.0 */;
      for (Integer j=iter2; j<=last_; j++) { aux += M1[j]*v[j];}
      // save p_iter1 in the unusued part of D_ and compute p'v
      aux1 += (D_[iter1] = beta*aux) /* *1.0 */;
      // other cols
      for (Integer i=iter2; i<=last_; i++)
      {
        // Computation of p_i = beta * P_ v using the lower part of P_
        aux = M1[i] /* *1.0 */;
        for (Integer j=iter2; j<=i;    j++)  { aux += P_(i,j)*v[j];}
        for (Integer j=i+1;   j<=last_; j++) { aux += P_(j,i)*v[j];}
        // save p_i in the unusued part of D_ and compute p'v
        aux1 += (D_[i] = beta*aux) * v[i];
      }
      // 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 D_)
      aux = (D_[iter1] + (aux1 *= (beta/2.0)));
      M1[iter1] += 2.0*aux;
      // update lower part: all rows
      // compute q_i and update the lower part of P_
      for (Integer i=iter2; i<=last_; i++)
      {
        // get v_i
        Real v_i = v[i];
        // get q_i and save it in D_i=q_i = p_i + <p,v> * beta * v_i/2
        Real q_i = (D_[i] += aux1 * v_i);
        // Compute P_ + u q' + q u',
        // update the row i, first_ element
        M1[i] += q_i /* *1.0 */+ v_i * aux;
        // update the row i: all cols under the diagonal
        for (Integer j=iter2; j<=i; j++)
          P_(i,j) += v[j] * q_i + v_i * D_[j];
      }
    }
  }
  // Last col: F[last_] = 0.0;
  D_[last_] = P_(last_,last_);
}

Here is the call graph for this function:

void STK::EigenvaluesSymmetric::compHouse ( ) [private]

compute the Householder matrix and P

Definition at line 315 of file STK_EigenvaluesSymmetric.cpp.

References beta(), first_, last_, and P_.

Referenced by run().

{
  // compute P_
  // iter0 is the column of the Householder vector
  // iter is the current column to compute
  for ( Integer iter0=last_-1, iter=last_, iter1=last_+1
      ; iter0>=first_
      ; iter0--, iter--, iter1--)
  {
    // reference on the Householder vector
    Vector v(P_, Range(iter, last_), iter0);
    // Get Beta
    Real beta = v[iter];
    if (beta)
    {
      // Compute Col iter -> P_ e_{iter}= e_{iter}+ beta v
      P_(iter, iter) = 1.0 + beta /* *1.0 */;
      for (Integer i=iter1; i<=last_; i++)
      { P_(i, iter) = beta * v[i];}

      // Update the other Cols
      for (Integer i=iter1; i<=last_; i++)
      {
        // compute beta*<v, P_^i>
        Real aux = 0.0;
        for (Integer j=iter1; j<=last_; j++)
        { aux += P_(j, i) * v[j]; }
        aux *= beta;
        // first_ row (iter)
        P_(iter, i) = aux;
        // other rows
        for (Integer j=iter1; j<=last_; j++)
        { P_(j, i) += aux * v[j];}
      }
    }
    else // beta = 0, nothing to do
    { P_(iter,iter)=1.0;
      for (Integer j=iter1; j<=last_; j++)
      { P_(iter, j) =0.0; P_(j, iter) = 0.0;}
    }
  }
  // first_ row and first_ col
  P_(first_, first_) = 1.0;
  for (Integer j=first_+1; j<=last_; j++)
  { P_(first_, j) = 0.0; P_(j, first_) = 0.0;}
}

Here is the call graph for this function:

void STK::EigenvaluesSymmetric::diagonalize ( ) [private]

computing the diagonalization of D_ and F_

Definition at line 363 of file STK_EigenvaluesSymmetric.cpp.

References STK::abs(), D_, error_, F_, first_, STK::ITContainer1D< TYPE, TContainer1D >::last(), last_, MAXITER, norm(), STK::norm(), P_, STK::rightGivens(), STK::sign(), STK::sum(), STK::IArray1DBase< TYPE, PTRELT, TArray1D >::swap(), and STK::ITArray2DBase< TYPE, PTRCOL, TArray2D >::swapCols().

Referenced by run().

{
  // Diagonalisation of P_
  for (Integer iend=last_; iend>=first_+1; iend--)
  {
    Integer iter;
    for (iter=0; iter<MAXITER; iter++) // 30 iterations max
    {
      // check cv for the last element
      Real sum = abs(D_[iend]) + abs(D_[iend-1]);
      // if the first element of the small subdiagonal
      // is 0. we stop the QR iterations and increment iend
      if (abs(F_[iend-1])+sum == sum)
      { F_[iend-1] = 0.0 ; break;}
      // look for a single small subdiagonal element to split the matrix
      Integer ibeg = iend-1;
      while (ibeg>first_)
      {
        ibeg--;
        // if a subdiagonal is zero, we get a sub matrix unreduced
        sum = abs(D_[ibeg])+abs(D_[ibeg+1]);
        if (abs(F_[ibeg])+sum == sum) { F_[ibeg] = 0.; ibeg++; break;}
      };
      // QR rotations between the rows/cols ibeg et iend
      // Computation of the Wilkinson's shift
      Real aux = (D_[iend-1] - D_[iend])/(2.0 * F_[iend-1]);
      // initialisation of the matrix of rotation
      // y is the current F_[k-1],
      Real y = D_[ibeg]-D_[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
      Integer 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   = D_[k] + DeltaDk;      // compute current D[k]
        Real aux1 = 2.0 * cosinus * Fk + sinus * (Dk - D_[k1]);
        // compute DeltaD(k+1)
        D_[k]     = Dk - (DeltaDk =  sinus * aux1);  // compute D_[k]
        y         = cosinus*aux1 - Fk;    // temporary F_[k]
        Fk        = cosinus * F_[k1];     // temporary F_[k+1]
        z         = -sinus * F_[k1];      // temporary z
        // update P_
        rightGivens(P_, k1, k, cosinus, sinus);
      }
      // k=iend if z!=0 and k<iend if z==0
      D_[k] += DeltaDk ; F_[k-1] = y;
      // restore F[ibeg-1]
      F_[ibeg-1] = 0.;
    } // iter
    if (iter == MAXITER) { error_ = true;}
    // We have to sort the eigenvalues : we use a basic strategy
    Real z = D_[iend];        // current value
    for (Integer i=iend+1; i<=D_.last(); i++)
    { if (D_[i] > z)           // if the ith eigenvalue is greater
      { D_.swap(i-1, i);       // swap the cols
        P_.swapCols(i-1, i);
      }
      else break;
    } // end sort
  } // iend
  // sort first value
  Real z = D_[first_];        // current value
  for (Integer i=first_+1; i<=D_.last(); i++)
  { if (D_[i] > z)                // if the ith eigenvalue is greater
    {
      D_.swap(i-1, i);       // swap the cols
      P_.swapCols(i-1, i);
    }
    else break;
  } // end sort
}

Here is the call graph for this function:

void STK::EigenvaluesSymmetric::compEstimates ( ) [private]

compute rank, norm and determinant.

Definition at line 448 of file STK_EigenvaluesSymmetric.cpp.

References STK::abs(), D_, det_, first_, last_, norm_, and rank_.

Referenced by run().

{
  // compute 2-norm_
  norm_ = D_[first_];
  // trivial case
  if (norm_ < Arithmetic<Real>::epsilon())
  {
    rank_ = 0;
    det_ = 1;
    // compute determinant
    for (Integer i = first_; i<= last_; i++)
    { det_ *= D_[i];}
    return;
  }
  // compute tolerance
  Real tol = Arithmetic<Real>::epsilon() * norm_;
  det_ = 1;
  // compute rank_and determinant
  for (Integer i = first_; i<= last_; i++)
  {
    det_ *= D_[i];
    if (abs(D_[i])> tol ) rank_++;
  }
}

Here is the call graph for this function:

Member Data Documentation

MatrixSquare STK::EigenvaluesSymmetric::P_ [protected]

P_ Square matrix or the eigenvectors.

P_ is initialized using the matrix passe to the constructor. The initialization can be done by reference, in this case A will be overwritten in the diagonalization process.

Definition at line 146 of file STK_EigenvaluesSymmetric.h.

Referenced by compHouse(), diagonalize(), ginv(), operator=(), rotation(), run(), and tridiagonalize().