A unified approach to smoothing formulas

Automatica, Vol. 12, pp. 147-157. Pergamon Press, 1976. Printed in Great Britain

A Unified Approach to Smoothing Formulas*t LENNART

LJUNG:~ and THOMAS

Summary--Based on various approaches, several different solutions to the smoothing problem have been given. The relationships between these solutions are not immediate, although they solve the same problem. Making use of a certain framework from scattering theory, we derive two families of solutions, with equations evolving forwards and backwards in time, respectively. Within these families three major previous approaches are obtained as special cases, and their relationships are clarified. The set of solutions also contains as a fourth special case a (new) backwards analog of the innovations solution. The Mayne-Fraser two-filter formula belongs to the set of backwards solutions, and within this framework certain difficulties with its interpretation can be resolved.

• (t) = .~(t; ~-, PO'), 2(0) :~o(t, ~') = 2 0 ; ~-, O, O) ~,(~" I t, II, xo°) estimate of x(~-) based on the assumption P(~-) = II, ~0") = x, °

~,(~ I t) = ~,(~ I t, P(-r), :~('r)) Script letters refer to the extended systems. Subscript B refers to the backwards models.

I. Introduction TH~ L[~EArt least-squares smoothing problem has been very widely studied, see, e.g. the references in the recent surveys [1, 2], and several different formulas and approaches have been used. The various solutions must, of course, all be the same on abstract grounds since they solve the same problem, but somewhat surprisingly it is not always easy to establish this equivalence by direct means. For example, Lainlotis [3, 4] has obtained certain interesting identities by comparing three smoothing formulas--the innovations formula of Kailath and Frost [5], the 'forward-backward' two-filter formula of Mayne [6] and Fraser [7] and a so-called 'partition' formula of Lainiotis [3]. One such identity, in the standard Kalman filter notation defined in Section II, is

Nomenclature ~t s

the left-hand point of the considered time interval the right-hand point of the considered time interval a running time variable. Mainly used in equations that are naturally solved from s = ~" cr a running time variable. Mainly used in equations that are naturally solved from tr = t, i.e. backwards in time. tr runs forwards in time, so with the lefthand side of the ODEs written as - d/do, the fighthand side shows the stability properties of the algorithm as it is solved backwards P(r) true initial covariance at ~II assumed initial covariance at r P(t, r ; H) solution to (7) with P(~-, 1-; II) = II

¢o(~', t) = eb-l(~", t),

(1)

where COO',t) is an observability matrix of a 'forward' filter,

e(t, z) = P(t, r; P(r)) Po(t, ~) = e(t, "r; O) P,(r [ t; II) smoothed error covariance, given by, e.g. (12)

~o0", t) =

f'• SoT(S, z ) H T (s) H(s) So(S, "f) ds,

d So(S, T) = [F(s)-Po(s, ~) Hr(s) H(s)] S0(s, t), ds

with P(~-) = II P,(r I t) = Ps(r I t; e0")) eb( ~" I t) = P,0" I t, co) S(t, s; ~', H) transition matrix for F - P ( t , r ; I I ) H T H with S(s, s; ~', II) = I S(t, s) = S(t, s; ~', P(~')); So(t, ~') = S(t, ~'; ~', 0) uZ(~',s;t, II) transition matrix for F + G G T I I * - X ( s ) P~(s [ t, II) H T H where II*(s) obeys (80) with II*(~-) = II and with UZ(s, s; t, I-I) = I ~(~-) true initial estimate at ~" xl ° assumed initial estimate at ~- for the filter x, ° assumed initial estimate at ~- for the smoother 2(t; ~-, II, xl °) estimate of x(t), based on the assumption e(~-) II, 2(~-) x~° =

KAILATH:~

So(y, *) = I,

d po(s ' z) = F(s) Po(s, "r)+Po(s, "f) FT(s)--Po(s, "r) ds x HT(s) H(s) Po(s, "c) + G(s) GT(s), Po('C, T) = O, and Pb-a(cr, t) is the solution of a 'backwards' Riccati equation, d

--'dG Pb-l((7' t) = Pb-l((r, t) F((7)+ FT(o") Pb-l(cr, t)

=

-- Pb-l((7, t) G(cr) GT(cr) eb-l((r, t)

* Received 9 May 1975; revised 8 September 1975. The original version of this paper was not presented at any I F A C meeting. It was recommended for publication in revised form by associate editor B. D. O. Anderson. ~"This work was supported in part by the Air Force Office of Scientific Research, A F Systems Command, under contract A F 44-620-69-C-0101, the Joint Services Electronics Program under contract N-00014-67-A-01120044 and the Industrial Affiliates Program at the Information Systems Laboratory, Stanford University. Information Systems Laboratory, Department of Electrical Engineering, Stanford University, Stanford, CA 94305, U.S.A. Ljung is now with the Department of Automatic Control, Lurid Institute of Technology, S-220 07 Lurid, Sweden.

+HT(cr) H(cr),

Pb-l(t, t) = O.

No simple algebraic proof of this identity seems to be available in the estimation literature, and the reader can verify that the usual methods, e.g. showing that both sides obey the same linear differential equation and are equal at some point, would be somewhat laborious. The same remark applies to several other identities in [4]. Furthermore, it is natural to wonder about the particular boundary conditions Po(~-, ~-) = 0, Pb-X(t, t) = 0 and to ask if there is not a form of the identity (1) for different boundary conditions. Again this question seems to be difficult to resolve by current techniques. The quantity Pb-l('r, t) arises in the Mayne-Fraser formula [6, 7] and the subscript b is used because Pb(~', t) 147

148

LENNART LJUNG a n d THOMAS KAILATH

is interpreted as the error-covariance matrix of a 'backwards' Kalman filter. However, as discussed in detail in Section IX, this interpretation is not easy to establish and in fact turns out to be a somewhat accidental by-product of the 'end' condition Pb(t, t) = oO, a choice whose rationale in the literature is, in our opinion, not very satisfactory. And again, we can ask whether there is not a version of the Mayne-Fraser formula, cf. (16), (17), that would hold for finite Pb(t, t) ? In this paper we shall show how a certain algebraic structure drawn from scattering theory, ~t la Redheffer [8, 9], can be combined with certain stochastic interpretations (in fact, one used by Zachrisson [10] to reduce smoothing problems to filtering problems) to give a general framework in which not only can the particular questions raised above be satisfactorily answered, but several new results in smoothing, and also other problems [1 I-13], can be obtained. Outline of the paper. In Section II we define the problem and describe the three solutions in [3-7] mentioned above. In Section III we give a very brief statement of some relevant results from scattering theory--a calculus for a semi-group of operators with a certain composition rule, Redheffer's star product. Further discussion of the background and of the simple and useful physical significance of this calculus can be found in [8-9]. In Section IV we provide a stochastic interpretation of certain scattering matrices encountered in Section III via a method first used by Zachrisson [10] to reduce the smoothing problem to an extended filtering problem. In Section V we shall show how the results of Sections III and IV immediately suggest a unified approach to the three known formulas of Section II for the covariance matrix of the smoothed error. Furthermore, the results of Sections III and IV suggest the possibility of several new formulas, especially some that appear to be quite difficult to obtain, or even to guess, with current techniques. Thus, in Section IV, we present the backwards analog of the innovations smoothing formula, using certain results on backwards Markov models obtained in [11]. A general feature of scattering formulations is the ease with which they lend themselves to studying the effects of modified initial conditions, as is well known in circuit theory applications, see, e.g. [14-15]. In Section VII we exploit this fact to give families of forwards and backwards smoothing formulas, in which the basic results are first obtained for certain 'nice' initial conditions and then corrected to correspond to the true initial conditions. The results in Sections V-VII are all for the error-covariance matrices. In Sections VIII and IX we shall give the corresponding formulas for the smoothed and filtered estimates. An elegant and far-reaching identity arises naturally in this analysis and for completeness we give a direct algebraic proof of it in Appendix I. An overview of the results and some indications of further directions are given in the concluding Section X, which the reader might read directly after this Introduction. The paper unavoidably contains a large number of symbols and to aid the reader some notations and notational conventions are summarized in the nomenclature.

We shall usually suppress the time argument in F, G and H. We also assume a priori knowledge, which may, for example, be the result of previous filtering, of the mean value and variance of X(T), E x(~) = .fir), E (xCz)-~(r)) (x(r)-~(~)) T = e(r).

(3)

Given observed data ~,~ = {y(s), T ~
(4)

f* )t(~', t) = Jr* ~ T(s, T) HT[y(s)-- H .~(s)] ds,

(5)

where

and .~(s) is the Kalman filter estimate, obeying d ~(s) = [ F - P ( s ) H T HI ~(s)+P(s) H T y(s), .~(~') = ~?(T), (6) d p(s) = FP(s) + P(s) FT + GG T - P(s) H T HP(s), P(~-) = P(~-).

(7)

The smoothed error covariance is given by Ps(T I t) = P(r) -P('r) ¢(T, t) P(r)

(8)

with 0(r, t) =

T(S, "r) H a: H~(s, ~') ds.

(9)

~b(s, T) is the transition matrix for the Kalman filter, d ~ ( s , ~') = [ F - P ( s ) H T HI $(s, ~'), ~0", T) = I. (10) From these formulas, several others can also readily be deduced, e.g. those of Ranch et aL [16] and those of Meditch [17]. A discussion of such relationships can be found in [5]. II.3 Two-filter approach. [6, 7] and is given by

This solution is treated in

~(~ It) = Ps(~')[z0", t)+P(~) -1 ~0")]

(1 l)

e.0" [ t) = [Pb(T, t)-l+P(7")-l] -I,

(12)

with II. Some previous solutions to the smoothing problem ILl Problem formulation. The smoothing and filtering problems to be treated in this paper are the following: Consider a signal process z(.), corrupted by additive white noise v('), and with a state space description 5c(t) = F(t) x(t)+G(t) u(t), z(t) = H ( t ) x(t),

t>~'r, ) j

d

= - p b - 1 F - - F T p b - l - - H T H + p b -x GGTpb -1, Pb(t, 0 -1 = O, (13) d

(2)

P~(cr, 0 -1

z(cr, t) = -- [ F - GG T Pb(~r, t)-l] T Z(O, t)--HTy(a),

y(t) = z(t)+v(t).

Z(t,t) = 0.

where, with T denoting transpose and E expectation, E u(t) u(s) ~ = I ~(t--s), E u(t) v(s) T - O.

E v(t) v(s) T = I ~ ( t - s ) ,

(14)

Notice that z(~, t) can be written as z(~, t) = P~-I(a, t) 2b(o, t),

(15)

A u n i f i e d a p p r o a c h to s m o o t h i n g f o r m u l a s where

d gb(~, t) = (Fq-Pb H T H) gb(tr, t)--Pb(tr, t) H T y(tr), d~ gb(t, t) arbitrary,

(16)

and d -~ff eb(<7, t) ---- FPb(tr, t)-t-eb(cr, t) F T - GG T -I-Pb(<7,t)

XH T HPb(cr, t),

Pb(t, t) = cO.

(17)

We may note that (16) and (17) are usually regarded as defining a backwards Kalman filter, solved from ¢r = t. Note that the end condition in (17) can be treated by taking the limit of the solution Pb(Cr, t) as Pb(t, t) -->oo. Then (13) and (14) are the equivalent information filter version of this Kalman filter. While this gives a nice intuitive picture of the smoothing solution as the combination of a forward and a backward estimate, the usual backwards-filter interpretation of (16) leaves unsettled some difficult questions, as we shall explain in more detail later in Section IX. II.4 Partitioned approach. This approach [3, 4], gives a different form of the solution to the smoothing, as well as to the filtering problem, by first considering certain special solutions, conditioned with respect to x(~-). The results are &(,r I t) = P,(,r I t) [ao(,r, t)+P(,r) -1 g0)] P+(,r I t) = [~0(-r, t) + P ( z ) - q -x,

09)

A0(±r,t) -- f~ ¢o T (s, ,r) H r [ y ( s ) - H go(S, ,r)] ds,

III. Some simple facts from scattering theory Scattering theory is a way of analyzing how certain propagating entities, such as neutrons or electromagnetic waves, travel through a medium and interact with it. The analysis is based on certain operators called transmission and reflection operators and the interesting fact is that Riccati equations arise naturally in this analysis, see, for example, [8, 9]. Since Riccati equations also arise in leastsquares estimation problems, we have the possibility of interpreting certain estimation equations in terms of a scattering model. The advantages of such an identification will soon be made clear. We shall only need some simple facts from the theory, which can be developed in a self-contained way [8, 9]. Consider entities u+(z) and u_(s) incident upon an 'obstacle', say a layer extending from T to s of some medium, and let u_(r) and u+(s) be emergent entities. We assume there exist four operators a, r, p, o~such that

u_(r)

ro(,r, s)

then it is easy to see by algebraic manipulation or by tracing rays through the figure that

u-('r) ] = so(,r, t) u+('r) ] , u+(t) u_(t) where S0(-r, t) is obtained by the decomposition rule

t

T

$0(±r, t) = ['¢0 (s, .r) n r riCer(S, ,r) ds,

(21)

So(T, t) = S0(r, s)* SO(s, t)

Jr

d

u_(s)

u_(s)U+(t)] =so(s,t) u+(s)

(20)

where

~o(z, s)

The matrix of the four operators will be written S0(-r, s) and called the scattering matrix. If we now have another layer extending from s to t, so that as illustrated in Fig. 1

(18)

with

149

(30a)

with

go(S, ,r) = (F-Po H T H) go(S, ,r)+Po H T y(s), go(-r, .r) = 0,

(22)

c

d

=

d

~s Po(s, ,r) = FPo + Po FT + GGT-- Po H T HPo, Po(-r,-r) = O, (23) and

C

D

[ A(I-bC) -la

B+ Ab(I-Cb) -1 D ] d . (30b) C+ d C ( l - bC) -1 a ( l - Cb) -a D

where I is the identity operator. This is known as Redbeffer's 'star-product' composition rule.

d

d--s•°(s' z) = [F-Po(s, .r) H T HI Co(S, ,r), ¢0(,r, ,r) = t.

(24)

ro(s,i~~/ ) po(S,f)

Furthermore, g(t) = g0(t)+~o(t, -r) gs(-r ] t),

(25)

e(t) = Vo(t, ,r)+ Co(t, ,r) e+(z I t) CoT(t, ,r).

(26)

The significance of this solution is that estimates based on a certain assumption of initial conditions, namely P(,r) = 0, can be modified to those for the true initial conditions without having to reprocess the data. This result will in the present paper be extended to any assumed initial condition, not necessarily P(,r) = 0. Womble [18] has pointed out that (26) can be used to improve the numerical properties oftbe solution of the Riccati equation. By comparing (19) with (12), Lainiotis [3, 4], obtained the interesting equation 0o(,r, t) = Pb-X(,r, t)

(27)

mentioned in the Introduction. Other relationships follow similarly, e.g. [ref. 4; equation (20)]

¢(t, ,r) P(z) = Co(t, ,r) P+(,r I t).

r/A

FI~. 1. To illustrate propagation through two adjacent layers (shown apart for clarity). Now we may ignore the underlying physical problem of propagating entities and layers and concentrate on the fact that what we have is a semi-group of operators closed under a particular composition rule. We can develop a calculus of such operators by examining the effect of infinitesimal changes. Thus suppose that

So(s, s+As) = So(s, s)+M(s) As+o(As), where

(28)

We shall see later the natural family of identities that includes the above results and several others, Sections V-VII, and shall also describe simple direct proofs.

$

&(s, s) = X

and

M ( s ) = [ f(s ) g(s) ] h(s) e(s) '

(31)

150

LENNART LJUNG a n d THOMAS KAILATH

so that

which leads to

M(s) = lim So(s, s + A s ) - I Aa-*o

d do. So(°" t)

As

is the il~nitesirnal generator of S. Now note that

r ao(o., t) [f(o.)+g(o.) ro(o., t)]

a0(o., t)g(o.)

"1

|

So(r, s+ As) = So(r, s)*So(s, s+As)

(32)

= o(As) + So(s, r) + As

I

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

= | |

.............

L ~(r, s) h(s)ao(% s)

t)+ro(o., t)f(o.)

. . . . .

I'e(o.)+ro(o.> t)

+ ro(o., t) g(o') ro(o., t)

r [f(s) + Po(% s) h(s)] i g(s) +f(s) po(r, s) + po(r, s)] × | ......

h(o.)+e(o.) ro(o.,

x g(o.)] x %(o., t)

(33)

So(t, t) = L

oto(r,s) [e(s)+ h(s) Po(% s)] 1

J ,

(39)

Note that there are certain interesting 'dualities' between (34) and (39). However, our main point now is that l'[*So(o., t) does not obey the same backwards equation, (39), as So(o-, t). The proper equation can nevertheless be calculated as follows, where for later application, we also explicitly let 1"[ be a function of a. We can write

From this we readily obtain the differential equations

d

d-~So(r, s) =

S ( a - d o . , t; fl) = f l ( a - d a ) * s o ( a - d o ' , t)

[f(s) + po(r, s) h(s)]

= l"I(o.-do.)*So(a-do., o.)*So(o., t)

ao(r, s) g(s)+f(s) po(Z,s)

= n(o.-do.)*[l+M(o.) do.+o.(do.)]*

+ po(r, s) e(s)

x n-l(o'),]['l(o-),so(o-, t)

+ po(r, s) h(s) x ...............................................

= [I+ J / H O ) do.]*S(o., t; 1"]),

po(r, s)

~(r, s) h(s) ao(r, s)

where ..¢t'ri(o.) is defined implicitly by

ao(r, s) [e(s) + h(s) p0(r, s)]

So(r, r) = 1 =

l'l (0. - d a ) * [ / + M(o.) do. + o(do.)] ,1"I-1(o.)

(10) 0

= I+MdU(o.) do.+o(do.).

(34)

I

Note that the reflection coefficient Po obeys a Riccati equation. Furthermore, Reid [19], showed that (34) can itself be rewritten as a Riccati equation if the columns of So are interchanged. That is, if

~ o ( s , r ) = ( p°(r's)

(40)

, ................................

[The reader may readily verify that the inverse 1"i-a in the star-product algebra coincides with the 'usual' inverse.] It follows by comparing (38) and (40) that S(o., t; 17[)obeys a differential equation in o. that has the same structure as (39), where however the entries f, g, h and e of M(o.) are replaced by the corresponding entries of..Ct'U(o.),denoted by

"tli'[(o.) =

a°(r's) ) = S o ( r , s ) j '

o~(r, s) ro(r, s)

J

=

=

(o.) /

o

+

'

(35)

o o°]

rl(o.) =

II(o.)] II(o.) '

d H(s) = f(s) II(s) + II(s) e(s) + g(s) (36)

circuit theory, a particular advantage of the scattering formulation is the ease with which initial conditions can be changed. Thus note that

So(o'-do., t) = So(a-do', o.).so(a, t) (38)

(42a)

in which case the entries of d / ~ ( . ) are readily found to be

ftl(a) = - e(o.) - II-1 (0.) g, gB(o.) = h(o.), I en(o.) = - f ( o . ) - g H-l(o.), ha(o.) g(o.). )

(42b)

Let us also remark that w i t h J as in (15) and ~(t, 0., 1"i) & [l'][(a)*So(o., t ) ] J we can write the backwards equation as

(37)

obeys the same differential equations as So, but with S(r, r ; f][) = /"I[. Hence #~(s, ~-; 1"][) = S(r, s; fl),,¢ obeys (36) with . ~ ( r , r ; 1"I)= l"][J. This fact leads to some useful identities which will be discussed later in Section V. The solution So(r, t) is a very useful one because its 'backwards' derivative is much easier to find than that of S(r, t; 1~). Thus note that we can write

= [I+ M(o,)do. + o(dcr)]*So(~r, t)

[II(o.) H(o.)

where

Change of initial conditions. As it is well known in

S(r, s; 1"]) & 1].S0(r, s)

[ f~(o.) ga(o') ] hn(o.) eB(o.) "

It is important to note that the exact form of .-~'rl(') will depend critically upon our choice for 1"I. If 17[ = L then the backwards equation is just (39) and dt'rI = M. Another important case to be encountered below is when

+ [ g(S)o oO]+~°(s'r)[ h(s)o o°] ~°(s'r)' ~o(r, r) = J .

(41)

[0 e 0
x

x

[oo}[oo] 0 fB(o')

+

[°° 1

0 hB(a)

+ ~(t, cr; ill)

P(t, o.; 1"1), ~(t, t; 1"1) = l ' l ( t ) J .

0 gB(o-)

(42c)

This concludes our presentation of the simple facts we shall need from scattering theory. It is clear that the above

A unified approach to smoothing formulas results can be derived without any reference to starproducts or scattering theory, as done, for example, by Reid [19], but we have found the above language to be mnemonic, suggestive and simple.

151

Hence the framework of Section III gives us several possible ways of solving for P(t, 7.), which we shall examine in the next three sections.

V. Solving for the covariance matrix IV. The filtering and smoothing problem as an extended

filtering problem In this section we shall provide a stochastic interpretation of the matrix S(% t, 1"I) of the previous section, cf. (35-37). This interpretation will relate the scattering problem to the estimation problem and will be the vehicle for many later discussions. Considering x(7.) as an unknown parameter, we can use the approach of the extended Kalman filter to estimate it. This is in fact the approach suggested by Zachrisson in [10] to solve the smoothing problem. Let ~'(s, 7.) =

(*') x(7.)

' 7. <~s ~ t,

(43)

From Section III it follows that So(7.,t) obeys (34) (equals (36)) with S0(7., 7.) = I or (39) with So(t, t) = L In addition S(7., t) = ~'*So(7., t) obeys (34) with S(7., 7.) = ~ . Hence (52) immediately suggests three ways of calculating ~(t, 7.), which we shall find will yield the three smoothing formulas of Section II. There is also a natural fourth way of calculating ~(t, 7.) from (52) and this will lead to certain new formulas to be discussed in Section VI.

1. First approach (forward innovations). Find S(7, t) by solving (34) with initial condition S(% 7 ) = g~. This is perhaps the most natural approach, and it corresponds to solving the Riccati equation (47) with the true initial condition. The solution is (cf. (34))

where x(.) is defined by equations (2). Then we can write

~ ~(s, 7.) = ~ ~ ( s , 7.) + flu(s),

(44)

d ply(t, 7.) = [F_Pxl(t ' 7.) H T H] Ply(t, 7.), dt P12(7., 7.) = e(r),

y(s) = ~a~(s, 7.) + v(s),

(o)

where ~'=

0

0

, if=

0

d

, ~=

(H

(54)

Ptl(t, 7") = FPlt + P u FT + G G T - P l l HT HP11, Pl1(7., 7.) = e(7.),

0). (45)

(55)

d

The filtered estimate of .~(t, 7.), given {y(s), 7.
d ~(s, 7.) = [ ~ - - ~ ( s , 7.) ~ T ~ ] ~(s, 7.) + ~ ( s , 7.) ~ ds

~t P22(t, 7.) = -P~x(t, 7.) H T HPI~(t, 7.), P~dT., 7.) = P(7.),

pzl(t, 7.) = pl,T(t, 7.).

y(s), (46)

(56) (57)

Now compare (54)-(57) with (7) and (10) which yields

where ~(s, 7.) obeys the Riccati equation

d ~(s, 7.)

elx(t, 7.) = P(t),

(58)

Ply(t, 7.) = r~(t, 7.) P(T).

(59)

Hence from (56) = .~'.~(S, 7.)+.~(S, 7.).v#-T+ ~ T

~(S, 7.),~T . ~ .~(S, 7.)

(47)

P22(t, 7.) = P,(7. I t)

with initial conditions

= P(7.) -P(7.) (P(r)

~(7.,7.) =

~(7.)

'

•~(r, 7.) =

P(7.)

(48) .~(t) ] x,(~" I t)

so that we here have formulas for both the filtered estimate 0~(t) and the smoothed estimate £,(7.1t). Therefore we should also get formulas for the smoothing error by examining the Riccati equation (47). Denote

2. Second approach (partitioned formulas). Find S(7., t) by solving first for S0(7., t) using (34) with So(r, r) = I, and then forming S(7., t) = ~*So(7., t) using the star-product. This gives, with ( P.2°(t, 7.) P11°(t, 7.) ) So(% t) =

d p12°(t, 7.)

E

P22°(t, 7.) P2I°(t, 7.)

P2~(t, z)

P~(t,

7.)

Since •~(t, z) = E [~(t, 7.)-~'(t, 7.)] (~(t, 7.)-~'(t, 7.))T,

(50)

it follows immediately that

e ~ ( t , 7.) = P(t),

e ~ ( t , 7.) = e , ( r I t).

(51)

•~(t, ~r) = [~'*S0(7., t ) ] / ~ SO', t ) /

PllO H T HPlsO,

Ptt2(7., ~') = 0,

(62)

F(s) -- HT(S) H(s)

P22°(7., 7.) = 0.

(63)

NOW compare with (23) and (24), which yields P11°(t, 7.) = Po(t, 7.),

(64)

Pl~°(t, 7.) = ~o(t, 7.).

(65)

Next (63) gives

G(s) G(s)~ ) FT(s)

-dt Pzz°(t' T) = -P21°(t, 7.) H T HPla°(t, 7.),

(52)

as defined in the previous section, if we take

(

d plxo(t ' 7.) = Fpno+plx o FT + GGT dt

(61)

d

Furthermore, by comparing (47), (48) with formulas (35)(37) from the previous section on scattering theory, we see that we can write

M(s) =

P1~°(7., 7.) = /,

(49)

"

,

[F-- Plt°( t, 7.) H T H] P12°(t, T),

( P11(t, 7.) P12(t, 7.) ) ~(t, 7.) =

(60)

which is the expression (8) for the smoothed covariance naturally obtained in the innovations approach.

P(7.)

Note now that ~'(t, 7.) =

T(S, 7.) H T H ~(s, 7.) ds P(7.),

P(7.) )

(53) "

7.)

E

I~0T(S, "r) H T H~bo(S, 7.) ds = ~o(~', t), dr

(66)

152

LENNART LJUNG a n d THOMAS KAILATH

where the second equality is (21). Finally, from S(r,t)=

=

( S(t,'r)P(r)

P(t)

e,o" I t)

P(t) ST(t, r)

other backwards in time. However, we do not have the backwards counterpart of the forwards innovations formulas, e.g. (60). The main stumbling block to obtaining this result is the absence of results on backwards Markovian state models equivalent to a given forwards Markovian state model. However, insights from scattering theory have recently enabled us to obtain such a model, using which we can complete the results of this section.

)

~ * S o ( % t)

we obtain [cf. (30b)], P,(r [ t) = P(r)--P0-) 00(% t) [I+P(~-) tVo(~-, t)]-x P0") (67)

VI. The backwards counterpart of the innovations smoothing formula

where the second equation follows via the so-called ABCD lemma

The innovations formulas for the covariances (54)-(60) were obtained by solving for S(% t) = .~*S0(~-, t) using the forward equations (33). The backwards counterpart is naturally obtained by solving for this quantity using the backwards equations (39), with f, g, h and e given by the entries of ..¢t'ii(o), defined in (41). Before doing this, let us note that we can be somewhat more general and allow the 'boundary' condition ~' to be a function of the left end point. That is, we can consider

= [00(% t) +P0")-I] -a,

(A + BCD) -1 = A-x _ A-a B(C-a + DA-1 B)-XDA-a. (68) (Take A = P0-) -a, B = 00-I , C = D = L) Furthermore, P(t) = Po(t, r)+ So(t, r)P(T) (I+ ~o(% t)P0")) -1 SoT(t, ~-) = Po(t, ~') + So(t, ~') P,0" I t) SoT (t, ~)

(69)

and so that

S(t, T) P(T) = So(t, T) (I+ P(r) ¢0(% t)) -x P('r) = So(t, ~') e,(~" I t).

as o -->t, (70)

Clearly this solution corresponds to the partitioned approach described in Section II.4, and (67) is (19), (69) is (26), while (70) is (28). 3. Third approach (two-filter formula with P(t, t) = oo). The third approach is to solve for S0(T, t) backwards using (39) with So(t , t) = I (i.e. (42c) with 1"[ = I). This gives

as a -~ r,

rl(s) =

(71)

d dtr Pn°(t' o) = Pas°(t, o) GGT Pzt°(t, tr), Pn°(t, t) = 0, d do P~s°(t' o) = ( - H T H + F T P2s°+P2s ° F+Pss ° GGT Pss°), e22°(t, t) = O.

(73)

(74)

where the second equality is just the previously mentioned identity (27). Moreover, we recognize the state transition matrix for the filter (14) as P~°(t, r) = ~o(t, z). Finally, Pn°(t, z) = Po(t, r) =

P12°(t, s) GG T Psx°(t, s) ds. (75)

(YI(s) H(s)

II(s)) II(s) '

(79)

where II(r) = P(r).

(80) Then, it follows from (42) that the entries of .-~'ii(o) are given by fn = FnT(o) = -- ( F + GG T II-I(O))T,

gB = -- H T H, hB = GGT,

[ ]

(81)

Consequently, S ( o , t ) = l"l(o).S0(o, t) obeys (39), with fB, ga, h8 and eB given by (81) and with S(t, t) = l'l(t), cf. also (42c). More explicitly, we have d do Pl~(t' o) = - Pl~(t, ~) [F T + II-a(o) GGT + H T x HPaz(t, o)], Paz(t, t) = II(t), d - - - ~ Pea(t, o) = GGT - [F+ GG T II-l(o)] Pz~(t, o)

Using (74) in (67) clearly gives P,0"I t) = (Pb-l(r, t)+P(~-)-x) -1

17[(r) = ~ .

eB = - (F + GGT H-l(o)).

Next compare (73) with (13), which gives Pss°(t, r) = --Pb-l(% t)) ( = --d~o(t, ~')),

S(t, t) = I'l(t),

d II(s) = F II(s) + II(s) F T + GG T, (72)

(78)

Different choices of the function 1"][(o)will lead to different families of solutions. The extra freedom obtained by allowing n to depend on tr turns out to be essential in interpreting the results in terms of a backwards Markovian model for the process {y(s)} defined by (2). The choice that allows this interpretation is the following. Let f][(o) obey

d do Pxz°(t' 0") = P12°(t, o) [F+ GG T Pz~(t, o)], Pas°(t, t) = I,

S(o, t) = l"J[(o),S0(o, t)

(76)

which is the formula (12) used in the two-filter solution. Note incidentally that the two different expressions for P~2°(t, -r) give the identity, which we have used elsewhere [20] S0(t, "r) = ~0T(t, -r). (77) We see, therefore, that the three apparently different approaches described in Section II can be regarded as being induced by three different ways of solving for ~(t, ~'), all of which are natural from the scattering formulation. Moreover, the story does not end here, since there are more ways of solving for ~(t, ~), which do not correspond to any presently known smoothing formula. In particular, we see from the scattering framework that the partitioning and two-filter approaches are natural counterparts of each other, one forwards in time and the

(82)

-- Pa2(t, o) (F + GGT II-a(o)) T -P22(t, tr) H T HP2~(t , o),

P~z(t, t) = II(t), d - - ~ P u ( t , o) = P12(t, o) H T HP2t(t, o),

(83)

Plx(t, t) = II(t).

(84)

If we compare (82)-(84) with the equations defining a Kalman filter, say (54)-(56), we find that Psi(t, o) is the filtered error covariance and P~z(t, o) the closed loop transition matrix of the Kalman filter for the system - d"'~xB(o) = FR(o') xn(o') + G uB(a),

yn(o) = Hxn(o) + vB(o),

(85)

A unified approach to smoothing formulas evolving backwards in time from tr = t to tr = ~'. Now we can readily calculate the covariance function E yn(s)yB(t) = Rn(s, t), under the assumption

where d H*(s) = F I I * ( s ) + II*(s) F T + and

=

l"I*(t)

and E us(a) xsT(t) = 0,

a < t.

It turns out, see [11], that this covariance function is the same as that of the process {y(s)} defined by (2) It also follows that E xB(s)xB(t) T = E x(s)x(t) T. Hence (85) defines a state representation of y(.) such that the state is a Markov process in decreasing time: It is a 'backwards model' for y(.). Filtering for this backwards process corresponds to smoothing for the forward one, and hence the smoothed error covariance P,(r I t) = P~2(t, r) obeys a Riccati equation backwards in time. These interpretations of (82)--(84) in terms of the backwards model (85) clearly show that the approach is a backwards analog of the innovations of Section V, approach 1. VII. Families of smoothing formulas The four solutions described in Sections V and VI are quite natural from the scattering formulation of Section III. The partitioned approach, as well as the MayneFraser approach, correspond to assuming a certain initial condition, ~(T, ~-) = ,t, and then modifying the result to the true initial condition ~ ( % ~ - ) = ~ . However, in principle, there is no specific reason to choose just,,¢ as the assumed initial condition. Since the star product is a powerful tool for handling changes in initial conditions, we can easily extend the partitioned approach to any assumed initial condition n : First solve forwards for

S(T, t; 1"I) = r[*s0(% t)

(86)

using (34) ( = 36) with S(% r ; rl) = n for any rI, and then take SO', t) = SO', t; ~ ) = ~ , r [ - x * s ( % t; 1~),

Remark. I f n =

( II II

I I ) , then 1,i_l does not exist. H

II-X-P(r) - 1 ) 0

I

'

(88)

since, cf. (30b),

0

I

II II

--

P22(a, t) = Pb(a, t).

(94)

Consequently, solving (89)-(92) backwards with H = oo corresponds to the 'backward Kalman filter' (16), (17) while approach 3 gives the equivalent information filter (13), (14). To summarize this discussion, we have described two families of solutions to determine the filtered and smoothed error covariances as well as the transition matrices for the corresponding filter. A family of forwards solutions is defined by (86)-(88). The special case r][ = ~ corresponds to the innovations approach, while I"[ = I gives the partitioned formulas. The family of backwards solutions (89)-(93) has the special case II = II*(s) = o% which corresponds to the Mayne-Fraser solution in Kalman filter form. This can, maybe more conveniently, also be written using the equivalent information filter formulation by solving for Pb-1~ and P C 1 rather than for ~ and Pb, which corresponds to Section V, solution 3. It should also be noted that our approach in fact gives algebraic proofs of equalities like (67), (69), (70), (74) ( = 1) and (77). VIII. Formulas for the estimates, forwards solutions The solution to the filtering and smoothing problem is given by (46)-(48). We have in the previous section discussed various ways of determining the solution to (47), P(t, "r). They are essentially obtained by assuming a certain initial condition, solving either forwards or backwards in time, and modifying the result to the true initial condition. The obvious question now is whether we can do the same with the estimates. Suppose that we solve (46)-(48) with initial conditions:

XIo ) xsO '

,~(,r, ,r) = ~"~= ( ]]1 I~3), 1-/2 II~

(89) Denote the result by

rq,

~(t, ~-; 1"I,~r0) =

x?) ]

~'(r I t; fI, ~ ) l"

Then, it is readily seen from (46) that

~f(t, ~-; fI, ~ ) =

[F+ GGTH*(~)-x] T,

gn(a) = - H T H, hB(a) = GGT, en(~) =fnT(~),

(93)

P12®(a, t) is the transition matrix for the filter (16) and

Y~(T, T) = ~.o =

where S(% t; 1"I)is found by solving (39) ( = 42c) backwards with S(t, t, fl) = l"[*(t) (90) and fB(~) = (FB*(~)) T =

H*(t)) II*(t) "

S(r, t, oo) = ( Pl~®(t, T) Pn®(t, T) ), Pi2®(t, r) P~x®(t,"r)

P(T) P(r)

[ O II-1-P-a('r) ],S(r,t;l-[), I 1

(II*(t) II*(t)

Consequently, we have a forwards and a backwards family of solutions, each parametrized by n . Examining (89)-(93) with II = 0% we find by letting II(s) -~ co that, with

Similarly, we can obtain a family of backwards solutions as follows: we use the formula

S(% t) =

=

(87)

However, ~.1"I -~ can then be replaced by

(I

GGT, H*(r) = II (92)

E un(s) uzT(t) 8(s--t), E on(s) vnT(t) 8(s-t), E us(s) vnT(t)-O, E xn(t) xn(t) T = II(t) =

153

¢(t,

(91)t

I

nl) x,0 +

.;.,

nl) P(s, . ; n,) .T

s)

ds1 1

II2 f~(s, ~'; % HI) H T [y(s)-- HYc(s, ~'; H1, x I0)] ds+ x, 0

"["We may remind the reader of our convention of suppressing the time argument in F, G, H.

'

l

(95)

154

LENNART LJUNG a n d THOMAS KAILATH

where P(s, r ; lla) is the solution of the Riccati equation (7) with P(r, r; Hi) = II~ and q~(t, s; r, II1) is the transition matrix of [F-P(t, r ; IIt) H ~r H] with ~(s, s; r, Ha)

Z

Remark. As seen from (95) the solution depends only on II~ and 172 due to the special structure of the extended problem. From the estimate (95) it is possible to transform to the true initial conditions Y'0") and ~ , and to any other initial conditions, see [13]. The result is

'"(1, T)= ~ 0 P(T) I I ~ - I - - P ( T ) f j : ( S , T ) T H T ~ \

× H¢(s, r) ds (P('r) -- I]1) II2 -1 /

/ I\

(96)

IX. Formulas for the estimates, backwards solutions The fact that ~(t, ~-) can be solved for using backwards equations, i.e. differential equations in ~-, with initial condition in r = t, suggests that the estimates also can be obtained in this way. The two-filter formula described in Section II.3 clearly has such an interpretation. However, it is not immediately clear why the initial (i.e. final) condition has to be chosen to be infinity for the backward Kalman filter. Neither is it obvious how the fact that u(s) and x(t), s ~ t are correlated in the state model (2) can be overlooked. In fact, x(s), s<~t produced by (2) backwards in time is no Markovian representation. The solution to these problems lies in the backwards model discussed in Section VI. The model (85) where lI(s) obeys (80) gives an XB(a) that is a Markov representation backwards in time, and such that xB(t) and uz~(cr) are uncorrelated for o < t. This is true for any initial condition on H0") in (80). Consequently, a backwards Kalman filter can be constructed for (85) and the smoothing problem can be transformed into a filtering problem. To show the symmetry between forward and backward time, we can differentiate the extended state (43) with respect to r and obtain

l

4 ~(t, a) = ~B(O) ~(t, o) + ~ UB(a), f

- x?) + (~0")- x, °) /

y(o)

where ¢(t, ~') = ¢(t, "r; ~-, P0"))-

where

The proof consists of straightforward calculations. The innovations formula (4)-(5) clearly is obtained with 17[ = ~

defined in (48),

~°(7) = [20") 20")] ~.

Considernowthechoicei~=l]o=

( 0I

01 )

"~B

and ~o

= 0. This gives flt= ~o(t, s)Po(s, r) H w

a~n.qt'(t,o) + VB(O),

(0

,~B(s) = =

o

)

0 -(F+GG ~c17-a(s)) (0

(99)

)

(0) '

fib =

G

H),

(100)

where II(s) obeys (80). This approach clearly generates a family of backwards solutions: Choose any n, take the corresponding backwards model (85) with FB = F~* as in (91) where II*(s) obeys (92) with II*(r) = II. Then construct the Kalman filter for this model

x y(s) ds ~'(t, z, r, flo, 0) A ~'0(t, r) ........................................

do;. d~ °2 B(t, a; H, 0) = (~~*(a) - ~( t, a)

j~t= ¢oT(S, ~') HT (y(s)

-- H ~0(s)) [ =

Y(t,t) = 0, (101) where ~(t, a) obeys the Riccati equation (cf. (42c)) (97)

where ~o(t, s) = ~(t, s; ~-, 0), Po(s,"r) = P(s, ~', O) and 2o(t) and AoO-,t) are defined by (22) and (20) respectively, and 1 q~(t, r) P0") .....................................................................

~!t':'~!"~!r"! -P('r)J~'= ~(s, r) T H T H~(s, +

~(r)

Using the identities (8) (or

(60))

x

~(r)

do

~) + P(t, o) ~'n*(a) T + fin ~B 'j! (102a)

- ~ ( t , a) ~ B T ~ B ~(t, or) with

II*(t) II*(t) ] II*(t)

II*(t)

"

(102b)

[Equation (102) is, of course, (39) with (90), (91).] Introduce

x ds P0")

x~C0(t,~')+

- --d ~(t, a) = ~'B*(O) ~(t,

d~(t, t) =

:-:(,,

]

II, O)+ ~(t, o) .~B T y(a),

ds

2o(t) ] A0(r, t) '

o

a~etOTB .9~B) ~ B(t, cr;

( ,'~BS(t[ O; II, 0) ) ds

• (98)

~n(t, a; 17, 0) =

~n(O, t; H, 0)

"

Then the estimate at a = r can be written II*(t) f:~T(S, t; t, II) HT(y(s)

--P(r)ff=:(s,,) T H T "Cfi(s,t)-t-I= Ps('r,t)P('r) -1,

--H 2B(s, t; II, 0)) ds

and (70),

~s(t, r ; II, O) = ¢(t, 0 co-) = ¢0(t, *) P,(T I t),

we see that the first block row of (98) is (25), while the second block row is (18). Consequently, this choice of initial conditions gives the partitioned solution to the filtering and smoothing problem.

.................................................... ~0", s; t, 1]) P,(s I t, II) H T x y(s) ds (103)

A unified approach to smoothing formulas where ~b0-, s; t, 1-I) is the transition matrix for

F+ GG T H * - l ( s ) - P , ( s [ t, H) H T H with

~(s,s;t, II) = I

and

P,(s [ t, II)

is the solution to the backwards Riccati equation (83) with P,(t I t, H) = II*(t). From stochastic arguments it is clear that £~B(t, 1-; II, 0) is the estimate ~(t, 0-; II, 0) from the forward model and forward Kalman filter given by (95) with II, = H a = H. We shall, however, give an independent algebraic proof of this fact.

We may note that the a priori estimate can naturally be thought of as the result of filtering all data prior to s = ~'. The fact that the a priori data are given at the 'left-hand' point % as well as the fact that the original model (2) is Markov in increasing time, gives rise to the slight asymmetry in the forwards and backwards solutions. Clearly, with a priori data in t, the analogy would be complete. Let us collect x0") and x(t) into a vector

~ ' ( t ' ¢ ) = ( x ( t")x O ' )

=

H

~(t,r) =

P,(s I t, II) H T y(s) ds

9~ (s, ~'; % I'i)H T ( y ( s ) - H 2(s, ¢;

(108)

We then have the following possibilities to find the estimate

Lemma 9.1. A useful identity. f ~ o ' , s; t, rl)

155

( ) 2s(r [ t)

:

1. Consider the forward model II, 0))

ds (104)

~(s,r)=

0 0

o

and

~"~
H*(t) fi'g:(s, t; t, II) H T (y(s)--H ~B(S, t; I I ,

y(s) = (H O) ~(s, "0 + v(s).

0)) ds

) (109)

(t, s; % II) P(s, % II) 1-1T y(s) ds.

=

(105)

Choose any initial condition (110)

The proof of this lemma is given in Appendix I. By rearranging terms it is readily seen that (104) is equivalent to ~b0", s; t,

H~

1-Ix

x, °

and determine the Kalman filter estimate for (109) [defined by (46)-(47) and (95)]

H) P,(s I t, H) = II[~T(s, ~'; % 11)

~0-, t; It, ~0).

_f~=, ~T(~, z; ~', II) H T x H~(tr, s; % H) de P(s, "r; II)] (106) and similarly for (105).

Then calculate ~(~', t) using (96). The corresponding covariances are calculated using (86), (87). This is the family of forward solutions discussed in Section VIII. 2. Choose an initial condition

Remarks. The identity (106) is in fact a variant of the Krein-Siegert-Bellman resolvent identity for the Fredholm equations, see, for example, [21]. By differentiation of (106), and its counterpart, certain 'Levinson-type' equations can be derived also for nonstationary processes. This subject is further discussed in [12]. The conclusion of this section is that ~'(t, ~', 1"I, 0) can be computed for any

n = II

II

'

~.0 = 0

(lll)

and determine II*(s) as the solution of dlI*(s) = FH*+II*FT+GGT; ds

II*0") = H.

(112)

F o r m the corresponding backwards model

>o using the backwards filter (101)-(103) and then be adjusted to the true initial conditions by formula (96). Clearly the two-fdter solution described in Section II.3 is obtained when II -~ oo. It is also possible to include an assumption about x/° = x, ° = x* ~ 0. However, since the backward model (85) is constructed for a zero-mean process {y(t)}, we must proceed as follows: Form the modified, zero-mean corrected output

yM(G)

--~

y(o)-- H(aF(tr, ~) x*,

(107)

where Se(o, ~') is the transition matrix for ~', and let it be the input to the backward Kalman filter. Then add x* to the obtained result in z. Another approach is, of course, to use y(tr) straightforwardly in the filter and add x* to the final result, using (96).

X. Conclusions We have been studying the problem of finding the best linear least-squares estimates of x0") and x(t), I- < t, based on a priori estimates ~1-) with error covariance P(r) at time ~" and on observations {y(s); ,r<~s<~t}, where y and x are related according to (2).

~(t, o) =

(00

o

- (F+ GGTII*-I(~))

+( ; )uz(tr),

"r<~o<~t,

y(o) = (H 0) ~r(t, o) + v(o)

113) and determine the Kalman filter estimate for (113) with initial conditions in

( II*(t) o = t,

~(t,t)=

II*(t) )

II*(t) II*(t) '

£~B(t,t)=O"

It is given by (101), (102) and (103). Then calculate £~(% t) from ~B(~', t, II, 0) using (96). The covariances are obtained from (89)-(93). This is the family of backwards solutions discussed in Section IX. Note that in this case, i.e. 1I = II 1 = II~ and ~f0 = 0, the second block row of (96) can be written, using (60) and

(76)

~8(r ] t) = P,(~" I t) {p(¢)-1 g(¢)+ [II-1 +pb-10., t)] X XB0", t; 1-I, 0)}.

156

LENNART LJUNG a n d THOMAS KAILATH

This expression shows more clearly how we have obtained a generalization of the Mayne-Fraser formula (11), (15). Within these two families we recognize as special cases the presently known innovations approach, the two-filter approach and the partitioned approach corresponding to specific choices of II and ~o. The relationships between these seemingly different approaches have been clarified, and our discussion, in fact, contains algebraic analysis of their connections, as well as independent derivations of each of them, using only the Kalman filter results. Moreover, we believe that our approach gives a deeper understanding of the backward formulas. The fact that the Mayne-Fraser backward filter seems to be based on the forward model is illusory and is directly connected with the assumption of infinite initial covariance. For any other initial condition the correlation between u(s) and x(t) has to be decoupled using the true backwards model. Of course, we should ask what the advantage is of such a wealth of solutions. Apart from the satisfaction of a general solution, we may point out the following practical aspects. In the first place, the choice between forwards and backwards solutions is usually dictated by the kind of smoothing problem. For the fixed point smoothing problem, .~(r I t) is estimated for increasing t. Then forwards solutions, i.e. recursions in t, are most natural. For the fixed interval smoothing problem, t is fixed and our point of interest, ~-, is varying. Then recursions in ~-, i.e. backwards solutions, are the most natural ones. The possibility of choice between different initial conditions may be useful to improve the numerical properties of the filters, see, e.g. [18]. Also, as discussed in [13], for time-invariant systems this choice of II can lead to fast algorithms like the Chandrasekhar algorithms [22], for any initial conditions, whereas, as pointed out in [22], the rank reduction basically takes place only for very special initial conditions. Clearly the discussion of this paper can also be made for discrete-time systems. The results are analogous, and it turns out that the star-product is a very useful tool for handling the more complex expressions for filtering and smoothing in discrete time [23]. The explicit discussion is deferred to a subsequent paper.

References [1] J. S. MEDITCH: A survey of data smoothing for linear and nonlinear dynamic systems. Automatica 9, 151162 (1973). [2] T. KAILATH: Supplement to "A survey of data smoothing". Automatica 11, 109--111 (1975). [3] D. G. LAINmTIS: Optimal linear smoothing; continuous data case. Int. J. Control 17, 921-930 (1973). [4] D. G. L~a~OTIS: Partitioned estimation algorithms, II: Linear estimation. Information Sci. 7 317-340 (1974). [5] T. KAILATHand P. FROST: An innovations approach to least-squares estimation, Part II: Linear smoothing in additive white noise. IEEE Trans. Autom. Control AC-13, 655-660 (1968). [6] D. Q. MAYBE: A solution of the smoothing problem for linear dynamic systems. Automatica 4, 73-92 (1966). [7] D. C. FRASER: A New Technique for Optimal Smoothing of Data. Sc.D. Dissertation, MIT, Cambridge, MA (1967). [8] R. REDHEFFERl On the relation of transmission-line theory to scattering and transfer. J. Math. Phys. XLI, 1-41 (1962). [9] R. REDHEFFER: Difference equations and functional equations in transmission-line theory. In Modern Mathematics for the Engineer, Chapter 12, Second series. (E. F. BECKm,mACH, ed.), McGraw-Hill, New York (1961). [10] L. ZACHmSSON:On optimal smoothing of continuoustime Kalman processes. Information Sci. 1, 143-172 (1969).

[11 a] L. LJUNG,T. KAILATHand B. FRIEDLANDERi Scattering theory and linear least-squares estimation, Part I: Continuous-time problems. Proc. IEEE64, (1), (1976). [llb] L. LJUNG and T. KAILATH: Backwards Markovian models for second-order stochastic processes. IEEE Trans. Inf. Theory (Submitted). [12] T. KAILATH and L. LJUNG: A scattering theory framework for fast least-squares algorithms. Proc. 4th lnt. Symposium on Multiv. Analyt., Dayton, Ohio, June (1975). North Holland (1976). [13] L. LJUNG and T, KAILATH: Formulas for efficient change of initial conditions in least-squares estimation. (Submitted for publication.) [14] H. J. CARLIN and A. B. GIORDANO: Network Theory. Prentice-Hall, Englewood Cliffs, N.J. (1964). [15] H. J. CARLIN: The scattering matrix in network theory. IRE Trans. PGCT CT-3, 88-97 (1956). (Special Issue on the Scattering Matrix.) [16] H. E. RAUCH, F. TUNG and C. T. STRIEaEL:Maximumlikelihood estimates of linear dynamic systems. AIAA J. 3, 1445-1450 (1965). [17] J. S. MEDITC8: On optimal linear smoothing theory. J. Information Control 10, 598-615 (1967). [18] M. E. WOMBLE: The Linear Quadratic Gaussian Problem with Ill-conditioned Riccati Matrices. Ph.D. Dissertation, MIT, Cambridge, MA (1972). [19] W. T. REID: Riccati Differential Equations. Academic Press, New York (1972). [20] T. KAILATHand L. LJUNG: Asymptotic behaviour of constant-coefficient Riccati differential equations. IEEE Trans. Autom. ControL (Submitted). [21] T. KAILATH:A view of three decades of linear filtering theory. IEEE Trans. Inform. Theory IT-20, 146-181 (1974). [22] T. KAILATH: Some new algorithms for recursive estimation in constant linear systems. 1EEE Trans. lnform. Theory IT-19, 750-760 (1973). [23] B. FRIEDLANDER,T. KAILATHand L. LJUNG: Scattering theory and linear least-squares estimation, Part II: Discrete-time problems. J. Franklin Inst. 301, 71-82 (1976).

Appendix L An algebraic proof of lemma 9.1 Consider (104). Let, for convenience, the argument 1-1 be suppressed. First, we have for the second term in the right-hand side t

T

Js q[ (S, "1";"r) H T H Yc(s,"r) ds = fs=r f:=r ¢ T ( s ' ' r ; ~ ' ) H T H ¢ ( s ' a ; ~ ' ) P ( a ) H T y ( a ) d a da.

J~

Hence the coefficient for HTy(s) in the integrand of the right-hand side of (104) is R0", s, t) ~ II

T(s, "r; r ) -

t~T(tr,r ; ~') H T H =8

x ~(a, s; "r) do P(s, "r)].

(114)

Equation (104) is therefore shown if (114) is equal to

L(z, s, t) ~- ~(7, s, t ) P , ( s l t ) .

(115)

To show this take first t = s. This gives R(~', t, t) = II~bT(t, ~'; q')

(116)

L0", t, t) = ~b0",t; t) H*(t).

(117)

and But (116) is the (2, 2) element of S(~', t; l'I) = 1]*(T)*So(r, t) solved for forwards, while (117) is the expression for the same element solved backwards. Hence R(~-, t, t) = L(~', t, t).

(118)

A unified approach to smoothing formulas Moreover, d d--t R(T, s, t) = - I I ~T(t, T; T) H T H ~(t, s; ¢)P(s, ¢).

(119) To find (d/dt) L(¢, s, t) we first note that ~('r, s; t) = ~('r, t; t) q~-l(s, t; t)

157

which is a well-known formula for changing initial conditions in Riccati equations, and hence d ~-~L('r, s, t) -- - II ~T(t, "r; "r) H T H ~(t, s; T) P(s, m')

x {[II*-'(s)-e-'(s,

~')] e,(s I t)+P-'(s, ~)

x ~-l(t, s; ~) ~(t, s, s) II*(s)}.

= II~T(t, r; ~') ~-~(t, s; s) (n*(s))-~. Hence, using the forwards expressions (33) for (d/dt) $(t, m-) and (d/dt)Ps(¢ [ t)

It now only remains to show that the expression within brackets is identically I. This is readily established by differentiation with respect to t, and making use of (120). Hence

d'tdL(T, s, t) = II{~T(t, ¢; m') [ F T - H T He(t, ~)] ~-T X (t, S; S) II*-l(s)-- ~T(t, "r; r) x [ F T - H T HP(t, s)] $-~r(t, s; s) II*-l(s)} xP,(s [ t)-- II~T(t, r ; T) ~T(t, s; s) II-*(s) x II*(s)~T(t, s; s) n T n q~(t, s, s)II*(s) = II~T(t, ~'; ~) H T H{[P(t, s ) - e ( t , m')]

~tLO ", s, t) =

d

R(r, s, t)

and since L(¢, t, t) = R0", t, t), (104) follows. Equation (105) is, of course, analogous. This algebraic proof is rather long, and the reader may find it more profitable to study the 'obvious' stochastic proof noted in Section IX. The identities (104)-(106) turn out to be useful in many calculations.

x ~-T(t, s; s)II-*(s) P,(s [ t)+q~(t, s; s) x

II*(s)}.

Now P(t, s)--P(t, T) = --~(t, s; r) [P(s, r ) - - II*(s)] ~T(t, s; s) (120)

Note added in proof: After the submission of this paper, we have been aware that results, similar to those of Section V, have been derived also by G. S. Sidhu and U. B. Desai in 'A method for orthogonal directions'. (Submitted to S I A M J. Control.)