Logo ROOT  
Reference Guide
Matrix Linear Algebra

The ROOT Matrix Linear Algebra package.

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.

Matrix classes

ROOT provides the following matrix classes, among others:

For 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 arbitrary 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)).

Matrix properties

A matrix has five properties, which are all set in the constructor:

Accessing matrix properties

Use one of the following methods to access the information about the relevant matrix property:

Setting matrix properties

Use one of the following methods to set a matrix property:

Creating and filling a matrix

Use one of the following constructors to create a matrix:

Use one of the following methods to fill a matrix:

Example

A Hilbert matrix is created by copying an array.

TMatrixD h(5,5);
TArrayD data(25);
for (Int_t = 0; i < 25; i++) {
const Int_t ir = i/5;
const Int_t ic = i%5;
data[i] = 1./(ir+ic);
}
h.SetMatrixArray(data.GetArray());

You can also assign the data array to the matrix without actually copying it.

TMatrixD h; h.Use(5,5,data.GetArray());
h.Invert();

The array data now contains the inverted matrix.

Now a unit matrix in sparse format is created.

TMatrixDSparse unit1(5,5);
TArrayI row(5),col(5);
for (Int_t i = 0; i < 5; i++) row[i] = col[i] = i;
TArrayD data(5); data.Reset(1.);
unit1.SetMatrixArray(5,row.GetArray(),col.GetArray(),data.GetArray());
TMatrixDSparse unit2(5,5);
unit2.SetSparseIndex(5);
unit2.SetRowIndexArray(row.GetArray());
unit2.SetColIndexArray(col.GetArray());
unit2.SetMatrixArray(data.GetArray());

Inverting a matrix

TMatrixD a(...);
a.Invert();

– or –

TMatrixD b(kInvert,a);

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.

Using decomposition classes for inverting

You can also use the following decomposition classes (see → Matrix decompositions) for inverting a matrix:

Name Matrix Comment
TDecompLU General
TDecompQRH General
TDecompSVD General Can manipulate singular matrix.
TDecompBK Symmetric
TDecompChol Symmetric Matrix should also be positive definite.
TDecompSparse Sparse

If the required matrix type is general, you also can handle symmetric matrices.

Example

This example shows how to check whether the matrix is singular before attempting to invert it.

TDecompLU lu(a);
TMatrixD b;
if (!lu.Decompose()) {
cout << "Decomposition failed, matrix singular ?" << endl;
cout << "condition number = " << = a.GetCondition() << endl;
} else {
lu.Invert(b);
}

Matrix operators and methods

The matrix/vector operations are classified according to BLAS (basic linear algebra subroutines) levels.

Arithmetic operations between matrices

Description Format Comment
Element C=A+B overwrites A
Wise sum A+=B
Add (A,alpha,B)
TMatrixD(A,TMatrixD::kPlus,B)
A = A + α B constructor
Element wise subtraction C=A-B A-=B
TMatrixD(A,TMatrixD::kMinus,B)
overwrites A
constructor
Matrix multiplication C=A*B
A*=B
C.Mult(A,B)
TMatrixD(A,TMatrixD::kMult,B)
TMatrixD(A, TMatrixD(A, TMatrixD::kTransposeMult,B)
TMatrixD(A, TMatrixD::kMultTranspose,B)
overwrites A


constructor of A·B
constructor of AT·B
constructor of A·BT
Element wise multiplication ElementMult(A,B) A(i,j)*= B(i,j)
Element wise division ElementDiv(A,B) A(i,j)/= B(i,j)

Arithmetic operations between matrices and real numbers

Description Format Comment
Element wise sum C=r+A C=A+r A+=r overwrites A
Element wise subtraction C=r-A C=A-r A-=r overwrites A
Matrix multiplication C=r*A C=A*r A*=r overwrites A

Comparison and Boolean operations

Comparison between two matrices

Description Output Descriptions
A == B Bool_t equal to
A != B matrix not equal
A > B matrix greater than
A >= B matrix greater than or equal to
A < B matrix smaller than
A <= B matrix smaller than or equal to
AreCompatible(A,B) Bool_t compare matrix properties
Compare(A,B) Bool_t return summary of comparison
VerifyMatrixIdentity(A,B,verb, maxDev) check matrix identity within maxDev tolerance

Comparison between matrix and real number

Format Output Description
A == r Bool_t equal to
A != r Bool_t not equal
A > r Bool_t greater than
A >= r Bool_t greater than or equal to
A < r Bool_t smaller than
A <= r Bool_t smaller than or equal to
VerifyMatrixValue(A,r,verb, maxDev) Bool_t compare matrix value with r within maxDev tolerance
A.RowNorm() Double_t norm induced by the infinity vector norm
A.NormInf() Double_t
A.ColNorm() Double_t norm induced by the 1 vector norm
A.Norm1() Double_t
A.E2Norm() Double_t square of the Euclidean norm
A.NonZeros() Int_t
A.Sum() Double_t number of elements unequal zero
A.Min() Double_t
A.Max() Double_t
A.NormByColumn (v,"D") TMatrixD
A.NormByRow (v,"D") TMatrixD

Matrix views

With the following matrix view classes, you can access the matrix elements:

Matrix view operators

For the matrix view classes TMatrixDRow, TMatrixDColumn and TMatrixDDiag, the necessary assignment operators are available to interact with the vector class TVectorD.
The sub matrix view classes TMatrixDSub has links to the matrix classes TMatrixD and TMatrixDSym.

The next table summarizes how to access the individual matrix elements in the matrix view classes.

Format Comment
TMatrixDRow(A,i)(j) TMatrixDRow(A,i)[j] element Aij
TMatrixDColumn(A,j)(i) TMatrixDColumn(A,j)[i] element Aij
TMatrixDDiag(A(i) TMatrixDDiag(A[i] element Aij
TMatrixDSub(A(i) TMatrixDSub(A,rl,rh,cl,ch)(i,j) element Aij
element Arl+i,cl+j

Matrix decompositions

There are the following classes available for matrix decompositions:

Matrix Eigen analysis

With the TMatrixDEigen and TMatrixDSymEigen classes, you can compute eigenvalues and eigenvectors for general dense and symmetric real matrices.

Additional Notes

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 is efficient.

Since TMatrixT et al. are fully integrated in ROOT, they of course can be stored in a ROOT database.

How to efficiently use this package

1. Never return complex objects (matrices or vectors)

Danger: For example, when the following snippet:

TMatrixD foo(int n)
{
TMatrixD foom(n,n); fill_in(foom); return foom;
}
TMatrixD m = foo(5);

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. Some 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 single call to the constructor.

2. Use "two-address instructions"

void TMatrixD::operator += (const TMatrixD &B);

as much as possible. That is, to add two matrices, it's much more efficient to write

A += B;

than

TMatrixD C = A + B;

(if both operand should be preserved, TMatrixD C = A; C += B; is still better).

3. Use glorified constructors when returning of an object seems inevitable:

TMatrixD A(TMatrixD::kTransposed,B);
TMatrixD C(A,TMatrixD::kTransposeMult,B);

like in the following snippet (from $ROOTSYS/test/vmatrix.cxx) that verifies that for an orthogonal matrix T, T'T = TT' = E.

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

4. Accessing row/col/diagonal of a matrix without much fuss

(and without moving a lot of stuff around):

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;

Note, constructing of, say, TMatrixDDiag does not involve any copying of any elements of the source matrix.

5. It's possible (and encouraged) to use "nested" functions

For example, creating of a Hilbert matrix can be done as follows:

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

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 within 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:

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

Validation code $ROOTSYS/test/vmatrix.cxx and vvector.cxx contain a few more examples of that kind.

6. Lazy matrices:

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:

TMatrixD haar = THaarMatrixD(5);

THaarMatrixD() is a class, 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!

Acknowledgements

  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
  2. Chris R. Birchenhall, We adapted his idea of the implementation for the decomposition classes instead of our messy installation of matrix inversion His installation of matrix condition number, using an iterative scheme using the Hage algorithm is worth looking at ! Chris has a nice writeup (matdoc.ps) on his matrix classes at ftp://ftp.mcc.ac.uk/pub/matclass/
  3. Mark Fischler and Steven Haywood of CLHEP They did the slave labor of spelling out all sub-determinants for Cramer inversion of (4x4),(5x5) and (6x6) matrices The stack storage for small matrices was also taken from them
  4. Roldan Pozo of TNT (http://math.nist.gov/tnt/) He converted the EISPACK routines for the eigen-vector analysis to C++ . We started with his implementation
  5. Siegmund Brandt (http://siux00.physik.uni-siegen.de/~brandt/datan We adapted his (very-well) documented SVD routines

Classes

class  TDecompBase
 Decomposition Base class. More...
 
class  TDecompBK
 The Bunch-Kaufman diagonal pivoting method decomposes a real symmetric matrix A using. More...
 
class  TDecompChol
 Cholesky Decomposition class. More...
 
class  TDecompLU
 LU Decomposition class. More...
 
class  TDecompQRH
 QR Decomposition class. More...
 
class  TDecompSparse
 Sparse Symmetric Decomposition class. More...
 
class  TDecompSVD
 Single Value Decomposition class. More...
 
class  TMatrixDEigen
 TMatrixDEigen. More...
 
class  TMatrixDSymEigen
 TMatrixDSymEigen. More...
 
class  TMatrixT< Element >
 TMatrixT. More...
 
class  TMatrixTBase< Element >
 TMatrixTBase. More...
 
class  TMatrixTLazy< Element >
 Templates of Lazy Matrix classes. More...
 
class  TMatrixTSparse< Element >
 TMatrixTSparse. More...
 
class  TMatrixTSym< Element >
 TMatrixTSym. More...
 
class  TVectorT< Element >
 TVectorT. More...