PREPRINT

1

SVD-based Kalman Filter Derivative Computation

arXiv:1612.04777v2 [cs.SY] 6 Apr 2017

J.V. Tsyganova and M.V. Kulikova

Abstract— Recursive adaptive filtering methods are often used for
solving the problem of simultaneous state and parameters estimation
arising in many areas of research. The gradient-based schemes for
adaptive Kalman filtering (KF) require the corresponding filter sensitivity
computations. The standard approach is based on the direct differentiation of the KF equations. The shortcoming of this strategy is a numerical
instability of the conventional KF (and its derivatives) with respect to
roundoff errors. For decades, special attention has been paid in the KF
community for designing efficient filter implementations that improve
robustness of the estimator against roundoff. The most popular and
beneficial techniques are found in the class of square-root (SR) or UD
factorization-based methods. They imply the Cholesky decomposition of
the corresponding error covariance matrix. Another important matrix
factorization method is the singular value decomposition (SVD) and,
hence, further encouraging KF algorithms might be found under this
approach. Meanwhile, the filter sensitivity computation heavily relies
on the use of matrix differential calculus. Previous works on the
robust KF derivative computation have produced the SR- and UD-based
methodologies. Alternatively, in this paper we design the SVD-based
approach. The solution is expressed in terms of the SVD-based KF
covariance quantities and their derivatives (with respect to unknown
system parameters). The results of numerical experiments illustrate
that although the newly-developed SVD-based method is algebraically
equivalent to the conventional approach and the previously derived SRand UD-based strategies, it outperforms the mentioned techniques for
estimation accuracy in ill-conditioned situations.
Index Terms— Kalman filter, filter sensitivity equations, SVD factorization, array algorithms.

I. I NTRODUCTION
The problem of filter sensitivities evaluation plays a key role
in many areas of research; for instance, in state estimation and
parameter identification realm [1], [2], in the field of optimal input
design [3], [4], in information theory for computing the Fisher
information matrix [5]–[7] etc. In this paper we explore linear
discrete-time stochastic systems where the associated Kalman filter
(KF) is used for estimating the unknown dynamic states. Therefore,
the standard approach for computing the filter sensitivities (with
respect to unknown system parameters) is a direct differentiation of
the KF equations. This conventional methodology is comprehensively
studied in [3], [8], [9]. The shortcoming of this strategy is a numerical
instability of the conventional KF (and its derivatives) with respect
to roundoff errors discussed in [10], [11]. Due to this fact, special
attention has been paid in the KF community for designing robust
KF implementation methods. The most popular techniques belong
to the class of square-root (SR) or UD factorization-based methods;
see [12]–[15] and many others. These algorithms imply the Cholesky
decomposition and its modification for the corresponding covariance
matrix factorization [13], [16], [17]. We may note that the Cholesky
decomposition exists and is unique when the symmetric matrix to
be decomposed is positive definite [18]. If it is a positive semidefinite, then the Cholesky decomposition still exists, however, it
Manuscript received ??; revised ??. The second author thanks the support of
Portuguese National Fund (Fundação para a Ciência e a Tecnologia) within
the scope of project UID/Multi/04621/2013.
The first author is with Ulyanovsk State University, Str. L. Tolstoy 42,
432017 Ulyanovsk, Russian Federation. The second author is with CEMAT
(Center for Computational and Stochastic Mathematics), Instituto Superior
Técnico, Universidade de Lisboa, Av. Rovisco Pais 1, 1049-001 LISBOA,
Portugal; Emails: TsyganovaJV@gmail.com; maria.kulikova@ist.utl.pt

is not unique [19]. Further encouraging KF implementation methods
might be found with the use of singular value decomposition (SVD).
Some evidences of better estimation quality obtained under the SVDbased approach exist in the field of nonlinear filtering; for instance,
see discussion in [20]–[22] and others. For linear filtering problem
examined in this paper, the first SVD-based KF was, to the best of
our knowledge, designed in [23]. Our recent analysis exposes that the
mentioned SVD-based filter can be further improved for enhancing its
numerical robustness. This result is comprehensively studied in [24],
where some new stable SVD-based KF implementations are designed.
Despite the existence of inherently more stable SR-, UD- and SVDbased KF variants, the problem of robust filter derivative computation
is seldom addressed in practice because of its complicated matter.
The solution to the mentioned problem heavily relies on the use
of matrix differential calculus. The first SR-based information-type
algorithm for the KF derivative computations belongs to Bierman et
al. and was appeared in 1990; see [25]. Alternatively, the SR-based
covariance-type method was proposed in [26] as well as the UDbased scheme designed in [27]. Later on, a general “differentiated”
SR-based methodology was designed for both orthogonal and Jorthogonal transformations involved in the filtering equations (and
their derivatives) in [28]–[30]. Alternatively, in this technical note we
develop the SVD-based approach for the KF derivative computation.
We show that the new technique is algebraically equivalent to the
conventional “differentiated” KF, but it improves the robustness
against roundoff errors as well as the existing “differentiated” SR- and
UD-based methodologies. However, motivated by the results obtained
in nonlinear filtering realm, we expect that the newly-designed SVDbased method outperforms the previously derived algorithms while
solving the parameters estimation problem, especially when the error
covariance is ill-conditioned.

II. F ILTER SENSITIVITY EQUATIONS : CONVENTIONAL APPROACH
Consider the state-space equations
xk = F (θ)xk−1 + B(θ)uk−1 + G(θ)wk−1 ,

k ≥ 1,

zk = H(θ)xk + vk , vk ∼ N (0, R(θ)) , wk ∼ N (0, Ω(θ))

(1)
(2)

where zk ∈ Rm , uk ∈ Rd , xk ∈ Rn and θ ∈ Rp are, respectively,
the vectors of available measurements, the known deterministic
control input, the unknown dynamic state and the unknown system
parameters that need to be estimated from the available experimental
data, {z1 , . . . , zN }. The process and the measurement noises are
independent Gaussian zero-mean white-noise processes that are also
independent from the initial state x0 ∼ N (x̄0 , Π0 (θ)). The covariances are assumed to be Ω(θ) ≥ 0, R(θ) > 0 and Π0 (θ) ≥ 0.
Equations (1), (2) represent a set of the state-space models (SSMs).
Each of them corresponds to a particular system parameter value. This
means that for any fixed value of θ, say θ̂∗ , the system matrices are
known, i.e. there is no uncertainty in model (1), (2). For simplicity,
throughout the paper we write F etc. instead of F (θ̂∗ ) etc. when
evaluating at the fixed point θ̂∗ . The associated KF yields the linear
minimum least-square estimate of the unknown dynamic state that

PREPRINT

2

can be recursively computed via the equations [16, Theorem 9.2.1]:
ek = zk − H x̂k|k−1 ,
Kp,k = F Pk|k−1 H

T

x̂1|0 = x̄0 ,

−1
Re,k
,

k ≥ 1,

(3)
T

Re,k = R + HPk|k−1 H ,

x̂k+1|k = F x̂k|k−1 + Buk + Kp,k ek

(4)
(5)

where {ek } are innovations of the discrete-time KF. The important
property of the
 KF for Gaussian SSMs is ek ∼ N (0, Re,k ). The
Pk|k−1 = E (xk − x̂k|k−1 )(xk − x̂k|k−1 )T is the one-step ahead
predicted error covariance matrix computed as follows:
T
Pk+1|k = F Pk|k−1 F T + GΩGT − Kp,k Re,k Kp,k
, P1|0 = Π0 . (6)

The conventional approach for deriving the related sensitivity
model is based on differentiation of the corresponding filtering
equations. Let A(θ) ∈ Rm×n , B(θ) ∈ Rn×q be matrices, which
entries are differentiable functions of the parameter vector θ ∈ Rp .
The m×n matrix ∂i A = ∂A/∂θi implies the partial derivative of the
A with respect
Pto the i-th component of θ, i = 1, . . . p. The m × n
matrix dA = pi=1 (∂i A) · (dθi ) is the differential form of first-order
derivatives of A(θ). Taking into account the matrix product rule of
differentiation [31, p. 955]: d (AB)
= (dA) B + A (dB), and the

fact dI = 0, we derive d A−1 = −A−1 (dA) A−1 for any square
and invertible matrix A (it is also known as the Jacobi’s formula);
see also [8, p. 546]. Using these differentiation rules, the necessary
differentials of (3)-(6) can be written as follows [8], [9]:


dek = − (dH) x̂k|k−1 + H dx̂k|k−1 ,
(7)

dx̂k+1|k = (dF ) x̂k|k−1 + F dx̂k|k−1 + (dB) uk
+ (dKp,k ) ek + Kp,k (dek ),
T

(8)



−1
Re,k

+ F dPk|k−1 H
dKp,k = (dF ) Pk|k−1 H


T
−1
+ F Pk|k−1 dH
Re,k

T

−1
Re,k

−1
−1
,
(dRe,k ) Re,k
− F Pk|k−1 H T Re,k
T



(9)
T

dRe,k = dR + (dH) Pk|k−1 H + H dPk|k−1 H


(10)
+ HPk|k−1 dH T ,
 T
T
dPk+1|k = (dF ) Pk|k−1 F + F dPk|k−1 F


+ F Pk|k−1 dF T + (dG) ΩGT + G (dΩ) GT


T
+ GΩ dGT − (dKp,k ) Re,k Kp,k


T
T
− Kp,k (dRe,k ) Kp,k
− Kp,k Re,k dKp,k
.
(11)

In deriving the equations above we take into account that dzk = 0
and duk = 0, because the observations zk and the control input uk do
not depend on the parameters (i.e. their realizations are independent
of variations in θ) and therefore have a differential equal to zero.
We may also note that except for the scalar factor dθi , ∂i A is a
special case of dA, so that to obtain partial-derivative forms from
differential forms, we only have to everywhere replace operator d (·)
with ∂i (·) for i = 1, . . . p [8, p. 546]. Hence, from (7) – (11) we
obtain a set of p vector equations, known as the filter sensitivity
equations, for computing ∂i x̂k+1|k , i = 1, . . . p, and a set of p
matrix equations, known as the Riccati-type sensitivity equations,
for computing ∂i Pk+1|k , i = 1, . . . p. This approach for the KF
sensitivity model derivation is called the “differentiated KF”. Its main
drawback is a numerical instability of the conventional KF (3) – (6)
and inherently its derivative (7) – (11) with respect to roundoff errors.
The goal of this paper is to design a robust methodology for
updating the “differentiated” KF equations above in terms of SVD
factors (and their derivatives) of the error covariance matrices Pk|k−1
instead of using the full matrices Pk|k−1 (and their derivatives).

III. SVD FACTORIZATION - BASED K ALMAN FILTERING
To the best of our knowledge, the first SVD-based KF was by Wang
et al. and appeared in 1992; see Eqs (17), (22), (23) in [23, pp. 12251226]. Our recent research shows that although that implementation
is inherently more stable than the KF (3) – (6), it is still sensitive
to roundoff and poorly treats ill-conditioned problems. The cited
analysis exposes that the SVD-based filter can be further improved for
enhancing its numerical robustness. This result is comprehensively
studied in [24], where new stable SVD-based KF implementations
are designed. The readers are referred to the cited paper for the
detailed derivations, numerical stability discussion and proofs. Here,
we briefly outline the principle steps for construction of the most
advanced SVD-based KF variant. Next, we extend it to a stable filter
sensitivities computation, which is the main purpose of this study.
Consider the SVD factorization [32, Theorem 2.8.1]: suppose A ∈
Cm×n , rank A = r. There exist positive numbers σ1 ≥ . . . σr > 0
and unitary matrices W ∈ Cm×m and V ∈ Cn×n such that


S 0
∗
∈ Cm×n , S = diag{σ1 , . . . , σr }
A = W ΣV , Σ =
0 0

where V ∗ is the conjugate transpose of V .
The diagonal entries of Σ are known as the singular values of A.
The non-zero σi (i = 1, . . . , r) are the square roots of the non-zero
eigenvalues of both A∗ A and AA∗ .
If A is a square matrix such that A∗ A = AA∗ , then the A can be
diagonalized using a basis of eigenvectors according to the spectral
theorem, i.e. it can be factorized as follows: A = QDQ∗ where Q
is a unitary matrix and D is a diagonal matrix, respectively. If A
is also positive semi-definite, then the spectral decomposition above,
A = QDQ∗ , is also a SVD factorization, i.e. the diagonal matrix
D contains the singular values of A. For the SSMs examined in this
paper, the initial error covariance Π0 ∈ Rn is a symmetric positive
semi-definite matrix and, hence, the spectral decomposition implies
Π0 = QΠ0 DΠ0 QTΠ0 where QΠ0 and DΠ0 are the orthogonal and
diagonal matrices, respectively. It is also a SVD factorization, i.e.
the factor DΠ0 contains the singular values of Π0 .
Now, we are ready to present the SVD-based KF implementation
developed recently in [24]. Instead of conventional recursion (3)-(6)
1/2
for Pk|k−1 , we update only their SVD factors, {QPk|k−1 , DPk|k−1 },
at each iteration step of the filter as shown below.
I NITIAL S TEP (k = 0). Apply the SVD factorization for the initial
error covariance matrix Π0 = QΠ0 DΠ0 QTΠ0 and, additionally, for
the process and measurement noise covariances: Ω = QΩ DΩ QTΩ
and R = QR DR QTR , respectively. Set the initial values as follows:
1/2
1/2
QP1|0 = QΠ0 , DP1|0 = DΠ0 and x̂1|0 = x̄0 .
M EASUREMENT U PDATE (k = 1, . . . , N ). Build the pre-arrays
from the filter quantities that are currently available and, then, apply
the SVD factorizations in order to obtain the corresponding SVD
factors of the updated filter quantities as follows:
#
"
#
"
1/2
1/2
DR QTR
DRe,k
(1)
QTRe,k , (12)
=
W
1/2
MU
DPk|k−1 QTPk|k−1 H T
0
|
{z
} |
{z
}
Pre−array
Post−array SVD factors


K̄k = QPk|k−1 DPk|k−1 QTPk|k−1 H T QRe,k , (13)
#
"
#
" 1/2
1/2
DPk|k−1 QTPk|k−1 (I − Kk H)T
DPk|k
(2)
QTPk|k (14)
= WM U
1/2
0
DR QTR KkT
{z
} |
{z
}
|
Pre−array

Post−array SVD factors

(1)
−1
where we denote Kk = K̄k DR
QTRe,k . The matrices WM U ∈
e,k
(2)
m×m
(n+m)×(n+m)
(m+n)×(m+n)
and WM U ∈ R
,
R
, QRe,k ∈ R
QPk|k ∈ Rn×n are the orthogonal matrices of the corresponding

PREPRINT

3

1/2

1/2

SVD factorizations in (12), (14). Next, DRe,k ∈ Rm×m and DPk|k ∈
Rn×n are diagonal matrices with square roots of the singular values
of Re,k and Pk|k , respectively.
It can be easily seen that the required SVD factors of the inno1/2
vation covariance Re,k , i.e. {QRe,k , DRe,k }, and a posteriori error
1/2

covariance matrix Pk|k , i.e. {QPk|k , DPk|k }, are directly read-off
from the post-array factors in (12) and (14), respectively. Finally,
find a posteriori estimate x̂k|k through equations

−1
x̂k|k = x̂k|k−1 + K̄k DR
ē , ēk = QTRe,k zk − H x̂k|k−1 . (15)
e,k k
T IME U PDATE (k = 1, . . . , N ). Build the pre-array and apply the
SVD factorization to obtain a priori error covariance SVD factors
1/2
{QPk+1|k , DPk+1|k } as follows:
" 1/2
#
DPk|k QTPk|k F T
1/2

DΩ QTΩ GT
|
{z
Pre−array

}

= WT U
|

"

1/2
DPk+1|k

0
{z

#

QTPk+1|k

Post−array SVD factors

and find a priori estimate x̂k+1|k as follows:

x̂k+1|k = F x̂k|k + Buk .

(16)

}

the derivative of the pre-array, A′θ , the following formulas calculate
the corresponding derivatives of the post-arrays:
 ′
h
i
S
Σ′θ = θ ,
Sθ′ = diag WT A′θ V
,
(19)
0
s×s
h
i
V′θ = V L̄T2 − L̄2
(20)



where WT A′θ V s×s denotes the main s × s block of the matrix
product WT A′θ V, and L̄2 is a strictly lower triangular matrix, which
entries are computed as follows:
(l̄2 )ij =

ūji σj + ¯
lij σi
,
σi2 − σj2

i = 2, . . . , s, j = 1, . . . , i − 1.

(21)

In equation above, the quantities ūji and l̄ji denote the entries of
matrices L̄ and Ū , respectively. The L̄,Ū are strictly
 lower and upper
triangular parts of the matrix product WT A′θ V s×s , respectively.
Proof: By differentiating (18) with respect to θ, we obtain
A′θ = W′θ ΣVT + WΣ′θ VT + WΣ (VT )′θ .

(22)
T

(17)

The SVD-based KF implementation above is formulated in twostage form. Meanwhile, following [15], the conventional KF (3) –
(6) is expressed in the so-called “condensed” form. Nevertheless,
these KF variants are algebraically equivalent. It is easy to prove
if we take into account the SVD factorization A = WΣVT and
the properties of orthogonal matrices. Indeed, for each pre-array to
be decomposed we have AT A = (VΣWT )(WΣVT ) = VΣ2 VT .
Next, by comparing both sides of the obtained matrix equations, we
come to the corresponding SVD-based KF formulas. The detailed
derivation can be found in [24].
IV. F ILTER SENSITIVITY EQUATIONS : SVD- BASED APPROACH
To begin constructing the “differentiated” SVD-based method for
computing the filter sensitivities, we pay attention to the underlying
SVD-based filter and remark that it is formulated in the so-called
array form. This makes the modern KF algorithms better suited to
parallel implementation and to very large scale integration (VLSI)
implementation as mentioned in [15]. Each iteration of the SVDbased filter examined has the following pattern: given a pre-array A ∈
R(k+s)×s , compute the post-array SVD factors W ∈ R(k+s)×(k+s) ,
Σ ∈ R(k+s)×s and V ∈ Rs×s by means of the SVD factorization
 
S
T
, S = diag{σ1 , . . . , σs }
(18)
A=WΣV , Σ=
0
where the matrix A is of full column rank, i.e. rank A = s; the W,
V are orthogonal matrices and S is a diagonal matrix with singular
values of the pre-array A.
The goal of our study is to develop the method that naturally
extends formula (18) on the post-array SVD factors’ derivative computation. More precisely, the computational procedure is expected to
utilize the pre-array A and its derivative A′θ for reproducing the SVD
post-arrays {W, Σ, V} together with their derivatives {W′θ , Σ′θ , V′θ }.
To achieve our goal, we prove the result presented below. We also
bear in mind that the SVD post-array factor W is of no interest in the
presented SVD-based KF for performing the next step of the filter
recursion and, hence, the quantity W′θ is not required to be computed.
Lemma 1: Consider the SVD factorization in (18). Let entries of
the pre-array A(θ) be known differentiable functions of a scalar
parameter θ. We assume that σi (θ) 6= σj (θ), j 6= i, for all θ. Given

Having applied a right-multiplier V and a left-multiplier W to
equation (22), we have
h
i
h
i
WT A′θ V = WT W′θ Σ + Σ′θ + Σ (VT )′θ V .
(23)
In deriving the equation above we take into account the properties of
any orthogonal matrix Q, i.e. QQT = QT Q = I.
It is also easy to show that for any orthogonal matrix Q the product
Q′θ QT is a skew symmetric matrix. Indeed, by differentiating
′ both
sides of QQT = I with respect to θ, we get Q′θ QT + Q QT θ = 0,
T
or in the equivalent form Q′θ QT = − Q′θ QT . The latter implies
′
T
that the matrix Qθ Q is skew symmetric.
For the sake of simplicity we introduce the following notations:
Υ = WT W′θ and Λ = VT V′θ . As discussed above, the matrices
Υ ∈ R(k+s)×(k+s) and Λ ∈ Rs×s are skew symmetric, because W
and V are orthogonal matrices. Hence, we have ΛT = −Λ. Taking
into account this fact, we obtain the following partitioning of the
matrix form of equation (23):
"
# 

   ′  
WT A′θ V s×s
S
[Υ]s×s [Υ]s×k S
S
 T ′ 
+ θ −
Λ.
=
0
0
[Υ]k×s [Υ]k×k 0
W Aθ V k×s
From the equation above, we derive the formula for the main s × s
block of the matrix product WT A′θ V
h
i
= [Υ]s×s S + Sθ′ − SΛ.
(24)
WT A′θ V
s×s

Hence, the diagonal matrix Sθ′ obeys the equation
i
h
− [Υ]s×s S + SΛ.
Sθ′ = WT A′θ V
s×s

(25)

Now, let us discuss formula (25) in details. Recall the matrices Υ
and Λ are skew symmetric matrices and, hence, their diagonal entries
are equal to zero. The multiplication of any skew symmetric matrix
by a diagonal matrix does not change the matrix structure, i.e. the
diagonal entries of the matrix products [Υ]s×s S and SΛ are equal to
zero as well. Meanwhile, the matrix WT A′θ V s×s is a full matrix
and contains a diagonal part. Hence, from equation (25) we conclude
that diagonal matrix Sθ′ is, in fact, a diagonal part of the main s × s
block of the matrix product WT A′θ V. This completes the proof of
formulas in equation (19).
Finally, we need to compute V′θ where Λ = VT V′θ . Since V is an
orthogonal matrix, we obtain V′θ = VΛ. Next, any skew symmetric
matrix can be presented as a difference of a strictly lower triangular

PREPRINT

4

matrix and its transpose. Hence, the skew symmetric matrices [Υ]s×s
and Λ can be represented as follows:
[Υ]s×s = L̄T1 − L̄1

Λ = L̄T2 − L̄2

(26)

where L̄1 and L̄2 are strictly lower triangular
matrices.


Next, we split the matrix product WT A′θ V s×s into strictly
lower
diagonal and strictly upper triangular parts, i.e.
 T triangular,

W A′θ V s×s = L̄ + D + Ū . It was proved above that Sθ′ = D.
Taking into account this fact, the substitution of both formulas in (26)
into (25) yields
h
i
h
i
T
T
D
=
L̄
+
D
+
Ū
.
(27)
S
+
S
L̄
−
L̄
−
L̄
−
L̄
2
1
2
1
|{z} |
{z
}
| {z }
| {z }
′
′V
Sθ
T
W
A
[
θ ]s×s
Λ
[Υ]s×s

Hence, we obtain

L̄ + Ū =

[L̄T1

− L̄1 ]S −

S[L̄T2

− L̄2 ].

(28)

In (28), the L̄, L̄1 , L̄2 are strictly lower triangular matrices, the
Ū is a strictly upper triangular matrix and S is a diagonal. Hence,
equation (28) implies

Ū = L̄T1 S − S L̄T2 ,
L̄ = −L̄1 S + S L̄2 .

A LGORITHM 2. D IFFERENTIATED SVD- BASED KF(x̄0 , Π0 )
Initial Step (k = 0).
1
Ω = QΩ DΩ QTΩ and R = QR DR QTR ✄ SVD factorization
✄ SVD factorization
2
Π0 = QΠ0 DΠ0 QTΠ0
1/2
1/2
3
Set QP1|0 = QΠ0 , DP1|0 = DΠ0 and x̂1|0 = x̄0 .

7
8
9

10
11
12
13
14
15
16

It can be solved with respect to entries of L̄2 as follows:
(l̄2 )ij =

ūji sj + l̄ij si
, i = 2, . . . , s, j = 1, . . . , i − 1.
s2i − s2j

The formula above is exactly equation (21). Having computed the
entries (l̄2 )ij we can form the matrix Λ = L̄T2 − L̄2 in (26) and,
then, compute the derivative V′θ = VΛ. This completes the proof
of (20) and Lemma 1.
Remark 1: The assumption of singular values of A(θ) being
distinct for all values of parameter θ is necessary for avoiding the
division by zero in formula (21). In future, if possible, we will intend
for relaxing this restriction, which reduces the practical applicability
of the proposed method.
For readers’ convenience, Algorithm 1 provides a pseudocode for
the computational scheme derived in Lemma 1.
A LGORITHM 1. D IFFERENTIATED SVD(A, A′θ )
Input: A, A′θ
✄ Pre-array and its derivative
1 Apply SVD from (18) to the pre-array A. Save W, S, V.
2 Compute the matrix product WT A′θV.

3 Extract the main s × s block M = WT A′θ V s×s .
4 M = L̄ + D + Ū . ✄ Split into strictly lower triangular, diagonal
✄ and strictly upper triangular parts
5 Given L̄, Ū and
S, compute
the lower triangular L̄2 by (21).


6 Find V′θ = V L̄T2 − L̄2 .
T
7 Find Sθ′ = D. Hence, Σ′θ = [Sθ′ | 0] .
′
′
Output: Σ, V and Σθ , Vθ ✄ Post-arrays and their derivative
The theoretical result presented in Lemma 1 can be further applied
to the SVD factorization-based KF discussed in Section III. The
obtained computational scheme is summarized in Algorithm 2 and
shown in the form of a pseudocode. The new “differentiated” SVDbased KF extends the underlying SVD-based filter on the derivative
computation (with respect to unknown system parameters) for updating the corresponding filter sensitivities equations. The method can
be used for replacing the conventional “differentiated KF” approach
discussed in Section II by inherently more stable approach, which
is preferable for practical implementation. Finally, we would like
to remark that any “differentiated” filtering technique consists of
two parts: i) the underlying KF variant, and ii) its “differentiated”
extension used for the filter sensitivities computation.

1/2

1/2

Set ∂i QP1|0 = ∂i QΠ0 , ∂i DP1|0 = ∂i DΠ0 , ∂i x̂1|0 = 0.
Measurement Update: (k = 1, . . . , N ).
5
Build pre-array A in (12) and its derivatives ∂i A, i = 1, p.
6
[Σ,
n V, ∂i Σ, ∂i V]o← Differentiated SVD(A, ∂i A).
4

17

18
19
20
21
22
23

1/2

1/2

DRe,k , ∂i DRe,k ← read-off from Σ, ∂i Σ (i = 1, p).
QRe,k , ∂i QRe,k ← read-off from V, ∂i V (i = 1, p).
−1
Find K̄k from (13) and Kk = K̄k DR
QTRe,k .
e,k


∂i K̄k = ∂i QPk|k−1 DPk|k−1 QTPk|k−1 H T QRe,k , i = 1, p.


Build pre-array A in (14) and its derivatives ∂i A, i = 1, p.
[Σ,
n V, ∂i Σ, ∂i V]o← Differentiated SVD(A, ∂i A).
1/2
1/2
D
, ∂i DPk|k ← read-off from Σ, ∂i Σ (i = 1, p).
n Pk|k
o
QPk|k , ∂i QPk|k ← read-off from V, ∂i V (i = 1, p).

Find a posteriori
estimate
x̂k|k and ēk from (15).



T
∂i ēk = ∂i QRe,k zk − H x̂k|k−1


− QTRe,k (∂i H) x̂k|k−1 + H ∂i x̂k|k−1 , i = 1, p.
 −1
−1
∂i x̂k|k = ∂i x̂k|k−1 + ∂i K̄k DRe,k ēk + K̄k DR
(∂i ēk )
e,k
 −1
−1
− K̄k DRe,k ∂i DRe,k DRe,k ēk , i = 1, p.
Time Update: (k = 1, . . . , N ).
Build pre-array A in (16) and its derivatives ∂i A, i = 1, p.
[Σ,
Differentiated SVD(A, ∂i A).
n V, ∂i Σ, ∂i V] ← o
1/2

1/2

D
, ∂i DPk+1|k ← read-off from Σ, ∂i Σ (i = 1, p).
n Pk+1|k
o
QPk+1|k , ∂i QPk+1|k ← read-off from V, ∂i V (i = 1, p).

Find a priori estimate x̂k+1|k from equation
(17).

∂i x̂k+1|k = (∂i F ) x̂k|k + F ∂i x̂k|k + (∂i B) uk , i = 1, p.
End.

At the same manner, one can naturally augment any existing SVDbased KF variant (see, for instance, the algorithms in [23], [24]) or
potentially new SVD-based KF implementation on the corresponding
filter sensitivities computation.
Finally, taking into account the properties of orthogonal matrices,
it is not difficult to show that the negative log likelihood function
(LF) given as [33]:
N


o
1 Xn
−1
L θ, Z1N = c0 +
ln (det Re,k ) + eTk Re,k
ek
2 k=1

can be rewritten in terms of the SVD filter variables QRe,k , DRe,k
and ēk appeared in equations (12) – (17) as follows:
N
o



1 Xn
−1
ēk
ln det DRe,k + ēTk DR
(29)
L θ, Z1N = c0 +
e,k
2 k=1

where Z1N = {z1 , . . . , zN } is N -step measurement history and c0 is
ln(2π).
a constant value where c0 = Nm
2
Taking into account that the matrix DRe,k is diagonal and using the
Jacobi’s formula, d A−1 = −A−1 (dA) A−1 , from (29) we obtain
the expression for the log LF gradient evaluation in terms of the SVD
filter variables and their derivatives computed in the newly-developed
Algorithm 2 (for each i = 1, . . . , p):
N




 −1 i
1 Xn h
−1
∂i L θ, Z1N =
tr ∂i DRe,k DR
+2 ∂i ēTk DR
ēk
e,k
e,k
2
k=1
o
 −1
−1
ē .
−ēTk DR
∂i DRe,k DR
(30)
e,k k
e,k

PREPRINT

5

V. N UMERICAL EXPERIMENTS
By using simple test problem, we would like to demonstrate
thoroughly each step of the method summarized in Algorithm 1.
Example 1: Given pre-array A(θ) and its derivative A′θ


−2θ
sin(θ)
2

2θ
θ 


2
sin
(θ)
1/3
θ3 
A(θ) = 
,


θ 2θ2 − 1 
cos2 (θ) θ3 + θ2

compute the corresponding SVD post-arrays Σ, V and their derivative
Σ′θ , V′θ at the point θ̂ = 0.5.
Table 1 illustrates each step of the computational scheme in
Algorithm 1. To assess the accuracy of computations, we compute
′
′
l∞ = (AT A)θ̂=0.5 − (VΣ2 VT )θ̂=0.5
. This quantity should be
∞
small. Indeed, taking into account the properties of diagonal and
orthogonal matrices, from (18) we have AT A = VΣT WT WΣVT =
VΣ2 VT . Hence, the derivatives of both sides of the last formula
should coincide as well. In our numerical experiment we obtain
l∞ = 1.99 · 10−15 . This justifies the correctness of computations via
Algorithm 1 and confirms the theoretical derivations in Lemma 1.
TABLE I
A LGORITHM 1 ILLUSTRATIVE CALCULATIONS FOR E XAMPLE 1

Input

Pre-array derivative:

Line 1.

 −0.6070

W=
Σ=



0.8776
1.0000
0.2500 
2.0000
1.7500



" 1.7061

0#
 0.9967 0.0811 
0 0.8185
0
0 , V = −0.0811 0.9967
0
0
0  0

2.2959 −1.6522
1.1584 0.5691
M =  0.5177 −0.1427 .
−0.2470 −2.2944
−1.8181 −0.8517

Compute

Line 3.

Extract [M ]2×2 =

Line 4.

Split [M ]2×2 =

Line 5.

Compute L̄2 =

Line 6.

V′θ

Output

" −1.0000

0.4848 0.1556 0.2057 0.5745
0.5723 0.4035 0.0539 −0.5533 0.4478
0.1323 0.0735 0.9579 0.1059 −0.2197 
0.3159 −0.5593 0.0946 0.4337 0.6247
0.4321 0.5329 −0.2152 0.6724 −0.1756

Line 2.

Line 7.

0.4794 #
1.0000 0.2500
0.2298 0.0417
0.5000 −0.5000
0.7702 
0.3750
−2.0000
2.0000
A′θ θ̂=0.5 =  0.8415
1.0000
−0.8415

Pre-array: A|θ̂=0.5 =

θ̂=0.5

Σ′θ θ̂=0.5

=
=

Post-arrays:

Post-arrays’





 2.2959

1.1584

0
0
1.1584 0
0
0
0.8348 0

−1.6522
0.5691



  2.2959 0   0 −1.6522 
+
0
0.5691 + 0
0



 0.0677 −0.8321 
0.8321

" 2.2959
0
0
0
0

0.0677

0
0.5691
0
0
0

#

" 1.7061

#
0
0
0.8185
0
0
Σ|θ̂=0.5 =
0
0
0
h 0
i
0.9967 −0.0811
V|θ̂=0.5 = −0.0811 −0.9967
derivative: Σ′θ θ̂=0.5 and V′θ θ̂=0.5

(Lines 6,7)

Next, we wish to demonstrate how the novel method for the filter
sensitivities evaluation (Algorithm 2) works in practice. For that,
we consider the parameter estimation problem where the gradientbased optimization method is applied for finding the optimal value
of unknown system parameters. We test the conventional “differentiated” KF (Eqs (3) – (11) in Section II) and the previously
derived SR- and UD-based “differentiated” KF variants from [26]
and [27], respectively, against the new “differentiated” SVD-based

KF (Algorithm 2). As discussed in Section IV, all “differentiated”
methods consist of two parts and, hence, they compute the Log
LF and its gradient simultaneously. These values are utilized by a
gradient-based optimization method for maximizing the log LF with
respect to system parameters. Our library of codes is implemented in
MATLAB where we use the built-in optimization method fminunc.
Example 2: Consider a linearized version of the in-track motion
dynamic when a satellite travels in a circular orbit [34]:




0 0 0 0
1 1 0.5
0.5


0 1
1
1 
 xk−1 + wk−1 , Ω = 0 0 0 0 
xk = 


0 0
0 0 0 0
1
0
0 0 0 q1
0 0
0
0.606


 2 2

θ δ
0
1 1 1
1
xk + vk ,
R=
zk =
0
θ2 δ 2
1 1 1 1+δ
where q1 = 0.63 · 10−2 , x0 ∼ N (0, θ2 I4 ) and θ is the unknown
system parameter that needs to be estimated. In contrast to [34], we
wish to test both well-conditioned and ill-conditioned situations. For
that, following [17], we simulate the roundoff by parameter δ. It is
assumed to be δ 2 < ǫroundof f , but δ > ǫroundof f where ǫroundof f
denotes the unit roundoff error1 . When δ → ǫroundof f , i.e. the
machine precision limit, the problem above becomes ill-conditioned.
By varying the ill-conditioning parameter δ, we are able to explore
some numerical insights of each method assessed.
The numerical experiment is organized as follows. For each fixed
value of ill-conditioning parameter δ, the SSM in Example 2 is
simulated for θ∗ = 5 to generate N = 100 measurements. Next,
the unknown system parameter θ is estimated from the available
experimental data, Z1N = {z1 , . . . , zN }, by using gradient-based
adaptive KF techniques examined, i.e. by four “differentiated” KF
methods mentioned earlier in this Section. For a fair comparison,
each “differentiated” algorithm utilizes the same data Z1N and the
same initial value for the optimization method, θ̂(0) = 1. Next, the
obtained optimal estimate θ̂∗ is compared with the “true” value of
θ∗ = 5 for assessing the estimation quality of each method. We repeat
the experiment M = 100 times and calculate a posterior mean of the
estimate, the root mean squared error (RMSE) and the mean absolute
percentage error (MAPE) over 100 Monte Carlo runs.
Having carefully analyzed the obtained numerical results summarized in Table 2, we make a few important conclusions. First, all
“differentiated” KF variants work equally well when δ is about 10−1
and 10−2 , i.e. when the problem is not ill-conditioned. This confirms that all “differentiated” techniques are algebraically equivalent.
Second, among all methods examined, the conventional approach
(“differentiated” KF) shows the worst performance. It degrades faster
than any other algorithms when δ → ǫroundof f . Furthermore, the
line in Table 2 means that MATLAB can not even run the algorithm.
Third, we analyze the outcomes obtained by other methods tested
and observe that the UD- and SVD-based “differentiated” techniques
produce a better estimation quality than the SR-based counterpart.
This conclusion is reasonable if we recall that in this paper we do
not explore the filtering algorithms, but their differential form for
the KF sensitivities computation. Any existing “differentiated” SR1/2
based scheme requires the triangular matrix inversion Re,k that is
a square-root factor of the innovation covariance Re,k ; see Eq (6)
in [26]. In contrast, the UD- and SVD-based “differentiated” methods
involve the inversion of only diagonal matrix DRe,k ; see (30) and
Eq (8) in [27]. Finally, we observe that the new SVD-based approach
slightly outperforms the UD-based counterpart when δ → ǫroundof f .
1 Computer roundoff for floating-point arithmetic is 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

6

TABLE II
E FFECT OF ROUNDOFF ERRORS IN ILL - CONDITIONED TEST PROBLEMS IN E XAMPLE 2; EXACT θ ∗ = 5, 100 M ONTE C ARLO RUNS

δ
10−1
10−2
10−3
10−4
10−5
10−6
10−7
10−8
10−9
10−10

“differentiated” KF
Mean
RMSE
MAPE%
5.0046
4.9649
5.2764
8.8812
0.2803
-0.1315
−
−
−
−

0.2485
0.2784
0.7027
4.1440
8.0217
7.2403
−
−
−
−

3.8829
4.2892
9.7757
77.623
>100%
>100%
−
−
−
−

“differentiated” SR-based KF
Mean
RMSE
MAPE%

“differentiated” UD-based KF
Mean
RMSE
MAPE%

“differentiated” SVD-based KF
Mean
RMSE
MAPE%

5.0046
4.9649
5.0083
4.9879
4.9509
4.9310
4.9298
−
−
−

5.0046
4.9649
5.0083
4.9879
4.9508
4.9323
4.9268
5.0437
6.0119
6.7496

5.0046
4.9649
5.0083
4.9879
4.9509
5.0288
4.9249
5.0493
5.9738
6.7021

0.2485
0.2784
0.3555
0.3715
0.3352
1.0362
0.3658
−
−
−

3.8829
4.2883
5.7217
5.8595
5.6154
6.8368
5.8586
−
−
−

In summary, the previously derived UD- and the new SVDbased techniques provide the best estimation quality when solving
parameter estimation problem by the gradient-based adaptive filtering
methodology. This creates a strong background for their practical use.
In our ill-conditioned test example, the new SVD-based approach
even slightly outperforms the UD-based counterpart.
R EFERENCES
[1] R. K. Mehra, “Approaches to adaptive filtering,” IEEE Trans. Automat.
Contr., vol. 17, no. 5, pp. 693–698, Oct. 1972.
[2] G. Bastin and M. R. Gevers, “Stable adaptive observers for nonlinear
time-varying systems,” IEEE Trans. Automat. Contr., vol. 33, no. 7, pp.
650–658, Jul. 1988.
[3] R. K. Mehra, “Optimal input signals for parameter estimation in dynamic
systems – survey and new results,” IEEE Trans. Automat. Contr., vol.
AC-19, no. 6, pp. 753–768, Dec. 1974.
[4] N. K. Gupta and R. K. Mehra, “Computational aspects of maximum
likelihood estimation and reduction in sensitivity function calculations,”
IEEE Trans. Automat. Contr., vol. AC-19, no. 6, pp. 774–783, Dec. 1974.
[5] P. A. Zadrozny and S. Mittnik, “Kalman-filtering methods for computing
information matrices for time-invariant, periodic, and generally timevarying VARMA models and samples,” Computers & Mathematics with
Applications, vol. 28, no. 4, pp. 107–119, 1994.
[6] A. Klein and G. Mélard, “Computation of the Fisher information
matrix for time series models,” Journal of computational and applied
mathematics, vol. 64, no. 1, pp. 57–68, 1995.
[7] A. Klein, G. Mélard, and T. Zahaf, “Construction of the exact Fisher
information matrix of Gaussian time series models by means of matrix
differential rules,” Linear Algebra and its Applications, vol. 321, no. 1,
pp. 209–232, 2000.
[8] P. A. Zadrozny, “Analytic derivatives for estimation of linear dynamic
models,” Computers & Mathematics with Applications, vol. 18, no. 6,
pp. 539–553, 1989.
[9] A. Klein and H. Neudecker, “A direct derivation of the exact Fisher
information matrix of Gaussian vector state space models,” Linear
Algebra and its Applications, vol. 321, no. 1, pp. 233–238, 2000.
[10] M. Verhaegen and P. Van Dooren, “Numerical aspects of different
kalman filter implementations,” IEEE Trans. Automat. Contr., vol. AC31, no. 10, pp. 907–917, Oct. 1986.
[11] M. Verhaegen, “Round-off error propagation in four generallyapplicable, recursive, least-squares estimation schemes,” Automatica,
vol. 25, no. 3, pp. 437–444, 1989.
[12] P. G. Kaminski, A. E. Bryson, and S. F. Schmidt, “Discrete square-root
filtering: a survey of current techniques,” IEEE Trans. Automat. Contr.,
vol. AC-16, no. 6, pp. 727–735, Dec. 1971.
[13] G. J. Bierman, Factorization Methods For Discrete Sequential Estimation. New York: Academic Press, 1977.
[14] A. H. Sayed and T. Kailath, “Extended Chandrasekhar recursion,” IEEE
Trans. Automat. Contr., vol. AC-39, no. 3, pp. 619–622, Mar. 1994.
[15] P. Park and T. Kailath, “New square-root algorithms for Kalman
filtering,” IEEE Trans. Automat. Contr., vol. 40, no. 5, pp. 895–899,
May 1995.
[16] T. Kailath, A. H. Sayed, and B. Hassibi, Linear Estimation. New Jersey:
Prentice Hall, 2000.

0.2485
0.2784
0.3555
0.3715
0.3353
1.0333
0.3562
0.3757
1.2179
2.6030

3.8829
4.2883
5.7217
5.8596
5.6162
6.8265
5.6883
6.0712
20.762
49.405

0.2485
0.2784
0.3555
0.3715
0.3352
0.3138
0.3507
0.3790
1.1853
2.5286

3.8829
4.2883
5.7217
5.8597
5.6150
4.8826
5.5674
6.0946
20.106
49.252

[17] M. Grewal and A. Andrews, Kalman filtering: theory and practice. New
Jersey: Prentice Hall, 2001.
[18] G. H. Golub and C. F. Van Loan, Matrix computations. Baltimore,
Maryland: Johns Hopkins University Press, 1983.
[19] N. J. Higham, “Analysis of the Cholesky decomposition of a semidefinite matrix,” University of Manchester, Tech. Rep. MIMS EPrint:
2008.56, 1990.
[20] X. Zhang, W. Hu, Z. Zhao, Y. Wang, and Q. Wei, “SVD based Kalman
particle filter for robust visual tracking,” in Proceedings of the 19th
International Conference on Pattern Recognition. Tampa, FL, USA:
IEEE, 2008, pp. 1–4.
[21] O. Straka, J. Dunik, M. Simandl, and J. Havlik, “Aspects and comparison
of matrix decompositions in unscented Kalman filter,” in Proceedings of
the IEEE American Control Conference (ACC), 2013, pp. 3075–3080.
[22] Q. Zhang, X. Meng, S. Zhang, and Y. Wang, “Singular Value
Decomposition-based Robust Cubature Kalman Filtering for an Integrated GPS/SINS Navigation System,” Journal of Navigation, vol. 68,
no. 3, pp. 549–562, 2014.
[23] L. Wang, G. Libert, and P. Manneback, “Kalman Filter Algorithm
based on Singular Value Decomposition,” in Proceedings of the 31st
Conference on Decision and Control. Tuczon, AZ, USA: IEEE, 1992,
pp. 1224–1229.
[24] M. V. Kulikova and J. V. Tsyganova, “Improved discrete-time
Kalman filtering within singular value decomposition,” IET Control
Theory & Applications, 2016, (in progress). [Online]. Available:
https://arxiv.org/abs/1611.03686
[25] G. J. Bierman, M. R. Belzer, J. S. Vandergraft, and D. W. Porter,
“Maximum likelihood estimation using square root information filters,”
IEEE Trans. Automat. Contr., vol. 35, no. 12, pp. 1293–1298, Dec. 1990.
[26] M. V. Kulikova, “Likelihood gradient evaluation using square-root
covariance filters,” IEEE Trans. Automat. Contr., vol. 54, no. 3, pp. 646–
651, Mar. 2009.
[27] J. V. Tsyganova and M. V. Kulikova, “State sensitivity evaluation within
UD based array covariance filter,” IEEE Trans. Automat. Contr., vol. 58,
no. 11, pp. 2944–2950, Nov. 2013.
[28] M. V. Kulikova and A. Pacheco, “Kalman filter sensitivity evaluation
with orthogonal and J-orthogonal transformations,” IEEE Trans. Automat. Contr., vol. 58, no. 7, pp. 1798–1804, Jul. 2013.
[29] M. V. Kulikova and J. V. Tsyganova, “Constructing numerically stable
Kalman filter-based algorithms for gradient-based adaptive filtering,” International Journal of Adaptive Control and Signal Processing, vol. 29,
no. 11, pp. 1411–1426, 2015.
[30] ——, “A unified square-root approach for the score and Fisher information matrix computation in linear dynamic systems,” Mathematics and
Computers in Simulation, vol. 119, pp. 128–141, 2016.
[31] H. Neudecker, “Some theorems on matrix differentiation with special
reference to Kronecker matrix products,” Journal of the American
Statistical Association, vol. 64, no. 327, pp. 953–963, 1969.
[32] E. E. Tyrtyshnikov, A brief introduction to numerical analysis. Springer
Science & Business Media: Springer, 2012.
[33] F. C. Schweppe, “Evaluation of likelihood functions for Gaussian
signals,” IEEE Trans. Inf. Theory, vol. IT-11, pp. 61–70, Jan. 1965.
[34] H. E. Rauch, C. T. Striebel, and F. Tung, “Maximum likelihood estimates
of linear dynamic systems,” AIAA journal, vol. 3, no. 8, pp. 1445–1450,
1965.

