STK++ 1.0
Public Member Functions | Static Public Member Functions | Protected Attributes | Private Member Functions | Private Attributes

STK::Eigenvalues Class Reference
The Algebra project

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

#include <STK_Eigenvalues.h>

List of all members.

Public Member Functions

 Eigenvalues (const MatrixSquare &A=MatrixSquare(), const bool &ref=false)
 Eigenvalues (const Eigenvalues &S)
virtual ~Eigenvalues ()
Eigenvalues & operator= (const Eigenvalues &S)
Range range () const
 Range of P_.
Real matrixNorm () const
 norm of the matrix
Integer matrixRank () const
 rank of the matrix
bool error () const
 if an error occur during diagsym()
MatrixSquare const & rotation () const
 get rotation matrix
Vector const & eigenvalues () const
MatrixSquare * ginv ()
void ginv (MatrixSquare &res)

Static Public Member Functions

static Real tridiag (MatrixSquare &P, Point &D, Point &F)
 compute the tri-diagonalization of the matrix P
static void compHouse (MatrixSquare &P)
 compute the Householder matrix P
static bool diag (MatrixSquare &P, Point &D, Point &F, const bool &withP=true)
 compute the tri-diagonalization of the tri-diagonal matrix

Protected Attributes

MatrixSquare P_
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
bool error_
 Everything OK during computation ?

Private Member Functions

void diagsym ()
 Diagonalization algorithm of P_.
void tridiag ()
 compute the tri-diagonalization of P_
void compHouse ()
 compute the Householder matrix and P
void diag ()
 computing the diagonalization of D_ and F_

Private Attributes

Vector F_
 Values of the sub-diagonal.

Detailed Description

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

The decomposition of a symmetric matrix require

  • Input: Symmetric matrix A(n,n)
  • Output:
    1. P Matrix (n,n).
    2. D Point (n)
    3. A = PDP' The matrix A can be copied or overwritten by the class.

Definition at line 56 of file STK_Eigenvalues.h.

Constructor & Destructor Documentation

home iovleff Developpement workspace stkpp projects Algebra src STK_Eigenvalues cpp home iovleff Developpement workspace stkpp projects Algebra src STK_Eigenvalues cpp STK::Eigenvalues::Eigenvalues ( const MatrixSquare &  A = MatrixSquare(),
const bool &  ref = false 
)

Default constructor.

Parameters:
AThe symmetric matrix to decompose.
reftrue if A is overwritten.

Definition at line 51 of file STK_Eigenvalues.cpp.

{
  if (A.empty()) return;
  // diagonalize
  diagsym();
}

/* Copy constructor                                                                                                                                */
STK::Eigenvalues::Eigenvalues ( const Eigenvalues &  S)

Copy constructor.

Parameters:
Sthe EigenValue to copy

Definition at line 66 of file STK_Eigenvalues.cpp.

{ ;}

/* Operator = : overwrite the Eigenvalues with S.
STK::Eigenvalues::~Eigenvalues ( ) [virtual]

virtual destructor

Definition at line 92 of file STK_Eigenvalues.cpp.

Member Function Documentation

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

Operator = : overwrite the Eigenvalues with S.

Parameters:
SEigenvalues to copy
Returns:
a reference on this

Definition at line 78 of file STK_Eigenvalues.cpp.

Range STK::Eigenvalues::range ( ) const [inline]

Range of P_.

Definition at line 105 of file STK_Eigenvalues.h.

References first_, and last_.

{ return Range(first_, last_);}
Real STK::Eigenvalues::matrixNorm ( ) const [inline]

norm of the matrix

Definition at line 107 of file STK_Eigenvalues.h.

References norm_.

{ return norm_;}
Integer STK::Eigenvalues::matrixRank ( ) const [inline]

rank of the matrix

Definition at line 109 of file STK_Eigenvalues.h.

References rank_.

{ return rank_;}
bool STK::Eigenvalues::error ( ) const [inline]

if an error occur during diagsym()

Definition at line 111 of file STK_Eigenvalues.h.

References error_.

Referenced by diag().

{ return error_;}
MatrixSquare const& STK::Eigenvalues::rotation ( ) const [inline]

get rotation matrix

Definition at line 113 of file STK_Eigenvalues.h.

References P_.

Referenced by STK::LocalVariance::computeAxis().

{ return P_;}
Vector const& STK::Eigenvalues::eigenvalues ( ) const [inline]

get Diagonal

Definition at line 115 of file STK_Eigenvalues.h.

References D_.

Referenced by STK::LocalVariance::computeAxis().

{ return D_;}
MatrixSquare * STK::Eigenvalues::ginv ( )

Compute the generalized inverse of the diagonalized 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_Eigenvalues.cpp.

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

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

  {
    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;
}

/*
void STK::Eigenvalues::ginv ( MatrixSquare &  res)

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

Parameters:
resthe generalized inverse of the matrix P_.

Definition at line 129 of file STK_Eigenvalues.cpp.

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

  {
    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;
    }
  }
}
Real STK::Eigenvalues::tridiag ( MatrixSquare &  P,
Point &  D,
Point &  F 
) [static]

compute the tri-diagonalization of the matrix P

Parameters:
Pthe matrix to tri-diagonalize, overwritten by the Householder vectors.
Dthe Vector of the diagonal
Fthe Vector of the sub-diagonal
Returns:
the estimated norm of the Matrix

Definition at line 437 of file STK_Eigenvalues.cpp.

  {
    // 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 P1(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 * M v using the lower part of M
      Real aux = P1[iter1] /* *1.0 */;
      for (Integer j=iter2; j<=last; j++) { aux += P1[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 * M v using the lower part of M
        aux = P1[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)));
      P1[iter1] += 2.0*aux;
      // update lower part: all rows
      // compute q_i and update the lower part of M
      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 M + u q' + q u',
        // update the row i, first element
        P1[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];
      }
    }
    // Estimation of the norm
    norm = max(abs(D[iter])+abs(F[iter]), norm);
  }
  // Last col: F[last] = 0.0;
  D[last] = P(last,last);
  // return norm estimated
  return sqrt(norm);
}

// diagonalize D_ and F_
void STK::Eigenvalues::compHouse ( MatrixSquare &  P) [static]

compute the Householder matrix P

Parameters:
Pthe matrix with the Householder vectors

Definition at line 382 of file STK_Eigenvalues.cpp.

References beta(), and STK::Funct::P.

  {
    // 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;}
}
bool STK::Eigenvalues::diag ( MatrixSquare &  P,
Point &  D,
Point &  F,
const bool &  withP = true 
) [static]

compute the tri-diagonalization of the tri-diagonal matrix

Parameters:
Pthe Householder Matrix
Dthe Vector of the diagonal
Fthe Vector of the sub-diagonal
withPtrue if P is needed
Returns:
true if the iteration terminate normaly, false otherwise.

Definition at line 519 of file STK_Eigenvalues.cpp.

References STK::abs(), error(), STK::norm(), STK::rightGivens(), STK::sign(), and STK::sum().

{
  // error detected in the iterations
  bool error =false;
  // dimensions
  Integer first = D.first(), last =D.last(); 
  // Diagonalisation of P_
  for (Integer iend=last; iend>=first+1; iend--)
  {
    Integer iter;
    for (iter=0; iter<30; iter++) // 30 iterations max
    {
      // check Cv fo 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, 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 = 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_
        if (withP) 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 == 30) { error = true;}
  } // iend
  return error;
}

} // Namespace STK
void STK::Eigenvalues::diagsym ( ) [private]

Diagonalization algorithm of P_.

Definition at line 154 of file STK_Eigenvalues.cpp.

void STK::Eigenvalues::tridiag ( ) [private]

compute the tri-diagonalization of P_

Definition at line 168 of file STK_Eigenvalues.cpp.

  {
    // 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_);
}

// Compute P_ from the Householder rotations
void STK::Eigenvalues::compHouse ( ) [private]

compute the Householder matrix and P

Definition at line 239 of file STK_Eigenvalues.cpp.

References beta(), last_, and P_.

  {
    // 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;}
}

// diagonalize D_ and F_
void STK::Eigenvalues::diag ( ) [private]

computing the diagonalization of D_ and F_

Definition at line 287 of file STK_Eigenvalues.cpp.

References STK::abs(), D_, error_, F_, first_, STK::IContainer1D::last(), STK::norm(), P_, STK::rightGivens(), STK::sign(), STK::sum(), STK::IArray1DBase< TYPE, PTRELT, TArray1D >::swap(), and STK::IArray2DBase< TYPE, PTRCOL, TArrayHo, TArrayVe, TArray2D >::swapCols().

  {
    Integer iter;
    for (iter=0; iter<30; 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, 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 = 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 == 30) { 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

  // compute norm_
  norm_ = D_[first_];
  // compute tolerance
  Real tol = Arithmetic<Real>::epsilon() * norm_;
  // compute rank_
  for (Integer i = first_; i<= last_; i++)
  {
    if (abs(D_[i])> tol ) rank_++;
  }
}

// Compute P_ from the Householder rotations

Member Data Documentation

Vector STK::Eigenvalues::F_ [private]

Values of the sub-diagonal.

Definition at line 60 of file STK_Eigenvalues.h.

Referenced by diag().

MatrixSquare STK::Eigenvalues::P_ [protected]

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

Definition at line 68 of file STK_Eigenvalues.h.

Referenced by compHouse(), diag(), ginv(), and rotation().

Vector STK::Eigenvalues::D_ [protected]

Array of the eigenvalues.

Definition at line 70 of file STK_Eigenvalues.h.

Referenced by diag(), eigenvalues(), and ginv().

Integer STK::Eigenvalues::first_ [protected]

first row/col of P_

Definition at line 72 of file STK_Eigenvalues.h.

Referenced by diag(), ginv(), and range().

Integer STK::Eigenvalues::last_ [protected]

last row/col of P_

Definition at line 74 of file STK_Eigenvalues.h.

Referenced by compHouse(), ginv(), and range().

Real STK::Eigenvalues::norm_ [protected]

norm of the matrix

Definition at line 76 of file STK_Eigenvalues.h.

Referenced by matrixNorm().

Integer STK::Eigenvalues::rank_ [protected]

rank of the matrix

Definition at line 78 of file STK_Eigenvalues.h.

Referenced by matrixRank().

bool STK::Eigenvalues::error_ [protected]

Everything OK during computation ?

Definition at line 80 of file STK_Eigenvalues.h.

Referenced by diag(), and error().

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