STK_Eigenvalues.cpp

Go to the documentation of this file.

/*--------------------------------------------------------------------*/
/*     Copyright (C) 2004-2007  Serge Iovleff

    This program is free software; you can redistribute it and/or modify
    it under the terms of the GNU Lesser General Public License as
    published by the Free Software Foundation; either version 2 of the
    License, or (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU Lesser General Public License for more details.

    You should have received a copy of the GNU Lesser General Public
    License along with this program; if not, write to the
    Free Software Foundation, Inc.,
    59 Temple Place,
    Suite 330,
    Boston, MA 02111-1307
    USA

    Contact : Serge.Iovleff@stkpp.org
*/

/*
 * Project:  Algebra
 * Purpose:  Implementation of the EigenValue class
 * Author:   Serge Iovleff, serge.iovleff@stkpp.org
 *
 **/

#include "../include/STK_Vector.h"
#include "../include/STK_Householder.h"
#include "../include/STK_Givens.h"
#include "../include/STK_LinAlgebra1D.h"
#include "../include/STK_TExpAlgebra.h"

#include "../include/STK_Eigenvalues.h"

namespace STK
{
/* Constructors.
  *  A is the matrix to decompose.
  **/
Eigenvalues::Eigenvalues( MatrixSquare const& A, const bool &ref)
                        : P_(A, ref)
                        , D_(A.rangeVe())
                        , first_(A.firstCol())
                        , last_(A.lastCol())
                        , norm_(0)
                        , rank_(0)
                        , error_(false)
{
  if (A.empty()) return;
  // diagonalize
  diagsym();
}

/* Copy constructor                                                                                                                                */
Eigenvalues::Eigenvalues( const Eigenvalues &S)
                        : P_(S.P_)
                        , D_(S.D_)
                        , first_(S.first_)
                        , last_(S.last_)
                        , norm_(S.norm_)
                        , rank_(S.rank_)
                        , error_(S.error_)
{ ;}

/* Operator = : overwrite the Eigenvalues with S.
  **/
Eigenvalues& Eigenvalues::operator=(const Eigenvalues &S)
{
  P_     = S.P_;        // Matrix P
  D_     = S.D_;        // Singular values
  first_ = S.first_;    // first value
  last_  = S.last_;     // last value
  norm_  = S.norm_;     // rank of the matrix
  rank_  = S.rank_;     // rank of the matrix
  error_ = S.error_;    // Everything OK during computation ?

  return *this;
}

/* destructor  */
Eigenvalues::~Eigenvalues()
{ ;}

/*
 * Compute the generalized inverse of the diagonalized matrix.
 * The result is allocated dynamically and is not liberated by this
 * class.
 * @return a pointer on the generalized inverse.
 */
MatrixSquare* Eigenvalues::ginv()
{
  // 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;
}

/*
 * Compute the generalized inverse of the diagonalized matrix P_ and put
 * the result in res.
 * @param res the generalized inverse of the matrix P_.
 */
void Eigenvalues::ginv(MatrixSquare& res)
{
  // 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;
    }
  }
}


/* Main methods.  */
/* Compute diagonalization of the symmetric matrix                     */
void Eigenvalues::diagsym()
{
  // tridiagonalize P_
  tridiag();
  // compute P_
  compHouse();  
  // Diagonalize
  diag(); 
}

/* tridiagonalisation of the symetric matrix P_. Only the lower
 * part of P_ used. P_ is overwritten with the Householder vectors.
 * D_ contains the diagonal. F_ contains the subdiagonal.
 **/
void Eigenvalues::tridiag()
{
  // 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_);
}

// Compute P_ from the Householder rotations
void Eigenvalues::compHouse()
{
  // 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;}
}

// diagonalize D_ and F_
void Eigenvalues::diag()
{
  // Diagonalisation of P_
  for (Integer iend=last_; iend>=first_+1; iend--)
  {
    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
void Eigenvalues::compHouse(MatrixSquare &P)
{
  // compute the number of iteration
  const Integer first = P.firstRow();
  const Integer last  = P.lastRow();
  // 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;}
}


/* tridiagonalisation of the symetric matrix P. Only the lower
 * part of P used. P is overwritten with the Householder vectors.
 * D contains the diagonal. F contains the subdiagonal.
 **/
Real Eigenvalues::tridiag( MatrixSquare &P, Point &D, Point& F)
{
  // norm of the matrix P
  Real norm  = 0.0;
  // compute the number of iteration
  Integer first = P.firstRow();
  Integer last  = P.lastRow();
  // Diagonal values
  D.resize(P.rangeVe());
  // 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 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_
bool Eigenvalues::diag( MatrixSquare &P, Point &D, Point &F
                      , const bool &withP)
{
  // 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