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. *


ROOT page - Class index - Class Hierarchy - Top of the page

This page has been automatically generated. If you have any comments or suggestions about the page layout send a mail to ROOT support, or contact the developers with any questions or problems regarding ROOT.