ZOOM ZOOM

The ZOOM CovMatrices package

The ZOOM CovMatrices Package is intended as a lean and mean set of classes representing The goal is to provide classes suitable for high-speed tracking computation.

Classes


CovMatrix Class 
SPD Matrices  1-Column Matrices  2-Column Matrices 
CovMatrix2 
Vector2 
BiVector2 
CovMatrix3 
Vector3 
BiVector3 
CovMatrix4 
Vector4 
BiVector4 
CovMatrix5 
Vector5 
BiVector5 
CovMatrix6 
Vector6 
BiVector6 
CovMatrix 
BiVector 

* std::vector < double > is used to represent a 1-column matrix of arbitrary length.

Constructors

A CovMatrixN may be constructed supplying its (N(N+1)/2) distinct elements in row (or equivalently column) order. The usual copy constructor is present, and the default constructor constructs a unit matrix.

A Vector may be constructed from from its elements, or from an array conaining them.

A BiVector may be constructed from two vectors, or from its elements in the order of first one vector than the other.

Operations

Assuming that M is a CovMatrix of some size, and v and w are respectively 1 and 2-column matrices of that same size, the following operations are supported, and are fast:
  M = M1 + M2;
  bool OK = M.invert();
  double d = M.determinant();
  Vector5 eigenvals; Vector[5] eigenvecs[5]; // here we assume M is 5x5
  M.diagonalize(eigenvals, eigenvecs);
  M.trace();
  cout << M; // nicely formatted output
  M = M1 - M2; // See note
  M*v;
  M*w;
  v*M;
  w*M;
  v+v; v-v;
  w+w; w-w;
Note that unless user knows M2 is small, the answer to M = M1 - M2 may not be positive definite. Nonetheless this operation is supplied because there are tracking algorithms that need to do this as one step, and which know that M2 is supposed to be small.

Element access

The user can do M (i, j) to get individual elements (indices are 0-based).

All elements of the upper part of the matrix, in row order, may be dumped into an array of doubles.

  CovMatrix5 M;
  double elements[15];
  M.dump(elements);

Speed

The following tables compare the speed of CovMatrix operations, with the analogous operations in the CLHEP package (CLHEP 1.8) and in LinearAlgebra (LinearAlgebra 2.1). Note that the routines in CLHEP beyond 1.7 already have had considerable optimization work, so improvements in speed are very meaningful. These improvements stem from explointing the positive definite properties of the matrix, and from tweaking until the optimum algorithm was found.

The bottom line is that for 5x5 and 6x6 matrices, inversion proceeds almost twice as fast using the CovMatrices package as it does using CLHEP Matrices.

Performing matrix addition timings table  (more condensed form of the same table: only the minimum values are reported)
Performing matrix subtraction timings table  (more condensed form of the same table: only the minimum values are reported)
Performing matrix to matrix multiplication timings table  (more condensed form of the same table: only the minimum values are reported)
  (Note: matrix*matrix can't be performed in CovMatrices)
Performing matrix to vector multiplication timings table  (more condensed form of the same table: only the minimum values are reported)
Performing matrix to bivector multiplication timings table (more condensed form of the same table: only the minimum values are reported)
  (Note: since matrix*bivector can't be performed in LinearAlgebra and CLHEP, matrix*matrix with a 2-column matrix as the 2nd matrix was done)
Performing vector to vector multiplication timings table  (more condensed form of the same table: only the minimum values are reported)
  (Note: in order to perform vector*vector in CLHEP, we have to transpose one of the two vectors. However, when we transpose a HepVector we get back a HepMatrix. Thus, we end up performing vector*matrix multiplication and not vector*vector.)
Performing matrix inversion timings table  (more condensed form of the same table: only the minimum values are reported)
Performing random matrix inversion timings table  (more condensed form of the same table: only the minimum values are reported)
Calculating matrix determinant timings table  (more condensed form of the same table: only the minimum values are reported)
Calculating matrix trace timings table  (more condensed form of the same table: only the minimum values are reported)
 

Correctness Validation

Since the intent is that codes which may already be debugged consider moving to use this package for speed, absolute confidence in correctness is a must.
These test programs verify that all operations have been tested and behave properly.

The following pages compare the results of various matrix operations in CovMatrices with the same operations in LinearAlgebra and CLHEP.
The link on each matrix operation leads to a comparison summary for the matrix operation, whereas the link on each package leads to the actual test.

Matrix additionLinearAlgebra, CLHEP, CLHEP (Symmetric matrices), CovMatrices
Matrix subtractionLinearAlgebra, CLHEP, CLHEP (Symmetric matrices), CovMatrices
Matrix * Matrix multiplication: LinearAlgebra, CLHEP, CLHEP (Symmetric matrices), CovMatrices (dot product sums) (*)
Matrix * Vector multiplication: LinearAlgebra, CLHEP, CLHEP (Symmetric matrices), CovMatrices
Vector * Matrix multiplication: LinearAlgebra, CLHEP (**), CLHEP (Symmetric matrices) (**), CovMatrices
(Column)Vector * (Row)Vector multiplication: LinearAlgebra, CLHEP (***), CovMatrices (dot product) (*)
(Row)Vector * (Column)Vector multiplication: LinearAlgebra, CLHEP (***), CovMatrices (dot product) (*)
Matrix traceLinearAlgebra, CLHEP, CLHEP (Symmetric matrices), CovMatrices
Matrix determinantLinearAlgebra, CLHEP, CLHEP (Symmetric matrices), CovMatrices
Matrix inversionLinearAlgebra, CLHEP, CLHEP (Symmetric matrices), CovMatrices

(*) Dot products do not test the ability of the package itself to perform the clculations, they don't use/test a method native to the package; they are provided for symmetry's sake.
(**) In order to perform vector*matrix in CLHEP we have to transpose the vector. However, when we transpose a HepVector we get back a HepMatrix. Thus, we end up performing matrix*matrix multiplication and not vector*matrix.
(***) In order to perform vector*vector in CLHEP we have to transpose one of the two vectors. However, when we transpose a HepVector we get back a HepMatrix. Thus, we end up performing vector*matrix multiplication and not vector*vector.
 
 


ZOOM Home Page - Fermilab at Work - Fermilab Home


David SachsNick Macks,
Last modified: September 16, 2002