Updateable Cholesky decomposition: A=LDL'. More...

#include <vnl_ldl_cholesky.h>

Public Types
enum	Operation { quiet, verbose, estimate_condition }
	Modes of computation. See constructor for details. More...

Public Member Functions
	vnl_ldl_cholesky (vnl_matrix< double > const &M, Operation mode=verbose)
	Make cholesky decomposition of M optionally computing the reciprocal condition number. More...

	~vnl_ldl_cholesky ()=default

vnl_vector< double >	solve (vnl_vector< double > const &b) const
	Solve LS problem M x = b. More...

void	solve (vnl_vector< double > const &b, vnl_vector< double > *x) const
	Solve LS problem M x = b. More...

void	solve_lx (vnl_vector< double > &y)
	Solve equation of form Lx=y (in-place). More...

double	determinant () const
	Compute determinant. More...

void	rank1_update (const vnl_vector< double > &v)
	Compute rank-1 update, ie the decomposition of (M+v.v'). More...

void	update (const vnl_matrix< double > &W)
	Multi-rank update, ie the decomposition of (M+W.W'). More...

vnl_matrix< double >	inverse () const
	Compute inverse. Not efficient. More...

const vnl_matrix< double > &	lower_triangle () const
	Return lower-triangular factor. More...

vnl_matrix< double >	upper_triangle () const
	Return upper-triangular factor. More...

const vnl_vector< double > &	diagonal () const
	Return elements of diagonal matrix D in LDL'. More...

double	xt_m_inv_x (const vnl_vector< double > &x) const
	Efficient computation of x' * inv(M) * x. More...

double	xt_m_x (const vnl_vector< double > &x) const
	Efficient computation of x' * M * x. More...

int	rank_deficiency () const
	A Success/failure flag. More...

double	rcond () const
	Return reciprocal condition number (smallest/largest singular values). More...

vnl_vector< double > &	nullvector ()
	Return computed nullvector. More...

vnl_vector< double > const &	nullvector () const

Protected Attributes
vnl_matrix< double >	L_
	Lower triangular matrix. More...

vnl_vector< double >	d_
	Elements of diagonal matrix. More...

double	rcond_
	1/(condition number). More...

long	num_dims_rank_def_

vnl_vector< double >	nullvector_

Private Member Functions
	vnl_ldl_cholesky (vnl_ldl_cholesky const &that)=delete
	Copy constructor - privatised to avoid it being used. More...

vnl_ldl_cholesky &	operator= (vnl_ldl_cholesky const &that)=delete
	Assignment operator - privatised to avoid it being used. More...

void	inplace_solve (double *x) const
	Solve Mx=b, overwriting input vector with the solution. More...

Detailed Description

Updateable Cholesky decomposition: A=LDL'.

A class to hold the Cholesky decomposition of a positive definite symmetric matrix, using the form A=LDL', where L is lower triangular with ones on the leading diagonal, and D is a diagonal matrix. This differs from vnl_cholesky, which decomposes as A=LL', without a constraint on the diagonal. The advantage of the LDL' form is that it can be efficiently updated. Thus given the decomposition of A, we can compute the decomposition of (A+v.v') in O(n^2) time.

This can be used to solve linear systems, compute determinants and inverses.

To check that the decomposition can be used safely for solving a linear equation it is wise to construct with mode==estimate_condition and check that rcond()>sqrt(machine precision). If this is not the case it might be a good idea to use vnl_svd instead.

Definition at line 29 of file vnl_ldl_cholesky.h.

Member Enumeration Documentation

◆ Operation

enum vnl_ldl_cholesky::Operation

Modes of computation. See constructor for details.

Enumerator
quiet
verbose
estimate_condition

Definition at line 33 of file vnl_ldl_cholesky.h.

Constructor & Destructor Documentation

◆ vnl_ldl_cholesky() [1/2]

vnl_ldl_cholesky::vnl_ldl_cholesky	(	vnl_matrix< double > const &	M,
		Operation	mode = `verbose`
	)

Make cholesky decomposition of M optionally computing the reciprocal condition number.

Cholesky decomposition.

Make cholesky decomposition of M optionally computing the reciprocal condition number. If mode is estimate_condition, the condition number and an approximate nullspace are estimated, at a cost of a factor of (1 + 18/n). Here's a table of 1 + 18/n:

 n:              3      5     10     50    100    500   1000
 slowdown:     7.0    4.6    2.8    1.4   1.18   1.04   1.02

Sqrt of elements of diagonal matrix.

Definition at line 25 of file vnl_ldl_cholesky.cxx.

◆ ~vnl_ldl_cholesky()

vnl_ldl_cholesky::~vnl_ldl_cholesky ( )

default

◆ vnl_ldl_cholesky() [2/2]

vnl_ldl_cholesky::vnl_ldl_cholesky ( vnl_ldl_cholesky const & that )

privatedelete

Copy constructor - privatised to avoid it being used.

Member Function Documentation

◆ determinant()

double vnl_ldl_cholesky::determinant ( ) const

Compute determinant.

Definition at line 177 of file vnl_ldl_cholesky.cxx.

◆ diagonal()

const vnl_vector<double>& vnl_ldl_cholesky::diagonal ( ) const

inline

Return elements of diagonal matrix D in LDL'.

Definition at line 78 of file vnl_ldl_cholesky.h.

◆ inplace_solve()

void vnl_ldl_cholesky::inplace_solve ( double * x ) const

private

Solve Mx=b, overwriting input vector with the solution.

x points to beginning of an n-element vector containing b On exit, x[i] filled with solution vector.

Definition at line 98 of file vnl_ldl_cholesky.cxx.

◆ inverse()

vnl_matrix< double > vnl_ldl_cholesky::inverse ( ) const

Compute inverse. Not efficient.

Note that you rarely need the inverse - backsubstitution is faster and less prone to rounding errors.

Definition at line 249 of file vnl_ldl_cholesky.cxx.

◆ lower_triangle()

const vnl_matrix<double>& vnl_ldl_cholesky::lower_triangle ( ) const

inline

Return lower-triangular factor.

Definition at line 72 of file vnl_ldl_cholesky.h.

◆ nullvector() [1/2]

vnl_vector<double>& vnl_ldl_cholesky::nullvector ( )

inline

Return computed nullvector.

Not calculated unless Operation mode at construction was estimate_condition.

Definition at line 99 of file vnl_ldl_cholesky.h.

◆ nullvector() [2/2]

vnl_vector<double> const& vnl_ldl_cholesky::nullvector ( ) const

inline

Definition at line 100 of file vnl_ldl_cholesky.h.

◆ operator=()

vnl_ldl_cholesky& vnl_ldl_cholesky::operator= ( vnl_ldl_cholesky const & that )

privatedelete

Assignment operator - privatised to avoid it being used.

◆ rank1_update()

void vnl_ldl_cholesky::rank1_update ( const vnl_vector< double > & v )

Compute rank-1 update, ie the decomposition of (M+v.v').

If the initial state is the decomposition of M, then L and D are updated so that on exit LDL'=M+v.v'

Uses the algorithm given by Davis and Hager in "Multiple-Rank Modifications of a Sparse Cholesky Factorization",2001.

Definition at line 191 of file vnl_ldl_cholesky.cxx.

◆ rank_deficiency()

int vnl_ldl_cholesky::rank_deficiency ( ) const

inline

A Success/failure flag.

Definition at line 89 of file vnl_ldl_cholesky.h.

◆ rcond()

double vnl_ldl_cholesky::rcond ( ) const

inline

Return reciprocal condition number (smallest/largest singular values).

As long as rcond()>sqrt(precision) the decomposition can be used for solving equations safely. Not calculated unless Operation mode at construction was estimate_condition.

Definition at line 95 of file vnl_ldl_cholesky.h.

◆ solve() [1/2]

vnl_vector< double > vnl_ldl_cholesky::solve ( vnl_vector< double > const & b ) const

Solve LS problem M x = b.

Solve least squares problem M x = b.

Definition at line 167 of file vnl_ldl_cholesky.cxx.

◆ solve() [2/2]

void vnl_ldl_cholesky::solve	(	vnl_vector< double > const &	b,
		vnl_vector< double > *	xp
	)		const

Solve LS problem M x = b.

Solve least squares problem M x = b.

The right-hand-side std::vector x may be b, which will give a fractional increase in speed.

Definition at line 158 of file vnl_ldl_cholesky.cxx.

◆ solve_lx()

void vnl_ldl_cholesky::solve_lx ( vnl_vector< double > & x )

Solve equation of form Lx=y (in-place).

Solve Lx=y (in-place).

x is overwritten with solution

Definition at line 88 of file vnl_ldl_cholesky.cxx.

◆ update()

void vnl_ldl_cholesky::update ( const vnl_matrix< double > & W0 )

Multi-rank update, ie the decomposition of (M+W.W').

If the initial state is the decomposition of M, then L and D are updated so that on exit LDL'=M+W.W'

Definition at line 216 of file vnl_ldl_cholesky.cxx.

◆ upper_triangle()

vnl_matrix<double> vnl_ldl_cholesky::upper_triangle ( ) const

inline

Return upper-triangular factor.

Definition at line 75 of file vnl_ldl_cholesky.h.

◆ xt_m_inv_x()

double vnl_ldl_cholesky::xt_m_inv_x ( const vnl_vector< double > & x ) const

Efficient computation of x' * inv(M) * x.

Useful when M is a covariance matrix!

Useful when M is a covariance matrix! Solves Ly=x for y, then returns sum y[i]*y[i]/d[i]

Definition at line 121 of file vnl_ldl_cholesky.cxx.

◆ xt_m_x()

double vnl_ldl_cholesky::xt_m_x ( const vnl_vector< double > & x ) const

Efficient computation of x' * M * x.

Twice as fast as explicitly computing x' * M * x

Definition at line 138 of file vnl_ldl_cholesky.cxx.

Member Data Documentation

◆ d_

vnl_vector<double> vnl_ldl_cholesky::d_

protected

Elements of diagonal matrix.

Definition at line 109 of file vnl_ldl_cholesky.h.

◆ L_

vnl_matrix<double> vnl_ldl_cholesky::L_

protected

Lower triangular matrix.

Definition at line 106 of file vnl_ldl_cholesky.h.

◆ nullvector_

vnl_vector<double> vnl_ldl_cholesky::nullvector_

protected

Definition at line 114 of file vnl_ldl_cholesky.h.

◆ num_dims_rank_def_

long vnl_ldl_cholesky::num_dims_rank_def_

protected

Definition at line 113 of file vnl_ldl_cholesky.h.

◆ rcond_

double vnl_ldl_cholesky::rcond_

protected

1/(condition number).

Definition at line 112 of file vnl_ldl_cholesky.h.

The documentation for this class was generated from the following files:

core/vnl/algo/vnl_ldl_cholesky.h
core/vnl/algo/vnl_ldl_cholesky.cxx

Public Types

Public Member Functions

Protected Attributes

Private Member Functions

Detailed Description

Member Enumeration Documentation

◆ Operation

Constructor & Destructor Documentation

◆ vnl_ldl_cholesky() [1/2]

◆ ~vnl_ldl_cholesky()

◆ vnl_ldl_cholesky() [2/2]

Member Function Documentation

◆ determinant()

◆ diagonal()

◆ inplace_solve()

◆ inverse()

◆ lower_triangle()

◆ nullvector() [1/2]

◆ nullvector() [2/2]

◆ operator=()

◆ rank1_update()

◆ rank_deficiency()

◆ rcond()

◆ solve() [1/2]

◆ solve() [2/2]

◆ solve_lx()

◆ update()

◆ upper_triangle()

◆ xt_m_inv_x()

◆ xt_m_x()

Member Data Documentation

◆ d_

◆ L_

◆ nullvector_

◆ num_dims_rank_def_

◆ rcond_