Analysis of Costate Discretizations in Parameter Estimation for Linear Evolution Equations.
A widely used approach to parameter identification is the output least-squares formulation. Numerical methods for solving the resulting minimization problem almost invariably require the computation of the gradient of the output least-squares functional. When the identification problem involves time-dependent distributed parameter systems (or approximations thereof), numerical evaluation of the gradient can be extremely time consuming. The costate method can greatly reduce the cost of computing these gradients. However, questions have been raised concerning the accuracy and convergence of costate approximations, even when the numerical methods being used are known to converge rapidly on the forward problem.  In this paper it is shown that the use of time-marching schemes that yield high-order accuracy on the forward problem does not necessarily lead to high-order accurate costate approximations. In fact, in some cases these approximations do not converge at all. However, under certain circumstances, rapidly converging gradient approximations do result because of rapid weak-star-type convergence of the costate approximations. These issues are treated both theoretically and numerically.
Introduction
. In this paper we analyze temporal discretizations of the costate
method for computing gradients in the output least-squares approach to parameter esti-
mation. This analysis applies to initial value problems of the form
Here A(q) is a bounded linear operator on a Hilbert space H, and the inner product on H
is denoted h\Delta; \Deltai H . The dot over u indicates differentiation with respect to t. We will refer
to u as the state variable and to (1.1) as the state equation. We assume q lies in set QAD of
admissible parameters contained in a normed linear "parameter" space Q. Throughout the
paper we assume that the map q 7! A(q) is Gateaux differentiable in the operator norm.
In applications of interest, (1.1) is a finite dimensional (i.e., dim(H) ! 1) approximation
of a time-dependent PDE. An important example is the diffusion equation
@t
Department of Mathematical Sciences, Montana State University, Bozeman MT 59717. Research was
supported in part by the NSF under Grant DMS-9106609 and by the Center for Interfacial Microbial Process
Engineering at Montana State University, an NSF-sponsored Engineering Research Center, and the Center's
Industrial Associates.
y Center for Applied Mathematical Sciences, University of Southern California, Los Angeles, CA 90089-
1113. Address as of July 1, 1992: Department of Mathematics, Texas A & M University, College Station,
77843-3368. Research was supported in part by the Air Force Office of Scientific Research under grants
AFOSR-90-0091 and AFOSR-86-0126, and by the NSF under grants DMS-8818530 and DMS-8704169.
In these situations, dim(H) can be arbitrarily large.
In general we assume that the solution u lies in the "state space"
which is a Hilbert space with inner product
As in Banks and Kunish [3], we assume the existence of an "observation space" Z, which is
a Hilbert space with inner product hh \Delta; \Delta ii Z , and an "observation operator"
Given an observation z of u, one wishes to estimate the parameter q. In the output least-squares
approach to parameter estimation, q 2 QAD is selected to minimize the functional
where u(q) is the solution to (1.1). Here the scalar ff is a positive regularization parameter
and N 2 (q) is a regularization functional, whose purpose is to stabilize the minimization.
Usually, N (q) is a norm or seminorm on Q.
Computational methods to minimize T ff (q) typically require gradients or gradient ap-
proximations. For example, assuming a discretization of the parameter of the form
one often (see, for example, [3, Section V.6]) approximates the components of the gradient
using finite differences, e.g.,
for - a small scalar. Note that in the limit as - ! 0, the i th component of rT ff (q) is the
directional (i.e., Gateaux) derivative of T ff (q) in the direction / i . When dim(H) and n q are
large, gradient approximations based directly on (1.5) are extremely expensive, requiring
evaluations of the functional T ff , and hence, n q solutions of the evolution equation
(1.1).
An attractive alternative is the costate approach [3, Section V.5], [7], which is described
in its continuous form in Section 2. In sharp contrast to the (n q + 1) state equation solutions
needed for the finite difference method (1.5), the costate method requires the solution of
only two evolution equations, the state equation and the "costate" or "adjoint" equation,
followed by n q H-inner products.
Unfortunately, computational experience with adjoint approximations in parameter estimation
for certain evolution equations has been puzzlingly disappointing. For example, the
authors of [4] concluded that this approach was unsatisfactory because extremely
difficult to obtain accurate search directions with gradients computed in this manner". Also,
in [2] authors the suggest costate method for a certain class of damped elastic systems and
explain how it may facilitate the efficient use of gradient-based methods such as conjugate
gradients or BFGS. They also report having used it in some of their numerical experiments.
However, in none of the specific examples upon which they report did they use this method.
Subsequent private communication with two of the authors revealed that they encountered
unexpected difficulty with the costate approach. Later commenting on this, in [1] Banks
states that experienced so much difficulty with the costate based gradient methods
that we abandoned them
Similar difficulties have been encountered in optimal control settings. In [5], which dealt
with computational methods for control systems governed by delay differential equations,
the authors observed that certain spline approximations to the solution of operator Ricatti
equations failed to converge strongly. In a subsequent work [6] it was carefully shown that
the underlying difficulty was lack of strong convergence of the adjoint approximations.
In this paper we focus on the temporal discretizations of the costate system in least-squares
parameter estimation. Since the costate method of computing rT ff (q) requires the
solution of two evolution equations plus n q inner products, two important components of
any numerical scheme for approximating the gradient are (a) the scheme for approximating
the solution of the evolution equations (which we assume is a time-marching scheme), and
(b) the numerical quadrature scheme. If both of these are, say, - th -order accurate, then
one would reasonably hope to attain - th -order convergence of the gradient approximations.
However, if the observation operator C involves pointwise evaluation in time, then certain
subtleties arise and one may observe unexpectedly poor convergence of the gradient approxi-
mations. Perhaps the most striking examples which we present below involve the fourth-order
Runge-Kutta (RK4) time-marching scheme. We show that if RK4 is used in conjunction with
Simpon's rule (which is fourth-order) for numerical quadrature, then the gradient approximations
fail altogether to converge. We show that in fact these RK4/Simpson's approximations
converge with second-order accuracy to the (3=2) times the true gradient! And more positively
but perhaps equally surprisingly, we prove that RK4 together with the second-order
accurate "trapezoidal" quadrature rule yields fourth-order convergence!
The underlying reason for these strange phenomena is that when C involves pointwise
evaluation in time, the costate equation is an evolution equation with Dirac-delta functions in
the forcing term. For this reason, the costate approximate solutions do not converge strongly.
However, in some cases they do exhibit high-order convergence in the weak   topology of the
dual space of C - .
These considerations are similar in spirit to the work reported in [6]. In the introduction
of that paper, the authors state "We feel that many distributed parameter control
systems are such that might lead to numerical difficulties'' when used
in optimization schemes, and "we hope that the reader will be motivated to think about
similar problems for more complex distributed parameter systems". While that paper addresses
the optimal control problem, it is nonetheless relevant to us to the extent that output
least-squares parameter estimation is mathematically similar to optimal control.
Banks [1] has suggested that non-convergence of costate approximations can also occur
in parameter estimation because of a lack of convergence of the spatial discretization. For
example, in systems such as (1.2), one must discretize with respect to each component of x as
well as with respect to t. Then a result needed for the costate approximations is essentially
convergence the dual semigroups. These considerations suggest directions for future
work, but we do not pursue them here.
In Section 2 we discuss the costate method for the continuous problem, apart from
any temporal approximations. In Section 3 we introduce a fairly wide class of time-marching
schemes for (1.1) and a corresponding costate approximation scheme. We then analyze these
approximations in some detail, making certain assumptions about the properties of the time-marching
scheme, the observation operator, and the various operator discretizations. Section
contains some numerical examples which illustrate the results of the analysis. In section
5, several alternative approaches which avoid the difficulties mentioned above are presented.
Conclusions of this work are discussed in Section 6.
2. The Continuous Costate Approach. As stated in the introduction, typical implementations
of the output least-squares approach to parameter estimation require gradients
of the functional T ff (q) given in (1.4). Most of the effort in obtaining this gradient lies in
computing the gradient of the least squares functional
where solves the state equation
Since the gradient of T (q) can be obtained by computing directional derivatives, we
focus throughout on the computation of the directional derivative of T at q in the direction
p. This is given by
The residual is defined by
and hence,
The solution to (2.2) is given in terms of the "solution operator" S(q) by
Z te A(q)s
In terms of S(q), we have
Thus
Taking adjoints of C and S(q), one obtains
where y is the solution to the "costate equation"
The operator S   (q) is given, for
[S
e \GammaA   (q)s
The differential equation which y satisfies is
Since this equation has a final condition instead of an initial condition, it is useful to note
that it is equivalent to
y:
Here J is the "time reversal" operator on H defined by
The costate approach to computing rT (q) then consists of the following four steps.
The Continuous Costate Method:
(i) Solve the state equation (2.2) for the state variable
(ii) Compute the residual
(iii) Solve the costate equation (2.9) for the costate variable y.
(iv) Compute the directional derivatives of T in directions p (c.f. equation
(1.5), where
Note that the this approach requires the solution of only two evolution equations, the
state equation (2.2) and costate equation (2.9), as opposed to the (n q +1) evolution equations
of (1.5). Moreover, if A is linear in q (as is the case with the approximations of (1.2) for
example), then
so that the implementation of step (iv) is straightforward. Finally, we note that to complete
the computation of rT ff (q) it is necessary to compute the gradient rN 2 (q) of the
regularization functional, but this is usually trivial.
3. Analysis of Discrete Approximations.
3.1. The Discrete Approximations. Usually in solving (2.2) numerically, one uses
some sort of finite-difference or time-marching scheme (tms). In this section, we pose and
analyze a rather natural algorithm for computing costate approximations based on a given
tms and the Continuous Costate Method described above.
For a given n, let
denote a specified mesh on [0; t F ], and define h k
k . For later use in the discussions
of asymptotic rates of convergence, we define
Also, let u
k is the approximation to the state variable u at t n
k which is
obtained using a particular tms with this mesh. If we denote ff(t n
k )g by f n , then we may
express this approximation by
This approximation is then used in a minimization scheme for T ff (q) in the parameter
estimation problem. Corresponding to the continuous least-squares functional in (2.1), define
Here z n , Z n , and C n are discretizations of z, Z, and C which are discussed in detail below. If
the minimization scheme requires directional derivatives (i.e., gradients), then a seemingly
natural approach ([3, Section V.5]) which we call the "discretized costate approximation",
based on an obvious discretization of the continuous costate method, suggests itself:
The Discretized Costate Approximation:
(i) Compute
(ii) Compute the residual r
(iii) Compute y n by
~
Here, J n is the approximation of J (c.f. (2.10)) given by
and (C   ) n is an approximation to the adjoint of C
(iv) Approximate directional derivatives in directions p by
Here, hh \Delta; \Delta ii H n is a discretization of the continuous inner product (1.3).
The convergence of g
depends on factors such as the convergence of the
time marching scheme tms on the forward problem, the convergence of approximations C n
to the observation operator C, and convergence of the discrete inner products hh \Delta; \Delta ii H n and
to the continuous inner products hh \Delta; \Delta ii H and hh \Delta; \Delta ii Z , respectively. To carry out a
convergence analysis, we first carefully state the discrete version of the problem.
O
H:
For
and for f
k=0 and
k=0 in H n , the inner product is
where fw n
k=0 is a given sequence of weights. The purpose of these weights is to facilitate
definition of hh \Delta; \Delta ii H n in such a way that it well-approximates hh \Delta; \Delta ii H . Accordingly, we
assume, essentially without loss of generality, that there exists an independent of
n, for which
This forms the basis for the following lemma:
Lemma 1. If a sequence ff n g 2 H n has the property that for some constant - ? 0,
0-k-n
kf n
then
kf n kH
We also denote by hh \Delta; \Delta ii E n the inner product which, loosely speaking, uses the standard
Euclidean inner product for the time component. Specifically,
This inner product will be used below in adjoint computations. The two inner products are
related by
where W n is the diagonal operator on H n defined by
We assume that the time marching scheme tms is of the form
where S n (q) is defined by recursion, for
for some bounded operator B k (q). Expressed in terms of components,
See Section 4 below for specific examples.
We make the following assumptions on the tms:
(A1) The tms is stable [10], i.e., there exists
Y
The terms in this product are ordered from left to right with decreasing indices.
(A2) For some constant - ?
as in (A2), for any f 2 C - ([0; t F ]; H) we have
As f(t ds +O(jhj -
(A4) In the state equation (2.2), the forcing function f lies in C - ([0; t F ]; H).
The following theorem is then readily established.
Theorem 1. If (A1)-(A4) hold, then
using (A2)-(A4), (2.2), and
(3.12) one obtains
From this one can use (A1) and induction to show that
0-k-n
The theorem now follows from Lemma 1.
3.2. The Rate of Convergence of the Gradients. We now address the convergence
of the (state) directional derivative approximations . The results obtained in this
subsection will later be used in proving convergence of the discrete costate approximations
We assume the existence of operators
Also, note that C . We make the following assumptions on the relationships
between and the inner product approximations.
(A5) For
We now consider the convergence of the directional derivatives. From (2.5) and (3.1),
where
If (A6) holds, the rate of convergence of e 1 is directly determined by the rate of convergence
of solutions to the forward problem, while that of e 3 is assured under (A5). The rate
convergence of e 2 is determined by how fast
Taking directional derivatives in the component form of the recursion (3.12),
On the other hand,
As f(t n
so taking directional derivatives gives
Comparing (3.19) and (3.20), one can show the desired convergence provided the directional
derivative operator preserves convergence rates and smoothness in t. More precisely, we
assume
R h ds +O(jhj -
We then obtain the following analogue of Theorem 1.
Theorem 2. If assumptions (A1)-(A4) and (A2')-(A4') hold, then
Proof: For notational convenience, we suppress dependence on q. First we show that each
k kH is bounded independently of k or n. Applying (3.19) recursively,
where
From assumption (A1),
But for each j, by assumptions (A2')-(A3') and Theorem 1,
Z h je As f(t
The integral term in the above inequality is bounded independently of j as a consequence of
assumption (A4').
Subtracting (3.20) from (3.19),
As f(t n
ds
From the boundedness of the
k 's, the result of Theorem 1, and the assumptions, all the
terms on the right side of this equation are O(jhj -+1 ), except the first term. Thus this
equation has the form (3.13), and the same argument used in the proof of Theorem 1 apply.
From (3.15) one obtains the following result.
Theorem 3. Under assumptions (A1)-(A6) and (A2')-(A4'),
Proof: Referring to equation (3.16), from Schwartz's Inequality and the definitions of r,
r n , and z n ,
The first term within the above parentheses is O(jhj - ) by Theorem 1, while the second is
(A6). The kffi p u n kH n are bounded by Theorem 2. Consequently, je
Similarly, from equation (3.17), Theorem 2, and assumption (A6),
Finally, from (3.18) and (A5), je 3
3.3. Continuous-time versus Discrete-time observations. In the following sub-section
we analyse the convergence of the costate approximations. In this subsection, we
make some preliminary considerations toward that end, leading to a specification of two
distinct classes of observation operators C.
For brevity of notation, we set
From (2.11) and (3.4) one obtains
As a consequence of Theorem 1, the boundedness of ffi p A, and the fact that operators ffi p A
and P n commute, we have
If the sequence fky n kH ng is bounded, then Schwartz's Inequality and (3.24) imply that
ng is bounded if jE Hence, the rate of convergence of the
costate approximations to the directional derivatives depends on the rate at which
From (2.9) we see that
where the operator has a representation
Z te A(q)   s
Similarly, from (3.2)
are the operators which are obtained if A is replaced by A   (q) in the formation
of S n and R n ; c.f. (3.10)-(3.11). The operator (C   is an approximation
to the adjoint of C Z. For now we select
the adjoint of the operator different choice might be made on the basis of
the discussion in Section 5 below.
Since the costate vector y is given in terms of an evolution equation with C   r as the
source term (and a similar statement is true for the approximations), we find it convenient
to specify C further. This facilitates analysis of the convergence of jE 2 j. It commonly
happens in applications that C inherits a type of tensor-product structure from the state
space In particular, we assume the existence of a bounded "spatial
observation operator"
where Z is a Hilbert space related to Z as descibed below. We distinguish two cases which
are of practical importance.
Observations continuous in time: In this first case, the observation operator C is "con-
tinuous in time". It is defined in terms of a space Z and a bounded operator
which acts at each t by
In this case,
O
and is defined by
In addition, we define Q n and the inner product approximation in a manner analogous
to (3.5) and (3.6):
and
Note that in this case, the validity of assumption (A5) is determined by the quadrature
weights w n
k , which also determine the accuracy by which hh \Delta; \Delta ii H n approximates
Assumption (A6) is trivially satisfied in this case since C n P
From (3.28) and the definition of the adjoint,
The residuals r and r n appearing here are given by
Observations discrete in time: In the second case we assume that the observation operator
C consists of some spatial observations taken at m discrete points -
in time. Accordingly, we assume a discrete observation space Z of the form
O
Z:
For simplicity we assume that, given the set of temporal observation points f- i g,
the time-marching grid t n is always chosen so that f- i g ae t n , so that there is always
an injective map
ng
such that for f 2 C((0; t F
In this case, since the observation space is already discrete, we define
the identity on Z;
(3.
d
In this case, assumption (A5) is always true since Z Also, (A6) is always
true since, as in the continuous-time observation case, C n P C. (This is true
since we have assumed that f- i g ae t n .)
From (3.29) and the definition of the adjoint, we find that C   : Z 7! H is given, for
where ffi(\Delta) denotes the Dirac delta function. The discrete analogue of (3.32) is
where ffi -;k denotes the Kronecker delta function for integer pairs.
3.4. The Rate of Convergence of the Costate Approximations. As discussed in
the previous subsection, immediately following (3.23), it suffices to consider the convergence
of the term jE 2 j to zero. We address the simple case, that of continuous time observations,
first. In this case, - th -order convergence of E 2 results if the observed data z is smooth
enough.
Theorem 4. If (A1)-(A6) and (A2')-(A4') hold and if the data z lies in C - ([0; t F ]; Z),
then
Proof: From (3.25) and (3.25),
Note that J , C, and C   each preserve smoothness with respect to t. From the smoothness
of z and the smoothness of u (which follows from (A4)), JC   r 2 C -+1 ([0; t F ]; H). Applying
Theorem 1 with JC   r in place of the forcing function f and A   in place of A, and noting
that
the last term in (3.34) is O(jhj - ). Also, by (A1) and (A3), kS n
k is bounded independently
of n. From (3.35) and (3.36), it suffices to show
But
)), so this follows from Theorem 1 and the
boundedness of C and C   .
With the result of this theorem and an assumption that the quadrature weights are
chosen correctly (c.f. (3.6)), we obtain costate convergence:
Corollary 1. If the hypotheses of Theorem 4 hold and if
for any g 1 and g 2 in C - ([0; t f ]; H), then
Proof: Theorem 4 implies the sequence fky n kH ng is bounded. From the discussion immediately
following (3.23), it is suffices to show that
The first term is O(jhj - ) by Theorem 4. Since y and OE are sufficiently smooth, the second
term is O(jhj - ) by (3.37).
We next address the case of discrete time observations. In this case, C   r is given by
a linear combination of Dirac delta functions - (3.32). Since this is the forcing term in
the costate equation, the question of how fast (if at all) g
rather delicate. In
particular, the smoothness of C   r played a crucial role in the proof of Theorem 4. But in
the present case, C   r is not smooth at all; it only exists as a distribution. Thus the proof of
Theorem 4 (and hence, Corollary 1) is not valid in this case.
On the other hand, the convergence in (3.33) is stronger than what is actually needed.
It is essentially a statement of strong, or pointwise, convergence in H. But from (3.23), we
see that a sort of weak   convergence would be sufficient. Rewriting E 2 as
and noting that OE 2 C -+1 ([0; t F ]; H), we need only establish - th -order weak   convergence of
n y to y in the dual space of C -+1 ([0; t F ]; H).
There are seemingly natural situations in which this weak   convergence does not apply,
and a number of subtleties arise in our exploration of these matters. These are discussed
shortly. There is, however, one class of problems which we fully analyze here. It covers a
variety of second-order methods with uniform meshes t n .
Theorem 5. Assume that
(a) The assumptions (A1)-(A4) hold;
(b) In (A2)-(A4) we have
(c) The observation points - i all lie in the interior of (0; t F ).
(d) The partition t n is uniform, in the sense that h
(e) The quadrature weights w k in equation (3.6) are given by w
weights correspond to the "trapezoid
rule" for numerical quadrature);
(f) The operators B k in terms of which S n is defined (cf. (3.11)) are the same for
all k; there is an operator B for which
(g) The operator R n may be written, for
(Note from (3.11) that the last component, [R no role in the tms.
Our purpose in defining it this way here will become clear below.)
(h) The operators B, R 0 and R 1 are all rational functions of A(q).
Then
Proof: As in Corollary 1, it suffices to show This leads us to the study of
From (3.26), we see that
We shall take adjoints in the inner product on the right side of this equation, and the main
effort in this proof lies in the subsequent simplification. In particular, our goal is to show
that
We now proceed with the details. First note from (3.9) that
But since [P n OE] 0, from condition (e) above we obtain
Referring to (3.42) and (3.44) and noting that J n J n equals the identity, we obtain
J
J
(R n
denotes the adjoint with respect to the E n -inner product defined in (3.9). We have
used here the fact that J n is self-adjoint with respect to hh \Delta; \Delta ii E n .
Next we assert that
(R n
From (3.11), (3.39), and (3.40), the operators S n and R n have the following block Toeplitz
matrix representations:
I
.
and
R
.
The blocks are operators on H. The operators S n
and R n
have an identical block Toeplitz
form, except B, R 0 , and R 1 are replaced by their H-adjoints. From consideration of these
block matrix representations, the action of J n , and the E n -adjoint, one can verify that (3.46)
and (3.47) hold. Thus (3.45) simplifies to
Thus the question arises as to how much S n R n and R n S n differ. Since B, R 0 and R 1
are rational functions of A(q), they commute with each other. This and computations based
directly on (3.11) and (3.40) reveal that, for arbitrary
[S (R
for
(R
(R
for
These results may also be obtained directly from the matrix representations (3.48) and (3.49).
Thus the commutator
Equation (3.51) holds for any g 2 H n . In the case of present interest, namely (3.50), the
role of g is played by P n OE, and [P n OE] Thus from (3.50) and (3.51), we have
But the observation points - i are all less than t F and yet occur at mesh points - see the
discussion leading up to (3.30). Thus [(C n )   r n and the last term in this equation
drops out. Also, since [S n R n P n OE] we obtain from condition (e)
Combining this with (3.42) yields (3.43).
Theorem 1. And by reasoning similar to that in the proof of
Theorem 1, since OE is smooth, we also have kS Consequently,
the right side of (3.43) may be replaced by hh Z to within second-order. Using
this, the definition of E 2 in (3.23) and the fact that hh we find that
Now the result follows from assumption (A6).
A result similar to Theorem 5 can be obtained for certain fourth-order methods (includ-
ing Runge-Kutta) on a uniform mesh, provided that the forcing function f is zero at
and provided that appropriate quadrature weights are chosen. Surprisingly, these weights
correspond to second order quadrature rather than fourth order.
Theorem 6. Assume that the forcing term f satisfies and that
(a) Same as Theorem 5, hypothesis (a).
(b) In (A2)-(A4) we have
(c) Same as Theorem 5, hypothesis (c).
(d) Same as Theorem 5, hypothesis (d).
Same as Theorem 5, hypothesis (e).
(f) Same as Theorem 5, hypothesis (f).
(g) The operator R n may be written, for
Same as Theorem 5, hypothesis (h).
Then,
Proof: For given g 2 H n , [R n g] is defined for all k ? 0 but [R n g] 0 is unspecified. Define
f
R n by
Note that f
(because (A3) is satisfied by R n and not necessarily by f
R n ), but that
by (c). Therefore we may replace R n by f
R n in (3.26). Since the
block matrix representation of f
R n is Toeplitz, by following reasoning similar to that leading
to (3.50) in the proof of Theorem 5 we obtain
R
As in the proof of Theorem 5, we now consider the extent to which f
R n and S n fail to
commute. Our goal here is to show that
R
Then we may invoke reasoning similar to that in and following (3.52). With
f
R
we have
From (3.11), (3.54) and the fact that after lengthy but straightforward
computations that for
Also, from the definitions of S n , R n and f
R n and from (A3), we find that for 1 - k - n,
Z he As OE(t n
ds
Substituting these into (3.56) yields, for
\Gammah
Z he As OE(t n
ds +O(h 4 )
The integral term in this equation is approximated to within fourth-order accuracy if
we replace OE(t) with its second-order Taylor polynomial -(t) about
arguments used to show fourth-order convergence of Simpson's rule):h
Z he As OE(t n
Z he As -(t n
ds +O(h 4
Now, -(t) can be evaluated at t n
so that we can use assumption (A3)
to make the substitutionh
Z he As OE(t n
by hypothesis, so that u 0 which means that
Thush
Z he As OE(t n
h
Also, OE(t n) and OE(t n) can be replaced by -(t n), and -(t n), respectively, to within third-order
accuracy in (3.58). This yields
But by (g) and the fact that (which is true by (A2)).
Thus we arrive at
and so we have shown that for 1
R
Using this, (3.55) and reasoning similar to that in (3.53) we obtain
It is also natural to ask what happens if we use a fourth-order tms with a fourth-order
quadrature scheme, such as Simpson's rule. We show that the gradient approximations thus
computed fail to converge, although their directions in the space Q converges. In particular,
the RK4/Simpson costate approximations converge quadratically to (3/2) times the true
gradient.
Corollary 2. Assume that all of the hypotheses of Theorem 6 hold except (e), and
assume further that
(e') The quadrature weights are those arising from Simpson's rule, so that
(i) The limit n ! 1 is taken by repeated doubling of n. Precisely, there is an
integer n 0 such that each n equals n 0 2 l for some positive integer l.
Then,
By arguments similar to those used to obtain (3.55), we find that
R
R n is as in (3.54). Define
T rap is the weight matrix corresponding to the trapezoidal rule:
We shall show that k~v \Gamma vkH which implies that
R
We show by induction that k~v \Gamma vkH We note first that ~ v
assume that for some k - 1 we have ~ v We note from (3.11) that for any
Here, fw n
k g are the Simpson weights - the diagonal entries of W n
Simp . Thus,
But from the definition of W n
Simp , the smoothness of OE and the fact that
follows that the second parenthesized term above is equal to h( 2hOE(t n
(A2) we have Using these observations and the definition of W n
and dividing by 2, we thus find that
By the inductive hypothesis, v
so we may replace v k\Gamma1 by ~
v k\Gamma1 on the
right side of this equation. The resulting expression is then, by definition, the formula for
~
. Thus we have
Since it is also true that
we obtain ~ v holds. We may
then use the reasoning leading to (3.59) to obtain
Next we turn our attention to (C n )   r n g. Referring to (3.3) we see that g depends
linearly on f1=w k g. We show that if the Simpson weights are replaced by the trapezoidal
weights, then the result - call it ~
(3=2)~g. Indeed, this follows easily from
hypothesis (i) by which we are guaranteed that -(i) is even for all so that the
Simpson weights w k satisfy w m. But (2=3)h is simply (3/2)
times the -(i) th trapezoidal weight. Thus the right side of (3.62) is equal, to within O(h 2 ), to
(3=2) times the right side of (3.55). Thus the result of Theorem (6) applies, and we obtain
the desired result.
4. Numerical Examples. To illustrate the results of the previous section we focus on
a spline-based Galerkin approximation (with fixed dimensionality) of the diffusion equation
(1.2) in one space dimension and t directional derivatives for the corresponding
least-squares functional (2.1). Accordingly we
1), and for fixed N - 1 take as
basis functions fb i g N
the usual piecewise linear "hat" functions on a uniform mesh. With
the nodal points x i
1), these are given, for 1 - i - N , by
We seek approximations u of the solution of the diffusion equation (1.2) as linear combinations
of these basis functions, i.e.,
For fixed t we denote the coefficient vector fu i (t)g simply by u(t). An N-dimensional linear
system of the form (2.2) is obtained for u from the standard Galerkin method by requiring
that u(x;
Z 1-
Z 1/
\Gammaq(x)
@x
dx
dx:
In our implementation, the integrals on the left side of this this equation are encoded exactly,
and the ones on the right side are approximated using the trapezoid rule.
We report results below for six numerical experiments. In each of these, 3. Thus
H is the subspace of L 2 (0; 1) consisting of continuous functions on [0; 1] which equal zero at
the endpoints and which are piecewise linear with interior nodes at x
One of the experiments is an example of "continuous time observations", and the
rest are based on "discrete time observations". In each of these cases, the spatial observation
operator C, as in (3.27), is the same. It is defined as pointwise evaluation of functions in H
at the points
The function q(x) and the perturbations p(x) are represented as linear splines on the
described above. Perturbations p are elements of H; specifically, they are zero
at the endpoints and can vary at any of the x 3.
Thus the q-dependent linear evolution equation is given by M -
and
.
The "data" z are generated by carrying out the approximations as outlined above, with
approximated by its linear spline interpolate. The source term f(x; t) is chosen so that the
solution of u(x; t) of (1.2) is u(x; . Then the directional derivatives (2.3)
and their approximations (3.4) are investigated based on q(x) j 1.
Three of the experiments are based on non-uniform meshes t n . To explain how these
meshes are generated, it is sufficient that we explain how the vectors fh k g are generated.
To highlight the dependence on n here we denote these vectors by h (n) , and the individual
h-values by h (n)
k . To generate non-uniform h (n) we begin with a given n and create a uniform
mesh by setting h (n)
=n. Then with a given "mesh ratio" r satisfying 1, the
"mesh refinement scheme" for generating h (2n) from h (n) is
n. The choice refined meshes which are each uniform. This is
clearly not true for any other choice of r.
We examine the convergence properties of costate approximations using four different
time-marching schemes in (3.2). For the state equation (2.2), these methods are:
CN. The trapezoidal, or "Crank-Nicholson", method
RK2e. The explicit second-order Runge-Kutta method
RK2i. The implicit second-order Runge-Kutta method
RK4. The explicit fourth-order Runge-Kutta method
where
In order to express the RK4 method in a form compatible with hypothesis
(g) of Theorem 6, we first obtain a fourth-order accurate approximation of
f k+1=2 . This we do by evaluating at (t n
k+1 )=2 the cubic polynomial which
interpolates
2. For a uniform mesh, this yields
We carried out a rather extensive numerical investigation, using the "MATLAB" software
package [9] on a Sun Sparcstation 2 and a DECstation 5000/200. In all cases, we
compare the costate approximation g
to centered-difference approximations
These are computed by evaluation of the expression under the limit sign in (2.3) with T
replaced by T n , and with
ffl. The value of ffl (called "eps" in MATLAB) is the machine-dependent
constant defined as the smallest integer-power of 2 for which "1
On the Sparcstation 2, so that numbers are represented up to sixteen
decimal places. The choice of
ffl in the finite difference computation reflects a desire to
strike a balance between discretization error and roundoff error. By heuristics in [8, pages
31,32], we expect that (ffi p T n ) FD thus computed should agree with the correct value of
up to roughly half of the sixteen decimal places.
By the result of Theorem 3, the true gradient of the discrete functional
with full accuracy to ffi p T (q) (under the assumptions (A5') and (A6'), of course). Since
to within eight digits or so, the triangle inequality implies that the
rate of convergence of g
correctly indicated by the rate of convergence of g
to long as the relative difference of these two quantities is larger than about
For each of the six experiments, then, we tabulate the relative error
1-i-m
for the four methods as jhj varies. From these results we then estimate the actual rate of
convergence -
- based on the assumption that E
Case 1. Observations are continuous in time, and the meshes are non-uniform with refinement
1=3. The trapezoidal quadrature rule is used for the H n inner
product for the three second-order methods, and Simpson's rule is used for the RK4
method.
6.67e-3 1.6e-6 2.6e-4 2.0e-4 2.9e-6
4.44e-3 5.2e-7 7.8e-5 6.8e-5 2.1e-7
- 2.7 2.7 2.7 6.4
Case 2. This experiment is motivated by Theorem 5. We take observations discrete in time
on uniform meshes. The quadrature weights, as in (3.6), are those arising from the
trapezoidal rule. The results are presented in the following table.
1.0e-2 8.8e-10 8.7e-5 6.9e-5 3.9e-6
The three Runge-Kutta methods appear to be converging at a second-order rate,
as predicted by Theorem 5. Crank-Nicholson costate approximations agree with
the forward-difference approximations to within nine decimal places - roughly the
expected accuracy of (ffi p T n ) FD . In fact, for uniform meshes, g
computed by
the Crank-Nicholson method exactly equals (neglecting roundoff error). We
outline a proof of this assertion: The recursion formula (3.12) for
takes the form
with
Then,
Using these in (3.19) leads to
Referring to and (3.22) and (4.2) we see that this last equation is just the recursion
But by reasoning similar to that which led to (3.43), the right side of this equation
can be rewritten as
which by (3.4), (3.22) and (3.26) is the same as g
the numerator in (4.1) is non-zero only to the extent that This
explains the Crank-Nicholson column in the table above.
Case 3. This example is the same as Case 2 but with non-uniform meshes, generated with
1=3. In this case we have no theoretical basis on which to expect convergence,
and indeed the table below suggests that convergence is not obtained.
6.7e-3 1.6e-3 2.7e-2 1.9e-2 9.8e-3
4.4e-3 1.6e-3 2.4e-2 2.0e-2 4.2e-1
Case 4. We conjecture that in some cases, even though the meshes t n may be non-uniform,
high-order convergence can still be achieved provided that the meshes are "locally
uniform" in the sense that for a given n, all of the h k 's ``near'' a given observation
point - i are constant. More precisely, there is a positive integer oe and a real ~ h for
which t n
-(i) is from (3.30).
In this example we repeat the experiment reported under Case 3, but took steps to
ensure that the refined meshes were "locally uniform" with
6.7e-3 5.8e-5 6.2e-3 2.2e-3 9.8e-3
4.4e-3 2.5e-5 2.1e-3 1.2e-3 4.2e-1
The table above suggests that high-order convergence can sometimes be obtained
on "locally uniform" meshes. However, we do not pursue this idea further.
Case 5. Next we present an example illustrating Theorem 6. The forcing term f in the
state equation is changed so that The RK4 method is used together
with the quadrature weight matrix W T rap arising from the trapezoidal quadrature
rule. These weights appear in Theorem 6, hypothesis (d). The mesh refinement is
uniform.
1.0e-2 6.5e-8
4.3
Case 6. Finally, we illustrate Corollary 2. We repeat experiment shown in Case 5 but
using Simpson's rule instead of the trapezoidal rule, and we compare g
1.0e-2 3.1e-3
2.2
5. Alternate Approaches. Throughout the preceding sections we have discussed the
use of the "discretized costate approximation" and have attempted to illustrate the delicacy
of that procedure when pointwise (in time) observations are involved. There are at least
three other possible approaches which also merit consideration when pointwise observations
are involved, which we now briefly discuss.
The strategy pursued in Sections 2 and 3 was to take the discretization of the adjoint
system; that is, to first derive the costate method for the continuous problem, and then
to discretize the resulting equations and integrals. An alternative to this would be to use
the adjoint of the discrete system. While this approach is straightforward when applied to
boundary value problems, certain complications may arise in the implementation for evolution
equations, particularly when time marching schemes are used for temporal discretization.
In particular, from the definition (3.1) of the discrete least-squares functional we have
From the expression for components of ffi p u n in (3.19) and the definition (3.11) of S n (q), we
see that
where fl n (q; p) is given by
Thus
The advantage of this approach is that it will yield exact values for the directional derivatives
of T n . Hence in numerical attempts to minimizeT n for any fixed n it would be ideal. However,
it involves the use of fl n (q; p) as given in (5.1), which may require significant amounts of
additional mathematical effort, code complexity and computational execution time. These
disadvantages may be negated in whole or in part for certain special cases if the expression
for can be simplified. Such simplification occurred with the Crank-Nicholson method
for uniform time meshes, as described in Case 2 of Section 4. However, we do not pursue
this approach further here.
Another alternative [3, Section V, (5.30)-(5.32)]) is to retain the strategy of discretizing
the continuous adjoint system, but to transform the costate equation (2.9) so as to increase
the smoothness of the source term. This can be done by the introduction of two new variables
/ and - given by
Z t[J(C   r)](s) ds
and
The idea is that since / is piecewise constant (as opposed to J(C   )r which is a linear combination
of Dirac delta functions), it should be easier to obtain accurate numerical solutions
of y by first approximating - as given here.
We briefly explored this idea numerically. We repeated the relevant numerical examples,
which are "Case 2" and "Case 3" of Section 4. In both cases, all four of the tms's were
used to approximate - instead of y, and then we set y In each of these runs,
we observed linear convergence of the directional derivatives; e.g., -
- was approximately 1
in each case. Obviously, for uniform time meshes (as in Case 2) this represents a loss of
accuracy, while for nonuniform meshes (Case 3) it is a definite improvement.
One other alternative is to consider modifying step (iii) of the Discretized Costate Approximation
along the following lines. Instead of attempting to integrate -
y over the whole
interval (0; t F ) in "one shot" with the tms, one should treat each subinterval (-
rately, using the fact that J n (C   is zero in the interior of these subintervals and that y
has jump discontinuities at - i . Correspondingly, the quadrature in (3.4) must be modified.
Although the implementation is more complicated, analysis along the lines of the proof of
Corollary 1 above then applies. Consequently, for sufficiently smooth state-variable forcing
functions f , full accuracy of the tms and the quadrature will yield full accuracy of the costate
approximations.
6. Summary and Conclusions. While the costate method offers potentially dramatic
time savings when computing gradients in parameter estimation problems for time-dependent
systems, we have indicated that care must be taken when using this approach. We have shown
both analytically and numerically that the use of convergent time-marching and quadrature
schemes in seemingly reasonable combinations can yield surprisingly poor results.
The difficulty lies in the fact that the forcing term of the costate equation may be very
non-smooth, so that the assumptions under which the tms and quadrature scheme converge
are violated. In the description of the "Discretized Costate Approximation" at the beginning
of Section 3, we assumed that the tms is used in step (iii) as a "black box" without regard
to the smoothness (or lack thereof) of J n (C
The costate implementation which we have analyzed involves discretization of the (con-
tinuous) costate equation using standard time marching schemes. Due to its simplicity, this
approach is very desirable when the resulting gradient approximation has the same order of
accuracy as the particular time marching scheme. We presented conditions which guarantee
this. Some of these conditions are straightforward. For instance, several popular second order
time marching schemes yield second order gradient approximations provided a uniform time
mesh is used in combination with (second order accurate) trapazoidal quadrature. Other
conditions are less obvious (c.f., Theorem 6). We also demonstrate that certain "reasonable"
implementations yield gradient approximations which fail to converge at all, e.g., fourth order
time marching on a uniform time mesh combined with (fourth order) trapa-
zoidal quadrature. Finally, several alternative (though less simple) costate implementations
were presented.



--R

"Computational issues in parameter estimation and feedback control problems for partial differential equations"
"Methods for the identification of material parameters in distributed models for flexible structures"
Estimation Techniques for Distributed Parameter Systems
"Numerical schemes for the estimation of functional parameters in distributed models for mixing mechanisms in lake and sea sediment cores"
"A spline-based technique for computing Riccati operator and feedback controls in regulator problems for delay equations"
"On nonconvergence of adjoint semigroups for control systems with delays"
de la Non-Linearite' D'Une Equation Parabolique Quasi- lineare
Numerical Methods for Unconstrained Optimization and Nonlinear Equations

Difference Methods for Initial Value Problems
--TR
costate method;evolution equations;parameter estimation