Detailed Description
class SRI encapsulates all the information associated with the solution of a set of simultaneous linear equations. It is used in least squares estimation (linear and linearized) and is the basis of the preferred implementation of Kalman filtering. An SRI consists of just three things: (1) 'R', the 'information matrix', which is an upper triangular matrix of dimension N, equal to the inverse of the square root (or Cholesky decomposition) of the solution covariance matrix, (2) 'Z', the 'SRI state vector' of length N (parallels the components of R), (not to be confused with the regular state vector X), and (3) 'names', a Namelist used to label the elements of R and Z (parallels and labels rows and columns of R and elements of Z). A Namelist is part of class SRI because the manipulations of SRI (see functions below) requires a consistent way of manipulating the different individual elements of R and Z, in addition it allows the user to attach 'human-readable' labels to the elements of the state vector, which is useful in adding, dropping and bumping states, and it makes printed results more readable (see the LabeledMatrix class in Namelist.hpp).
The set of simultaneous equations represented by an SRI is R * X = Z, where X is the (unknown) state vector (the conventional solution vector) also of dimension N. The state X is solved for as X = inverse(R) * Z, and the covariance matrix of the state X is equal to transpose(inverse(R))*inverse(R).
Least squares estimation via SRI is very simple and efficient; it uses the Householder transformation to convert the problem to upper triangular form, and then uses very efficient algorithms to invert the information matrix to find the solution and its covariance. The usual matrix equation is H * X = D, where H is the 'design matrix' or the 'partials matrix', of dimension M x N, X is the (unknown) solution vector of length N, and D is the 'data' or 'measurement' vector of length M. In the least squares 'update' of the SRI, this set of information {H,D} is concatenated with the existing SRI {R,Z} to form an (N+M x N+1) matrix Q which has R in the upper left, Z upper right, H lower left and D lower right. This extended matrix is then subjected to a Householder transformation (see class Matrix), which will put (at least the first N columns of) Q into upper triangular form. The result is a new, updated SRI (R and Z) in the place of the old, while in place of D are residuals of fit corresponding to the measurements in D (the H part of Q is trashed). This result, in fact (see the reference), produces an updated SRI which gives precisely the usual least squares solution for the combined 'a priori SRI + new data' problem. This algorithm is called a 'measurement update' of the SRI.
It is most enlightening to think of the SRI and this process in terms of 'information'. The SRI contains all the 'information' which has come from updates that have been made to it using (H,D) pairs. Initially, the SRI is all zeros, which corresponds to 'no information'. This overcomes one serious problem with conventional least squares and the Kalman algorithm, namely that a 'zero information' starting value cannot be correctly expressed, because in that case the covariance matrix is singular and the state vector is indeterminate; in the SRI method this is perfectly consistent - the covariance matrix is singular because the information matrix (R) is zero, and thus the state is entirely indeterminate. As new 'information' (in the form of data D and partials matrix H pairs) is added to the SRI (via the Householder algorithm), the 'information' stored in R and Z is increased and they become non-zero. (By the way note that the number of rows in the {H,D} information is arbitrary - information can be added in 'batches' - M large - or one - M=1 - piece at a time.) When there is enough information, R becomes non-singular, and so can be inverted and the solution and covariance can be computed. As the amount of information becomes large, elements of R become large, and thus elements of the covariance (think of covariance as a measure of uncertainty - the opposite or inverse of information) become small.
The structure of the SRI method allows some powerful techniques to be used in manipulating, combining and separating state elements and the information associated with them in SRIs. For example, if the measurement updates have failed to increase the information about one particular state element, then that element, and its information, may be removed from the problem by deleting that element's row and column of R, and its element of Z (and then re-triangularizing the SRI). In general, any subset of an SRI may be separated, or the SRI split (see the routine of that name below - note the caveats) into two separate SRIs. For another example, SRI allows the information of a each state element to be selectively reduced or even zeroed, simply by multiplying the corresponding elements of R and Z by a factor; in Kalman filtering this is called a 'Q bump' of the element and is very important in some filtering applications. There are methods (see below) consistently to merge (operator+()), split, and permute elements of, SRIs.
Kalman filtering is an important application of SRI methods (actually it is called 'square root information filtering' or SRIF - technically the term 'Kalman filter algorithm' is reserved for the classical algorithm just as Kalman presented it, in terms of a state vector and its covariance matrix). The measurment update described above (which is actually just linear least squares) is half of the SRIF (Kalman filter) - there is a 'time update' that propagates the SRI (and thus the state and covariance) forward in time using the dynamical model of the filter. These are algebraically equivalent to the classical Kalman algorithm, but are more efficient and numerically stable (actually the Kalman algorithm has been shown to be numerically unstable!). There are even SRI smoothing algorithms, corresponding to Kalman smoothers, which consist of a 'backwards' filter, implemented by applying a 'smoother update' to the SRI at each point in reverse order.
Ref: Bierman, G.J. "Factorization Methods for Discrete Sequential Estimation," Academic Press, 1977.
#include <SRI.hpp>
Public Member Functions | |
| void | addAPriori (const Matrix< double > &Cov, const Vector< double > &X) |
| void | addAPrioriInformation (const Matrix< double > &ICov, const Vector< double > &X) |
| SRI & | append (const SRI &S) |
| void | getConditionNumber (double &small, double &big) const |
| std::string | getName (const unsigned int in) const |
| Namelist | getNames () const |
| Matrix< double > | getR () const |
| void | getState (Vector< double > &X, int *ptrSingularIndex=NULL) const |
| void | getStateAndCovariance (Vector< double > &X, Matrix< double > &C, double *ptrSmall=NULL, double *ptrBig=NULL) const |
| Vector< double > | getZ () const |
| unsigned int | index (std::string &name) |
| void | measurementUpdate (Matrix< double > &Partials, Vector< double > &Data) |
| void | measurementUpdate (SparseMatrix< double > &Partials, Vector< double > &Data) |
| void | merge (const SRI &S) |
| SRI & | operator+= (const Namelist &NL) |
| SRI & | operator+= (const SRI &S) |
| SRI & | operator= (const SRI &right) |
| operator= More... | |
| void | permute (const Namelist &NL) |
| void | Qbump (const unsigned int in, double q=0.0) |
| void | reshape (const Namelist &NL) |
| void | retriangularize (const Matrix< double > &A) |
| void | retriangularize (Matrix< double > RR, Vector< double > ZZ) |
| void | setFromCovState (const Matrix< double > &Cov, const Vector< double > &State, const Namelist &NL) |
| bool | setName (const unsigned int in, const std::string &label) |
| void | shift (const Vector< double > &X0) |
| void | shiftZ (const Vector< double > &Z0) |
| unsigned int | size () const |
| void | split (const Namelist &NL, SRI &S) |
| SRI () | |
| empty constructor More... | |
| SRI (const Matrix< double > &R, const Vector< double > &Z, const Namelist &NL) | |
| SRI (const Namelist &NL) | |
| SRI (const SRI &s) | |
| copy constructor More... | |
| SRI (const unsigned int N) | |
| void | stateFix (const std::string &name, double value, double sigma, bool restore) |
| void | stateFix (const unsigned int index, double value, double sigma, bool restore) |
| void | stateFixAndRemove (const Namelist &drops, const Vector< double > &values) |
| void | stateFixAndRemove (const unsigned int index, double value) |
| void | transform (const Matrix< double > &invT, const Namelist &NL) |
| void | zeroAll (const unsigned int n=0) |
| void | zeroOne (const unsigned int n) |
| void | zeroState () |
| Zero out (set all elements to zero) the state (Vector Z) only. More... | |
Protected Attributes | |
| Namelist | names |
| Namelist parallel to R and Z, labelling the elements of the state vector. More... | |
| Matrix< double > | R |
| Information matrix, an upper triangular (square) matrix. More... | |
| Vector< double > | Z |
| SRI state vector, of length equal to the dimension (row and col) of R. More... | |
Friends | |
| SRI | operator+ (const SRI &Sleft, const SRI &Sright) |
| std::ostream & | operator<< (std::ostream &s, const SRI &) |
| output operator More... | |
Constructor & Destructor Documentation
◆ SRI() [1/5]
◆ SRI() [2/5]
| gnsstk::SRI::SRI | ( | const unsigned int | N | ) |
◆ SRI() [3/5]
◆ SRI() [4/5]
explicit constructor - throw if the dimensions are inconsistent. User is responsible for ensuring the input is self-consistent.
- Exceptions
-
MatrixException if dimensions are not consistent
◆ SRI() [5/5]
Member Function Documentation
◆ addAPriori()
Add a priori or constraint information in the form of an ordinary state vector and covariance matrix. The matrix must be non-singular.
- Parameters
-
Cov Covariance matrix of same dimension as this SRIFilter X State vector of same dimension as this SRIFilter
- Exceptions
-
MatrixException if input is invalid: dimensions are wrong or Cov is singular.
◆ addAPrioriInformation()
| void gnsstk::SRI::addAPrioriInformation | ( | const Matrix< double > & | ICov, |
| const Vector< double > & | X | ||
| ) |
Add a priori or constraint information in the form of an information matrix (inverse covariance) and ordinary state. ICov must be non-singular.
- Parameters
-
ICov Inverse covariance matrix of same dimension as this SRIFilter X State vector of same dimension as this SRIFilter
- Exceptions
-
MatrixException if input is invalid: dimensions are wrong.
◆ append()
append an SRI onto this SRI. Similar to opertor+= but simpler; input SRI is simply appended, first using operator+=(Namelist), then filling the new portions of R and Z, all without final Householder transformation of result. Do not allow a name that is already present to be added: throw.
- Parameters
-
S input SRI to be appended
- Returns
- appended SRI
- Exceptions
-
SomeUnknownException if a name is repeated MatrixException VectorException
◆ getConditionNumber()
| void gnsstk::SRI::getConditionNumber | ( | double & | small, |
| double & | big | ||
| ) | const |
Compute the condition number, or rather the largest and smallest eigenvalues of the SRI matrix R (the condition number is the ratio of the largest and smallest eigenvalues). Note that the condition number of the covariance matrix would be the square of the condition number of R.
- Parameters
-
small smallest eigenvalue big largest eigenvalue, condition = big/small
- Exceptions
-
MatrixException
◆ getName()
|
inline |
◆ getNames()
|
inline |
◆ getR()
|
inline |
◆ getState()
Compute the state X without computing the covariance matrix C. R*X=Z so X=inverse(R)*Z; in this routine the state is computed explicitly, without forming the inverse of R, using the fact that R is upper triangular. NB. The matrix is singular if and only if one or more of the diagonal elements is zero; in the case the routine will still have valid entries in the state vector for index greater than the largest index with zero diagonal
- Parameters
-
X State vector (output) ptrSingularIndex if ptr is non-null, on output *ptr will be the largest index of singularity
- Exceptions
-
MatrixException if R is singular.
◆ getStateAndCovariance()
| void gnsstk::SRI::getStateAndCovariance | ( | Vector< double > & | X, |
| Matrix< double > & | C, | ||
| double * | ptrSmall = NULL, |
||
| double * | ptrBig = NULL |
||
| ) | const |
Compute the state X and the covariance matrix C of the state, where C = transpose(inverse(R))*inverse(R) and X = inverse(R) * Z. Optional pointer arguments will return smallest and largest eigenvalues of the R matrix, which is a measure of singularity. NB this is the most efficient way to invert the SRI equation.
- Parameters
-
X State vector (output) C Covariance of the state vector (output) ptrSmall Pointer to double, on output *ptrSmall set to smallest eigenvalue of R ptrBig Pointer to double, on output *ptrBig set to largest eigenvalue of R
- Exceptions
-
MatrixException if R is singular. VectorException
◆ getZ()
|
inline |
◆ index()
|
inline |
◆ measurementUpdate() [1/2]
◆ measurementUpdate() [2/2]
|
inline |
◆ merge()
merge an SRI into this one. NB names may be reordered in the result. NB this is just operator+=()
- Parameters
-
S SRI to be merged into this
- Exceptions
-
MatrixException VectorException
◆ operator+=() [1/2]
◆ operator+=() [2/2]
◆ operator=()
◆ permute()
◆ Qbump()
| void gnsstk::SRI::Qbump | ( | const unsigned int | in, |
| double | q = 0.0 |
||
| ) |
Decrease the information in this SRI for, or 'Q bump', the element with the input index. This means that the uncertainty and the state element given by the index are divided by the input factor q; the default input is zero, which means zero out the information (q = infinite). A Q bump by factor q is equivalent to 'de-weighting' the element by q. No effect if in is out of range.
- Parameters
-
in index to bump q factor by which to 'de-weight' the element(in), default 0.0 which implies all information removed.
- Exceptions
-
MatrixException VectorException
◆ reshape()
◆ retriangularize() [1/2]
Retriangularize the SRI, when it has been modified to a non-UT matrix (e.g. by transform()). Given the matrix A=[R||Z], apply HH transforms to retriangularize it and pull out new R and Z. NB caller must modify names, if necessary
- Parameters
-
A Matrix<double> which is [R || Z] to be retriangularizied.
- Exceptions
-
MatrixException if dimensions are wrong.
◆ retriangularize() [2/2]
Retriangularize the SRI, that is assuming R has been modified to a non-UT matrix (e.g. by transform()). Given RR and ZZ, apply HH transforms to retriangularize, and store as R,Z. NB caller must modify names, if necessary
- Parameters
-
[in] RR the modified (non-UT) R [in] ZZ the (potentially) modified Z
- Exceptions
-
MatrixException if dimensions are wrong.
◆ setFromCovState()
◆ setName()
◆ shift()
◆ shiftZ()
◆ size()
|
inline |
◆ split()
split an SRI into two others, this one matching the input Namelist, the other containing whatever is left. The input Namelist must be a non-trivial subset of this->names; throw MatrixException if it is not. NB. Interpreting the results of a split() and merge (operator+()) operations should be done very carefully; remember that the SRI contains both solution and noise, and that the results of these operations are not always as expected, particularly note that split() and operator+() are usually NOT reversible. NB output SRI S will be singular.
- Parameters
-
NL new namelist to be given to this object S this SRI after the NL part has been removed
- Exceptions
-
MatrixException VectorException
◆ stateFix() [1/2]
| void gnsstk::SRI::stateFix | ( | const std::string & | name, |
| double | value, | ||
| double | sigma, | ||
| bool | restore | ||
| ) |
Fix one state element (with the given name) to a given value, and set the information for that element (== 1/sigma) to a given value. No effect if name is not found
- Parameters
-
name string labeling the state in Namelist names value to which the state element is fixed sigma (1/information) assigned to the element restore if true, permute back to the original form after fixing
- Exceptions
-
Exception
◆ stateFix() [2/2]
| void gnsstk::SRI::stateFix | ( | const unsigned int | index, |
| double | value, | ||
| double | sigma, | ||
| bool | restore | ||
| ) |
Fix one state element (at the given index) to a given value, and set the information for that element (== 1/sigma) to a given value. No effect if index is out of range.
- Parameters
-
index of the element to fix value to which the state element is fixed sigma (1/information) assigned to the element restore if true, permute back to the original form after fixing
- Exceptions
-
Exception
◆ stateFixAndRemove() [1/2]
Vector version of stateFixAndRemove, with Namelist identifying the states. Fix the given state elements to the input value, and collapse the SRI by removing those elements. No effect if name is not found.
- Parameters
-
drops Namelist of states to fix values vector parallel to drops of values to which elements are fixed
- Exceptions
-
MatrixException VectorException
◆ stateFixAndRemove() [2/2]
| void gnsstk::SRI::stateFixAndRemove | ( | const unsigned int | index, |
| double | value | ||
| ) |
Fix the state element with the input index to the input value, and collapse the SRI by removing that element. No effect if index is out of range.
- Parameters
-
index of the element to fix value to which the state element is fixed
- Exceptions
-
MatrixException VectorException
◆ transform()
Transform the state by the transformation matrix T; i.e. X -> T*X; this is done by right multiplying R by inverse(T), which is the input. Thus R -> R*inverse(T), so Z -> R*inverse(T)*T*X = Z. [R|Z] -> [R*invT|Z]. NB Input is the inverse of the transformation.
- Parameters
-
invT Matrix<double> inverse of the transformation T : X->T*X NL Namelist of the transformed SRI, SRI.names is set to this
- Exceptions
-
MatrixException if input dimensions are wrong.
◆ zeroAll()
| void gnsstk::SRI::zeroAll | ( | const unsigned int | n = 0 | ) |
◆ zeroOne()
| void gnsstk::SRI::zeroOne | ( | const unsigned int | n | ) |
◆ zeroState()
|
inline |
Friends And Related Function Documentation
◆ operator+
merge two SRIs to produce a third. ? should this be operator&() ?
- Exceptions
-
MatrixException VectorException
◆ operator<<
Member Data Documentation
◆ names
|
protected |
◆ R
|
protected |
◆ Z
|
protected |
The documentation for this class was generated from the following files: