PREPRINT

1

Likelihood Gradient Evaluation Using Square-Root
Covariance Filters

arXiv:1605.06654v1 [] 21 May 2016

M.V. Kulikova
Abstract— Using the array form of numerically stable square-root
implementation methods for Kalman filtering formulas, we construct
a new square-root algorithm for the log-likelihood gradient (score)
evaluation. This avoids the use of the conventional Kalman filter with
its inherent numerical instabilities and improves the robustness of
computations against roundoff errors. The new algorithm is developed
in terms of covariance quantities and based on the ”condensed form” of
the array square-root filter.
Index Terms— identification, maximum likelihood estimation, gradient
methods, Kalman filtering, numerical stability.

I. I NTRODUCTION
Consider the discrete-time linear stochastic system
xk

=

Fk xk−1 + Gk wk ,

(1)

zk

=

Hk xk + vk ,

(2)

n

k = 1, . . . , N,

m

where xk ∈ R and zk ∈ R are, respectively, the state and the measurement vectors; k is a discrete time, i.e. xk means x(tk ). The noises
wk ∈ Rq , vk ∈ Rm and the initial state x0 ∼ N (x̄0 , Π0 ) are taken
from mutually independent Gaussian distributions with zero mean and
covariance matrices Qk and Rk , respectively, i.e. wk ∼ N (0, Qk ),
vk ∼ N (0, Rk ). Additionally, system (1), (2) is parameterized by a
vector of unknown system parameters θ ∈ Rp , which needs to be
estimated. This means that the entries of the matrices Fk , Gk , Hk ,
Qk , Rk and Π0 are functions of θ ∈ Rp . However, for the sake of
simplicity we will suppress the corresponding notations below, i.e
instead of Fk (θ), Gk (θ), Hk (θ), Qk (θ), Rk (θ) and Π0 (θ) we will
write Fk , Gk , Hk , Qk , Rk and Π0 .
Solving the parameter estimation problem by the method of maximum likelihood requires the maximization of the likelihood function
(LF) with respect to unknown system parameters. It is often done by
using a gradient approach where the computation of the likelihood
gradient (LG) is necessary. For the state-space system (1), (2) the
negative Log LF is given as [1]:
N


o
1 X nm
−1
Lθ Z1N =
ln(2π) + ln (det Re,k ) + eTk Re,k
ek
2
2
k=1

Z1N

where
= [z1 , . . . , zN ] is N -step measurement history and ek are
the innovations, generated by the discrete-time Kalman filter (KF),
with zero mean and covariance matrix Re,k . They are ek = zk −
Hk x̂k|k−1 and Re,k = Hk Pk|k−1 HkT + Rk , respectively. The KF
defines the one-step ahead predicted state estimate x̂k|k−1 and the
one-step predicted error covariance matrix Pk|k−1 .
Straight forward differentiation of the KF equations is a direct
approach to the Log LG evaluation, known as a ”score”. This leads
to a set of p vector equations, known as the filter sensitivity equations,
for computing ∂ x̂k|k−1 /∂θ, and a set of p matrix equations, known
as the Riccati-type sensitivity equations, for computing ∂Pk|k−1 /∂θ.
Consequently, the main disadvantage of the standard approach
is the problem of numerical instability of the conventional KF, i.e
divergence due to the lack of reliability of the numerical algorithm.
Solution of the matrix Riccati equation is a major cause of numerical difficulties in the conventional KF implementation, from the
Manuscript received September 14, 2007; revised May 9, 2008.
M.V. Kulikova is with the School of Computational and Applied
Mathematics, University of the Witwatersrand, South Africa; Email:
Maria.Kulikova@wits.ac.za.

standpoint of computational load as well as from the standpoint of
computational errors [2].
The alternative approach can be found in, so-called, squareroot filtering algorithms. It is well known that numerical solution
of the Riccati equation tends to be more robust against roundoff
errors if Cholesky factors or modified Cholesky factors (such as
the U T DU -algorithms [3]) of the covariance matrix are used as the
dependent variables. The resulting KF implementation methods are
called square-root filters (SRF). They are now generally preferred
for practical use [2], [4], [5]. For more insights about numerical
properties of different KF implementation methods we refer to the
celebrated paper of Verhaegen and Van Dooren [6].
Increasingly, the preferred form for algorithms in many fields
is now the array form [7]. Several useful SRF algorithms for KF
formulas formulated in the array form have been recently proposed
by Park and Kailath [8]. For this implementations the reliability of
the filter estimates is expected to be better because of the use of
numerically stable orthogonal transformations for each recursion step.
Apart from numerical advantages, array SRF algorithms appear to be
better suited to parallel and to very large scale integration (VLSI)
implementations [8], [9].
The development of numerically stable implementation methods
for KF formulas has led to the hope that the Log LG (with respect
to unknown system parameters) might be computed more accurately.
For this problem, a number of questions arise:
• Is it possible to extend reliable array SRF algorithms to the case
of the Log LG evaluation?
• If such methods exist, will they inherit the advantages from
the source filtering implementations? In particular, will they
improve the robustness of the computations against roundoff
errors compared to the conventional KF technique? The question
about suitability for parallel implementation is beyond the scope
of this paper.
The first attempt to answer these questions belongs to Bierman
et al. [10]. The authors used the square-root information filter,
developed by Dyer and McReynolds [11] and later extended by
Bierman [3], as a source filter implementation and constructed the
method for score evaluation. The algorithm was developed in the
form of measurement and time updates. However, the accuracy of
the proposed method has not been investigated.
In contrast to the main result of [10], we focus on the dual class
of KF implementation methods (that is the class of covariance-type
methods) and discuss the efficient Log LG evaluation in square-root
covariance filters. More precisely, we consider the array form of the
square-root covariance filter eSRCF introduced in [8]. The purpose
of this paper is to design the method for the Log LG evaluation in
terms of the square-root covariance variables, i.e. in terms of the
quantities that appear naturally in the eSRCF. This avoids the use
of the conventional KF with its inherent numerical instabilities and
gives us an opportunity to improve the robustness of the Log LG
computation against roundoff errors.
II. E XTENDED S QUARE -ROOT C OVARIANCE F ILTER
To achieve our goal, we are first going to present the extended
square-root covariance filter (eSRCF), proposed in [8], and second,
we will derive the expression for the Log LG evaluation in terms of
the variables that are generated by the eSRCF implementation.
Notations to be used: For the sake of simplicity, we denote the
one-step predicted state estimate as x̂k and the one-step predicted
error covariance matrix as Pk . We use Cholesky decomposition of
T /2 1/2
1/2
the form Pk = Pk Pk , where Pk is an upper triangular matrix.
1/2
1/2
1/2
Similarly, we define Rk , Qk , Re,k . For convenience we will write

PREPRINT

2



1/2

Rk


 1/2
Qk  Pk HkT

0

1/2
Re,k

= 
0
0

0
1/2

Pk

FkT

1/2

Qk GTk
T
K̄p,k
1/2
Pk+1

0

−T /2

−Rk

−T /2

Pk

zk

x̂k

0
−ēk
−T /2
Pk+1 x̂k+1
γk

1/2

∂θi Rk


1/2
∂θi Pk HkT
0

Xi
Ni
Bi

A−1/2 = (A1/2 )−1 , A−T /2 = (A−1/2 )T and ∂θi A implies the
partial derivative of the matrix A with respect to the ith component
of θ, i.e ∂A/∂θi .
In this paper, we deal with the ”condensed form”1 of the eS−T /2
1/2
x̄0 , recurRCF [8]: Assume that Rk > 0. Given Π0 and Π0
1/2
−T /2
sively update Pk and Pk
x̂k as follows:


1/2
−T /2
Rk
0
−Rk
zk


−T /2
Qk  Pk1/2 HkT Pk1/2 FkT
Pk
x̂k 
1/2 T
0
0
Qk Gk


1/2
T
−ēk
Re,k K̄p,k


1/2
−T /2
(3)
= 0
Pk+1 Pk+1 x̂k+1 
γk
0
0
where Qk is any orthogonal rotation that upper-triangularizes the
first two (block) columns of the matrix on the left-hand side of (3);
−1/2
−T /2
K̄p,k = Fk Pk HkT Re,k and ēk = Re,k ek .
One can easily obtain the expression for the negative Log LF in
terms of the eSRCF variables:
N


o


1 X nm
1/2
Lθ Z1N =
ln(2π) + 2 ln det Re,k + ēTk ēk . (4)
2 k=1 2

Let θ = [θ1 , . . . , θp ] denote the vector of parameters with respect
to which the likelihood function is to be differentiated. Then from (4),
we have
N 
h 
i 1

 X
h
i
1/2
+ ∂θi ēTk ēk
∂θi ln det Re,k
∂θi Lθ Z1N =
. (5)
2
k=1
1/2
Re,k

Taking into account that the matrix
is upper triangular, we
derive
"
"m
#
#
m
h 
i
X
X  jj 
1
1/2
jj
∂θi ln det Re,k =∂θi
=
ln re,k
· ∂θi re,k
jj
re,k
j=1
j=1
i
h
−1/2
1/2
i = 1, . . . , p (6)
=tr Re,k · ∂θi Re,k ,

jj
where the re,k
, j = 1, . . . , m denote the diagonal elements of the
1/2
matrix Re,k .
Substitution of (6) into (5) yields the result that we are looking for
N n
i
o

h
 X
−1/2
1/2
∂θi Lθ Z1N =
tr Re,k · ∂θi Re,k + ēTk · ∂θi ēk . (7)
k=1

Ultimately, our problem is to compute Log LG (7) by using the
eSRCF equation (3). Before we come to the main result of this paper,
there are a few points to be considered. As can be seen from (7),
1/2
the elements ēk and Re,k involved in the Log LG evaluation are
obtained from the underlying filtering algorithm directly, i.e. from (3).
No additional computations are needed. Hence, our aim is to explain
1 The

0

”condensed form” of filtering algorithms refers to the case when implementation method for the KF formulas is not divided into the measurement
and time updates.

Yi
Vi
Ki





1/2
∂θi Pk FkT


1/2
∂θi Qk GTk



−T /2
zk
∂θi −Rk


−T /2
x̂k
∂θi Pk
0

Mi

Wi  .
Ti






(8)

1/2

how the last two terms in the Log LG expression, ∂θi ēk and ∂θi Re,k ,
can be computed using quantities available from eSRCF (3).
III. S UGGESTED S QUARE -ROOT M ETHOD FOR S CORE
E VALUATION
We can now prove the following result.
Theorem 1: Let the entries of the matrices Fk , Gk , Hk , Qk , Rk ,
Π0 describing the linear discrete-time stochastic system (1), (2) be
differentiable functions of a parameter θ ∈ Rp . Then in order to
compute the Log LF and its gradient (with respect to unknown system
parameter θ) the eSRCF, which is used to filter the data, needs to be
extended as follows. Assume that R
k > 0. Given
 the initial values
−T /2
1/2
−T /2
1/2
x̄0 , recursively update
x̄0 and ∂θi Π0 , ∂θi Π0
Π0 , Π0


−T /2
1/2
−T /2
1/2
x̂k as follows:
Pk , Pk
x̂k and ∂θi Pk , ∂θi Pk

I. Replace the eSRCF equation (3) by (8) where Qk is any orthogonal rotation that upper-triangularizes the first two (block)
columns of the matrix on the left-hand side of (8).
II. Having computed the elements of the right-hand side matrix
in (8), calculate for each θi :
#
#
"
"
1/2
T
T
h
i R1/2 K̄p,k
∂θi Re,k ∂θi K̄p,k
e,k
T
,
= L̄i + Di + Ūi
1/2
1/2
0 Pk+1
0
∂θi Pk+1
(9)
"

−∂θi ēk


−T /2

#
h
i
 = L̄T − L̄
i
i

−ēk
−T /2



Pk+1 x̂k+1
∂θi Pk+1 x̂k+1
#−T  
" 1/2
 
T
Re,k K̄p,k
Mi
Bi
(10)
γ
+
+
k
1/2
Wi
Ki
0 Pk+1

where L̄i , Di and Ūi are strictly lower triangular, diagonal and
strictly upper triangular parts of the following matrix product:
#
 " −1/2

−1/2 T
−1/2
Xi Yi Re,k −Re,k K̄p,k Pk+1
= L̄i + Di + Ūi .
−1/2
Ni Vi
0
Pk+1
(11)
1/2
1/2
III. Having determined ēk , Re,k and ∂θi ēk , ∂θi Re,k compute Log
LF (4) and Log LG (7).
Proof: As discussed earlier, the main difficulty in score evalu1/2
ation (7) is to define ∂θi Re,k and ∂θi ēk from the underlying filter,
i.e. from (3). We divide the proof into two parts, first proving (9) for
1/2
the ∂θi Re,k evaluation and then validating (10) for ∂θi ēk .
1/2
Part I. Our goal is to express ∂θi Re,k in terms of the variables
that appear naturally in the eSRCF implementation. First, we can
note that the eSRCF transformation in (3) has a form
QA = B
where A is a rectangular matrix, and Q is an orthogonal transformation that block upper-triangularizes B. If matrix A is square and

PREPRINT

3

daij
invertible, then given the matrix of derivatives A′θ =
we can
dθ
compute Bθ′ as follows [10]:
h
i
Bθ′ = LT + D + U B
(12)

where L, D and U are, respectively, strictly lower triangular, diagonal
and strictly upper triangular parts of the matrix QA′θ B −1 .
However, this idea cannot be applied to the eSRCF because the
matrix to be triangularized, i.e. the first two (block) columns of the
matrix on the left-hand side of (3), is not square and, hence, not invertible. By using the pseudoinversion (Moore-Penrose inversion) we
avoid this obstacle and generalize the scheme of computations (12)
to the case of eSRCF (3).
To begin constructing the method for score evaluation, we augment
the matrix to be triangularized by q columns of zeros. Hence, we
obtain

 

1/2
1/2
T
Rk
0
0
0
Re,k K̄p,k

 

1/2
Qk  Pk1/2 HkT Pk1/2 FkT 0  =  0
Pk+1 0  . (13)
1/2 T
0
0
0
0
Qk Gk 0

The matrices in (13) have dimensions (m + n + q) × (m + n + q).
For the sake of simplicity, we denote the left-hand side and the righthand side matrices of (13) as Ak and Bk , respectively. Then, by
differentiating (13) with respect to the components of θ, we obtain
∂θi Qk · Ak + Qk · ∂θi Ak = ∂θi Bk .

(14)

Multiplication both sides of (14) by the pseudoinverse matrix Bk+
yields

∂θi Bk · Bk+ =∂θi Qk Ak Bk+ + Qk · ∂θi Ak · Bk+


=∂θi Qk QTk Bk Bk+ + (Qk · ∂θi Ak ) Bk+ . (15)
One can easily obtain the explicit expression for Bk+ :

 −1/2
−1/2 T
−1/2
Re,k −Re,k K̄p,k
Pk+1 0


−1/2
Bk+ =  0
Pk+1
0 .

h
i
∂θi Qk · QTk

m+n

= L̄Ti − L̄i .

(19)

Formulas (19) and (18) are, in fact, equations (9) and (11) of the
proposed method for score evaluation. The theorem is half proved.
Part II. We need to verify (10). By differentiating the last equation
of the eSRCF with respect to the components of θ
 


−T /2
−ēk
−Rk
zk
−T /2


Qk  P −T /2 x̂k  = Pk+1 x̂k+1 
k
γk
0
we obtain


−∂θi ēk



−T /2
∂θi Pk+1 x̂k+1
∂θi γk





−T /2
−Rk
zk



T
 = ∂θi Qk · Qk · Qk  Pk−T /2 x̂k 


0 

−T /2
−∂θi Rk
zk




+Qk 
 ∂θi Pk−T /2 x̂k  . (20)


0

Next, we replace the last term in (20) with the quantities already
computed and collected in the right-hand side matrix of (8). Furthermore, it is useful to note that the element ∂θi γk is of no interest here.
These two steps give us
"
#


h
i
−∂θi ēk
−ēk

 = ∂ Q · QT
−T /2
−T /2
θi k
k
∂θi Pk+1 x̂k+1
m+n Pk+1 x̂k+1
 
h
irow: 1:m+n
Mi
(21)
+ ∂θi Qk · QTk
γk +
col: last q
Wi
irow: 1:m+n
h
where ∂θi Qk · QTk
stands for the (m + n) × q matrix
col: last q

(16)

0
0
0
By using (8), we replace Qk · ∂θi Ak in (15) by the quantities
already computed. Then, taking into account (16), we derive the
equation for the (m + n) × (m + n) main block of the matrix Bk :
# " 1/2 T #−1
"
1/2
T
h
i
Re,k K̄p,k
∂θi Re,k ∂θi K̄p,k
T
=
∂
Q
×
Q
θ
k
k
i
1/2
1/2
m+n
0 Pk+1
0
∂θi Pk+1
#−1
"


1/2
T
Xi Yi Re,k K̄p,k
(17)
+
1/2
Ni Vi
0 Pk+1


where ∂θi Qk · QTk m+n denotes the (m+n)×(m+n) main block
of the matrix ∂θi Qk · QTk .
As discussed in [10], the matrix ∂θi Qk · QTk is skew symmetric
and, hence, can be represented in the form L̄T − L̄ where L̄ is strictly
lower triangular.
Now, let us consider matrix equation (17). As can be seen, the
matrix on the left-hand side of (17) is block upper
 triangular. Thus,
the strictly lower triangular part of the matrix ∂θi Qk · QTk m+n
must exactly cancel the strictly lower triangular part of the second
term on the right-hand side of (17). In other words, if
#−1
 " 1/2

T
Xi Yi Re,k K̄p,k
= L̄i + Di + Ūi ,
1/2
Ni Vi
0 Pk+1
then

Substitution of (18) into (17) leads to the result
#
"
#
"
1/2
T
T
i R1/2 K̄p,k
h
∂θi Re,k ∂θi K̄p,k
e,k
T
.
= L̄i + Di + Ūi
1/2
1/2
0 Pk+1
0
∂θi Pk+1

(18)

composed of the entries that are located at the intersections of the
last q columns with the first m + n rows of ∂θi Qk · QTk .
Taking into account (18), from the equation above we obtain
"
#

h
i
−∂θi ēk
−ēk

 = L̄T − L̄
i
−T
/2
−T /2
i
Pk+1 x̂k+1
∂θi Pk+1 x̂k+1
 
irow: 1:m+n
h
Mi
(22)
γk +
+ ∂θi Qk · QTk
col: last q
Wi
where L̄i is strictly lower triangular part of the matrix in (11).
Since ∂θi Qk · QTk is skew symmetric, we can write down
h
irow: 1:m+n
irow: last q T
h
T
T
∂θi Qk · Qk
= − ∂θi Qk · Qk
(23)
col: last q

irow:

h

col: 1:m+n

last q

stands for the q × (m + n) matrix
where ∂θi Qk · QTk
col: 1:m+n
composed of the entries that are located at the intersections of the
last q rows with the first (m + n) columns of ∂θi Qk · QTk .
To evaluate the right-hand side of (23), we return to (15) and write
it in the matrix form:
+
  1/2

1/2
T
T
Re,k K̄p,k
0
∂θi Re,k ∂θi K̄p,k
0



T
1/2
1/2  
 0
0 P
0 = ∂θi Qk · Qk
∂θ P
0
i

0

k+1

k+1

0

0
0
+

  1/2
T
Xi Yi 0
Re,k K̄p,k
0


+ Ni Vi 0  0 P 1/2 0 . (24)
k+1
Bi Ki 0
0
0
0

0

0

PREPRINT

4

(26)

Final substitution of (26) into (22) validates (10) of the proposed
method for the Log LG evaluation. This completes the proof.
Remark 1: The method for score evaluation introduced above has
been derived from the eSRCF implementation. As a consequence, the
proposed method is of covariance-type.
Remark 2: The new square-root algorithm for score evaluation
naturally extends the eSRCF filter and, hence, consists of two parts.
They are the ”filtered” and ”differentiated” parts. This structure
allows the Log LF and its gradient to be computed simultaneously.
Thus, the method is ideal for simultaneous state estimation and
parameter identification.
Remark 3: In the KF formulation of the Log LG evaluation, it is
necessary to run the ”differentiated” KF for each of the parameters
θi to be estimated. As in [10], in the eSRCF formulation this ”bank”
of filters is replaced with the augmented arrays to which orthogonal
transformations are applied.

254

Negative Log Likelihood Function

By substituting (25) into (23), we obtain
" 1/2 T #−T  
h
irow: 1:m+n
Re,k K̄p,k
Bi
.
∂θi Qk · QTk
=
1/2
col: last q
K
i
0 Pk+1

256

252
250
248
246
244
242
240
238

0

First, we would like to check our theoretical derivations. To do so,
we apply the square-root algorithm introduced in Theorem 1 to the
following simple test problem.
Example 1: Consider the special case of the system (1), (2) being

 
  
dk
0
1
∆t
wk , zk = [1 0] xk + vk
x
+
=
xk =
k−1
−∆t/τ
sk
1
0 e
where wk ∼ N (0, I2 ), vk ∼ N (0, I1 ), In denotes the n × n identity
matrix and τ is a parameter which needs to be estimated.
In our simulation experiment, we compute the negative Log LF and
its gradient by the proposed square-root method and, then, compare
the results to those produced by the conventional KF approach. The
outcomes of this experiments are illustrated by Fig. 1 and 2.
As can be seen from Fig. 2, all algorithms for score evaluation
produce exactly the same result and give the same zero point that
further coincides with the minimum point of the negative Log LF (see
Fig. 1). All these evidences substantiate the theoretical derivations of
Section III.
Next, we wish to answer the second question posed in this paper:
does the algorithm for score evaluation derived from numerically
stable square-root implementation method improve the robustness
of computations against roundoff errors? The previously obtained
results (Example 1) indicate that both methods, i.e. the conventional
KF technique and the new square-root algorithm, produce exactly
the same answer for the Log LF and Log LG evaluation. However,
numerically they no longer agree. We are now going to explore the
accuracy of the numerical algorithms.
To begin designing the ill-conditioned test problem we, first, stress
the type of the proposed method. As discussed in Remark 1, the
new square-root algorithm belongs to the class of covariance-type
methods. From Verhaegen and Van Dooren’s celebrated paper [6], we
know that the condition number of the innovation covariance matrix
K(Re,k ) is the key parameter determining the numerical behavior of

5

8

10

15

20

τ (in sec.)

25

30

35

Fig. 1. The negative Log LF computed by the eSRCF and the conventional
KF for Example 1
2

0

−2

−4

−6

−8

−10

−12

Conventional KF technique
Suggested square−root method, the "differentiated" part
−14

IV. N UMERICAL R ESULTS

Conventional KF technique
Suggested square−root method, the "filtered" part

236

Likelihood Gradient

As can be seen, the last (block) row of the left-hand side matrix
in (24) is zero. Thus, the last (block) row of the matrix ∂θi Qk · QTk
must exactly cancel the last (block) row of the second term in (24):
#−1
" 1/2
T
h
irow: last q
Re,k K̄p,k
T
∂θi Qk · Qk
. (25)
= − [Bi Ki ]
1/2
col: 1:m+n
0 Pk+1

0

5

8

10

15

20

τ (in sec.)

25

30

35

Fig. 2. The Log LG computed by the proposed square-root method and the
conventional KF for Example 1

the covariance algorithms. Taking into account these two important
facts, we construct the following ill-conditioned test problem.
Example 2: Consider the problem with the measurement sensitivity matrix


1 1 1
Hk =
and Fk = I3 , Gk = 0, Qk = I1 , Rk = δ 2 θI2
1 1 1+δ
with x0 ∼ N (0, θI3 ), where θ is an unknown system parameter. To
simulate roundoff we assume that δ 2 < ǫroundof f , but δ > ǫroundof f
where ǫroundof f denotes the unit roundoff error2 .
When θ = 1, Example 2 coincides with well-known ill-conditioned
filtering problem (see, for instance, [2]) and demonstrates how
a problem that is well-conditioned, as posed, can be made illconditioned by the filter implementation. The difficulty to be explored
is in matrix inversion. As can be seen, although rank H = 2, the
matrix Re,1 is singular in machine precision that yields the failure
of the conventional KF implementation. We introduced an unknown
system parameter θ making sure that the same problem is applied to
the matrix (Re,1 )′θ for each value of θ. Thus, both parts of the method
for score evaluation, that are the ”filtered” and ”differentiated” parts,
fail after processing the first measurement. From the discussion above
we understand that Example 2 demonstrates the difficulty only for
the covariance-type methods.
Our simulation experiments presented below are organized as
follows. All methods were implemented in the same precision (64-bit
2 Computer roundoff for floating-point arithmetic is often characterized by a
single parameter ǫroundof f , defined in different sources as the largest number
such that either 1 + ǫroundof f = 1 or 1 + ǫroundof f /2 = 1 in machine
precision.

PREPRINT

5

TABLE I
E FFECT OF ROUNDOFF ERRORS ON THE COMPUTED SOLUTIONS FOR THE SET OF TEST PROBLEMS FROM E XAMPLE 2

K(Re,1 )

10−2

103

10−4

107

10−6
10−8
10−9
10−10

1011
1015
1016
∞

Conventional KF technique
∆P1
10−13

1·
5 · 10−10
2 · 10−6
3 · 10−3
3 · 10−1
NaN

∆P1′
10−10

1·
9 · 10−4
2 · 10−1
2 · 10−1
NaN
NaN

∆LogLG

10−13

10−13

2·
4 · 10−9
2 · 10−5
3 · 10−1
4 · 100
NaN

floating point) in MatLab where the unit roundoff error is 2−53 ≈
1.11 · 10−16 . The MatLab function eps is twice the unit roundoff
error and δ = eps2/3 satisfies the conditions δ 2 < ǫroundof f and
δ > ǫroundof f from Example 2. We provide the computations for
for different values of δ, say δ ∈ [10−9 eps2/3 , 109 eps2/3 ]. This
means that we consider a set of test problems from Example 2. The
unknown system parameter θ is fixed, say θ = 2. The exact answers
are produced by the Symbolic Math Toolbox of MatLab.
Experiment 1: In this experiment we are going to use the performance profile technique to compare the conventional KF approach
for score evaluation with the square-root algorithm introduced in this
paper. The performance profile method was developed by Dolan and
Moré [12] to answer a common question in scientific computing:
how to compare several competing methods on a set of test problem.
Now, it can be found in textbooks (see, for instance, [13]).
In our simulation experiments we consider a set A of n = 2
algorithms, mentioned above. The performance measure, ta (p), is a
measure of accuracy. More precisely, ta (p) is the maximum absolute
error in Log LG computed for 7 different values of δ. Thus, we
consider a set P of m = 7 test problems from Example 2;
δ ∈ [10−2 , 10−3 , 10−4 , 10−5 , 10−6 , 10−7 , 10−8 ]. According to the
performance profile technique, we compute the performance ratio
rp,a =

ta (p)
≥ 1,
min{tσ (p) : σ ∈ A}

which is the performance of algorithm a on problem p divided
by the best performance of all the methods (we mean a particular
implementation method for score evaluation) on this problem. The
performance profile of algorithm a is the function
1
× number of p ∈ P such that rp,a ≤ µ,
m
which is monotonically increasing. Thus, φa (µ) is the probability
that the performance of algorithm a is within a factor µ of the best
performance over all implementations on the given set of problems.
The results of this experiment are illustrated by Fig. 3. For each
method, µ is plotted against the performance profile φa (µ), for µ ∈
[0, 3]. We are now going to explain Fig. 3.
Let us consider the left-hand side of Fig. 3, where µ = 1. We
can say that the new square-root algorithm proposed in this paper is
the most accurate implementation on ≈ 71% of the problems, with
the conventional KF being accurate on 30% of the problems. Next,
we consider the middle of the plot, looking where the curve first hit
probability 1. We conclude that the suggested square-root method is
within a factor µ ≈ 1.3 of being the most accurate implementation on
every test problem. However, the conventional KF approach for score
evaluation will never manage all 7 problems (as δ → ǫroundof f , the
machine precision limit, the test problems become ill-conditioned).
We need to increase µ to ≈ 2.7 to be able to say that for ≈ 58% of
the test problems the conventional KF provides an accurate Log LG
evaluation within a factor µ ≈ 2.7.
φa (µ) =

Suggested square-root method

∆LogLF

∆P1
10−15

1·
1 · 10−9
6 · 10−6
2 · 10−2
4 · 100
NaN

4·
4 · 10−13
3 · 10−11
3 · 10−10
2 · 10−8
2 · 10−7

∆P1′

∆LogLF

∆LogLG

7 · 10−16
7 · 10−14
1 · 10−11
2 · 10−10
7 · 10−9
1 · 10−8

1 · 10−13
6 · 10−10
9 · 10−6
2 · 10−1
1 · 100
2 · 104

9 · 10−14
7 · 10−10
4 · 10−6
9 · 10−3
5 · 101
2 · 104

1
0.9
0.8
0.7

a

δ

Performance profile φ (µ)

Problem conditioning

0.6
0.5
0.4
0.3
0.2
0.1
0
1

Suggested square−root algorithm
Conventional KF technique

1.2

1.4

1.6

1.8

2
2.2
Factor, µ

2.4

2.6

2.8

3

Fig. 3.
Performance profiles of the methods for score evaluation: the
conventional KF implementation and the new square-root algorithm proposed
in this paper, – on the set of test problems from for Example 2.

Thus, the performance profiles clearly indicate that on the set of
the test problems from Example 2 the new square-root algorithm
derived in this paper provides more accurate evaluation of the Log
LG compared with the conventional KF approach.
Experiment 2: In this experiment we use the conventional KF
technique and the proposed square-root method to compute the
maximum absolute error in Log LF, denoted as ∆LogLF , and its
gradient, denoted as ∆LogLG. The results of this experiment are
summarized in Table I. We also present the maximum absolute error
among elements in matrices P1 and (P1 )′θ (denoted as ∆P1 and ∆P1′ ,
respectively) to explore the numerical behavior of the ”filtered” and
”differentiated” parts of the methods for score evaluation.
As can be seen from Table I, the square-root implementation of
the Riccati-type sensitivity equation degrades more slowly than the
conventional Riccati-type sensitivity recursion as δ → ǫroundof f ,
the machine precision limit (see columns denoted as ∆P1′ ). For
instance, the ”filtered” (columns ∆P1 ) and ”differentiated” (columns
∆P1′ ) parts of the proposed square-root method for score evaluation
maintain about 7 and 8 digits of accuracy, respectively, at δ = 10−9 .
The conventional KF technique provides essentially no correct digits
in both computed solutions. Besides, it seems that the roundoff errors
tend to accumulate and degrade the accuracies of the Log LF and Log
LG faster than the accuracies of ∆P1 and ∆P1′ . Indeed, for the same
δ = 10−9 we obtain no correct digits in the computed solutions for
all methods. In MatLab, the term ’NaN’ stands for ’Not a Number’
that actually means the failure of the numerical algorithm.
Remark 4: The results of Experiment 2 indicate that the new
square-root algorithm provides more accurate computation of the
sensitivity matrix (Pk )′θ compared to the conventional KF. Hence,
it can be successfully used in all applications where this quantity is
required.

PREPRINT

6

V. C ONCLUDING R EMARKS
In this paper, a numerically stable square-root implementation
method for KF formulas, the eSRCF, has been extended in order
to compute the Log LG for linear discrete-time stochastic systems.
The preliminary analysis indicates that the new algorithm for score
evaluation provides more accurate computations compared with the
conventional KF approach. The new result can be used for efficient
calculations in sensitivity analysis and in gradient-search optimization algorithms for the maximum likelihood estimation of unknown
system parameters.
As an extension of the eSRCF, the new method for score evaluation
is expected to inherit its benefits. However, the question about
suitability for parallel implementation is still open.
It can be mentioned that another approach to construct numerically
stable implementation method for score evaluation is to use the U D
filter [3]. Being the modification of the square-root implementations,
the U D-type algorithms improve the robustness of computations
against roundoff errors, but compared with SRF, the U D filter reduces
the computational cost (see [3], [6], [5]). As mentioned in [10] and
as far as this author knows, it is still not known how to use the U D
filter to compute the score.
R EFERENCES
[1] F.C. Schweppe, ”Evaluation of likelihood functions for Gaussian signals”, IEEE Trans. Inf. Theory, vol. IT-11, pp. 61–70, Jan. 1965.
[2] M.S. Grewal and A.P. Andrews, Kalman Filtering: Theory and Practice
Using MATLAB, 2nd ed., Ed. New York: John Wiley & Sons, 2001.
[3] G.J. Bierman, Factorization Methods for Discrete Sequential Estimation.
Ed. New York: Academic Press, 1977.
[4] P.G. Kaminski, A.E. Bryson and S.F. Schmidt, ”Discrete square-root
filtering: a survey of current techniques”, IEEE Trans. Autom. Control,
vol. AC-16, pp. 727–735, Dec. 1971.
[5] T. Kailath, A.H. Sayed and B. Hassibi, Linear Estimation, Ed. New
Jersey: Prentice Hall, 2000.
[6] M. Verhaegen and P.V. Dooren, ”Numerical aspects of different Kalman
filter implementations”, IEEE Trans. Autom. Control, vol. AC-31,
pp. 907–917, Oct. 1986.
[7] T. Kailath, ”Array algorithms for structured matrices”, presented at the
conference of the International Linear Algebra Society, Winnipeg, USA,
June 3-6, 1998.
[8] P. Park and T. Kailath, ”New square-root algorithms for Kalman filtering”, IEEE Trans. Autom. Control, vol. 40, pp. 895–899, May 1995.
[9] E.K.B. Lee and S. Maykin, ”Parallel implementation of the extended
square-root covariance filter for tracking applications”, IEEE Trans.
Parallel Distrib. Syst., vol. 4, pp. 446–457, April 1993.
[10] G.J. Bierman, M.R. Belzer, J.S. Vandergraft and D.W. Porter, ”Maximum
likelihood estimation using square root information filters”, IEEE Trans.
Autom. Control, vol. 35, pp. 1293–1298, Dec. 1990.
[11] P. Dyer and S. McReynolds, ”Extensions of square root filtering to
include process noise”, J. Opt. Theory Appl., vol. 3, pp. 444–459,
June 1969.
[12] E.D. Dolan and J.J. Moré, ”Benchmarking optimization software with
performance profiles”, Math. Programming, vol. 91:201-213, pp. 320–
348, 2002.
[13] D.J. Higham and N.J. Higham, MatLab Guide, 2nd ed., Ed. Philadelphia:
SIAM, 2005.

