STK::Svd< Array > Class Template Reference

The class Svd compute the Singular Value Decomposition of a Array with the Golub-Reinsch Algorithm. More...

#include <STK_Svd.h>

Inheritance diagram for STK::Svd< Array >:

Public Types
typedef ISvd< Svd< Array > >	Base
typedef hidden::Traits< Array >::Col	ColVector
typedef hidden::Traits< Array >::Row	RowVector

Public Member Functions
	Svd (Array const &A, bool ref=false, bool withU=true, bool withV=true)
	Default constructor.
template<class OtherArray >
	Svd (ArrayBase< OtherArray > const &A, bool withU=true, bool withV=true)
	constructor with other kind of array/expression
	Svd (const Svd &S)
	Copy Constructor.
virtual	~Svd ()
	destructor.
Svd &	operator= (const Svd &S)
	Operator = : overwrite the Svd with S.
bool	runImpl ()
	run the Svd
Public Member Functions inherited from STK::ISvd< Svd< Array > >
Type	det () const
Type	trace () const
Type	norm () const
int	rank () const
ArrayU const &	U () const
ArrayV const &	V () const
ArrayD const &	D () const
virtual bool	run ()
	implement the run method
void	setData (OtherArray const &A, bool withU=true, bool withV=true)
	Set a new data set to ISvd class.
OtherArray &	ginv (OtherArray &res) const
	Compute the generalized inverse of the matrix and put the result in res.
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.

Static Public Member Functions
static bool	diag (ArrayDiagonalX &D, VectorX &F, Array &U, ArraySquareX &V, bool withU=true, bool withV=true, Real const &tol=Arithmetic< Real >::epsilon())
	Computing the diagonalization of a bi-diagonal matrix.
static void	rightEliminate (ArrayDiagonalX &D, VectorX &F, int const &nrow, ArraySquareX &V, bool withV=true, Real const &tol=Arithmetic< Real >::epsilon())
	right eliminate the element on the subdiagonal of the row nrow

Private Member Functions
bool	computeSvd ()
	Svd main steps.
void	compU ()
	Compute U (if withU_ is true)
void	compV ()
	Compute V (if withV_ is true)

Private Attributes
VectorX	F_
	Values of the Sub-diagonal.

Additional Inherited Members
Protected Types inherited from STK::ISvd< Svd< Array > >
typedef hidden::AlgebraTraits< Svd< Array > >::ArrayU	ArrayU
typedef hidden::AlgebraTraits< Svd< Array > >::ArrayD	ArrayD
typedef hidden::AlgebraTraits< Svd< Array > >::ArrayV	ArrayV
typedef ArrayU::Type	Type
Protected Member Functions inherited from STK::ISvd< Svd< Array > >
	ISvd (ArrayU const &A, bool ref, bool withU=true, bool withV=true)
	Default constructor.
	ISvd (ArrayBase< OtherDerived > const &A, bool withU=true, bool withV=true)
	constructor with other kind of array/expression
	ISvd (ISvd const &S)
	Copy Constructor.
virtual	~ISvd ()
	destructor.
ISvd &	operator= (const ISvd &S)
	Operator = : overwrite the ISvd with S.
virtual void	finalize ()
	Finalize any operations that have to be done after the computation of the decomposition.
int	nrowU () const
int	ncolU () const
int	nrowD () const
int	ncolD () const
int	nrowV () const
int	ncolV () const
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::ISvd< Svd< Array > >
ArrayU	U_
	U_ matrix.
ArrayV	V_
	V_ matrix.
ArrayD	D_
	Diagonal array of the singular values.
bool	withU_
	Compute U_ ?
bool	withV_
	Compute V_ ?
Type	norm_
	trace norm
int	rank_
	rank
Type	trace_
	trace norm
Type	det_
	determinant
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 Array>
class STK::Svd< Array >

The class Svd compute the Singular Value Decomposition of a Array with the Golub-Reinsch Algorithm.

The method take as:

• input: a matrix A(nrow,ncol)
• output:
  1. U Array (nrow,ncol).
  2. D diagonal matrix (min(norw,ncol))
  3. V Array (ncol,ncol). and perform the decomposition:
• A = UDV' (transpose V). U can have more columns than A, and it is possible to compute some (all) vectors of Ker(A).

See also

SK::ISvd, STK::lapack::Svd

Definition at line 116 of file STK_Svd.h.

Member Typedef Documentation

Base

template<class Array >

typedef ISvd< Svd<Array> > STK::Svd< Array >::Base

Definition at line 119 of file STK_Svd.h.

ColVector

template<class Array >

typedef hidden::Traits<Array>::Col STK::Svd< Array >::ColVector

Definition at line 120 of file STK_Svd.h.

RowVector

template<class Array >

typedef hidden::Traits<Array>::Row STK::Svd< Array >::RowVector

Definition at line 121 of file STK_Svd.h.

Constructor & Destructor Documentation

Svd() [1/3]

template<class Array >

STK::Svd< Array >::Svd ( Array const &  A, bool  ref = false, bool  withU = true, bool  withV = true  )

inline

Default constructor.

Parameters

A	the matrix to decompose.
ref	if true, U_ is a reference of A.
withU	if true, we save the left housolder transforms in U_.
withV	if true, we save the right housolder transforms in V_.

Definition at line 139 of file STK_Svd.h.

              : Base(A, ref, withU, withV)
    {}

Svd() [2/3]

template<class Array >

template<class OtherArray >

STK::Svd< Array >::Svd ( ArrayBase< OtherArray > const &  A, bool  withU = true, bool  withV = true  )

inline

constructor with other kind of array/expression

Parameters

A	the matrix/expression to decompose.
withU	if true save the left housolder transforms in U_.
withV	if true save the right housolder transforms in V_.

Definition at line 148 of file STK_Svd.h.

              : Base(A, withU, withV) {}

Svd() [3/3]

template<class Array >

STK::Svd< Array >::Svd

(

const Svd< Array > &

S

)

Copy Constructor.

Parameters

S	the Svd to copy

Definition at line 208 of file STK_Svd.h.

: Base(S), F_(S.F_) {}

~Svd()

template<class Array >

virtual STK::Svd< Array >::~Svd ( )

inlinevirtual

destructor.

Definition at line 155 of file STK_Svd.h.

{}

Member Function Documentation

compU()

template<class Array >

void STK::Svd< Array >::compU ( )

private

Compute U (if withU_ is true)

Definition at line 386 of file STK_Svd.h.

{
  int niter = D_.size();            // Number of iterations
  int ncol  = std::min(nrowU(), ncolU()); // number of non zero cols of U_
 
  // initialization of the remaining cols of U_ to 0.0
  // put 0 to unused cols
  U_.col(Range(ncol+1, ncolU(), 0)) = 0.0;
  // Computation of U_
  for (int iter=niter, iter1=niter+1; iter>=1; iter--, iter1--)
  {
    // ref of the column iter
    ColVector X(U_, Range(iter1,nrowU(), 0), iter);
    // ref of the row iter
    RowVector P(U_, Range(iter,ncolU(), 0), iter);
    // Get beta and test
    Real beta = P[iter];
    if (beta)
    {
      // update the column iter
      P[iter] = 1.0 + beta;
      // Updating the cols iter+1 to ncolU_
      for (int j=iter1; j<=niter; j++)
      { // product of U_iter by U_j
        Real aux;
        ColVector Y(U_, Range(iter1, nrowU(), 0), j); // ref on the column j
        // U_j = aux = beta * X'Y
        P[j] = (aux = dot( X, Y) *beta);
        // U^j += aux * U^iter
        for (int i= iter1; i <= nrowU(); ++i)
        {
          Y[i] += X[i] * aux;
        }
      }
      // compute the vector v
      X *= beta;
    }
    else // U^iter = identity
    {
      P[iter] = 1.0;
      X = 0.0;
    }
    // update the column iter
    U_.col(Range(1,iter-1, 0), iter) = 0.0;
  }
}

References STK::dot().

computeSvd()

template<class Array >

bool STK::Svd< Array >::computeSvd ( )

private

Svd main steps.

Definition at line 235 of file STK_Svd.h.

{
  // if the container is empty, there is nothing to do
  if (U_.empty())
  { rank_ = 0;
    norm_ = 0.0;
    return true;
  }
  int beginRow = U_.beginRows(), beginCol = U_.beginCols();
  // if U_ is just a copy of A, translate begin to 1
  // if U_ is a ref on A, this can generate an error
  U_.shift(1,1);
  // Bidiagonalize (U_)
  norm_ = bidiag(U_, D_, F_);
  // right householder vectors are in upper part of U_
  // We need to create V_ before rightEliminate
  if (withV_) { compV();}
  // rightEliminate last element of F_ if any
  if (nrowU() < ncolV())
  { rightEliminate(D_, F_, nrowU(), V_, withV_, norm_);}
  // If (U_) is not needed, we can destroy the storage
  if (withU_) { compU();}
  // Diagonalize
  bool error = diag(D_, F_, U_, V_, withU_, withV_, norm_);
  // The sub diagonal is now zero
  F_.resize(0,0);
  U_.shift(beginRow, beginCol);
  D_.shift(beginCol);
  V_.shift(beginCol);
  return error;
}

References STK::IRunnerBase::error().

compV()

template<class Array >

void STK::Svd< Array >::compV ( )

private

Compute V (if withV_ is true)

Definition at line 324 of file STK_Svd.h.

{
  // Construction of V_
  V_.resize(U_.cols());
  // Number of right Householder rotations
  int  niter = (ncolV()>nrowU()) ? (nrowU()) : (ncolV()-1);
  // initialization of the remaining rows and cols of V_ to Identity
  for (int iter=niter+2; iter<=ncolV(); iter++)
  {
    VectorX W(V_, V_.cols(), iter);
    W       = 0.0;
    W[iter] = 1.0;
  }
 
  Range range1(niter+1, ncolV(), 0), range2(niter+2, ncolV(), 0);
  for ( int iter0=niter, iter1=niter+1, iter2=niter+2; iter0>=1
      ; iter0--, iter1--, iter2--
      , range1.decFirst(1), range2.decFirst(1)
      )
  {
    // get beta and test
    Real beta = U_(iter0, iter1);
    if (beta)
    {
      // ref on the row iter1 of V_
      PointX  Vrow1(V_, range1, iter1);
      // diagonal element
      Vrow1[iter1] = 1.0+beta;
      // ref on the column iter1
      VectorX Vcol1(V_, range2, iter1);
      // get the Householder vector
      Vcol1 = RowVector(U_, range2, iter0);
      // Apply housholder to next cols
      for (int j=iter2; j<=ncolV(); j++)
      {
        Real aux;
        // ref on the column j
        VectorX Vcolj( V_, range2, j);
        // update column j
        Vrow1[j] = (aux = dot(Vcol1, Vcolj) * beta);
        for (int i= iter2; i <= ncolV(); ++i)
        { Vcolj[i]   += Vcol1[i] * aux;}
      }
      // compute the Householder vector
      Vcol1 *= beta;
    }
    else // nothing to do
    {
      V_(range2, iter1) = 0.0;
      V_(iter1, iter1)  = 1.0;
      V_(iter1, range2) = 0.0;
    }
  }
  // First column and rows
  V_(1,1) =1.0;
  V_(Range(2,ncolV(), 0),1) =0.0;
  V_(1,Range(2,ncolV(), 0)) =0.0;
}

References STK::dot().

diag()

template<class Array >

bool STK::Svd< Array >::diag ( ArrayDiagonalX &  D, VectorX &  F, Array &  U, ArraySquareX &  V, bool  withU = true, bool  withV = true, Real const &  tol = Arithmetic<Real>::epsilon()  )

static

Computing the diagonalization of a bi-diagonal matrix.

Parameters

D	the diagonal of the matrix
F	the subdiagonal of the matrix
U	a left orthogonal Array
withU	true if we want to update U
V	a right orthogonal Square Array
withV	true if we want to update V
tol	the tolerance to use

Definition at line 503 of file STK_Svd.h.

{
  // result of the diag process
  bool error = false;
  // Diagonalization of A : Reduction of la matrice bidiagonale
  for (int end=D.lastIdx(); end>=D.begin(); --end)
  { // 30 iter max
    int iter;
    for (iter=1; iter<=MAX_ITER; iter++)
    { // if  the last element of the subdiagonale is 0.0
      // stop the iterations
      int beg;
      if (std::abs(F[end-1])+tol == tol)  { F[end-1] = 0.0; break;}
      // now F[end-1] !=0
      // if  D[end] == 0, we can annulate F[end-1]
      // with rotations of the columns.
      if (std::abs(D[end])+tol == tol)
      {
        D[end]   = 0.0;
        rightEliminate(D, F, end-1, V, withV, tol);
        break; // Last element of the subdiagonal is 0
      }
      // now D[end] != 0 and F[end-1] != 0
      // look for the greatest matrix such that all the elements
      // of the diagonal and subdiagonal are not zeros
      for (beg = end-1; beg>D.begin(); --beg)
      {
        if ((std::abs(D[beg])+tol == tol)||(std::abs(F[beg])+tol == tol))
        break;
      }
      // now F[beg-1]!=0
      // if D[beg] == 0 and F[beg] != 0,
      // we can eliminate the element F[beg]
      // with rotations of the rows
      if ((std::abs(D[beg])+tol == tol) && (std::abs(F[beg])+tol != tol))
      {
        D[beg] = 0.0;
        leftEliminate(D, F, beg, U, withU, tol);
      }
 
      // Si F[beg]==0, on augmente beg
      if (std::abs(F[beg])+tol == tol) { F[beg] = 0.0; beg++;}
 
      // On peut commencer les rotations QR entre les lignes beg et end
      // Shift computation
      // easy shift : commented
      // Real aux = norm(D[end],F[end-1]);
      // Real y   = (D[beg]+aux)*(D[beg]-aux);
      // Real z   = D[beg]*F[beg];
      // Wilkinson shift : look at the doc
      Real dd1 = D[end-1];      // d_1
      Real dd2 = D[end];        // d_2
      Real ff1 = F[end-2];      // f_1
      Real ff2 = F[end-1];      // f_2
      // g
      Real aux = ( (dd1+dd2)*(dd1-dd2)
                 + (ff1-ff2)*(ff1+ff2))/(2*dd1*ff2);
      // A - mu
      Real d1 = D[beg];
      Real f1 = F[beg];
      Real y  = (d1+dd2)*(d1-dd2)
              + ff2*(dd1/(aux+sign(aux,norm(1.0,aux)))- ff2);
      Real z  = d1*f1;
      // chase no null element
      int k, k1;
      for (k=beg, k1 = beg+1; k<end; ++k, ++k1)
      { // Rotation colonnes (k,k+1)
        // Input :  d1 contient D[k],
        //          f1 contient F[k],
        //          y  contient F[k-1]
        // Output : d1 contient F[k],
        //          d2 contient D[k+1],
        //          y  contient D[k]
        Real cosinus=1., sinus=0.;
        Real d2 = D[k1];
        F[k-1] = (aux = norm(y,z));                            // F[k-1]
       // arbitrary rotation if y = z = 0.0
        if (aux)
          y  = (cosinus = y/aux) * d1 - (sinus = -z/aux) * f1; // D[k]
        else
          y  = cosinus * d1 - sinus * f1;                      // D[k]
 
        z      = -sinus * d2;                                  // z
        d1     =  sinus * d1 + cosinus * f1;                   // F[k]
        d2    *=  cosinus;                                     // D[k+1]
        // Update V
        if (withV)
          rightGivens(V, k1, k, cosinus, sinus);
        // avoid underflow
        // Rotation lignes (k,k+1)
        // Input :  d1 contient F[k],
        //          d2 contient D[k+1],
        //          y  contient D[k]
        // Output : d1 contient D[k+1],
        //          f1 contient F[k+1],
        //          y  contient F[k]
        f1   = F[k1];
        D[k] = (aux = norm(y,z));                              // D[k]
        // arbitrary rotation if y = z = 0.0
        if (aux)
          y  = (cosinus = y/aux) * d1 - (sinus = -z/aux) * d2; // F[k]
        else
          y   = cosinus * d1 - sinus * d2;                     // F[k]
 
        z   = -sinus * f1;                                    // z
        d1  = sinus *d1 + cosinus * d2;                       // D[k+1]
        f1 *= cosinus;                                        // F[k+1]
        // Update U
        if (withU)
          rightGivens(U, k1, k, cosinus, sinus);
      } // end of the QR updating iteration
      D[end]   = d1;
      F[end-1] = y;
      F[beg-1] = 0.0;  // F[beg-1] is overwritten, we have to set 0.0
    } // iter
    // too many iterations
    if (iter >= 30) { error = true;}
    // positive singular value only
    if (D[end]< 0.0)
    {
      // change sign of the singular value
      D[end]= -D[end];
      // change sign of the column end of V
      if (withV)
      {
        for (int i= V.beginRows(); i <= V.lastIdxRows(); ++i)
        { V(i,end) = -V(i,end);}
      }
    }
 
    // We have to sort the singular value : we use a basic strategy
    Real z = D[end];        // current value
    for (int i=end+1; i<=D.lastIdx(); i++)
    { if (D[i]> z)                // if the ith singular value is greater
      { D.swap(i-1, i);       // swap the cols
        if (withU) U.swapCols(i-1, i);
        if (withV) V.swapCols(i-1, i);
      }
      else break;
    } // end sort
  } // boucle end
  return error;
}

References STK::IArray2D< Derived >::beginRows(), d1, d2, STK::IRunnerBase::error(), STK::IArray2D< Derived >::lastIdxRows(), MAX_ITER, STK::norm(), STK::rightGivens(), STK::sign(), STK::IArray2D< Derived >::swap(), and STK::IArray2D< Derived >::swapCols().

operator=()

template<class Array >

Svd< Array > & STK::Svd< Array >::operator=

(

const Svd< Array > &

S

)

Operator = : overwrite the Svd with S.

Parameters

S	the Svd to copy

Definition at line 212 of file STK_Svd.h.

{
  U_ = S.U_;
  V_ = S.V_;
  D_ = S.D_;
  withU_ = S.withU_;
  withV_ = S.withV_;
  norm_ = S.norm_;
  rank_ = S.rank_;
  return *this;
}

rightEliminate()

template<class Array >

void STK::Svd< Array >::rightEliminate ( ArrayDiagonalX &  D, VectorX &  F, int const &  nrow, ArraySquareX &  V, bool  withV = true, Real const &  tol = Arithmetic<Real>::epsilon()  )

static

right eliminate the element on the subdiagonal of the row nrow

Parameters

D	the diagonal of the matrix
F	the subdiagonal of the matrix
nrow	the number of the row were we want to rightEliminate
V	a right orthogonal Square Array
withV	true if we want to update V
tol	the tolerance to use

Definition at line 436 of file STK_Svd.h.

{
  // the element to eliminate
  Real z = F[nrow];
  // if the element is not 0.0
  if (std::abs(z)+tol != tol)
  {
    // column of the element to eliminate
    int ncol1 = nrow+1;
    // begin the Givens rotations
    for (int k=nrow, k1=nrow-1; k>=1 ; k--, k1--)
    {
      // compute and apply Givens rotation to the rows (k, k+1)
      Real aux, sinus, cosinus;
      Real y = D[k];
      D[k]   = (aux     = norm(y,z));
      z      = (sinus   = -z/aux) * F[k1];
      F[k1] *= (cosinus =  y/aux);
      // Update V_
      if (withV)
        rightGivens(V, ncol1, k, cosinus, sinus);
      // if 0.0 we can break now
      if (std::abs(z)+tol == tol) break;
    }
  }
  // the element is now 0
  F[nrow] = 0.0;        // is 0.0
}

References STK::norm(), and STK::rightGivens().

runImpl()

template<class Array >

bool STK::Svd< Array >::runImpl

(

)

run the Svd

Definition at line 226 of file STK_Svd.h.

{
  if (!computeSvd()) return false;
  return true;
}

Member Data Documentation

F_

template<class Array >

VectorX STK::Svd< Array >::F_

private

Values of the Sub-diagonal.

Definition at line 196 of file STK_Svd.h.

---

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

• STK_Svd.h