 |
ROOT
Reference Guide |
|
| |
ROOT
Reference Guide |
|
The ROOT linear algebra package provides a complete environment in ROOT to perform matrix calculations such as matrix-vector and matrix-matrix multiplications and other linear algebra calculations like equation solving and eigenvalue decompositions.
The present package implements all the basic algorithms dealing with vectors, matrices, matrix columns, rows, diagonals, etc. In addition eigen-Vector analysis and several matrix decomposition have been added (LU,QRH,Cholesky,Bunch-Kaufman and SVD) . The decompositions are used in matrix inversion, equation solving.
ROOT provides the following matrix classes, among others:
TMatrixDBaseTMatrixFTMatrixFSymTVectorFTMatrixDTMatrixDSymTMatrixDSparseTDecompBaseTDecompCholFor a dense matrix, elements are arranged in memory in a ROW-wise fashion . For (n x m) matrices where n*m <=kSizeMax (=25 currently) storage space is available on the stack, thus avoiding expensive allocation/deallocation of heap space . However, this introduces of course kSizeMax overhead for each matrix object . If this is an issue recompile with a new appropriate value (>=0) for kSizeMax
Sparse matrices are also stored in row-wise fashion but additional row/column information is stored, see TMatrixTSparse source for additional details .
Another way to assign and store matrix data is through Use see for instance stressLinear.cxx file .
Unless otherwise specified, matrix and vector indices always start with 0, spanning up to the specified limit-1. However, there are constructors to which one can specify aribtrary lower and upper bounds, e.g. TMatrixD m(1,10,1,5) defines a matrix that ranges from 1..10, 1..5 (a(1,1)..a(10,5)).
A matrix has five properties, which are all set in the constructor:
precision precision is float (i.e. single precision), use the TMatrixF class family. If the precision is double, use the TMatrixD class family.typegeneral (TMatrixD), symmetric (TMatrixDSym) or sparse (TMatrixDSparse).sizeindexsparse mapUse one of the following methods to access the information about the relevant matrix property:
Int_t GetRowLwb(): Row lower-bound index.Int_t GetRowUpb(): Row upper-bound index.Int_t GetNrows(): Number of rows.Int_t GetColLwb(): Column lower-bound index.Int_t GetColUpb(): Column upper-bound index.Int_t GetNcols(): Number of columns.Int_t GetNoElements(): Number of elements, for a dense matrix this equals: fNrows x fNcols.Double_t GetTol(): Tolerance number that is used in decomposition operations.Int_t *GetRowIndexArray(): For sparse matrices, access to the row index of fNrows+1 entries.Int_t *GetColIndexArray(): For sparse matrices, access to the column index of fNelems entries.Use one of the following methods to set a matrix property:
SetTol (Double_t tol)ResizeTo (Int_t nrows,Int_t ncols, Int_t nr_nonzeros=-1)nrows x ncols. Index will start at 0.ResizeTo(Int_t row_lwb,Int_t row_upb, Int_t col_lwb,Int_t col_upb, Int_t nr_nonzeros=-1)row_lwb:row_upb x col_lwb:col_upb.SetRowIndexArray (Int_t *data)fNrows+1 entries column lower-bound index.SetColIndexArray (Int_t *data)fNelems entries.SetSparseIndex (Int_t nelems new)nelems_new elements and copies (if exists) at most nelems_new matrix elements over to the new structure.SetSparseIndex (const TMatrixDBase &a)a.SetSparseIndexAB (const TMatrixDSparse &a, const TMatrixDSparse &b)a and b.Use one of the following constructors to create a matrix:
TMatrixD(Int_t nrows,Int_t ncols)TMatrixD(Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb)TMatrixD(Int_t nrows,Int_t ncols,const Double_t *data, Option_t option= "")TMatrixD(Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb, const Double_t *data,Option_t *option="")TMatrixDSym(Int_t nrows)TMatrixDSym(Int_t row_lwb,Int_t row_upb)TMatrixDSym(Int_t nrows,const Double_t *data,Option_t *option="")TMatrixDSym(Int_t row_lwb,Int_t row_upb, const Double_t *data, Option_t *opt="")TMatrixDSparse(Int_t nrows,Int_t ncols)TMatrixDSparse(Int_t row_lwb,Int_t row_upb,Int_t col_lwb, Int_t col_upb)TMatrixDSparse(Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb, Int_t nr_nonzeros,Int_t *row,Int_t *col,Double_t *data)Use one of the following methods to fill a matrix:
SetMatrixArray(const Double_t*data,Option_t*option="")option="F", the array fills the matrix column-wise else row-wise. This option is implemented for TMatrixD and TMatrixDSym. It is expected that the array data contains at least fNelems entries.SetMatrixArray(Int_t nr,Int_t *irow,Int_t *icol,Double_t *data)nr entries with row index, column index and data entry. Only the entries with non-zero data value are inserted.operator(), operator[]i,j) is found in the sparse index table it will be entered, which involves some memory management. Therefore, before invoking this method in a loop set the index table first through a call to the SetSparseIndex() method.SetSub(Int_t row_lwb,Int_t col_lwb,const TMatrixDBase &source)row_lwb,col_lwb) can be both, dense or sparse.Use()Use() method.Use(TMatrixD &a)Use(Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,Double_t *data)Use(Int_t nrows,Int_t ncols,Double_t *data)Use(TMatrixDSym &a)Use(Int_t nrows,Double_t *data)Use(Int_t row_lwb,Int_t row_upb,Double_t *data)Use(TMatrixDSparse &a)Use(Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,Int_t nr_no nzeros, Int_t *pRowIndex,Int_t *pColIndex,Double_t *pData)Use(Int_t nrows,Int_t ncols,Int_t nr_nonzeros,Int_t *pRowIndex,Int_t *pColIndex,Double_t *pData)Example
A Hilbert matrix is created by copying an array.
You can also assign the data array to the matrix without actually copying it.
The array data now contains the inverted matrix.
Now a unit matrix in sparse format is created.
Invert(Double_t &det=0) function to invert a matrix:– or –
Both methods are available for general and symmetric matrices.
For matrices whose size is less than or equal to 6x6, the InvertFast(Double_t &det=0) function is available. Here the Cramer algorithm will be applied, which is faster but less accurate.
You can also use the following decomposition classes (see → [Matrix decompositions](#matrix-decompositions")) for inverting a matrix: <table width="100%"> <tr> <th scope="col">Name</th> <th scope="col">Matrix type</th> <th scope="col">Comment</th> </tr> <tr> <td>TDecompLU</td> <td>General</td> <td></td> </tr> <tr> <td>TDecompQRH</td> <td>General</td> <td></td> </tr> <tr> <td>TDecompSVD</td> <td>General</td> <td>Can manipulate singular matrix.</td> </tr> <tr> <td>TDecompBK</td> <td>symmetric</td> <td></td> </tr> <tr> <td>TDecompChol</td> <td>Symmetric</td> <td>Matrix should also be positive definite.</td> </tr> <tr> <td>TDecompSparse</td> <td>Sparse</td> <td></td> </tr> </table> If the required matrix type is general, you also can handle symmetric matrices. <em><strong>Example</strong></em> This example shows how to check whether the matrix is singular before attempting to invert it. @code TDecompLU lu(a); TMatrixD b; if (!lu.Decompose()) { cout << "Decomposition failed, matrix singular ?" << endl; cout << "condition number = " << = a.GetCondition() << endl; } else { lu.Invert(b); } @endcode <h3>Matrix operators and methods</h3> The matrix/vector operations are classified according to BLAS (basic linear algebra subroutines) levels. <h4>Arithmetic operations between matrices</h4> <table width="100%"> <tr> <th scope="col">Description</th> <th scope="col">Format</th> <th scope="col">Comment</th> </tr> <tr> <td>Element</td> <td>C=A+B</td> <td>overwrites A</td> </tr> <tr> <td>Wise sum</td> <td>A+=B<br>Add (A,alpha,B) <br>TMatrixD(A,TMatrixD::kPlus,B)</td> <td>A = A + α B constructor</td> </tr> <tr> <td>Element wise substraction</td> <td>C=A-B A-=B<br> TMatrixD(A,TMatrixD::kMinus,B)</td> <td>overwrites A<br> constructor</td> </tr> <tr> <td>Matrix multiplication</td> <td>C=A*B<br> A*=B<br> C.Mult(A,B)<br>TMatrixD(A,TMatrixD::kMult,B)<br>TMatrixD(A, TMatrixD(A, TMatrixD::kTransposeMult,B)<br>TMatrixD(A, TMatrixD::kMultTranspose,B)</td> <td>overwrites A<br> <br> <br>constructor of A.B<br>constructor of A<sup>T</sup> .B<br>constructor of A.B<sup>T</sup></td> </tr> <tr> <td>Element wise multiplication</td> <td>ElementMult(A,B)</td> <td>A(i,j)*= B(i,j)</td> </tr> <tr> <td>Element wise division</td> <td>ElementDiv(A,B)</td> <td>A(i,j)/= B(i,j)</td> </tr> </table> <h4>Arithmetic operations between matrices and real numbers</h4> <table width="100%"> <tr> <th scope="col">Description</th> <th scope="col">Format</th> <th scope="col">Comment</th> </tr> <tr> <td>Element wise sum</td> <td>C=r+A C=A+r A+=r</td> <td>overwrites A</td> </tr> <tr> <td>Element wise subtraction</td> <td>C=r-A C=A-r A-=r</td> <td>overwrites A</td> </tr> <tr> <td>Matrix multiplication</td> <td>C=r*A C=A*r A*=r</td> <td>overwrites A</td> </tr> </table> <h4>Comparison and Boolean operations</h4> <strong>Comparison between two matrices</strong> <table width="100%"> <tr> <th scope="col">Description</th> <th scope="col">Output</th> <th scope="col">Descriptions</th> </tr> <tr> <td>A == B</td> <td>Bool_t</td> <td>equal to</td> </tr> <tr> <td>A != B</td> <td>matrix</td> <td>not equal</td> </tr> <tr> <td>A > B</td> <td>matrix</td> <td>greater than</td> </tr> <tr> <td>A >= B</td> <td>matrix</td> <td>greater than or equal to</td> </tr> <tr> <td>A < B</td> <td>matrix</td> <td>smaller than</td> </tr> <tr> <td>A <= B</td> <td>matrix</td> <td>smaller than or equal to</td> </tr> <tr> <td>AreCompatible(A,B)</td> <td>Bool_t</td> <td>compare matrix properties</td> </tr> <tr> <td>Compare(A,B)</td> <td>Bool_t</td> <td>return summary of comparison</td> </tr> <tr> <td>VerifyMatrixIdentity(A,B,verb, maxDev)</td> <td> </td> <td>check matrix identity within maxDev tolerance</td> </tr> </table> <strong>Comparison between matrix and real number</strong> <table width="100%"> <tr> <th scope="col">Format</th> <th scope="col">Output</th> <th scope="col">Description</th> </tr> <tr> <td>A == r</td> <td>Bool_t</td> <td>equal to</td> </tr> <tr> <td>A != r</td> <td>Bool_t</td> <td>not equal</td> </tr> <tr> <td>A > r</td> <td>Bool_t</td> <td>greater than</td> </tr> <tr> <td>A >= r</td> <td>Bool_t</td> <td>greater than or equal to</td> </tr> <tr> <td>A < r</td> <td>Bool_t</td> <td>smaller than</td> </tr> <tr> <td>A <= r</td> <td>Bool_t</td> <td>smaller than or equal to</td> </tr> <tr> <td>VerifyMatrixValue(A,r,verb, maxDev)</td> <td>Bool_t</td> <td>compare matrix value with r within maxDev tolerance</td> </tr> <tr> <td>A.RowNorm()</td> <td>Double_t</td> <td>norm induced by the infinity vector norm</td> </tr> <tr> <td>A.NormInf()</td> <td>Double_t</td> <td> </td> </tr> <tr> <td>A.ColNorm()</td> <td>Double_t</td> <td>norm induced by the 1 vector norm</td> </tr> <tr> <td>A.Norm1()</td> <td>Double_t</td> <td> </td> </tr> <tr> <td>A.E2Norm()</td> <td>Double_t</td> <td>square of the Euclidean norm</td> </tr> <tr> <td>A.NonZeros()</td> <td>Int_t</td> <td> </td> </tr> <tr> <td>A.Sum()</td> <td>Double_t</td> <td>number of elements unequal zero</td> </tr> <tr> <td>A.Min()</td> <td>Double_t</td> <td> </td> </tr> <tr> <td>A.Max()</td> <td>Double_t</td> <td> </td> </tr> <tr> <td>A.NormByColumn (v,"D")</td> <td>TMatrixD</td> <td> </td> </tr> <tr> <td>A.NormByRow (v,"D")</td> <td>TMatrixD</td> <td> </td> </tr> </table> <h3>Matrix views</h3> With the following matrix view classes, you can access the matrix elements: - <tt>TMatrixDRow</tt> - <tt>TMatrixDColumn</tt> - <tt>TMatrixDDiag</tt> - <tt>TMatrixDSub</tt> <h4>Matrix view operators</h4> For the matrix view classes <tt>TMatrixDRow</tt>, <tt>TMatrixDColumn</tt> and <tt>TMatrixDDiag</tt>, the necessary assignment operators are available to interact with the vector class <tt>TVectorD</tt>.<br>The sub matrix view classes <tt>TMatrixDSub</tt> has links to the matrix classes <tt>TMatrixD</tt> and <tt>TMatrixDSym.</tt> The next table summarizes how to access the individual matrix elements in the matrix view classes. <table width="100%"> <tr> <th scope="col">Format</th> <th scope="col">Comment</th> </tr> <tr> <td>TMatrixDRow(A,i)(j) TMatrixDRow(A,i)[j]</td> <td>element A<sub>ij</sub></td> </tr> <tr> <td>TMatrixDColumn(A,j)(i) TMatrixDColumn(A,j)[i]</td> <td>element A<sub>ij</sub></td> </tr> <tr> <td>TMatrixDDiag(A(i) TMatrixDDiag(A[i]</td> <td>element A<sub>ij</sub></td> </tr> <tr> <td>TMatrixDSub(A(i) TMatrixDSub(A,rl,rh,cl,ch)(i,j)</td> <td>element A<sub>ij</sub><br>element A<sub>rl+i,cl+j</sub></td> </tr> </table> <h4>Matrix decompositions</h4> There are the following classes available for matrix decompositions: - {% include ref class="TDecompLU" %}: Decomposes a general <tt>n x n</tt> matrix <tt>A</tt> into <tt>P A = L U</tt>. - {% include ref class="TDecompBK" %}: The Bunch-Kaufman diagonal pivoting method decomposes a real symmetric matrix <tt>A</tt>. - {% include ref class="TDecompChol" %}: The Cholesky decomposition class, which decomposes a symmetric, positive definite matrix <tt>A = U^T * U</tt> where <tt>U</tt> is a upper triangular matrix. - {% include ref class="TDecompQRH" %}: QR decomposition class. - {% include ref class="TDecompSVD" %}: Single value decomposition class. - {% include ref class="TDecompSparse" %}: Sparse symmetric decomposition class. <h3>Matrix Eigen analysis</h3> With the <tt>TMatrixDEigen</tt> and <tt>TMatrixDSymEigen</tt> classes, you can compute eigenvalues and eigenvectors for general dense and symmetric real matrices. <h2>Additonal Notes</h2> The present package provides all facilities to completely AVOID returning matrices. Use "TMatrixD A(TMatrixD::kTransposed,B);" and other fancy constructors as much as possible. If one really needs to return a matrix, return a TMatrixTLazy object instead. The conversion is completely transparent to the end user, e.g. "TMatrixT m = THaarMatrixT(5);" and <em>is</em> efficient. Since TMatrixT et al. are fully integrated in ROOT, they of course can be stored in a ROOT database. <h3>How to efficiently use this package</h3> <h4>1. Never return complex objects (matrices or vectors)</h4> Danger: For example, when the following snippet: @code TMatrixD foo(int n) { TMatrixD foom(n,n); fill_in(foom); return foom; } TMatrixD m = foo(5); @endcode runs, it constructs matrix foo:foom, copies it onto stack as a return value and destroys foo:foom. Return value (a matrix) from foo() is then copied over to m (via a copy constructor), and the return value is destroyed. So, the matrix constructor is called 3 times and the destructor 2 times. For big matrices, the cost of multiple constructing/copying/destroying of objects may be very large. <em>Some</em> optimized compilers can cut down on 1 copying/destroying, but still it leaves at least two calls to the constructor. Note, TMatrixDLazy (see below) can construct TMatrixD m "inplace", with only a <em>single</em> call to the constructor. <h4>2. Use "two-address instructions"</h4> @code "void TMatrixD::operator += (const TMatrixD &B);" @endcode as much as possible. That is, to add two matrices, it's much more efficient to write @code A += B; @endcode than @code TMatrixD C = A + B; @endcode (if both operand should be preserved, <tt>TMatrixD C = A; C += B;</tt> is still better). <h4>3. Use glorified constructors when returning of an object seems inevitable:</h4> @code "TMatrixD A(TMatrixD::kTransposed,B);" "TMatrixD C(A,TMatrixD::kTransposeMult,B);" @endcode like in the following snippet (from $ROOTSYS/test/vmatrix.cxx) that verifies that for an orthogonal matrix T, T'T = TT' = E. @code TMatrixD haar = THaarMatrixD(5); TMatrixD unit(TMatrixD::kUnit,haar); TMatrixD haar_t(TMatrixD::kTransposed,haar); TMatrixD hth(haar,TMatrixD::kTransposeMult,haar); TMatrixD hht(haar,TMatrixD::kMult,haar_t); TMatrixD hht1 = haar; hht1 *= haar_t; VerifyMatrixIdentity(unit,hth); VerifyMatrixIdentity(unit,hht); VerifyMatrixIdentity(unit,hht1); @endcode <h4>4. Accessing row/col/diagonal of a matrix without much fuss</h4> (and without moving a lot of stuff around): @code TMatrixD m(n,n); TVectorD v(n); TMatrixDDiag(m) += 4; v = TMatrixDRow(m,0); TMatrixDColumn m1(m,1); m1(2) = 3; // the same as m(2,1)=3; @endcode Note, constructing of, say, TMatrixDDiag does <em>not</em> involve any copying of any elements of the source matrix. <h4>5. It's possible (and encouraged) to use "nested" functions</h4> For example, creating of a Hilbert matrix can be done as follows: @code void foo(const TMatrixD &m) { TMatrixD m1(TMatrixD::kZero,m); struct MakeHilbert : public TElementPosActionD { void Operation(Double_t &element) { element = 1./(fI+fJ-1); } }; m1.Apply(MakeHilbert()); } @endcode of course, using a special method THilbertMatrixD() is still more optimal, but not by a whole lot. And that's right, class MakeHilbert is declared <em>within</em> a function and local to that function. It means one can define another MakeHilbert class (within another function or outside of any function, that is, in the global scope), and it still will be OK. Note, this currently is not yet supported by the interpreter CINT. Another example is applying of a simple function to each matrix element: @code void foo(TMatrixD &m,TMatrixD &m1) { typedef double (*dfunc_t)(double); class ApplyFunction : public TElementActionD { dfunc_t fFunc; void Operation(Double_t &element) { element=fFunc(element); } public: ApplyFunction(dfunc_t func):fFunc(func) {} }; ApplyFunction x(TMath::Sin); m.Apply(x); } @endcode Validation code $ROOTSYS/test/vmatrix.cxx and vvector.cxx contain a few more examples of that kind. <h4>6. Lazy matrices:</h4> instead of returning an object return a "recipe" how to make it. The full matrix would be rolled out only when and where it's needed: @code TMatrixD haar = THaarMatrixD(5); @endcode THaarMatrixD() is a <em>class</em>, not a simple function. However similar this looks to a returning of an object (see note #1 above), it's dramatically different. THaarMatrixD() constructs a TMatrixDLazy, an object of just a few bytes long. A special "TMatrixD(const TMatrixDLazy &recipe)" constructor follows the recipe and makes the matrix haar() right in place. No matrix element is moved whatsoever! <h3>Acknowledgements</h3> 1. Oleg E. Kiselyov First implementations were based on the his code . We have diverged quite a bit since then but the ideas/code for lazy matrix and "nested function" are 100% his . You can see him and his code in action at http://okmij.org/ftp
