Numerical Approximation of an SQP-Type Method for Parameter
This paper deals with the numerical approximation of the Levenberg--Marquardt SQP (LMSQP) method for parameter identification problems, which has been presented and analyzed in [M. Burger and W. Mhlhuber, Inverse Problems, pp. 943--969]. It is shown that a Galerkin-type discretization leads to a convergent approximation and that the indefinite system arising from the Karush--Kuhn--Tucker (KKT) system is well-posed.In addition, we present a multilevel version of the Levenberg--Marquardt method and discuss the simultaneous solution of the discretized KKT system by preconditioned iteration methods for indefinite problems. From a discussion of the numerical effort we conclude that these approaches may lead to a considerable speed-up with respect to standard iterative regularization methods that eliminate the underlying state equation. The numerical efficiency of the LMSQP method is confirmed by numerical examples.
Introduction
Parameter identication denotes the procedure of determining unknown parameters appearing
in an underlying state equation (usually a partial dierential equation), from indirect
measurements related to the solution of this equation. Such problems frequently appear in
many applications, where mathematical models of physical, chemical, biological or economical
processes are used (cf. e.g. [1, 12, 16] and the references there).
Since such problems are ill-posed in general, i.e., the parameter to be reconstructed does
not depend on the observation in a stable way, regularization methods have to be used in
order to compute a stable approximation of the parameter in presence of data noise. Due to
the ill-posedness of the identication problem, the numerical approximation of such problems
is not a simple task. The standard approach that can be found in literature is based on an
a-priori elimination of the state equation, and an application of a discretized regularization
method to the resulting operator equation involving the parameter-to-output map, which is
the operator mapping the parameter to the corresponding observation. The main part in the
This work has been supported by the Austrian National Science Foundation FWF under project grants F
13/08 and F 13/09.

INTRODUCTION

evaluation of this map is the solution of the underlying state equation for given parameter,
which is numerically realized by standard discretizations such as nite elements.
This approach, in particular combined with iterative regularization methods (cf. [17] for an
overview), has been applied with success even to rather complicated parameter identication
problems (cf. e.g. [9, 24, 25]). However, since this methods need a high number of direct solves
(i.e., solutions of the state equation), ne discretizations of the parameter yield a considerable
computational eort, which results in high CPU-times or even in the impossibility to use ne
discretizations. Another drawback of this approach is that the discretizations of state and
parameter are rather independent, which makes the numerical analysis extremely di-cult.
Therefore, fundamentally dierent methods for the solution of parameter identication
problems have been investigated recently, whose common idea is to avoid the a-priori elimination
of the state equation (cf. [10, 20, 26]). The aim of this paper is to discuss the numerical
approximation of an iterative regularization method based on the idea of sequential quadratic
programming (cf. [10]). We investigate Galerkin-type discretizations in the product space
for parameter, state variable and a corresponding Lagrangian variables, which leads to a sequence
of well-posed indenite systems. With this approach we are able to show convergence
of the numerical approximation both for the quadratic programming problem arising in each
iteration step and for the overall minimization procedure.
The general setup in this paper is as follows: we assume that we are given a noisy measurement
z - satisfying
where the exact data satisfy
z := E^u; (1.2)
Our aim is to identify the parameter q 2 Q ad  Q (where
Q ad is a closed subset of Q with non empty interior) in the underlying equation
is a continuous nonlinear operator with
In this setup X, X  , Q and Z are Hilbert spaces, and X  can be identied with the dual of
X. Finally, we assume that e is continuously Frechet-dierentiable on X  Q and that the
partial derivative e u 2 L(X; X  ) is self-adjoint and satises the coercivity condition
he u (u; q)v; vi   e kvk 2
for some  e
The above setup is typical for a partial dierential equation of elliptic type, which is
also the main type of application we have in mind. We want to mention that the innite-
dimensional analysis carried out in the preceding paper [10] was not restricted to elliptic
problems, but only assumed well-posedness of the state equation for given parameter. How-
ever, since the numerical approximation techniques for elliptic problems dier from the ones
for parabolic or hyperbolic problems (cf. e.g. [32] for an overview), one cannot expect a successful
unied approach to corresponding parameter identication problems. For this reason
we start with an investigation of the elliptic case in this paper, but we want to mention that
the numerical identication of parameters in transient equations or even mixed systems of
equations is an important and challenging problem for future research.
In [10], it has been mentioned that the parameter identication problem in the above setup
is an ill-posed inverse problem and we have proposed the following iterative regularization
method based on the idea of sequential quadratic programming:
Programming Method).
Let
be a given initial value and let ( k ) k2N be a bounded sequence of positive real numbers.
The Levenberg-Marquardt sequential quadratic programming (LMSQP) method consists of
the iteration procedure
where is the minimizer of the quadratic programming problem2
subject to the linear constraint
The iteration procedure is stopped as soon as
with appropriately chosen  > 1.
Due to the results of [10], the LMSQP-method is a convergent regularization method,
in particular the quadratic programming problems of the form (1.7), (1.8), which have to
be solved in each iteration step, are well-posed. Our aim in this paper is to investigate
the numerical approximation of the LMSQP-method by a Galerkin-type approach. We shall
show below that this leads to an indenite system in each iteration step, whose solution is an
approximation of optimal order to the solution of (1.7), (1.8). Moreover, we show that the
reconstructions obtained with the discretized LMSQP-method converge to a solution of the
parameter identication problem as the noise level and the discretization size tend to zero, if
an appropriate stopping rule is used, which relates the residual to the noise level and some
measures for the discretization.
Moreover, we shall discuss the solution of the discretized Karush-Kuhn Tucker system,
which is an indenite linear system to be solved for the discretized equivalents of state, parameter
and Lagrangian variable. The standard approaches to the solution of such discretized
problems arising from partial dierential equations are reduced SQP-methods, where state and
Lagrangian variable are eliminated a-priorily. We recall the basic properties of the reduced
SQP-approach, but we mainly focus on the iterative solution of the whole system with appropriate
preconditioning. This promising approach has been employed recently for parameter
identication (cf. [20, 26]) and optimal control problems (cf. [2, 3, 4, 5]) with good numerical
results, in particular with respect to e-ciency.
The paper is organized as follows: in Section 2 we investigate the numerical approximation
of the LMSQP-method by a Galerkin-type approach and discuss the well-posedness,
stability and approximation properties of the discretized Karush-Kuhn-Tucker (KKT) system;
the convergence of the discretized solutions is shown in Section 3. Some further numerical
methods and the implementation of the outer iteration, i.e., the SQP-iteration under the
assumption that we are able to solve the quadratic optimization problems arising in each
step of the LMSQP method, are examined in Section 4. We discuss the correct scaling of
variables, globalization strategies as well as a multi-level approach, which leads to a further
speed-up of the method. Section 5 is devoted to the inner iteration, i.e., the numerical solution
of the discretized KKT-system. Some basic properties of this symmetric indenite
problem are studied, as well as its iterative solution with appropriate preconditioning. As a
rst application we investigate the identication of a potential in an elliptic boundary-value
problem, where we can give quantitative error estimates in terms of the discretization sizes.
Some numerical experiments related to this identication problems are presented in Section
7, before we nally conclude and give an outlook to further interesting problems related to
this topic in Section 8.
Discretization Techniques
In the following we investigate the discretization of the LMSQP-method by a Galerkin ap-
proach. First of all, we assume that we have discretized data z -; 2 Z   Z of the form
z
is the orthogonal projector onto the nite-dimensional subspace Z  . Note
that we can give an error estimate for z -; using (1.1) and kR
Now let X h  X, Q h  Q be nite-dimensional subspaces of X and Q, with the corresponding
orthogonal projectors
. Then we can discretize the
LMSQP-Method as follows:
be as above and let
be a given initial value. Moreover, let ( k ) k2N be a bounded sequence of positive real numbers.
The Galerkin Levenberg-Marquardt sequential quadratic programming (GLMSQP) method
consists of the iteration procedure
where is the minimizer of the quadratic programming problem2 kR  (Eu z - )k 2
subject to the linear constraint
2.1 The Discretized Karush-Kuhn-Tucker System 5
Note that the constraint (2.5) can be rewritten in operator form as
~
to be solved for (u; with the notation
and P
is the adjoint of P h . Under the assumption (1.5), we obtain that
for all v 2 X h , i.e., the discrete bilinear form associated with the operator P
coercive
on X h . This implies by the Lax-Milgram theorem, that (2.6) is uniquely solvable with respect
to u for given q 2 Q h . Consequently, in an analogous way to the proof of Proposition 2.1 in
[10] we may show the following result on the well-posedness of the quadratic programming
problem that has to be solved in each step of Method 2:
Proposition 2.1. Let e be continuously Frechet-dierentiable, let (1.5) hold and let  k > 0.
Then the quadratic programming problem (2.4), (2.5) has a unique solution
which is also the only local minimum.
2.1 The Discretized Karush-Kuhn-Tucker System
In [10], the Karush-Kuhn-Tucker system for the innite-dimensional version of the LMSQP-
method has been derived and analyzed in the framework of linear saddle point problems.
Now we will discuss the discretized analogue of this system, namely the rst-order optimality
conditions for the quadratic programming problem (2.4), (2.5).
The Lagrangian of (2.4), (2.5) is given by
for
are equal to the identity on X h and Q h ,
respectively, we can rewrite the Lagrangian as
~
with the operators K k and L k dened by (2.7), (2.8). The KKT-system can now be deduced
by computing the partial derivatives of the Lagrangian with respect to u, q and , i.e.,
solves the linear saddle-point problem@ P
~
~
~
~
As in [10], we dene the symmetric bilinear form a
a
and the bilinear form b
Moreover, we use the right-hand sides
Then the KKT-system (2.12) can be interpreted as the Galerkin approximation of an indenite
variational problem, i.e., (u; q; is the solution of
a
In an analogous way to the proof of Theorem 2.3 in [10] we can show that the bilinear
form a satises the kernel-ellipticity condition on X h Q h , i.e., there exists a constant  a > 0
such that
a
and that b satises the LBB-condition
sup
for some  b > 0. This implies the following well-posedness result (cf. [7, 8]) for the discretized
problem (2.17), (2.18):
Theorem 2.2. Let e be continuously Frechet-dierentiable, let (1.5) hold and let  k > 0.
Then the indenite system (2.17), (2.18) has a unique solution (u; q;
which depends continuously on the right-hand sides f k and g
k .
Since the constants  a and  b are the same as in the corresponding innite-dimensional
conditions in X  Q, they are in particular independent of the discrete subspaces X h and
This allows to deduce an approximation result for the solutions of (2.17), (2.18) to the
solution of the innite-dimensional KKT-System, given in variational
form as
a k (u; q;
with a k given by
as above and g k dened by
2.1 The Discretized Karush-Kuhn-Tucker System 7
Theorem 2.3. Suppose that the assumptions of Theorem 2.2 are satised and let
denote the unique solution of (2.17), (2.18). Then there exists a constant c > 0 independent
of X h and Q h such that
where (u; q; ) denotes the unique solution of (2.19), (2.20) and
Proof. First, let (~u h ; ~
the solution of (2.17), (2.18) with a
k replaced by a k ,
k . Then Theorem 2.1 in [8] implies the existence of a constant c 1 > 0 (independent of X h
and Q h ) such that
Moreover, the stable dependence of the solutions of (2.17), (2.18) on the right-hand side
implies the existence of c 2 > 0 with

sup
hg
ja

sup
(R
(R
and with the triangle inequality we may conclude (2.23).
Theorem 2.3 provides an error estimate for the solutions of the discretized saddle-point
problem (2.17), (2.18), consisting of two parts corresponding to the numerical approximation
in the image space Z and in the pre-image spaces X and Q. An obvious estimate for the rst
term is
;h  inf
~
which possibly does not lead to a quantitative estimate, since there is no additional information
on the smoothness of the noisy data. An alternative estimate is
~
The inmum of can usually be estimated more easily, since the exact data z are
smoother due to the fact that ^ u is the solution of the state equation for some parameter ^
q.
E.g., if the state equation is of elliptic type with solution ^
is the embedding operator, and R  results from a standard nite element discretization on a
grid with neness , then we have at least
Another important observation is that the last term vanishes if the discrete spaces Z  and
are equal, which can be achieved in some applications.
The second term in (2.23) shows that the Galerkin approximation of the KKT-system is
of optimal order in X h Q h X h ; it can be estimated by standard methods for nite element
discretizations; quantitative estimates can be obtained using the regularity of the iterates.
This part depends of course strongly on the specic application.
3 Convergence Analysis
In this section we will analyze the Galerkin LMSQP method with respect to convergence, i.e.,
the convergence of the reconstruction obtained with an appropriate stopping rule as the noise
level and the measure for the discretization neness tend to zero. With
identify the innite-dimensional case, i.e., We assume that
the discrete subspaces satisfy
If we denote by e k and f k the error terms
we can rewrite the Karush-Kuhn-Tucker system (2.12) as@ P
~
~
~
~
~
~
r kA (3.2)
where the r k denotes the remainder
As in [10], we require a condition on the nonlinearity, which is summarized in the following:
Assumption 1. Let (1.5) be satised and dene the remainder r(u; q) by
Then we assume that there exists a constant
kEe
and that there exists a solution (^u;
of the parameter identication
problem.
If we dene the discretization measures  h ,  h by
where
kEe
and   by
Then for all (u;
holds, where
Remark 1. If X h , Z  and Q h are standard nite-element spaces on some triangulations, then
h ,   and  h can be estimated by the approximation error of these elements. In particular,
if the discretization parameter (i.e., the maximal size of a triangle) tends to zero and if the
triangulation is regular, one can guarantee that  h ,   and  h tend to zero (cf. [32] for further
details).
For the choice of the stopping index we use a numerical version of (1.9), which involves
the discretization measures dened above:
For an appropriate choice of  , this allows us to prove the following monotonicity property of
the iterates:
Lemma 3.1. Let Assumption 1 be fullled, let the noise be bounded by (1.1), and assume
that
In addition,  k is chosen such that  k   k 1 for all k 2 N and that
s
and the stopping index k  is chosen according to the generalized discrepancy principle (3.11)
with
then and the estimates
and
hold for all k < k  .
Proof. Assume that q
we deduce the identity
The noise bound (1.1) implies that
and using the Cauchy-Schwarz inequality together with (3.9) we obtain the estimate
(3.16) follows from dividing (3.15) by  k and the fact that
s
By induction we can now show that q k 2 B  (q 0 ) for k < k  and  satisfying (3.14).
In an analogous way to the proof of Lemma 3.2 in [10] we can prove the following statement
on the niteness of the stopping index k  if - > 0:
Lemma 3.2. Under the assumptions of Lemma 3.1, the discrepancy principle (3.11) yields
a nite stopping index k  if
and  is chosen according to (3.14).
One observes that in the above estimates, the term - ;h now plays the same role as the noise
level - in the innite-dimensional setup. Therefore it is also possible to prove convergence
as in the same way as convergence in the innite-dimensional case for - ! 0 (cf.
[10, Theorem 3.5]). Consequently, we do not give the detailed convergence proof, but refer to
[10] for further details on the technique of the proof. We only recall the basic assumptions
on e and give the nal convergence result, where we use the notation (u -;;h
k ) for the
iteration according to (2.12) with initial value (P h
discretization
parameter h and .
Assumption 2. In addition to Assumption 1, assume that e is of the form
with continuously Frechet-dierentiable (nonlinear) operators A
, such that
Moreover, we assume that A and N satisfy the nonlinearity conditions
kEe
and
kEe
for some positive constants
2 and
3 .
Theorem 3.3 (Convergence). Let Assumption 2 and (3.12) be fullled with ,  suciently
small, and let the noise be bounded by (1.1). Moreover, let  k be chosen such that
and that (3.13) is satised. If the perturbed iteration is stopped with
according to the generalized discrepancy principle (3.11) with
(uniformly bounded in h and ) satisfying (3.14), then
(q -;;h
where (u; q) is a solution of (1.3) with
Proof. Analogous to the proof of Theorem 3.5 in [10].
4 Numerical Realization of the SQP-Iteration
In the following we want to discuss some numerical methods and variants for the 'outer
iteration', i.e., the Galerkin LMSQP algorithm under the assumption that we are able to
solve the discretized KKT-system numerically. The 'inner iteration', namely the numerical
solution of the indenite system (2.12) will be investigated in Section 5.
4.1 Scaling of State Variable, Parameter and Lagrangian Variable
The performance of an iteration algorithm often depends crucially on the way the problem
is formulated. Scaling is a well-known technique for reformulating an optimization problem
whose main objective is twofold: on the one hand all the variables should be of similar
magnitude, on the other hand also the value of the derivatives should all be of similar size. In
unconstrained optimization, a problem should be rescaled in such a way, that changes of the
iterate in one direction do not result in by far larger changes of the value of the objective than
changes in another direction. In constrained optimization the above statements are also true
for each constraint. Additionally the set of constraints should be well balanced with respect
to each other such that each constraint has equal weight. Furthermore the set of constraints
should be balanced with respect to the objective. As scaling is of high practical importance
for any optimization problem, many aspects can be found in monographs on optimization (cf.
e.g. [19, 30]).
We want to consider only the last aspect in this context, i.e., the scaling of the state
constraint with respect to the objective which is also of high importance for achieving fast
convergence of the outer iteration. For the inner iteration, the aspect of scaling can be
included in the construction of a good preconditioner. The outer iteration of an SQP method
tries to attain two goals at the same time: feasibility of the iterate with respect to the state
constraint and optimality of the iterate with respect to the objective. One aspect dominating
the other results usually in bad convergence properties: If the feasibility aspect dominates,
only very small changes of the iterate are possible in order to ensure "almost" feasibility. If
the optimality aspect dominates, any violation of the state constraint is reduced too slowly.
For the LMSQP method in the form of (2.4),(2.5) it turned out that in many situations
the feasibility aspect is strongly dominating. Using line search methods for globalization (see
also Subsection 4.2) this results usually in step lengths much smaller than one. Replacing the
state constraint by a preconditioned state constraint leads to a better balanced formulation
and to much faster convergence. Furthermore a step length parameter equal to one is accepted
in almost all steps. Another aspect of this kind of rescaling is treated in Subsection 6.2.
4.2 Globalization Strategies
The LMSQP method is a variant of Newton's method and therefore only locally convergent
(see also the analysis in Section 3). For this reason, globalization strategies, such as trust
region methods or line search strategies (which are the two most popular classes of globalization
techniques in optimization), are needed.
The basic idea of trust region methods is to add an additional constraint on the maximal
increment to the quadratic optimization problem for the correction step of the current iterate,
i.e. instead of (2.4), (2.5) one would solve (2.4), (2.5) and k(u u k ; q q k )k   k with  k chosen
appropriately. Trust region methods have been successfully applied to PDE constrained
optimization problems (see e.g. [14, 38]), often using a reduced SQP approach. We want to
mention that a similar eect as with trust-region methods could be reached in principle by
controlling the penalty parameter  k , which also restricts the step size and produced good
numerical results (see Example 7.1). For a comprehensive overview of trust region methods
we refer to Conn et. al. [13].
In the code used for two-dimensional problems (cf. [10] and Example 7.2), we use a line
search algorithm for globalization. In contrast to trust region methods, the calculation of the
increment is split into two phases: rst of all, a search direction is determined, and secondly
the estimation of a step length parameter indicating how far into the search direction one
should go. For the computation of the search direction we solve the optimization problem
and (2.5). In order to determine the step length we cannot use the objective itself as a
criterion (as in unconstrained optimization), but have to use a merit function which balances
the minimization of the objective with the feasibility with respect to the state constraint.
Applied to a discretized optimization problem of the form
~
subject to an equation constraint of the form
possible choices are the l 1 -merit function
4.3 Nested Multi-Level Optimization Techniques 13
and its variants, and the augmented Lagrangian
where  is an estimate of the Lagrangian variable corresponding to the discretized equation
constraint. Both merit functions are exact in the sense that for  su-ciently large, minimizers
of the original constrained optimization problem also minimize the merit function.
A crucial property in the design of a merit function is that it should accept step length
one close to a solution in order to preserve the quadratic convergence of the SQP method.
The augmented Lagrangian works well, as long as the estimate for the Lagrangian multiplier
is accurate enough, whereas the l 1 -merit function sometimes suers from the so-called
"Marathos-eect", i.e. it does not accept unit step length and therefore causes a slow-down
of the convergence. A strategy to overcome this di-culty using a second order correction
can be found in [30]; nevertheless, it performed very well in our numerical experiments (see
Example 7.2).
4.3 Nested Multi-Level Optimization Techniques
Important tools for the e-cient numerical approximation of innite-dimensional optimization
problems are multi-level optimization methods. In the nested multi-level setup, one starts the
optimization procedure at a coarse level X
, where the iteration procedure can be
carried out e-ciently. If an appropriate stopping rule is satised, one interpolates the state
and parameter obtained in this way to a ner level X h 2
(for serving now
as a starting value on this level. This procedure is repeated until the nest level is reached.
Usually, nested space are used in this approach, i.e., X h 1
(for
which leads to simple interpolation operators. Since one cannot choose the discretization of
the data arbitrarily in general , we consider only the case of xed  here, but a multi-level
approach in  can be realized in an analogous way, if necessary.
Nested multi-level methods outperform standard discretization techniques in many cases
(cf. e.g. [21, 22, 29]); usually a considerable number of iterations is needed on the coarse level
only , where the numerical eort per iteration is very low. On the nest levels, the stopping
rule is often satised already after one iteration step and so the overall eort is less than for a
direct discretization on the nest level. For the Galerkin LMSQP method, we can formulate
a multi-level algorithm as follows:
Algorithm 4.1 (Nested Multi-Level Galerkin LMSQP). Given a decreasing sequence
'=1;:::;L with nested spaces X h '  X h '+1 , Q h '  Q h '+1 (e.g. h non-increasing
sequence  ' satisfying (3.14), the nested multi-level Galerkin LMSQP method consists
of the following iterative procedure:
1.
2. Perform the Galerkin LMSQP method until the stopping criterion (3.11) is satised
with stopping index k  (').
3. If stop the iteration, else prolongate the iterate (u '
k ) to the ner level
which results in a new starting value (u '+1
and go to step 2.
The analysis in Section 3 shows that for  '
, the estimate
holds, where  ' is the error corresponding to the interpolation of the iterates from level ' 1
to level ', i.e.,
kR  Ee '
Z
kf '
This monotonicity estimate corresponds very well to the intuition that only few iterations are
needed on the ne levels, in particular if  '
k is decreasing, which leads to
kR  Ee '
Z   '
For a ne level with small , we can expect that
and the second term  '
can be expected to be negligible. I.e., the stopping
rule at level ' is probably satised with k
Under typical conditions, where X h ' and Q h ' correspond to standard nite-element spaces
on dierent renement levels of an initial triangulation of a
domain
one can show that at
least  consequently
ch
ch 1
for some constant c 2 R+ , where
Together with the above estimate one can show with a standard proof technique that the
converges to a solution (u; q) of the parameter identication problem for
5 Numerical Solution of the KKT-System
In the following we will discuss the numerical solution of the discretized KKT-system (2.12)
for xed iteration number k. We have seen above that the Galerkin-type approximation (2.12)
of the original KKT-system is stable and convergent, now we discuss some of its structural
properties, which are important for the application of iterative solution methods and for the
construction of preconditioners.
5.1 The System Matrix M 15
Choosing bases
of the nite-dimensional subspaces X h and Q h , we may represent
via
with coordinate vectors V;  In order to transform (2.12) into a linear
system for V , S and , we dene the matrices
and the vectors
This allows us to rewrite the discretized KKT-system (with penalty parameter
respectively as
where M is the matrix in (5.6) and
The structural properties of M and its sub-matrices will be examined in the following section.
5.1 The System Matrix M
Due to the well-posedness result on the discretized KKT-system (2.12) (cf. Theorem 2.2), we
may conclude that the system matrix M is regular. In order to obtain further insight into
the structure of M , we investigate the properties of the sub-matrices G, H, K and L:
Proposition 5.1. The matrices K 2 R mm and H 2 R nn are symmetric positive denite,
and the matrix G 2 R mm is symmetric positive semi-denite. If, in addition, the operator
is injective on X h , then G is regular, too.
Proof. Let u and q be as in (5.2), then there exist constants c 1 (h) and c 2 (h) such that
where j:j denotes the euclidean norm in R n and R m , respectively. Thus, we have
and
Moreover, the identity
implies that G is positive semi-denite and regular under the assumption that E is injective
on X h . The symmetry of the matrices G, H and K can be veried in a similar way, using
the symmetry of scalar products and the self-adjointness of the operator K k .
The matrix L 2 R mn is di-cult to analyze, it is neither symmetric nor regular in general
(in particular if n 6= m). However, some fundamental properties of M (such as its regularity)
rely rather on G, H and K than on L. Moreover, the classical splitting of a symmetric
saddle-point problem as@ G 0 K T
I
I
where H  := H and C is the Schur-complement
is possible if we only know that G and H are regular. In particular, we may conclude that
M has n +m positive and m negative eigenvalues.
5.2 Reduced SQP Approaches
The basic idea of reduced SQP-methods is the a-priori elimination of the equality constraint,
which can be written in matrix form as
which is equivalent to an elimination of V and  in (5.6).
Due to Proposition 5.1, K is a regular, symmetric matrix and thus, we may compute
which yields after some calculations the n  n-system
with
The reduced SQP-approach seems of particular interest if n  m, which is a frequently
used discretization strategy for parameter identication and optimal control problems (cf. e.g.
[35, 36, 37]). The original matrix M is an indenite matrix of size (2m+n) (2m+n), while
5.3 Simultaneous Solution of the KKT-System 17
the reduced system matrix M r in (5.12) is of size n  n. However, M r is not a sparse matrix
even if all the sub-matrices of M are sparse, since it involves the inverse of K. Moreover,
the evaluation of M r is more expensive than the evaluation of the original system matrix M ,
since it involves the solution of two systems of the form
with dierent right-hand sides g, while for the evaluation of M only direct evaluations of K are
needed, which are very cheap for typical nite element discretization of the state constraint.
In practice, one usually tries to compensate this disadvantage of reduced SQP-methods by
using a Broyden-type update for the reduced system matrix instead of the exact matrix M r ,
which leads to e-cient optimization algorithms for small n.
5.3 Simultaneous Solution of the KKT-System
Recently, the simultaneous solution of KKT-systems by iterative methods has been investi-
gated, in particular in connection with optimal control problems (cf. [2, 4, 5, 20]). Compared
to the reduced SQP-approach, a simultaneous solution strategy has the obvious advantage
that the allocation and evaluation of the system matrix M is much cheaper than of M r . The
pay-o is that M is indenite and larger than M r , which might cause additional eort. How-
ever, the main eort in the reduced SQP-approach is related to the evaluation or assembly of
the system matrix M r , respectively, and therefore a simultaneous solution of the KKT-system
can result in a tremendous speed-up of the SQP-method, in particular for ne discretizations.
At a rst glance, it seems rather straight-forward to solve (5.7) by a standard iterative
method for indenite systems such as inexact Uzawa methods (cf. [6, 15]) or Krylov-subspace
methods such as GMRES (cf. [34]), MINRES (cf. [31]) and QMR (cf. [18]). However, in
the case of large-scale problems, we have to expect a large condition number (note that
is usually small and that M is singular for and a complicated eigenvalue pattern of
the matrix M , which might cause iterative methods to diverge or to need a high number of
iterations. Therefore, an appropriate preconditioning technique seems necessary for any of
the methods. We do not go into details here, but refer to the forthcoming paper [11] for a
discussion of preconditioners.
In the following we distinguish two types of solvers that seem appropriate for the solution
of the indenite system (5.7) and discuss their basic properties with respect to the special
structure of M .
Inexact Uzawa Iterations
Inexact Uzawa methods and similar iteration procedures have been developed for the solution
of the classical Stokes system and similar problems (cf. [32] for an overview). The classical
Uzawa method is just a gradient method for the dual of the corresponding Lagrange functional,
the inexact Uzawa method can be interpreted as a preconditioned version (cf. [32]). Following
the exposition by Zulehner [39], we can write an inexact Uzawa method for a system of the
form (5.6) as
A
followed by
A is a preconditioner for the diagonal matrix
A :=
C is a preconditioner for the Schur-complement C dened by (5.8). In terms of (5.7) we
can write the inexact Uzawa iteration as
M is a preconditioner for the system matrix, given by
.
A convergence analysis of this method is available only in the case when A is a regular
matrix (cf. [6, 39]), which means that we have to assume that G is regular. The latter
is true e.g. if the data z represent distributed data for the state, i.e., E is an embedding
operator. In this case, the structure of A is rather simple and it is not a di-cult task to
construct a preconditioner, even exact preconditioning seems possible (note that G is just a
mass matrix for a typical nite element discretization). Since the matrices G and H do not
change during the SQP-iteration we may even compute decompositions in a preprocessing
step. The construction of a preconditioner for the Schur-complement C is more di-cult and
must take into account the specic nature of the underlying state equation.
Krylov-Subspace Methods
The Krylov-subspace methods GMRES and QMR are variants of the CG-algorithm that are
applicable to indenite problems, too. The basic idea of such methods is a defect minimization
in the Krylov-subspace
generated by X 1 , in the k-th iteration step. Since preconditioned CG-methods are probably
the most successful class of iteration methods for positive denite systems, such methods
seem very attractive also in the indenite case, although additional di-culties may arise (cf.
e.g. [34]).
The convergence analysis in [34] and [18] shows that the error bounds obtained for both
methods are essentially the same, and mainly dependent on the eigenvalue distribution and
the condition number of the system matrix M . Therefore, appropriate preconditioning is
again of high importance, in this case also with the possibility that G is singular. We refer
to [11] for a detailed discussion of this problem.
As a rst application we investigate the identication of the potential q in the elliptic boundary
value problem
in
@
from a state observation in L
2 which is a well-studied problem in literature (cf. e.g. [33]).
In [10], it has been shown that in the setup (d denotes the space dimension)
the operators
satisfy all assumptions needed for the convergence analysis of the LMSQP-method. Now we
shall study a concrete nite-element discretization of the KKT-system and the derivation of
estimates for the numerical errors   ,  h and  h .
6.1 Error Estimates for the Discretized KKT-System
In this case we can write the whole KKT-system in classical form as
in
in
in
again with homogenous Dirichlet boundary conditions upon v and  on @ where L d is a
dimension-dependent dierential operator of order 2d corresponding the norm in H
7 e.g.,
we have
supplemented by homogenous boundary conditions up to order d 1. If f 2 L
and
3 a standard elliptic regularity argument shows that ^
0(for all k 2 N. In the same way we can show that  k 2 H
and s k 2 H
7 This
additional regularity can be employed to derive standard error estimates for nite-element
discretizations of the KKT-system (2.12).
If we use piecewise linear nite elements on regular triangulations T  and T h for the
discretization spaces Z  and X h , where  and h represent the neness of the grids, then a
classical approximation result for nite elements (cf. [32, p.96]) implies that
Of course, one could also use piecewise constant elements on T  , which would yield
However, in practical applications a higher-order approximation in  is often desirable, since
can be signicantly larger than a reasonable choice of h. A canonical approximation of the
parameter q is a nite element space of order greater or equal d on a regular triangulation T ~ h
assumptions on the exact solution ^
q one can obtain quantitative estimates
for  h in terms ~
h. At a rst glance it seems surprising that one needs a-priori assumptions
on the parameter, but not on the state in order to derive error estimates. However, due to
the ill-posedness of the identication problem with respect to the parameter q, such a-priori
knowledge seems to be necessary. The approximation of the state corresponds rather to the
approximation of the underlying elliptic state equation, which is well-posed with respect to
u and yields further regularity. We nally want to mention that according to the theory
developed above, one could choose T ~ h independent of T h , but this would cause unnecessary
complications in the implementation of the method.
We note that alternatively one can use the space
for d  3, which yields
An appropriate discretization strategy is e.g. to choose Q h as the space of piecewise constant
elements on an underlying grid T ~ h . The advantage of this approach is that elements of order
greater than one, which are necessary for
(d  2), can be avoided.
6.2 Structure of the System Matrix
For the potential identication problem, some parts of the system matrix M in (5.6) are well-
understood. First of all, G is an L 2 -mass matrix and it is positive denite if the triangulations
T  and T h coincide, which we assume in the following. The eigenvalues of G are then all of
order h d . The matrix H is the stiness matrix for the dierential operator L d , with minimal
eigenvalue of order h d and maximal eigenvalue of order h d .
The matrix K is the sum of a stiness matrix for the Laplacian and a weighted mass
matrix (with weight q k in the L 2 -scalar product), where one can expect the rst part in
this sum to be dominating. Thus, the stiness matrix ^
K for the Laplacian will be a good
preconditioner for K. The maximal and minimal eigenvalues of K and
K are of order h d 2
and h d , respectively. The remaining part in the system matrix, namely the matrix L, is
di-cult to understand, since its elements are weighted L 2 -scalar products of basis functions
of dierent nite element spaces. However, the spectral norm of L can be estimated, it is of
order ~ h d .
The construction of preconditioners for G and H is well-investigated, even exact preconditioning
seems to be applicable. For K it seems reasonable to use a preconditioner ^
K for the
Laplacian, e.g. a multi-grid preconditioner. With preconditioning for K, the system matrix
can be transformed to
~
with the corresponding Schur-complement
~
K is an appropriate preconditioner for K, then we can estimate the minimal eigenvalue by
min ( ~
C)   min
and the maximal eigenvalue by
O
Hence, the condition number of ~
C is independent of h, but only depends on  and ~ h
h . One
observes that the condition number is decreasing as ~ h tends to h from above (note that
usually ~ h  h). For the Uzawa iteration, one can choose the preconditioner ^
C in this case as
a multiple of ^
K 1 or even of G 1 . If ~ h  h, the Uzawa iteration seems not to
be optimal, in this case one can apply either a reduced SQP-approach or use Krylov-subspace
methods with dierent preconditioning strategies. For the details on the latter we refer to
[11].
7 Numerical Experiments
In order to test our theoretical results, we numerically solve some model problems, which have
already been investigated with respect to the convergence behavior of the LMSQP-method in
[10].
Example 7.1. Our rst example is the identication of the potential q in (6.1), (6.2) from a
state observation u 2 L
2(with
The exact potential is given by
which is an element of
3 This problem was implemented in the software-system
MATLAB.
The data are generated by solving the state equation on a ne grid and subsequent interpolation
to a coarser grid; the noise is an additive high-frequency perturbation. We used
uniform grids with m nodes for the discretization of the state u and the Lagrange-parameter
and n nodes for the parameter q, i.e., . The parameters
k are chosen according to  which lead to convergence of the
method even for starting value q  0.
The KKT-system (5.6) is solved by the QMR method, using an Uzawa-type preconditioner
as described in Section 6.2, with ^
K a preconditioner for the Laplacian and ^
. The convergence results for the overall LMSQP-method have been shown in [10] and
compared to a Levenberg-Marquardt method following the feasible path. It turned out that
both methods lead to almost the same iteration sequence q k . In particular, the number of
iterations needed until the stopping rule is satised, is the same for both methods. Now
we compare the numerical e-ciency of the LMSQP-method with feasible path approaches,
namely the Levenberg-Marquardt method (LM) on the feasible path (with the same Galerkin
discretization as for LMSQP and solution of the Gauss-Newton system by a preconditioned
CG-method) and a Broyden-type variant of the Levenberg-Marquardt method (cf. [23] for
further details).
22 7 NUMERICAL EXPERIMENTS

Table

1: CPU-time (in seconds) needed for the LMSQP-method, the LM-method and a
Broyden-type variant of the LM-method.
For this sake we choose dierent discretization levels (xed during the iteration) and
measure the CPU-time needed for the LMSQP-method, until the stopping rule is satised (for
xed noise level -). From the results shown in Table 1 one observes that the LMSQP-method
with simultaneous solution of the KKT-system outperforms the feasible-path approaches for
all dierent discretizations. Since the LMSQP and the LM-method need the same number of
outer iterations, the dierence in the numerical eort is caused by the fact that the eort for
the evaluation of the system matrix in the LM-method is signicantly higher than evaluation
and preconditioning of the system matrix in the simultaneous LMSQP-method. Obviously,
the gain in the numerical eort for the evaluation of the system matrix increases with the
number of discretization points, which explains the extremely large CPU-time for the LM-
method at the nest discretization level
is much faster than the LM-method, which is again caused by the fact that the evaluation of
the system matrix can be carried out e-ciently. However, the number of iterations needed
for the Broyden-type variant is much larger than for the other two methods, which use the
full information about the derivatives.
Finally, we investigate the spectral condition of the system matrix M as well as of the
matrix ~
M dened by (6.13), where we use a preconditioner for the Laplacian as ^
K. >From the
left picture in Figure 1, which shows the condition number as a function of the discretization
size h (in logarithmic scale) for xed one observes that the condition number of M
grows quadratically with h 1 , while the condition number of ~
M is much smaller and almost
independent of h. The second part of Figure 1 shows a plot of the condition numbers vs.
the parameter  in doubly logarithmic scale, from which it seems that the growth of the
condition number as  ! 0 is slower for ~
M than for the original matrix M . In both cases, the
condition number seems to be a convex function of , which has a unique minimum at some
. However, this value  is rather large and values of  that are signicantly larger than
are not of interest for our purpose, since they would cause a tremendous slow-down of the
outer iteration. Therefore we can focus our attention to the case  < , where the condition
number increases in a monotonically with  1 .
Example 7.2. Our second numerical example is the identication of the conductivity q 2
log(cond(M))
Condition Number vs. Discretization Size
Original
Preconditioned State
log b
log(cond(M))
Condition Number vs. b
Original
Preconditioned State

Figure

1: Plot of the spectral condition of the matrix M vs. the discretization size h (in
logarithmic scale, left) and vs. the parameter  (in doubly logarithmic scale, right). The
solid line shows the condition number of the original matrix M , the dashed line of the matrix
~
M with preconditioned state equation.
in
in
on
@
from a state observation u 2 L
2 The
domain
is a ball in R 2 with missing rst quadrant,
i.e., in radial
coordinates
The exact parameter to be reconstructed is ^
q  1, the right-hand side in (7.1) is given by
r
with
The corresponding solution of the state equation is ^
3r). The data are generated
using the exact solution ^ u perturbed by uniformly distributed random noise. For the discretization
we used triangular nite elements with piecewise quadratic shape functions for
the state u and the Lagrange parameter  and piecewise constant shape functions for the
parameter q. The results were calculated using the nite element code FEPP [27], developed
at the Department for Analysis and Computational Mathematics of the University of Linz.
We want to mention that this identication problem is quite challenging not only due to
the complicated geometry, but also due to the fact that q is not identiable along a level
line in the interior, where u attains an extremum. This does not destroy the theoretical
identiability results, because it is a set of Lebesgue-measure zero, but it can be expected to
create numerical di-culties.
Results for exact data can be found in Table 2. The good performance of the method
with respect to both, CPU time and number of outer iterations can be observed clearly.
Especially for problems with ne discretizations of the parameter q, this method can still be
realized e-ciently, while classical approaches do not yield results in reasonable time. A plot
Level dim q dim u avg QMR it SQP it time

Table

2: CPU-time and number of inner (QMR) and outer (SQP) iterations for exact data

Figure

2: Parameter distribution for exact data at level 4
of the parameter q can be found in Figure 2, from which one observes that the parameter is
reconstructed very well except in a neighborhood of the level curve 0g.
Additional speed-up can be gained using a multi-level approach as described in Subsection
4.3. We used nested spaces for q and u by subdividing each triangular element into four smaller
elements, when rening the mesh. Table 3 presents results for this approach. It can be seen
that on ne discretization levels one SQP step is su-cient for fullling the stopping criterion,
which corresponds very well to the theoretical predictions made in Section 4.3. A comparison
of the results to the ones in Table 2 shows that for xed discretization level, the solution of
the identication problem on level 5 is only slightly faster than the identication of q on level
6 (with about the fourfold number of parameters) using a multi-level approach. A plot of
Level dim q dim u avg QMR it SQP it time acc. time

Table

3: CPU-time per level, accumulated time and number of inner (QMR) and outer (SQP)
iterations for exact data using a nested multi-level approach
Figure

3: Parameter distribution for exact data at level 4 using a nested multi-level approach
the parameter can be found in Figure 3. Here the approximation of the parameter in the
area where it can not identied is by far better than in the classical approach using only one
discretization level (compare Figure 2).
8 Conclusions and Outlook
We have developed a framework for Galerkin-type approximations of the LMSQP-method
for parameter identication problems in elliptic partial dierential equations and we have
discussed the implementation of the Galerkin LMSQP-method with iterative solution of the
KKT-system. The numerical results show that the resulting iteration method clearly outperforms
state-of-the-art methods for iterative regularization and provides a tool for the e-cient
solution of identication problems with ne discretizations. Moreover, we have developed a
multi-level version of the Galerkin-LMSQP method, which yields a further speed-up.
The crucial point for the possibility to obtain an e-cient implementation of the LMSQP-
method is the preconditioning of the KKT-system, which is then solved iteratively as an
indenite problem in the product space for state, parameter and Lagrangian variable. The
construction of such preconditioners is not a simple task and has not been discussed in detail in
the present paper, but will be investigated in [11], where dierent preconditioning techniques
will be compared.
Other numerical aspects to be investigated in future research are adaptive discretization
strategies and fast parallel solvers based on domain-decomposition techniques. The adaptive
discretization of optimal control problems, which is a closely related subject, has been
discussed by Becker et al. [3]; possibly the ideas of this work can be carried over to identi-
cation problems, too. The parallel solution of optimal control problems has been investigated
by Lions and Pironneau [28] in the case of quadratic problems; recently Biros and Ghattas
[4, 5] performed a numerical study of a parallel solver with an SQP-method for the outer and
preconditioned Krylov-subspace methods for the inner iteration. Many of their ideas seem
to be applicable also for parameter identication problems that are solved with the LMSQP-
method, which rises the hope that e-cient parallel versions of the LMSQP-method can be
designed also for large-scale identication problems such as impedance tomography.
Finally, we want to recall that the framework of this problem does not apply to transient
problems of parabolic or hyperbolic type. Since numerical methods for dierent types of
partial dierential equations have many type-specic features in general, it is not surprising
that also the numerical treatment of parameter identication problems should depend on
the type of the underlying state equation. However, it seems possible to construct e-cient
and convergent discretized methods at least in the case of parabolic equations, which is an
important task for future research.

Acknowledgments

The authors thank Dr. Walter Zulehner (University of Linz) and Dr. Joachim Schoberl (cur-
rently Texas A & M University) for useful and stimulating discussions on the preconditioning
of the indenite system (5.6).



--R

Estimation Techniques for Distributed Parameter Systems
Preconditioners for Karush-Kuhn-Tucker matrices arising in the optimal control of distributed systems

Parallel Lagrange-Newton-Krylov-Schur methods for PDE-constrai- ned optimization problems
Parallel Lagrange-Newton-Krylov-Schur methods for PDE-constrai- ned optimization problems
Analysis of the inexact Uzawa algorithm for saddle point problems
On the existence
Mixed and Hybrid Finite Element Methods
Iterative regularization of a parameter identi


Inverse Problems in Partial Di


Inexact and preconditioned Uzawa algorithms for saddle point problems
Inverse Problems in Di
Convergence rate results for iterative methods for solving non-linear ill-posed problems
QMR: a quasi-minimal residual method for non-Hermitian linear systems

Preconditioned all-at-once methods for large
The numerical solution of a control problem governed by a phase

On Broyden's method for the regularization of nonlinear ill-posed prob- lems
A projection-regularized Newton method for nonlinear ill-posed problems with application to parameter identi cation problems with nite element discretization



Sur le controle parallele des systemes distribues


Solution of sparse inde
Numerical approximation of partial di
Determination of a source term in the linear di
GMRES: a generalized minimal residual algorithm for solving non-symmetric linear systems
Control applications of reduced SQP methods
Partially reduced SQP methods for large-scale nonlinear optimization problems
Solving discretized optimization problems by partially reduced SQP methods
Global convergence of trust-region interior-point algorithms for in nite-dimensional nonconvex minimization subject to pointwise bounds
Analysis of iterative methods for saddle point problems: a uni
--TR
parameter identification;iterative regularization;galerkin methods;sequential quadratic programming;indefinite systems