template<class F>
class ROOT::Math::CholeskyDecompGenDim< F >

class to compute the Cholesky decomposition of a matrix

class to compute the Cholesky decomposition of a symmetric positive definite matrix when the dimensionality of the problem is not known at compile time

provides routines to check if the decomposition succeeded (i.e. if matrix is positive definite and non-singular), to solve a linear system for the given matrix and to obtain its inverse

the actual functionality is implemented in templated helper classes which have specializations for dimensions N = 1 to 6 to achieve a gain in speed for common matrix sizes

usage example:

// let m be a symmetric positive definite SMatrix (use type float
// for internal computations, matrix size is 4x4)
CholeskyDecompGenDim<float> decomp(4, m);
// check if the decomposition succeeded
if (!decomp) {
  std::cerr << "decomposition failed!" << std::endl;
} else {
  // let rhs be a vector; we seek a vector x such that m * x = rhs
  decomp.Solve(rhs);
  // rhs now contains the solution we are looking for
 
  // obtain the inverse of m, put it into m itself
  decomp.Invert(m);
}

Definition at line 342 of file CholeskyDecomp.h.

Public Member Functions
template<class M >
	CholeskyDecompGenDim (unsigned N, const M &m, std::vector< F > wksp={})
	perform a Cholesky decomposition

template<typename G >
	CholeskyDecompGenDim (unsigned N, G *m, std::vector< F > wksp={})
	perform a Cholesky decomposition

template<typename G >
bool	getL (G *m) const
	obtain the decomposed matrix L

template<class M >
bool	getL (M &m) const
	obtain the decomposed matrix L

template<typename G >
bool	getLi (G *m) const
	obtain the inverse of the decomposed matrix L

template<class M >
bool	getLi (M &m) const
	obtain the inverse of the decomposed matrix L

template<typename G >
bool	Invert (G *m) const
	place the inverse into m

template<class M >
bool	Invert (M &m) const
	place the inverse into m

bool	ok () const
	returns true if decomposition was successful

	operator bool () const
	returns true if decomposition was successful

std::vector< F >	releaseWorkspace ()
	release the workspace of the decomposition for reuse

template<class V >
bool	Solve (V &rhs) const
	solves a linear system for the given right hand side

Private Attributes
std::vector< F >	fL
	lower triangular matrix L

unsigned	fN = 0
	dimensionality dimensionality of the problem

bool	fOk = false
	flag indicating a successful decomposition

#include <Math/CholeskyDecomp.h>

Constructor & Destructor Documentation

◆ CholeskyDecompGenDim() [1/2]

template<class F >

template<class M >

ROOT::Math::CholeskyDecompGenDim< F >::CholeskyDecompGenDim	(	unsigned	N,
		const M &	m,
		std::vector< F >	wksp = {} )

inline

perform a Cholesky decomposition

perform a Cholesky decomposition of a symmetric positive definite matrix m

this is the constructor to uses with an SMatrix (and objects that behave like an SMatrix in terms of using operator()(int i, int j) for access to elements)

The optional wksp parameter that can be used to avoid (re-)allocations on repeated decompositions of the same size (by passing a workspace of size N * (N + 1), or reusing one returned from releaseWorkspace by a previous decomposition).

Definition at line 369 of file CholeskyDecomp.h.

◆ CholeskyDecompGenDim() [2/2]

template<class F >

template<typename G >

ROOT::Math::CholeskyDecompGenDim< F >::CholeskyDecompGenDim	(	unsigned	N,
		G *	m,
		std::vector< F >	wksp = {} )

inline

perform a Cholesky decomposition

perform a Cholesky decomposition of a symmetric positive definite matrix m

this is the constructor to use in special applications where plain arrays are used

NOTE: the matrix is given in packed representation, matrix element m(i,j) (j <= i) is supposed to be in array element (i * (i + 1)) / 2 + j

The optional wksp parameter that can be used to avoid (re-)allocations on repeated decompositions of the same size (by passing a workspace of size N * (N + 1), or reusing one returned from releaseWorkspace by a previous decomposition).

Definition at line 394 of file CholeskyDecomp.h.

Member Function Documentation

◆ getL() [1/2]

template<class F >

template<typename G >

bool ROOT::Math::CholeskyDecompGenDim< F >::getL ( G * m ) const

inline

obtain the decomposed matrix L

Returns: if the decomposition was successful

NOTE: the matrix is given in packed representation, matrix element m(i,j) (j <= i) is supposed to be in array element (i * (i + 1)) / 2 + j

Definition at line 542 of file CholeskyDecomp.h.

◆ getL() [2/2]

template<class F >

template<class M >

bool ROOT::Math::CholeskyDecompGenDim< F >::getL ( M & m ) const

inline

obtain the decomposed matrix L

This is the method to use with a plain array.

Returns: if the decomposition was successful

Definition at line 514 of file CholeskyDecomp.h.

◆ getLi() [1/2]

template<class F >

template<typename G >

bool ROOT::Math::CholeskyDecompGenDim< F >::getLi ( G * m ) const

inline

obtain the inverse of the decomposed matrix L

Returns: if the decomposition was successful

NOTE: the matrix is given in packed representation, matrix element m(j,i) (j <= i) is supposed to be in array element (i * (i + 1)) / 2 + j

Definition at line 598 of file CholeskyDecomp.h.

◆ getLi() [2/2]

template<class F >

template<class M >

bool ROOT::Math::CholeskyDecompGenDim< F >::getLi ( M & m ) const

inline

obtain the inverse of the decomposed matrix L

This is the method to use with a plain array.

Returns: if the decomposition was successful

Definition at line 564 of file CholeskyDecomp.h.

◆ Invert() [1/2]

template<class F >

template<typename G >

bool ROOT::Math::CholeskyDecompGenDim< F >::Invert ( G * m ) const

inline

place the inverse into m

This is the method to use with a plain array.

Returns: if the decomposition was successful

NOTE: the matrix is given in packed representation, matrix element m(i,j) (j <= i) is supposed to be in array element (i * (i + 1)) / 2 + j

Definition at line 495 of file CholeskyDecomp.h.

◆ Invert() [2/2]

template<class F >

template<class M >

bool ROOT::Math::CholeskyDecompGenDim< F >::Invert ( M & m ) const

inline

place the inverse into m

This is the method to use with an SMatrix.

Returns: if the decomposition was successful

Definition at line 474 of file CholeskyDecomp.h.

◆ ok()

template<class F >

bool ROOT::Math::CholeskyDecompGenDim< F >::ok ( ) const

inline

returns true if decomposition was successful

Returns: true if decomposition was successful

Definition at line 445 of file CholeskyDecomp.h.

◆ operator bool()

template<class F >

ROOT::Math::CholeskyDecompGenDim< F >::operator bool ( ) const

inline

returns true if decomposition was successful

Returns: true if decomposition was successful

Definition at line 448 of file CholeskyDecomp.h.

◆ releaseWorkspace()

template<class F >

std::vector< F > ROOT::Math::CholeskyDecompGenDim< F >::releaseWorkspace ( )

inline

release the workspace of the decomposition for reuse

If you use CholeskyDecompGenDim repeatedly with the same size N, it is possible to avoid repeatedly allocating the internal workspace. The constructors take an optional wksp argument, and this routine can be called to get the workspace out of a decomposition when you are done with it.

Please note that once you call releaseWorkspace, you cannot do further operations on the decomposition object (and this is indicated by having it return false when you call ok()).

Here is an example snippet of (pseudo-)code which illustrates the use of releaseWorkspace for inverting the covariance matrices of a number of items (which must all be of the same size):

std::vector<float> wksp;
for (const auto& item: items) {
    CholeskyDecompGenDim decomp(item.cov().size(), item.cov(),
                                std::move(wksp));
    if (decomp.ok()) {
        decomp.Invert(item.invcov());
    }
    wksp = std::move(decomp.releaseWorkspace());
}

In the above code snippet, only the first CholeskyDecompGenDim constructor call will allocate memory. Subsequent calls in the loop will reuse the buffer from the first iteration.

Definition at line 436 of file CholeskyDecomp.h.

◆ Solve()

template<class F >

template<class V >

bool ROOT::Math::CholeskyDecompGenDim< F >::Solve ( V & rhs ) const

inline

solves a linear system for the given right hand side

Note that you can use both SVector classes and plain arrays for rhs. (Make sure that the sizes match!). It will work with any vector implementing the operator [i]

Returns: if the decomposition was successful

Definition at line 459 of file CholeskyDecomp.h.

Member Data Documentation

◆ fL

template<class F >

std::vector<F> ROOT::Math::CholeskyDecompGenDim< F >::fL

private

lower triangular matrix L

lower triangular matrix L, packed storage, with diagonal elements pre-inverted

Definition at line 350 of file CholeskyDecomp.h.

◆ fN

template<class F >

unsigned ROOT::Math::CholeskyDecompGenDim< F >::fN = 0

private

dimensionality dimensionality of the problem

Definition at line 346 of file CholeskyDecomp.h.

◆ fOk

template<class F >

bool ROOT::Math::CholeskyDecompGenDim< F >::fOk = false

private

flag indicating a successful decomposition

Definition at line 352 of file CholeskyDecomp.h.

math/smatrix/inc/Math/CholeskyDecomp.h

Public Member Functions

Private Attributes

Constructor & Destructor Documentation

◆ CholeskyDecompGenDim() [1/2]

◆ CholeskyDecompGenDim() [2/2]

Member Function Documentation

◆ getL() [1/2]

◆ getL() [2/2]

◆ getLi() [1/2]

◆ getLi() [2/2]

◆ Invert() [1/2]

◆ Invert() [2/2]

◆ ok()

◆ operator bool()

◆ releaseWorkspace()

◆ Solve()

Member Data Documentation

◆ fL

◆ fN

◆ fOk