domain: info math

text:
C Duval

P Erhard

Ch Faure

J Ch

Gilbert

Application of the automatic differentiation tool Odysse ´e to a system of thermohydraulic equations 1

The applicability of automatic differentiation on a set of partial differential equations governing thermohydraulic phenomena in heat exchangers is examined. More specifically, the challenge is to differentiate Thyc-1D, a 1-D mockup of the 3-D code Thyc implementing these equations, by means of the automatic differentiator Odysse ´e with as few manual interventions as possible. The program to differentiate contains 23 subroutines, including linear solvers and black-box functions, whose code is not available.

Introduction

Automatic differentiation is a set of techniques that are aimed at differentiating functions from a program that computes its values at arbitrary points (see for example [START_REF] Corliss | Automatic Differentiation of Algorithms: Theory, Implementation, and Application[END_REF][START_REF] Gilbert | La diffe ´rentiation automatique de fonctions repre ´sente ´es par des programmes[END_REF][START_REF] Rostaing-Schmidt | Diffe ´rentiation automatique : application a `un proble `me d'optimisation en me ´te ´orologie[END_REF]). The method differs from finite differences in that the value of the derivatives are computed exactly (up to rounding errors) in a much more efficient way. This is particularly true when gradients (set of partial derivatives) are to be computed. In this case, the reverse mode of automatic differentiation can be viewed as a method for generating adjoint codes automatically.

The aim of this work is to differentiate a one-dimensional mockup Thyc-1D of the 3-D code Thyc from EDF. This code implements a set of five partial differential equations governing thermohydraulic phenomena in heat exchangers. Thyc-1D contains 23 subroutines, including linear solvers and black-box functions, whose code is not available. This is achieved by improving and using the automatic differentiator Odysse ´e [START_REF] Rostaing | Automatic differentiation in Odysse ´e[END_REF], developed at INRIA. The challenge is to differentiate the program with as few manual interventions as possible, in such a way that the resulting code will be efficient with respect to the CPU time, sufficiently modest on its memory needs, and will provide accurate derivatives.

Description of the thermohydraulic application

More concretely, we aim at computing the sensitivity of a partial differential equation system answer to a set of input parameters. In fact, even if an industrial computer code is validated on a large experimental data basis, it could provide bad results particularly when it is used far from the limits of its applicability domain. Sensitivity evaluation provides information on influent parameters according to this calculation response. These parameters could be boundary or initial conditions, closure laws, : : : An uncertainty may be evaluated on calculation results by using system answers and their derivatives given by the derivated computer code. In our study, we aim at knowing closure laws influence in order to perform the most sensitive ones. The chosen model is a mock-up of an industrial computer code for thermohydraulics in bundles, Thyc-1D developed at EDF/Direction des Etudes et Recherches. This code is used to study a single or two-phase flow in reactor cores, steam generators and condensers. The differentiated system consists of: three conservation equations for the mixture (mass, momentum and energy), one conservation equation for the vapor mass and a last one for relative momentum between liquid and vapor phases. The system is shown below. @ @t + d i v u = 0 @ u @t + div( u2 + c(1 ; c)u 2 r )) = ;rp + g ; Mw @ h @t + d i v ( uh + Lc(1 ; c)ur) = @p @t + urp + q @ c @t + d i v ( c(u + ( 1 ; c)ur)) = @ur @t + ( urr)u + (((1 ; 2c)ur + u)r)ur

The main variables of the previous system are: enthalpy h, pressure p, mass flow rate of the mixture u, vapor quality c, and relative velocity between phases ur.

Many laws close this system: friction pressure loss tensor (Mw), volumetric power (q), interfacial drag coefficient (beta terms), interfacial mass transfer term ( ), : : : This last closure law is not well known and system answers sensitivities to this parameters are very important.

Having the derivated code of this model is interesting to study one of the many variables sensitivity ur, the relative velocity between vapor and liquid phases according to two parameters included in the interfacial drag coefficient (cd and qsi) and the following ones: relaxation time tau included in the previously mentioned interfacial mass exchange closure law, and power exchanged with the fluid q.

The studied flow is a liquid-vapor mixture. Entering in a heating tube at its lower part, the fluid gets out with a high vapor ratio at its top.

The automatic differentiation tool

Odysse ´e Odysse ´e can differentiate a program as a whole, detecting active variables, those whose value depends on the variables with respect to which the function is differentiated. Both the direct and reverse modes of differentiation are implemented. The first one is appropriate for computing directional derivatives, the second is adapted to the efficient computation of gradients. In Odysse ´e, the adjoint code generated by the reverse mode has the same structure as the original code.

Analysis of the code

The process of differentiation consists of four phases:

1. preparation of the code, 2. interprocedural analysis of the program, 3. differentiation of the algebraic expressions, 4. simplification of the code.

During the first phase, the system modifies the original units and puts them in a suitable form for differentiation. It rewrites the nondifferentiable operators abs, min into if-statements, which are tractable by Odysse ´e. The function statements are also in-lined into the code. The third transformation applied in this phase comes from theoretical reasons: if the code is formed of binary expressions, then the complexity of the calculation of the function and its derivative is less than four or five times the one of the function alone (this depends on the differentiation mode, see Section 3.2). We have implemented a process that splits expressions into equivalent binary expressions by introducing intermediate variables (see [START_REF] Faure | Splitting of Algebraic Expressions for Automatic Differentiation[END_REF]).

During the second phase, Odysse ´e makes an interprocedural analysis of the program, which results in a dependency graph between the input/output variables of each unit. These dependencies are stored in a data basis, which is filled by the differentiator for the available subroutines. For the subroutines whose code is not supplied (black-box subroutines of Thyc-1D or user-supplied derivatives), the user has to fill first the data basis, by declaring what are the arguments of the subroutine and which are the input or/and output variables. This information is enough for Odysse ´e to make a consistent analysis of the overall program. At that point, the system is able to propagate the active variables from the head-unit to all the other units. Thus the active inputs of each subroutine are known and Odysse ´e can differentiate them independently.

In the third phase, the unit by unit differentiation of the program is done, according to the following two rules. First, each subroutine of the program is differentiated with respect to its input variables, and not with respect to the input variables of the head-unit. Secondly, the differentiation is maximal in the sense that a subroutine is differentiated with respect to the maximum set of input variables appearing in every call statement. This allows Odysse ´e to generate only one subroutine for each unit of the original program. The methods used to differentiate a subroutine are described in the next section.

The fourth phase aims at simplifying the resulting code. The algebraic expressions are simplified in the same way as in Computer Algebra Systems except that the arguments of sums and products cannot be reordered modulo associativity and commutativity. The dead code is suppressed, which is very important in the reverse mode of differentiation, because the system generates to many saving instructions.

Two modes of differentiation

Automatic differentiation is based on two principles. The first one is that a program can be seen as a composition of functions, the second one is the chain rule.

In its simplest form a program is formed of a finite sequence of assignment statements. At each of these statements, one can associate a function leaving unchanged all the variables of the program except the one modified by the statement. Then, the whole program can be viewed a composition of these functions. Now, suppose that there are K statements in the program. Then the function f that it computes is the composition of K element functions: f = fK f2 f1. Denote by J k the Jacobian matrix of f k computed at (f k;1 f k;2 : : : f1)(x). From the chain rule, the Jacobian matrix of f at x is J = JK J2J1

and the transposed Jacobian matrix of f is of course

J > = J > 1 J > 2 J > K : Section Title 2 C. Duval et al.
The direct mode of differentiation consists in calculating the directional derivative J uwith the first formula, while the reverse mode of differentiation consists in computing J > v with the second formula. If f is defined from R n to R, one sees that the reverse mode can compute the n partial derivatives forming the gradient of f by a sequence of K matrix-vector products (take v = 1 above). Given a program P (an example is the subroutine sample given in Figure 1) and some input or independent variables, Odysse ´e can generate two kinds of codes: the tangent code T P and the cotangent code T P. The tangent code computes the product of the Jacobian matrix by a vector u (direct mode)

and the cotangent code computes the product of the transposed Jacobian matrix by a vector v (reverse mode). It is important to note that in the reverse mode, it is necessary either to store or to recalculate the intermediate values (f1(x), f2(f1(x)), : : :

in the example above) to be able to compute the transposed Jacobian matrix. 1), a new variable vtl (with suffix tl) is introduced for each variable v in the original code. It represents its directional derivative. The tangent code is obtained by inserting one assignment statement before each assignment statement of the original code. Its goal is to compute the directional derivative of the modified variable. The control structure of the code is preserved.

In the case of the subroutine sampletl, the directional derivatives of the input variables x and y are supposed given in xtl and ytl, and the directional derivative of the output variable z is placed in ztl. After execution, one has ztl := @z @x xtl + @z @y ytl:

SUBROUTINE sampletl (x, y, z, xtl, ytl, ztl) PARAMETER (n = 3) DIMENSION x(n), y(n), xtl(n), ytl(n) CALL norme2tl (x, xx, xtl, xxtl) CALL norme2tl (y, yy, ytl, yytl) IF (xx.LE.yy) THEN ztl = -yytl+xxtl z = -yy+xx ELSE ztl = yytl-xxtl z = yy-xx END IF RETURN END SUBROUTINE norme2tl (x, y, xtl, ytl) PARAMETER (n = 3) DIMENSION x(n), xtl(n) ytl = 0. y = 0. DO i = 1, n ytl = ytl+2*xtl(i)*x(i) y = y+x(i)**2 END DO RETURN END
Each subroutine of the cotangent code T P generated by Odysse ´e is formed of two parts (see Figure 3 for the cotangent code samplead associated with the subroutine sample of Figure 1). In the first part, the original code is executed with appropriate savings. In the second part, the adjoint variables vad (with suffix ad) associated with the original variables v are computed by transposition of tangent code, as explained above. As a result, the execution flow of the second part is reversed with respect to the one of the original code. For example, the do-loop DO i=1,n of the subroutine norme2 becomes DO i=n,1,-1 in norme2ad.

In the case of the subroutine samplead, the adjoint variables xad and yad are computed from the adjoint variables xad, yad, zad, and from the value of the original variables x, y, and z. After execution, one has xad := xad + @z @x zad yad := yad + @z @y zad zad := 0:

(n = 3) DIMENSION x(n), xad(n) DO i = n, 1, -1 xad(i) = xad(i)+yad*(2*x(i)) END DO RETURN END

Limitations and future work

The current limitations of Odysse ´e are the following:

common blocs with different partitions and equivalence statements are not accepted, -an argument of a subroutine cannot be a subroutine name, -only subroutines can be differentiated (not the functions or programs), -the gotos are not accepted in reverse mode. Odysse ´e generates the derivatives of all the variables of the code that depend on the input variables. In the present version, Odysse ´e does not take advantage of a selection of output variables to reduce the execution time.

Using Odysse ´e on Thyc-1D

As described in Section 2, we study how the variables cd, qsi, tau, and puisvol = q influence ur = ur. Therefore the program Thyc-1D must be differentiated with respect to these four input variables. More precisely, we want to know how the maximum (in time and space) of ur is influenced by these variables. Since this maximum occurs at the time step 306 and the node 23 in space, it is the program that computes urmax = ur(306,23) that is differentiated.

Treatment of black-box subroutines and linear solvers

Two difficulties have been encountered in differentiating Thyc -1D. First, the program has black-box subroutines, whose code is not available. Secondly, Thyc-1D contains linear solvers implementing iterative algorithms (Jacobi's and Gauss-Seidel's methods) with stopping tests, so that the number of iterations is not known beforehand and depends on the value given to the independent variables. The derivatives of the black-box subroutines have been hand-coded, using finite differences to approximate the derivatives.

For the linear solvers, we have chosen to supply the associated subroutines to Odysse ´e instead of letting it generate the code. This simplifies the code generation (differentiating iterative processes is a rather sophisticated operation, see [START_REF] Ch | Automatic differentiation and iterative processes[END_REF]) and makes the resulting code very efficient.

To be more specific, suppose that the original subroutine solves the linear system in y Ay = b

where A is a nonsingular matrix and b is a vector. Then in direct mode, the associated subroutine has to solve Ay = b to find y and the linear system A _ y = _ b ; _ Ay to find _ y. Here _ v denotes the directional derivatives of the variable v. In reverse mode, the adjoint variables A, b, and y are obtained by

> < > :

A > z = y (to solve in z) b := b + z A := A ; zy > y := 0: More details can be found in [START_REF] Eyssette | Applicabilite ´de la diffe ´rentiation automatique a `un syste `me d'e ´quations aux de ´rive ´es partielles re ´gissant les phe ´nome `nes thermohydrauliques dans un tube chauffant[END_REF].

Description of user-differentiated subroutines

The subroutines that must not be processed by Odysse ´e (here, the black-box subroutines and the linear solvers) have to be described in a file that is read by the system before it analyses the program. This description specifies what are Section Title 4 C. Duval et al.

the dummy arguments and whether they are input or output variables.

An appropriate language is used for this task. For example, we have described the black-box function tbalps by: setdummys tbalps (p1 p2 p3 p4 p5 tbalps) setinout tbalps (p2 p3) () (tbalps) () and the linear solver tdma by:

setdummys tdma (c a d b x imax) setinout tdma (c a d b x) () (c a d b x) ()
This information is enough for Odysse ´e to make a consistent analysis of the overall program.

Towards derivatives

As described in Section 3, Odysse ´e can derive a set of subroutines, whose call graph forms a tree, the root of which is called the head-unit. To give this structure to the set of units to derive in Thyc-1D, we had to build the head-unit, called princp sub, and had to specify the input variables of this subroutine.

Then the global analysis of the set of units starting with princp sub was performed by Odysse ´e, which was able to differentiate the overall program either in direct mode to generate the tangent code princp subtl or in reverse mode to generate the cotangent code princp subad. The handderivated subroutines (for the black-box functions and linear solvers) fit the required format, so that all the subroutines could be compiled and linked together.

Then, we had to write the main programs princptl and princpad that compute the four partial derivatives of urmax with respect to the four input variables, either in direct mode (by calling princp subtl) or in reverse mode (by calling princp subad).

In direct mode, princp subtl has to be called four times to get the four derivatives. The part of the code computing and printing the derivatives with respect to cd and qsi is: cdtl = 1. qsitl = 0. tautl = 0. puisvoltl = 0. call princp_subtl (npas) print *, "durmaxdcd = ", urmaxtl cdtl = 0. qsitl = 1. tautl = 0. puisvoltl = 0. call princp_subtl (npas) print *, "durmaxdqsi = ", urmaxtl ... In reverse mode, princp subad has to be called only once to get the four derivatives. The part of the code computing and printing the four derivatives is: cdad = 0. qsiad = 0. tauad = 0. puisvolad = 0. urmaxad = 1. call princp_subad (npas) print *, "durmaxdcd = ", cdad print *, "durmaxdqsi = ", qsiad print *, "durmaxdtau = ", tauad print *, "durmaxdpuisvol = ", puisvolad

Numerical results

We have compared the results obtained by the tangent and cotangent codes generated by Odysse ´e with the derivatives approximated by finite differences. For the latter, we have chosen the step-size giving the largest number of correct significant digits (this requires trying many step-sizes). Comparisons have been made in single and double precision.

Here, we present results obtained for the derivatives of all the 25 space components of ur, not only the 23rd giving the maximum of ur, with respect to the four inputs. This requires running princp subtl four times (in direct mode) and princp subad twenty five times (in reverse mode).

Comparison between tangent and cotangent codes

The differentiation of the black-box subroutines by finite differences can be done in different ways. In order to make the tangent and cotangent codes more similar, we have chosen to compute in both modes the Jacobian matrix of the function realized by each of these subroutines. In direct mode there is a cheaper way of doing, which consists in computing only the directional derivative of the function. Figures 4 and6 show the values of the 25 components of the derivatives of ur with respect to cd and qsi, as computed by the tangent code. We do not show those computed by the cotangent code because they are similar. Table 1 shows the values of the derivatives of urmax = ur(306,23) computed with the two modes of differentiation. We see that these derivatives have between 8 and 12 significant digits in common.

Comparisons with finite differences

The step-size for the finite differences has been chosen to make the derivatives as precise as possible. This optimal stepsize has been obtained in the following way. The approximate derivatives of urmax has been computed for several stepsizes and the one that gives the greatest number of unvariant significant digits has been chosen. For example, for the derivative with respect to cd, we have tried ten different step-sizes. The results are given in Table 2. In this case, we see that the optimal step-size is 5 1 0 ;6 and that the correct derivative is likely to be ;2:6877707 : : : . In Table 3, we show the approximate derivatives obtained by finite differences and those computed by Odysse ´e. In Section Title 5 C. Duval et al.

- the latter case, we only display the most significant digits in common to the tangent and cotangent derivatives. The last column shows the maximum number of correct significant digits that Odysse ´e can yield with respect to "optimal" finite differences.

In Figures 5 and7, we have compared the values of the derivatives computed by finite differences with optimal stepsize, with those computed by Odysse ´e. The curves show the relative differences between two methods (the vertical axis has a logarithmic scale): one compares the tangent code to finite differences (dashed line), another compares the cotangent code to finite differences (dotted line) and the last one Step-size Approximate derivative compares the tangent and cotangent codes (plain line). The lowest curve indicates the couple of methods that most agree.

In each case, it is formed of the tangent and cotangent codes of Odysse ´e. The same tests have been performed in single precision. The derivatives of urmax obtained by Odysse ´e are given in Table 4.

As in double precision, we have chosen the "optimal" stepsize for computing approximate derivatives by finite differences and have compared them with the one obtained by Odysse ´e. The results are shown in Table 5 Table 3. If we trust the first significant digits obtained by finite differences in double precision, we see that with the optimal step-size, only the last digit may be erroneous in the second column of Table 5. We also see that the digits in common to the direct and reverse mode of Odysse ´e are not necessary correct: up to 2 digits can be erroneous. This implies that the number in the column 'Max gain' in Table 3 are not the actual number of digits that Odysse ´e can yield with respect to finite differences. It is an upper bound for this number. On the other hand, always based on the digits obtained in double precision, we can give the number of digits that Odysse ´e yields in single precision, with respect to finite differences. They are given in the last column of Table 5.

Finite differences

Odysse

Efficiency in memory-space and CPU time

In order to compare the efficiency of the codes generated by Odysse ´e, we have generated the derivatives of Thyc-1D with respect to 1 input variable (cd), 2 input variables (cd, qsi), 3 input variables (cd, qsi, tau), and 4 input variables (cd, qsi, tau, puisvol). Let T (P ) be the execution time of the original code, T (T P ) be the execution time of the tangent code, T (T P) be the execution time of the cotangent code, and n be the number of input variables. Extrapolating the results in Table 6 leads to the following ratios: T (T P ) T (P ) 2:3 n and T (T P) T (P ) 4:6:

These numbers are in accordance with the theoretical bounds, which are 4n for the first ratio and 5 for the second (see for example [START_REF] Morgenstern | How to compute fast a function and all its derivatives, a variation on the theorem of Baur-Strassen[END_REF][START_REF] Griewank | On automatic differentiation[END_REF][START_REF] Gilbert | La diffe ´rentiation automatique de fonctions repre ´sente ´es par des programmes[END_REF]).

In terms of the size of the code, the tangent code takes 35% more memory space whereas the cotangent code takes about twice that of the original code. As for the size of the data space, it is approximately doubled in the tangent code and multiplied by 25 in the cotangent code. One must notice that the only memory management strategy implemented in the current version of Odysse ´e saves all the modified variables. Other strategies are being studied for the next version of Odysse ´e.

Conclusion

These tests show that, in double precision, 7 or 8 correct significant digits can be obtained by "optimal" finite differences. With respect to this, the derivatives obtained by Odysse ´e are likely to have between 1 and 3 more correct significant digits. We believe that the accuracy of the derivatives obtained by Odysse ´e have been limited by the need to differentiate by finite differences frequently used black-box functions.

As for the CPU time efficiency, it fits the theory: the evaluation of the tangent code takes about 2:3 times the time to run the original code, and the evaluation of the cotangent code takes around 4:6 times the time to evaluate the original function. The size of the memory used by the tangent code is twice as large as the memory used by the original code, while the cotangent code used about 25 times more memory than the original code. It is likely that this last number could be reduced by using more sophisticated strategies for determining the information to save during the computation of the original code. This optimization will become necessary when automatic differentiation of the industrial code Thyc Section Title 7 C. Duval et al.

will be undertaken, which is next stage of the EDF-INRIA collaboration.

Figure 1 .

 1 Figure 1. Sample code

Figure 2 .

 2 Figure 2. Tangent code of sample generated by Odysse ´e, for the input variables x and y

Figure 3 .

 3 Figure 3. Cotangent (or adjoint) code of sample generated by Odysse ´e, for the input variables x and y

Figure 4 .Figure 5 .

 45 Figure 4. Derivatives in direct mode with respect to cd

Figure 6 .Figure 7 .

 67 Figure 6. Derivatives in direct mode with respect to qsi

Table 1 .

 1 Derivatives of ur(306,23) by Odysse ´e in double precision

Table 2 .

 2 Trial step-sizes for the derivative of urmax with respect to cd.

	5:E ; 11 ;2:687743361207140E + 0 0 5:E ; 10 ;2:687798428269161E + 0 0 5:E ; 09 ;2:687771161191677E + 0 0 5:E ; 08 ;2:687771023524022E + 0 0 5:E ; 07 ;2:687770758846852E + 0 0 5:E ; 06 ;2:687770744280727E + 0 0 5:E ; 05 ;2:687770756817365E + 0 0 5:E ; 04 ;2:687771680981221E + 0 0 5:E ; 03 ;2:687864088674541E + 0 0 5:E ; 02 ;2:697140333223182E + 0 0

  . To give an interpretation to these results, let us compare them to those in

	cd qsi tau pui.	Finite differences ;2:6877707D + 0 0 ;2:68777074766D + 0 0 +4 Odysse ´e Max gain ;1:1283833D ; 02 ;1:1283833581 D ; 02 +3 2:395347 D ; 01 2:3953473 D ; 01 +1 8:4189731D ; 09 8:4189731 D ; 09 +0
	Section Title	6	C. Duval et al.

Table 3 .

 3 Comparison of derivatives in double precision

	cd qsi tau puisvol	Tangent code ;2:68763E + 0 0 ;2:68763E + 0 0 Cotangent code ;1:12781E ; 02 ;1:12782E ; 02 2:39515E ; 01 2:39520E ; 01 8:41967E ; 09 8:41967E ; 09

Table 4 .

 4 Derivatives of ur(306,23) by Odysse ´e in single precision

Table 5 .

 5 Comparisons of derivatives in single precision

Table 6 .

 6 Table 6 shows the results: the CPU column gives the CPU time in seconds, the text column gives the text length in bytes, and the bss column gives the size of the initialized data in bytes. Execution time and memory space

			Tangent code			Cotangent code
		CPU	text	bss	CPU	text	bss
	0	5:72 270336 105608	5:93 270336	105608
	1 14:08 417792 266680 27:61 573440 2540712
	2 26:91 417792 266712 27:30 573440 2540912
	3 38:25 417792 267592 27:20 581632 2540544
	4 55:21 417792 267616 27:31 581632 2540608

This work was initiated and financed in part by EDF-DER.

EDF-DER, De ´partement Physique des Re ´acteurs, 1, avenue du Ge ´ne ´ral de Gaulle, F-92141 Clamart Cedex, France.

INRIA Sophia Antipolis, 2004 route des Lucioles, BP 93, F-06902 Sophia-Antipolis Cedex, France.

INRIA Rocquencourt, BP 105, F-78153 Le Chesnay Cedex, France.


Jean Charles Gilbert

Piecewise line-search techniques for constrained minimization by quasi-Newton algorithms

Keywords: constrained optimization, nonlinear programming, piecewise line-search, quasi-Newton BFGS updates, reduced curvature condition, reduced Hessian method, SQP algorithm, Wolfe's line-search

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

The paper is organized as follows. In Section 1, we give a succinct presentation of quasi-Newton algorithms for solving unconstrained and equality constrained minimization problems. We recall what is the process to maintain the positive definiteness of the updated matrices in unconstrained optimization and the difficulties to extend this approach to constrained problems.

In section 2, we discuss the different choices of pairs (γ k , δ k ) that have been proposed to update the matrices in the quasi-Newton versions of the reduced Hessian and SQP algorithms. We also sketch the principles of the piecewise line-search (PLS), further developed in Sections 4 and 5.

In Section 3, we present a new argument, based on geometrical considerations, that can be helpful in determining the form of the pairs (γ k , δ k ) to use for updating the matrices in reduced quasi-Newton algorithms. It is argued that, when γ k is the change in the reduced gradient, the bases tangent to the constraint manifolds should satisfy a so-called "zero Lie bracket bases" property, which is equivalent to the existence of local parametrizations of these manifolds, whose derivatives give the bases. If this property holds for the tangent bases computed by partitioning the Jacobian matrix of the constraints, it is usually not satisfied when the tangent bases are orthonormal.

Section 4 deals with the PLS technique for minimizing functions subject to equality constraints by reduced quasi-Newton algorithms. The technique is viewed as an explicit Euler discretization of a smooth curve, along which a so called "reduced curvature condition" can be satisfied, while decreasing a merit function. A new update criterion is introduced, with which a strong q-superlinear convergence result can be proved.

In Section 5, it is shown how the PLS technique can be used for the quasi-Newton version of the SQP algorithm. A strategy to maintain positive definite the matrices M k approximating of the Hessian of the augmented Lagrangian is presented. It is based on two observations. The first one is that the penalty parameter can be set to ensure d ⊤ M k d > 0 when d is not in the tangent space to the constraint manifold. The second one is that the PLS can ensure the positivity of d ⊤ M k d when d is in this tangent space. The update scheme combines these two ideas without having to switch between them.

We end the paper with a conclusion section listing a few open problems.

1 Quasi-Newton algorithms

In this section, we present quasi-Newton algorithms for solving unconstrained (Section 1.1) and equality constrained (Section 1.2) minimization problems. We restrict our presentation to line-search methods and to algorithms having the potential to converge one step q-superlinearly (or one step q-quadratically, when second derivatives are computed). We emphasize the similarities between the approaches.

The Wolfe line-search in unconstrained optimization

Consider the problem of minimizing a smooth function f : R n → R, which we write: min f (x).

(1.1)

Descent direction algorithms solve problem (1.1) by generating a sequence of approximations x k ∈ R n of a solution x * of (1.1). These iterates are obtained by the recurrence

x k+1 = x k + α k d k ,
where d k is a descent direction of f at x k , which means that f ′ (x k ) • d k < 0, and α k > 0 is a step-size found by a technique called line-search. The aim of this technique is, in particular, to force the decrease of f at each iteration. Frequently, this is controlled by requiring α k to satisfy the inequality

f (x k + α k d k ) ≤ f (x k ) + ω 1 α k ∇f (x k ) ⊤ d k , (1.2)
where ω 1 ∈ (0, 1) is a constant and ∇f (x k ) denotes the gradient of f at x k associated with the Euclidean (or ℓ 2 ) inner product. Since the last term in (1.2) is negative, f decreases at each iteration. The type of direction d k characterizes the minimization algorithm. In quasi-Newton methods, this direction takes the form

d k = -M -1 k ∇f (x k ),
where M k is a nonsingular symmetric matrix approximating the Hessian of f at x k .

In practice, it is suitable to maintain M k positive definite, for two reasons that we will meet again when we will consider algorithms for constrained problems. First, this is natural since the Hessian of f is positive definite at a strong minimum, i.e., a solution of (1.1) satisfying the second order sufficient conditions of optimality. Secondly, this makes d k a descent direction of f at x k : f ′ (x k ) • d k = -∇f (x k ) ⊤ M -1 k ∇f (x k ) < 0, when x k is not stationary.

The way of maintaining the positive definiteness of M k is the central theme of this paper. In unconstrained optimization, it is obtained by a nice combination of ideas. First M k is updated by a formula that allows this property to be sustained from one iteration to the next one. An example, which is very efficient in practice, is the BFGS formula:

M k+1 = M k - M k δ k δ ⊤ k M k δ ⊤ k M k δ k + γ k γ ⊤ k γ ⊤ k δ k . (1.3)
In (1.3), γ k and δ k are two vectors of R n ; γ k gives the change in the gradient of f and δ k is the step between two successive iterates:

γ k = ∇f (x k+1 ) -∇f (x k ) and δ k = x k+1 -x k = α k d k .
Since M k+1 defined by (1.3) satisfies the so-called quasi-Newton equation

γ k = M k+1 δ k , (1.4)
it incorporates information from the Hessian of f . On the other hand, taking the inner product of the two sides of (1.4) with

δ k gives γ ⊤ k δ k = δ ⊤ k M k+1 δ k . Hence, if M k+1 is positive definite, one has γ ⊤ k δ k > 0. (1.5)
An interesting property of formula (1.3), which will play a key role in this paper, is that the converse is true: when M k is positive definite and (1.5) holds, then M k+1 is also positive definite. This implies that satisfying (1.5) is crucial for the success of the algorithm.

In unconstrained optimization, (1.5) is realized by the line-search. The so-called Wolfe line-search [START_REF] Wolfe | Convergence conditions for ascent methods[END_REF][START_REF] Wolfe | Convergence conditions for ascent methods II: some corrections[END_REF] is the appropriate concept. It finds a step-size α k > 0 satisfying (1.2) and

∇f (x k + α k d k ) ⊤ d k ≥ ω 2 ∇f (x k ) ⊤ d k , (1.6)
where ω 1 < ω 2 < 1. These two conditions (1.2) and (1.6) can be satisfied simultaneously under mild assumptions and there are specific algorithms that can realize this; see [START_REF] Fletcher | Practical Methods of Optimization[END_REF][START_REF]A view of line-searches[END_REF][START_REF] Moré | Line search algorithms with guaranteed sufficient decrease[END_REF][START_REF] Armand | A line-search technique with sufficient decrease and curvature conditions[END_REF]. By subtracting ∇f (x k ) ⊤ d k from both sides of (1.6), we see that (1.6) implies (1.5), so that we finally have the series of implications:

(1.6) =⇒ γ ⊤ k δ k > 0 =⇒ M k+1 is positive definite =⇒ d k+1 is a descent direction.
All of this is described in detail in textbooks or review articles such as [START_REF] Gill | Practical Optimization[END_REF][START_REF] Dennis | Numerical Methods for Unconstrained Optimization and Nonlinear Equations[END_REF][START_REF] Fletcher | Practical Methods of Optimization[END_REF][START_REF] Powell | A review of algorithms for nonlinear equations and unconstrained optimization[END_REF][START_REF] Dennis | A view of unconstrained optimization[END_REF][START_REF]Theory of algorithms for unconstrained optimization[END_REF][START_REF] Fletcher | An overview of unconstrained optimization[END_REF].

If satisfying inequality (1.5) is mandatory to maintain the positive definiteness of the updated matrices, the way of satisfying it in constrained optimization is still a subject of debate and research. In this paper, we describe a step-size determination process that implies (1.5), so generalizing the Wolfe line-search. Note that, for constrained problems, γ k is no longer the change in the gradient of f . Also, the meaning of γ k and its size change according to the algorithm considered.

Two classes of quasi-Newton algorithms for constrained problems

Consider now the problem of minimizing a smooth function f : R n → R subject to equality constraints given by a smooth function c : R n → R m (m < n):

min f (x) c(x) = 0. (1.7)
We consider problems with equality constraints only. Algorithms dealing with inequalities are often of two types: active set methods (see [START_REF] Gill | Practical Optimization[END_REF][START_REF] Fletcher | Practical Methods of Optimization[END_REF]) or interior point methods (see [START_REF] Byrd | A trust region method based on interior point techniques for nonlinear programming[END_REF][START_REF] Coleman | An interior trust region approach for nonlinear minimization subject to bounds[END_REF][START_REF] Dennis | Trust-region interior-point SQP algorithms for a class of nonlinear programming problems[END_REF]). These algorithms may have to solve subproblems like (1.7), so that algorithms for solving (1.7) are also interesting in that case.

The optimality conditions at a solution x * of (1.7) can be written: there exists

λ * ∈ R m such that ∇ℓ(x * , λ * ) = 0 c(x * ) = 0, (1.8)
where ℓ(x, λ) = f (x) + λ ⊤ c(x) is the Lagrangian of the problem and ∇ denotes the gradient with respect to x. The existence of λ * is ensured when the m × n Jacobian matrix of the constraint at the solution

A(x * ) = c ′ (x * ) is surjective. We assume that A(x) = c ′ (x) is surjective for all x ∈ R n .
Remark 1.1. Below, more hypotheses will be made on the problem data. Some of them are very restrictive when they are assumed to hold for all x ∈ R n (for example, the existence of a continuous field of bases tangent to the constraint manifold). It would be cleaner to specify the open set Ω containing x * on which these hypotheses hold. However, the case when Ω = R n is more complicated to describe, since the line-search has to control the step-size such that the next iterate is also in Ω. To avoid unnecessary complications, which would hide the essence of the algorithms, we will not do so, knowing that simple modifications can take into account the case when Ω = R n .

The SQP algorithm

The Sequential Quadratic Programming (SQP) method determines the search direction d k at an iterate (x k , λ k ) by solving the quadratic program in the variable

d ∈ R n : min ∇f (x k ) ⊤ d + 1 2 d ⊤ M k d c(x k ) + A(x k )d = 0, (1.9)
where M k is the Hessian of the Lagrangian evaluated at (x k , λ k ) or an approximation to it (see [START_REF] Fletcher | Practical Methods of Optimization[END_REF][START_REF] Gill | Constrained nonlinear programming[END_REF][START_REF] Boggs | Sequential quadratic programming[END_REF]). It also makes sense to force M k to be an approximation of the Hessian of the augmented Lagrangian

ℓ τ (x, λ) = ℓ(x, λ) + τ 2 c(x) 2 2 ,
where • 2 denotes the ℓ 2 -norm (see [START_REF] Han | Superlinearly convergent variable metric algorithms for general nonlinear programming problems[END_REF][START_REF] Tapia | Diagonalized multiplier methods and quasi-Newton methods for constrained optimization[END_REF][START_REF] Tapia | On secant updates for use in general constrained optimization[END_REF]). This is the point of view that we take in this paper since, for τ sufficiently large, ∇ 2 ℓ τ (x * , λ * ) is positive definite at a strong solution of (1.7) (a solution satisfying the strong sufficient conditions of optimality), so that it is reasonable to update positive definite approximations M k .

As in unconstrained problems, another reason to maintain M k positive definite is that the solution d k of (1.9) can be a descent direction of the merit function

Θ σ (x) = f (x) + σ c(x) P , (1.10)
where σ > 0 is a penalty parameter and

• P denotes an arbitrary norm on R m . Indeed Θ ′ σ (x k ; d k ) = -d ⊤ k M k d k + (λ QP k ) ⊤ c(x k ) -σ c(x k ) P , where λ QP
k is the multiplier associated with the constraints in (1.9). Hence Θ

′ σ (x k ; d k ) is negative if M k is positive definite, σ > λ QP k D , (1.11)
and x k is not stationary for (1.7). We have used the dual norm

u D = sup v P =1 u ⊤ v.
For the globalization of the SQP algorithm, i.e., for forcing its convergence from a remote starting point x 1 , the following property is important. When σ > λ * D , a strong solution of (1.7) is a strict local minimum of Θ σ : Θ σ is said to be an exact penalty function. In that case, Θ σ can be used to measure the progress to optimality and to determine the step-size α k along the descent direction d k . The natural extension of inequality (1.2) is then

Θ σ k (x k + α k d k ) ≤ Θ σ k (x k ) + ω 1 α k Θ ′ σ k (x k ; d k ).
(1.12)

In this inequality we have allowed the penalty parameter σ k to change at each iteration to satisfy (1.11). To ensure convergence, the step-size cannot be too small and, for example, can be determined by backtracking [START_REF] Armijo | Minimization of functions having Lipschitz continuous first partial derivatives[END_REF]. Then, the next iterate is set to

x k+1 = x k + α k d k .

The reduced Hessian algorithm

The reduced Hessian algorithm is particularly useful when n-m ≪ n, because, as we will see, its quasi-Newton version, also called reduced quasi-Newton algorithm, needs to update matrices of order nm only, while the SQP method requires updating matrices of order n.

An iteration of the reduced Hessian algorithm consists of two stages; see [START_REF] Coleman | Nonlinear programming via an exact penalty function: asymptotic analysis[END_REF][START_REF] Hoyer | Variants of the reduced Newton method for nonlinear equality constrained optimization problems[END_REF][START_REF] Ch | On the local and global convergence of a reduced quasi-Newton method[END_REF][START_REF] Byrd | On the convergence of constrained optimization methods with accurate Hessian information on a subspace[END_REF]. Starting an iteration at x k ∈ R n , a tangent step t k is first computed by solving the quadratic program in the variable t ∈ R n :

min ∇f (x k ) ⊤ t + 1 2 t ⊤ M k t A(x k )t = 0, (1.13)
where M k is again the Hessian of the Lagrangian or an approximation to it. This problem is similar to (1.9), except that t k is now in the null space of A(x k ) which is also the tangent space at x k to the manifold

M x k = {x ∈ R n : c(x) = c(x k )}.
Therefore, only the part of M k acting on this tangent space is visible by the algorithm (M k is not used in the second stage of the algorithm). It is that part that will have to be updated in the quasi-Newton version of the method.

Next, a Newton-like direction r k for solving the second equation in (1.8) is determined. For this, one uses a right inverse

A -(x k ) of the Jacobian A(x k ) (hence A(x k )A -(x k ) = I m ∈ R m×m ) and one computes r k = -A -(x k )c(x k + t k ). (1.14)
Finally, in the local algorithm (without line-search), x k+1 is set to

x k + t k + r k . Now, let Z -(x k
) be a basis of the null space N (A(x k )) of A(x k ), i.e., an n×(n-m) injective matrix such that A(x k )Z -(x k ) = 0. Then, the solution of (1.13) can be written

t k = -Z -(x k )G -1 k g(x k ), (1.15)
where

G k = Z -(x k ) ⊤ M k Z -(x k ) is or approximates the reduced Hessian of the La- grangian G(x, λ) = Z -(x) ⊤ ∇ 2 ℓ(x, λ)Z -(x) and g(x) = Z -(x) ⊤ ∇f (x)
is the reduced gradient of f . Note that for any choice of right inverse A -(x) of A(x) and basis Z -(x) of N (A(x)), there exists an (nm) × n matrix Z(x), uniquely determined by

A -(x)A(x) + Z -(x)Z(x) = I n .
(1.16)

We deduce from this that

Z(x)Z -(x) = I n-m and Z(x)A -(x) = O (n-m)×m , so that Z -(x) is a right inverse of Z(x).
When G k is a quasi-Newton approximation of the reduced Hessian of the Lagrangian, it is interesting to maintain this matrix positive definite for two reasons similar to those met in Section 1.1 for unconstrained problems and in this section for the SQP algorithm. First, the reduced Hessian of the Lagrangian is positive definite at a strong solution of (1.7). Secondly, when G k is positive definite and g(x k ) = 0, the directional derivative

Θ ′ σ k (x k ; t k ) = -g(x k ) ⊤ G -1 k g(x k ), is negative, making t k a descent direction of Θ σ k at x k .
Then, the globalization of the reduced quasi-Newton algorithm can be done along the following lines [START_REF] Ch | On the local and global convergence of a reduced quasi-Newton method[END_REF][START_REF] Ch | Superlinear convergence of a reduced BFGS method with piecewise line-search and update criterion[END_REF], similar to those followed above. Since the direction r k is computed by evaluating functions at two different points, x k and x k + t k , there is no guarantee that t k + r k will be a descent direction of Θ σ k at x k . On the other hand, introducing

r 0 k = -A -(x k )c(x k ), (1.17)
we observe that

Θ ′ σ k (x k ; t k + r 0 k ) = -g(x k ) ⊤ G -1 k g(x k ) + λ(x k ) ⊤ c(x k ) -σ k c(x k ) P ,
where λ(x) is the multiplier estimate

λ(x) = -A -(x) ⊤ ∇f (x). (1.18)
We deduce from this that

t k + r 0 k is a descent direction of Θ σ k provided σ k > λ(x k ) D
and x k is not stationary for (1.7). Therefore, if σ k is sufficiently large, one can force the decrease of Θ σ k by determining a step-size α k > 0 along the quadratic curve joining x k and x k + t k + r k , and tangent to

t k + r 0 k at x k , α → p 0 k (α) = x k + α(t k + r 0 k ) + α 2 (r k -r 0 k ).
The step-size will be determined, for example by backtracking, so that

Θ σ k (p 0 k (α k )) ≤ Θ σ k (x k ) + ω 1 α k Θ ′ σ k (x k ; t k + r 0 k ). (1.19)
The next iterate is then set to

x k+1 = p 0 k (α k ). 1.3 Realizing γ ⊤ k δ k > 0
If we compare the two algorithms of the previous section, with the unconstrained algorithm of Section 1.1, we see that what is missing in the algorithms for constrained problems is the counterpart of condition (1.6), which ensures the inequality γ ⊤ k δ k > 0 and by that the positive definiteness of the updated matrices.

In constrained optimization, realizing γ ⊤ k δ k > 0 is a much harder task, at least if one tries to mimic the Wolfe line-search used in unconstrained optimization. Remark, however, that once this condition is realized and is compatible with the inequalities (1.12) or (1.19), the algorithms are quite simple to implement, since they amount to do what is sketched in Section 1.2. In Sections 4 and 5, we dissect the content of the black boxes, introduced in [START_REF] Ch | On the realization of the Wolfe conditions in reduced quasi-Newton methods for equality constrained optimization[END_REF][START_REF] Armand | A piecewise line-search technique for maintaining the positive definiteness of the updated matrices in the SQP method[END_REF], that realize the condition γ ⊤ k δ k > 0, by means of a step-size determination process.

Regardless of the meaning of γ k and δ k , there is always the possibility to get γ ⊤ k δ k > 0 by using the modification of the BFGS formula suggested by Powell [START_REF] Powell | A fast algorithm for nonlinearly constrained optimization calculations[END_REF]. In that approach, the vector γ k is modified to γk = θγ k + (1θ)M k δ k , where M k stands for the matrix to update and θ is the number in (0, 1], the closest to 1, such that the inequality γ⊤

k δ k ≥ η δ ⊤ k M k δ k is satisfied (the constant η ∈ (0, 1
) is usually set to 0.1 or 0.2 [START_REF] Powell | Algorithms for nonlinear constraints that use Lagrangian func-Mathematical Programming[END_REF][START_REF] Powell | The performance of two subroutines for constrained optimization on some difficult test problems[END_REF]). Then γk is used instead of γ k in the BFGS formula (1.3). Powell's correction of γ k is certainly the most widely used technique in practice. Its success is due to its appealing simplicity and its usually good numerical efficiency. We believe, however, that this technique is not satisfactory for at least three reasons: (i) it is hard to figure out what part of the problem data is taken into account by this modification, (ii) the asymptotic r-superlinear rate of convergence that can be proved in theory [START_REF] Powell | The convergence of variable metric methods for nonlinearly constrained optimization calculations[END_REF] is not as strong as one could expect; the stronger q-superlinear convergence that is obtained by the BFGS algorithm for unconstrained problems would be more satisfactory, and (iii) the technique can yield artificial illconditioning of the updated matrices, deteriorating the numerical efficiency of the algorithms (see [46, p. 125] and [START_REF] Powell | The performance of two subroutines for constrained optimization on some difficult test problems[END_REF]). These facts have motivated further studies (see also [START_REF] Coleman | Partitioned quasi-Newton methods for nonlinear constrained optimization[END_REF][START_REF] Gurwitz | Local convergence of a two-piece update of a projected Hessian matrix[END_REF][START_REF] Biegler | A reduced Hessian method for large-scale constrained optimization[END_REF][START_REF] Coleman | A quasi-Newton L 2 -penalty method for minimization subject to nonlinear equality constraints[END_REF]).

Choices of γ k and δ k for constrained problems

Let us now consider how to choose the vectors γ k and δ k in the algorithms described in Section 1.2. If we denote by M * the matrix at the optimal point that is approximated by M k , what must guide us is the necessity to have

γ k -M * δ k δ k → 0, when k → ∞. (2.1)
This estimate is required by the asymptotic analysis. Observe that this rule is satisfied in unconstrained optimization, since M k approximates ∇ 2 f (x k ) and, when the Hessian of f is Lipschitz continuous and

x k → x * , one has γ k = ∇ 2 f (x k )δ k + O( δ k 2 ) = ∇ 2 f (x * )δ k + o( δ k ), when k → ∞.
In this section, we limit the discussion on the choices of γ k and δ k for algorithms without line-search, since the line-search techniques described in Section 1.2 will be modified anyway. Our aim is to give insight into what γ k and δ k should be near the solution and to make explicit the difficulties that are encountered.

Vectors γ k and δ k in reduced quasi-Newton algorithms

Let us consider first the choice of γ k and δ k in reduced quasi-Newton algorithms. We observe that, since the reduced gradient is given by g

(x) = Z -(x) ⊤ ∇f (x) = Z -(x) ⊤ ∇ℓ(x, λ * ), its derivative at a stationary point x * is g ′ (x * ) = Z -(x * ) ⊤ L * , (2.2)
where L * = ∇ 2 ℓ(x * , λ * ) is the Hessian of the Lagrangian with respect to x at (x * , λ * ).

Since the updated matrices G k must approach the reduced Hessian of the Lagrangian at the optimal point,

G * = Z -(x * ) ⊤ L * Z -(x * ) = g ′ (x * )Z -(x *
), the rule (2.1) and formula (2.2) suggest to take

γ k = g(x k + t k ) -g(x k ) δ k = Z(x k )t k .
(2.3) Indeed, with this choice

t k = Z -(x k )δ k and, when g ′ is Lipschitz continuous, we have γ k = Z -(x * ) ⊤ L * t k + o( t k ) = G * δ k + o( δ k ), when k → ∞. (2.4)
Another possibility is to take

γ k = Z -(x k ) ⊤ ∇ℓ(x k + t k , λ k ) -∇ℓ(x k , λ k ) δ k = Z(x k )t k , (2.5)
where λ k is λ(x k ) or the multiplier associated with the constraints in (1.13). If f ′′ and c ′′ are Lipschitz continuous, we obtain the same estimate (2.4) when (x k , λ k ) converges to (x * , λ * ).

Other formulas are proposed in [START_REF] Coleman | On the local convergence of a quasi-Newton method for the nonlinear programming problem[END_REF][START_REF] Nocedal | Projected Hessian updating algorithms for nonlinearly constrained optimization[END_REF][START_REF] Byrd | An analysis of reduced Hessian methods for constrained optimization[END_REF]. Close to a solution, they are all equivalent, so that an asymptotic analysis cannot distinguish between them. Now, their global efficiency can be very different. In Section 3 below, we bring a new argument that can help in selecting one of these formulas according to the choice of tangent bases Z -.

The choices of γ k and δ k mentioned above are safe, since they yield (2.1), but they can be very expensive in computing time. They require indeed linearizing functions (i.e., computing ∇f , A and sometimes Z -) at the intermediate point x k + t k , while this linearization is not necessary in the local method of Section 1.2. Therefore, researchers have suggested to avoid this linearization by taking [START_REF] Nocedal | Projected Hessian updating algorithms for nonlinearly constrained optimization[END_REF], even in a global framework [START_REF] Ch | Mise à jour de la métrique dans les méthodes de quasi-Newton réduites en optimisation avec contraintes d'égalité[END_REF][START_REF] Ch | Superlinear convergence of a reduced BFGS method with piecewise line-search and update criterion[END_REF], but in practice the numerical results are sometimes disappointing. Now, the criterion (2.7) is rather crude because there is no need to have a small transversal displacement r k to update the matrix. This is particularly clear when formula (2.3) is used for γ k . We see that what is necessary for a safe update is to have g(x k+1 ) close to g(x k + t k ), which can occur even if x k+1 and x k + t k are far from each other, provided these two points are close the same reduced gradient manifold. In other words, and this formulation of the remark is also valid when the step-size differs from one, the useful information to update G k is the tangent first order approximation of the vector joining x k and the intersection of the constraint manifold and the reduced gradient manifold. To our knowledge, there is no update criterion based on this observation. We believe that this topic deserves more attention. Now that the asymptotic criterion (2.1) has led us to appropriate formulas for γ k and δ k , one can ask whether a line-search can help in getting γ ⊤ k δ k > 0. The answer is negative in general, even if a search is made along the tangent direction t k . An example with (2.3) is given in [START_REF] Ch | Maintaining the positive definiteness of the matrices in reduced secant methods for equality constrained optimization[END_REF]. Therefore, more sophisticated step-size determination techniques have to be introduced. This is the matter of Section 4.

γ k = g(x k+1 ) -g(x k ) or γ k = Z -(x k ) ⊤ (∇ℓ(x k+1 , λ k ) -∇ℓ(x k , λ k )) (2.

Vectors γ k and δ k in the SQP algorithm

Let us now consider the case of the SQP algorithm and let us use the same criterion (2.1) for determining γ k and δ k .

Since M k may approximate the Hessian of the augmented Lagrangian, a first possibility is to use

γ k = ∇ℓ τ (x k+1 , λ QP k ) -∇ℓ τ (x k , λ QP k ) δ k = x k+1 -x k .
where λ QP k is the dual solution of (1.9) (see [START_REF] Han | Superlinearly convergent variable metric algorithms for general nonlinear programming problems[END_REF][START_REF] Tapia | Diagonalized multiplier methods and quasi-Newton methods for constrained optimization[END_REF][START_REF] Glad | Properties of updating methods for the multipliers in augmented Lagrangians[END_REF]). This approach has however serious practical difficulties: (i) a priori knowledge of the threshold value of τ making positive definite the Hessian of ℓ τ at the solution is generally unavailable, (ii) large values of τ present severe numerical problems (examples are given in [START_REF] Tapia | Diagonalized multiplier methods and quasi-Newton methods for constrained optimization[END_REF][START_REF] Nocedal | Projected Hessian updating algorithms for nonlinearly constrained optimization[END_REF]), and (iii) far from the solution, there may be no value of τ and no step-size along d k that make the inner product γ ⊤ k δ k positive. Some of the inconvenients of the preceding formula can be overcome by taking advantage of the structure of the Hessian of the augmented Lagrangian at the solution:

∇ 2 ℓ τ (x * , λ * ) = L * + τ A(x * ) ⊤ A(x * ). (2.8)
This suggests to take [START_REF] Tapia | On secant updates for use in general constrained optimization[END_REF]

γ k = ∇ℓ(x k+1 , λ QP k ) -∇ℓ(x k , λ QP k ) + τ A(x k ) ⊤ A(x k )δ k δ k = x k+1 -x k .
(2.9)

The inner product of γ k and δ k is then

γ ⊤ k δ k = ∇ℓ(x k+1 , λ QP k ) -∇ℓ(x k , λ QP k ) ⊤ δ k + τ A(x k )δ k 2 ,
so that, one can get γ ⊤ k δ k > 0 by taking τ sufficiently large, as long as A(x k )δ k = 0. It is clear that this strategy does not work when A(x k )δ k is numerically zero and γ ⊤ k δ k is negative. In [START_REF] Byrd | An SQP augmented Lagrangian BFGS algorithm for constrained optimization[END_REF], Byrd, Tapia, and Zhang propose the back-up strategy that consists in replacing A(x k ) ⊤ A(x k )δ k by δ k in formula (2.9) when A(x k )δ k is small and γ ⊤ k δ k is not sufficiently positive. Then the positivity of γ ⊤ k δ k can be obtained as before by taking τ sufficient large. Numerical experiment in [START_REF] Byrd | An SQP augmented Lagrangian BFGS algorithm for constrained optimization[END_REF] has shown that this approach is numerically competitive with Powell's correction of the BFGS formula. They also proved that the convergence of the sequence {x k } implies its r-superlinear convergence and even its q-superlinear convergence if the penalty parameter τ is eventually maintained fixed and sufficiently large. This nice result is not completely satisfactory because the threshold value for τ giving the q-superlinear result is usually unknown.

An interesting aspect of the approach in [START_REF] Byrd | An SQP augmented Lagrangian BFGS algorithm for constrained optimization[END_REF] is to give an update rule for the penalty parameter τ . Clearly, the technique appropriately deals with the transversal component of the matrix M k (its action in the range space of A ⊤ ) by setting the parameter τ , but it needs a backup strategy for its longitudinal component (its action in the null space of A). In particular, when the constraints are linear and the iterates are feasible and away from the solution, the back-up strategy may often be active, which is not very desirable. It is clear that in this particular case Wolfe's line-search would overcome the difficulty. In Section 5, we will present a technique combining the use of τ for dealing with the transversal part of M k and a piecewise line-search that takes care of the longitudinal part of M k .

Principle of the new approach

The piecewise line-search (PLS) techniques presented in Sections 4 and 5 are based on the following principle. First, it is observed that a so-called reduced curvature condition can be satisfied along a curve that is defined as the solution of an ordinary differential equation. This condition implies the positivity of γ ⊤ k δ k , in the same way as (1.6) implies (1.5). Computing a step-size by moving along this curve would be computationally too expensive in general. Also, a piecewise linear approximation of the curve is introduced, using an explicit Euler discretization. At each point of discretization, the reduced curvature condition is tested. If it holds the PLS is interrupted; otherwise, the search is pursued along a new direction.

Geometrical considerations

When the pair (γ k , δ k ) used to update the matrix G k in the reduced quasi-Newton algorithm is defined by (2.3), one has

γ k = 1 0 g ′ (x k + αt k )Z -(x k ) dα δ k .
Therefore, (γ k , δ k ) collects information from a matrix close to g ′ (x)Z -(x). This section addresses the question of the symmetry of this matrix. This is a very desirable property. Indeed, the pair (γ k , δ k ) is used to update the symmetric matrix G k , so that it is better when (γ k , δ k ) does not contain useless information from an unsymmetric matrix.

The matrix g ′ (x)Z -(x) is clearly symmetric at the solution, because by formula

(2.2), g ′ (x * )Z -(x * ) = Z -(x * ) ⊤ L * Z -(x * )
is the reduced Hessian of the Lagrangian. But away from x * , this matrix may not be symmetric. An equivalent condition to the symmetry of g ′ (x)Z -(x) is given in Proposition 3.1 below. It uses the following definitions. For further details on the geometrical concepts used in this section, we refer the reader to [START_REF] Spivak | A Comprehensive Introduction to Differential Geometry[END_REF][START_REF] Boothby | An Introduction to Differentiable Manifolds and Differential Geometry[END_REF][START_REF] Conlon | Differentiable Manifolds -A first Course[END_REF][START_REF] Carmo | Riemannian Geometry[END_REF].

A vector field on R n is a smooth map X : R n → R n . The Lie bracket of two vector fields X and Y on R n is the vector field on R n , denoted by [X, Y ] and defined by

[X, Y ](x) = Y ′ (x) • X(x) -X ′ (x) • Y (x).
In this formula, Y ′ (x) • X(x) denotes the usual directional derivative at x of the function Y in the direction X(x). It is not difficult to show that if X and Y are tangent to a submanifold of R n in an open neighborhood of a point x, then [X, Y ] is also tangent to this submanifold in this neighborhood. Suppose that the function c : R n → R m defining the constraints in (1.7) is smooth and has surjective Jacobian matrices A(x), for all x ∈ R n . Then, the set

M x := {y ∈ R n : c(y) = c(x)} is a submanifold of R n . Let be given a smooth map x ∈ R n → Z -(x) ∈ R n×(n-m)
such that the columns of Z -(x) span the null space of A(x) (see Remark 1.1). We denote by Z - ,k the vector field defined by the kth column of Z -. These vector fields are tangent to the manifolds M x , so that their Lie brackets [Z - ,k , Z - ,l ] are also tangent to these manifolds. In particular, g ′ (x)Z -(x) is symmetric when g(x) = 0.

Proof. Denoting by {e k } the canonical basis of R n-m , we have

e ⊤ k g ′ (x)Z -(x)e l = g ′ k (x) • Z - ,l (x). Using g k (x) = Z - ,k (x) ⊤ ∇f (x), e ⊤ k g ′ (x)Z -(x)e l = ((Z - ,k ) ′ (x) • Z - ,l (x)) ⊤ ∇f (x) + Z - ,k (x) ⊤ ∇ 2 f (x)Z - ,l (x).
Therefore g ′ (x)Z -(x) is symmetric if and only if (3.1) holds. When g(x) = 0, ∇f (x) is perpendicular to the tangent space to M x at x. Since [Z - ,k , Z - ,l ](x) belongs to that space, (3.1) holds and g ′ (x)Z -(x) is symmetric. ✷

Since at a non stationary point, ∇f is arbitrary, the only way to be sure to have ∇f (x) perpendicular to [Z - ,k , Z - ,l ](x) is to choose tangent bases Z -with zero Lie bracket columns. As shown in Proposition 3.2, this last property is equivalent to the existence of local parametrizations of the manifolds M x , whose derivatives give Z -. Before making this statement precise, we recall the definition of a parametrization.

A parametrization of the manifold M x around y ∈ M x is a map ψ : U → M x defined on an open set U of R n-m such that y ∈ ψ(U ), ψ : U → ψ(U ) is a homeomorphism when ψ(U ) is endowed with the topology induced from that of R n , and i • ψ (i denotes the canonical injection from M x to R n ) is smooth and has injective derivatives (i • ψ) ′ (u) for all u ∈ U . Proposition 3.2. Consider the vector fields on M x given by the columns of the tangent bases Z -. The Lie brackets [Z - ,k , Z - ,l ](y) = 0 for k, l ∈ {1, . . . , n -m} and all y ∈ M x in a neighborhood of x if and only if there exists a parametrization ψ : U → M x of M x around x such that for all u ∈ U , (i

• ψ) ′ (u) = Z -(ψ(u)) (i denotes the canonical injection from M x to R n ).
It is easy to see that the if-part of the proposition holds. For this, assume the existence of a parametrization ψ with (i • ψ) ′ (u) = Z -(ψ(u)) and let k, l ∈ {1, . . . , n -m}. Then, for u ∈ U and y = ψ(u), one has

(i • ψ) ′ (u) • e k = Z - ,k (ψ(u)), (i • ψ) ′′ (u) • (e k , e l ) = (Z - ,k ) ′ (y) • Z - ,l (y). Therefore [Z - ,k , Z - ,l ](y) = 0 by the symmetry of (v 1 , v 2 ) → (i • ψ) ′′ (u) • (v 1 , v 2 )
. The only-if-part is more involved. For a proof, see [START_REF] Conlon | Differentiable Manifolds -A first Course[END_REF]Theorem 4.3.1].

For the interpretation of the search algorithm described in Section 4, we will need one of the two equivalent statements of Proposition 3.2, which we quote in the form of an assumption.

Assumption 3.3 (zero Lie bracket bases).

The tangent bases Z -are such that for all x ∈ R n there exists a parametrization ψ : U → M x around x such that (i • ψ) ′ (u) = Z -(ψ(u)) for all u ∈ U (i denotes the canonical injection from M x to R n ). Equivalently, the Lie brackets [Z - ,k , Z - ,l ](x) = 0 for k, l ∈ {1, . . . , n -m} and x ∈ R n .

Combining Propositions 3.1 and 3.2, we see that when Assumption 3.3 holds, the matrix g ′ (x)Z -(x) is symmetric for all x ∈ R n . In fact this statement is also a consequence of the following simple calculation.

Proposition 3.4. Suppose that ψ : U ⊂ R n-m → M x is a parametrization of M x around x such that (i • ψ) ′ (u) = Z -(ψ(u)) for all u ∈ U . Then, for u ∈ U and y = ψ(u), one has g ′ (y)Z -(y) = ∇ 2 (f • ψ)(u).
Proof. By definition of g, we have for u ∈ U , y = ψ(u), and

v 1 ∈ R n-m (f • ψ) ′ (u) • v 1 = f ′ (ψ(u)) • Z -(ψ(u))v 1 = g(ψ(u)) ⊤ v 1 .
Differentiating again in u in the direction

v 2 ∈ R n-m gives (f • ψ) ′′ (u) • (v 1 , v 2 ) = v ⊤ 1 g ′ (y)Z -(y)v 2 .
Hence the result. ✷

The question that arises now is whether the usual ways of computing tangent bases give rise to matrices Z -satisfying Assumption 3.3. The answer is positive for the tangent bases obtained by partitioning the Jacobian A and is usually negative for orthonormal bases, including those obtained by the QR factorization of A ⊤ . Let us consider these two cases successively.

Since A(x) is surjective, it has m linearly independent columns, which, for simplicity, will be supposed to be the first ones. Then, the Jacobian of c can written

A(x) = ( B(x) N (x) ) ,
where B(x) is an order m nonsingular matrix. Clearly, the columns of the matrix

Z -(x) = -B(x) -1 N (x) I n-m (3.2)
form a basis of the null space of A(x). The situation is completely different when the columns of Z -(x) are chosen to be orthonormal, i.e., Z -(x) ⊤ Z -(x) = I n-m for all x, no matter how this matrix is computed. To explain this, we use an argument from Riemannian geometry (see for example [START_REF] Carmo | Riemannian Geometry[END_REF]). Let x be a point fixed in R n . With the bases Z -, one can associate a Riemannian structure on the manifold M x , by defining the inner product of two tangent vectors X = Z -(y)u and Y = Z -(y)v at y ∈ M x by g(X, Y ) = u ⊤ v. If Assumption 3.3 holds, it gives parametrizations ψ, which are local isometries between R n-m and (M x , g). Indeed, ψ * g is the inner product of R n-m :

(ψ * g)(u, v) := g(ψ * u, ψ * v) := g(Z -(y)u, Z -(y)v) = u ⊤ v,
where ψ * and ψ * denote the tangent and cotangent maps associated with ψ. A consequence of this is that the Riemannian curvature of (M x , g) is zero, as the one of R n-m . Now, if the columns of the matrices Z -(y) are orthonormal, the Riemannian structure of (M x , g) is also the one induced by R n on M x considered as a submanifold of R n , say (M x , g E ). Indeed, for this induced Riemannian structure, the inner product of X = Z -(y)u and Y = Z -(y)v is the one in R n :

g E (X, Y ) = X ⊤ Y = u ⊤ Z -(y) ⊤ Z -(y)v = u ⊤ v = g(X, Y ),
by the orthonormality of the columns of Z -(y). Therefore, always under Assumption 3.3, we have shown that the Riemannian curvature of (M x , g E ) is zero (it is the same as the one of (M x , g)). This property is satisfied by very particular submanifolds of R n . For example, it occurs when M x is an affine subspace of R n or a cylinder (the product S 1 × E of the sphere of dimension one S 1 and an affine subspace E). Therefore, when Z -has orthonormal columns, Assumption 3.3 is rarely satisfied and the matrix g ′ (x)Z -(x) is not necessarily symmetric.

The conclusion of this discussion is that when the pairs (γ k , δ k ) used to update the matrix G k are given by ( 2.3), one should use tangent bases Z -satisfying Assumption 3.3, in order to have g ′ (x)Z -(x) symmetric. In particular, formula (3.2) is suitable, but orthonormal bases are usually inadequate.

Finally, note that the question of the symmetry of g ′ (x)Z -(x) does not arise when one uses the pairs (γ k , δ k ) given by formula (2.5), since these pairs collect information from a matrix close to Z -(y) ⊤ ∇ 2 ℓ(y, λ)Z -(y), which is always symmetric.

We believe that the comments given in this section could explain the much better efficiency of the reduced quasi-Newton algorithm tested in [START_REF] Ch | On the realization of the Wolfe conditions in reduced quasi-Newton methods for equality constrained optimization[END_REF], which uses the pairs (γ k , δ k ) given by (2.3), when the bases satisfy the zero Lie bracket bases assumption 3.3. Further numerical tests should be necessary to confirm this impression.

PLS for reduced quasi-Newton algorithms 4.1 Longitudinal search

Suppose that Assumption 3.3 holds. Then, the parametrizations ψ can be used to interpret algorithms in R n as piecewise linear approximations of the image by ψ of algorithms defined in the reduced space R n-m . For example, let

ψ k : U → M x k be a local parametrization of the submanifold M x k around x k = ψ k (u k ), u k ∈ U , satisfying (i • ψ k ) ′ (u) = Z -(ψ k (u)), for all u ∈ U . A quasi-Newton algorithm in R n-m , whose search path has the form α → u k + αδ k (δ k = -G -1 k g(x k )), in R n-m , is transformed by ψ k , in the curve α → pk (α), solution of the differential equation p′ k (α) = Z -(p k (α))δ k pk (0) = x k .
The first order approximation of the path α → pk (α) is α → x k + αt k , where t k is the tangent step (1.15) of the reduced Hessian algorithm.

The interest of this interpretation is that we know that under mild assumptions the Wolfe conditions (1.2) and (1.6) can be satisfied in the linear space R n-m on the function (f • ψ k ) along the direction δ k . This implies that one can find a step-size α k > 0 such that (0 < ω 1 < ω 2 < 1):

(f • ψ k )(u k + α k δ k ) ≤ (f • ψ k )(u k ) + ω 1 α k (f • ψ k ) ′ (u k ) • δ k , (f • ψ k ) ′ (u k + α k δ k ) • δ k ≥ ω 2 (f • ψ k ) ′ (u k ) • δ k . (4.1)
Using the merit function Θ σ k defined by (1.10) and the properties of the parametrization ψ k , these conditions can be rewritten (note that c(p k (α)) remains constant):

Θ σ k (p k (α k )) ≤ Θ σ k (x k ) + ω 1 α k g(x k ) ⊤ δ k , (4.2)
g(p k (α k )) ⊤ δ k ≥ ω 2 g(x k ) ⊤ δ k . (4.3)
The last inequality is called the reduced curvature condition. This inequality is very interesting because it implies that (g(p k (α k ))g(x k )) ⊤ δ k > 0, which is the desired inequality γ ⊤ k δ k > 0, provided γ k is the change in the reduced gradient between the points x k and pk (α k ).

Let us stress the fact that the step-size α k satisfying (4.2) and (4.3) exists without Assumption 3.3, as this can be observed by using the same argument as in unconstrained optimization (assuming the boundedness of Θ σ k from below on the manifold M x k ). This assumption was only used to view the path α → pk (α) as the image by ψ k of a straight line in the reduced space.

Two operations are necessary to derive an implementable algorithm from the path pk . First a discretization of pk must be introduced. This one can be done such that a sufficient decrease of the merit function Θ σ k and the reduced curvature condition hold also along the discretized path. Next, a restoration step r k (see (1.14)) has to be added to complete the iteration, in order to force the decrease of the norm of the constraints. This approach was developed in [START_REF] Ch | Maintaining the positive definiteness of the matrices in reduced secant methods for equality constrained optimization[END_REF].

Longitudinal and transversal search

For some applications, the algorithm outlined in the last paragraph may be too expensive to use. Indeed, even if a single step is performed in the longitudinal part of the algorithm (formed by the discretization of the path pk ), the functions have to be linearized twice per iteration: at x k and x k + α k t k . In fact, one can view this algorithm as a globalization of the local reduced quasi-Newton algorithm, in which the matrix G k is updated with the pair (γ k , δ k ) defined by (2.3), requiring the additional linearization at the intermediate point

x k + t k .
A way of avoiding this additional linearization is to add a transversal component (in the range space of A -) to the derivative of pk and to define a new search path pk by p′

k (α) = Z -(p k (α))δ k -A -(p k (α))c(p k (α)) pk (0) = x k , (4.4)
where

δ k = -G -1 k g(x k
), as before. The first order approximation of this path is α → x k + α(t k + r 0 k ), where r 0 k is defined by (1.17). This is a simplified version of the reduced Hessian algorithm, r 0 k being used instead of r k . Under the assumptions g ′ A -≡ 0 and g ′ Z -is constant on the reduced gradient manifolds, one has g(p k (α)) = g(p k (α)), as long as both paths pk and pk exist [START_REF] Ch | On the realization of the Wolfe conditions in reduced quasi-Newton methods for equality constrained optimization[END_REF]. In this case, the same interpretation holds for a search along pk and pk : realizing

g(p k (α k )) ⊤ δ k ≥ ω 2 g(x k ) ⊤ δ k is equivalent to realizing (4.
1) in the reduced space. These conditions are used to give an interpretation to the search along pk in terms of a search along a straight line. Now they are not necessary to be able to find an adequate step-size along pk , as shown by the following proposition [START_REF] Ch | On the realization of the Wolfe conditions in reduced quasi-Newton methods for equality constrained optimization[END_REF]. Proposition 4.1. Suppose that the path α → pk (α) defined by (4.4) exists for a sufficiently large step-size α ≥ 0, that Θ σ k is bounded from below along this path, that σ k ≥ λ(p k (α)) D whenever pk (α) exists, and that ω 2 ∈ (0, 1). Then, the inequalities

Θ σ k (p k (α)) ≤ Θ σ k (x k ), g(p k (α)) ⊤ Z(x k )t k ≥ ω 2 g(x k ) ⊤ Z(x k )t k (4.5)
are satisfied for some α > 0.

Again, the reduced curvature condition (4.5) is very attractive, since after discretization of the search path, it provides the desired inequality γ ⊤ k δ k > 0, this time with γ k = g(x k+1 )g(x k ), hence without having to linearize the functions at the intermediate point x k + t k .

Like for the longitudinal search, it remains to discretize the path pk to obtain an implementable algorithm. Let us denote by

α 0 k = 0 < α 1 k < • • • < α i k k = α k the discretization step-
sizes, which are not given a priori but computed as explained below, and by x i k the points approximating pk (α i k ) (with

x 0 k = x k and x i k k = x k+1
). An explicit Euler approximation gives

x i+1 k = x i k + (α i+1 k -α i k )d i k , i = 1, . . . , i k -1,
where

d i k = Z -(x i k )δ k -A -(x i k )c(x i k ). For i = 1, . . . , i k -1, the step-size α i+1 k is determined such that Θ σ k decreases sufficiently: Θ σ k (x i k + (α i+1 k -α i k )d i k ) ≤ Θ σ k (x i k ) + ω 1 (α i+1 k -α i k )Θ ′ σ k (x i k ; d i k ).
The first intermediate point x 1 k is determined in a slightly different way, so that, when the unit step-size is accepted, x 1 k = x k + t k + r k and the local method is recovered, allowing the q-superlinear convergence of the algorithm; see [START_REF] Ch | Superlinear convergence of a reduced BFGS method with piecewise line-search and update criterion[END_REF] for the details. Once α i+1 k has been computed, the reduced curvature condition is tested at At the intermediate points x i k , the functions defining the problem have to be linearized. This may look expensive, but this impression have to be appreciated in view of the following observations. 1. Even when intermediate linearizations occur, the displacement along d i k helps in decreasing Θ σ k , so that a progress to the solution is done and this displacement is not useless. 2. The PLS is only used when the update criterion is verified (see Section 2.1), because there is no need to have γ ⊤ k δ k > 0 when the matrix G k is not updated. If an update criterion of the form (2.7) is used, the transversal part of the displacement is small, so that the search path pk is close to the manifold M x k . As for unconstrained optimization, we have observed that in practice the unit step-size is then accepted most of the time, so that no intermediate points are necessary. This fact is corroborated by the asymptotic result below (Theorem 4.2, outcome (ii)). To conclude this section, we mention a superlinear convergence result obtained with the algorithm described above, equipped with an update criterion that we describe now. This one is of the form (2.7) with a sequence {µ k } ruled by the algorithm itself (not a sequence given a priori as in [START_REF] Nocedal | Projected Hessian updating algorithms for nonlinearly constrained optimization[END_REF]). The idea behind this criterion is the following. The sequence µ k in (2.7) need not be decreased too rapidly, otherwise the matrix G k would be rarely updated, which would prevent superlinear convergence. This suggests not changing µ k , i.e., taking µ k+1 = µ k , when the update criterion does not hold at iteration k. This strategy looks safe in particular in the case when the update criterion is not satisfied after a given iteration, say k 0 . Then µ k is constant for k ≥ k 0 and t k = O( r k ). It can be shown that the latter estimate readily implies the superlinear convergence. More difficult is the case when G k is both updated infinitely often and kept unchanged infinitely often. For dealing with this situation, the following update criterion is suitable

x i+1 k : g(x i+1 k ) ⊤ δ k ≥ ω 2 g(x k ) ⊤ δ k . If it
c(x k ) + c(x k + t k ) ≤ µ e 1 k⊖2 δ k . (4.6)
In (4.6), µ > 0 is an arbitrary constant and e

1 k = x 1 k -x k is the step from x k to the first intermediate iterate x 1 k . The index k ⊖ 2 is used to express that µ k = µ e 1 k⊖2
will not change when there is no matrix update. More precisely, k ⊖ 1 is the greatest index less than k at which an update occurred. Hence, (k + 1) ⊖ 1 (the greatest index less than k + 1 at which an update occurred) is the same index as k ⊖ 1 if there is no matrix update at iteration k. For technical reasons, it is not the index k ⊖ 1 that must be used in (4.6) but k ⊖ 2 := (k ⊖ 1) ⊖ 1 (the greatest index less than k ⊖ 1 at which an update occurred). Suppose also that the reduced quasi-Newton algorithm with BFGS matrix updates, PLS, and the update criterion (4.6) described above generates a sequence of points {x i k } k≥1,0≤i≤i k -1 in Ω different from x * and converging to x * in the sense that

max 0≤i≤i k -1 x i k -x * → 0, when k → ∞.
Then, the following properties hold: (i) the sequences of matrices {G k } k≥1 and {G -1 k } k≥1 are bounded; (ii) the ideal step-size is accepted eventually: i k = 1 and α k = 1 for k large; (iii) the sequence {x k } k≥1 converges q-superlinearly in two steps to x * ; (iv) the sequence {x k + t k } k≥1 converges q-superlinearly to x * . This result given in [START_REF] Ch | Superlinear convergence of a reduced BFGS method with piecewise line-search and update criterion[END_REF] is, to our knowledge, the first extension of a similar result proved by Powell [START_REF] Powell | Some global convergence properties of a variable metric algorithm for minimization without exact line searches[END_REF] for the BFGS algorithm in unconstrained optimization, at least with hypotheses as weak as those listed above. This is not an obvious result. For example, it may occur that the matrices G k are not updated infinitely often, although the sequence {x k + t k } still converges superlinearly. This theorem shows that q-superlinear convergence with quasi-Newton methods for minimizing equality constrained problems is possible with a single linearization of the constraints per iteration (for a result with two linearizations per iteration, see [START_REF] Byrd | An analysis of reduced Hessian methods for constrained optimization[END_REF]) and without knowing the threshold value of the penalty parameter τ making the Hessian of the augmented Lagrangian positive definite (a superlinear result with this hypothesis for the SQP method is proved in [START_REF] Byrd | An SQP augmented Lagrangian BFGS algorithm for constrained optimization[END_REF]).

PLS for quasi-Newton-SQP algorithms

We have seen in Section 1.3 that, when (γ k , δ k ) ∈ R 2n is given by (2.9), the penalty parameter τ can be used to ensure the positive definiteness of M k+1 (by means of γ ⊤ k δ k > 0), the matrix approximating the Hessian of the augmented Lagrangian, provided δ k = x k+1x k is not in the tangent space to the constraint manifold. A back-up strategy is proposed in [START_REF] Byrd | An SQP augmented Lagrangian BFGS algorithm for constrained optimization[END_REF] to deal with the case when A(x k )δ k is relatively near zero. In this case, δ k is almost tangent to the constraint manifold and we have seen, with the longitudinal search of Section 4, that the PLS is appropriate to find a longitudinal displacement that can ensure the positive definiteness of matrices approximating Z -(x * ) ⊤ L * Z -(x * ). Here, we follow [START_REF] Armand | A piecewise line-search technique for maintaining the positive definiteness of the updated matrices in the SQP method[END_REF] and show how to combine the two ideas: controlling the "transversal part" of M k+1 by τ and its "longitudinal part" by the PLS algorithm.

The description of this approach requires some preliminaries. First, the presentation is simplified if we introduce the right inverse A - k (x) of A(x) associated with the quadratic program (1.9). Assume that the matrix M k ∈ R n×n is positive definite in the null space of A(x). Then, the quadratic program

min 1 2 d ⊤ M k d c + A(x)d = 0.
(5.1) has a unique solution, which is given by the optimality conditions of (5.1): for some

λ ∈ R m , M k d + A(x) ⊤ λ = 0 A(x)d = -c. (5.2)
This shows that the solution d is a linear function of the vector c ∈ R m . We denote by -A - k (x) the n × m matrix representing this linear function (it depends on M k , hence the index k in A - k (x)). Then, the solution of (5.1) can be written -A - k (x)c. Substituting this quantity in the second equation of (5.2) shows that A - k (x) is a right inverse of A(x). Suppose now that a tangent basis Z -(x) is given and that we denote by Z k (x) the unique (nm) × n matrix satisfying

A - k (x)A(x) + Z -(x) Z k (x) = I n .
Multiplying the first equation of (5.2) by Z -(x) ⊤ leads to

Z -(x) ⊤ M k A - k (x) = 0,
so that the solution of (1.9) can be written

d k = Z -(x k ) Z k (x k )d k -A - k (x k )c(x k ), where Z k (x k )d k = -(Z -(x k ) ⊤ M k Z -(x k )) -1 g(x k ).
It is also useful to compute the derivative of the multiplier estimate λ defined by (1.18). One has λ

(x) = -A -(x) ⊤ ∇f (x) = -A -(x) ⊤ ∇ℓ(x, λ * ) + λ * , so that λ ′ (x * ) = -A -(x * ) ⊤ L * . (5.3)
We are now ready to present the approach of [START_REF] Armand | A piecewise line-search technique for maintaining the positive definiteness of the updated matrices in the SQP method[END_REF]. Since M k+1 has to approximate the Hessian of the augmented Lagrangian, (γ k , δ k ) should verify

γ k ≃ L * δ k + τ A(x k ) ⊤ A(x k )δ k ,
at the first order (see (2.1) and (2.8)). The first term of the sum in the right hand side is approximated, in (2.9), by the change in the gradient of the Lagrangian. Here, we split this term in two parts, for a reason that will be clear below. Using (1.16), (2.2), and (5.3):

L * δ k = Z(x k ) ⊤ Z -(x k ) ⊤ L * δ k + A(x k ) ⊤ A -(x k ) ⊤ L * δ k ≃ Z(x k ) ⊤ (g ′ (x * ) • δ k ) -A(x k ) ⊤ (λ ′ (x * ) • δ k ) ≃ Z(x k ) ⊤ (g(x k+1 )-g(x k )) -A(x k ) ⊤ (λ(x k+1 )-λ(x k )), provided δ k ≃ x k+1 -x k .
This approximate computation leads us to the following formula:

γ k = Z(x k ) ⊤ (g(x k+1 )-g(x k )) -A(x k ) ⊤ (λ(x k+1 )-λ(x k )) + τ A(x k ) ⊤ A(x k )δ k ,
where τ has to be adjusted. The inner product of γ k with δ k is

γ ⊤ k δ k = (g(x k+1 )-g(x k )) ⊤ Z(x k )δ k -(λ(x k+1 )-λ(x k )) ⊤ A(x k )δ k + τ A(x k )δ k 2 2 .
Like for formula (2.9), one can make γ ⊤ k δ k positive by taking τ sufficiently large, provided A(x k )δ k = 0. Recall that, when δ k is tangent, the approach in [START_REF] Byrd | An SQP augmented Lagrangian BFGS algorithm for constrained optimization[END_REF] requires a back-up strategy. Here, when A(x k )δ k = 0,

γ ⊤ k δ k = (g(x k+1 ) -g(x k )) ⊤ Z(x k )δ k , (5.4)
so that the positivity of γ ⊤ k δ k can be obtained by choosing the next iterate such that g(x k+1 ) ⊤ Z(x k )δ k > g(x k ) ⊤ Z(x k )δ k . This reminds us of the reduced curvature inequality (4.5) satisfied by the PLS in reduced quasi-Newton method. Here, this inequality may not be feasible, because Z(x k )δ k may not be a reduced descent direction, meaning that g(x k ) ⊤ Z(x k )δ k may not be negative when A(x k )δ k = 0. Therefore, we ask instead to realize the following reduced curvature condition

g(x k+1 ) ⊤ Z k (x k )d k ≥ ω 2 g(x k ) ⊤ Z k (x k )d k , where d k is the SQP direction. Note that Z(x k )d k = Z k (x k )d k when A(x k )d k = 0, so that γ ⊤ k δ k > 0 when A(x k )δ k = 0 and δ k is parallel to d k .
The following result is similar to proposition 4.1.

Proposition 5.1. Suppose that the path α → p k (α) defined by

p ′ k (α) = Z -(p k (α)) Z k (x k )d k -A - k (p k (α))c(p k (α)) p k (0) = x k
exists for a sufficiently large step-size α ≥ 0. Suppose also that Θ σ k is bounded from below along this path, that σ k ≥ λ QP k (p k (α)) D whenever p k (α) exists, that M k is positive definite, and that ω 2 ∈ (0, 1). Then, the inequalities

Θ σ k (p k (α)) ≤ Θ σ k (x k ) g(p k (α)) ⊤ Z k (x k )d k ≥ ω 2 g(x k ) ⊤ Z k (x k )d k
are satisfied for some α > 0.

Like for the reduced quasi-Newton algorithm, an explicit Euler discretization of the path p k is introduced:

α 0 k = 0 < α 1 k < • • • < α i k k = α k
are the discretization step-sizes determined as explained below and the points x i k are approximations of p k (α i k ) (with x 0 k = x k and x i k k = x k+1 ) given by

x i+1 k = x i k + (α i+1 k -α i k )d i k , i = 1, . . . , i k -1,
where

d i k = Z -(x i k ) Z k (x k )d k -A -(x i k )c(x i k ).
For i = 0, . . . , i k -1, the step-size α i+1 k is determined such that Θ σ k decreases sufficiently (σ i k may need to be adapted before the determination of α i+1 k , in order to make d i k a descent direction):

Θ σ i k (x i k + (α i+1 k -α i k )d i k ) ≤ Θ σ i k (x i k ) + ω 1 (α i+1 k -α i k )Θ ′ σ i k (x i k ; d i k ).
Once α i+1 k has been computed, the reduced curvature condition is tested at x i+1 k :

g(x i+1 k ) ⊤ Z k (x k )d k ≥ ω 2 g(x k ) ⊤ Z k (x k )d k .
(5.5)

If it holds the search is completed, otherwise the search is pursued along the new intermediate direction d i+1 k . It can be shown that this PLS algorithm terminates in a finite number of trials.

Let us describe more precisely how δ k is computed. On the one hand, we have said that it is necessary to take δ k ≃ x k+1x k . On the other hand, one needs to control the positivity of γ ⊤ k δ k , even when A(x k )d k = 0. In our algorithm, γ ⊤ k δ k is given by (5.4) and the positivity of this inner product is not guaranteed by the reduced curvature condition (5.5) when δ k = x k+1x k , since this vector may not be parallel to d k . Therefore, we prefer to take the following approximation of x k+1 -x k :

δ k = α k Z -(x k ) Z k (x k )d k -α A k A - k (x k )c(x k ), (5.6)
where

α A k = i k -1 i=0 (α i+1 k -α i k )e -α i k .
This formula aims at taking into account the fact that the value of c at x i k is used in the search directions d i k , while c(x k ) is used in δ k . It is based on the observation that along the path p k , we have c(p k (α)) = e -α c(x k ), which gives after discretization: c(x i k ) ≃ e -α i k c(x k ). Let us now check that the form (5.6) of δ k is appropriate. Suppose that A(x k )δ k = 0. Then, c(x k ) = 0, δ k = α k Z -(x k ) Z k (x k )d k , and

γ ⊤ k δ k = (g(x k+1 )-g(x k )) ⊤ Z(x k )δ k = α k (g(x k+1 )-g(x k )) ⊤ Z k (x k )d k > 0,
by the reduced curvature condition (5.5), which is satisfied for i = i k -1 (in which case x i+1 k = x k+1 ). The conclusion of this discussion is that for any k ≥ 1, one can find a (finite) τ ≥ 0 such that γ ⊤ k δ k > 0, either because A(x k )δ k = 0 or because A(x k )δ k = 0 and γ ⊤ k δ k > 0 by the reduced curvature condition (5.5). More details, including a global convergence result and numerical experiment, are given in [START_REF] Armand | A piecewise line-search technique for maintaining the positive definiteness of the updated matrices in the SQP method[END_REF].

Conclusion

By way of conclusion, we list some open problems and questions raised in this paper.

First, it would be very useful to have an update criterion based on the difference between g(x k+1 ) and g(x k +t k ) instead of the criterion (2.7) based on the comparison between the norm of the tangent step t k and the norm of the restoration step r k , as explained in Section 1. [START_REF] Armijo | Minimization of functions having Lipschitz continuous first partial derivatives[END_REF].

Another problem that would deserve more attention is to give a geometrical interpretation to the piecewise line-search introduced in Section 4 and 5, in the case when Assumption 3.3 does not hold. The PLS can still find a step-size satisfying reduced curvature conditions but its interpretation as a usual search along a straight line in the reduced space no longer holds.

As shown by the discussion in Section 3, γ k given by formula (2.5) or the second formula in (2.6) may be more appropriate, in particular when the tangent bases are chosen orthonormal. In this case, can we still use the PLS approach presented here for a γ k given by (2.3) or the first formula in (2.6)?

The q-superlinear convergence of the BFGS version of the SQP algorithm is still an open problem, for hypotheses similar to those of Theorem 4.2. Can a technique ensuring γ ⊤ k δ k > 0 at each iteration, such as the PLS technique, be helpful in getting this result?

Proposition 3 . 1 .

 31 The matrix g ′ (x)Z -(x) is symmetric if and only if ∇f (x) ⊥ [Z - ,k , Z - ,l ](x), for all k, l ∈ {1, . . . , n -m}.(3.1)

Proposition 3 . 5 .

 35 The tangent bases given by (3.2) satisfy Assumption 3.3. In particular, for these bases, g ′ (x)Z -(x) is a symmetric matrix.Proof. Let us partition the components of y ∈ R n in (ξ, u), with ξ ∈ R m and u ∈ R n-m . Then, by the implicit function theorem, the nonsingularity of B(x) implies the existence an implicit function u ∈ U → ξ(u) ∈ R m , where U is an open set of R n-m , such that c(ξ(u), u) = c(x) for all u ∈ U . The parametrization satisfying Assumption 3.3 is ψ(u) = (ξ(u), u). Then, by proposition 3.4, g ′ (x)Z -(x) is symmetric. ✷

  holds the search is completed, otherwise the search is pursued along the new intermediate direction d i+1 k . It can be shown that this PLS algorithm terminates in a finite number of trials.

Theorem 4 . 2 .

 42 Let (x * , λ * ) be a primal-dual solution of problem (1.7) such that c(x * ) = 0, g(x * ) = 0, and Z -(x * ) ⊤ L * Z -(x * ) is positive definite. Suppose that there is an open convex neighborhood Ω of x * , such that: (a) f and c are twice continuously differentiable on Ω; (b) the Jacobian matrix A(x) of c is surjective on Ω and the map x → (A -(x), Z -(x)) is Lipschitz continuous on Ω; (c) g is differentiable on Ω with Lipschitz continuous derivative.

  [START_REF] Boothby | An Introduction to Differentiable Manifolds and Differential Geometry[END_REF] and to update the matrix G k only when the update with this new γ k looks safe. The difficulty with these values of γ k is that the estimate (2.1) no longer holds. To measure the appropriateness of an update, a possibility is to compare the length of the tangent step t k to the length of the restoration step r k . An update can occur if the update criterionr k ≤ µ k t k (2.7)holds[START_REF] Nocedal | Projected Hessian updating algorithms for nonlinearly constrained optimization[END_REF][START_REF] Ch | Sur quelques problèmes d'identification et d'optimisation rencontrés en physique des plasmas[END_REF]. The sequence {µ k } in (2.7) has to converge to zero if G k is updated infinitely often. The aim of (2.7) is to select the iterations where the cheap γ k (given by (2.6)) and the safe γ k (given by(2.3) or (2.5)) are similar because x k + t k and x k+1 are closer and closer to each other relatively to the distance between x k and x k + t k . It is not difficult to see that when (2.7) holds with µ k → 0, γ k given by one of the formulas in (2.6) and δ k = Z(x k )t k satisfy (2.1). The update criterion (2.7) works well in theory

† INRIA Rocquencourt, BP 105, 78153 Le Chesnay Cedex (France) ; e-mail : Jean-Charles. Gilbert@inria.fr.