TDecompSVD

class TDecompSVD : public TDecompBase

    protected:

                    static Bool_t Bidiagonalize(TMatrixD& v, TMatrixD& u, TVectorD& sDiag, TVectorD& oDiag)
                      static void Diag_1(TMatrixD& v, TVectorD& sDiag, TVectorD& oDiag, Int_t k)
                      static void Diag_2(TVectorD& sDiag, TVectorD& oDiag, Int_t k, Int_t l)
                      static void Diag_3(TMatrixD& v, TMatrixD& u, TVectorD& sDiag, TVectorD& oDiag, Int_t k, Int_t l)
                    static Bool_t Diagonalize(TMatrixD& v, TMatrixD& u, TVectorD& sDiag, TVectorD& oDiag)
      virtual const TMatrixDBase& GetDecompMatrix() const
                      static void SortSingular(TMatrixD& v, TMatrixD& u, TVectorD& sDiag)

    public:

                       TDecompSVD()
                       TDecompSVD(Int_t nrows, Int_t ncols)
                       TDecompSVD(Int_t row_lwb, Int_t row_upb, Int_t col_lwb, Int_t col_upb)
                       TDecompSVD(const TMatrixD& m, Double_t tol = 0.0)
                       TDecompSVD(const TDecompSVD& another)
               virtual ~TDecompSVD()
        static TClass* Class()
      virtual Double_t Condition()
        virtual Bool_t Decompose()
          virtual void Det(Double_t& d1, Double_t& d2)
        const TMatrixD GetMatrix() const
         virtual Int_t GetNcols() const
         virtual Int_t GetNrows() const
       const TVectorD& GetSig() const
       const TMatrixD& GetU() const
       const TMatrixD& GetV() const
                  void Invert(TMatrixD& inv)
              TMatrixD Invert()
       virtual TClass* IsA() const
           TDecompSVD& operator=(const TDecompSVD& source)
          virtual void Print(Option_t* opt = "") const
          virtual void SetMatrix(const TMatrixD& a)
          virtual void ShowMembers(TMemberInspector& insp, char* parent)
        virtual Bool_t Solve(TVectorD& b)
      virtual TVectorD Solve(const TVectorD& b, Bool_t& ok)
        virtual Bool_t Solve(TMatrixDColumn& b)
          virtual void Streamer(TBuffer& b)
                  void StreamerNVirtual(TBuffer& b)
        virtual Bool_t TransSolve(TVectorD& b)
      virtual TVectorD TransSolve(const TVectorD& b, Bool_t& ok)
        virtual Bool_t TransSolve(TMatrixDColumn& b)

Data Members

    protected:

      TMatrixD fU    orthogonal matrix
      TMatrixD fV    orthogonal matrix
      TVectorD fSig  diagonal of diagonal matrix

    public:

      static const enum TDecompSVD:: kWorkMax  

Class Description

                                                                       
 Single Value Decomposition class                                      
                                                                       
 For an (m x n) matrix A with m >= n, the singular value decomposition 
 is                                                                    
 an (m x m) orthogonal matrix fU, an (m x n) diagonal matrix fS, and   
 an (n x n) orthogonal matrix fV so that A = U*S*V'.                   
                                                                       
 If the row/column index of A starts at (rowLwb,colLwb) then the       
 decomposed matrices/vectors start at :                                
  fU   : (rowLwb,colLwb)                                               
  fV   : (colLwb,colLwb)                                               
  fSig : (colLwb)                                                      
                                                                       
 The diagonal matrix fS is stored in the singular values vector fSig . 
 The singular values, fSig[k] = S[k][k], are ordered so that           
 fSig[0] >= fSig[1] >= ... >= fSig[n-1].                               
                                                                       
 The singular value decompostion always exists, so the decomposition   
 will (as long as m >=n) never fail. If m < n, the user should add     
 sufficient zero rows to A , so that m == n                            
                                                                       
 Here fTol is used to set the threshold on the minimum allowed value   
 of the singular values:                                               
  min_singular = fTol*max(fSig[i])                                     
                                                                       

 Bidiagonalize the (m x n) - matrix a (stored in v) through a series of Householder
 transformations applied to the left (Q^T) and to the right (H) of a ,
 so that A = Q . C . H^T with matrix C bidiagonal. Q and H are orthogonal matrices .

 Output:
   v     - (n x n) - matrix H in the (n x n) part of v
   u     - (m x m) - matrix Q^T
   sDiag - diagonal of the (m x n) C
   oDiag - off-diagonal elements of matrix C

  Test code for the output:
    const Int_t nRow = v.GetNrows();
    const Int_t nCol = v.GetNcols();
    TMatrixD H(v); H.ResizeTo(nCol,nCol);
    TMatrixD E1(nCol,nCol); E1.UnitMatrix();
    TMatrixD Ht(TMatrixDBase::kTransposed,H);
    Bool_t ok = kTRUE;
    ok &= VerifyMatrixIdentity(Ht * H,E1,kTRUE,1.0e-13);
    ok &= VerifyMatrixIdentity(H * Ht,E1,kTRUE,1.0e-13);
    TMatrixD E2(nRow,nRow); E2.UnitMatrix();
    TMatrixD Qt(u);
    TMatrixD Q(TMatrixDBase::kTransposed,Qt);
    ok &= VerifyMatrixIdentity(Q * Qt,E2,kTRUE,1.0e-13);
    TMatrixD C(nRow,nCol);
    TMatrixDDiag(C) = sDiag;
    for (Int_t i = 0; i < nCol-1; i++)
      C(i,i+1) = oDiag(i+1);
    TMatrixD A = Q*C*Ht;
    ok &= VerifyMatrixIdentity(A,a,kTRUE,1.0e-13);

 Diagonalizes in an iterativ fashion the bidiagonal matrix C as described through
 sDiag and oDiag, so that S' = U'^T . C . V' is diagonal. U' and V' are orthogonal
 matrices .

 Output:
   v     - (n x n) - matrix H . V' in the (n x n) part of v
   u     - (m x m) - matrix U'^T . Q^T
   sDiag - diagonal of the (m x n) S'

   return convergence flag:  0 -> no convergence
                             1 -> convergence

  Test code for the output:
    const Int_t nRow = v.GetNrows();
    const Int_t nCol = v.GetNcols();
    TMatrixD tmp = v; tmp.ResizeTo(nCol,nCol);
    TMatrixD Vprime  = Ht*tmp;
    TMatrixD Vprimet(TMatrixDBase::kTransposed,Vprime);
    TMatrixD Uprimet = u*Q;
    TMatrixD Uprime(TMatrixDBase::kTransposed,Uprimet);
    TMatrixD Sprime(nRow,nCol);
    TMatrixDDiag(Sprime) = sDiag;
    ok &= VerifyMatrixIdentity(Uprimet * C * Vprime,Sprime,kTRUE,1.0e-13);
    ok &= VerifyMatrixIdentity(Q*Uprime * Sprime * Vprimet * Ht,a,kTRUE,1.0e-13);

 Perform a permutation transformation on the diagonal matrix S', so that
 matrix S'' = U''^T . S' . V''  has diagonal elements ordered such that they
 do not increase.

 Output:
   v     - (n x n) - matrix H . V' . V'' in the (n x n) part of v
   u     - (m x m) - matrix U''^T . U'^T . Q^T
   sDiag - diagonal of the (m x n) S''

 Reconstruct the original matrix using the decomposition parts

 Solve Ax=b assuming the SVD form of A is stored . Solution returned in b.
 If A is of size (m x n), input vector b should be of size (m), however,
 the solution, returned in b, will be in the first (n) elements .

 For m > n , x  is the least-squares solution of min(A . x - b)

 Solve Ax=b assuming the SVD form of A is stored . Solution returned in the
 matrix column cb b.
 If A is of size (m x n), input vector b should be of size (m), however,
 the solution, returned in b, will be in the first (n) elements .

 For m > n , x  is the least-squares solution of min(A . x - b)

 Solve A^T x=b assuming the SVD form of A is stored . Solution returned in b.

 For a matrix A(m,n), its inverse A_inv is defined as A * A_inv = A_inv * A = unit
 The user should always supply a matrix of size (m x m) !
 If m > n , only the (n x m) part of the returned (pseudo inverse) matrix
 should be used .

 For a matrix A(m,n), its inverse A_inv is defined as A * A_inv = A_inv * A = unit
 (n x m) Ainv is returned .

Inline Functions

                       void ~TDecompSVD()
        const TMatrixDBase& GetDecompMatrix() const
                 TDecompSVD TDecompSVD(const TDecompSVD& another)
                      Int_t GetNrows() const
                      Int_t GetNcols() const
            const TMatrixD& GetU() const
            const TMatrixD& GetV() const
            const TVectorD& GetSig() const
                     Bool_t Solve(TMatrixDColumn& b)
                     Bool_t TransSolve(TMatrixDColumn& b)
                TDecompSVD& operator=(const TDecompSVD& source)
                    TClass* Class()
                    TClass* IsA() const
                       void ShowMembers(TMemberInspector& insp, char* parent)
                       void Streamer(TBuffer& b)
                       void StreamerNVirtual(TBuffer& b)