NAG Library Routine Document
g13ebf (multi_kalman_sqrt_invar)
1
Purpose
g13ebf performs a combined measurement and time update of one iteration of the timeinvariant Kalman filter using a square root covariance filter.
2
Specification
Fortran Interface
Subroutine g13ebf ( 
transf, n, m, l, a, lds, b, stq, q, ldq, c, ldm, r, s, k, h, u, tol, iwk, wk, ifail) 
Integer, Intent (In)  ::  n, m, l, lds, ldq, ldm  Integer, Intent (Inout)  ::  ifail  Integer, Intent (Out)  ::  iwk(m)  Real (Kind=nag_wp), Intent (In)  ::  q(ldq,*), r(ldc,m), tol  Real (Kind=nag_wp), Intent (Inout)  ::  a(lds,n), b(lds,l), c(ldc,n), s(lds,n), k(lds,m), h(ldc,m), u(lds,*)  Real (Kind=nag_wp), Intent (Out)  ::  wk((n+m)*(n+m+l))  Logical, Intent (In)  ::  stq  Character (1), Intent (In)  ::  transf 

C Header Interface
#include nagmk26.h
void 
g13ebf_ (const char *transf, const Integer *n, const Integer *m, const Integer *l, double a[], const Integer *lds, double b[], const logical *stq, const double q[], const Integer *ldq, double c[], const Integer *ldc, const double r[], double s[], double k[], double h[], double u[], const double *tol, Integer iwk[], double wk[], Integer *ifail, const Charlen length_transf) 

3
Description
The Kalman filter arises from the state space model given by
where
${X}_{i}$ is the state vector of length
$n$ at time
$i$,
${Y}_{i}$ is the observation vector of length
$m$ at time
$i$ and
${W}_{i}$ of length
$l$ and
${V}_{i}$ of length
$m$ are the independent state noise and measurement noise respectively. The matrices
$A,B$ and
$C$ are time invariant.
The estimate of
${X}_{i}$ given observations
${Y}_{1}$ to
${Y}_{i1}$ is denoted by
${\hat{X}}_{i\mid i1}$ with state covariance matrix
$\mathrm{Var}\left({\hat{X}}_{i\mid i1}\right)={P}_{i\mid i1}={S}_{i}{S}_{i}^{\mathrm{T}}$ while the estimate of
${X}_{i}$ given observations
${Y}_{1}$ to
${Y}_{i}$ is denoted by
${\hat{X}}_{i\mid i}$ with covariance matrix
$\mathrm{Var}\left({\hat{X}}_{i\mid i}\right)={P}_{i\mid i}$. The update of the estimate,
${\hat{X}}_{i\mid i1}$, from time
$i$ to time
$\left(i+1\right)$ is computed in two stages. First, the measurementupdate is given by
where
${K}_{i}={P}_{i\mid i}{C}^{\mathrm{T}}{\left[C{P}_{i\mid i}{C}^{\mathrm{T}}+{R}_{i}\right]}^{1}$ is the Kalman gain matrix. The second stage is the timeupdate for
$X$, which is given by
where
${D}_{i}{U}_{i}$ represents any deterministic control used.
The square root covariance filter algorithm provides a stable method for computing the Kalman gain matrix and the state covariance matrix. The algorithm can be summarised as
where
$U$ is an orthogonal transformation triangularizing the lefthand prearray to produce the righthand postarray. The triangularization is carried out via Householder transformations exploiting the zero pattern of the prearray. The relationship between the Kalman gain matrix
${K}_{i}$ and
${G}_{i}$ is given by
In order to exploit the invariant parts of the model to simplify the computation of
$U$ the results for the transformed state space
${U}^{*}X$ are computed where
${U}^{*}$ is the transformation that reduces the matrix pair
$\left(A,C\right)$ to lower observer Hessenberg form. That is, the matrix
${U}^{*}$ is computed such that the compound matrix
is a lower trapezoidal matrix. Further the matrix
$B$ is transformed to
${U}^{*}B$. These transformations need only be computed once at the start of a series, and
g13ebf will, optionally, compute them.
g13ebf returns transformed matrices
${U}^{*}A{U}^{*T}$,
${U}^{*}B$,
$C{U}^{*T}$ and
${U}^{*}A{K}_{i}$, the Cholesky factor of the updated transformed state covariance matrix
${S}_{i+1}^{*}$ (where
${U}^{*}{P}_{i+1\mid i}{U}^{*T}={S}_{i+1}^{*}{S}_{i+1}^{*T}$) and the matrix
${H}_{i}^{1/2}$, valid for both transformed and original models, which is used in the computation of the likelihood for the model. Note that the covariance matrices
${Q}_{i}$ and
${R}_{i}$ can be timevarying.
4
References
Vanbegin M, van Dooren P and Verhaegen M H G (1989) Algorithm 675: FORTRAN subroutines for computing the square root covariance filter and square root information filter in dense or Hessenberg forms ACM Trans. Math. Software 15 243–256
Verhaegen M H G and van Dooren P (1986) Numerical aspects of different Kalman filter implementations IEEE Trans. Auto. Contr. AC31 907–917
5
Arguments
 1: $\mathbf{transf}$ – Character(1)Input

On entry: indicates whether to transform the input matrix pair
$\left(A,C\right)$ to lower observer Hessenberg form. The transformation will only be required on the first call to
g13ebf.
 ${\mathbf{transf}}=\text{'T'}$
 The matrices in arrays a and c are transformed to lower observer Hessenberg form and the matrices in b and s are transformed as described in Section 3.
 ${\mathbf{transf}}=\text{'H'}$
 The matrices in arrays a, c and b should be as returned from a previous call to g13ebf with ${\mathbf{transf}}=\text{'T'}$.
Constraint:
${\mathbf{transf}}=\text{'T'}$ or $\text{'H'}$.
 2: $\mathbf{n}$ – IntegerInput

On entry: $n$, the size of the state vector.
Constraint:
${\mathbf{n}}\ge 1$.
 3: $\mathbf{m}$ – IntegerInput

On entry: $m$, the size of the observation vector.
Constraint:
${\mathbf{m}}\ge 1$.
 4: $\mathbf{l}$ – IntegerInput

On entry: $l$, the dimension of the state noise.
Constraint:
${\mathbf{l}}\ge 1$.
 5: $\mathbf{a}\left({\mathbf{lds}},{\mathbf{n}}\right)$ – Real (Kind=nag_wp) arrayInput/Output

On entry: if
${\mathbf{transf}}=\text{'T'}$, the state transition matrix,
$A$.
If ${\mathbf{transf}}=\text{'H'}$, the transformed matrix as returned by a previous call to g13ebf with ${\mathbf{transf}}=\text{'T'}$.
On exit: if
${\mathbf{transf}}=\text{'T'}$, the transformed matrix,
${U}^{*}A{U}^{*T}$, otherwise
a is unchanged.
 6: $\mathbf{lds}$ – IntegerInput

On entry: the first dimension of the arrays
a,
b,
s,
k and
u as declared in the (sub)program from which
g13ebf is called.
Constraint:
${\mathbf{lds}}\ge {\mathbf{n}}$.
 7: $\mathbf{b}\left({\mathbf{lds}},{\mathbf{l}}\right)$ – Real (Kind=nag_wp) arrayInput/Output

On entry: if
${\mathbf{transf}}=\text{'T'}$, the noise coefficient matrix
$B$.
If ${\mathbf{transf}}=\text{'H'}$, the transformed matrix as returned by a previous call to g13ebf with ${\mathbf{transf}}=\text{'T'}$.
On exit: if
${\mathbf{transf}}=\text{'T'}$, the transformed matrix,
${U}^{*}B$, otherwise
b is unchanged.
 8: $\mathbf{stq}$ – LogicalInput

On entry: if
${\mathbf{stq}}=\mathrm{.TRUE.}$, the state noise covariance matrix
${Q}_{i}$ is assumed to be the identity matrix. Otherwise the lower triangular Cholesky factor,
${Q}_{i}^{1/2}$, must be provided in
q.
 9: $\mathbf{q}\left({\mathbf{ldq}},*\right)$ – Real (Kind=nag_wp) arrayInput

Note: the second dimension of the array
q
must be at least
${\mathbf{l}}$ if
${\mathbf{stq}}=\mathrm{.FALSE.}$.
On entry: if
${\mathbf{stq}}=\mathrm{.FALSE.}$,
q must contain the lower triangular Cholesky factor of the state noise covariance matrix,
${Q}_{i}^{1/2}$. Otherwise
q is not referenced.
 10: $\mathbf{ldq}$ – IntegerInput

On entry: the first dimension of the array
q as declared in the (sub)program from which
g13ebf is called.
Constraints:
 if ${\mathbf{stq}}=\mathrm{.FALSE.}$, ${\mathbf{ldq}}\ge {\mathbf{l}}$;
 otherwise ${\mathbf{ldq}}\ge 1$.
 11: $\mathbf{c}\left({\mathbf{ldm}},{\mathbf{n}}\right)$ – Real (Kind=nag_wp) arrayInput/Output

On entry: if
${\mathbf{transf}}=\text{'T'}$, the measurement coefficient matrix,
$C$.
If ${\mathbf{transf}}=\text{'H'}$, the transformed matrix as returned by a previous call to g13ebf with ${\mathbf{transf}}=\text{'T'}$.
On exit: if
${\mathbf{transf}}=\text{'T'}$, the transformed matrix,
$C{U}^{*T}$, otherwise
c is unchanged.
 12: $\mathbf{ldm}$ – IntegerInput

On entry: the first dimension of the arrays
c,
r and
h as declared in the (sub)program from which
g13ebf is called.
Constraint:
${\mathbf{ldm}}\ge {\mathbf{m}}$.
 13: $\mathbf{r}\left({\mathbf{ldm}},{\mathbf{m}}\right)$ – Real (Kind=nag_wp) arrayInput

On entry: the lower triangular Cholesky factor of the measurement noise covariance matrix ${R}_{i}^{1/2}$.
 14: $\mathbf{s}\left({\mathbf{lds}},{\mathbf{n}}\right)$ – Real (Kind=nag_wp) arrayInput/Output

On entry: if
${\mathbf{transf}}=\text{'T'}$ the lower triangular Cholesky factor of the state covariance matrix,
${S}_{i}$.
If ${\mathbf{transf}}=\text{'H'}$ the lower triangular Cholesky factor of the covariance matrix of the transformed state vector ${S}_{i}^{*}$ as returned from a previous call to g13ebf with ${\mathbf{transf}}=\text{'T'}$.
On exit: the lower triangular Cholesky factor of the transformed state covariance matrix, ${S}_{i+1}^{*}$.
 15: $\mathbf{k}\left({\mathbf{lds}},{\mathbf{m}}\right)$ – Real (Kind=nag_wp) arrayOutput

On exit: the Kalman gain matrix for the transformed state vector premultiplied by the state transformed transition matrix, ${U}^{*}A{K}_{i}$.
 16: $\mathbf{h}\left({\mathbf{ldm}},{\mathbf{m}}\right)$ – Real (Kind=nag_wp) arrayOutput

On exit: the lower triangular matrix ${H}_{i}^{1/2}$.
 17: $\mathbf{u}\left({\mathbf{lds}},*\right)$ – Real (Kind=nag_wp) arrayOutput

Note: the second dimension of the array
u
must be at least
${\mathbf{n}}$ if
${\mathbf{transf}}=\text{'T'}$.
On exit: if
${\mathbf{transf}}=\text{'T'}$ the
$n$ by
$n$ transformation matrix
${U}^{*}$, otherwise
u is not referenced.
 18: $\mathbf{tol}$ – Real (Kind=nag_wp)Input

On entry: the tolerance used to test for the singularity of
${H}_{i}^{1/2}$. If
$0.0\le {\mathbf{tol}}<{m}^{2}\times \mathit{machineprecision}$, then
${m}^{2}\times \mathit{machineprecision}$ is used instead. The inverse of the condition number of
${H}^{1/2}$ is estimated by a call to
f07tgf (dtrcon). If this estimate is less than
tol then
${H}^{1/2}$ is assumed to be singular.
Suggested value:
${\mathbf{tol}}=0.0$.
Constraint:
${\mathbf{tol}}\ge 0.0$.
 19: $\mathbf{iwk}\left({\mathbf{m}}\right)$ – Integer arrayWorkspace
 20: $\mathbf{wk}\left(\left({\mathbf{n}}+{\mathbf{m}}\right)\times \left({\mathbf{n}}+{\mathbf{m}}+{\mathbf{l}}\right)\right)$ – Real (Kind=nag_wp) arrayWorkspace

 21: $\mathbf{ifail}$ – IntegerInput/Output

On entry:
ifail must be set to
$0$,
$1\text{ or}1$. If you are unfamiliar with this argument you should refer to
Section 3.4 in How to Use the NAG Library and its Documentation for details.
For environments where it might be inappropriate to halt program execution when an error is detected, the value
$1\text{ or}1$ is recommended. If the output of error messages is undesirable, then the value
$1$ is recommended. Otherwise, if you are not familiar with this argument, the recommended value is
$0$.
When the value $\mathbf{1}\text{ or}\mathbf{1}$ is used it is essential to test the value of ifail on exit.
On exit:
${\mathbf{ifail}}={\mathbf{0}}$ unless the routine detects an error or a warning has been flagged (see
Section 6).
6
Error Indicators and Warnings
If on entry
${\mathbf{ifail}}=0$ or
$1$, explanatory error messages are output on the current error message unit (as defined by
x04aaf).
Errors or warnings detected by the routine:
 ${\mathbf{ifail}}=99$
An unexpected error has been triggered by this routine. Please
contact
NAG.
See
Section 3.9 in How to Use the NAG Library and its Documentation for further information.
 ${\mathbf{ifail}}=399$
Your licence key may have expired or may not have been installed correctly.
See
Section 3.8 in How to Use the NAG Library and its Documentation for further information.
 ${\mathbf{ifail}}=999$
Dynamic memory allocation failed.
See
Section 3.7 in How to Use the NAG Library and its Documentation for further information.
7
Accuracy
The use of the square root algorithm improves the stability of the computations as compared with the direct coding of the Kalman filter. The accuracy will depend on the model.
8
Parallelism and Performance
g13ebf is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.
g13ebf makes calls to BLAS and/or LAPACK routines, which may be threaded within the vendor library used by this implementation. Consult the documentation for the vendor library for further information.
Please consult the
X06 Chapter Introduction for information on how to control and interrogate the OpenMP environment used within this routine. Please also consult the
Users' Note for your implementation for any additional implementationspecific information.
For models with timevarying
$A,B$ and
$C$,
g13eaf can be used.
The initial estimate of the transformed state vector can be computed from the estimate of the original state vector
${\hat{X}}_{1\mid 0}$, say, by premultiplying it by
${U}^{*}$ as returned by
g13ebf with
${\mathbf{transf}}=\text{'T'}$; that is,
${\hat{X}}_{1\mid 0}^{*}={U}^{*}{\hat{X}}_{1\mid 0}$. The estimate of the transformed state vector
${\hat{X}}_{i+1\mid i}^{*}$ can be computed from the previous value
${\hat{X}}_{i\mid i1}^{*}$ by
where
are the independent onestep prediction residuals for both the transformed and original model. The estimate of the original state vector can be computed from the transformed state vector as
${U}^{*T}{\hat{X}}_{1+1\mid i}^{*}$. The required matrixvector multiplications can be performed by
f06paf (dgemv).
If
${W}_{i}$ and
${V}_{i}$ are independent multivariate Normal variates then the loglikelihood for observations
$i=1,2,\dots ,t$ is given by
where
$\kappa $ is a constant.
The Cholesky factors of the covariance matrices can be computed using
f07fdf (dpotrf).
Note that the model
can be specified either with
b set to the identity matrix and
${\mathbf{stq}}=\mathrm{.FALSE.}$ and the matrix
${Q}^{1/2}$ input in
q or with
${\mathbf{stq}}=\mathrm{.TRUE.}$ and
b set to
${Q}^{1/2}$.
The algorithm requires
$\frac{1}{6}{n}^{3}+{n}^{2}\left(\frac{3}{2}m+l\right)+2n{m}^{2}+\frac{2}{3}{p}^{3}$ operations and is backward stable (see
Verhaegen and van Dooren (1986)). The transformation to lower observer Hessenberg form requires
$\mathit{O}\left(\left(n+m\right){n}^{2}\right)$ operations.
10
Example
This example first inputs the number of updates to be computed and the problem sizes. The initial state vector and the Cholesky factor of the state covariance matrix are input followed by the model matrices
$A,B,C,{R}^{1/2}$ and optionally
${Q}^{1/2}$ (the Cholesky factors of the covariance matrices being input). At the first update the matrices are transformed using the
${\mathbf{transf}}=\text{'T'}$ option and the initial value of the state vector is transformed. At each update the observed values are input and the residuals are computed and printed and the estimate of the transformed state vector,
${\hat{U}}^{*}{X}_{i\mid i1}$, and the deviance are updated. The deviance is
$2\times \text{loglikelihood}$ ignoring the constant. After the final update the estimate of the state vector is computed from the transformed state vector and the state covariance matrix is computed from
s and these are printed along with the value of the deviance.
The data is for a twodimensional time series to which a VARMA
$\left(1,1\right)$ has been fitted. For the specification of a VARMA model as a state space model see the
G13 Chapter Introduction. The means of the two series are included as additional states that do not change over time. The initial value of
$P$,
${P}_{0}$, is the solution to
10.1
Program Text
Program Text (g13ebfe.f90)
10.2
Program Data
Program Data (g13ebfe.d)
10.3
Program Results
Program Results (g13ebfe.r)