hide long namesshow long names
hide short namesshow short names
Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

NAG Toolbox: nag_tsa_multi_kalman_sqrt_invar (g13eb)

Purpose

nag_tsa_multi_kalman_sqrt_invar (g13eb) performs a combined measurement and time update of one iteration of the time-invariant Kalman filter using a square root covariance filter.

Syntax

[a, b, c, s, k, h, u, ifail] = g13eb(transf, a, b, stq, q, c, r, s, 'n', n, 'm', m, 'l', l, 'tol', tol)
[a, b, c, s, k, h, u, ifail] = nag_tsa_multi_kalman_sqrt_invar(transf, a, b, stq, q, c, r, s, 'n', n, 'm', m, 'l', l, 'tol', tol)

Description

The Kalman filter arises from the state space model given by
Xi + 1 = AXi + BWi, Var(Wi) = Qi
Yi = CXi + Vi, Var(Vi) = Ri
Xi+1=AXi+BWi, Var(Wi)=Qi Yi=CXi+Vi, Var(Vi)=Ri
where XiXi is the state vector of length nn at time ii, YiYi is the observation vector of length mm at time ii and WiWi of length ll and ViVi of length mm are the independent state noise and measurement noise respectively. The matrices A,BA,B and CC are time invariant.
The estimate of XiXi given observations Y1Y1 to Yi1Yi-1 is denoted by ii1X^ii-1 with state covariance matrix Var(ii1) = Pii1 = SiSiTVar(X^ii-1)=Pii-1=SiSiT while the estimate of XiXi given observations Y1Y1 to YiYi is denoted by iiX^ii with covariance matrix Var(ii) = PiiVar(X^ii)=Pii. The update of the estimate, ii1X^ii-1, from time ii to time (i + 1)(i+1) is computed in two stages. First, the measurement-update is given by
ii = ii1 + Ki[YiCii1]
X^ii=X^ii-1+Ki[Yi-CX^ii-1]
(1)
where Ki = PiiCT[CPiiCT + Ri]1Ki=PiiCT[CPiiCT+Ri] -1 is the Kalman gain matrix. The second stage is the time-update for XX, which is given by
i + 1i = Aii + DiUi
X^i+1i=AX^ii+DiUi
(2)
where DiUiDiUi 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
  Ri1 / 2 0 CSi 0 BQi1 / 2 ASi  
U =
  Hi1 / 2 0 0 Gi Si + 1 0  
Ri1/2 0 CSi 0 BQi1/2 ASi U= Hi1/2 0 0 Gi Si+1 0
where UU is an orthogonal transformation triangularizing the left-hand pre-array to produce the right-hand post-array. The triangularization is carried out via Householder transformations exploiting the zero pattern of the pre-array. The relationship between the Kalman gain matrix KiKi and GiGi is given by
AKi = Gi (Hi1 / 2)1.
AKi=Gi (Hi1/2) -1.
In order to exploit the invariant parts of the model to simplify the computation of UU the results for the transformed state space U*XU*X are computed where U*U* is the transformation that reduces the matrix pair (A,C)(A,C) to lower observer Hessenberg form. That is, the matrix U*U* is computed such that the compound matrix
[ CU * T U*AU * T ]
[ CU*T U*AU*T ]
is a lower trapezoidal matrix. Further the matrix BB is transformed to U*BU*B. These transformations need only be computed once at the start of a series, and nag_tsa_multi_kalman_sqrt_invar (g13eb) will, optionally, compute them. nag_tsa_multi_kalman_sqrt_invar (g13eb) returns transformed matrices U*AU * TU*AU*T, U*BU*B, CU * TCU*T and U*AKiU*AKi, the Cholesky factor of the updated transformed state covariance matrix Si + 1 * Si+1* (where U*Pi + 1iU * T = Si + 1 * Si + 1 * TU*Pi+1iU*T=Si+1*Si+1 *T) and the matrix Hi1 / 2Hi1/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 QiQi and RiRi can be time-varying.

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. AC-31 907–917

Parameters

Compulsory Input Parameters

1:     transf – string (length ≥ 1)
Indicates whether to transform the input matrix pair (A,C)(A,C) to lower observer Hessenberg form. The transformation will only be required on the first call to nag_tsa_multi_kalman_sqrt_invar (g13eb).
transf = 'T'transf='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 [Description].
transf = 'H'transf='H'
The matrices in arrays a, c and b should be as returned from a previous call to nag_tsa_multi_kalman_sqrt_invar (g13eb) with transf = 'T'transf='T'.
Constraint: transf = 'T'transf='T' or 'H''H'.
2:     a(lds,n) – double array
lds, the first dimension of the array, must satisfy the constraint ldsnldsn.
If transf = 'T'transf='T', the state transition matrix, AA.
If transf = 'H'transf='H', the transformed matrix as returned by a previous call to nag_tsa_multi_kalman_sqrt_invar (g13eb) with transf = 'T'transf='T'.
3:     b(lds,l) – double array
lds, the first dimension of the array, must satisfy the constraint ldsnldsn.
If transf = 'T'transf='T', the noise coefficient matrix BB.
If transf = 'H'transf='H', the transformed matrix as returned by a previous call to nag_tsa_multi_kalman_sqrt_invar (g13eb) with transf = 'T'transf='T'.
4:     stq – logical scalar
If stq = truestq=true, the state noise covariance matrix QiQi is assumed to be the identity matrix. Otherwise the lower triangular Cholesky factor, Qi1 / 2Qi1/2, must be provided in q.
5:     q(ldq, : :) – double array
The first dimension, ldq, of the array q must satisfy
  • if stq = falsestq=false, ldqlldql;
  • otherwise ldq1ldq1.
The second dimension of the array must be at least ll if stq = falsestq=false and at least 11 if stq = truestq=true
If stq = falsestq=false, q must contain the lower triangular Cholesky factor of the state noise covariance matrix, Qi1 / 2Qi1/2. Otherwise q is not referenced.
6:     c(ldm,n) – double array
ldm, the first dimension of the array, must satisfy the constraint ldmmldmm.
If transf = 'T'transf='T', the measurement coefficient matrix, CC.
If transf = 'H'transf='H', the transformed matrix as returned by a previous call to nag_tsa_multi_kalman_sqrt_invar (g13eb) with transf = 'T'transf='T'.
7:     r(ldm,m) – double array
ldm, the first dimension of the array, must satisfy the constraint ldmmldmm.
The lower triangular Cholesky factor of the measurement noise covariance matrix Ri1 / 2Ri1/2.
8:     s(lds,n) – double array
lds, the first dimension of the array, must satisfy the constraint ldsnldsn.
If transf = 'T'transf='T' the lower triangular Cholesky factor of the state covariance matrix, SiSi.
If transf = 'H'transf='H' the lower triangular Cholesky factor of the covariance matrix of the transformed state vector Si * Si* as returned from a previous call to nag_tsa_multi_kalman_sqrt_invar (g13eb) with transf = 'T'transf='T'.

Optional Input Parameters

1:     n – int64int32nag_int scalar
Default: The first dimension of the arrays a, b, s and the second dimension of the arrays a, c, s. (An error is raised if these dimensions are not equal.)
nn, the size of the state vector.
Constraint: n1n1.
2:     m – int64int32nag_int scalar
Default: The first dimension of the arrays c, r and the second dimension of the array r. (An error is raised if these dimensions are not equal.)
mm, the size of the observation vector.
Constraint: m1m1.
3:     l – int64int32nag_int scalar
Default: The second dimension of the array b.
ll, the dimension of the state noise.
Constraint: l1l1.
4:     tol – double scalar
The tolerance used to test for the singularity of Hi1 / 2Hi1/2. If 0.0tol < m2 × machine precision0.0tol<m2×machine precision, then m2 × machine precisionm2×machine precision is used instead. The inverse of the condition number of H1 / 2H1/2 is estimated by a call to nag_lapack_dtrcon (f07tg). If this estimate is less than tol then H1 / 2H1/2 is assumed to be singular.
Default: 0.00.0
Constraint: tol0.0tol0.0.

Input Parameters Omitted from the MATLAB Interface

lds ldq ldm iwk wk

Output Parameters

1:     a(lds,n) – double array
ldsnldsn.
If transf = 'T'transf='T', the transformed matrix, U*AU * TU*AU*T, otherwise a is unchanged.
2:     b(lds,l) – double array
ldsnldsn.
If transf = 'T'transf='T', the transformed matrix, U*BU*B, otherwise b is unchanged.
3:     c(ldm,n) – double array
ldmmldmm.
If transf = 'T'transf='T', the transformed matrix, CU * TCU*T, otherwise c is unchanged.
4:     s(lds,n) – double array
ldsnldsn.
The lower triangular Cholesky factor of the transformed state covariance matrix, Si + 1 * Si+1*.
5:     k(lds,m) – double array
ldsnldsn.
The Kalman gain matrix for the transformed state vector premultiplied by the state transformed transition matrix, U*AKiU*AKi.
6:     h(ldm,m) – double array
ldmmldmm.
The lower triangular matrix Hi1 / 2Hi1/2.
7:     u(lds, : :) – double array
The first dimension of the array u will be nn
The second dimension of the array will be nn if transf = 'T'transf='T', and at least 11 otherwise
ldsnldsn.
If transf = 'T'transf='T' the nn by nn transformation matrix U*U*, otherwise u is not referenced.
8:     ifail – int64int32nag_int scalar
ifail = 0ifail=0 unless the function detects an error (see [Error Indicators and Warnings]).

Error Indicators and Warnings

Errors or warnings detected by the function:
  ifail = 1ifail=1
On entry,transf'T'transf'T' or 'H''H',
orn < 1n<1,
orm < 1m<1,
orl < 1l<1,
orlds < nlds<n,
orldm < mldm<m,
orstq = truestq=true and ldq < 1ldq<1,
orstq = falsestq=false and ldq < lldq<l,
ortol < 0.0tol<0.0.
  ifail = 2ifail=2
The matrix Hi1 / 2Hi1/2 is singular.

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.

Further Comments

For models with time-varying A,BA,B and CC, nag_tsa_multi_kalman_sqrt_var (g13ea) can be used.
If WiWi and ViVi are independent multivariate Normal variates then the log-likelihood for observations i = 1,2,,ti=1,2,,t is given by
t t
l(θ) = κ(1/2) l n (det(Hi)) (1/2)(YiCiXii1)THi1(YiCiXii1)
i = 1 i = 1
l(θ) = κ - 12 i=1 t l n ( det(Hi) ) - 12 i=1 t ( Yi - Ci Xii-1 )T Hi-1 ( Yi - Ci X ii-1 )
where κκ is a constant.
The Cholesky factors of the covariance matrices can be computed using nag_lapack_dpotrf (f07fd).
Note that the model
Xi + 1 = AXi + Wi, Var(Wi) = Qi
Yi = CXi + Vi, Var(Vi) = Ri
Xi+1=AXi+Wi, Var(Wi)=Qi Yi=CXi+Vi, Var(Vi)=Ri
can be specified either with b set to the identity matrix and stq = falsestq=false and the matrix Q1 / 2Q1/2 input in q or with stq = truestq=true and b set to Q1 / 2Q1/2.
The algorithm requires (1/6)n3 + n2((3/2)m + l) + 2nm2 + (2/3)p316n3+n2(32m+l)+2nm2+23p3 operations and is backward stable (see Verhaegen and van Dooren (1986)). The transformation to lower observer Hessenberg form requires O((n + m)n2)O((n+m)n2) operations.

Example

function nag_tsa_multi_kalman_sqrt_invar_example
transf = 'T';
a = [0.607, -0.033, 1, 0, 0, 0;
     0, 0.543, 0, 1, 0, 0;
     0, 0, 0, 0, 0, 0;
     0, 0, 0, 0, 0, 0;
     0, 0, 0, 0, 1, 0;
     0, 0, 0, 0, 0, 1];
b = [1, 0;
     0, 1;
     0.543, 0.125;
     0.134, 0.026;
     0, 0;
     0, 0];
stq = false;
q = [1.612, 0;
     0.347, 2.282];
c = [1, 0, 0, 0, 1, 0;
     0, 1, 0, 0, 0, 1];
r = [0, 0;
     0, 0];
s = [2.8648, 0, 0, 0, 0, 0;
     0.7191, 2.729, 0, 0, 0, 0;
     0.5169, 0.2194, 0.781, 0, 0, 0;
     0.1266, 0.0449, 0.1899, 0.0098, 0, 0;
     0, 0, 0, 0, 0, 0;
     0, 0, 0, 0, 0, 0];
[aOut, bOut, cOut, sOut, k, h, u, ifail] = ...
    nag_tsa_multi_kalman_sqrt_invar(transf, a, b, stq, q, c, r, s)
 

aOut =

    0.8035   -0.0165    0.7341         0         0         0
         0    0.7715    0.0051    0.7431         0         0
    0.0526    0.0096   -0.1245   -0.0103    0.2587         0
   -0.0004    0.0702    0.0062   -0.1339    0.0207   -0.2917
    0.1887    0.0325   -0.4466   -0.0336    0.9273    0.0072
    0.0157   -0.2154   -0.0536    0.4132    0.0072    0.9060


bOut =

   -0.7071         0
         0   -0.7071
   -0.3338   -0.1045
   -0.1252    0.1934
    0.8279    0.0858
    0.0152   -0.6787


cOut =

   -1.4142         0         0         0         0         0
         0   -1.4142         0         0         0         0


sOut =

    1.2666         0         0         0         0         0
    0.2794    1.6137         0         0         0         0
    0.4511    0.2352   -0.3882         0         0         0
    0.1037   -0.4422   -0.0912   -0.0000         0         0
   -1.4634   -0.1949    0.1105    0.0000   -0.0000         0
    0.2262    1.5486   -0.0302   -0.0000   -0.0000    0.0000


k =

   -0.5425   -0.0335
   -0.0283   -0.3956
    0.1459    0.0179
    0.0077    0.1215
    0.5230    0.0611
    0.0163   -0.3726


h =

    2.8648         0
    0.7191    2.7290


u =

   -0.7071         0         0         0   -0.7071         0
         0   -0.7071         0         0         0   -0.7071
    0.1893    0.0159   -0.9632         0   -0.1893   -0.0159
   -0.0013    0.2173    0.0067   -0.9516    0.0013   -0.2173
    0.6790    0.0516    0.2685    0.0236   -0.6790   -0.0516
    0.0563   -0.6707    0.0000   -0.3065   -0.0563    0.6707


ifail =

                    0


function g13eb_example
transf = 'T';
a = [0.607, -0.033, 1, 0, 0, 0;
     0, 0.543, 0, 1, 0, 0;
     0, 0, 0, 0, 0, 0;
     0, 0, 0, 0, 0, 0;
     0, 0, 0, 0, 1, 0;
     0, 0, 0, 0, 0, 1];
b = [1, 0;
     0, 1;
     0.543, 0.125;
     0.134, 0.026;
     0, 0;
     0, 0];
stq = false;
q = [1.612, 0;
     0.347, 2.282];
c = [1, 0, 0, 0, 1, 0;
     0, 1, 0, 0, 0, 1];
r = [0, 0;
     0, 0];
s = [2.8648, 0, 0, 0, 0, 0;
     0.7191, 2.729, 0, 0, 0, 0;
     0.5169, 0.2194, 0.781, 0, 0, 0;
     0.1266, 0.0449, 0.1899, 0.0098, 0, 0;
     0, 0, 0, 0, 0, 0;
     0, 0, 0, 0, 0, 0];
[aOut, bOut, cOut, sOut, k, h, u, ifail] = g13eb(transf, a, b, stq, q, c, r, s)
 

aOut =

    0.8035   -0.0165    0.7341         0         0         0
         0    0.7715    0.0051    0.7431         0         0
    0.0526    0.0096   -0.1245   -0.0103    0.2587         0
   -0.0004    0.0702    0.0062   -0.1339    0.0207   -0.2917
    0.1887    0.0325   -0.4466   -0.0336    0.9273    0.0072
    0.0157   -0.2154   -0.0536    0.4132    0.0072    0.9060


bOut =

   -0.7071         0
         0   -0.7071
   -0.3338   -0.1045
   -0.1252    0.1934
    0.8279    0.0858
    0.0152   -0.6787


cOut =

   -1.4142         0         0         0         0         0
         0   -1.4142         0         0         0         0


sOut =

    1.2666         0         0         0         0         0
    0.2794    1.6137         0         0         0         0
    0.4511    0.2352   -0.3882         0         0         0
    0.1037   -0.4422   -0.0912   -0.0000         0         0
   -1.4634   -0.1949    0.1105    0.0000   -0.0000         0
    0.2262    1.5486   -0.0302   -0.0000   -0.0000    0.0000


k =

   -0.5425   -0.0335
   -0.0283   -0.3956
    0.1459    0.0179
    0.0077    0.1215
    0.5230    0.0611
    0.0163   -0.3726


h =

    2.8648         0
    0.7191    2.7290


u =

   -0.7071         0         0         0   -0.7071         0
         0   -0.7071         0         0         0   -0.7071
    0.1893    0.0159   -0.9632         0   -0.1893   -0.0159
   -0.0013    0.2173    0.0067   -0.9516    0.0013   -0.2173
    0.6790    0.0516    0.2685    0.0236   -0.6790   -0.0516
    0.0563   -0.6707    0.0000   -0.3065   -0.0563    0.6707


ifail =

                    0



PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

© The Numerical Algorithms Group Ltd, Oxford, UK. 2009–2013