library: libMatrix #include "TDecompSVD.h"

# TDecompSVD

class description - source file - inheritance tree (.pdf)

## class TDecompSVD : public TDecompBase

| Inheritance Chart:
`TObject`
| <- | `TDecompBase`
| <- | `TDecompSVD` |

        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])

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 &a,Double_t tol)

TDecompSVD(const TDecompSVD &another): TDecompBase(another)

Bool_t Decompose()

Bool_t Bidiagonalize(TMatrixD &v,TMatrixD &u,TVectorD &sDiag,TVectorD &oDiag)

     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);

Bool_t Diagonalize(TMatrixD &v,TMatrixD &u,TVectorD &sDiag,TVectorD &oDiag)

     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);

void Diag_1(TMatrixD &v,TVectorD &sDiag,TVectorD &oDiag,Int_t k)

void Diag_2(TVectorD &sDiag,TVectorD &oDiag,Int_t k,Int_t l)

void Diag_3(TMatrixD &v,TMatrixD &u,TVectorD &sDiag,TVectorD &oDiag,Int_t k,Int_t l)

void SortSingular(TMatrixD &v,TMatrixD &u,TVectorD &sDiag)

     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''

const TMatrixD GetMatrix()

     Reconstruct the original matrix using the decomposition parts

void SetMatrix(const TMatrixD &a)

Bool_t Solve(TVectorD &b)

     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)

Bool_t Solve(TMatrixDColumn &cb)

     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)

Bool_t TransSolve(TVectorD &b)

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

Bool_t TransSolve(TMatrixDColumn &cb)

Double_t Condition()

void Det(Double_t &d1,Double_t &d2)

void Invert(TMatrixD &inv)

     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 .

TMatrixD Invert()

     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 .

void Print(Option_t *opt) const

### 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)

Last update: root/matrix:$Name: $:$Id: TDecompSVD.cxx,v 1.22 2005/02/15 16:17:09 brun Exp $
Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. *