ROOT logo
ROOT » MATH » MATRIX » TDecompSVD

class TDecompSVD: public TDecompBase


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


Function Members (Methods)

public:
TDecompSVD()
TDecompSVD(const TDecompSVD& another)
TDecompSVD(Int_t nrows, Int_t ncols)
TDecompSVD(const TMatrixD& m, Double_t tol = 0.0)
TDecompSVD(Int_t row_lwb, Int_t row_upb, Int_t col_lwb, Int_t col_upb)
virtual~TDecompSVD()
voidTObject::AbstractMethod(const char* method) const
virtual voidTObject::AppendPad(Option_t* option = "")
virtual voidTObject::Browse(TBrowser* b)
static TClass*Class()
virtual const char*TObject::ClassName() const
virtual voidTObject::Clear(Option_t* = "")
virtual TObject*TObject::Clone(const char* newname = "") const
virtual Int_tTObject::Compare(const TObject* obj) const
virtual Double_tCondition()
virtual voidTObject::Copy(TObject& object) const
virtual Bool_tDecompose()
virtual voidTObject::Delete(Option_t* option = "")MENU
virtual voidDet(Double_t& d1, Double_t& d2)
virtual Int_tTObject::DistancetoPrimitive(Int_t px, Int_t py)
virtual voidTObject::Draw(Option_t* option = "")
virtual voidTObject::DrawClass() constMENU
virtual TObject*TObject::DrawClone(Option_t* option = "") constMENU
virtual voidTObject::Dump() constMENU
virtual voidTObject::Error(const char* method, const char* msgfmt) const
virtual voidTObject::Execute(const char* method, const char* params, Int_t* error = 0)
virtual voidTObject::Execute(TMethod* method, TObjArray* params, Int_t* error = 0)
virtual voidTObject::ExecuteEvent(Int_t event, Int_t px, Int_t py)
virtual voidTObject::Fatal(const char* method, const char* msgfmt) const
virtual TObject*TObject::FindObject(const char* name) const
virtual TObject*TObject::FindObject(const TObject* obj) const
Int_tTDecompBase::GetColLwb() const
Double_tTDecompBase::GetCondition() const
Double_tTDecompBase::GetDet1() const
Double_tTDecompBase::GetDet2() const
virtual Option_t*TObject::GetDrawOption() const
static Long_tTObject::GetDtorOnly()
virtual const char*TObject::GetIconName() const
const TMatrixDGetMatrix()
virtual const char*TObject::GetName() const
virtual Int_tGetNcols() const
virtual Int_tGetNrows() const
virtual char*TObject::GetObjectInfo(Int_t px, Int_t py) const
static Bool_tTObject::GetObjectStat()
virtual Option_t*TObject::GetOption() const
Int_tTDecompBase::GetRowLwb() const
const TVectorD&GetSig()
virtual const char*TObject::GetTitle() const
Double_tTDecompBase::GetTol() const
const TMatrixD&GetU()
virtual UInt_tTObject::GetUniqueID() const
const TMatrixD&GetV()
virtual Bool_tTObject::HandleTimer(TTimer* timer)
virtual ULong_tTObject::Hash() const
virtual voidTObject::Info(const char* method, const char* msgfmt) const
virtual Bool_tTObject::InheritsFrom(const char* classname) const
virtual Bool_tTObject::InheritsFrom(const TClass* cl) const
virtual voidTObject::Inspect() constMENU
TMatrixDInvert()
Bool_tInvert(TMatrixD& inv)
TMatrixDInvert(Bool_t& status)
voidTObject::InvertBit(UInt_t f)
virtual TClass*IsA() const
virtual Bool_tTObject::IsEqual(const TObject* obj) const
virtual Bool_tTObject::IsFolder() const
Bool_tTObject::IsOnHeap() const
virtual Bool_tTObject::IsSortable() const
Bool_tTObject::IsZombie() const
virtual voidTObject::ls(Option_t* option = "") const
voidTObject::MayNotUse(const char* method) const
virtual Bool_tTDecompBase::MultiSolve(TMatrixD& B)
virtual Bool_tTObject::Notify()
static voidTObject::operator delete(void* ptr)
static voidTObject::operator delete(void* ptr, void* vp)
static voidTObject::operator delete[](void* ptr)
static voidTObject::operator delete[](void* ptr, void* vp)
void*TObject::operator new(size_t sz)
void*TObject::operator new(size_t sz, void* vp)
void*TObject::operator new[](size_t sz)
void*TObject::operator new[](size_t sz, void* vp)
TDecompSVD&operator=(const TDecompSVD& source)
virtual voidTObject::Paint(Option_t* option = "")
virtual voidTObject::Pop()
virtual voidPrint(Option_t* opt = "") constMENU
virtual Int_tTObject::Read(const char* name)
virtual voidTObject::RecursiveRemove(TObject* obj)
voidTObject::ResetBit(UInt_t f)
virtual voidTObject::SaveAs(const char* filename = "", Option_t* option = "") constMENU
virtual voidTObject::SavePrimitive(basic_ostream<char,char_traits<char> >& out, Option_t* option = "")
voidTObject::SetBit(UInt_t f)
voidTObject::SetBit(UInt_t f, Bool_t set)
virtual voidTObject::SetDrawOption(Option_t* option = "")MENU
static voidTObject::SetDtorOnly(void* obj)
virtual voidSetMatrix(const TMatrixD& a)
static voidTObject::SetObjectStat(Bool_t stat)
Double_tTDecompBase::SetTol(Double_t newTol)
virtual voidTObject::SetUniqueID(UInt_t uid)
virtual voidShowMembers(TMemberInspector& insp, char* parent)
virtual Bool_tSolve(TVectorD& b)
virtual Bool_tSolve(TMatrixDColumn& b)
virtual TVectorDSolve(const TVectorD& b, Bool_t& ok)
virtual voidStreamer(TBuffer& b)
voidStreamerNVirtual(TBuffer& b)
virtual voidTObject::SysError(const char* method, const char* msgfmt) const
Bool_tTObject::TestBit(UInt_t f) const
Int_tTObject::TestBits(UInt_t f) const
virtual Bool_tTransSolve(TVectorD& b)
virtual Bool_tTransSolve(TMatrixDColumn& b)
virtual TVectorDTransSolve(const TVectorD& b, Bool_t& ok)
virtual voidTObject::UseCurrentStyle()
virtual voidTObject::Warning(const char* method, const char* msgfmt) const
virtual Int_tTObject::Write(const char* name = 0, Int_t option = 0, Int_t bufsize = 0)
virtual Int_tTObject::Write(const char* name = 0, Int_t option = 0, Int_t bufsize = 0) const
protected:
static Bool_tBidiagonalize(TMatrixD& v, TMatrixD& u, TVectorD& sDiag, TVectorD& oDiag)
static voidDiag_1(TMatrixD& v, TVectorD& sDiag, TVectorD& oDiag, Int_t k)
static voidDiag_2(TVectorD& sDiag, TVectorD& oDiag, Int_t k, Int_t l)
static voidDiag_3(TMatrixD& v, TMatrixD& u, TVectorD& sDiag, TVectorD& oDiag, Int_t k, Int_t l)
static Bool_tDiagonalize(TMatrixD& v, TMatrixD& u, TVectorD& sDiag, TVectorD& oDiag)
static voidTDecompBase::DiagProd(const TVectorD& diag, Double_t tol, Double_t& d1, Double_t& d2)
virtual voidTObject::DoError(int level, const char* location, const char* fmt, va_list va) const
virtual const TMatrixDBase&GetDecompMatrix() const
Int_tTDecompBase::Hager(Double_t& est, Int_t iter = 5)
voidTObject::MakeZombie()
voidTDecompBase::ResetStatus()
static voidSortSingular(TMatrixD& v, TMatrixD& u, TVectorD& sDiag)

Data Members

public:
enum { kWorkMax
};
enum TDecompBase::EMatrixDecompStat { kInit
kPatternSet
kValuesSet
kMatrixSet
kDecomposed
kDetermined
kCondition
kSingular
};
enum TDecompBase::[unnamed] { kWorkMax
};
enum TObject::EStatusBits { kCanDelete
kMustCleanup
kObjInCanvas
kIsReferenced
kHasUUID
kCannotPick
kNoContextMenu
kInvalidObject
};
enum TObject::[unnamed] { kIsOnHeap
kNotDeleted
kZombie
kBitMask
kSingleKey
kOverwrite
kWriteDelete
};
protected:
Int_tTDecompBase::fColLwbColumn lower bound of decomposed matrix
Double_tTDecompBase::fConditionmatrix condition number
Double_tTDecompBase::fDet1determinant mantissa
Double_tTDecompBase::fDet2determinant exponent for powers of 2
Int_tTDecompBase::fRowLwbRow lower bound of decomposed matrix
TVectorDfSigdiagonal of diagonal matrix
Double_tTDecompBase::fTolsqrt(epsilon); epsilon is smallest number number so that 1+epsilon > 1
TMatrixDfUorthogonal matrix
TMatrixDfVorthogonal matrix

Class Charts

Inheritance Inherited Members Includes Libraries
Class Charts

Function documentation

TDecompSVD(Int_t nrows,Int_t ncols)
 Constructor for (nrows x ncols) matrix
TDecompSVD(Int_t row_lwb, Int_t row_upb, Int_t col_lwb, Int_t col_upb)
 Constructor for ([row_lwb..row_upb] x [col_lwb..col_upb]) matrix
TDecompSVD(const TMatrixD &a,Double_t tol)
 Constructor for general matrix A .
TDecompSVD(const TDecompSVD& another)
 Copy constructor
Bool_t Decompose()
 SVD decomposition of matrix
 If the decomposition succeeds, bit kDecomposed is set , otherwise kSingular
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 iterative 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)
 Step 1 in the matrix diagonalization
void Diag_2(TVectorD& sDiag, TVectorD& oDiag, Int_t k, Int_t l)
 Step 2 in the matrix diagonalization
void Diag_3(TMatrixD& v, TMatrixD& u, TVectorD& sDiag, TVectorD& oDiag, Int_t k, Int_t l)
 Step 3 in the matrix diagonalization
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)
 Set matrix to be decomposed
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)
 Solve A^T x=b assuming the SVD form of A is stored . Solution returned in b.
Double_t Condition()
 Matrix condition number
void Det(Double_t& d1, Double_t& d2)
 Matrix determinant det = d1*TMath::Power(2.,d2)
Bool_t 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(Bool_t &status)
 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
 Print class members
TDecompSVD & operator=(const TDecompSVD& source)
 Assignment operator
const TMatrixDBase & GetDecompMatrix() const
{ return fU; }
TDecompSVD()
{}
virtual ~TDecompSVD()
{}
Int_t GetNrows() const
{ return fU.GetNrows(); }
Int_t GetNcols() const
{ return fV.GetNcols(); }
const TMatrixD & GetU()
const TMatrixD & GetV()
const TVectorD & GetSig()
Bool_t Solve( TVectorD &b)
Bool_t TransSolve( TVectorD &b)
Bool_t Invert(TMatrixD &inv)