STK_ISvd.h

Go to the documentation of this file.

/*--------------------------------------------------------------------*/
/*     Copyright (C) 2004-2016  Serge Iovleff, Université Lille 1, Inria
 
    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 : S..._Dot_I..._At_stkpp_Dot_org (see copyright for ...)
*/
 
/*
 * Project:  stkpp::Algebra
 * Purpose:  Define The Interface ISvd Class.
 * Author:   Serge Iovleff, S..._Dot_I..._At_stkpp_Dot_org (see copyright for ...)
 *
 **/
 
#ifndef STK_ISVD_H
#define STK_ISVD_H
 
#include <STKernel.h>
 
#include "STK_Algebra_Util.h"
 
namespace STK
{
template<class Derived>
class ISvd: public IRunnerBase, public IRecursiveTemplate<Derived>
{
  protected:
    typedef typename hidden::AlgebraTraits<Derived>::ArrayU ArrayU;
    typedef typename hidden::AlgebraTraits<Derived>::ArrayD ArrayD;
    typedef typename hidden::AlgebraTraits<Derived>::ArrayV ArrayV;
    typedef typename ArrayU::Type Type;
 
    inline ISvd( ArrayU const& A, bool ref, bool withU = true, bool withV = true)
               : U_(A, ref), V_(), D_()
               , withU_(withU), withV_(withV), norm_(0), rank_(0), trace_(0), det_(0)
    {}
    template<class OtherDerived>
    inline ISvd( ArrayBase<OtherDerived> const& A, bool withU = true, bool withV = true)
               : U_(A), V_(), D_()
               , withU_(withU), withV_(withV), norm_(0), rank_(0), trace_(0), det_(0)
    {}
    inline ISvd( ISvd const& S)
               : U_(S.U_, S.U_.isRef()), V_(S.V_), D_(S.D_)
               , withU_(S.withU_), withV_(S.withV_)
               , norm_(S.norm_), rank_(S.rank_), trace_(S.trace_), det_(S.det_)
    {}
    inline virtual ~ISvd() {}
    ISvd& operator=(const ISvd &S)
    {
      U_ = S.U_; V_ = S.V_;  D_ = S.D_;
      withU_ =  S.withU_; withV_ = S.withV_;
      norm_ = S.norm_; rank_ = S.rank_; trace_ = S.trace_; det_ = S.det_;
      return *this;
    }
    virtual void finalize()
    {
      // compute sign and trace
      trace_ = 0.;
      int s = 1;
      for (int i=D_.begin(); i< D_.end(); ++i )
      {
        Type value = D_[i];
        s *= sign(value);
        trace_ += value;
      }
 
      // compute norm_ (2-norm) and determinant
      norm_ = D_.maxElt();
      det_  = 0;
      if (D_.abs().minElt() > 0)
      { det_ = s * std::exp(D_.abs().log().sum());}
 
      // compute norm_ rank_
      rank_ = D_.size(); // full rank
      Type tol = norm_ * Arithmetic<Type>::epsilon();
      for (int i=D_.begin(); i< D_.end(); ++i )
      { if (std::abs(D_[i]) < tol) { rank_--;}}
 
    }
 
  public:
    // getters
    inline Type det() const { return det_;}
    inline Type trace()  const { return trace_;}
    inline Type norm()  const { return norm_;}
    inline int rank()  const { return rank_;}
    inline ArrayU const& U() const { return U_;}
    inline ArrayV const& V() const { return V_;}
    inline ArrayD const&  D() const { return D_;}
 
    virtual bool run()
    {
      if (U_.empty()) { return true;}
      // compute Svd decomposition
      if (!this->asDerived().runImpl()) return false;
      finalize();
      return true;
    }
    template<class OtherArray>
    void setData( OtherArray const& A, bool withU = true, bool withV = true)
    {
      U_ = A;           // Copy A in U_
      withU_ = withU;   // copy withU_ value
      withV_ = withV;   // copy withV_ value
      V_.resize(0,0), D_.resize(0);
    }
    template<class OtherArray>
    OtherArray& ginv(OtherArray& res) const;
 
  protected:
    ArrayU U_;
    ArrayV V_;
    ArrayD D_;
    bool withU_;
    bool withV_;
    Type norm_;
    int rank_;
    Type trace_;
    Type det_;
 
    inline int nrowU() const { return U_.sizeRows();}
    inline int ncolU() const { return U_.sizeCols();}
    inline int nrowD() const { return D_.sizeRows();}
    inline int ncolD() const { return D_.sizeCols();}
    inline int nrowV() const { return V_.sizeRows();}
    inline int ncolV() const { return V_.sizeCols();}
};
 
/* Compute the generalized inverse of the matrix and put the result in res.
 *  @param res array with the result
 *  @return the result
 */
template<class Derived>
template<class OtherArray>
OtherArray& ISvd<Derived>::ginv(OtherArray& res) const
{
  Type tol = Arithmetic<Type>::epsilon() * norm_;
  if(tol==0) { tol = Arithmetic<Type>::epsilon();}
  // compute D
  if (U_.cols() != D_.range())
  {
    Array2DDiagonal<Real> diag(U_.cols(), 0.0);
    for (int i= D_.begin(); i< D_.end(); ++i) { diag[i]=D_[i];}
    res = V_ * diag.safeInverse(tol) * U_.transpose();
  }
  else
  {
    res = V_ * D_.diagonalize().safeInverse(tol) * U_.transpose();
  }
  // compute and return UD^{-1}V'
  return (res);
}
 
} // namespace STK
 
#endif /* STK_ISVD_H */