# A piecewise line-search technique for maintaining the positive definiteness of the matrices in the SQP method

Paul Armand, Jean Charles Gilbert

<figure><img src="image_1.png" /></figure>

## To cite this version:

Paul Armand, Jean Charles Gilbert. A piecewise line-search technique for maintaining the positive definiteness of the matrices in the SQP method. Computational Optimization and Applications, 2000, 16, pp.121-158. 10.1023/A:1008701324575 . inria-00074070

HAL Id: inria-00074070

https://inria.hal.science/inria-00074070

Submitted on 24 May 2006

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
A piecewise line-search technique
for maintaining the positive definiteness
of the matrices in the SQP method

by

Paul ARMAND† and Jean Charles GILBERT‡

Research Report n° 2615

July 17, 1995

ABSTRACT

A technique for maintaining the positive definiteness of the matrices in the quasi-Newton version of the SQP algorithm is proposed. In our algorithm, approximations of the Hessian of the augmented Lagrangian are updated. The positive definiteness of these matrices in the space tangent to the constraint manifold is ensured by a piecewise line-search technique, while their positive definiteness in a decoupled complementary subspace is obtained by setting the augmentation parameter. The combination of these two ideas makes the new approach more robust in our experiment with respect to existing approaches.

**Key words:** BFGS formula, equality constrained optimization, piecewise line-search, quasi-Newton algorithm, successive quadratic programming.

**Abbreviated title:** Piecewise line-search for SQP.

**AMS Subject Classification:** primary: 65K05; secondary: 49M37, 90C30.

<footnote>†LACO - Département de Mathématiques - 123, av. A. Thomas - 87060 Limoges Cedex (France).</footnote>

<footnote>‡Projet PROMATH, INRIA Rocquencourt, BP 105, 78153 Le Chesnay Cedex (France).</footnote>
Une méthode de recherche linéaire brisée
pour maintenir la définie positivité
des matrices dans la méthode PQS

par

Paul ARMAND† et Jean Charles GILBERT‡

Rapport de Recherche n° 2615
17 juillet 1995

RÉSUMÉ

Une technique pour maintenir la définie positivité des matrices dans la version quasi-newtonienne de la PQS est proposée. Dans notre algorithme, ce sont des approximations du hessien du lagrangien augmenté qui sont mises à jour. Leur définie positivité dans le plan tangent à la variété des contraintes est assurée par une technique de recherche linéaire brisée, tandis que leur définie positivité dans un espace supplémentaire découplé du premier est obtenue par le réglage du paramètre d'augmentation. La combinaison de ces deux idées rend la nouvelle approche plus robuste dans nos expériences, quand on la compare à d'autres approches existantes.

<footnote>†LACO - Département de Mathématiques - 123, av. A. Thomas - 87060 Limoges Cedex (France).</footnote>

<footnote>‡Projet PROMATH, INRIA Rocquencourt, BP 105, 78153 Le Chesnay Cedex (France).</footnote>
# 1 Introduction

We consider with a numerical point of view the nonlinear equality constrained optimization problem

$$ \begin{array}{ll} \min & f(x) \\ \text{subject to} & c(x) = 0, \quad x \in \Omega, \end{array} \qquad (1.1) $$

where $\Omega$ is an open set of $\mathbb{R}^n$ and the two functions $f: \Omega \to \mathbb{R}$ and $c: \Omega \to \mathbb{R}^m$ ($1 \le m < n$) are sufficiently smooth. Since $\Omega$ is supposed open, this set does not define general constraints. It is just the set on which good properties of $f$ and $c$ hold. For example, we assume that the $m \times n$ Jacobian matrix of $c$ at $x$, denoted by $A(x) = \nabla c(x)^\top$, has full rank on $\Omega$.

The Lagrangian function associated to problem (1.1) is defined on $\Omega \times \mathbb{R}^m$ by

$$ \ell(x, \lambda) = f(x) + \lambda^\top c(x). \qquad (1.2) $$

The vector $\lambda$ is called the *Lagrange multiplier*. The first order optimality conditions of problem (1.1) at a local solution $x_*$ with associate multiplier $\lambda_*$ can be written

$$ \nabla f_* + A_*^\top \lambda_* = 0 \quad \text{and} \quad c_* = 0. \qquad (1.3) $$

Throughout the paper, the notation $f_* = f(x_*)$, $\nabla f_* = \nabla f(x_*)$, $A_* = A(x_*)$, etc, is used and should be unambiguous from the context. It is also assumed that the second order sufficient conditions of optimality hold at a local solution of problem (1.1). Using the notation $L_* = \nabla^2_{xx}\ell(x_*, \lambda_*)$, this can be written

$$ h^\top L_* h > 0, \quad \text{for all } h \neq 0 \text{ such that } A_* h = 0. \qquad (1.4) $$

The sequential quadratic programming (SQP) algorithm is a Newton-like method in $(x, \lambda)$ applied to the first order optimality conditions (1.3) (see for example Fletcher [8] or the recent survey paper [1]). The $k$th iteration of the algorithm can be described as follows. Given an iterate pair $(x_k, \lambda_k) \in \Omega \times \mathbb{R}^m$, the following quadratic subproblem in $d$ is solved:

$$ \begin{array}{ll} \min & \nabla f_k^\top d + \frac{1}{2} d^\top M_k d \\ \text{s.t.} & c_k + A_k d = 0. \end{array} \qquad (1.5) $$

We adopt the notation $f_k = f(x_k)$, $\nabla f_k = \nabla f(x_k)$, $c_k = c(x_k)$, $A_k = A(x_k)$, etc. In (1.5), it is suitable to take for $M_k$ the Hessian of the Lagrangian or an approximation to it. Let us denote by $(d_k, \lambda_k^{\text{QP}})$ a primal-dual solution of (1.5), i.e., a solution of its optimality conditions

$$ \begin{array}{l} \nabla f_k + M_k d_k + A_k^\top \lambda_k^{\text{QP}} = 0 \\ c_k + A_k d_k = 0. \end{array} \qquad (1.6) $$

The link between (1.3) and (1.5) is that, when $M_k$ is the Hessian of the Lagrangian, $(d_k, \lambda_k^{\text{QP}} - \lambda_k)$ is the Newton step for the system (1.3) at $(x_k, \lambda_k)$.

The convergence of this algorithm from remote starting points is often obtained by using $d_k$ as search direction, along which a stepsize $\alpha_k > 0$ is chosen. The stepsize is adjusted such that the next iterate

$$ x_{k+1} = x_k + \alpha_k d_k, $$
reduces sufficiently the value of some merit function.

In the quasi-Newton version of the method, $M_k$ approximates some Hessian and is commonly forced to be positive definite to ensure descent properties. This can be achieved by using the BFGS formula for instance: for some vectors $\gamma_k$ and $\delta_k$ in $\mathbb{R}^n$,

$$M_{k+1} = M_k - \frac{M_k \delta_k \delta_k^\top M_k}{\delta_k^\top M_k \delta_k} + \frac{\gamma_k \gamma_k^\top}{\gamma_k^\top \delta_k}. \quad (1.7)$$

With this formula, it is well known that the positive definiteness is sustained from $M_k$ to $M_{k+1}$ if and only if the following curvature condition holds:

$$\gamma_k^\top \delta_k > 0. \quad (1.8)$$

When $M_k$ is taken as an approximation of the Hessian of the Lagrangian, it makes sense to take for $\gamma_k$ in (1.7) the vector

$$\gamma_k^\ell = \nabla_x \ell(x_{k+1}, \lambda) - \nabla_x \ell(x_k, \lambda), \quad (1.9)$$

where $\lambda$ is a some multipier, usually $\lambda_k^{\text{QP}}$. However, the lack of positive definiteness of the Hessian of the Lagrangian function at $(x_*, \lambda_*)$ makes this approach difficult. Indeed, with this choice of $\gamma_k$, the curvature condition may never be realized for any displacement $\delta_k = x_{k+1} - x_k$ along $d_k$, because the Lagrangian function may have negative curvature along this direction, even close to the solution.

The idea of modifying the vector $\gamma_k^\ell$ to force satisfaction of the curvature condition goes back at least to Powell [20], who suggested to set $\gamma_k$ to a convex combination of $\gamma_k^\ell$ and $M_k \delta_k$:

$$\gamma_k^P = \theta \gamma_k^\ell + (1 - \theta) M_k \delta_k, \quad (1.10)$$

where $\theta$ is the number in (0, 1], the closest to 1, such that the inequality

$$\gamma_k^{\text{P}\top} \delta_k \geq \eta \delta_k^{\text{T}} M_k \delta_k$$

is satisfied. The constant $\eta$ is set to 0.2 in [20] and to 0.1 in [22]. Powell's correction of $\gamma_k^\ell$ is certainly the most widely used technique in practice. Its success is due to its appealing simplicity and its usually good numerical performance. The fact that it may encounter difficulties partly motivates this study (see [22] or [23, p. 125]). Another motivation is that the best known result obtained so far on the speed of convergence with Powell's correction (namely, the r-superlinear convergence, see [21]) is not as good as one can reasonably expect, which is the q-superlinear convergence.

Another modification of $\gamma_k^\ell$ is to take for vector $\gamma_k$ an approximation of the change in the gradient of the augmented Lagrangian, which is the function (1.2) with the augmentation term $\frac{\pi}{2}\|c(x)\|^2$ ($\| \cdot \|$ denotes the $\ell_2$-norm). This idea, proposed by Tapia [28], has roots in the work of Han [14] and Tapia [27] and was refined later by Byrd, Tapia, and Zhang [4]. In this approach, the matrix $M_k$ is viewed as an approximation of the Hessian of the augmented Lagrangian function and $\gamma_k$ is set to

$$\gamma_k^S = \gamma_k^\ell + r_k A_{k+1}^\top A_{k+1} \delta_k, \quad (1.11)$$
where the augmentation parameter $r_k$ is the smallest nonnegative number satisfying

$$ \gamma_k^{\mathrm{T}} \delta_k \geq \max\{|\gamma_k^{\mathrm{T}} \delta_k|, \nu_{\mathrm{BTZ}} \|A_{k+1} \delta_k\|^2\}. \quad (1.12) $$

The positive constant $\nu_{\mathrm{BTZ}}$ is set to 0.01 in [4]. It is clear that this strategy does not work when $A_{k+1} \delta_k = 0$ and $\gamma_k^{\mathrm{T}} \delta_k \le 0$. In this case, the authors propose the following “back-up strategy”. When (1.12) does not hold with $r_k = 0$ and

$$ \|A_{k+1} \delta_k\| < \min\{\beta_{\mathrm{BTZ}}, \|\delta_k\|\} \|\delta_k\| \quad (1.13) $$

where $\beta_{\mathrm{BTZ}}$ is a small positive number (the value 0.01 is proposed in [4]), then the vector $A_{k+1}^\top A_{k+1} \delta_k$ in (1.11) is replaced by $\delta_k$ and $r_k$ is set such that (1.12) is satisfied. Numerical experiment in [4] shows that Byrd, Tapia, and Zhang’s strategy is numerically competitive with Powell’s correction.

Let us mention that there are techniques to maintain the positive definiteness of the updated matrices that do not modify the vector $\gamma_k^\ell$: see Coleman and Fenyes [6].

The present paper is mainly motivated by the desire to realize the curvature condition (1.8). To carry out this task, we propose a new update scheme that is based on two principles. On the one hand, as in the work of Byrd, Tapia, and Zhang, we take the point of view to force $M_k$ to be an approximation of the Hessian of the augmented Lagrangian. This is achieved by taking for $\gamma_k$ the following value:

$$ \gamma_k = \tilde{\gamma}_k + r_k A_k^\top A_k \delta_k, \quad \text{with} \quad \tilde{\gamma}_k \approx L_* \delta_k. $$

The precise form of $\gamma_k$ will be given in Section 3.1. As we have seen, the strategy of forcing the positivity of $\gamma_k^\ell \delta_k$ by tuning $r_k$ is not appropriate when $A_k \delta_k$ is too close to zero and $\tilde{\gamma}_k^\ell \delta_k$ is negative. Instead of using a back-up strategy as Byrd, Tapia, and Zhang, we observe that when $A_k \delta_k = 0$ and $\delta_k$ is parallel to the SQP direction, this direction is tangent to the constraint manifold. In this case, the piecewise line-search (PLS) technique introduced by Gilbert [11, 12] in the framework of reduced quasi-Newton methods, is well adapted to making the term $\tilde{\gamma}_k^\ell \delta_k$ positive. Our algorithm is based on a combination of these two ideas.

The PLS technique aims at generalizing what is done in unconstrained optimization, where the curvature condition is fulfilled by a line-search algorithm realizing the Wolfe conditions (see for instance [16, Chap. II, § 3.3]). When constraints are present, the problem is more difficult, since the curvature condition may never hold along a straight line. However, there is a path defined by a particular differential equation along which such condition can be realized. The technique consists then in following a piecewise linear approximation of this “guiding path”, each discretization point being successively chosen by means of an Armijo line-search procedure, until a Wolfe point is found. In the present case, this technique allows us to choose a suitable value $r_k$ such that the curvature condition is satisfied, even when $A_k \delta_k = 0$.

The paper is organized as follows. In Section 2, we make precise our notation, the form of the SQP direction, and our choice of merit function. Section 3 presents our approach to satisfy the curvature condition (1.8) and outlines the PLS technique. Section 4 shows the finite termination of the search algorithm. The overall minimization algorithm is given in Section 5 and its convergence is proved. Section 6 gives
more details on implementation issues and Section 7 relates numerical tests. The paper terminates with a conclusion section.

## 2 Background material and notation

Let us first introduce two decompositions of $\mathbb{R}^n$ that will be useful throughout the paper. Each of them decomposes the variable space in two complementary subspaces and is characterized by a triplet $(Z^{-}(x), A^{-}(x), Z(x))$. The columns of the matrices $Z^{-}(x)$ and $A^{-}(x)$ span the two complementary subspaces and $Z(x)$ is deduced from $Z^{-}(x), A^{-}(x)$, and $A(x)$.

In the first decomposition, $Z^{-}(x)$ and $A^{-}(x)$ are typically given by the user. These operators and $Z(x)$ have to satisfy the following properties.

* $Z^{-}(x)$ is an $n \times (n-m)$ matrix, whose columns form a basis of the null space $\mathcal{N}(A(x))$ of $A(x)$:

$$A(x)Z^{-}(x) = O_{m \times (n-m)}. \qquad (2.1)$$

* $A^{-}(x)$ is an $n \times m$ right inverse of $A(x)$:

$$A(x)A^{-}(x) = I_{m}. \qquad (2.2)$$

In particular, the columns of $A^{-}(x)$ form a basis of a subspace complementary to $\mathcal{N}(A(x))$.

* $Z(x)$ is the unique $(n-m) \times n$ matrix such that

$$Z(x)Z^{-}(x) = I_{n-m} \quad \text{and} \quad Z(x)A^{-}(x) = 0_{(n-m) \times m}. \qquad (2.3)$$

From these properties, we can deduce the following identity:

$$A^{-}(x)A(x) + Z^{-}(x)Z(x) = I_{n}. \qquad (2.4)$$

For a motivation of this choice of notation and for practical examples of operators $A^{-}(x)$ and $Z^{-}(x)$, see Gabay [9].

From these operators can be introduced the notions of reduced gradient and Lagrange multiplier estimate. The reduced gradient of $f$ at $x$ is defined by

$$g(x) = Z^{-}(x)^{\top}\nabla f(x). \qquad (2.5)$$

Using (2.1), we have $g(x) = Z^{-}(x)^{\top}\nabla_{x}\ell(x, \lambda_{*})$, so that

$$\nabla g_{*}^{\top} = Z_{*}^{-\top}L_{*}. \qquad (2.6)$$

The first equation in (1.3) and (2.2) imply that $\lambda_{*} = -A_{*}^{-\top}\nabla f_{*}$. Therefore, we can take as Lagrange multiplier estimate, the vector

$$\lambda(x) = -A^{-}(x)^{\top}\nabla f(x). \qquad (2.7)$$
As previously, using (2.2) we have $\lambda(x) = -A^{-}(x)^{\top}\nabla_{x}\ell(x, \lambda_{*}) + \lambda_{*}$, so that

$$ \nabla\lambda_{*}^{\top} = -A_{*}^{-\top}L_{*}. \qquad (2.8) $$

The second useful decomposition of $\mathbb{R}^n$ differs from the first one by the choice of the subspace complementary to $\mathcal{N}(A(x))$. It comes from the form of the solution of the quadratic subproblem (1.5) and therefore it depends only on the problem data. Let $M_k$ be the current approximate Hessian with the property that $Z^{-}(x)^{\top}M_kZ^{-}(x)$ is positive definite, and define

$$ H_k(x) = (Z^{-}(x)^{\top}M_kZ^{-}(x))^{-1}. $$

Let $x$ be a point in $\Omega$ and consider the quadratic subproblem in $d$:

$$ \begin{array}{ll} \min & \nabla f(x)^{\top} d + \frac{1}{2} d^{\top} M_k d \\ \text{s.t.} & c(x) + A(x) d = 0. \end{array} \qquad (2.9) $$

Let us denote by $(d_k^{\text{QP}}(x), \lambda_k^{\text{QP}}(x))$ the primal-dual solution of (2.9). Using the first decomposition of $\mathbb{R}^n$ at $x$, it is not difficult to see that the primal solution can be written (see also Gabay [10])

$$ d_k^{\text{QP}}(x) = -Z^{-}(x)H_k(x)g(x) - (I - Z^{-}(x)H_k(x)Z^{-}(x)^{\top}M_k)A^{-}(x)c(x). \qquad (2.10) $$

Using (2.1) and (2.2), we find that the factor of $c(x)$ above satisfies

$$ A(x) [(I - Z^{-}(x)H_k(x)Z^{-}(x)^{\top}M_k) A^{-}(x)] = I_m. $$

Hence, the product of matrices inside the square brackets forms a right inverse of $A(x)$, which is denoted by $\hat{A}_k^{-}(x)$. Defining

$$ \hat{Z}_k(x) = H_k(x)Z^{-}(x)^{\top}M_k, \qquad (2.11) $$

we have

$$ \hat{A}_k^{-}(x) = (I - Z^{-}(x)\hat{Z}_k(x)) A^{-}(x), \qquad (2.12) $$

and thus (2.10) can be rewritten

$$ d_k^{\text{QP}}(x) = -Z^{-}(x)H_k(x)g(x) - \hat{A}_k^{-}(x)c(x). \qquad (2.13) $$

We have built a triplet $(Z^{-}(x), \hat{A}_k^{-}(x), \hat{Z}_k(x))$ satisfying conditions (2.1), (2.2), and (2.3), hence defining suitably a second decomposition of $\mathbb{R}^n$. In particular, we have

$$ \hat{A}_k^{-}(x)A(x) + Z^{-}(x)\hat{Z}_k(x) = I_n. \qquad (2.14) $$

Note that despite $A^{-}(x)$ and $Z^{-}(x)$ are used in formula (2.12), the operator $\hat{A}_k^{-}(x)$ does not depend on the choice of right inverse and tangent basis. Indeed, $-\hat{A}_k^{-}(x)c(x)$ is also defined as the solution of the quadratic subproblem (2.9) in which $\nabla f(x)$ is set to zero (see also the proof of Lemma 4.3 below).
In order to simplify the notation, we denote by $\hat{Z}_k$ and $\hat{A}_k^-$ the matrices $\hat{Z}_k(x_k)$ and $\hat{A}_k^-(x_k)$. With this convention, the direction $d_k$, solution of the quadratic subproblem (1.5) can be written

$$d_k = -Z_k^- H_k g_k - \hat{A}_k^- c_k. \qquad (2.15)$$

The vector $\hat{Z}_k d_k = -H_k g_k$ is called the reduced tangent direction.

Using $\hat{Z}_k(x) \hat{A}_k^-(x) = 0$ and the nonsingularity of $H_k(x)$, we have from (2.11) the following useful identity

$$Z^{-}(x)^{\top} M_k \hat{A}_{k}^{-}(x) = 0_{(n-m) \times m}. \qquad (2.16)$$

In particular, if $L_*$ is used in place of $M_k$ in the previous equality and if $x = x_*$, we obtain

$$Z_*^{-\top} L_* \hat{A}_*^- = 0. \qquad (2.17)$$

With (2.6), this shows that the columns of $\hat{A}_*^-$ form a basis of the space tangent to the reduced gradient manifold $\{g=0\}$ at $x_*$. Therefore, from (2.15), we see that the SQP direction $d_k$ has a longitudinal component $-Z_k^- H_k g_k$, tangent to the manifold $\{c=c_k\}$, and a transversal component $-\hat{A}_k^- c_k$, which tends to be tangent to the manifold $\{g=g_k\}$ when the pair $(x_k, Z_k^{-\top} M_k)$ is close to $(x_*, Z_*^{-\top} L_*)$.

In this paper, the globalization of the SQP method follows the approach of Bonnans [2]. We take as merit function the nondifferentiable augmented Lagrangian

$$\Theta_{\mu, \sigma}(x) = f(x) + \mu^{\top} c(x) + \sigma \|c(x)\|_{p}, \qquad (2.18)$$

in which $\mu \in \mathbb{R}^m$, $\sigma$ is a positive number, and $\|\cdot\|_p$ is an arbitrary (primal) norm on $\mathbb{R}^m$. This norm may differ from the $l_2$-norm and it is not squared in $\Theta_{\mu, \sigma}$. We denote by $\|\cdot\|_D$ the dual norm associated to $\|\cdot\|_p$ with respect to the Euclidean scalar product:

$$\|u\|_{D} = \sup_{\|v\|_{P}=1} u^{\top} v.$$

The penalty function $\Theta_{\mu, \sigma}$ is convenient for globalizing the SQP method for at least two reasons. On the one hand, the penalization is exact, provided the exactness condition

$$\|\mu - \lambda_*\|_D < \sigma \qquad (2.19)$$

holds (see for example Han and Mangasarian [15], and Bonnans [2]). On the other hand, the Armijo inequality using this function accepts the unit stepsize asymptotically, under some natural conditions (this is analyzed in Section 5, see also [2]).

We recall that $(\psi \circ \phi)$ has directional derivatives at a point $x$, if $\psi$ is Lipschitz continuous in a neighborhood of $\phi(x)$ and has directional derivatives at $\phi(x)$, and if $\phi$ has directional derivatives at $x$. Furthermore, $(\psi \circ \phi)'(x; h) = \psi'(\phi(x); \phi'(x; h))$. In particular, due to its convexity, a norm has the properties of function $\psi$ above, and since $f$ and $c$ are supposed smooth, $\Theta_{\mu, \sigma}$ has directional derivatives.

We conclude this section by giving formulae for the directional derivatives of $\Theta_{\mu, \sigma}$ and by giving conditions for having descent directions. Let $d$ be a vector of $\mathbb{R}^n$
satisfying the linear constraints $c(x) + A(x)d = 0$. The directional derivative of $\Theta_{\mu,\sigma}$ at $x$ in the direction $d$ is given by

$$ \Theta'_{\mu,\sigma}(x; d) = \nabla f(x)^{\top} d - \mu^{\top} c(x) - \sigma \|c(x)\|_{p} \quad (2.20) $$

(for the differentiation of the term with the norm, use the very definition of directional derivative, see for example [12]). For any multiplier $\lambda$, we then have

$$ \Theta'_{\mu,\sigma}(x; d) = \nabla_x \ell(x, \lambda)^{\top} d + (\lambda - \mu)^{\top} c(x) - \sigma \|c(x)\|_{p}. \quad (2.21) $$

Therefore, if $d$ is a descent direction of the Lagrangian function at $(x, \lambda)$, in the sense that $\nabla_x \ell(x, \lambda)^{\top} d < 0$ (which in particular holds for the direction $d_k$ when $(x, \lambda) = (x_k, \lambda_k^{QP})$), then $d$ is also a descent direction of $\Theta_{\mu,\sigma}$ at $x$ provided the descent condition

$$ \|\lambda - \mu\|_{D} \leq \sigma \quad (2.22) $$

holds (compare with the exactness condition (2.19)).

# 3 The approach

This section describes our quasi-Newton version of the SQP algorithm in a global framework. The aim we pursue is to develop a consistent way of updating the positive definite matrix $M_k$, using convenient vectors $\gamma_k$ and $\delta_k$.

## 3.1 Computation of $\gamma_k$

As we said in the introduction, we take the point of view to force $M_k$ to be an approximation of the Hessian of the augmented Lagrangian. This is equivalent to considering the problem

$$ \begin{array}{ll} \min & f(x) + \frac{r}{2} \|c(x)\|^2 \\ \text{s.t.} & c(x) = 0, \quad x \in \Omega, \end{array} $$

for some $r \ge 0$. This problem has the same solutions as problem (1.1) and has a Lagrangian whose Hessian at $(x_*, \lambda_*)$ is

$$ L_{*}^{r} = L_{*} + r A_{*}^{\top} A_{*}. $$

It is well known that when (1.4) holds, $L_*^r$ is positive definite when $r$ is larger than some threshold. Therefore, it makes sense to force $M_k$ to approach $L_*^r$ for some sufficiently large $r$ and to keep its positive definiteness.

For this purpose, we would like to have for some $r_k \ge 0$:

$$ \gamma_k \approx L_*^{r_k} \delta_k \approx L_* \delta_k + r_k A_k^\top A_k \delta_k. $$

Using successively (2.4), (2.6), and (2.8), we get

$$ \begin{aligned} L_* \delta_k &= Z_k^\top Z_k^{-\top} L_* \delta_k + A_k^\top A_k^{-\top} L_* \delta_k \\ &\approx Z_k^\top \nabla g_*^\top \delta_k - A_k^\top \nabla \lambda_*^\top \delta_k \\ &\approx Z_k^\top (g_{k+1} - g_k) - A_k^\top (\lambda_{k+1} - \lambda_k), \end{aligned} \quad (3.1) $$
provided

$$ \delta_k \approx x_{k+1} - x_k. $$

This approximate computation motivates our choice of $\gamma_k$, which is

$$ \gamma_k = Z_k^\top (g_{k+1} - g_k) - A_k^\top (\lambda_{k+1} - \lambda_k) + r_k A_k^\top A_k \delta_k. \quad (3.2) $$

This choice is very close to the value given by formula (1.11), which is used by Byrd, Tapia, and Zhang (BTZ for short). The main difference is that $\gamma_k^\ell$ is split in two terms for reasons that are discussed now. For this, let us look at the form of the scalar product $\gamma_k^\top \delta_k$, which we want to have positive:

$$ \gamma_k^\top \delta_k = (g_{k+1} - g_k)^\top Z_k \delta_k - (\lambda_{k+1} - \lambda_k)^\top A_k \delta_k + r_k \|A_k \delta_k\|^2. \quad (3.3) $$

When $A_k \delta_k \neq 0$, it is clear that the curvature condition (1.8) can be satisfied by choosing $r_k$ sufficiently large. Remember that when $A_k \delta_k$ is close to zero, the BTZ approach needs a back-up strategy. For our form of $\gamma_k$, $A_k \delta_k = 0$ implies that

$$ \gamma_k^\top \delta_k = (g_{k+1} - g_k)^\top Z_k \delta_k. $$

A possible way of satisfying the curvature condition in this case would be to choose the next iterate $x_{k+1}$ such that $g_{k+1}^\top Z_k \delta_k > g_k^\top Z_k \delta_k$. We believe, however, that this may not be possible at iteration where $A_k \delta_k \neq 0$, because $Z_k \delta_k$ may not be a reduced descent direction (meaning that $g_k^\top Z_k \delta_k$ may not be negative). Now when $A_k \delta_k = 0$ and $\delta_k$ is parallel to the SQP direction $d_k$, we have $c_k = 0$ and, from (2.15), $d_k$ reduces to $d_k = -Z_k^- H_k g_k = Z_k^- \hat{Z}_k d_k$, which implies that $Z_k d_k = \hat{Z}_k d_k$. Therefore, by forcing the inequality

$$ g_{k+1}^\top \hat{Z}_k d_k > g_k^\top \hat{Z}_k d_k, $$

the curvature condition can be fulfilled when $A_k \delta_k = 0$. The important point is that, as we shall see, it is always possible to realize this inequality, even when $c_k \neq 0$, because $\hat{Z}_k d_k = -H_k g_k$ is a reduced descent direction ($g_k^\top \hat{Z}_k d_k < 0$). The piecewise line-search (PLS) technique introduced for reduced quasi-Newton methods in [11] and extended in [12] is designed for realizing this inequality.

One can view our update scheme of $M_k$ as follows. Multiplying both sides of $L_*^{r_k}$, the matrix $M_k$ has to approach, by the left hand side of (2.14) evaluated at $x = x_*$ and using the decoupling identity (2.17), we obtain

$$ L_*^{r_k} = \hat{Z}_*^\top Z_*^{-\top} L_* Z_*^{-\top} \hat{Z}_* + A_*^\top (\hat{A}_*^{-\top} L_* \hat{A}_*^{-\top} + r_k I) A_* $$

We see that, since the cross-term $\hat{A}_*^{-\top} L_*^{r_k} Z_*^{-\top}$ vanishes, the positive definiteness of $L_*^{r_k}$ can be obtained from that of $Z_*^{-\top} L_* Z_*^{-\top}$ (recall assumption (1.4)) and by forcing that of $(\hat{A}_*^{-\top} L_* \hat{A}_*^{-\top} + r_k I)$ with a sufficiently large $r_k$. This decomposition in longitudinal and transversal components is met again in the tools used to form $\gamma_k$. The PLS takes care of the case when $\delta_k$ is tangent to the constraint manifold by finding a suitable point $x_{k+1}$, while the setting of $r_k$ is helpful for a displacement $\delta_k$ in the range space of $A_k^-$.
## 3.2 Guiding path

From the discussion above, a central point of our algorithm is to find the next iterate $x_{k+1}$ in order to get, in particular, the following reduced Wolfe condition

$$g_{k+1}^{\top} \hat{Z}_k d_k \geq \omega_2 g_k^{\top} \hat{Z}_k d_k, \qquad (3.4)$$

for some constant $\omega_2 \in (0, 1)$.

Contrary to the unconstrained case, condition (3.4) may fail, whatever point $x_{k+1}$ is taken along the SQP direction $d_k$. On the other hand, Proposition 3.2 below shows that along the path $p_k$ defined by the following differential equation

$$\left\{ \begin{array}{l} p'_{k}(\xi) = Z^{-}(p_{k}(\xi))\hat{Z}_{k}d_{k} - \hat{A}_{k}^{-}(p_{k}(\xi))c(p_{k}(\xi)) \\ p_{k}(0) = x_{k}, \end{array} \right. \qquad (3.5)$$

one can find a stepsize $\xi_k$, such that the merit function $\Theta_{\mu,\sigma}$ decreases and the reduced Wolfe condition (3.4) holds:

$$\Theta_{\mu,\sigma}(p_k(\xi_k)) \leq \Theta_{\mu,\sigma}(x_k) \quad \text{and} \quad g(p_k(\xi_k))^{\top} \hat{Z}_k d_k \geq \omega_2 g_k^{\top} \hat{Z}_k d_k. \qquad (3.6)$$

Note that the reduced tangent component of $p'_{k}(\xi)$ keeps the constant value $\hat{Z}_k d_k$ along the path. This is further motivated in [12].

In the proof of Proposition 3.2, we will need the following lemma. We say that a function $\phi$ is locally Lipschitz continuous on a set $X$ if any point of $X$ has a neighborhood on which $\phi$ is Lipschitz continuous.

**Lemma 3.1.** Let $\alpha > 0$ and $\phi : [0, \alpha] \to \Omega$ be a continuous function having right derivatives on $(0, \alpha)$. Suppose that $f$ and $c$ are locally Lipschitz continuous on $\Omega$ and have directional derivatives on $\phi((0, \alpha))$. Then there exists $\bar{\alpha} \in (0, \alpha)$ such that

$$\Theta_{\mu,\sigma}(\phi(\alpha)) - \Theta_{\mu,\sigma}(\phi(0)) \leq \alpha \Theta'_{\mu,\sigma}(\phi(\bar{\alpha}); \phi'(\bar{\alpha}; 1)).$$

**Proof.** Since $c$ is locally Lipschitz continuous on $\Omega$, so is $\|c(\cdot)\|_P$. Furthermore, by the hypotheses on $c$ and the convexity of the norm, $\|c(\cdot)\|_P$ has directional derivatives on $\phi((0, \alpha))$. Therefore, with the hypotheses, we deduce that $\Theta_{\mu,\sigma}$ is locally Lipschitz continuous on $\Omega$ and has directional derivatives on $\phi((0, \alpha))$. Now with the properties of $\phi$, we see that $\Theta_{\mu,\sigma} \circ \phi$ has right derivatives on $(0, \alpha)$ and that for any $\bar{\alpha} \in (0, \alpha)$:

$$(\Theta_{\mu,\sigma} \circ \phi)'(\bar{\alpha}; 1) = \Theta'_{\mu,\sigma}(\phi(\bar{\alpha}); \phi'(\bar{\alpha}; 1)).$$

On the other hand, the function $\Theta_{\mu,\sigma} \circ \phi$ is continuous on $[0, \alpha]$ and, since it has right derivatives on $(0, \alpha)$, there exists $\bar{\alpha} \in (0, \alpha)$ such that

$$(\Theta_{\mu,\sigma} \circ \phi)(\alpha) - (\Theta_{\mu,\sigma} \circ \phi)(0) \leq \alpha (\Theta_{\mu,\sigma} \circ \phi)'(\bar{\alpha}; 1).$$

(see for instance Schwartz [26, Chap. III, § 2, Remarque 11]).

Combining this inequality with the preceding equality gives the result. $\square$
**Proposition 3.2.** Suppose that the path $\xi \mapsto p_k(\xi)$ defined by (3.5) exists for sufficiently large stepsize $\xi \ge 0$. Suppose also that $f$ and $c$ are continuously differentiable, that $\Theta_{\mu,\sigma}$ is bounded from below along the path $p_k$, that $\|\lambda_k^{\text{QP}}(p_k(\xi)) - \mu\|_D \le \sigma$ whenever $p_k(\xi)$ exists, that $M_k$ is positive definite, and that $\omega_2 \in (0, 1)$. Then, the inequalities in (3.6) are satisfied for some stepsize $\xi_k > 0$.

**Proof.** To lighten the notation in the proof, we denote by $(d(\xi), \lambda(\xi)) = (d_k^{\text{QP}}(p_k(\xi)), \lambda_k^{\text{QP}}(p_k(\xi)))$ the primal-dual solution of the quadratic subproblem

$$ \begin{array}{ll} \min & \nabla f(p_k(\xi))^{\top} d + \frac{1}{2} d^{\top} M_k d \\ \text{s.t.} & c(p_k(\xi)) + A(p_k(\xi))d = 0. \end{array} $$

The first order optimality conditions give

$$ \nabla_x \ell(p_k(\xi), \lambda(\xi)) = -M_k d(\xi), $$

and $d(\xi)$ can be written (see (2.13))

$$ d(\xi) = -Z^{-}(p_{k}(\xi))H_{k}(p_{k}(\xi))g(p_{k}(\xi)) - \hat{A}_{k}^{-}(p_{k}(\xi))c(p_{k}(\xi)). $$

Let us show that, when the second inequality in (3.6) or reduced Wolfe condition does not hold for $\xi_k = \xi$, then

$$ \Theta'_{\mu,\sigma}(p_k(\xi); p'_k(\xi)) < -\omega_2 g_k^{\top} H_k g_k. \quad (3.7) $$

Using successively (2.21), the hypothesis $\|\lambda(\xi) - \mu\|_D \le \sigma$, the optimality condition above, the form of $d(\xi)$ and $p'_k(\xi)$, the identity (2.16), and the positive definiteness of $M_k$, we get

$$ \begin{align*} \Theta'_{\mu,\sigma}(p_k(\xi); p'_k(\xi)) & = \nabla_x \ell(p_k(\xi), \lambda(\xi))^{\top} p'_k(\xi) + (\lambda(\xi) - \mu)^{\top} c(p_k(\xi)) - \sigma \|c(p_k(\xi))\|_P \\ & \le \nabla_x \ell(p_k(\xi), \lambda(\xi))^{\top} p'_k(\xi) \\ & = -d(\xi)^{\top} M_k p'_k(\xi) \\ & = -g(p_k(\xi))^{\top} H_k g_k - c(p_k(\xi))^{\top} \hat{A}_k^{-} (p_k(\xi))^{\top} M_k \hat{A}_k^{-} (p_k(\xi)) c(p_k(\xi)) \\ & \le -g(p_k(\xi))^{\top} H_k g_k. \end{align*} $$

Therefore, when the reduced Wolfe condition does not hold, we have (3.7).

On the other hand, we see by using Lemma 3.1 with $\phi = p_k$ that, as long as the path $p_k$ exists, for any $\xi > 0$, one can find $\bar{\xi} \in (0, \xi)$ such that

$$ \Theta_{\mu,\sigma}(p_k(\xi)) - \Theta_{\mu,\sigma}(x_k) \le \xi \Theta'_{\mu,\sigma}(p_k(\bar{\xi}); p'_k(\bar{\xi})). $$

Therefore, if the reduced Wolfe condition is never realized along the path $p_k$, we would have by (3.7)

$$ \Theta_{\mu,\sigma}(p_k(\xi)) - \Theta_{\mu,\sigma}(x_k) < -\xi \omega_2 g_k^{\top} H_k g_k, \quad (3.8) $$
which would imply the unboundedness of the merit function along this path and would contradict the hypotheses.

At the first stepsize $\xi_k > 0$ at which the reduced Wolfe condition is satisfied, inequality (3.8) shows that the merit function $\Theta_{\mu,\sigma}$ has decreased. This concludes the proof. $\square$

The inequality $||\lambda_k^{\text{QP}}(p_k(\xi)) - \mu||_D \le \sigma$ used as hypothesis in the previous proposition can be compared with the descent condition (2.22).

## 3.3 Outline of the PLS algorithm

The success of the path $p_k$ defined by (3.5) suggests to search for the next iterate $x_{k+1}$ along a discretized version of this path. Taking a precise discretization may not succeed and would be computationally expensive. Also, we propose to take as often as possible a unit stepsize along the directions obtained by an explicit Euler discretization of the differential equation (3.5). With this technique, the search path becomes piecewise linear. It is proved in Section 4 that the search along this path succeeds in a finite number of trials.

The piecewise line-search (PLS) algorithm generates intermediate points $x_{k,i}$, for $i = 0, \dots, i_k$, with $x_{k,0} = x_k$ and $x_{k,i_k} = x_{k+1}$. We adopt the notation $f_{k,i} = f(x_{k,i})$, $\nabla f_{k,i} = \nabla f(x_{k,i})$, $c_{k,i} = c(x_{k,i})$, $Z_{k,i} = Z^{-}(x_{k,i})$, $A_{k,i}^{-} = A^{-}(x_{k,i})$, $\hat{A}_{k,i}^{-} = \hat{A}_{k}^{-}(x_{k,i})$, and $\hat{Z}_{k,i} = \hat{Z}_{k}(x_{k,i})$. The iterations of the PLS algorithm, computing $x_{k,i+1}$ from $x_{k,i}$, are called inner iterations and their number is denoted by $i_k$.

The point $x_{k,i+1}$ is obtained from $x_{k,i}$ by

$$x_{k,i+1} = x_{k,i} + \alpha_{k,i} d_{k,i}, \qquad (3.9)$$

where the stepsize $\alpha_{k,i} > 0$ is determined along the direction

$$d_{k,i} = -Z_{k,i}^{-} H_{k} g_{k} - \hat{A}_{k,i}^{-} c_{k,i}. \qquad (3.10)$$

This direction is obtained by evaluating the right hand side of (3.5) at a discretization point $x_{k,i}$ of the path $p_k$. The stepsize is chosen such that the following two conditions are satisfied for $\alpha = \alpha_{k,i}$:

$$x_{k,i} + \alpha d_{k,i} \in \Omega, \qquad (3.11)$$

$$\Theta_{\mu_k, \sigma_{k,i}}(x_{k,i} + \alpha d_{k,i}) \le \Theta_{\mu_k, \sigma_{k,i}}(x_{k,i}) + \omega_1 \alpha \Theta'_{\mu_k, \sigma_{k,i}}(x_{k,i}; d_{k,i}). \qquad (3.12)$$

Condition (3.12) imposes a sufficient decrease of the merit function and will be called the Armijo condition.

At each inner iteration $i$, the penalty parameter $\sigma_{k,i}$ may need to be adapted so that $d_{k,i}$ is a descent direction of $\Theta_{\mu_k, \sigma_{k,i}}$ at $x_{k,i}$. An adaptation rule will be given in Section 4.

Next the reduced Wolfe condition

$$g(x_{k,i+1})^{\top} \hat{Z}_k d_k \ge \omega_2 g_k^{\top} \hat{Z}_k d_k \qquad (3.13)$$
is tested. If it holds, the PLS is completed and $i_k$ is set to $i+1$. Otherwise, the index $i$ is increased by one and the search is pursued along a new direction $d_{k,i}$.

From the description of the algorithm, we have

$$x_{k+1} = x_{k,i_k} = x_k + \sum_{i=0}^{i_k-1} \alpha_{k,i} d_{k,i}.$$

It is interesting to compare the PLS algorithm with a skipping rule strategy (no matrix update when $\gamma_k^\top \delta_k$ is not sufficiently positive). Indeed, the intermediate search directions $d_{k,i}$ are close to the SQP direction at $x_{k,i}$, with two differences however. First, the matrix $M_k$ is kept unchanged as long as the reduced Wolfe condition is not satisfied, which is similar to the skipping rule strategy. On the other hand, the reduced gradient used in these directions is also kept unchanged. As we will see (Theorem 4.4), this gives a chance to the matrix to be updated.

## 3.4 Computation of $\delta_k$

The choice of $\delta_k$ is governed by the necessity to have $\delta_k \simeq x_{k+1} - x_k$, as required by the discussion in Section 3.1, and the desire to control closely the positivity of $\gamma_k^\top \delta_k$ when $A_k \delta_k = 0$. We have already observed that when $A_k \delta_k = 0$,

$$\gamma_k^\top \delta_k = (g_{k+1} - g_k)^\top Z_k \delta_k,$$

so that $r_k$ cannot be used to get $\gamma_k^\top \delta_k > 0$.

Suppose that we choose $\delta_k = x_{k+1} - x_k$. Then, from the identity above, we have $\gamma_k^\top \delta_k = (g_{k+1} - g_k)^\top Z_k (x_{k+1} - x_k)$, and it is not clear how the reduced Wolfe condition (3.4) can assure the positivity of $\gamma_k^\top \delta_k$, since $x_{k+1} - x_k$ is not parallel to $d_k$. For this reason, we prefer to take for $\delta_k$ the following approximation of $x_{k+1} - x_k$:

$$\begin{aligned}
\delta_k &= -\alpha_k^z Z_k^- H_k g_k - \alpha_k^\Lambda \hat{A}_k^- c_k \quad (3.14) \\
&\approx \sum_{j=0}^{i_k-1} \alpha_{k,i} (-Z_{k,i}^- H_k g_k - \hat{A}_{k,i}^- c_{k,i}) \\
&= x_{k+1} - x_k.
\end{aligned}$$

In (3.14), the longitudinal stepsize $\alpha_k^z$ and the transversal stepsize $\alpha_k^\Lambda$ are defined by

$$\alpha_k^z = \sum_{i=0}^{i_k-1} \alpha_{k,i} \quad \text{and} \quad \alpha_k^\Lambda = \sum_{i=0}^{i_k-1} \alpha_{k,i} e^{-\xi_{k,i}}, \quad (3.15)$$

with $\xi_{k,i} = \sum_{j=0}^{i_k-1} \alpha_{k,j}$. The form of $\alpha_k^\Lambda$ aims at taking into account the fact that the value of $c$ at $x_{k,i}$ is used in the search directions $d_{k,i}$, while only $c_k$ is used in $\delta_k$. It is based on the observation that along the path $p_k$ defined by (3.5), we have $c(p_k(\xi)) = e^{-\xi} c_k$ (multiply both sides of (3.5) by $A_k(p_k(\xi))$ and integrate). After discretization: $c_{k,i} \simeq e^{-\xi_{k,i}} c_k$.
To check that our choice (3.14) for $\delta_k$ is appropriate, suppose that $A_k\delta_k = 0$. Then, we have $c_k = 0$, hence $\delta_k = \alpha_k^z Z_k^- \hat{Z}_k d_k$, and this allows us to write

$$ \begin{aligned} \gamma_k^\top \delta_k &= (g_{k+1} - g_k)^\top Z_k \delta_k \\ &= \alpha_k^z (g_{k+1} - g_k)^\top \hat{Z}_k d_k \\ &> 0, \end{aligned} $$

by the reduced Wolfe condition (3.4). By a continuity argument, one can claim that $(g_{k+1} - g_k)^\top Z_k \delta_k$ is also positive when $x_k$ is close to the constraint manifold, provided the stepsizes are determined by processes depending continuously on $x_k$ and (3.4) is realized with strict inequality (in this case, the number of inner iterations in the PLS algorithm does not change in the neighborhood of a point on the constraint manifold). In the algorithm below, we shall not impose strict inequality in (3.4), because we believe that this continuity argument is not important in practice.

The conclusion of this discussion is that for any $k \ge 1$, one can find a (finite) $r_k \ge 0$ such that $\gamma_k^\top \delta_k > 0$, either because $A_k \delta_k \ne 0$ or because $A_k \delta_k = 0$ and $\gamma_k^\top \delta_k > 0$ by the reduced Wolfe condition (3.4).

## 3.5 Computation of $r_k$

Formula (2.10) shows that only the part $Z^-(x)^\top M_k$ of $M_k$ plays a role in the determination of the SQP direction and that this direction is well defined provided $Z^-(x)^\top M_k Z^-(x)$ is nonsingular. In this case, because of (2.1), adding a positive multiple of $A_k^\top A_k$ to $M_k$ does not modify the SQP direction. In our case, $r_k$ is aimed at forcing $M_{k+1}$ to approach the Hessian of the augmented Lagrangian $L_* + r_k A_*^\top A_*$. But since $r_k$ intervenes nonlinearly in the matrix $M_{k+1}$ via the vector $\gamma_k$ and the BFGS formula (1.7), its value affects $d_{k+1}$. This discussion suggests, however, that the value of $r_k$ could be set from considerations based only on the matrix update.

In the algorithm below, we choose $r_k$ in order to minimize a measure of the conditioning of the matrix $M_{k+1}$. This measure is the function

$$ \omega(M) = \frac{\operatorname{tr} M}{\det M^{1/n}} $$

introduced by Dennis and Wolkowicz [7]. Interestingly enough, minimizing $\omega(M)$ on a subset $S$ of the set of positive definite matrices is equivalent to minimizing in $(\zeta, M) \in (0, +\infty) \times S$ the function $\psi(\zeta M)$, where $\psi(M) = \operatorname{tr}(M) - \ln \det(M)$ is the conditioning measure introduced by Byrd and Nocedal [3]. In both cases, one tries to find the matrix $M \in S$, in a certain sense the closest to the set $\{\zeta I : \zeta > 0\}$ of positive definite matrices with unit condition number.

The next proposition analyzes the problem of minimizing $\omega(M_{k+1})$ with respect to $r_k$.

**Proposition 3.3.** Let $\eta$, $\delta$, and $\tilde{\gamma}$ be vectors in $\mathbb{R}^n$ such that $\eta^\top \delta > 0$ and let $\gamma(r) = \tilde{\gamma} + r\eta$, where $r$ is a scalar parameter belonging to the interval $\mathcal{R} = \{r : \gamma(r)^\top \delta > 0\} =$
($-\tilde{\gamma}^\top \delta / \eta^\top \delta, +\infty$). Let $M$ be a positive definite matrix and let $M(r)$ be the matrix obtained by the BFGS formula, using $\gamma(r)$ and $\delta$:

$$M(r) = M - \frac{M \delta \delta^\top M}{\delta^\top M \delta} + \frac{\gamma(r) \gamma(r)^\top}{\gamma(r)^\top \delta}.$$

Then, the function $r \mapsto \omega(M(r))$ is uniquely minimized on $\mathbb{R}$ by

$$\bar{r} = \frac{\bar{t} - \tilde{\gamma}^{\top} \delta}{\eta^{\top} \delta}, \qquad (3.16)$$

where

$$\bar{t} = a + (a^2 + b)^{\frac{1}{2}}, \quad a = \frac{c_1}{2(n-1)c_2}, \quad b = \frac{n+1}{n-1} \frac{c_0}{c_2},$$

$$c_0 = \|\gamma(-\frac{\tilde{\gamma}^{\top}\delta}{\eta^{\top}\delta})\|^2, \quad c_1 = \text{tr} M - \frac{\|M\delta\|^2}{\delta^{\top} M \delta} + 2(\frac{\tilde{\gamma}^{\top}\eta}{\eta^{\top}\delta} - \|\eta\|^2 \frac{\tilde{\gamma}^{\top}\delta}{(\eta^{\top}\delta)^2}), \quad c_2 = \frac{\|\eta\|^2}{(\eta^{\top}\delta)^2}.$$

**Proof.** Let us show that $r \mapsto \omega(M(r))$ is pseudoconvex on $\mathbb{R}$, which in particular implies that any stationary point is a global minimizer (see [18]). By using

$$\text{tr} M(r) = \text{tr} M - \frac{\|M\delta\|^2}{\delta^{\top} M \delta} + \frac{\|\gamma(r)\|^2}{\gamma(r)^{\top} \delta}$$

and (see [19])

$$\det M(r) = \frac{\gamma(r)^{\top} \delta}{\delta^{\top} M \delta} \det M,$$

we have $\text{tr} M(r) = \varphi(\gamma(r)^{\top} \delta)$, where $\varphi(t) = \frac{c_0}{t} + c_1 + c_2t$ with the $c_i$ given in the statement of the proposition. Since $c_0$ is nonnegative, $\varphi$ is convex on $(0, +\infty)$ and, because $r \mapsto \gamma(r)^{\top} \delta$ is affine, $\text{tr} M(\cdot)$ is convex on $\mathbb{R}$. On the other hand, $M(r)$ is positive definite for $r \in \mathbb{R}$ and $r \mapsto \det M(r)$ is affine. Hence the function $\omega(M(\cdot))$ is the quotient of a positive convex function and a positive concave function and thus is pseudoconvex on $\mathbb{R}$ (see [18, p. 148]).

Now by using $\omega(M(r)) = (\frac{\delta^{\top} M \delta}{\det M})^{\frac{1}{n}} \tilde{\omega}(\gamma(r)^{\top} \delta)$, with $\tilde{\omega}(t) = t^{-\frac{1}{n}} \varphi(t)$, a straightforward calculation shows that $\bar{t}$ is the unique stationary point of $\tilde{\omega}(\cdot)$ on $(0, +\infty)$, hence $\bar{r}$ is the unique global minimum of $\omega(M(\cdot))$ on $\mathbb{R}$. $\square$

As we said above, setting $\eta = A_k^\top A_k \delta_k$, $\delta = \delta_k$, and $\tilde{\gamma} = \tilde{\gamma}_k = Z_k^\top (g_{k+1} - g_k) - A_k^\top (\lambda_{k+1} - \lambda_k)$ in Proposition 3.3, we take in our algorithm $r_k = \bar{r}$ given by formula (3.16). This value is not defined when $\eta_k^\top \delta_k = \|A_k \delta_k\|^2 = 0$. In this case, the discussion in Section 3.4 has shown that setting $r_k = 0$ is appropriate because the PLS alone ensures the positivity of the scalar product $\gamma_k^\top \delta_k$. From a numerical point of view, the difficulty is now to decide when $\|A_k \delta_k\|$ is sufficiently small so that $r_k$ can be set to zero in formula (3.2) while preserving the positivity of $\gamma_k^\top \delta_k$.

To appreciate the smallness of $\|A_k \delta_k\|$, we compare $\tilde{\gamma}_k^\top \delta_k^z$ to $\tilde{\gamma}_k^\top \delta_k^\wedge$, where $\delta_k^z$ and $\delta_k^\wedge$ are the longitudinal and transversal components of $\delta_k$, respectively:

$$\delta_k^z = -\alpha_k^z Z_k^- H_k g_k \quad \text{and} \quad \delta_k^\wedge = -\alpha_k^\wedge \hat{A}_k^- c_k.$$
Suppose that the reduced Wolfe condition (3.4) holds. Then $\tilde{\gamma}_k^\mathrm{T}\delta_k^z = \alpha_k^z(g_{k+1} - g_k)^\mathrm{T}$
$\hat{Z}_k d_k$ is positive and when

$$
\tilde{\gamma}_k^\mathrm{T} \delta_k^\mathrm{A} \geq -\beta \tilde{\gamma}_k^\mathrm{T} \delta_k^z, \quad (3.17)
$$

for some constant $\beta \in (0, 1)$, we have $\tilde{\gamma}_k^\mathrm{T}\delta_k = \tilde{\gamma}_k^\mathrm{T}\delta_k^z + \tilde{\gamma}_k^\mathrm{T}\delta_k^\mathrm{A} \ge (1-\beta)\tilde{\gamma}_k^\mathrm{T}\delta_k^z > 0$, so that $r_k$ can be set to zero. Therefore, we adopt the following rule for the update of $r_k$.

COMPUTATION OF $r_k$:

if (3.4) and (3.17) hold then $r_k = 0$,

else $r_k = \max(0, \bar{r}_k)$, where $\bar{r}_k$ is given by (3.16), in which $\eta = A_k^\mathrm{T} A_k \delta_k$,

$\delta = \delta_k$, and $\tilde{\gamma} = Z_k^\mathrm{T}(g_{k+1} - g_k) - A_k^\mathrm{T}(\lambda_{k+1} - \lambda_k)$.

In our numerical experiment, we set $\beta = 0.1$ in (3.17).

Since $\tilde{\gamma}_k$ aims at approximating $L_*\delta_k$ (see (3.1)), (3.17) may be seen as a way of comparing the curvature of the Lagrangian along $\delta_k^\mathrm{A}$ and $\delta_k^z$. The rule above for computing $r_k$ can be read as follows: if the curvature of the Lagrangian in the transversal direction is positive or not too negative (with respect to the longitudinal curvature), set $r_k$ to zero, otherwise use formula (3.16). One can also say that (3.17) is a way of measuring the smallness of $\|A_k\delta_k\|$, since when this quantity vanishes, $\delta_k^\mathrm{A} = 0$ and (3.17) readily holds.

4 The piecewise line-search

In this section, we make more precise the PLS algorithm outlined in Section 3.3, show
its well-posedness (Proposition 4.1), and prove its finite termination (Theorem 4.4).

4.1 Descent directions

A question we have not addressed so far is to know whether the *i*th inner search
direction $d_{k,i}$, given by (3.10) and used in the PLS algorithm, is a descent direction
of the merit function. The following proposition shows that this property holds when
the penalty parameter $\sigma_{k,i}$ in the merit function is larger than a threshold easy to
compute. For this, as in Proposition 3.2, the multiplier $\mu = \mu_k$ in the merit function is
compared to the multiplier $\lambda_k^{\mathrm{QP}}(x_{k,i})$ given by the quadratic program (2.9) at $x = x_{k,i}$.
The multiplier $\mu_k$ is indexed by $k$ because it will have to be modified at some iterations
of the overall algorithm below.

**Proposition 4.1.** Let $0 \le i < i_k$ be the index of an inner iteration of the PLS algorithm. Suppose that $x_k$ is not a stationary point, that $M_k$ is positive definite, and that

$$
\underline{\sigma}_k + \|\lambda_k^{\mathrm{QP}}(x_{k,i}) - \mu_k\|_D \le \sigma_{k,i}, \quad (4.1)
$$

where $\underline{\sigma}_k$ is a positive number. Then $d_{k,i}$ is a descent direction of $\Theta_{\mu_k, \sigma_{k,i}}$ at $x_{k,i}$, meaning that $\Theta'_{\mu_k, \sigma_{k,i}}(x_{k,i}; d_{k,i}) < 0$. For $i = 0$:

$$
\Theta'_{\mu_k, \sigma_{k,0}}(x_{k,0}; d_{k,0}) \le -d_k^\mathrm{T} M_k d_k - \underline{\sigma}_k \|c_k\|_P, \quad (4.2)
$$
while for $1 \le i < i_k$:

$$ \Theta'_{\mu_k, \sigma_{k,i}}(x_{k,i}; d_{k,i}) < -\omega_2 g_k^\top H_k g_k - c_{k,i}^\top \hat{A}_{k,i}^{-\top} M_k \hat{A}_{k,i}^{-} c_{k,i} - \underline{\sigma}_k \|c_{k,i}\|_p. \quad (4.3) $$

**Proof.** For $i = 0$, the search direction is $d_{k,0} = d_k$ and the optimality conditions of (2.9) give

$$ \nabla_x l(x_k, \lambda_k^{\text{QP}}(x_k))^{\top} d_{k,0} = -d_k^{\top} M_k d_k. $$

For $i = 1, \dots, i_k - 1$, we have by the optimality conditions of (2.9), formulae (2.13) and (3.10), identity (2.16), and the fact that the reduced Wolfe condition (3.13) is not satisfied at $x_{k,i}$:

$$ \begin{aligned} \nabla_x l(x_{k,i}, \lambda_k^{\text{QP}}(x_{k,i}))^{\top} d_{k,i} &= -d_k^{\text{QP}}(x_{k,i})^{\top} M_k d_{k,i} \\ &= -(-Z_{k,i}^{-} H_{k,i} g_{k,i} - \hat{A}_{k,i}^{-} c_{k,i})^{\top} M_k (-Z_{k,i}^{-} H_k g_k - \hat{A}_{k,i}^{-} c_{k,i}) \\ &= -g_{k,i}^{\top} H_k g_k - c_{k,i}^{\top} \hat{A}_{k,i}^{-\top} M_k \hat{A}_{k,i}^{-} c_{k,i} \\ &< -\omega_2 g_k^{\top} H_k g_k - c_{k,i}^{\top} \hat{A}_{k,i}^{-\top} M_k \hat{A}_{k,i}^{-} c_{k,i}. \end{aligned} $$

Then, from the estimates above, (2.21), and (4.1), we see that (4.2) and (4.3) hold. Since $x_k$ is not stationary, $d_{k,0} = d_k \ne 0$ and (4.2) shows that $d_{k,0}$ is a descent direction of $\Theta_{\mu_k, \sigma_{k,0}}$ at $x_k$. The strict inequality in (4.3) shows that for $i \ge 1$, $d_{k,i}$ is also a descent direction of $\Theta_{\mu_k, \sigma_{k,i}}$ at $x_{k,i}$. $\square$

The preceding result suggests the following rule for updating the penalty parameter $\sigma_{k,i}$. Let us denote by $\sigma_{k,0}$ the value of the penalty parameter at the beginning of the PLS. This value depends on the update of $\mu_k$ in the overall algorithm, which will be given in Section 5.

**UPDATE RULE OF $\sigma_{k,i}$:**

$$ \begin{array}{ll} \text{if } \underline{\sigma}_k + \|\lambda_k^{\text{QP}}(x_{k,i}) - \mu_k\|_D \le \sigma_{k,i-1} & \text{then } \sigma_{k,i} = \sigma_{k,i-1}, \\ \text{else } \sigma_{k,i} = \max(2\sigma_{k,i-1}, \underline{\sigma}_k + \|\lambda_k^{\text{QP}}(x_{k,i}) - \mu_k\|_D). & \end{array} $$

It follows that either

$$ \sigma_{k,i} = \sigma_{k,i-1}, \quad (4.4) $$

or

$$ \underline{\sigma}_k + \|\lambda_k^{\text{QP}}(x_{k,i}) - \mu_k\|_D > \sigma_{k,i-1} \quad \text{and} \quad \sigma_{k,i} \ge 2\sigma_{k,i-1}. \quad (4.5) $$

With this update rule, the search direction $d_{k,i}$ is a descent direction of $\Theta_{\mu_k, \sigma_{k,i}}$ at $x_{k,i}$. Then, by a standard argument, one can show that there is a stepsize $\alpha$ such that (3.11) and (3.12) hold. This shows that the PLS algorithm of Section 3.3 is well defined.
## 4.2 Finite termination

Before proving its finite termination, we give a precise description of the PLS algorithm. The algorithm starts at a point $x_k \in \Omega$ with a positive definite matrix $M_k$. It is assumed that the solution $(d_k, \lambda_k^{\text{QP}})$ of the quadratic program (1.5) is computed in the overall algorithm and that the penalty parameter $\sigma_k$ satisfies the descent condition

$$ \underline{\sigma}_k + \|\lambda_k^{\text{QP}} - \mu_k\|_D \le \sigma_k, $$

for some $\underline{\sigma}_k > 0$ and a multiplier estimate $\mu_k$ given by the overall algorithm. It is also supposed that two constants $\omega_1$ and $\omega_2$ are given in $(0, 1)$ and a constant $\rho$ is given in $(0, \frac{1}{2}]$.

PLS ALGORITHM:

0. Set $i = 0$, $x_{k,0} = x_k$, $d_{k,0} = d_k$, and $\sigma_{k,0} = \sigma_k$.

1. Find a stepsize $\alpha_{k,i}$ such that (3.11) and (3.12) hold for $\alpha = \alpha_{k,i}$. For this do the following:

1.0. Set $j = 0$ and $\alpha_{k,i,0} = 1$.

1.1. If (3.11) and (3.12) hold for $\alpha = \alpha_{k,i,j}$, set $\alpha_{k,i} = \alpha_{k,i,j}$, $x_{k,i+1} = x_{k,i} + \alpha_{k,i}d_{k,i}$, and go to Step 2.

1.2. Choose $\alpha_{k,i,j+1} \in [\rho\alpha_{k,i,j}, (1-\rho)\alpha_{k,i,j}]$,

1.3. Increase $j$ by 1 and go to Step 1.1.

2. If the reduced Wolfe condition (3.13) holds, set $i_k = i + 1$, $x_{k+1} = x_{k,i_k}$, and terminate.

3. Increase $i$ by 1, compute the multiplier estimate $\lambda_k^{\text{QP}}(x_{k,i})$ as the multiplier of problem (2.9) with $x = x_{k,i}$, compute $d_{k,i}$ by (3.10), update $\sigma_{k,i}$ according to the rule given in Section 4.1, and go to 1.

The behavior of the PLS algorithm is analyzed in Theorem 4.4, the proof of which uses the two lemmas below. We recall that a real-valued function $\phi$ is regular at $x$ (in the sense of Clarke [5, Definition 2.3.4]) if it has directional derivatives at $x$, and if for all $h$,

$$ \phi'(x; h) = \limsup_{\substack{x' \to x \\ t \to 0^+}} \frac{\phi(x' + th) - \phi(x')}{t}. $$

**Lemma 4.2.** Suppose that $f$ and $c$ are continuously differentiable on $\Omega$ and let $x \in \Omega$. Suppose also that $x_k \to x$ in $\Omega$, $d_k \to d$ in $\mathbb{R}^n$, and $\alpha_k \to 0$ in $\mathbb{R}_+$. Then, with the merit function $\Theta_{\mu,\sigma}$ defined by (2.18), we have

$$ \Theta'_{\mu,\sigma}(x; d) = \limsup_{k \to \infty} \frac{\Theta_{\mu,\sigma}(x_k + \alpha_k d_k) - \Theta_{\mu,\sigma}(x_k)}{\alpha_k}. $$

**Proof.** Since $\Theta_{\mu,\sigma}$ is Lipschitz continuous in a neighborhood of $x$ and $d_k \to d$, one can readily substitute $d_k$ by $d$ in the equality above, so that we have to prove that $\Theta_{\mu,\sigma}$ is regular. This is the case for $f(\cdot) + \mu^T c(\cdot)$, since this function is continuously
differentiable (use the corollary of Theorem 2.2.1 and Proposition 2.3.6 from [5]). For the regularity of the map $||c(\cdot)||_P$, use [5, Theorem 2.3.10] and the fact that the convexity of the norm implies its regularity [5, Theorem 2.3.6]. □

**Lemma 4.3.** Suppose that $M_k$ is positive definite, that $x \in \Omega \mapsto A(x)$ is bounded, and that the singular values of $A(x)$ are bounded away from zero on $\Omega$. Then, $x \in \Omega \mapsto \hat{A}_k^{-}(x)$ is continuous and bounded.

**Proof.** We have already seen that $u = -\hat{A}_k^{-}(x)c(x)$ is the solution of the quadratic subproblem (2.9), in which $\nabla f(x)$ is set to zero. Therefore there exists a multiplier $\lambda$ such that $(u, \lambda)$ is solution of the corresponding first order optimality conditions:

$$ M_k u + A(x)^\top \lambda = 0 \\ c(x) + A(x) u = 0. $$

Canceling $\lambda$ from these equations and observing that $c(x)$ is an arbitrary vector, the operator $\hat{A}_k^{-}(x)$ can be expressed as follows:

$$ \hat{A}_k^{-}(x) = M_k^{-1} A(x)^\top (A(x) M_k^{-1} A(x)^\top)^{-1}. $$

Then, the continuity and boundedness of $\hat{A}_k^{-}$ follow from the hypotheses. □

**Theorem 4.4.** Suppose that $f$ and $c$ are continuously differentiable on $\Omega$ and that $A(\cdot)$ and $Z^{-}(\cdot)$ are bounded on $\Omega$. Suppose also that $A(\cdot)$ has its singular values bounded away from zero on $\Omega$. If the PLS algorithm is applied from a point $x_k \in \Omega$ with a positive definite matrix $M_k$, then one of the following situations occurs.

(i) The number of iterations of the PLS algorithm is finite, in which case:

(a) $x_{k+1} \in \Omega$,

(b) at each inner iteration, the Armijo condition (3.12) holds,

(c) the reduced Wolfe condition (3.4) holds at $x_{k+1}$.

(ii) The algorithm builds a sequence $\{x_{k,i}\}_i$ in $\Omega$ and

(a) either $\limsup_{i \to \infty} \sigma_{k,i} = +\infty$, in which case $\limsup_{i \to \infty} ||\lambda_k^{\text{QP}}(x_{k,i})||_D = +\infty$,

(b) or $\sigma_{k,i} = \bar{\sigma}_k$ for large $i$, in which case either $\lim_{i \to \infty} \Theta_{\mu_k, \bar{\sigma}_k}(x_{k,i}) = -\infty$ or $\{x_{k,i}\}_i$ converges to a point on the boundary of $\Omega$.

**Proof.** Since $\Theta'_{\mu_k, \sigma_{k,i}}(x_{k,i}; d_{k,i}) < 0$ (Proposition 4.1), conditions (3.11) and (3.12) are satisfied for sufficiently small $\alpha$, and thus the PLS algorithm does not cycle in Step 1. It is also clear that when the number of inner iterations is finite, the conclusions of situation (i) occur.

Note that if $g_k = 0$, then $\hat{Z}_k d_k = -H_k g_k = 0$. In this case, the reduced Wolfe condition (3.13) is trivially satisfied and the algorithm terminates at the first stepsize $\alpha_{k,0}$ satisfying conditions (3.11) and (3.12).
Suppose now that (i) does not occur, then $g_k \neq 0$ and a sequence $\{x_{k,i}\}_i$ is built, such that for $i \geq 1$:

$$x_{k,i} \in \Omega, \qquad (4.6)$$

$$\Theta_{\mu_k, \sigma_{k,i}}(x_{k,i+1}) \leq \Theta_{\mu_k, \sigma_{k,i}}(x_{k,i}) + \omega_1 \alpha_{k,i} \Theta'_{\mu_k, \sigma_{k,i}}(x_{k,i}; d_{k,i}), \qquad (4.7)$$

and

$$g(x_{k,i})^{\top} \hat{Z}_k d_k < \omega_2 g_k^{\top} \hat{Z}_k d_k. \qquad (4.8)$$

Due to (4.4) and (4.5) it follows that either the sequence $\{\sigma_{k,i}\}_i$ is unbounded, which corresponds to conclusion (ii-a), or there exists $i_0$ such that $\sigma_{k,i} = \sigma_{k,i_0} = \bar{\sigma}_k$ for all $i \geq i_0$. It remains to show that in the latter case, the alternative in situation (ii-b) occurs. Up to the end of the proof, we simply denote by $\Theta$ the merit function $\Theta_{\mu_k, \bar{\sigma}_k}$. Let us prove situation (ii-b) by contradiction, assuming that the decreasing sequence $\{\Theta(x_{k,i})\}_i$ is bounded below and that $\{x_{k,i}\}$ does not converge to a point on the boundary of $\Omega$.

Inequality (4.7) implies

$$\Theta(x_{k,i+1}) \leq \Theta(x_{k,i_0}) + \omega_1 \sum_{l=i_0}^{i} \alpha_{k,l} \Theta'(x_{k,l}; d_{k,l}).$$

But, by the positive definiteness of $M_k$ and (4.3)

$$\Theta'(x_{k,i}; d_{k,i}) \leq -\omega_2 g_k^{\top} H_k g_k - \underline{\sigma}_k \|c_{k,i}\|_P.$$

The two latter inequalities, $\omega_1 > 0$, and our assumption on the boundedness from below of $\{\Theta(x_{k,i})\}_i$ imply the convergence of the series

$$\sum_{i \geq 0} \alpha_{k,i} \quad \text{and} \quad \sum_{i \geq 0} \alpha_{k,i} \|c_{k,i}\|_P. \qquad (4.9)$$

In particular, $\alpha_{k,i} \to 0$.

By definition of $x_{k,i}$ we have

$$x_{k,i+1} = x_k + \sum_{l=0}^{i} \left( -\alpha_{k,l} Z_{k,l}^{-\top} H_k g_k - \alpha_{k,l} \hat{A}_{k,l}^{-} c_{k,l} \right).$$

Since $Z^{-}(\cdot)$ and $\hat{A}_{k}^{-}(\cdot)$ are bounded on $\Omega$ (by hypothesis and Lemma 4.3), the convergence of the series (4.9) implies that the series above is absolutely convergent. It follows that the sequence $\{x_{k,i}\}_i$ converges to a point $\bar{x}_k$, which by our assumptions must be in $\Omega$. Therefore, using the continuity of $Z^{-}$, $c$, and $\hat{A}_{k}^{-}$ (Lemma 4.3)

$$d_{k,i} \to \bar{d}_k = -Z^{-}(\bar{x}_k) H_k g_k - \hat{A}_{k}^{-}(\bar{x}_k) c(\bar{x}_k).$$

In Step 1.0, the algorithm takes $\alpha_{k,i,0} = 1$ and we have $\alpha_{k,i} \to 0$. Therefore, for all large $i$, there must exists some index $j_i$ such that $\alpha_{k,i} \in [\rho\alpha_{k,i,j_i}, (1-\rho)\alpha_{k,i,j_i}]$. This means that either (3.11) or (3.12) is not verified for $\alpha = \alpha_{k,i,j_i}$. But for $i$ large,
condition (3.11) holds for $\alpha = \alpha_{k,i,j_i}$, because $x_{k,i} \to \bar{x}_k \in \Omega$, $\{d_{k,i}\}_i$ is bounded, and $\alpha_{k,i,j_i} \to 0$. Therefore, for all large $i$, it is the Armijo condition (3.12) that is not satisfied for $\alpha = \alpha_{k,i,j_i}$. This can be written

$$ \omega_1 \Theta'(x_{k,i}; d_{k,i}) < \frac{\Theta(x_{k,i} + \alpha_{k,i,j_i} d_{k,i}) - \Theta(x_{k,i})}{\alpha_{k,i,j_i}}. $$

When $i \to \infty$, the form of $\Theta'(x_{k,i}; d_{k,i})$ (see for example (2.20)) shows that the left hand side of the inequality above converges to $\omega_1 \Theta'(\bar{x}_k; \bar{d}_k)$. For the right hand side, we use Lemma 4.2, so that by taking the $\limsup_{i \to \infty}$ in the inequality, we obtain $\omega_1 \Theta'(\bar{x}_k; \bar{d}_k) \le \Theta'(\bar{x}_k; \bar{d}_k)$. Since $\omega_1 < 1$, this implies $\Theta'(\bar{x}_k; \bar{d}_k) \ge 0$.

On the other hand,

$$ \Theta'(x_{k,i}; d_{k,i}) \le -\omega_2 g_k^\top H_k g_k - \underline{\sigma}_k \|c_{k,i}\|_P. $$

Taking the limit in this inequality when $i \to \infty$ and recalling that $g_k \ne 0$, we obtain

$$ \Theta'(\bar{x}_k; \bar{d}_k) \le -\omega_2 g_k^\top H_k g_k - \underline{\sigma}_k \|c(\bar{x}_k)\|_P < 0, $$

a contradiction that concludes the proof. $\square$

# 5 Convergence results

In this section we give a global convergence result, assuming the boundedness of the generated matrices $M_k$ and their inverse. Despite this strong assumption, we believe that such a result is useful in that it shows that the different facets of the algorithm introduced in the previous sections can fit together.

The results given in this section deal with the behavior of the sequence of iterates $x_k$, so that it is implicitly assumed that a sequence $\{x_k\}$ is generated and therefore that the PLS algorithm has finite termination each time it is invoked.

Recall that the value of the penalty parameter $\sigma_k$ must be updated such that the descent condition

$$ \underline{\sigma}_k + \|\lambda_k^{\text{QP}} - \mu_k\|_D \le \sigma_k \qquad (5.1) $$

is satisfied. This corresponds to (4.1) with $i = 0$. In (5.1), $\underline{\sigma}_k$ is a positive number that is adapted at some iterations.

## 5.1 Admissibility of the unit stepsize

Admissibility of the unit stepsize by Armijo's condition on the penalty function $\Theta_{\mu_k, \sigma_k}$ is studied by Bonnans [2], but in a form that is not adapted to our algorithm. Proposition 5.1 gives a version suitable to us. Conditions for admissibility of the unit stepsize by the reduced Wolfe inequality are given in Proposition 5.2. These results are obtained by expanding $f$ and $c$ about the current iterate $x_k$. They are useful for determining how and when the multiplier $\mu_k$ and the penalty parameter $\sigma_k$ have to be adapted.
Proposition 5.1 requires that the multiplier estimate $\lambda_k^{QP}$ be used in the descent condition (5.1). It also requires that $\mu_k$ be sufficiently close to the optimal multiplier and that the penalty parameter be sufficiently small. In other words, near the solution, the merit function has to be sufficiently close to the Lagrangian function. This implies that, in the overall algorithm below, $\mu_k$ will have to be reset to $\lambda_k^{QP}$ and $\sigma_k$ will have to be decreased at some iterations.

**Proposition 5.1.** Suppose the $f$ and $c$ are twice continuously differentiable in a convex neighborhood of a local solution $x_*$ satisfying the second order sufficient condition of optimality (1.3)-(1.4). Suppose also that $x_k \to x_*$, $d_k \to 0$, $\omega_1 < 1/2$, the descent condition (5.1) holds, and

$$d_k^T(M_k - L_*^r)d_k = o(\|d_k\|^2), \qquad (5.2)$$

in which $r$ is a nonnegative scalar such that $L_*^r$ is positive definite. Then, there exists a constant $\varepsilon > 0$ such that when

$$\|\mu_k - \lambda_*\| \le \varepsilon \quad \text{and} \quad 0 \le \sigma_k \le \varepsilon,$$

and when $k$ is sufficiently large, the unit stepsize is accepted by the Armijo inequality:

$$\Theta_{\mu_k, \sigma_k}(x_k + d_k) \le \Theta_{\mu_k, \sigma_k}(x_k) + \omega_1 \Theta'_{\mu_k, \sigma_k}(x_k; d_k).$$

**Proof.** Since $d_k \to 0$, a second order expansion of $f(x_k + d_k)$ about $x_k$ gives with (1.6)

$$\begin{aligned} f(x_k + d_k) &= f_k + \nabla f_k^T d_k + \frac{1}{2} d_k^T \nabla^2 f(x_*) d_k + o(\|d_k\|^2) \\ &= f_k - d_k^T M_k d_k + (\lambda_k^{QP})^T c_k + \frac{1}{2} d_k^T \nabla^2 f(x_*) d_k + o(\|d_k\|^2). \end{aligned}$$

Similarly, for any component $c_{(i)}$ of $c$, we have with (1.6)

$$c_{(i)}(x_k + d_k) = \frac{1}{2} d_k^T \nabla^2 c_{(i)}(x_*) d_k + o(\|d_k\|^2).$$

Combining these two estimates and using (1.6), (2.21), and the hypotheses on $\mu_k$ and $\sigma_k$, we get

$$\begin{aligned} \Theta_{\mu_k, \sigma_k}(x_k + d_k) - \Theta_{\mu_k, \sigma_k}(x_k) - \omega_1 \Theta'_{\mu_k, \sigma_k}(x_k; d_k) &= -d_k^T M_k d_k + (\lambda_k^{QP} - \mu_k)^T c_k - \sigma_k \|c_k\|_P + \frac{1}{2} d_k^T L_* d_k \\ &\quad - \omega_1 \Theta'_{\mu_k, \sigma_k}(x_k; d_k) + (\|\mu_k - \lambda_*\| + \sigma_k) O(\|d_k\|^2) + o(\|d_k\|^2) \\ &\le (1 - \omega_1) \Theta'_{\mu_k, \sigma_k}(x_k; d_k) + \frac{1}{2} d_k^T L_*^r d_k - \frac{r}{2} \|c_k\|_2^2 \\ &\quad + \varepsilon O(\|d_k\|^2) + o(\|d_k\|^2). \end{aligned}$$
Now, splitting $(1-\omega_1)$ in $(\frac{1}{2}-\omega_1)+\frac{1}{2}$, using the fact that $\Theta'_{\mu_k, \sigma_k}(x_k; d_k) \le -d_k^T M_k d_k$ (see (4.2)), the nonnegativity of $r$, $\omega_1 < 1/2$, the positive definiteness of $L_*^r$ (which with (5.2) implies that $d_k^T M_k d_k \ge C' \|d_k\|^2$, for some constant $C' > 0$), and (5.2), we obtain

$$ \begin{aligned} & \Theta_{\mu_k, \sigma_k}(x_k + d_k) - \Theta_{\mu_k, \sigma_k}(x_k) - \omega_1 \Theta'_{\mu_k, \sigma_k}(x_k; d_k) \\ & \le \left(\frac{1}{2} - \omega_1\right) (-d_k^T M_k d_k) - \frac{1}{2} d_k^T (M_k - L_*^r) d_k + \varepsilon O(\|d_k\|^2) + o(\|d_k\|^2) \\ & \le -C \|d_k\|^2 + \varepsilon O(\|d_k\|^2) + o(\|d_k\|^2), \end{aligned} $$

for some constant $C > 0$. Since the last right hand side is negative when $k$ is large and $\varepsilon$ is sufficiently small, the proposition is proved. □

Proposition 5.1 suggests a way of updating the parameters $\mu_k$, $\sigma_k$, and $\underline{\sigma}_k$: $\mu_k$ should be close to $\lambda_k^{\text{QP}}$ and $\sigma_k$ should be kept small. The latter condition may require to decrease $\underline{\sigma}_k$.

In order to ensure convergence, we allow $\mu_k$ to change only when the iterates progress sufficiently to a local solution. This is measured by the following quantity:

$$ \varepsilon_k = \min(\|g_k\| + \|c_k\|_P, \varepsilon_{k-1}), \quad k \ge 0 \qquad (5.3) $$

($\varepsilon_{-1} = \|g_0\| + \|c_0\|_P$). It follows that the sequence $\{\varepsilon_k\}$ is nonincreasing and that $\lim_{k \to \infty} \varepsilon_k = 0$ if and only if $\liminf_{k \to \infty} \|g_k\| + \|c_k\|_P = 0$.

Now, suppose that a new iterate $x_{k+1}$ has been computed by the PLS algorithm. Recall that $\alpha_{k,0}$ is the first stepsize along the direction $d_k$ at which the Armijo condition (3.12) is satisfied. Our update rule for $\underline{\sigma}_k$ uses two constants $a_1 > 1$ and $a_2 > 1$. Let $k'$ be the index of the last iteration at which $\underline{\sigma}_k$ has been updated. Initially $k'$ is set to 0 in the overall algorithm.

UPDATE OF $\underline{\sigma}_k$:

if $\varepsilon_{k+1} \le \varepsilon_{k'}/a_1$ and $\alpha_{k,0} \ne 1$

then $k' = k + 1$ and $\underline{\sigma}_{k+1} = \underline{\sigma}_k/a_2$

else $\underline{\sigma}_{k+1} = \underline{\sigma}_k$.

In other words, $\underline{\sigma}_k$ is decreased when the iterates progress to a local solution, although the unit stepsize is not accepted along the direction $d_k$.

Recall that the value $\sigma_{k,i_{k-1}}$ used in the update rule below is the value of the penalty parameter at the end of the PLS. Let $a_3 > 1$ be another constant. The rule below updates an index $k''$, initially set to 0 in the overall algorithm. It is the index of the last iteration at which $\mu_k$ has been set to $\lambda_k^{\text{QP}}$.

UPDATE OF $\mu_k$ AND $\sigma_k$:

if $\varepsilon_{k+1} \le \varepsilon_{k''}/a_3$ then $k'' = k + 1$, $\mu_{k+1} = \lambda_{k+1}^{\text{QP}}$, and $\sigma_{k+1} = \underline{\sigma}_{k+1}$,
else $\mu_{k+1} = \mu_k$ and set $\sigma_{k+1}$ according to:

if $\underline{\sigma}_{k+1} + \|\lambda_{k+1}^{\text{QP}} - \mu_{k+1}\|_D \le \sigma_{k,i_k-1}$ then $\sigma_{k+1} = \sigma_{k,i_k-1}$,
else $\sigma_{k+1} = \max(2\sigma_{k,i_k-1}, \underline{\sigma}_{k+1} + \|\lambda_{k+1}^{\text{QP}} - \mu_{k+1}\|_D)$.

Note that, when $\mu_k$ is unchanged, the parameter $\sigma_k$ is updated by the same rule as $\sigma_{k,i}$ in Section 4.1. By using the notation $\sigma_{k+1,-1} = \sigma_{k,i_k-1}$, it follows that as long as $\mu_k$ is kept constant, the sequence $\sigma_{k,-1}$, $\sigma_{k,0}$ ($= \sigma_k$), $\sigma_{k,1}$, ..., $\sigma_{k+1,-1}$, $\sigma_{k+1,0}$, $\sigma_{k+1,1}$, ..., is nondecreasing and satisfies (4.4) or (4.5).

Let us consider now, the admissibility of the unit stepsize for the reduced Wolfe condition.

**Proposition 5.2.** Suppose that $g$ is continuously differentiable in a convex neighborhood of $x_*$ satisfying the optimality condition (1.3). Suppose also that $x_k \to x_*$, $d_k \to 0$, $\omega_2 > 0$, $\{M_k\}$ and $\{M_k^{-1}\}$ are bounded, $c_k = O(\|g_k\|)$, and

$$ (Z_k^{-\top} M_k - Z_*^{-\top} L_*) d_k = o(\|d_k\|). \qquad (5.4) $$

Then, when $k$ is sufficiently large, the unit stepsize is accepted by the reduced Wolfe condition:

$$ g(x_k + d_k)^{\top} \hat{Z}_k d_k \ge \omega_2 g(x_k)^{\top} \hat{Z}_k d_k. $$

**Proof.** From (1.6) we have $Z_k^{-\top} M_k d_k = -g_k$ and by the hypotheses on $M_k$ and $c_k$, $d_k = O(\|g_k\|)$. Then, using (2.6) and next (5.4), we obtain

$$ \begin{aligned} g(x_k + d_k) &= g_k + \nabla g(x_*)^{\top} d_k + o(\|d_k\|) \\ &= -(Z_k^{-\top} M_k - Z_*^{-\top} L_*) d_k + o(\|d_k\|) \\ &= o(\|g_k\|). \end{aligned} $$

Since $\hat{Z}_k d_k = -H_k g_k = O(\|g_k\|)$, we finally have

$$ g(x_k + d_k)^{\top} \hat{Z}_k d_k - \omega_2 g_k^{\top} \hat{Z}_k d_k = \omega_2 g_k^{\top} H_k g_k + o(\|g_k\|^2), $$

which is positive for $k$ large, by the uniform positive definiteness of $\{H_k\}$ and $\omega_2 > 0$. $\square$

Proposition 5.2 suggests to use a criterion for deciding when launching the PLS algorithm. Suppose indeed that, contrary to what is required by the hypotheses of this proposition $c_k \ne O(\|g_k\|)$ or equivalently that a subsequence of $\{\|g_k\|/\|c_k\|\}$ tends to zero. This means that the iterates $x_k$ approach $x_*$ tangently to the reduced gradient manifold $\{g = 0\}$. In this case, the longitudinal part of $M_k$ is less important and updating $M_k$ by adjusting the parameter $r_k$ only looks perfectly adequate. It is also feasible, since $c_k \ne 0$ for this subsequence. On the other hand, Proposition 5.2 tells us that using the PLS algorithm in this case may prevent the unit stepsize from being accepted at each iteration, which may prevent superlinear convergence from occurring. These two remarks suggest not launching a PLS in this case. Therefore, choosing some constant $K > 0$, we adopt the following criterion.
PLS CRITERION:

if $||c_k|| \le K||g_k||$ **then** call the PLS algorithm,
else perform only the first inner iteration of the PLS algorithm.

It follows that whenever $||c_k|| > K||g_k||$, the next iterate $x_{k+1}$ is set to $x_{k,1}$, which only satisfies condition (3.11) and the Armijo inequality (3.12), but not necessarily the reduced Wolfe condition (3.4).

Note finally that the conditions (5.2) and (5.4) on the updated matrix $M_k$, used in Propositions 5.1 and 5.2 are both satisfied when $(M_k - L_*^r)d_k = o(||d_k||)$, which is a reasonable condition to expect from the quasi-Newton theory.

## 5.2 Global convergence

We are now in position to give a complete description of our algorithm.

OVERALL ALGORITHM

0. Choose some constants $a_i > 1$ ($i = 1, \dots, 3$) for the update of $\mu_k$, $\sigma_k$, and $\underline{\sigma}_k$; a constant $K > 1$ for the PLS criterion; constants $\omega_1 \in (0, \frac{1}{2})$, $\omega_2 \in (0, 1)$, and $\rho \in (0, \frac{1}{2}]$ for the PLS algorithm; and a constant $\beta \in (0, 1)$ for the update of $r_k$.

Choose a starting point $x_0 \in \Omega$ and an initial symmetric positive definite matrix $M_0 \in \mathbb{R}^{n \times n}$.

Set $k = k' = k'' = 0$ (the indices $k'$ and $k''$ are reset by the update rules of $\mu_k$, $\sigma_k$, and $\underline{\sigma}_k$).

Solve the SQP subproblem (1.5) (with $k = 0$) giving $(d_0, \lambda_0^{\text{QP}})$.

Choose $\underline{\sigma}_0 > 0$, set $\mu_0 = \lambda_0^{\text{QP}}$ and $\sigma_0 = \underline{\sigma}_0$.

1. if $||c_k|| \le K||g_k||$ **then** call the PLS algorithm,
else perform only the first inner iteration of the PLS algorithm.
This gives a new iterate $x_{k+1}$.

2. Compute $\gamma_k$ and $\delta_k$ by formula (3.2) and (3.14), where $r_k$ is obtained as described in Section 3.5, and update $M_{k+1}$ by the BFGS formula (1.7).

3. Solve the SQP subproblem giving $(d_{k+1}, \lambda_{k+1}^{\text{QP}})$.

4. Update $\mu_{k+1}$, $\sigma_{k+1}$, and $\underline{\sigma}_{k+1}$ by the rules given in Section 5.1.

5. Increase $k$ by 1 and go to 1.

Below, we denote by $\text{dist}(x, \Omega^c)$ the Euclidean distance between a point $x$ and the complementary set of $\Omega$.

**Theorem 5.3.** Suppose that $\Omega$ is convex, that $f$ and $c$ are differentiable on $\Omega$ with Lipschitz continuous derivatives, and that $Z^{-}(\cdot)$ is bounded on $\Omega$. If the overall algorithm above generates a sequence $\{x_k\}$, using a bounded sequence of matrices $\{M_k\}$ with bounded inverses, then one of the following situations occurs.

(i) The algorithm converges in the sense that

$$ \liminf_{k \to \infty} (||g_k|| + ||c_k||_P) = 0. $$
(ii) There exists $k_0$ such that $\mu_k = \mu$, for all $k \ge k_0$ and

(a) either the set $\{\sigma_{k,i} : k \ge 0, 0 \le i < i_k\}$ is unbounded, implying that the set $\{\lambda_k^{\text{QP}}(x_{k,i}) : k \ge 0, 0 \le i < i_k\}$ is also unbounded,

(b) or there exists $k_1 \ge k_0$ such that $\sigma_{k,i} = \sigma$, for all $k \ge k_1$ and $0 \le i < i_k$, in which case $\Theta_{\mu,\sigma}(x_k) \to -\infty$ or $\liminf_{k \to \infty} \text{dist}(x_k, \Omega^c) = 0$.

**Proof.** To prove the first part, suppose that conclusion (i) does not occur. By the definition (5.3) of $\varepsilon_k$, this is equivalent to $\lim_{k \to \infty} \varepsilon_k > 0$. The update rules of $\mu_k$, $\sigma_k$ and $\underline{\sigma}_k$, given in Section 5.1, and (4.4) and (4.5) imply that there exists an index $k_0$, such that for all $k \ge k_0$, $\mu_k = \mu$, $\underline{\sigma}_k = \underline{\sigma}$ and the sequence $\sigma_{k,0}(= \sigma_k)$, $\sigma_{k,1}$, ..., $\sigma_{k,i_k-1}$, $\sigma_{k+1,0}(= \sigma_{k+1})$, $\sigma_{k+1,1}$, ..., is nondecreasing. This sequence is either unbounded, in which case conclusion (ii-a) follows, or there exists $k_1 \ge k_0$ such that $\sigma_{k,i} = \sigma$ for all $k \ge k_1$ and $0 \le i < i_k$.

It remains to prove that in the latter case the alternative given in (ii-b) holds. This is done by contradiction, assuming that $\Theta_{\mu,\sigma}(x_k)$ is bounded from below and that $\liminf_{k \to \infty} \text{dist}(x_k, \Omega^c) > 0$.

Since the Armijo inequality (3.12) holds at each inner iteration of the PLS algorithm (conclusion (i-b) of Theorem 4.4), we have

$$ \Theta_{\mu,\sigma}(x_{k+1}) \le \Theta_{\mu,\sigma}(x_k) + \omega_1 \sum_{i=0}^{i_k-1} \alpha_{k,i} \Theta'_{\mu,\sigma}(x_{k,i}; d_{k,i}). $$

Note that this inequality holds even if the PLS criterion is not satisfied (Step 1 of the overall algorithm). Recall that in our notation, $x_{k,0} = x_k$ and $d_{k,0} = d_k$. By using the fact that $d_{k,i}$ is a descent direction of $\Theta_{\mu,\sigma}$ at $x_{k,i}$, the previous inequality implies

$$ \Theta_{\mu,\sigma}(x_{k+1}) \le \Theta_{\mu,\sigma}(x_k) + \omega_1 \alpha_{k,0} \Theta'_{\mu,\sigma}(x_k; d_k). $$

From inequality (4.2) we have $\Theta'_{\mu,\sigma}(x_k; d_k) \le -d_k^T M_k d_k - \underline{\sigma} \|c_k\|_P$, so that the inequality above implies

$$ 0 \le \omega_1 \alpha_{k,0} (d_k^T M_k d_k + \underline{\sigma} \|c_k\|_P) \le \Theta_{\mu,\sigma}(x_k) - \Theta_{\mu,\sigma}(x_{k+1}). $$

Summing over $k$ and using the boundedness assumption on $\Theta_{\mu,\sigma}(x_k)$, we deduce the convergence of the series

$$ \sum_{k \ge 0} \alpha_{k,0} d_k^T M_k d_k < +\infty \quad \text{and} \quad \sum_{k \ge 0} \alpha_{k,0} \|c_k\|_P < +\infty. \qquad (5.5) $$

If $\liminf \alpha_{k,0} > 0$, then the convergence of these series and the boundedness of $\{M_k^{-1}\}$ would imply that $d_k \to 0$ and $c_k \to 0$, and since $g_k = -Z_k^{-T} M_k d_k$ (see (2.16)) we would have $(\|g_k\| + \|c_k\|_P) \to 0$, in contradiction with our initial assumption. Thus, a subsequence $\{\alpha_{k,0}\}_{k \in K}$ converges to 0. This means that either condition (3.11) or the Armijo condition (3.12) is not accepted for $\alpha = \alpha_{k,0,j_k} \le \rho^{-1}\alpha_{k,0}$ (Step 1 of the PLS algorithm). This can be written

$$ x_k + \alpha_{k,0,j_k} d_k \notin \Omega \qquad (5.6) $$
or

$$ \Theta_{\mu,\sigma}(x_k + \alpha_{k,0,j_k} d_k) > \Theta_{\mu,\sigma}(x_k) + \omega_1 \alpha_{k,0,j_k} \Theta'_{\mu,\sigma}(x_k; d_k). \quad (5.7) $$

Let us show that (5.6) does not hold for large $k$. Using the convergence of the first series in (5.5), the boundedness of {$M_k^{-1}$}, and $\alpha_{k,0} \le 1$, we have $\alpha_{k,0}||d_k|| \to 0$. But for $k \in \mathcal{K}$, $\alpha_{k,0,j_k} \le \rho^{-1}\alpha_{k,0}$, so that $\alpha_{k,0,j_k}||d_k|| \to 0$, which with (5.6) implies that $\liminf dist(x_k, \Omega^c) \to 0$, in contradiction with our assumptions.

Now, we show that (5.7) leads to a contradiction, which will prove the first part of the theorem. Expanding $f$ and $c$ at the first order and using the Lipschitz continuity of $\nabla f$ and $\nabla c$, we obtain the following estimates, when $\alpha \in (0,1]$ and $x_k + \alpha d_k \in \Omega$:

$$ f(x_k + \alpha d_k) \le f_k + \alpha \nabla f_k^\top d_k + C \alpha^2 ||d_k||^2, $$

$$ \mu^\top c(x_k + \alpha d_k) \le (1 - \alpha) \mu^\top c_k + C \alpha^2 ||d_k||^2, $$

and

$$ ||c(x_k + \alpha d_k)||_P \le (1 - \alpha) ||c_k||_P + C \alpha^2 ||d_k||^2, $$

where $C$ denotes a constant independent of $k$. Then, we deduce

$$ \begin{aligned} \Theta_{\mu,\sigma}(x_k + \alpha d_k) &= f(x_k + \alpha d_k) + \mu^\top c(x_k + \alpha d_k) + \sigma ||c(x_k + \alpha d_k)||_P \\ &\le f_k + \alpha \nabla f_k^\top d_k + (1 - \alpha) \mu^\top c_k + (1 - \alpha) \sigma ||c_k||_P + C \alpha^2 ||d_k||^2 \\ &\le \Theta_{\mu,\sigma}(x_k) + \alpha \Theta'_{\mu,\sigma}(x_k; d_k) + C \alpha^2 ||d_k||^2. \end{aligned} $$

The last inequality and (5.7) imply

$$ 0 < (1 - \omega_1) \Theta'_{\mu,\sigma}(x_k; d_k) + C \alpha_{k,0,j_k} ||d_k||^2. $$

Finally, using the inequality $\Theta'_{\mu,\sigma}(x_k; d_k) \le -d_k^\top M_k d_k - \underline{\sigma} ||c_k||_P$, the boundedness of {$M_k^{-1}$}, and next $\alpha_{k,0,j_k} \to 0$, when $k \to \infty$ in $\mathcal{K}$, we obtain a contradiction

$$ 0 < -(1 - \omega_1) \underline{\sigma} ||c_k||_P - C ||d_k||^2 \le 0. $$

This contradiction concludes the proof. □

We conclude this section by a result specifying the conditions under which the unit step-size is accepted asymptotically in the overall algorithm.

**Proposition 5.4.** Suppose that $f$ and $c$ are twice continuously differentiable and $g$ is continuously differentiable on a convex neighborhood of a point $x_*$ satisfying (1.3) and (1.4). If the overall algorithm generates a sequence {$x_k$} converging to $x_*$, such that $d_k \to 0$ and $(M_k - L_*^r)d_k = o(||d_k||)$ for some $r \ge 0$ such that $L_*^r$ is positive definite, then, for sufficiently large $k$, there is only one inner iteration in the PLS algorithm (i.e., $i_k = 1$) and the unit stepsize is accepted by the Armijo inequality (i.e., $\alpha_{k,0} = 1$).

**Proof.** Since $x_k$ converges to a local solution of the problem, then $g_k \to 0$ and $c_k \to 0$, and thus $\varepsilon_k \to 0$.
Suppose now that $\alpha_{k,0} \neq 1$ for a subsequence of iterates. Then, the update rule of $\underline{\sigma}_k$ implies the convergence of the sequence $\{\underline{\sigma}_k\}$ to 0 (it is here that the index $k'$ is useful). In the same way, the update rule of $\mu_k$ and $\sigma_k$ implies $\mu_k \to \lambda_*$ and $\sigma_k \to 0$ (usefulness of the index $k''$). It follows from Proposition 5.1 that $\alpha_{k,0} = 1$ for large $k$, contradicting our initial assumption.

Also $i_k = 1$ for large $k$, either because the PLS criterion does not hold or because of Proposition 5.2. □

# 6 Implementation issues

In this section, we discuss some issues related to the implementation of the algorithm described in the previous sections.

## 6.1 Using a QR factorization of $A^T$

In our experiment, the matrices $A^-(x)$ and $Z^-(x)$ described in Section 2 are obtained from a QR factorization of $A(x)^T$:

$$A(x)^T = (Y^-(x) \quad Z^-(x)) \begin{pmatrix} R(x) \\ 0 \end{pmatrix} = Y^-(x)R(x),$$

where $Y^-(x)$ and $Z^-(x)$ are respectively $n \times m$ and $n \times (n-m)$ matrices, such that $(Y^-(x) \quad Z^-(x))$ is orthogonal, and $R(x)$ is an order $m$ upper triangular matrix. Clearly, the columns of $Z^-(x)$ span the null space of $A(x)$, as desired. We choose as right inverse of $A(x)$ the Moore-Penrose pseudo-inverse $A(x)^T(A(x)A(x)^T)^{-1}$, which can be computed by

$$A^-(x) = Y^-(x)R(x)^{-T}.$$

With this choice for $Z^-$ and $A^-$, it follows that $Z(x) = Z^-(x)^T$.

## 6.2 Weighted augmentation

Byrd, Tapia and Zhang [4, p. 216] emphasize that, due to the augmentation term $A_k^T A_k$ in the vector $\gamma_k$, badly scaled constraints may have some negative effects on the next updated matrix. These authors prefer using a weighted augmentation term $A_k^T W_k A_k$, where $W_k$ is an order $m$ weighting matrix. They suggest to use $W_k = (A_k A_k^T)^{-1}$, because with the notation of Section 6.1

$$A_k^T W_k A_k = A_k^T (A_k A_k^T)^{-1} A_k = Y_k^{-} Y_k^{-T},$$

which is a well conditioned matrix.

It is clear that the same technique can be adopted in our algorithm. Therefore, in our experiment we replace the augmentation term $A_k^T A_k \delta_k$ by $Y_k^{-} Y_k^{-T} \delta_k$.
## 6.3 Scaling the tangential direction

It is shown in [12] that the number of inner iterations in the PLS algorithm can be reduced by using a scaling factor $\tau_{k,i} > 0$ on the longitudinal component of the inner search direction $d_{k,i}$. This leads to redefine the direction $d_{k,i}$ as

$$d_{k,i} = -\tau_{k,i} Z_{k,i}^{-} H_k g_k - \hat{A}_{k,i}^{-} c_{k,i}.$$

For $i = 0$, we set $\tau_{k,0} = 1$, so that when a unit stepsize is accepted ($i_k = 1$ and $\alpha_{k,0} = 1$), a plain SQP step is taken and superlinear convergence of the overall algorithm may occur.

With this change, the vector $\delta_k$ is still given by (3.14), but the longitudinal stepsize $\alpha_k^z$ has to be computed by

$$\alpha_k^z = \sum_{i=0}^{i_k-1} \alpha_{k,i} \tau_{k,i}.$$

It is not difficult to extend the finite termination result of Section 4.2 (Theorem 4.4), when $d_{k,i}$ is given as above, provided $\tau_{k,i}$ is maintained in a fixed interval: $0 < \underline{\tau}_k \le \tau_{k,i} \le \overline{\tau}_k$ for all $i \ge 0$. The global convergence result (Theorem 5.3) is not affected by the scaling factors $\tau_{k,i}$, since $\tau_{k,0} = 1$ and only the progress obtained by first inner iteration of the PLS is used in the convergence proof.

## 6.4 Speeded-up PLS technique

There is another way of speeding up the PLS algorithm that is useful in practice to reduce the number of inner iterations (this is discussed in [12]). It consists in resetting the right hand side of the Wolfe inequality (3.13) to $\omega_2 g_{k,i}^\top \hat{Z}_k d_k$, when the current iterate $x_{k,i}$ makes $g_{k,i}^\top \hat{Z}_k d_k$ more negative than at all the previous inner iterations. Synthetically, it consists in replacing the reduced Wolfe condition (3.4) by the less demanding inequality

$$g_{k+1}^\top \hat{Z}_k d_k \ge \omega_2 \min_{0 \le i < i_k} g_{k,i}^\top \hat{Z}_k d_k. \quad (6.1)$$

Since this inequality is more rapidly satisfied than (3.4), the finite termination result (Theorem 4.4) is still valid. Also, the global convergence theorem still holds, since the wolfe condition (3.4) plays no role in its proof.

Let us denote by $l_k$ the greatest index $i$ realizing the minimum in the right hand side of (6.1). When this condition is used in place of (3.4) and $l_k \ne 0$, the vectors $\gamma_k$ and $\delta_k$ have to be modified. Aiming at having $\delta_k \approx x_{k+1} - x_{k,l_k}$, the same reasoning as in Sections 3.1 and 3.4 shows that it is appropriate to set (we take into account the suggestions given in Sections 6.2 and 6.3):

$$\gamma_k = Z_{k,l_k}^\top (g_{k+1} - g_{k,l_k}) - A_{k,l_k}^\top (\lambda_{k+1} - \lambda_{k,l_k}) + r_k Y_{k,l_k}^{-} Y_{k,l_k}^{-\top} \delta_k$$

$$\delta_k = -\alpha_{k,l_k}^z Z_{k,l_k}^{-} H_k g_k - \alpha_{k,l_k}^\Lambda \hat{A}_{k,l_k}^{-} c_{k,l_k},$$
where $\alpha_{k,l_k}^z$ and $\alpha_{k,l_k}^A$ are defined by

$$ \alpha_{k,l_k}^z = \sum_{i=l_k}^{i_k-1} \alpha_{k,i} \tau_{k,i} \quad \text{and} \quad \alpha_{k,l_k}^A = \sum_{i=l_k}^{i_k-1} \alpha_{k,i} e^{-\xi_{k,i}}, $$

with $\xi_{k,i} = \sum_{j=l_k}^{i-1} \alpha_{k,j}$. Note that when $l_k = 0$ and $\tau_{k,i} = 1$, we recover the preceding formulae (3.2), (3.14), and (3.15).

Observe finally that when $A_{k,l_k} \delta_k = 0$, then $c_{k,l_k} = 0$ and we have

$$ \begin{aligned} \gamma_k^T \delta_k &= (g_{k+1} - g_{k,l_k})^T Z_{k,l_k} \delta_k \\ &= \alpha_{k,l_k}^z (g_{k+1} - g_{k,l_k})^T \hat{Z}_k d_k \\ &> 0, \end{aligned} $$

by (6.1). Hence, one can always have $\gamma_k^T \delta_k > 0$ by adjusting the value of $r_k$.

# 7 Numerical experiment

Our experiments were performed on a Power Macintosh 6100/66 in Matlab (release 4.2c.1) with a machine epsilon of about $2 \times 10^{-16}$. We have taken the same list of test problems as in the paper of Byrd, Tapia, and Zhang [4]. In each case, the standard starting points (those given in [17, 25]) were used, except for problems 12, 316-322, 336 and 338, for which we have used $x_0 = 10^{-4}(1, \dots, 1)$, since the Jacobian matrix is rank deficient at the standard initial point $x_0 = 0$.

Three updating methods have been tested:

* **Powell**: Powell's corrections, in which, at each iteration, $\gamma_k$ is set to $\gamma_k^P$ given by formula (1.10), with $\theta$ calculated as described after this formula ($\eta = 0.1$);

* **BTZ**: the Byrd, Tapia, and Zhang approach, in which $\gamma_k$ is given by formula (1.11), the rules (1.12) and (1.13), with $\nu_{\text{BTZ}} = \beta_{\text{BTZ}} = 0.01$, and the weighted augmentation described in Section 6.2;

* **PLS**: the algorithm presented in this paper.

All the methods use the same merit function (2.18), with $\| \cdot \|_P = \| \cdot \|_1$, and the same technique to update the parameters $\mu$, $\sigma$, and $\underline{\sigma}$. The constants, used by the update rules (see Section 5.1), are set as follows: $a_1 = 1.0001$, $a_2 = 10$, $a_3 = 1.0001$, and $\underline{\sigma}_0 = \|\lambda_0\|_{\infty}/10$. The constant used in the Armijo inequality is set to $\omega_1 = 0.01$. For **Powell** and **BTZ** algorithms, each stepsize $\alpha_k$ is found by a backtracking line-search along $d_k$, as in Step 1 of the PLS algorithm. In Step 1.2 of the PLS algorithm, quadratic or cubic interpolation formulae are used and we set $\rho = 0.1$. The constant for the reduced Wolfe condition is $\omega_2 = 0.9$ and the one used in the PLS criterion is $K = 100\|c_0\|/\|g_0\|$. Initially, $M_0 = I$ and the pre-update scaling $\gamma_k^T \delta_k / \|\delta_k\|^2 I$ is done before the first update by the BFGS formula. The stopping criterion for all the methods is

$$ \|g_k\|_2 + \|c_k\|_2 \le \varepsilon_{\text{tol}} (\|g_0\|_2 + \|c_0\|_2), \quad \text{with} \quad \varepsilon_{\text{tol}} = 10^{-7}. $$
<table><caption>Table 1: Results.</caption><thead><tr><th rowspan="2">P</th><th rowspan="2">n:m</th><th colspan="3">Powell</th><th colspan="3">BTZ</th><th colspan="3">PLS</th></tr><tr><th>ng/nf</th><th>cr</th><th>κ<sub>2</sub></th><th>ng/nf</th><th>cr</th><th>κ<sub>2</sub></th><th>ng/nf</th><th>cr</th><th>κ<sub>2</sub></th><th>ii</th></tr></thead><tbody><tr><td>26</td><td>3:1</td><td>25/26</td><td>-</td><td>4.</td><td>25/26</td><td>-</td><td>4.</td><td>25/26</td><td>1</td><td>5.</td><td>-</td></tr><tr><td>43</td><td>4:2</td><td>11/12</td><td>1</td><td>1.</td><td>11/12</td><td>-</td><td>1.</td><td>9/10</td><td>-</td><td>0.</td><td>-</td></tr><tr><td>46</td><td>5:2</td><td>36/40</td><td>-</td><td>6.</td><td>36/40</td><td>-</td><td>6.</td><td>36/40</td><td>-</td><td>6.</td><td>-</td></tr><tr><td>47</td><td>5:3</td><td>29/33</td><td>7</td><td>6.</td><td>26/29</td><td>2</td><td>7.</td><td>28/30</td><td>-</td><td>4.</td><td>-</td></tr><tr><td>56</td><td>7:4</td><td>**</td><td>5</td><td>13.</td><td>10/11</td><td>2</td><td>2.</td><td>13/15</td><td>2</td><td>2.</td><td>-</td></tr><tr><td>60</td><td>3:1</td><td>9/10</td><td>-</td><td>1.</td><td>9/10</td><td>-</td><td>1.</td><td>10/11</td><td>-</td><td>1.</td><td>-</td></tr><tr><td>63</td><td>3:2</td><td>12/13</td><td>7</td><td>10.</td><td>9/9</td><td>2</td><td>2.</td><td>8/8</td><td>5</td><td>0.</td><td>-</td></tr><tr><td>66</td><td>3:2</td><td>8/8</td><td>-</td><td>1.</td><td>8/8</td><td>-</td><td>1.</td><td>8/8</td><td>-</td><td>2.</td><td>-</td></tr><tr><td>78</td><td>5:3</td><td>9/10</td><td>1</td><td>2.</td><td>9/10</td><td>1</td><td>1.</td><td>8/9</td><td>1</td><td>1.</td><td>-</td></tr><tr><td>79</td><td>5:3</td><td>12/13</td><td>1</td><td>2.</td><td>12/13</td><td>-</td><td>2.</td><td>12/13</td><td>-</td><td>2.</td><td>-</td></tr><tr><td>100</td><td>7:2</td><td>14/17</td><td>-</td><td>2.</td><td>14/17</td><td>-</td><td>2.</td><td>12/15</td><td>-</td><td>2.</td><td>-</td></tr><tr><td>373</td><td>9:6</td><td>*</td><td>7</td><td>20.</td><td>27/29</td><td>1</td><td>14.</td><td>16/18</td><td>-</td><td>9.</td><td>-</td></tr><tr><td>7</td><td>2:1</td><td>14/15</td><td>2</td><td>1.</td><td>18/19</td><td>2</td><td>1.</td><td>10/10</td><td>1</td><td>1.</td><td>1</td></tr><tr><td>12</td><td>2:1</td><td>22/22</td><td>1</td><td>0.</td><td>28/34</td><td>1</td><td>0.</td><td>20/20</td><td>2</td><td>0.</td><td>2</td></tr><tr><td>40</td><td>4:3</td><td>7/8</td><td>2</td><td>5.</td><td>9/9</td><td>1</td><td>8.</td><td>7/7</td><td>1</td><td>0.</td><td>1</td></tr><tr><td>61</td><td>3:2</td><td>**</td><td>17</td><td>16.</td><td>42/53</td><td>5</td><td>1.</td><td>18/18</td><td>7</td><td>1.</td><td>2</td></tr><tr><td>104</td><td>8:4</td><td>27/28</td><td>-</td><td>3.</td><td>27/28</td><td>-</td><td>3.</td><td>25/26</td><td>-</td><td>3.</td><td>1</td></tr><tr><td>318</td><td>2:1</td><td>32/45</td><td>15</td><td>0.</td><td>*</td><td>183</td><td>0.</td><td>25/25</td><td>9</td><td>0.</td><td>1</td></tr><tr><td>319</td><td>2:1</td><td>***</td><td>20</td><td>Inf</td><td>40/66</td><td>10</td><td>1.</td><td>25/26</td><td>12</td><td>1.</td><td>1</td></tr><tr><td>320</td><td>2:1</td><td>***</td><td>38</td><td>Inf</td><td>32/40</td><td>3</td><td>1.</td><td>22/22</td><td>12</td><td>1.</td><td>2</td></tr><tr><td>321</td><td>2:1</td><td>***</td><td>157</td><td>Inf</td><td>31/37</td><td>3</td><td>2.</td><td>29/35</td><td>16</td><td>1.</td><td>4</td></tr><tr><td>322</td><td>2:1</td><td>38/47</td><td>22</td><td>4.</td><td>28/33</td><td>3</td><td>4.</td><td>23/25</td><td>10</td><td>3.</td><td>3</td></tr><tr><td>335</td><td>3:2</td><td>24/32</td><td>5</td><td>8.</td><td>25/31</td><td>1</td><td>6.</td><td>23/29</td><td>7</td><td>2.</td><td>6</td></tr><tr><td>355</td><td>4:1</td><td>35/60</td><td>3</td><td>6.</td><td>24/34</td><td>2</td><td>2.</td><td>18/21</td><td>-</td><td>2.</td><td>2</td></tr><tr><td>6</td><td>2:1</td><td>11/15</td><td>2</td><td>3.</td><td>9/10</td><td>1</td><td>2.</td><td>12/12</td><td>1</td><td>1.</td><td>2</td></tr><tr><td>10</td><td>2:1</td><td>11/11</td><td>-</td><td>1.</td><td>11/11</td><td>-</td><td>1.</td><td>11/12</td><td>-</td><td>1.</td><td>2</td></tr><tr><td>29</td><td>3:1</td><td>13/17</td><td>1</td><td>4.</td><td>9/9</td><td>1</td><td>1.</td><td>12/12</td><td>-</td><td>1.</td><td>1</td></tr><tr><td>71</td><td>4:3</td><td>6/6</td><td>-</td><td>1.</td><td>6/6</td><td>-</td><td>1.</td><td>6/6</td><td>-</td><td>2.</td><td>1</td></tr><tr><td>77</td><td>5:2</td><td>19/20</td><td>-</td><td>1.</td><td>19/20</td><td>-</td><td>1.</td><td>19/21</td><td>2</td><td>1.</td><td>2</td></tr><tr><td>80</td><td>5:3</td><td>7/7</td><td>-</td><td>1.</td><td>7/7</td><td>-</td><td>1.</td><td>7/7</td><td>1</td><td>1.</td><td>1</td></tr><tr><td>81</td><td>5:3</td><td>10/11</td><td>3</td><td>9.</td><td>9/9</td><td>2</td><td>9.</td><td>11/11</td><td>3</td><td>2.</td><td>2</td></tr><tr><td>93</td><td>6:2</td><td>24/26</td><td>1</td><td>3.</td><td>26/28</td><td>-</td><td>3.</td><td>27/29</td><td>-</td><td>3.</td><td>1</td></tr><tr><td>106</td><td>8:6</td><td>**</td><td>7</td><td>12.</td><td>*</td><td>2</td><td>11.</td><td>*</td><td>2</td><td>0.</td><td>197</td></tr><tr><td>316</td><td>2:1</td><td>33/34</td><td>2</td><td>0.</td><td>34/42</td><td>4</td><td>0.</td><td>32/36</td><td>2</td><td>0.</td><td>3</td></tr><tr><td>317</td><td>2:1</td><td>29/33</td><td>10</td><td>3.</td><td>*</td><td>183</td><td>0.</td><td>30/37</td><td>7</td><td>0.</td><td>5</td></tr><tr><td>336</td><td>3:2</td><td>31/32</td><td>1</td><td>3.</td><td>34/60</td><td>4</td><td>5.</td><td>29/36</td><td>2</td><td>4.</td><td>6</td></tr><tr><td>338</td><td>3:2</td><td>35/37</td><td>24</td><td>11.</td><td>128/334</td><td>7</td><td>9.</td><td>31/40</td><td>10</td><td>5.</td><td>7</td></tr><tr><td>375</td><td>10:9</td><td>10/11</td><td>2</td><td>4.</td><td>99/268</td><td>7</td><td>3.</td><td>15/19</td><td>6</td><td>2.</td><td>4</td></tr><tr><td>11</td><td>2:1</td><td>8/9</td><td>-</td><td>1.</td><td>8/9</td><td>-</td><td>1.</td><td>10/13</td><td>-</td><td>1.</td><td>4</td></tr><tr><td>27</td><td>3:1</td><td>19/21</td><td>4</td><td>3.</td><td>18/20</td><td>-</td><td>3.</td><td>29/36</td><td>1</td><td>3.</td><td>4</td></tr><tr><td>39</td><td>4:2</td><td>13/13</td><td>-</td><td>2.</td><td>13/13</td><td>-</td><td>2.</td><td>20/26</td><td>-</td><td>2.</td><td>4</td></tr><tr><td>65</td><td>3:1</td><td>13/17</td><td>1</td><td>2.</td><td>11/11</td><td>1</td><td>1.</td><td>15/23</td><td>-</td><td>1.</td><td>4</td></tr><tr><td>72</td><td>4:2</td><td>21/21</td><td>4</td><td>1.</td><td>21/21</td><td>-</td><td>1.</td><td>26/26</td><td>-</td><td>1.</td><td>7</td></tr><tr><td>216</td><td>2:1</td><td>11/14</td><td>-</td><td>2.</td><td>11/14</td><td>-</td><td>2.</td><td>17/26</td><td>2</td><td>0.</td><td>3</td></tr><tr><td>219</td><td>4:2</td><td>15/16</td><td>-</td><td>2.</td><td>15/16</td><td>-</td><td>2.</td><td>28/38</td><td>8</td><td>2.</td><td>5</td></tr></tbody></table>
The results are presented in Table 1. The columns in this table are labeled as follows: P is the problem number given in [17, 25], n is the number of variables, m is the number of constraints, ng is the number of gradient calculations and constraint linearizations (the main computation cost), nf is the number of function evaluations, cr is the number of corrections of $\gamma_k$ (Powell's corrections in algorithm Powell or $r_k > 0$ in BTZ and PLS), $\kappa_2$ gives the logarithm in base 10 of the $\ell_2$ condition number of the matrix $M_k$ at the last iteration, and ii is the number of extra inner iterations in the PLS (the associated simulation cost is taken into account in the counters ng and nf). A symbol '*' in the table indicates that the algorithm failed to find the solution in less than 201 linearizations. Other failures are of two kinds: either the number of backtrackings in the line-search exceeds 10 (symbol '**') or the matrix $M_k$ is so badly conditioned that $d_k$ is not a descent direction (symbol '**').

This experiment deserves some comments and allows us to draw some conclusions. Note that we do not pretend that this numerical experiment reflects the average behavior of the tested algorithms. It was mainly done to see whether our algorithm could be implemented and to compare it with other techniques on a small range of problems. Furthermore, the dimensions of the problems are very small, which prevents us from drawing firm conclusions. With this in mind, one can however quote the following points.

1. When there is no extra inner iteration in the PLS algorithm (first part of the table), PLS performs as well as Powell and BTZ. In this case, the difference between PLS and the two other techniques only lies in the computation of the vector $\gamma_k$. This indicates that our choice (3.2) of vector $\gamma_k$ is relevant.

2. In the other cases, and except the cycling observed in problem 106, on which all the algorithms fails, the number of extra inner iterations used by PLS is always small.

3. The last three parts of the table are classified according to the efficiency of PLS compared with the two other techniques. The second part of the table gathers the cases where the PLS technique reduces the overall computational cost, by having the best counters ng and nf. For a large amount of the problems (third part of the table), the performance of PLS is more or less the same as Powell or BTZ. Finally, it is worth noting that the PLS technique does not always improve the efficiency of the algorithm (last part of the table).

4. The technique for determining the augmentation parameter $r_k$ in PLS by minimizing a measure of the conditioning of the updated matrix (Section 3.5) makes their condition number smaller in PLS than in the other algorithms.

5. Finally, regarding the usual counters (ng and nf), there is no clear winner or loser, although PLS looks more robust in that it fails less often than Powell or BTZ.
# 8 Conclusion

This paper proposes a technique for maintaining the positive definiteness of the updated matrices in the quasi-Newton version of the SQP algorithm for equality constrained optimization problems. The overall algorithm generates approximations of the Hessian of the augmented Lagrangian, whose positive definiteness is obtained via the realization of a reduced Wolfe condition and an adequate setting of the augmentation parameter. The globalization is obtained by means of a nondifferentiable augmented Lagrangian function as merit function. Finite termination of the search algorithm, global convergence and admissibility of the unit step size are proved.

What is new in the proposed approach is the use of a piecewise line-search algorithm, mimicing what is done for unconstrained problems. This algorithm takes care of the positive definiteness of the matrices in the space tangent to the constraints. We believe that this brings a conceptual improvement on the algorithm proposed by Byrd, Tapia, and Zhang [4], with which our method has similarities.

The numerical experiment indicates that the proposed algorithm can be more robust than existing approaches, although there is no clear improvement on the usual counters.

Now, not all the facets of the algorithm have been studied and it is allowed to think that a precise analysis of the possibility to have superlinear convergence of the algorithm (such as in the papers [4] and [13]) may improve the algorithm. We plan to undertake this study.

# Acknowledgment

We would like to thank Richard Byrd and Jorge Nocedal for an exciting discussion at an early stage of the design of the algorithm, which has certainly influence our point of view on the implementation of the PLS technique for the SQP method.

# References

[1] P.T. Boggs and J.W. Tolle (1995). Sequential quadratic programming. *Acta Numerica 1995*, 1–51.

[2] J.F. Bonnans (1989). Asymptotic admissibility of the unit stepsize in exact penalty methods. *SIAM Journal on Control and Optimization*, 27, 631–641.

[3] R.H. Byrd and J. Nocedal (1989). A tool for the analysis of quasi-Newton methods with application to unconstrained minimization. *SIAM Journal on Numerical Analysis*, 26, 727–739.

[4] R.H. Byrd, R.A. Tapia, Y. Zhang (1992). An SQP augmented Lagrangian BFGS algorithm for constrained optimization. *SIAM Journal on Optimization*, 2, 210–241.

[5] F.H. Clarke (1983). *Optimization and Nonsmooth Analysis*. John Wiley & Sons, New York.

[6] T.F. Coleman and P.A. Fenyes (1992). Partitioned quasi-Newton methods for nonlinear constrained optimization. *Mathematical Programming Study*, 53, 17–44.
[7] J.E. Dennis and H. Wolkowicz (1993). Sizing and least-change secant methods. *SIAM Journal on Numerical Analysis*, 30, 1291–1314.

[8] R. Fletcher (1987). *Practical Methods of Optimization* (second edition). John Wiley & Sons, Chichester.

[9] D. Gabay (1982). Minimizing a differentiable function over a differential manifold. *Journal of Optimization Theory and Applications*, 37, 177–219.

[10] D. Gabay (1982). Reduced quasi-Newton methods with feasibility improvement for nonlinearly constrained optimization. *Mathematical Programming Study*, 16, 18–44.

[11] J.Ch. Gilbert (1991). Maintaining the positive definiteness of the matrices in reduced secant methods for equality constrained optimization. *Mathematical Programming*, 50, 1–28.

[12] J.Ch. Gilbert (1993). On the realization of the Wolfe conditions in reduced quasi-Newton methods for equality constrained optimization. Research Report 2127, INRIA, BP 105, 78153 Le Chesnay, France.

[13] J.Ch. Gilbert (1993). Superlinear convergence of a reduced BFGS method with piecewise line-search and update criterion. Research Report 2140, INRIA, BP 105, 78153 Le Chesnay, France.

[14] S.-P. Han (1976). Superlinearly convergent variable metric algorithms for general nonlinear programming problems. *Mathematical Programming*, 11, 263–282.

[15] S.-P. Han and O.L. Mangasarian (1979). Exact penalty functions in nonlinear programming. *Mathematical Programming*, 17, 251–269.

[16] J.-B. Hiriart-Urruty and C. Lemaréchal (1993). *Convex Analysis and Minimization Algorithms*. Grundlehren der mathematischen Wissenschaften 305 and 306, Springer-Verlag.

[17] W. Hock and K. Schittkowski (1981). *Test examples for nonlinear programming codes*. Lecture Notes in Economics and Mathematical Systems 187, Springer-Verlag, Berlin.

[18] O.L. Mangasarian (1969). *Nonlinear Programming*. McGraw-Hill, New York.

[19] J.D. Pearson (1969). Variable metric methods of minimization. *The Computer Journal*, 12, 171–178.

[20] M.J.D. Powell (1978). A fast algorithm for nonlinearly constrained optimization calculations. In G.A. Watson, ed., *Numerical Analysis Dundee 1977*, Lecture Notes in Mathematics 630 (Springer-Verlag, Berlin), 144–157.

[21] M.J.D. Powell (1978). The convergence of variable metric methods for nonlinearly constrained optimization calculations. O.L. Mangasarian, R.R. Meyer, and S.M. Robinson (eds.). *Nonlinear Programming 3*, 27–63.

[22] M.J.D. Powell (1985). The performance of two subroutines for constrained optimization on some difficult tests problems. In P.T. Boggs, R.H. Boyd and R.B. Schnabel, eds., *Numerical optimization 1984*, (SIAM, Philadelphia), 160–177.

[23] M.J.D. Powell (1991). A view of nonlinear optimization. *History of Mathematical Programming*. A Collection of Personal Reminiscences, 119–125. J.K. Lenstra, A.H.G. Rinnoy Kan, and A. Schrijver (eds). CWI North-Holland, Amsterdam.

[24] R.T. Rockafellar (1970). *Convex Analysis*. Princeton University Press, Princeton, New Jersey.
[25] K. Schittkowski (1987). *More test examples for nonlinear programming codes*. Lecture Notes in Economics and Mathematical Systems 282, Springer-Verlag, Berlin.

[26] L. Schwartz (1992). *Analyse II: Calcul Différentiel et Équations Différentielles*. Hermann, Paris.

[27] R.A. Tapia (1977). Diagonalized multiplier methods and quasi-Newton methods for constrained optimization. *Journal of Optimization Theory and Applications*, 22, 135–194.

[28] R.A. Tapia (1988). On secant updates for use in general constrained optimization. *Mathematics of Computation*, 51, 181–202.