# NAG CL Interfaceg13ecc (kalman_​sqrt_​filt_​info_​var)

## 1Purpose

g13ecc performs a combined measurement and time update of one iteration of the time-varying Kalman filter. The method employed for this update is the square root information filter with the system matrices in their original form.

## 2Specification

 #include
 void g13ecc (Integer n, Integer m, Integer p, Nag_ab_input inp_ab, double t[], Integer tdt, const double ainv[], Integer tda, const double b[], Integer tdb, const double rinv[], Integer tdr, const double c[], Integer tdc, const double qinv[], Integer tdq, double x[], const double rinvy[], const double z[], double tol, NagError *fail)
The function may be called by the names: g13ecc, nag_tsa_kalman_sqrt_filt_info_var or nag_kalman_sqrt_filt_info_var.

## 3Description

For the state space system defined by
 $X i+1 = A i X i + B i W i var W i = Q i Y i = C i X i + V i var V i = R i$
the estimate of ${X}_{i}$ given observations ${Y}_{1}$ to ${Y}_{i-1}$ is denoted by ${\stackrel{^}{X}}_{i\mid i-1}$ with $\mathrm{var}\left({\stackrel{^}{X}}_{i\mid i-1}\right)={P}_{i\mid i-1}={S}_{i}{S}_{i}^{\mathrm{T}}$. The function performs one recursion of the square root information filter algorithm, summarised as follows:
 $U 1 Q i - 1 / 2 0 Q i - 1 / 2 w ¯ i S i -1 A i -1 B i S i -1 A i -1 S i -1 X ^ i∣i 0 R i+1 - 1 / 2 C i+1 R i+1 - 1 / 2 Y i+1 = F i+1 - 1 / 2 * * 0 S i+1 -1 ξ i + 1∣i + 1 0 0 E i+1 Pre-array Post-array$
where ${U}_{1}$ is an orthogonal transformation triangularizing the pre-array. The triangularization is done entirely via Householder transformations exploiting the zero pattern of the pre-array. The term ${\overline{w}}_{i}$ is the mean process noise and ${E}_{i+1}$ is the estimated error at instant $i+1$. The inverse of the state covariance matrix ${P}_{i\mid i}$ is factored as follows
 $P i∣i -1 = S i -1 T S i -1$
where ${P}_{i\mid i}={S}_{i}{S}_{i}^{\mathrm{T}}$ (${S}_{i}$ is lower triangular).
The new state filtered state estimate is computed via
 $X ^ i + 1∣i + 1 = S i+1 ξ i + 1∣i + 1$
The function returns ${S}_{i+1}^{-1}$ and ${\stackrel{^}{X}}_{i+1\mid i+1}$ (see the G13 Chapter Introduction for more information concerning the information filter).

## 4References

Anderson B D O and Moore J B (1979) Optimal Filtering Prentice–Hall
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

## 5Arguments

1: $\mathbf{n}$Integer Input
On entry: the actual state dimension, $n$, i.e., the order of the matrices ${S}_{i}$ and ${A}_{i}^{-1}$.
Constraint: ${\mathbf{n}}\ge 1$.
2: $\mathbf{m}$Integer Input
On entry: the actual input dimension, $m$, i.e., the order of the matrix ${Q}_{i}^{-1/2}$.
Constraint: ${\mathbf{m}}\ge 1$.
3: $\mathbf{p}$Integer Input
On entry: the actual output dimension, $p$, i.e., the order of the matrix ${R}_{i+1}^{-1/2}$.
Constraint: ${\mathbf{p}}\ge 1$.
4: $\mathbf{inp_ab}$Nag_ab_input Input
On entry: indicates how the matrix ${B}_{i}$ is to be passed to the function.
${\mathbf{inp_ab}}=\mathrm{Nag_ab_prod}$
Array b must contain the product ${A}_{i}^{-1}{B}_{i}$.
${\mathbf{inp_ab}}=\mathrm{Nag_ab_sep}$
Then array b must contain ${B}_{i}$.
5: $\mathbf{t}\left[{\mathbf{n}}×{\mathbf{tdt}}\right]$double Input/Output
Note: the $\left(i,j\right)$th element of the matrix $T$ is stored in ${\mathbf{t}}\left[\left(i-1\right)×{\mathbf{tdt}}+j-1\right]$.
On entry: the leading $n$ by $n$ upper triangular part of this array must contain ${S}_{i}^{-1}$ the square root of the inverse of the state covariance matrix ${P}_{i\mid i}$.
On exit: the leading $n$ by $n$ upper triangular part of this array contains ${S}_{i+1}^{-1}$, the square root of the inverse of the of the state covariance matrix ${P}_{i+1\mid i+1}$.
6: $\mathbf{tdt}$Integer Input
On entry: the stride separating matrix column elements in the array t.
Constraint: ${\mathbf{tdt}}\ge {\mathbf{n}}$.
7: $\mathbf{ainv}\left[{\mathbf{n}}×{\mathbf{tda}}\right]$const double Input
Note: the $\left(i,j\right)$th element of the matrix is stored in ${\mathbf{ainv}}\left[\left(i-1\right)×{\mathbf{tda}}+j-1\right]$.
On entry: the leading $n$ by $n$ part of this array must contain ${A}_{i}^{-1}$ the inverse of the state transition matrix.
8: $\mathbf{tda}$Integer Input
On entry: the stride separating matrix column elements in the array ainv.
Constraint: ${\mathbf{tda}}\ge {\mathbf{n}}$.
9: $\mathbf{b}\left[{\mathbf{n}}×{\mathbf{tdb}}\right]$const double Input
Note: the $\left(i,j\right)$th element of the matrix $B$ is stored in ${\mathbf{b}}\left[\left(i-1\right)×{\mathbf{tdb}}+j-1\right]$.
On entry: the leading $n$ by $m$ part of this array must contain ${B}_{i}$ (if ${\mathbf{inp_ab}}=\mathrm{Nag_ab_sep}$) or its product with ${A}_{i}^{-1}$ (if ${\mathbf{inp_ab}}=\mathrm{Nag_ab_prod}$).
10: $\mathbf{tdb}$Integer Input
On entry: the stride separating matrix column elements in the array b.
Constraint: ${\mathbf{tdb}}\ge {\mathbf{m}}$.
11: $\mathbf{rinv}\left[{\mathbf{p}}×{\mathbf{tdr}}\right]$const double Input
Note: the $\left(i,j\right)$th element of the matrix is stored in ${\mathbf{rinv}}\left[\left(i-1\right)×{\mathbf{tdr}}+j-1\right]$.
On entry: if the measurement noise covariance matrix is to be supplied separately from the output weight matrix, then the leading $p$ by $p$ upper triangular part of this array must contain ${R}_{i+1}^{-1/2}$, the right Cholesky factor of the inverse of the measurement noise covariance matrix. If this information is not to be input separately from the output weight matrix (see below) then the array rinv must be set to NULL.
12: $\mathbf{tdr}$Integer Input
On entry: the stride separating matrix column elements in the array rinv.
Constraint: ${\mathbf{tdr}}\ge {\mathbf{p}}$ if rinv is defined.
13: $\mathbf{c}\left[{\mathbf{p}}×{\mathbf{tdc}}\right]$const double Input
Note: the $\left(i,j\right)$th element of the matrix $C$ is stored in ${\mathbf{c}}\left[\left(i-1\right)×{\mathbf{tdc}}+j-1\right]$.
On entry: if the array rinv has been set to NULL then the leading $p$ by $n$ part of this array must contain ${C}_{i+1}$, the output weight matrix (or its product with ${R}_{i+1}^{-1/2}$) of the discrete system at instant $i+1$.
14: $\mathbf{tdc}$Integer Input
On entry: the stride separating matrix column elements in the array c.
Constraint: ${\mathbf{tdc}}\ge {\mathbf{n}}$.
15: $\mathbf{qinv}\left[{\mathbf{m}}×{\mathbf{tdq}}\right]$const double Input
Note: the $\left(i,j\right)$th element of the matrix is stored in ${\mathbf{qinv}}\left[\left(i-1\right)×{\mathbf{tdq}}+j-1\right]$.
On entry: the leading $m$ by $m$ upper triangular part of this array must contain ${Q}_{i}^{-1/2}$ the right Cholesky factor of the inverse of the process noise covariance matrix.
16: $\mathbf{tdq}$Integer Input
On entry: the stride separating matrix column elements in the array qinv.
Constraint: ${\mathbf{tdq}}\ge {\mathbf{m}}$.
17: $\mathbf{x}\left[{\mathbf{n}}\right]$double Input/Output
On entry: this array must contain the estimated state ${\stackrel{^}{X}}_{i\mid i}$
On exit: the estimated state ${\stackrel{^}{X}}_{i+1\mid i+1}$.
18: $\mathbf{rinvy}\left[{\mathbf{p}}\right]$const double Input
On entry: this array must contain ${R}_{i+1}^{-1/2}{Y}_{i+1}$, the product of the upper triangular matrix ${R}_{i+1}^{-1/2}$ and the measured output ${Y}_{i+1}$.
19: $\mathbf{z}\left[{\mathbf{m}}\right]$const double Input
On entry: this array must contain ${\overline{w}}_{i}$, the mean value of the state process noise.
20: $\mathbf{tol}$double Input
On entry: tol is used to test for near singularity of the matrix ${S}_{i+1}^{-1}$. If you set tol to be less than ${n}^{2}×\epsilon$ then the tolerance is taken as ${n}^{2}×\epsilon$, where $\epsilon$ is the machine precision.
21: $\mathbf{fail}$NagError * Input/Output
The NAG error argument (see Section 7 in the Introduction to the NAG Library CL Interface).

## 6Error Indicators and Warnings

NE_2_INT_ARG_LT
On entry, ${\mathbf{tda}}=〈\mathit{\text{value}}〉$ while ${\mathbf{n}}=〈\mathit{\text{value}}〉$. These arguments must satisfy ${\mathbf{tda}}\ge {\mathbf{n}}$.
On entry, ${\mathbf{tdb}}=〈\mathit{\text{value}}〉$ while ${\mathbf{m}}=〈\mathit{\text{value}}〉$. These arguments must satisfy ${\mathbf{tdb}}\ge {\mathbf{m}}$.
On entry, ${\mathbf{tdc}}=〈\mathit{\text{value}}〉$ while ${\mathbf{n}}=〈\mathit{\text{value}}〉$. These arguments must satisfy ${\mathbf{tdc}}\ge {\mathbf{n}}$.
On entry, ${\mathbf{tdq}}=〈\mathit{\text{value}}〉$ while ${\mathbf{m}}=〈\mathit{\text{value}}〉$. These arguments must satisfy ${\mathbf{tdq}}\ge {\mathbf{m}}$.
On entry, ${\mathbf{tdr}}=〈\mathit{\text{value}}〉$ while ${\mathbf{p}}=〈\mathit{\text{value}}〉$. These arguments must satisfy ${\mathbf{tdr}}\ge {\mathbf{p}}$.
On entry, ${\mathbf{tdt}}=〈\mathit{\text{value}}〉$ while ${\mathbf{n}}=〈\mathit{\text{value}}〉$. These arguments must satisfy ${\mathbf{tdt}}\ge {\mathbf{n}}$.
NE_ALLOC_FAIL
Dynamic memory allocation failed.
On entry, argument inp_ab had an illegal value.
NE_INT_ARG_LT
On entry, ${\mathbf{m}}=〈\mathit{\text{value}}〉$.
Constraint: ${\mathbf{m}}\ge 1$.
On entry, ${\mathbf{n}}=〈\mathit{\text{value}}〉$.
Constraint: ${\mathbf{n}}\ge 1$.
On entry, ${\mathbf{p}}=〈\mathit{\text{value}}〉$.
Constraint: ${\mathbf{p}}\ge 1$.
NE_MAT_SINGULAR
The matrix inverse(S) is singular.

## 7Accuracy

The use of the square root algorithm improves the stability of the computations.

## 8Parallelism and Performance

g13ecc 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 function. Please also consult the Users' Note for your implementation for any additional implementation-specific information.

The algorithm requires approximately $\frac{7}{6}{n}^{3}+{n}^{2}\left(\frac{7}{2}m+p\right)+n\left(\frac{1}{2}{p}^{2}+{m}^{2}\right)$ operations and is backward stable (see Verhaegen and van Dooren (1986)).

## 10Example

To apply three iterations of the Kalman filter (in square root information form) to the system $\left({A}_{i}^{-1},{A}_{i}^{-1}{B}_{i},{C}_{i+1}\right)$. The same data is used for all three iterative steps.

### 10.1Program Text

Program Text (g13ecce.c)

### 10.2Program Data

Program Data (g13ecce.d)

### 10.3Program Results

Program Results (g13ecce.r)