R. E. KALMAN
Research Institute for Advanced Study,2
Baltimore, Md.
Introduction
A New Approach to Linear Filtering
and Prediction Problems1
The classical filtering and prediction problem is re-examined using the Bode-
Shannon representation of random processes and the “state transition” method of
analysis of dynamic systems. New results are:
(1) The formulation and methods of solution of the problem apply without modifica-
tion to stationary and nonstationary statistics and to growing-memory and infinite-
memory filters.
(2) A nonlinear difference (or differential) equation is derived for the covariance
matrix of the optimal estimation error. From the solution of this equation the co-
efficients of the difference (or differential) equation of the optimal linear filter are ob-
tained without further calculations.
(3) The filtering problem is shown to be the dual of the noise-free regulator problem.
The new method developed here is applied to two well-known problems, confirming
and extending earlier results.
The discussion is largely self-contained and proceeds from first principles; basic
concepts of the theory of random processes are reviewed in the Appendix.
AN
IMPORTANT class of
theoretical and practical
problems in communication and control is of a statistical nature.
Such problems are: (i) Prediction of random signals; (ii) separa-
tion of random signals from random noise; (iii) detection of
signals of known form (pulses, sinusoids) in the presence of
random noise.
In his pioneering work, Wiener [1]3 showed that problems (i)
and (ii) lead to the so-called Wiener-Hopf integral equation; he
also gave a method (spectral factorization) for the solution of this
integral equation in the practically important special case of
stationary statistics and rational spectra.
Many extensions and generalizations followed Wiener’s basic
work. Zadeh and Ragazzini solved the finite-memory case [2].
Concurrently and independently of Bode and Shannon [3], they
also gave a simplified method [2] of solution. Booton discussed
the nonstationary Wiener-Hopf equation [4]. These results are
now in standard texts [5-6]. A somewhat different approach along
these main lines has been given recently by Darlington [7]. For
extensions to sampled signals, see, e.g., Franklin [8], Lees [9].
Another approach based on the eigenfunctions of the Wiener-
Hopf equation (which applies also to nonstationary problems
whereas the preceding methods in general don’t), has been
pioneered by Davis [10] and applied by many others, e.g.,
Shinbrot [11], Blum [12], Pugachev [13], Solodovnikov [14].
In all these works, the objective is to obtain the specification of
a linear dynamic system (Wiener filter) which accomplishes the
prediction, separation, or detection of a random signal.4
———
1 This research was supported in part by the U. S. Air Force Office of
Scientific Research under Contract AF 49 (638)-382.
2 7212 Bellona Ave.
3 Numbers in brackets designate References at end of paper.
4 Of course, in general these tasks may be done better by nonlinear
filters. At present, however, little or nothing is known about how to obtain
(both theoretically and practically) these nonlinear filters.
Contributed by the Instruments and Regulators Division and presented
at the Instruments and Regulators Conference, March 29– Apri1 2, 1959,
of THE AMERICAN SOCIETY OF MECHANICAL ENGINEERS.
NOTE: Statements and opinions advanced in papers are to be understood
as individual expressions of their authors and not those of the Society.
Manuscript received at ASME Headquarters, February 24, 1959. Paper
No. 59—IRD-11.
Present methods for solving the Wiener problem are subject to
a number of limitations which seriously curtail their practical
usefulness:
(1) The optimal filter is specified by its impulse response. It is
not a simple task to synthesize the filter from such data.
(2) Numerical determination of the optimal impulse response is
often quite involved and poorly suited to machine computation.
The situation gets rapidly worse with increasing complexity of
the problem.
(3) Important generalizations (e.g., growing-memory filters,
nonstationary prediction) require new derivations, frequently of
considerable difficulty to the nonspecialist.
(4) The mathematics of the derivations are not transparent.
Fundamental assumptions and their consequences tend to be
obscured.
This paper introduces a new look at this whole assemblage of
problems, sidestepping the difficulties just mentioned. The
following are the highlights of the paper:
(5) Optimal Estimates and Orthogonal Projections. The
Wiener problem is approached from the point of view of condi-
tional distributions and expectations. In this way, basic facts of
the Wiener theory are quickly obtained; the scope of the results
and the fundamental assumptions appear clearly. It is seen that all
statistical calculations and results are based on first and second
order averages; no other statistical data are needed. Thus
difficulty (4) is eliminated. This method is well known in
probability theory (see pp. 75–78 and 148–155 of Doob [15] and
pp. 455–464 of Loève [16]) but has not yet been used extensively
in engineering.
(6) Models for Random Processes. Following, in particular,
Bode and Shannon [3], arbitrary random signals are represented
(up to second order average statistical properties) as the output of
a linear dynamic system excited by independent or uncorrelated
random signals (“white noise”). This is a standard trick in the
engineering applications of the Wiener theory [2–7]. The
approach taken here differs from the conventional one only in the
way in which linear dynamic systems are described. We shall
emphasize the concepts of state and state transition; in other
words, linear systems will be specified by systems of first-order
difference (or differential) equations. This point of view is
Transactions of the ASME–Journal of Basic Engineering, 82 (Series D): 35-45. Copyright © 1960 by ASME
natural and also necessary in order to take advantage of the
simplifications mentioned under (5).
(7) Solution of the Wiener Problem. With the state-transition
method, a single derivation covers a large variety of problems:
growing and infinite memory filters, stationary and nonstationary
statistics, etc.; difficulty (3) disappears. Having guessed the
“state” of the estimation (i.e., filtering or prediction) problem
correctly, one is led to a nonlinear difference (or differential)
equation for the covariance matrix of the optimal estimation error.
This is vaguely analogous to the Wiener-Hopf equation. Solution
of the equation for the covariance matrix starts at the time t0 when
the first observation is taken; at each later time t the solution of
the equation represents the covariance of the optimal prediction
error given observations in the interval (t0, t). From the covariance
matrix at time t we obtain at once, without further calculations,
the coefficients (in general, time-varying) characterizing the
optimal linear filter.
(8) The Dual Problem. The new formulation of the Wiener
problem brings it into contact with the growing new theory of
control systems based on the “state” point of view [17–24]. It
turns out, surprisingly, that the Wiener problem is the dual of the
noise-free optimal regulator problem, which has been solved
previously by the author, using the state-transition method to great
advantage [18, 23, 24]. The mathematical background of the two
problems is identical—this has been suspected all along, but until
now the analogies have never been made explicit.
(9) Applications. The power of the new method is most ap-
parent in theoretical investigations and in numerical answers to
complex practical problems. In the latter case, it is best to resort to
machine computation. Examples of this type will be discussed
later. To provide some feel for applications, two standard
examples from nonstationary prediction are included; in these
cases the solution of the nonlinear difference equation mentioned
under (7) above can be obtained even in closed form.
For easy reference, the main results are displayed in the form of
theorems. Only Theorems 3 and 4 are original. The next section
and the Appendix serve mainly to review well-known material in
a form suitable for the present purposes.
Notation Conventions
Throughout the paper, we shall deal mainly with discrete (or
sampled) dynamic systems; in other words, signals will be ob-
served at equally spaced points in time (sampling instants). By
suitable choice of the time scale, the constant intervals between
successive sampling instants (sampling periods) may be chosen as
unity. Thus variables referring to time, such as t, t0, τ, T will
always be integers. The restriction to discrete dynamic systems is
not at all essential (at least from the engineering point of view);
by using the discreteness, however, we can keep the mathematics
rigorous and yet elementary. Vectors will be denoted by small
bold-face letters: a, b, ..., u, x, y, ... A vector or more precisely an
n-vector is a set of n numbers x1, ... xn; the xi are the co-ordinates
or components of the vector x.
Matrices will be denoted by capital bold-face letters: A, B, Q,
Φ, Ψ, …; they are m × n arrays of elements aij, bij, qij,... The
transpose (interchanging rows and columns) of a matrix will be
denoted by the prime. In manipulating formulas, it will be
convenient to regard a vector as a matrix with a single column.
Using the conventional definition of matrix multiplication, we
write the scalar product of two n-vectors x, y as
n
x'y = ∑
=
1
i
iyx
i
= y'x
The scalar product is clearly a scalar, i.e., not a vector, quantity.
Similarly, the quadratic form associated with the n × n matrix Q
is,
n
x'Qx = ∑
xqx
i
ij
j
=
1,
ji
We define the expression xy' where x' is an m-vector and y is an
n-vector to be the m × n matrix with elements xiyj.
We write E(x) = Ex for the expected value of the random vec-
tor x (see Appendix). It is usually convenient to omit the brackets
after E. This does not result in confusion in simple cases since
constants and the operator E commute. Thus Exy' = matrix with
elements E(xiyj); ExEy' = matrix with elements E(xi)E(yj).
For ease of reference, a list of the principal symbols used is
given below.
Optimal Estimates
t
t0
x1(t), x2(t)
y(t)
x1*(t1|t)
L
ε
time in general, present time.
time at which observations start.
basic random variables.
observed random variable.
optimal estimate of x1(t1) given y(t0), …, y(t).
loss function (non random function of its argument).
estimation error (random variable).
Orthogonal Projections
Y(t)
linear manifold generated by the random variables
y(t0), …, y(t).
x (t1|t) orthogonal projection of x(t1) on Y(t).
x~ (t1|t)
component of x(t1) orthogonal to Y(t).
Models for Random Processes
Φ(t + 1; t)
Q(t)
transition matrix
covariance of random excitation
Solution of the Wiener Problem
x(t) basic random variable.
y(t) observed random variable.
Y(t)
Z(t) linear manifold generated by y~ (t|t – 1).
x*(t1|t) optimal estimate of x(t1) given Y(t).
x~ (t1|t)
linear manifold generated by y(t0), …, y(t).
error in optimal estimate of x(t1) given Y(t).
Optimal Estimates
To have a concrete description or the type of problems to be
studied, consider the following situation. We are given signal
x1(t) and noise x2(t). Only the sum y(t) = x1(t) + x2(t) can be ob-
served. Suppose we have observed and know exactly the values
of y(t0), ..., y(t). What can we infer from this knowledge in regard
to the (unobservable) value of the signal at t = t1, where t1 may be
less than, equal to, or greater than t? If t1 < t, this is a data-
smoothing (interpolation) problem. If t1 = t, this is called
filtering. If t1 > t, we have a prediction problem. Since our treat-
ment will be general enough to include these and similar
problems, we shall use hereafter the collective term estimation.
As was pointed out by Wiener [1], the natural setting of the
estimation problem belongs to the realm of probability theory and
statistics. Thus signal, noise, and their sum will be random
variables, and consequently they may be regarded as random
processes. From the probabilistic description of the random
processes we can determine the probability with which a par-
ticular sample of the signal and noise will occur. For any given
set of measured values η(t0), ..., η(t) of the random variable y(t)
one can then also determine, in principle, the probability of
simultaneous occurrence of various values ξ1(t) of the random
variable x1(t1). This is the conditional probability distribution
function
Transactions of the ASME–Journal of Basic Engineering, 82 (Series D): 35-45. Copyright © 1960 by ASME
Pr[x1(t1) ≤ ξ1|y(t0) = η(t0), …, y(t) = η(t)] = F(ξ1) (1)
Evidently, F(ξ1) represents all the information which the meas-
urement of the random variables y(t0), ..., y(t) has conveyed about
the random variable x1(t1). Any statistical estimate of the random
variable x1(t1) will be some function of this distribution and
therefore a (nonrandom) function of the random variables y(t0), ...,
y(t). This statistical estimate is denoted by X1(t1|t), or by just X1(t1)
or X1 when the set of observed random variables or the time at
which the estimate is required are clear from context.
Suppose now that X1 is given as a fixed function of the random
variables y(t0), ..., y(t). Then X1 is itself a random variable and its
actual value is known whenever the actual values of y(t0), ..., y(t)
are known. In general, the actual value of X1(t1) will be different
from the (unknown) actual value of x1(t1). To arrive at a rational
way of determining X1, it is natural to assign a penalty or loss for
incorrect estimates. Clearly, the loss should be a (i) positive, (ii)
nondecreasing function of the estimation error ε = x1(t1) – X1(t1).
Thus we define a loss function by
L(0) = 0
L(ε2) ≥ L(ε1) ≥ 0 when ε2 ≥ ε1 ≥ 0 (2)
L(ε) = L(–ε)
Some common examples of loss functions are: L(ε) = aε2, aε4,
a|ε|, a[1 – exp(–ε2)], etc., where a is a positive constant.
One (but by no means the only) natural way of choosing the
random variable X1 is to require that this choice should minimize
the average loss or risk
E{L[x1(t1) – X1(t1)]} = E[E{L[x(t1) – X1(t1)]|y(t0), …, y(t)}] (3)
Since the first expectation on the right-hand side of (3) does not
depend on the choice of X1 but only on y(t0), ..., y(t), it is clear that
minimizing (3) is equivalent to minimizing
E{L[x1(t1) – X1(t1)]|y(t0), ..., y(t)} (4)
Under just slight additional assumptions, optimal estimates can be
characterized in a simple way.
Theorem 1. Assume that L is of type (2) and that the conditional
distribution function F(ξ) defined by (1) is:
(A) symmetric about the mean ξ :
F(ξ – ξ ) = 1 – F( ξ – ξ)
(B) convex for ξ ≤ ξ :
F(λξ1 + (1 – λ)ξ2) ≤ λF(ξ1) + (1 – λ)F(ξ2)
for all ξ1, ξ2 ≤ ξ and 0 ≤ λ ≤ 1
Then the random variable x1*(t1|t) which minimizes the average
loss (3) is the conditional expectation
x1*(t1|t) = E[x1(t1)|y(t0), …, y(t)] (5)
Proof: As pointed out recently by Sherman [25], this theorem
follows immediately from a well-known lemma in probability
theory.
Corollary. If the random processes {x1(t)}, {x2(t)}, and {y(t)}
are gaussian, Theorem 1 holds.
Proof: By Theorem 5, (A) (see Appendix), conditional distribu-
tions on a gaussian random process are gaussian. Hence the re-
quirements of Theorem 1 are always satisfied.
In the control system literature, this theorem appears some-
times in a form which is more restrictive in one way and more
general in another way:
Theorem l-a. If L(ε) = ε2, then Theorem 1 is true without as-
sumptions (A) and (B).
Proof: Expand the conditional expectation (4):
E[x1
2(t1)|y(t0), …, y(t)] – 2X1(t1)E[x1(t1)|y(t0), …, y(t)] + X1
2(t1)
and differentiate with respect to X1(t1). This is not a completely
rigorous argument; for a simple rigorous proof see Doob [15], pp.
77–78.
Remarks. (a) As far as the author is aware, it is not known what
is the most general class of random processes {x1(t)}, {x2(t)} for
which the conditional distribution function satisfies the re-
quirements of Theorem 1.
(b) Aside from the note of Sherman, Theorem 1 apparently has
never been stated explicitly in the control systems literature. In
fact, one finds many statements to the effect that loss functions of
the general type (2) cannot be conveniently handled mathe-
matically.
(c) In the sequel, we shall be dealing mainly with vector-
valued random variables. In that case, the estimation problem is
stated as: Given a vector-valued random process {x(t)} and ob-
served random variables y(t0), ..., y(t), where y(t) = Mx(t) (M
being a singular matrix; in other words, not all co-ordinates of
x(t) can be observed), find an estimate X(t1) which minimizes the
expected loss E[L(||x(t1) – X(t1)||)], || || being the norm of a
vector.
Theorem 1 remains true in the vector case also, provided we
re- quire that the conditional distribution function of the n co-
ordi- nates of the vector x(t1),
Pr[x1(t1) ≤ ξ1,…, xn(t1) ≤ ξn|y(t0), …, y(t)] = F(ξ1, …,ξn)
be symmetric with respect to the n variables ξ1 – ξ 1, …, ξn – ξ n
and convex in the region where all of these variables are
negative.
Orthogonal Projections
The explicit calculation of the optimal estimate as a function of
the observed variables is, in general, impossible. There is an
important exception: The processes {x1(t)}, {x2(t)} are gaussian.
On the other hand, if we attempt to get an optimal estimate
under the restriction L(ε) = ε2 and the additional requirement that
the estimate be a linear function of the observed random
variables, we get an estimate which is identical with the optimal
estimate in the gaussian case, without the assumption of linearity
or quadratic loss function. This shows that results obtainable by
linear estimation can be bettered by nonlinear estimation only
when (i) the random processes are nongaussian and even then (in
view of Theorem 5, (C)) only (ii) by considering at least third-
order probability distribution functions.
In the special cases just mentioned, the explicit solution of the
estimation problem is most easily understood with the help of a
geometric picture. This is the subject of the present section.
Consider the (real-valued) random variables y(t0), …, y(t). The
set of all linear combinations of these random variables with real
coefficients
)(
iya
i
(6)
t
∑
=
t
i
0
forms a vector space (linear manifold) which we denote by Y(t).
We regard, abstractly, any expression of the form (6) as “point”
or “vector” in Y(t); this use of the word “vector” should not be
confused, of course, with “vector-valued” random variables, etc.
Since we do not want to fix the value of t (i.e., the total number
of possible observations), Y(t) should be regarded as a finite-
dimensional subspace of the space of all possible observations.
Transactions of the ASME–Journal of Basic Engineering, 82 (Series D): 35-45. Copyright © 1960 by ASME
Given any two vectors u, v in Y(t) (i.e., random variables ex-
pressible in the form (6)), we say that u and v are orthogonal if
Euv = 0. Using the Schmidt orthogonalization procedure, as de-
scribed for instance by Doob [15], p. 151, or by Loève [16], p.
459, it is easy to select an orthonormal basis in Y(t). By this is
meant a set of vectors et0, …, et in Y(t) such that any vector in
Y(t) can be expressed as a unique linear combination of et0, …, et
and
Thus any vector x in Y(t). is given by
Eeiej = δi j = 1 if i = j
(i, j = t0, …, t) (7)
= 0 if i ≠ j
x
=
t
∑
=
t
i
0
iea
i
and so the coefficients ai can be immediately determined with the
aid of (7):
exE
j
=
E
t
∑
=
t
i
0
eea
i
i
=
j
t
∑
=
t
i
0
eEea
i
i
j
=
t
∑
=
t
i
0
a
i
=δ
ij
a
j
(8)
It follows further that any random variable x (not necessarily in
Y(t)) can be uniquely decomposed into two parts: a part x in Y(t)
and a part x~ orthogonal to Y(t) (i.e., orthogonal to every vector in
'Y(t)). In fact, we can write
=+=
x
~
x
x
t
∑
=
t
i
0
(
Exe
i
)
e
i
~
x
+
(9)
Thus x is uniquely determined by equation (9) and is obviously
a vector in Y(t). Therefore x~ is also uniquely determined; it
remains to check that it is orthogonal to Y(t):
−
(
xE
=
−
=
Exe
i
exE
i
~
exE
i
)
ex
i
Now the co-ordinates of x with respect to the basis et0, …, et, are
given either in the form
(as in (8)) or in the form Exei (as in
(9)). Since the co-ordinates are unique, Exei =
(i = t0, ..., t);
iexE~ = 0 and x~ is orthogonal to every base vector ei; and
hence
therefore to Y(t). We call x the orthogonal projection of x on
Y(t).
iexE
iexE
There is another way in which the orthogonal projection can be
characterized: x is that vector in Y(t) (i.e., that linear function of
the random variables y(t0), ..., y(t)) which minimizes the quad-
ratic loss function. In fact, if w is any other vector in Y(t), we
have
wxE
(
−
2
)
=
~(
xE
−+
wx
2
=
[(
xE
−
x
)
+
(
wx
−
2
)]
)
Since x~ is orthogonal to every vector in Y(t) and in particular
wx − we have
wxE
(
−
)
2
=
(
xE
2
−
x
)
+
wxE
(
−
)
2
≥
(
xE
2
−
x
)
(10)
This shows that, if w also minimizes the quadratic loss, we must
have
random
variables x and w are equal (except possibly for a set of events
whose probability is zero).
which means
− wxE
that
the
2 =
0
)
(
These results may be summarized as follows:
Theorem 2. Let {x(t)}, {y(t)} random processes with zero mean
(i.e., Ex(t) = Ey(t) = 0 for all t). We observe y(t0), …, y(t).
If either
(A) the random processes {x(t)}, {y(t)} are gaussian; or
(B) the optimal estimate is restricted to be a linear function of
the observed random variables and L(ε) = ε2;
to
then
x*(t1|t) = optimal estimate of x(t1) given y(t0), …, y(t)
= orthogonal projection x (t1|t) of x(t1) on Y(t). (11)
These results are well-known though not easily accessible in
the control systems literature. See Doob [15], pp. 75–78, or
Pugachev [26]. It is sometimes convenient to denote the
orthogonal projection by
(
tx
1
)|
t
≡
x
(*
t
1
)|
t
=
([ˆ
txE
1
)
|Y(t)]
The notation Eˆ is motivated by part (b) of the theorem: If the
stochastic processes in question are gaussian, then orthogonal
projection is actually identical with conditional expectation.
Proof. (A) This is a direct consequence of the remarks in con-
nection with (10).
(B) Since x(t), y(t) are random variables with zero mean, it is
clear from formula (9) that the orthogonal part x~ (t1|t) of x(t1)
with respect to the linear manifold Y(t) is also a random variable
with zero mean. Orthogonal random variables with zero mean are
uncorrelated; if they are also gaussian then (by Theorem 5 (B))
they are independent. Thus
0 = E x~ (t1|t)
= E[ x~ (t1|t)|y(t0), …, y(t)]
= E[x (t1) – x (t1|t)|y(t0), …, y(t)]
= E[x (t1)|y(t0), …, y(t)] – x (t1|t) = 0
Remarks. (d) A rigorous formulation of the contents of this
section as t → ∞ requires some elementary notions from the
theory of Hilbert space. See Doob [15] and Loève [16 ].
(e) The physical interpretation of Theorem 2 is largely a matter
of taste. If we are not worried about the assumption of gaus-
sianness, part (A) shows that the orthogonal projection is the op-
timal estimate for all reasonable loss functions. If we do worry
about gaussianness, even if we are resigned to consider only
linear estimates, we know that orthogonal projections are not the
optimal estimate for many reasonable loss functions. Since in
practice it is difficult to ascertain to what degree of approxima-
tion a random process of physical origin is gaussian, it is hard to
decide whether Theorem 2 has very broad or very limited sig-
nificance.
(f) Theorem 2 is immediately generalized for the case of
vector-valued random variables. In fact, we define the linear
manifold Y(t) generated by y(t0), ..., y(t) to be the set of all linear
combinations
t
m
∑∑
=
t
i
0
=
1
j
)(
iya
ij
j
of all m co-ordinates of each of the random vectors y(t0), …, y(t).
The rest of the story proceeds as before.
(g) Theorem 2 states in effect that the optimal estimate under
conditions (A) or (B) is a linear combination of all previous ob-
servations. In other words, the optimal estimate can be regarded
as the output of a linear filter, with the input being the actually
occurring values of the observable random variables; Theorem 2
gives a way of computing the impulse response of the optimal
filter. As pointed out before, knowledge of this impulse response
is not a complete solution of the problem; for this reason, no
explicit formulas for the calculation of the impulse response will
be given.
Models for Random Processes
In dealing with physical phenomena, it is not sufficient to give
an empirical description but one must have also some idea of the
underlying causes. Without being able to separate in some sense
causes and effects, i.e., without the assumption of causality, one
can hardly hope for useful results.
Transactions of the ASME–Journal of Basic Engineering, 82 (Series D): 35-45. Copyright © 1960 by ASME
It is a fairly generally accepted fact that primary macroscopic
sources of random phenomena are independent gaussian proc-
esses.5 A well-known example is the noise voltage produced in a
resistor due to thermal agitation. In most cases, observed random
phenomena are not describable by independent random variables.
The statistical dependence (correlation) between random signals
observed at different times is usually explained by the presence of
a dynamic system between the primary random source and the
observer. Thus a random function of time may be thought of as the
output of a dynamic system excited by an independent gaussian
random process.
An important property of gaussian random signals is that they
remain gaussian after passing through a linear system (Theorem 5
(A)). Assuming independent gaussian primary random sources, if
the observed random signal is also gaussian, we may assume that
the dynamic system between the observer and the primary source
is linear. This conclusion may be forced on us also because of
lack of detailed knowledge of the statistical properties of the
observed random signal: Given any random process with known
first and second-order averages, we can find a gaussian random
process with the same properties (Theorem 5 (C)). Thus gaussian
distributions and linear dynamics are natural, mutually plausible
assumptions particularly when the statistical data are scant.
How is a dynamic system (linear or nonlinear) described? The
fundamental concept is the notion of the state. By this is meant,
intuitively, some quantitative information (a set of numbers, a
function, etc.) which is the least amount of data one has to know
about the past behavior of the system in order to predict its future
behavior. The dynamics is then described in terms of state
transitions, i.e., one must specify how one state is transformed
into another as time passes.
A linear dynamic system may be described in general by the
vector differential equation
dx/dt = F(t)x + D(t)u(t)
y(t) = M(t)x(t)
and
(12)
u(t)
where x is an n-vector, the state of the system (the components xi
of x are called state variables); u(t) is an m-vector (m ≤ n)
representing the inputs to the system; F(t) and D(t) are n × n,
respectively, n × m matrices. If all coefficients of F(t), D(t), M(t)
are constants, we say that the dynamic system (12) is time-
invariant or stationary. Finally, y(t) is a p-vector denoting the
outputs of the system; M(t) is an n × p matrix; p ≤ n
∑
M(t)
.
The physical interpretation of (12) has been discussed in detail
elsewhere [18, 20, 23]. A look at the block diagram in Fig. 1 may
be helpful. This is not an ordinary but a matrix block diagram (as
revealed by the fat lines indicating signal flow). The integrator in
Fig 1. Matrix block diagram of the general linear continuous-dynamic
system
———
5 The probability distributions will be gaussian because macroscopic
random effects may be thought of as the superposition of very many
microscopic random effects; under very general conditions, such ag-
gregate effects tend to be gaussian, regardless of the statistical properties
of the microscopic effects. The assumption of independence in this context
is motivated by the fact that microscopic phenomena tend to take place
much more rapidly than macroscopic phenomena; thus primary random
sources would appear to be independent on a macroscopic time scale.
F (t)
D(t)
x(t)
x(t)
y(t)
Fig. 1 actually stands for n integrators such that the output of
each is a state variable; F(t) indicates how the outputs of the
integrators are fed back to the inputs of the integrators. Thus fij(t)
is the coefficient with which the output of the jth integrator is fed
back to the input of the ith integrator. It is not hard to relate this
formalism to more conventional methods of linear system
analysis.
If we assume that the system (12) is stationary and that u(t) is
constant during each sampling period, that is
u(t + τ) = u(t); 0 ≤ τ < 1, t = 0, 1, … (13)
then (12) can be readily transformed into the more convenient
discrete form.
x(t + 1) = Φ(1)x(t) + ∆(1)u(t); t = 0, 1, …
where [18, 20]
Φ(1) = exp F = ∑
∞
=0i
Fi/i! (F0 = unit matrix)
and
∆(1) = ( ∫ 1
0
exp Fτ dτ) D
u(t)
x(t)
y(t)
∆(t)
∑
M(t)
x(t + 1)
Φ (t + 1; t)
unit
delay
Fig 2. Matrix block diagram of the general linear discrete-dynamic
system
See Fig. 2. One could also express exp Fτ in closed form using
Laplace transform methods [18, 20, 22, 24]. If u(t) satisfies (13)
but the system (12) is nonstationary, we can write analogously
but of course now Φ(t + 1; t), ∆(t) cannot be expressed in gen-
eral in closed form. Equations of type (14) are encountered fre-
quently also in the study of complicated sampled-data systems
[22]. See Fig. 2
Φ(t + 1; t) is the transition matrix of the system (12) or (14).
The notation Φ(t2; t1) (t2, t1 = integers) indicates transition from
time t1 to time t2. Evidently Φ(t; t) = I = unit matrix. If the system
(12) is stationary then Φ(t + 1; t) = Φ(t + 1 – t) = Φ(1) = const.
Note also the product rule: Φ(t; s)Φ(s; r) = Φ(t; r) and the inverse
rule Φ–1(t; s) = Φ(s; t), where t, s, r are integers. In a stationary
system, Φ(t; τ) = exp F(t – τ).
x(t + 1) = Φ(t + 1; t) + ∆(t)u(t)
y(t) = M(t)x(t)
t = 0, 1, … (14)
As a result of the preceding discussion, we shall represent ran-
dom phenomena by the model
x(t + 1) = Φ(t + 1; t)x(t) + u(t) (15)
where {u(t)} is a vector-valued, independent, gaussian random
process, with zero mean, which is completely described by (in
view of Theorem 5 (C))
Eu(t) = 0 for all t;
Eu(t)u'(s) = 0 if t ≠ s
Eu(t)u'(t) = G(t).
Of course (Theorem 5 (A)), x(t) is then also a gaussian random
process with zero mean, but it is no longer independent. In fact, if
we consider (15) in the steady state (assuming it is a stable sys-
tem), in other words, if we neglect the initial state x(t0), then
Transactions of the ASME–Journal of Basic Engineering, 82 (Series D): 35-45. Copyright © 1960 by ASME
−
1t
x(t) = ∑
−∞=
Φ(t; r + 1)u(r).
Φ(t; r + 1)Q(r) Φ'(s; r + 1).
r
Therefore if t ≥ s we have
−
1s
Ex(t)x'(s) = ∑
r
−∞=
Thus if we assume a linear dynamic model and know the
statistical properties of the gaussian random excitation, it is easy
to find the corresponding statistical properties of the gaussian
random process {x(t)}.
In real life, however, the situation is usually reversed. One is
given the covariance matrix Ex(t)x'(s) (or rather, one attempts to
estimate the matrix from limited statistical data) and the problem
is to get (15) and the statistical properties of u(t). This is a subtle
and presently largely unsolved problem in experimentation and
data reduction. As in the vast majority of the engineering
literature on the Wiener problem, we shall find it convenient to
start with the model (15) and regard the problem of obtaining the
model itself as a separate question. To be sure, the two problems
should be optimized jointly if possible; the author is not aware,
however, of any study of the joint optimization problem.
In summary, the following assumptions are made about random
processes:
Physical random phenomena may be thought of as due to
primary random sources exciting dynamic systems. The primary
sources are assumed to be independent gaussian random
processes with zero mean; the dynamic systems will be linear. The
random processes are therefore described by models such as (15).
The question of how the numbers specifying the model are
obtained from experimental data will not be considered.
Solution of the Wiener problem
Let us now define the principal problem of the paper.
Problem I. Consider the dynamic model
x(t + 1) = Φ(t + 1; t)x(t) + u(t) (16)
y(t) = M(t)x(t) (17)
where u(t) is an independent gaussian random process of n-
vectors with zero mean, x(t) is an n-vector, y(t) is a p-vector (p ≤
n), Φ(t + 1; t), M(t) are n × n, resp. p × n, matrices whose
elements are nonrandom functions of time.
Given the observed values of y(t0), ..., y(t) find an estimate
x*(t1|t) of x(t1) which minimizes the expected loss. (See Fig. 2,
where ∆(t) = I.)
This problem includes as a special case the problems of filter-
ing, prediction, and data smoothing mentioned earlier. It in-
cludes also the problem of reconstructing all the state variables of
a linear dynamic system from noisy observations of some of the
state variables (p < n!).
From Theorem 2-a we know that the solution of Problem I is
simply the orthogonal projection of x(t1) on the linear manifold
Y(t) generated by the observed random variables. As remarked in
the Introduction, this is to be accomplished by means of a linear
(not necessarily stationary!) dynamic system of the general form
(14). With this in mind, we proceed as follows.
Assume that y(t0), ..., y(t – 1) have been measured, i.e., that Y(t
– 1) is known. Next, at time t, the random variable y(t) is
measured. As before let y~ (t|t – 1) be the component of y(t)
orthogonal to Y(t – 1). If y~ (t|t – 1) ≡ 0, which means that the
values of all components of this random vector are zero for almost
every possible event, then Y(t) is obviously the same as Y(t – 1)
and therefore the measurement of y(t) does not convey any addi-
tional information. This is not likely to happen in a physically
meaningful situation. In any case, y~ (t|t – 1) generates a linear
manifold (possibly 0) which we denote by Z(t). By definition,
Y(t – 1) and Z(t) taken together are the same manifold as Y(t),
and every vector in Z(t) is orthogonal to every vector in Y(t – 1).
Assuming by induction that x*(t1 – 1|t – 1) is known, we can
x*(t1|t) = Eˆ [x(t1)|Y(t)] = Eˆ [x(t1)|Y(t – 1)] + Eˆ [x(t1)|Z(t)]
= Φ(t + 1; t) x*(t1 – 1|t – 1) + Eˆ [u(t1 – 1)|Y(t – 1)]
+ Eˆ [x(t1)|Z(t)] (18)
write:
where the last line is obtained using (16).
Let t1 = t + s, where s is any integer. If s ≥ 0, then u(tl – 1) is
independent of Y(t – 1). This is because u(tl – 1) = u(t + s – 1) is
then independent of u(t – 2), u(t – 3), ... and therefore by (16–
17), independent of y(t0), ..., y(t – 1), hence independent of Y(t –
1). Since, for all t, u(t0) has zero mean by assumption, it follows
that u(tl – 1) (s ≥ 0) is orthogonal to Y(t – 1). Thus if s ≥ 0, the
second term on the right-hand side of (18) vanishes; if s < 0,
considerable complications result in evaluating this term. We
shall consider only the case tl ≥ t. Furthermore, it will suffice to
consider in detail only the case tl = t + 1 since the other cases can
be easily reduced to this one.
variable y~ (t |t – 1):
The last term in (18) must be a linear operation on the random
Eˆ [x(t + 1)|Z(t)] = ∆*(t) y~ (t|t – 1) (19)
where ∆*(t) is an n × p matrix, and the star refers to “optimal
filtering”.
The component of y(t) lying in Y(t – 1) is y (t|t – 1) =
M(t)x*(t|t – 1). Hence
y~ (t|t – 1) = y(t) – y (t|t – 1) = y(t) – M(t)x*(t|t – 1). (20)
Combining (18-20) (see Fig. 3) we obtain
x*(t + 1|t) = Φ*(t + 1; t)x*(t|t – 1) + ∆*(t)y(t) (21)
where
Φ*(t + 1; t) = Φ(t + 1; t) – ∆*(t)M(t) (22)
Thus optimal estimation is performed by a linear dynamic
system of the same form as (14). The state of the estimator is the
previous estimate, the input is the last measured value of the
observable random variable y(t) , the transition matrix is given by
(22). Notice that physical realization of the optimal filter requires
only (i) the model of the random process (ii) the operator ∆*(t).
The estimation error is also governed by a linear dynamic
system. In fact,
x~ (t + 1|t) = x(t + 1) – x*(t + 1|t)
= Φ(t + 1; t)x(t) + u(t) – Φ*(t + 1; t)x*(t|t – 1)
– ∆*(t)M(t)x(t)
y~ (t|t – 1)
y(t)
∑
∆*(t)
x*(t|t – 1)
∑
MODEL OF RANDOM PROCESS
x*(t + 1|t)
y (t|t – 1)
M(t)
x*(t + s|t)
Φ(t + s; t + 1)
unit
delay
Φ (t + 1; t)
– I
x*(t + 1|t – 1)
Fig. 3 Matrix block diagram of optimal filter
Transactions of the ASME–Journal of Basic Engineering, 82 (Series D): 35-45. Copyright © 1960 by ASME
where Q(t) = Eu(t)u'(t).
∆* (and thus also for Φ*). Since,
There remains the problem of obtaining an explicit formula for
= Φ*(t + 1; t) x~ (t|t – 1) + u(t)
(23)
Thus Φ* is also the transition matrix of the linear dynamic system
governing the error.
From (23) we obtain at once a recursion relation for the co-
variance matrix P*(t) of the optimal error x~ (t|t – 1). Noting that
u(t) is independent of x(t) and therefore of x~ (t|t – 1) we get
P*(t + 1) = E x~ (t + 1|t) x~ '(t + 1|t)
= Φ*(t + 1; t)E x~ (t|t – 1) x~ '(t|t – 1)Φ*' (t + 1; t) + Q(t)
= Φ*(t + 1; t)E x~ (t|t – 1) x~ '(t|t – 1)Φ'(t + 1; t) + Q(t)
= Φ*(t + 1; t)P*(t)Φ'(t + 1; t) + Q(t)
(24)
x~ (t + 1)|Z(t)) = x(t + 1) – Eˆ [x(t + 1)|Z(t)]
is orthogonal to y~ (t |t – 1), it follows that by (19) that
0 = E[x(t + 1) – ∆*(t) y~ (t|t – 1)] y~ '(t|t – 1)
= Ex(t + 1) y~ '(t|t – 1) – ∆*(t)E y~ (t|t – 1) y~ '(t|t – 1).
Noting that x (t + 1|t – 1) is orthogonal to Z(t), the definition of
P(t) given earlier, and (17), it follows further
0 = E x~ (t + 1|t – 1) y~ '(t|t – 1) – ∆*(t)M(t)P*(t)M'(t)
= E[Φ(t + 1; t) x~ (t|t – 1) + u(t|t – 1)] x~ '(t|t – 1)M'(t)
Finally, since u(t) is independent of x(t),
– ∆*(t)M(t)P*(t)M'(t).
0 = Φ(t + 1; t)P*(t)M'(t) – ∆*(t)M(t)P*(t)M'(t).
Now the matrix M(t)P*(t)M'(t)will be positive definite and hence
invertible whenever P*(t) is positive definite, provided that none
of the rows of M(t) are linearly dependent at any time, in other
words, that none of the observed scalar random variables yl(t), ...,
ym(t), is a linear combination of the others. Under these
circumstances we get finally:
∆*(t) = Φ(t + 1; t)P*(t)M'(t)[M(t)P*(t)M'(t)]–1 (25)
Since observations start at t0, x~ (t0|t0 – 1) = x(t0); to begin the
iterative evaluation of P*(t) by means of equation (24), we must
obviously specify P*(t0) = Ex(t0)x'(t0). Assuming this matrix is
positive definite, equation (25) then yields ∆*(t0); equation (22)
Φ*(t0 + 1; t0), and equation (24) P*(t0 + 1), completing the cycle.
If now Q(t) is positive definite, then all the P*(t) will be positive
definite and the requirements in deriving (25) will be satisfied at
each step.
Now we remove the restriction that t1 = t + 1. Since u(t) is
orthogonal to Y(t), we have
x*(t + 1|t) = Eˆ [Φ(t + 1; t)x(t) + u(t)|Y(t)] = Φ(t + 1; t)x*(t|t)
Hence if Φ(t + 1; t) has an inverse Φ(t; t + 1) (which is always the
case when Φ is the transition matrix of a dynamic system
describable by a differential equation) we have
x*(t|t) = Φ(t; t + 1)x*(t + 1|t)
If t1 ≥ t + 1, we first observe by repeated application of (16) that
x(t + s) = Φ(t + s; t + 1)x(t + 1)
−1
s
Φ(t + s; t + r)u(t + r) (s ≥ 1)
+ ∑
r
1
Since u(t + s – 1), …, u(t + 1) are all orthogonal to Y(t),
x*(t + s|t) = Eˆ [x(t + s)|Y(t)]
= Eˆ [Φ(t + s; t + 1)x(t + 1)|Y(t)]
= Φ(t + s; t + 1)x*(t + 1|t) (s ≥ 1)
If s < 0, the results are similar, but x*(t – s|t) will have (1 –
s)(n – p) co-ordinates.
The results of this section may be summarized as follows:
Theorem 3. (Solution of the Wiener Problem)
Consider Problem I. The optimal estimate x*(t + 1|t) of x(t +
1) given y(t0), ..., y(t) is generated by the linear dynamic system
x*(t + 1|t) = Φ*(t + 1; t)x*(t|t – 1) + ∆*(t)y(t) (21)
The estimation error is given by
x~ (t + 1|t) = Φ*(t + 1; t) x~ (t|t – 1) + u(t)
The covariance matrix of the estimation error is
cov x~ (t|t – 1) = E x~ (t|t – 1) x~ '(t|t – 1) = P*(t)
(23)
(26)
The expected quadratic loss is
n
∑
|(~
xE
t
2
i
i
t
−
=
1
)1
= trace P*(t)
(27)
The matrices ∆*(t), Φ*(t + 1; t), P*(t) are generated by the
recursion relations
∆*(t) = Φ(t + 1; t)P*(t)M'(t)[M(t)P*(t)M'(t)]–1
Φ*(t + 1; t) = Φ(t + 1; t) – ∆*(t)M(t)
P*(t + 1) = Φ*(t + 1; t)P*(t)Φ'(t + 1; t)
In order to carry out the iterations, one must specify the
covariance P*(t0) of x(t0) and the covariance Q(t) of u(t).
Finally, for any s ≥ 0, if Φ is invertible
x*(t + s|t) = Φ(t + s; t + 1)x*(t + 1)|t)
(28)
(29)
(30)
+ Q(t)
t ≥ t0
= Φ(t + s; t + 1)Φ*(t + 1; t)Φ(t; t + s – 1)
+ Φ(t + s; t + 1)∆*(t)y(t)
× x*(t + s – 1|t – 1)
(31)
so that the estimate x*(t + s|t) (s ≥ 0) is also given by a linear dy-
namic system of the type (21).
Remarks. (h) Eliminating ∆* and Φ* from (28–30), a nonlinear
difference equation is obtained for P*(t):
P*(t + 1) = Φ(t + 1; t){P*(t) – P*(t)M'(t)[M(t)P*(t)M'(t)]–1
× P*(t)M(t)}Φ'(t + 1; t) + Q(t) t ≥ t0 (32)
This equation is linear only if M(t) is invertible but then the
problem is trivial since all components of the random vector x(t)
are observable P*(t + 1) = Q(t). Observe that equation (32) plays
a role in the present theory analogous to that of the Wiener-Hopf
equation in the conventional theory.
Once P*(t) has been computed via (32) starting at t = t0, the
explicit specification of the optimal linear filter is immediately
available from formulas (29-30). Of course, the solution of
Equation (32), or of its differential-equation equivalent, is a much
simpler task than solution of the Wiener-Hopf equation.
(i) The results stated in Theorem 3 do not resolve completely
Problem I. Little has been said, for instance, about the physical
significance of the assumptions needed to obtain equation (25),
the convergence and stability of the nonlinear difference equa-
tion (32), the stability of the optimal filter (21), etc. This can
actually be done in a completely satisfactory way, but must be
left to a future paper. In this connection, the principal guide and
Transactions of the ASME–Journal of Basic Engineering, 82 (Series D): 35-45. Copyright © 1960 by ASME
tool turns out to be the duality theorem mentioned briefly in the
next section. See [29].
(j) By letting the sampling period (equal to one so far) ap-
proach zero, the method can be used to obtain the specification of
a differential equation for the optimal filter. To do this, i.e., to
pass from equation (14) to equation (12), requires computing the
logarithm F* of the matrix Φ*. But this can be done only if Φ* is
nonsingular—which is easily seen not to be the case. This is
because it is sufficient for the optimal filter to have n – p state
variables, rather than n, as the formalism of equation (22) would
seem to imply. By appropriate modifications, therefore, equation
(22) can be reduced to an equivalent set of only n – p equations
whose transition matrix is nonsingular. Details of this type will be
covered in later publications.
(k) The dynamic system (21) is, in general, nonstationary. This
is due to two things: (1) The time dependence of Φ(t + 1; t) and
M(t); (2) the fact that the estimation starts at t = t0 and improves as
more data are accumulated. If Φ, M are constants, it can be shown
that (21) becomes a stationary dynamic system in the limit t → ∞.
This is the case treated by the classical Wiener theory.
(l) It is noteworthy that the derivations given are not affected
by the nonstationarity of the model for x(t) or the finiteness of
available data. In fact, as far as the author is aware, the only
explicit recursion relations given before for the growing-memory
filter are due to Blum [12]. However, his results are much more
complicated than ours.
(m) By inspection of Fig. 3 we see that the optimal filter is a
feedback system, and that the signal after the first summer is
white noise since y~ (t|t – 1) is obviously an orthogonal random
process. This corresponds to some well-known results in Wiener
filtering, see, e.g., Smith [28], Chapter 6, Fig. 6–4. However, this
is apparently the first rigorous proof that every Wiener filter is
realizable by means of a feedback system. Moreover, it will be
shown in another paper that such a filter is always stable, under
very mild assumptions on the model (16–17). See [29].
The Dual Problem
Let us now consider another problem which is conceptually
very different from optimal estimation, namely, the noise-free
regulator problem. In the simplest cases, this is:
Problem II. Consider the dynamic system
x(t + 1) = Φˆ (t + 1; t)x(t) + Mˆ (t)u(t)
(33)
where x(t) is an n-vector, u(t) is an m-vector (m ≤ n), Φˆ , Mˆ are
n × n resp. n × m matrices whose elements are nonrandom func-
tions of time. Given any state x(t) at time t, we are to find a
sequence u(t), ..., u(T) of control vectors which minimizes the
performance index
+
1T
V[x(t)] = ∑
=τ
t
x'(τ)Q(τ)x(τ)
Where Qˆ (t) is a positive definite matrix whose elements are
nonrandom functions of time. See Fig. 2, where ∆ = Mˆ and M = I.
Probabilistic considerations play no part in Problem II; it is
implicitly assumed that every state variable can be measured
exactly at each instant t, t + 1, ..., T. It is customary to call T ≥ t
the terminal time (it may be infinity).
The first general solution of the noise-free regulator problem is
due to the author [18]. The main result is that the optimal control
vectors u*(t) are nonstationary linear functions of x(t). After a
change in notation, the formulas of the Appendix, Reference [18]
(see also Reference [23]) are as follows:
u*(t) = – ∆ˆ *(t)x(t)
(34)
Under optimal control as given by (34), the “closed-loop” equa-
tions for the system are (see Fig. 4)
x(t + 1) = Φˆ *(t + 1; t)x(t)
and the minimum performance index at time t is given by
V*[x(t)] = x'(t)P*(t – 1)x(t)
The matrices ∆ˆ *(t), Φˆ *(t + 1; t), Pˆ *(t) are determined by
'(t) Pˆ *(t) Mˆ (t)]–1 Mˆ '(t) Pˆ *(t) Φˆ (t + 1; t)
PHYSICAL SYSTEM TO BE CONTROLLED
the recursion relations:
∆ˆ *(t) = Mˆ[
Φˆ *(t + 1; t) = Φˆ (t + 1; t) – Mˆ (t) ∆ˆ *(t)
Pˆ *(t – 1) = Φˆ '(t + 1; t) Pˆ *(t) Φˆ *(t + 1; t)
Initially we must set Pˆ *(T) = Qˆ (T + 1).
Φˆ (t + 1; t)
unit
delay
x(t + 1)
u*(t)
∑
– I
(t)
Mˆ
x(t)
Fig. 4 Matrix block diagram of optimal controller
+ Qˆ (t)
t ≤ T
(35)
(36)
(37)
∆ˆ *(t)
bles.
Comparing equations (35–37) with (28–30) and Fig. 3 with
Fig. 4 we notice some interesting things which are expressed
precisely by
Theorem 4. (Duality Theorem) Problem I and Problem II are
duals of each other in the following sense:
Let τ ≥ 0. Replace every matrix X(t) = X(t0 + τ) in (28–30) by
Xˆ '(t) = Xˆ '(T – τ). Then One has (35–37). Conversely, replace
every matrix Xˆ (T – τ) in (35–37) by X'(t0 + τ). Then one has
(28–30).
Proof. Carry out the substitutions. For ease of reference, the
dualities between the two problems are given in detail in Table 1.
Table 1
1
2
Problem I
x(t)
(unobservable)
state
variables of random proc-
ess.
y(t) observed random varia-
Problem II
x(t) (observable) state varia-
to be
bles of plant
regulated.
u(t) control variables
t0 first observation.
3
4 Φ(t0 + τ +1; t0 + τ) transition
5 P*(t0 + τ) covariance of
T last control action.
Φˆ (T – τ + 1; T – τ) transi-
matrix.
tion matrix.
optimized estimation error. Pˆ *(T – τ) matrix of quad-
ratic form for performance
index under optimal regu-
lation.
∆ˆ *(T – τ) weighting of
state for optimal control.
Φˆ *(T – τ + 1; T – τ) transi-
tion matrix under optimal
regulation.
Mˆ (T – τ) effect of control
vectors on state.
Qˆ (T – τ) matrix of
quadratic form defining
error criterion.
6 ∆*(t0 + τ) weighting of ob-
servation for optimal esti-
mation.
7 Φ*(t0 + τ + 1; t0 + τ) transi-
tion matrix for optimal es-
timation error.
8 M(t0 + τ) effect of state on
9 Q(t0 + τ) covariance of ran-
Remarks. (n) The mathematical significance of the duality be-
tween Problem I and Problem II is that both problems reduce to
the solution of the Wiener-Hopf-like equation (32).
dom excitation.
observation.
(o) The physical significance of the duality is intriguing. Why
are observations and control dual quantities?
Transactions of the ASME–Journal of Basic Engineering, 82 (Series D): 35-45. Copyright © 1960 by ASME