# Global linear convergence of an augmented Lagrangian algorithm for solving convex quadratic optimization problems

Frédéric Delbos, Jean Charles Gilbert

► To cite this version:

Frédéric Delbos, Jean Charles Gilbert. Global linear convergence of an augmented Lagrangian algorithm for solving convex quadratic optimization problems. Journal of Convex Analysis, 2005, 12 (1), pp.25. inria-00071556

HAL Id: inria-00071556

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

Submitted on 23 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
INSTITUT NATIONAL DE RECHERCHE EN INFORMATIQUE ET EN AUTOMATIQUE

# Global linear convergence of an augmented Lagrangian algorithm for solving convex quadratic optimization problems

Frédéric DELBOS — J. Charles GILBERT

N° 5028

2 décembre 2003

THÈME 4

Rapport
de recherche
<figure><img src="image_2.png" /></figure>
# Global linear convergence of an augmented Lagrangian algorithm for solving convex quadratic optimization problems

Frédéric DELBOS*, J. Charles GILBERT†

Thème 4 — Simulation et optimisation
de systèmes complexes
Projet Estime

Rapport de recherche n° 5028 — 2 décembre 2003 — 27 pages

**Abstract:** We consider an augmented Lagrangian algorithm for minimizing a convex quadratic function subject to linear inequality constraints. Linear optimization is an important particular instance of this problem. We show that, provided the augmentation parameter is large enough, the constraint value converges *globally* linearly to zero. This property is proven by establishing first a global radial Lipschitz property of the reciprocal of the dual function subgradient. It is also a consequence of the proximal interpretation of the method. No strict complementarity assumption is needed. The result is illustrated by numerical experiments and algorithmic implications are discussed.

**Key-words:** Augmented Lagrangian – convex quadratic optimization – error bound – global linear convergence – linear constraints – proximal algorithm.

<footnote>* Institut Français du Pétrole, 1 & 4 avenue de Bois-Préau, 92852 Rueil-Malmaison, France; e-mail: Frederic.Delbos@ifp.fr.</footnote>

<footnote>† INRIA Rocquencourt, projet Estime, BP 105, 78153 Le Chesnay Cedex, France; e-mail: Jean-Charles.Gilbert@inria.fr.</footnote>
Convergence linéaire globale d'un algorithme de lagrangien augmenté pour la résolution d'un problème d'optimisation quadratique convexe

**Résumé** : Nous considérons un algorithme de lagrangien augmenté pour minimiser une fonction quadratique convexe sous contraintes linéaires. Les problèmes d'optimisation linéaire entrent dans ce cadre. Nous montrons que, lorsque le paramètre d'augmentation est suffisamment grand, la valeur des contraintes converge *globalement* linéairement vers zéro. Cette propriété est démontrée en établissant dans un premier temps une propriété de Lipschitz radiale globale de la fonction réciproque du sous-gradient de la fonction duale. C'est aussi une conséquence de l'interprétation proximale de l'algorithme. Cette propriété ne requière pas d'hypothèse de complémentarité stricte. Le résultat est illustré par des expériences numériques et ses implications algorithmiques sont brièvement discutées.

**Mots-clés** : Algorithme proximal – borne d’erreur – contraintes linéaires – convergence linéaire globale – lagrangien augmenté – optimisation quadratique convexe.
# 1 Introduction

Convex quadratic programs (QP) arise in their own right and as subproblems in some numerical algorithms for solving optimization problems. On the one hand, since no strictly convex assumption is made, the important class of linear optimization problems, with a zero quadratic term in their objective, enters this framework. On the other hand, the SQP algorithm decomposes a regularized constrained least-squares problem into a sequence of strictly convex QP's (see [10, 5] for an example in reflection tomography, which partly motivates this study; see [19, 2] for recent books describing the SQP algorithm). Finding efficient algorithms for solving this basic multi-faceted problem in all possible situations is an objective that has been pursued for decades (see for example the already 20 year old survey on quadratic programming in [20] and the monographs on interior point methods in [15, 7, 18, 30, 14, 29, 31]).

The method that we further explore in this paper fits into the class of dual approaches, since it is essentially the augmented Lagrangian (AL) method of Hestenes [11] and Powell [22] that is applied to a convex QP. This algorithm can be implemented in such a way that it does not require any matrix factorization. It is therefore appropriate when the problem is so large that such a factorization is impracticable or too much time consuming. This is a motivation for using the AL algorithm when the optimization problem deals with systems governed by partial differential equations [8]. In that case, however, a good preconditioner for the unavoidable CG iterations must be available.

The quadratic problem we consider in this paper is written

$$ \left\{ \begin{array}{l} \min_x \frac{1}{2} x^\top Q x + q^\top x \\ l \le C x \le u, \end{array} \right. \qquad (1.1) $$

where $Q$ is an $n \times n$ positive semi-definite matrix, $q \in \mathbb{R}^n$, $C$ is $m \times n$, and the $m$-dimensional vectors $l$ and $u$ satisfy $l < u$ and may have infinite components. With lower and upper bounds, problem (1.1) is close to what is actually implemented in numerical codes (see [6] for an example). We have not included equality constraints in (1.1) to make the presentation simple, but the results should not be altered by these constraints. Note that, since $Q$ may be zero, (1.1) also modelizes linear optimization.

The AL algorithm we study in this paper is defined on an equivalent form of (1.1) obtained by introducing an auxiliary variable $y \in \mathbb{R}^m$ [10]:

$$ \left\{ \begin{array}{l} \min_x \frac{1}{2} x^\top Q x + q^\top x \\ y = C x \\ l \le y \le u. \end{array} \right. \qquad (1.2) $$

The algorithm generates a sequence $\{\lambda_k\} \subset \mathbb{R}^m$ converging to some optimal multiplier associated with the equality constraint of (1.2). At each iteration, an auxiliary bound constrained QP has to be solved, so that the approach can be viewed as transforming (1.1) into a sequence of bound constrained convex quadratic subproblems. Two facts contribute to the possible success of this method. First, a bound constraint QP is much easier to solve than (1.1), which has general linear constraints (see [17, 9] and
the references therein). Second, because of its dual and constraint convergence, the AL
algorithm usually identifies the active constraints of (1.1) in a finite number of iterations.
Since often these constraints are also the active constraints of the subproblems close to
the solution, the combinatorial aspect of the bound constrained QP’s rapidly decreases
in intensity as the convergence progresses (and usually disappears after finitely many
AL iterations). This reasoning is valid for instance when $Q$ is positive definite and strict
complementarity holds at the solution.

The AL algorithm also generates primal iterates $(x_k, y_k) \in \mathbb{R}^n \times \mathbb{R}^m$ and is controlled
by the convergence of the constraint values to zero: if $\|y_k - Cx_k\|$ is less than a given
tolerance, optimality can be considered to be reached. The algorithm is also driven
by a so-called augmentation parameter $r_k$, whose role on the speed of convergence is
major. This paper essentially shows that, provided $r_k$ is larger than a certain positive
threshold, the convergence of the constraint norm to zero is *globally* linear, meaning that
*at each iteration* the constraint norm decreases by a factor uniformly less than one. This
property makes predictable the number of iterations to converge to a given precision
and offers a possibility to study the global iterative complexity of the algorithm.

The paper is organized as follows. In section 2, the AL algorithm under investigation
is presented with the appropriate level of details. In section 3, we give the tools from
convex analysis that are useful for the study of the method. We already expose some
of the properties of the algorithm. This section also gives a lemma of general interest,
which compares the distances from a point in the positive orthant to an affine space, on
the one hand, and to the polyhedron given by the intersection of this affine space and the
positive orthant, on the other hand. Section 4 deals with the global linear convergence
of the AL algorithm. It starts by showing a global error bound for the dual solution
set, in terms of the subgradient of the dual function. The global linear convergence is
then seen as a consequence of this property. We conclude in section 5 by relating some
numerical experiments on a seismic tomography problem and by a discussion.

Notation

We denote the Euclidean norm by $\| \cdot \|$. The distance associated with this norm is
denoted by “dist”, $B := \{ x : \|x\| \le 1 \}$ is the closed unit ball, and $\partial B := \{ x : \|x\| = 1 \}$
is the unit sphere. We note $\bar{\mathbb{R}} := \mathbb{R} \cup \{-\infty, +\infty\}$. The null space and range space of a
matrix $A$ are respectively denoted by $N(A)$ and $R(A)$. We write $A \succ 0$ [resp. $A \succ 0$]
to indicate that a symmetric matrix $A$ is positive semi-definite [resp. positive definite].
The nonnegative orthant of $\mathbb{R}^n$ is denoted by $\mathbb{R}_+^n := \{ x \in \mathbb{R}^n : x \ge 0 \}$.

Let $E$ be a finite dimensional Euclidean space. The indicator function of a set $S \subset E$ is denoted by $\mathcal{I}_S$ (this is the function that vanishes on $S$ and takes the value $+\infty$ outside $S$). The domain of a function $f : E \to \mathbb{R} \cup \{+\infty\}$ is defined by $\text{dom } f := \{x \in E : f(x) < +\infty\}$ and its epigraph by $\text{epi } f := \{(x, \alpha) \in E \times \mathbb{R} : f(x) \le \alpha\}$. As in [12], $\overline{\text{conv}}(E)$ is the set of functions $f : E \to \mathbb{R} \cup \{+\infty\}$ that are convex (epi $f$ is convex), proper (epi $f \ne \emptyset$), and closed (epi $f$ is closed). The subdifferential at $x \in E$ of a function $f \in \overline{\text{conv}}(E)$ is denoted by $\partial f(x)$. We denote by $N_C(x)$ the normal cone
at $x$ to a convex set $C \subset E$. The orthogonal projection of a point $x$ onto a nonempty closed convex set $C$ is denoted by $P_C(x)$.

## 2 An AL algorithm for solving the QP

We assume throughout that problem (1.2) has a solution and denote by $S_P$ the set of its solutions $(\bar{x}, \bar{y})$. The projections of $S_P$ onto $\mathbb{R}^n$ and $\mathbb{R}^m$ are respectively denoted by $S_P^x := \{\bar{x} \in \mathbb{R}^n : (\bar{x}, \bar{y}) \in S_P$ for some $\bar{y} \in \mathbb{R}^m\}$ and $S_P^y := \{\bar{y} \in \mathbb{R}^m : (\bar{x}, \bar{y}) \in S_P$ for some $\bar{x} \in \mathbb{R}^m\}$. Since the constraints of (1.2) are qualified, there exist optimal multipliers, which certainly implies that the affine subspace

$$ \Lambda := \{\lambda \in \mathbb{R}^m : C^\top \lambda \in q + R(Q)\} \quad (2.1) $$

is nonempty. Note that $\Lambda = \mathbb{R}^m$ if $Q > 0$.

The augmented Lagrangian is obtained by dualizing the equality constraint of (1.2). It is the function $\ell_r : (x, y, \lambda) \in \mathbb{R}^n \times \mathbb{R}^m \times \mathbb{R}^m \to \mathbb{R}$, defined by

$$ \ell_r(x, y, \lambda) = \frac{1}{2} x^\top Q x + q^\top x + \lambda^\top (y - Cx) + \frac{r}{2} \|y - Cx\|^2, \quad (2.2) $$

where $r \ge 0$ is called the augmentation parameter (see [1, 11, 22]).

We can now give a precise statement of the AL algorithm we study in this paper, which is basically the method of Hestenes [11] and Powell [22] applied to (1.2).

AL ALGORITHM for solving (1.1)

Initialization: choose $\lambda_0 \in \mathbb{R}^m$ and $r_0 > 0$.

Repeat for $k = 0, 1, 2, \dots$

1. Solve:

$$ \min_{x \atop l \le y \le u} \ell_{r_k}(x, y, \lambda_k). \quad (2.3) $$

Denote a solution by $(x_{k+1}, y_{k+1})$.

2. Update the multiplier

$$ \lambda_{k+1} = \lambda_k + r_k(y_{k+1} - Cx_{k+1}). \quad (2.4) $$

3. Stop if $y_{k+1} \approx Cx_{k+1}$.

4. Choose a new augmentation parameter: $r_{k+1} > 0$.

This algorithm deserves some comments.
1. Under the sole assumption that problem (1.1) has a solution, the QP in step (2.3) has also a solution. This fact is clarified in proposition 3.3. This solution is not necessary unique however.

2. Even though $(x_{k+1}, y_{k+1})$ is not uniquely determined as a solution to (2.3), $y_{k+1} - Cx_{k+1}$ is independent of that solution, so that the multipliers $\lambda_k$ are unambiguously determined.

3. The augmentation parameter $r_k$ can change from iteration to iteration, but the same value must be used in the AL minimized in step 1 and in the multiplier update in step 2. If the “step-size” in (2.4) is different from $r_k$ (with the aim at minimizing better the dual function, as in [21, section 4.2] for example), several properties of the AL algorithm may no longer hold, such as the finite identification of active constraints (in the presence of strict complementarity) and the global linear convergence of section 4.

4. The larger are the augmentation parameters $r_k$, the faster is the convergence. The only limitation on a large value for $r_k$ comes from the ill-conditioning that such a value induces in the AL and the resulting difficulty in solving (2.3). Actually, it is clear from the structure of the AL in (2.2) that a large $r$ gives priority to the restoration of the equality constraint, leaving aside the minimization of the Lagrangian (whose role is to provide optimality).

In comparison with an interior point method, which faces the combinatorial aspect of (1.1) by transforming the problem into a sequence of linear systems, the AL algorithm goes around this difficulty by transforming a general QP into a sequence of QP's with simple bounds, which are easier to solve. Indeed, a number of efficient algorithms are available for dealing with the bound constraints on the AL in step 1. A possibility would be to minimize first analytically $\ell_r$ in $y$ and then to minimize the resulting function in $x$. Unfortunately this function of $x$, which is the AL associated with the inequality constrained QP (1.1) [24, 26], has a combinatorial structure (it contains maxima) that is not easier to deal with than the direct numerical minimization of $\ell_r$ in $(x, y)$. In our code QPAL [6], used for the numerical experiments of section 5 and in [5], we have adapted to the $(x, y)$ structure of problem (1.2) an active set strategy together with the gradient projection algorithm and conjugate gradient iterations on the activated faces (see [17, 9] and the references therein).

As opposed to standard (non shifted) interior point methods, whose elementary linear systems have an exploding condition number, the AL algorithm does not require the penalty parameter $r_k$ to go to infinity. Actually, any sequence $\{r_k\}$ that remains bounded away from zero guarantees the convergence, even though large values speed it up [28]. Therefore, the bound constrained QP's in step 1 can be maintained reasonably well conditioned, keeping satisfactory the efficiency of a conjugate gradient based solver. This remark reinforces the viewpoint that considers the AL algorithm as a method suitable for large problems.
# 3 Convex analysis tools

## 3.1 Duality

As dual function associated with problem (1.2), we use the one obtained by dualizing its equality constraints. It is the function $\delta: \mathbb{R}^m \to \mathbb{R} \cup \{+\infty\}$ defined by

$$ \lambda \mapsto \delta(\lambda) := - \inf_{y \in [l, u]} \left( \frac{1}{2} x^\top Q x + q^\top x + \lambda^\top (y - Cx) \right). \quad (3.1) $$

Clearly $\delta \in \text{Conv}(\mathbb{R}^m)$ (it takes a finite value, for instance, when $\lambda$ is an optimal multiplier associated with the equality constraint of (1.2)).

For a given $\lambda \in \mathbb{R}^m$, $(x_\lambda, y_\lambda)$ is a solution to the Lagrange problem, the minimization problem in (3.1), if and only if $x_\lambda \in X_\lambda$ and $y_\lambda \in Y_\lambda$, where $X_\lambda$ is the affine space

$$ X_{\lambda} := \{x \in \mathbb{R}^n : Qx = C^{\top}\lambda - q\} \quad (3.2) $$

and $Y_\lambda$ is the Cartesian product of the following intervals

$$ (Y_{\lambda})_i = \begin{cases} [u_i, u_i] & \text{if } \lambda_i < 0 \\ [l_i, u_i] & \text{if } \lambda_i = 0 \\ [l_i, l_i] & \text{if } \lambda_i > 0. \end{cases} \quad (3.3) $$

These intervals, with their possible infinite bounds, have to be understood in a broad sense: for example, $[l_i, u_i]$ is the interval $]-\infty, u_i]$ if $l_i = -\infty$ and $u_i$ is finite, $[l_i, l_i]$ is the empty set if $l_i = -\infty$, etc. Since the Lagrange problem is always feasible and since a feasible convex quadratic problem has a solution if and only if it is bounded (see [2, theorem 17.1] for example), the domain of $\delta$ is the set of $\lambda$'s for which the Lagrange problem has a solution. Therefore $\text{dom} \delta$ can be written as the nonempty polyhedron

$$ \text{dom} \delta = \{\lambda \in \mathbb{R}^m : X_\lambda \neq \emptyset, Y_\lambda \neq \emptyset\} = \mathbb{R}^m_{l,u} \cap \Lambda, $$

where

$$ \mathbb{R}^m_{l,u} := \{\lambda \in \mathbb{R}^m : \lambda_i \le 0 \text{ if } l_i = -\infty, \lambda_i \ge 0 \text{ if } u_i = +\infty\}. $$

Observe finally that the multivalued function $\lambda \mapsto -Y_\lambda$ is monotone: for $\lambda$ and $\lambda' \in \mathbb{R}^m$, and for $y_\lambda \in Y_\lambda$ and $y_{\lambda'} \in Y_{\lambda'}$, there holds

$$ -(y_{\lambda'} - y_{\lambda})^{\top}(\lambda' - \lambda) \ge 0. \quad (3.4) $$

Let $Q^\dagger$ be the pseudo-inverse of $Q$ and take the notation

$$ H := CQ^\dagger C^\top \quad \text{and} \quad v := CQ^\dagger q. \quad (3.5) $$
Let $\lambda \in \text{dom } \delta$. Minimizing explicitly with respect to $x$ in (3.1) yields

$$ 
\begin{aligned}
\delta(\lambda) &= \sup_{y \in [l, u]} \left( \frac{1}{2} \lambda^\top H \lambda - (v+y)^\top \lambda + \frac{1}{2} q^\top Q^\dagger q \right) \\
&= \sup_{\substack{y=(y_i)_{i=1}^m \\ y_i = u_i \text{ if } \lambda_i < 0 \\ y_i = l_i \text{ if } \lambda_i > 0}} \left( \frac{1}{2} \lambda^\top H \lambda - (v+y)^\top \lambda + \frac{1}{2} q^\top Q^\dagger q \right).
\end{aligned} 
\quad (3.6) $$

On its domain, the dual function $\delta$ is therefore the maximum of a finite number of convex quadratic functions, which only differ by their slope at the origin (in particular, they have the same Hessian $H$).

**Lemma 3.1** The subdifferential of the dual function (3.1) is given at $\lambda \in \text{dom } \delta$ by

$$ \partial \delta(\lambda) = H\lambda - v - Y_{\lambda} + C(N(Q)). $$

**PROOF.** We write $\delta$ as the sum of three convex functions. Let $\tilde{l}$ and $\tilde{u} \in \mathbb{R}^m$ be chosen such that $\tilde{l} < \tilde{u}$, $\tilde{l}_i = l_i$ if $l_i$ is finite, and $\tilde{u}_i = u_i$ if $u_i$ is finite. Define $\tilde{Y}_\lambda$ by formula (3.3) with $l$ and $u$ respectively replaced by $\tilde{l}$ and $\tilde{u}$. Then the following finite value function $\tilde{\delta} \in \text{Conv}(\mathbb{R}^m)$ is identical to the right hand side of (3.6) on $\mathbb{R}_{l,u}^m$:

$$ \tilde{\delta}(\lambda) = \sup_{\substack{y=(y_i)_{i=1}^m \\ y_i = \tilde{l}_i \text{ or } \tilde{u}_i}} \left( \frac{1}{2} \lambda^\top H \lambda - (v+y)^\top \lambda + \frac{1}{2} q^\top Q^\dagger q \right). $$

Clearly $\delta = \tilde{\delta} + \mathcal{I}_{\mathbb{R}_{l,u}^m} + \mathcal{I}_\Lambda$, so that theorem 23.8 in [23] implies that for $\lambda \in \text{dom } \delta$:

$$ \partial \tilde{\delta}(\lambda) = \partial \tilde{\delta}(\lambda) + \partial \mathcal{I}_{\mathbb{R}_{l,u}^m}(\lambda) + \partial \mathcal{I}_\Lambda(\lambda). $$

Equality holds above because $\mathcal{I}_{\mathbb{R}_{l,u}^m}$ and $\mathcal{I}_\Lambda$ are polyhedral and because $\text{dom } \tilde{\delta} (= \mathbb{R}^m$, $\text{ri}$ denotes the relative interior), $\text{dom } \mathcal{I}_{\mathbb{R}_{l,u}^m} = \mathbb{R}_{l,u}^m$, and $\text{dom } \mathcal{I}_\Lambda = \Lambda$ have a point in common (one in $\text{dom } \delta \neq \emptyset$).

To compute $\partial \tilde{\delta}(\lambda)$, we use corollary VI.4.3.2 in [12]:

$$ \partial \tilde{\delta}(\lambda) = \text{conv} \left\{ H\lambda - v - y : y \in \tilde{Y}_{\lambda} \text{ and } (y_i = \tilde{l}_i \text{ or } \tilde{u}_i \text{ if } \lambda_i = 0) \right\} = H\lambda - v - \tilde{Y}_{\lambda}. $$

On the other hand, $\partial \mathcal{I}_{\mathbb{R}_{l,u}^m}(\lambda) = N_{\mathbb{R}_{l,u}^m}(\lambda)$, which is the set of vectors $\nu \in \mathbb{R}^m$ satisfying

$$ \left\{ 
\begin{array}{ll}
\nu_i \ge 0 & \text{when } \lambda_i = 0, l_i = -\infty, \text{ and } u_i \text{ is finite} \\
\nu_i \le 0 & \text{when } \lambda_i = 0, l_i \text{ is finite, and } u_i = +\infty \\
\nu_i \in \mathbb{R} & \text{when } \lambda_i = 0, l_i = -\infty, \text{ and } u_i = +\infty \\
\nu_i = 0 & \text{when } \lambda_i \ne 0.
\end{array}
\right. $$

We deduce from this computation that for $\lambda \in \text{dom } \delta$

$$ \tilde{Y}_{\lambda} - \partial \mathcal{I}_{\mathbb{R}_{l,u}^m}(\lambda) = Y_{\lambda}. $$
Finally $\partial \mathcal{I}_{\Lambda}(\lambda) = N_{\Lambda}(\lambda) = \{\mu \in \mathbb{R}^m : C^{\top}\mu \in R(Q)\}^{\perp} = C(N(Q))$. Adding the last three subdifferentials provides the formula of $\partial \delta(\lambda)$ given in the statement of the lemma. $\square$

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

Let us denote by $S_D$ the set of dual solutions:

$$S_D := \{\bar{\lambda} \in \mathbb{R}^m : 0 \in \partial\delta(\bar{\lambda})\}.$$

Not surprisingly, this is a convex polyhedron, which can be described in the standard form. It will be useful to make this form explicit and we do so in lemma 3.2 below. For this, we take a partition of $\{1, \dots, m\}$ into the index sets

$$ 
\begin{align} 
I_l &:= \{i : \bar{y}_i = l_i \text{ for all } (\bar{x}, \bar{y}) \in S_P\}, \nonumber \\
J &:= \{i : l_i < \bar{y}_i < u_i \text{ for some } (\bar{x}, \bar{y}) \in S_P\}, \tag{3.7} \\
I_u &:= \{i : \bar{y}_i = u_i \text{ for all } (\bar{x}, \bar{y}) \in S_P\}. \nonumber
\end{align} 
$$

We also introduce the orthant face $\mathcal{O}$ and the affine subspace $\mathcal{A}$

$$ 
\begin{align} 
\mathcal{O} &:= \{\lambda \in \mathbb{R}^m : \lambda_{I_l} \ge 0, \lambda_J = 0, \lambda_{I_u} \le 0\}, \nonumber \\
\mathcal{A} &:= \{\lambda \in \mathbb{R}^m : C^{\top}\lambda = Q\bar{x} + q\}. \tag{3.8}
\end{align} 
$$

In the definition of $\mathcal{A}$, $\bar{x}$ is an arbitrary primal solution. We have not made this dependence explicit in the symbol of the set since, as shown in the next lemma, $\mathcal{A}$ does not depend on the choice of $\bar{x} \in S_P^x$.

**Lemma 3.2** The set of dual solutions $S_D$ can be written as the intersection

$$S_D = \mathcal{O} \cap \mathcal{A}.$$

Furthermore, for any $\bar{\lambda} \in S_D$ and any $\bar{y} \in S_P^y$, we have $\bar{y} \in Y_{\bar{\lambda}}$ and $H\bar{\lambda} = v + \bar{y} + C\bar{u}$ for some $\bar{u} \in N(Q)$.

**PROOF.** Let $\ell(x, y, \lambda) = \frac{1}{2}x^\top Q x + q^\top x + \mathcal{I}_{[l,u]}(y) + \lambda^\top(y - Cx)$ be the Lagrangian function of the problem $\min_{(x,y)} \{\frac{1}{2}x^\top Q x + q^\top x + \mathcal{I}_{[l,u]}(y) : y = Cx\}$, which is equivalent to (1.2). Since the constraint of this problem is qualified, $\bar{\lambda} \in S_D$ if and only if $0 \in \partial_{(x,y)}\ell(\bar{x}, \bar{y}, \bar{\lambda})$, where $(\bar{x}, \bar{y})$ is an arbitrary primal solution. This can also be written $Q\bar{x} + q = C^{\top}\bar{\lambda}$ and $0 \in N_{[l,u]}(\bar{y}) + \bar{\lambda}$, which is equivalent to $\bar{\lambda} \in \mathcal{A}_{\bar{x}} \cap \mathcal{O}_{\bar{y}}$, where

$$ 
\begin{align*} 
\mathcal{A}_{\bar{x}} &:= \{\lambda \in \mathbb{R}^m : C^{\top}\lambda = Q\bar{x} + q\}, \\ 
\mathcal{O}_{\bar{y}} &:= \{\lambda \in \mathbb{R}^m : \lambda_i \ge 0 \text{ if } \bar{y}_i = l_i, \lambda_i = 0 \text{ if } l_i < \bar{y}_i < u_i, \lambda_i \le 0 \text{ if } \bar{y}_i = u_i\}. 
\end{align*} 
$$

By varying $(\bar{x}, \bar{y}) \in S_P$, we see that $\mathcal{A}_{\bar{x}} = \mathcal{A}$ is independent of the chosen primal solution $\bar{x} \in S_P^x$ and that $\bar{\lambda} \in \mathcal{O}_{\bar{y}} \cap \mathcal{O}_{\bar{y}}$.

For proving the second part of the lemma, take $\bar{\lambda} \in S_D$ and $(\bar{x}, \bar{y}) \in S_P$. We have shown that $\bar{\lambda} \in \mathcal{A}_{\bar{x}} \cap \mathcal{O}_{\bar{y}}$. Actually, $\bar{\lambda} \in \mathcal{O}_{\bar{y}}$ is equivalent to $\bar{y} \in Y_{\bar{\lambda}}$. By $\bar{\lambda} \in \mathcal{A}_{\bar{x}}$, we have that $C^{\top}\bar{\lambda} = Q\bar{x} + q$. Multiplying to the left both sides of this equation by $CQ^\dagger$ provides $H\bar{\lambda} = v + \bar{y} + C\bar{u}$, where $\bar{u} := (Q^\dagger Q - I)\bar{x} \in N(Q)$. $\square$

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

The fact observed in the proof above that the gradient of the criterion of the primal problem at a solution, here $Q\bar{x} + q$, is independent of the chosen solution is a property of general convex problems; see [16, 3].
## 3.2 Proximality

We will use the fundamental result of Rockafellar [25], according to which the AL algorithm of section 2 is the proximal algorithm on the dual function $\delta$. More precisely, the multiplier $\lambda_{k+1}$ computed in step 2 of the AL algorithm is also the unique solution to

$$ \inf_{\lambda \in \mathbb{R}^m} \left( \delta(\lambda) + \frac{1}{2r_k} \|\lambda - \lambda_k\|^2 \right). \qquad (3.9) $$

The same parameter $r_k > 0$ is used above and in (2.3). In addition, the optimal value of this problem is the opposite of the optimal value of problem (2.3). The optimality conditions of problem (3.9) can be written $0 \in \partial\delta(\lambda_{k+1}) + (\lambda_{k+1} - \lambda_k)/r_k$. Using (2.4), we see that:

$$ Cx_k - y_k \in \partial\delta(\lambda_k), \quad \forall k \ge 1. \qquad (3.10) $$

Note that, since $\lambda_{k+1}$ is uniquely determined as the solution to (3.9), this is also the case for $y_{k+1} - Cx_{k+1}$, even though $x_{k+1}$ and $y_{k+1}$ are not uniquely determined.

Let us now clarify the conditions ensuring that the augmented Lagrange problem (2.3) has a solution.

**Proposition 3.3** *The following three properties are equivalent:*

(i) $\text{dom } \delta \neq \emptyset$,

(ii) problem (1.1), with a possible finite shift of its finite bounds to make it feasible, has a solution,

(iii) for (some or any) $r_k > 0$ and $\lambda_k \in \mathbb{R}^m$, problem (2.3) has a solution.

**PROOF.** [(i) ⇒ (iii)] Fix $r_k > 0$ and $\lambda_k \in \mathbb{R}^m$ (not necessarily the $k$th iterate). Since $\text{dom } \delta \neq \emptyset$, the optimal value of (3.9) is finite, so that the optimal value of problem (2.3) is also finite. As a feasible bounded quadratic problem, (2.3) must have a solution [2, theorem 17.1].

[(iii) ⇒ (ii)] Let $\tilde{l}$ and $\tilde{u} \in \bar{\mathbb{R}}^m$ be such that $\tilde{l} < 0 < \tilde{u}$, $\tilde{l}_i = -\infty$ iff $l_i = -\infty$ and $\tilde{u}_i = +\infty$ iff $u_i = +\infty$ (these new bounds result from a finite shift of the finite bounds of (1.1) that makes this problem feasible) and assume that the problem $\min\{f(x) : \tilde{l} \le Cx \le \tilde{u}\}$, where $f(x) := (1/2)x^\top Q x + q^\top x$, has no solution. Then, there exists a sequence $\{x_j\}$ such that $Cx_j \in [\tilde{l}, \tilde{u}]$ and $f(x_j) \to -\infty$ when $j \to \infty$ (since a bounded feasible quadratic problem has a solution). Let $y_j := P_{[l,u]}(Cx_j)$ be the projection of $Cx_j$ on $[l, u]$. Then $\|y_j - Cx_j\| \le m^{1/2}\|y_j - Cx_j\|_{\infty} \le m^{1/2}\gamma$, where $\gamma := \max(\|\tilde{l} - l\|_{\infty}, \|\tilde{u} - u\|_{\infty})$ (these norms are taken on the finite components of $l$ and $u$), and $(x_j, y_j)$ is feasible for problem (2.3). On the other hand, for an arbitrary $r_k > 0$ and $\lambda_k \in \mathbb{R}^m$, $l_{r_k}(x_j, y_j, \lambda_k) \le f(x_j) + m^{1/2}\gamma\|\lambda_k\| + (r_k/2)m\gamma^2 \to -\infty$ when $j \to \infty$. Therefore problem (2.3) has no solution.

[(ii) ⇒ (i)] Introduce $\tilde{l}$, $\tilde{u}$, and $f$ as in the previous paragraph. By assumption, the problem $\min\{f(x) : \tilde{l} \le Cx \le \tilde{u}\}$ has a solution, $\tilde{x}$ say. Since its constraints are qualified, there exist $\tilde{\lambda}^l$ and $\tilde{\lambda}^u$ such that $Q\tilde{x} + q = C^\top(\tilde{\lambda}^l - \tilde{\lambda}^u)$, $\tilde{\lambda}^l \ge 0$, $\tilde{\lambda}^u \ge 0$, $\tilde{\lambda}_i^l = 0$
if $l_i = -\infty$, and $\tilde{\lambda}_i^u = 0$ if $u_i = +\infty$. It is easy to check that $\tilde{\lambda}^l - \tilde{\lambda}^u \in \mathbb{R}_{l,u}^m \cap \Lambda = \text{dom} \delta$. $\square$

If the original quadratic problem (1.1) has a solution, condition (ii) above holds (without having to shift the bounds), so that the augmented Lagrange problem (2.3) has a solution.

### 3.3 Projection onto a convex polyhedron

This section gives two lemmas related to the projection onto a convex polyhedron. The first lemma has a general interest. It compares the distance from a point $x$ in the positive orthant to a convex polyhedron $\mathcal{X}$ defined in the standard form and the distance from $x$ to the underlying affine space $\mathcal{A}$. It is claimed that the second distance is bounded below by a positive constant (independent of $x \ge 0$) times the first one. Of course, since $\mathcal{X} \subset \mathcal{A}$, $\text{dist}(x, \mathcal{A}) \le \text{dist}(x, \mathcal{X})$.

**Lemma 3.4** Let $A$ be an $m \times n$ matrix and $b \in \mathbb{R}^m$. Consider the affine subspace $\mathcal{A}$ and the convex polyhedron $\mathcal{X}$ defined by

$$ \mathcal{A} := \{x \in \mathbb{R}^n : Ax = b\} \quad \text{and} \quad \mathcal{X} := \{x \in \mathbb{R}^n : Ax = b, x \ge 0\}. $$

These sets are supposed to be nonempty. Then, there exists a constant $\gamma > 0$ such that

$$ \forall x \in \mathbb{R}_+^n, \quad \text{dist}(x, \mathcal{A}) \ge \gamma \text{dist}(x, \mathcal{X}). $$

**PROOF.** First stage: reformulation of the statement of the lemma. By using the triangle inequality, it is easy to see that the conclusion of the lemma is equivalent to claim the existence of a constant $\gamma' > 0$ such that for all $x \in \mathbb{R}_+^n$:

$$ \|x - P_{\mathcal{A}}(x)\| \ge \gamma' \|P_{\mathcal{A}}(x) - P_{\mathcal{X}}(x)\|, \quad (3.11) $$

This inequality suggests that a certain function (whose value is the left hand side of (3.11) divided by the factor of $\gamma'$ in the right hand side) has a positive slope. This is the strategy we follow to establish (3.11).

Let $x \in \mathbb{R}_+^n$ and let us simplify the notation by introducing $x^0 := P_{\mathcal{X}}(x)$ and $x^1 := P_{\mathcal{A}}(x)$. Observe that $x - x^1 = A^\top y$ for some $y \in \mathbb{R}^m$. On the other hand, one can suppose that $x^1 \ne x^0$, since otherwise (3.11) is trivially satisfied. Observe then that $x^1 - x^0$ is a nonzero direction, normal to $\mathcal{X}$ at $x^0$, which is in the null space of $\mathcal{A}$. As a result, an arbitrary point $x \in \mathbb{R}_+^n$ such that $P_{\mathcal{A}}(x) \ne P_{\mathcal{X}}(x)$, can be found from a point $x^0 \in \mathcal{X}$ having a unitary normal direction $d$ in the null space of $\mathcal{A}$, by adding a positive displacement along $d$, leading to $x_\alpha^1 := x^0 + \alpha d \in \mathcal{A}$ ($\alpha > 0$), and finally by adding a displacement $A^\top y$, normal to $N(\mathcal{A})$, such that $x = x_\alpha^1 + A^\top y \in \mathbb{R}_+^n$. For the points $x \in \mathbb{R}_+^n$, described in this manner, the inequality (3.11) can be written

$$ \|A^\top y\| \ge \gamma' \alpha. $$
For fixed $x^0$, $d$, and $\alpha$, it is sufficient to consider the points $x \ge 0$ giving the smallest value to $\|A^\top y\|$. This motivates the introduction of the function (its dependence on $x^0$ and $d$ is not mentioned to keep the notation light)

$$ \varphi : \alpha \mapsto \varphi(\alpha) := \inf_y \{ \|A^\top y\| : x^0 + \alpha d + A^\top y \ge 0 \}. \quad (3.12) $$

Therefore, we now have to show that there exists a constant $\gamma' > 0$, such that for all $x^0 \in \mathcal{X}$, all $d \in N_{\mathcal{X}}(x^0) \cap N(A) \cap \partial B$, and all $\alpha > 0$, there holds

$$ \varphi(\alpha) \ge \gamma'\alpha. $$

Observe that $\varphi(\alpha) \ge 0$, that $\varphi(0) = 0$, and that $\varphi(t\alpha) \le t\varphi(\alpha)$ when $\alpha \ge 0$ and $t \in ]0, 1]$. This last property follows from the fact that, when $\varphi(\alpha) < \infty$, the minimization problem in (3.12) has a solution, say $y_\alpha$; adding $(1-t)(x^0 + A^\top 0) \ge 0$ and $t(x_\alpha^1 + A^\top y_\alpha) \ge 0$ provides $x_{t\alpha}^1 + A^\top (ty_\alpha) \ge 0$, so that $\varphi(t\alpha) \le \|A^\top (ty_\alpha)\| = t\varphi(\alpha)$. This property implies that $\alpha \in \mathbb{R}_+ \mapsto \varphi(\alpha)/\alpha$ is nondecreasing. Therefore, we have reduced the problem to show that the right derivative of $\varphi$ at zero satisfies

$$ \varphi'(0; 1) \ge \gamma'. \quad (3.13) $$

where $\gamma' > 0$ is a constant independent of $x^0 \in \mathcal{X}$ and $d \in N_{\mathcal{X}}(x^0) \cap N(A) \cap \partial B$.

**Second stage: control of the decomposition of the normal directions.** Consider a point $x^0 \in \mathcal{X}$ having a unitary normal direction in the null space of $A$, say $d \in N_{x^0}(\mathcal{X}) \cap N(A) \cap \partial B$. Define $I := I(x^0) := \{i : x_i^0 = 0\}$ and $J := J(x^0) := \{i : x_i^0 > 0\}$. These directions $d$ are characterized by the conditions

$$ d = A^\top z - r, \quad r_I \ge 0, \quad r_J = 0, \quad Ad = 0, \quad \text{and} \quad \|d\| = 1, \quad (3.14) $$

for some vectors $z \in \mathbb{R}^m$ and $r \in \mathbb{R}^n$. The decomposition of $d$ in $A^\top z - r$ as above is not necessarily unique. It will be useful to identify a decomposition that provides the smallest value to $\|A^\top z\|$, which is therefore a solution to

$$ \left\{ \begin{array}{l} \min_{(z,r)} \frac{1}{2} \|A^\top z\|^2 \\ A^\top z - r = d \\ r_I \ge 0 \\ r_J = 0. \end{array} \right. $$

It is easy to show that this problem has a solution, which is characterized by (3.14) and:

$$ A(A^\top z - s) = 0, \quad s_I \ge 0, \quad \text{and} \quad s_I^\top r_I = 0, \quad (3.15) $$

for some vector $s \in \mathbb{R}^n$.

Let us show that

$$ \begin{array}{lll} \max_{x^0 \in \mathcal{X}} & \sup_{d \in N_{\mathcal{X}}(x^0)} & \min_{(z,r) \in \mathbb{R}^m \times \mathbb{R}^n} \\ & Ad=0 & A^\top z - r = d \\ & \|d\| = 1 & r_{I(x^0)} \ge 0 \\ & & r_{J(x^0)} = 0 \end{array} \quad \|A^\top z\| < +\infty. \quad (3.16) $$
We see on (3.14) that two points $x^0 \in \mathcal{X}$ having the same index set $I$ have the same normal cone. Therefore, the point $x^0 \in \mathcal{X}$ intervenes in (3.16) only through its index sets $I$ and $J$. Since there is a finite number of such sets, one can fix $x^0$, hence $I$ and $J$. Let us continue by contradiction, assuming that there exists a sequence $\{(d^k, z^k, r^k, s^k)\}$ such that $Ad^k = 0$, $\|d^k\| = 1$, $A^\top z^k - r^k = d^k$, $r_I^k \ge 0$, $r_J^k = 0$, $A(A^\top z^k - s^k) = 0$, $s_I^k \ge 0$, $(s_I^k)^\top r_I^k = 0$, and $\|A^\top z^k\| \to \infty$. Extracting a subsequence if necessary, it can be assumed that $A^\top z^k / \|A^\top z^k\| \to A^\top \bar{z}$, a vector of unit norm. Since $\{d^k\}$ is bounded, the identity $A^\top z^k - r^k = d^k$ shows that $r^k / \|A^\top z^k\|$ converges to $\bar{r} := A^\top \bar{z}$. Multiplying the identity $A(A^\top z^k - s^k) = 0$ by $\bar{z}$, one finds for sufficiently large $k$

$$0 = \bar{z}^\top A A^\top z^k - \bar{r}^\top s^k = \bar{z}^\top A A^\top z^k,$$

because, when $\bar{r}_i > 0$, then $i \in I$ and, for all sufficiently large $k$, $r_I^k > 0$, so that $s_I^k = 0$. Dividing the right hand side by $\|A^\top z^k\|$ and taking the limit, one would find $A^\top \bar{z} = 0$, which provides the expected contradiction.

*Third stage: lower bound for $\varphi'(0;1)$ and conclusion.* Let us introduce $v : \mathbb{R}^n \to \mathbb{R} \cup \{+\infty\}$, the value function of the problem

$$\left\{ \begin{array}{l} \inf_y \|A^\top y\| \\ x^0 + A^\top y \ge 0, \end{array} \right. \qquad (3.17)$$

which is the proper convex function defined by $v(p) := \inf\{\|A^\top y\| : x^0 + A^\top y \ge p\}$. Then, for fixed $x^0 \in \mathcal{X}$ and $d \in N_\mathcal{X}(x^0) \cap N(A) \cap \partial B$, $\varphi(\alpha)$ defined by (3.12) can be written $\varphi(\alpha) = v(-\alpha d)$. Therefore

$$\varphi'(0; 1) = v'(0; -d) \ge -g^\top d, \quad \forall g \in \partial v(0). \qquad (3.18)$$

As for the subdifferential $\partial v(0)$, it is formed of the optimal multipliers associated with the constraint of (3.17), which are the $g$-parts of the pairs $(g, u)$ satisfying

$$g \in u + N(A), \quad \|u\| \le 1, \quad g \ge 0, \quad \text{and} \quad (x^0)^\top g = 0. \qquad (3.19)$$

Let $d = A^\top z - r$ be a decomposition of $d$ satisfying (3.14)-(3.15). If $A^\top z = 0$, then $g := -\alpha d = \alpha r$ is a subgradient of $v$ at zero for any $\alpha \ge 0$ (the conditions (3.19) are satisfied with $u = 0$; recall that $d \in N(A)$), so that (3.18) shows that $\varphi'(0; 1) = +\infty$. If $A^\top z \ne 0$, then $g := (A^\top z - d) / \|A^\top z\| = r / \|A^\top z\|$ is a subgradient of $v$ at zero (the conditions (3.19) are satisfied with $u = A^\top z / \|A^\top z\|$). Then (3.18) shows that $\varphi'(0; 1) \ge 1 / \|A^\top z\|$. Since for these decompositions, stage 2 of the proof has shown that $A^\top z$ is bounded, (3.13) holds and, consequently, the result is proven. $\square$

As shown by the following example, lemma 3.4 no longer holds with all its generality when $\mathcal{X}$ is the intersection of an affine space $\mathcal{A}$ and an arbitrary closed convex cone.
**Example 3.5** Let us introduce the following closed convex cone $K := \{x \in \mathbb{R}^3 : x_2x_3 \ge x_1^2, x_2 \ge 0, x_3 \ge 0\}$, the $1 \times 3$ matrix $A := (0 \ 1 \ 0)$, and $b = 0 \in \mathbb{R}$. Define the affine space $\mathcal{A}$ and its intersection with $K$ by

$$ \mathcal{A} := \{x \in \mathbb{R}^3 : Ax = b\} = \{x \in \mathbb{R}^3 : x_2 = 0\}, \\ \mathcal{X} := K \cap \mathcal{A} = \{x \in \mathbb{R}^3 : x_1 = x_2 = 0, x_3 \ge 0\}. $$

Then the conclusion of Lemma 3.4 does not hold for these sets $\mathcal{A}$ and $\mathcal{X}$. To see this, fix $x_1 > 0$ and consider the points $x^t := (x_1, x_1^2/t, t)$ for $t \uparrow +\infty$. Clearly $x^t \in K$, $P_{\mathcal{A}}(x^t) = (x_1, 0, t)$, and $P_{\mathcal{X}}(x^t) = (0, 0, t)$. Therefore, $\|x^t - P_{\mathcal{A}}(x^t)\|/\|P_{\mathcal{A}}(x^t) - P_{\mathcal{X}}(x^t)\| = x_1/t$ is not bounded away from zero. $\square$

Actually, it will be useful below to have the following relaxed version of lemma 3.4. This one allows the projected point $x$ not to belong to $\mathbb{R}_+^n$. This point must however be sufficiently close to the positive orthant with respect to its distance to $\mathcal{X}$.

**Corollary 3.6** Assume the framework defined in the statement of lemma 3.4. Then, there exist two constants $\tau > 0$ and $\gamma > 0$ such that for all $x \in \mathbb{R}^n$,

$$ \mathrm{dist}(x, \mathbb{R}_+^n) \le \tau \mathrm{dist}(x, \mathcal{X}) \quad \implies \quad \mathrm{dist}(x, \mathcal{A}) \ge \gamma \mathrm{dist}(x, \mathcal{X}). $$

**PROOF.** Let $\gamma$ be the constant given by lemma 3.4 and set

$$ \tau := \frac{\gamma}{4(1 + \gamma)}. $$

Let $x$ be such that $\mathrm{dist}(x, \mathbb{R}_+^n) \le \tau \mathrm{dist}(x, \mathcal{X})$. To simplify the notation, let us define

$$ x^0 := P_{\mathcal{X}}(x), \quad x^1 := P_{\mathcal{A}}(x), \quad \text{and} \quad \bar{x} := P_{\mathbb{R}_+^n}(x), $$

Using several times the triangle inequality, lemma 3.4, the non-expansiveness of the projectors $P_{\mathcal{A}}$ and $P_{\mathcal{X}}$, and the definition of $\tau$, one can write

$$ \begin{aligned} \|x - x^1\| & \ge \|\bar{x} - P_{\mathcal{A}}(\bar{x})\| - \|P_{\mathcal{A}}(\bar{x}) - P_{\mathcal{A}}(x)\| - \|x - \bar{x}\| \\ & \ge \gamma \|\bar{x} - P_{\mathcal{X}}(\bar{x})\| - 2\|x - \bar{x}\| \\ & \ge \gamma \|\bar{x} - x^0\| - (2 + \gamma)\|x - \bar{x}\| \\ & \ge \gamma \|x - x^0\| - 2(1 + \gamma)\|x - \bar{x}\| \\ & \ge \gamma \|x - x^0\| - 2\tau(1 + \gamma)\|x - x^0\| \\ & = \frac{\gamma}{2} \|x - x^0\|. \end{aligned} $$

This is the expected inequality. $\square$

The following lemma will be also useful. If $I \subset \{1, \dots, n\}$, we note $I^c$ the complementary set of $I$ in $\{1, \dots, n\}$.
**Lemma 3.7** Let A be an $m \times n$ matrix, $b \in \mathbb{R}^m$, $I \subset \{1, \dots, n\}$, and $\varphi : \mathbb{R}^n \to \mathbb{R}$ be a convex differentiable function. If $\bar{x}$ is a solution to the problem

$$ \min \{\varphi(x) : Ax = b, x_I \ge 0, x_{Ic} = 0\}, $$

then there is a subset of indices $J \subset \{1, \dots, n\}$, containing $I$, such that $\bar{x}$ is also a solution to the problem

$$ \min \{\varphi(x) : Ax = b, x_J \ge 0, x_{Jc} \le 0\}. $$

**PROOF.** The constraints of the first problem are affine, hence qualified. Therefore, there exist vectors $y \in \mathbb{R}^m$ and $s \in \mathbb{R}^n$ such that

$$ \nabla \varphi(\bar{x}) + A^\top y + s = 0, \quad \bar{x}_I \ge 0, \quad s_I \le 0, \quad s_I^\top \bar{x}_I = 0, \quad \bar{x}_{Ic} = 0. $$

Define

$$ J := I \cup \{i \in I^c : s_i \le 0\}. $$

Then

$$ \nabla \varphi(\bar{x}) + A^\top y + s = 0, \quad \bar{x}_J \ge 0, \quad s_J \le 0, \quad s_J^\top \bar{x}_J = 0, \\ \bar{x}_{Jc} \le 0, \quad s_{Jc} \ge 0, \quad s_{Jc}^\top \bar{x}_{Jc} = 0. $$

By convexity, these conditions suffice to show that $\bar{x}$ is also a solution to the second problem. $\square$

# 4 Global linear convergence of the algorithm

The global linear convergence of the AL algorithm will be shown in section 4.2 to be a consequence of the radial Lipschitz continuity of the multifunction $\partial \delta^{-1}$, the reciprocal of the subgradient of the dual function (this argument is taken from [28]). The latter property is the subject of section 4.1.

## 4.1 Radial Lipschitz continuity of the subgradient reciprocal

Two normed spaces E and F being given, a multifunction $T : E \to F$ is said to be radially Lipschitz continuous at $x_0 \in E$ with module $L \ge 0$ if for all $x \in E$ and all $y \in T(x)$, there holds $\text{dist}(y, T(x_0)) \le L\|x - x_0\|$ ("dist" denotes here the distance associated with the norm of F). Consider the multifunction

$$ \partial \delta^{-1} : g \in \mathbb{R}^m \mapsto \{\lambda \in \mathbb{R}^m : g \in \partial \delta(\lambda)\} \subset \mathbb{R}^m, $$

where $\delta$ is the dual function defined in (3.1). Clearly $\partial \delta^{-1}(0) = S_D$, the set of dual solutions. Then $\partial \delta^{-1}$ is radially Lipschitz continuous at 0 with module $L \ge 0$, if

$$ \forall \lambda \in \mathbb{R}^m, \forall g \in \partial \delta(\lambda) : \text{dist}(\lambda, S_D) \le L \|g\|. \quad (4.1) $$
Such a property is sometimes called a *global error bound* for the dual solution set $S_D$ in terms of the dual function subgradient (see the review paper by Pang [21] and the contribution of Izmailov and Solodov [13]). In this section, we show that this property holds in a weaker form: $\lambda$ has to stay at a bounded distance from $S_D$ (the Lipschitz constant $L$ depends on this distance). Nevertheless, this property still has a global nature, since $\lambda$ is not required to be close to $S_D$ and $g$ is not required to be close to 0.

To show that this property is natural, consider first a quadratic problem with only equality constraints:

$$ \left\{ \begin{array}{l} \inf_x \frac{1}{2} x^\top Q x + q^\top x \\ \qquad C x = b. \end{array} \right. \qquad (4.2) $$

It is assumed that this problem is convex ($Q \succeq 0$) and has a solution. It is therefore feasible: $b \in R(C)$. Since the constraint is qualified, there exist optimal multipliers, which implies that the affine subspace $\Lambda$ defined in (2.1) is nonempty.

Using the pseudo-inverse $Q^\dagger$ of $Q$, the symmetric matrix $H \succeq 0$, and the vector $v$ defined in (3.5), the dual function $\delta$ associated with problem (4.2) can be written

$$ \delta(\lambda) = \begin{cases} \frac{1}{2}\lambda^\top H \lambda - (v+b)^\top \lambda + \frac{1}{2}q^\top Q^\dagger q & \text{for } \lambda \in \Lambda \\ +\infty & \text{otherwise.} \end{cases} \qquad (4.3) $$

A computation like in the proof of lemma 3.1 shows that

$$ \partial \delta(\lambda) = H\lambda - v - b + C(N(Q)), \quad \text{for } \lambda \in \Lambda. $$

Since $S_D$ is defined as the set of minimizers of $\delta$, one finds

$$ S_D = \{\bar{\lambda} \in \Lambda : H\bar{\lambda} \in v + b + C(N(Q))\}. $$

**Proposition 4.1** *Consider problem (4.2) with $Q \succeq 0$ and suppose that it has a solution. Then property (4.1) is satisfied by the dual function (4.3), with the Euclidean norm and a constant $L$ equal to the inverse of the smallest nonzero eigenvalue of $H$ ($L = 0$ if $H = 0$).*

**PROOF.** Since problem (4.2) has a solution and its constraint is affine, $S_D$ is nonempty (it is identical to the set of optimal multipliers). To prove (4.1), we only have to consider the dual variables $\lambda \in \Lambda$, since otherwise $\partial \delta(\lambda)$ is empty. Note also that we only have to consider the case when $H \neq 0$ since otherwise $S_D = \Lambda$ and (4.1) is trivially satisfied with $L = 0$.

Let $\lambda \in \Lambda$, $g \in \partial \delta(\lambda)$, and $\bar{\lambda}$ be the projection of $\lambda$ onto $S_D$. We have for some $u$ and $\bar{u} \in N(Q)$

$$ g = H\lambda - (v+b) + C\bar{u} \quad \text{and} \quad 0 = H\bar{\lambda} - (v+b) + C\bar{u}. $$

Subtracting these two identities and taking the scalar product with $(\lambda - \bar{\lambda})$ yield

$$ g^\top(\lambda - \bar{\lambda}) = (\lambda - \bar{\lambda})^\top H(\lambda - \bar{\lambda}) + (u - \bar{u})^\top C^\top(\lambda - \bar{\lambda}). $$
The last term vanishes, since $C^\top(\lambda - \bar{\lambda}) \in R(Q)$ and $Q(u - \bar{u}) = 0$. Now, since $\bar{\lambda} + N(H) \subset S_D$, $\lambda - \bar{\lambda} \in N(H)^\perp$. Therefore

$$g^\top(\lambda - \bar{\lambda}) \geq \left( \inf_{\substack{w \in N(H)^\perp \\ \|w\|=1}} w^\top H w \right) \|\lambda - \bar{\lambda}\|^2 = \frac{1}{L} \|\lambda - \bar{\lambda}\|^2.$$

Now (4.1) follows by using the Cauchy-Schwarz inequality on the left hand side. □

When $C$ is surjective and $Q \succ 0$, extending this result to the dual function associated with the strictly convex quadratic problem (1.2) is an easy exercise (then $H$ is positive definite and $L$ is the inverse of the smallest eigenvalue of $H$), but it presents little interest in practice. On the other hand, when $C$ is not surjective, property (4.1) cannot hold without being lightly weakened, as shown by the following example.

**Example 4.2** Consider indeed the special case with a single inactive constraint ($m = 1$, $C = 0$, and $l < 0 < u$) and a zero optimal value ($q = 0$). Then $\delta$ is the function

$$\delta(\lambda) = \begin{cases} -u\lambda & \text{if } \lambda \le 0 \\ -l\lambda & \text{if } \lambda > 0. \end{cases}$$

Clearly, (4.1) can hold only if $\lambda$ is not too far from the dual solution set: $|\lambda| \le L \min(-l, u)$. □

The analysis of the inequality constrained QP is more difficult since it has to cover simultaneously two extreme cases: the quadratic dual function of the equality constrained QP and the sharp dual function of the previous example.

In the case of an inequality constrained QP, the Lipschitz constant $L$ will be shown to depend on the gap $\Delta$ between $\bar{y}_J$ (for some $\bar{y} \in S_P^y$) and the bounds $l_J$ and $u_J$ (see (3.7) for the definition of the index set $J$). More precisely, $\Delta$ is defined by

$$\Delta := \sup_{\bar{y} \in S_P^y} \min \left( \min_{i \in I_l \cup J} (u_i - \bar{y}_i), \min_{i \in J \cup I_u} (\bar{y}_i - l_i) \right). \quad (4.4)$$

If $J = \emptyset$, either $I_l$ or $I_u \neq \emptyset$, and the fact that $l < u$ implies that $\Delta > 0$. If $J \neq \emptyset$, the convexity of $S_P$ implies that there is a $\bar{y} \in S_P^y$ such that $l_J < \bar{y}_J < u_J$, in which case also $\Delta > 0$. The dependence of $L$ on $\Delta$ is clearly visible in example 4.2: for $\lambda$ at a unit distance from the solution, we must have $L \ge 1/\min(-l, u) = 1/\Delta$. This lower bound on $L$ goes to infinity when $l$ or $u$ tends to zero, and it goes to zero when $l \to -\infty$ and $u \to +\infty$.

**Proposition 4.3** For any bounded set $B \subset \mathbb{R}^m$, there exists a constant $L$, such that

$$\forall \lambda \in S_D + B, \forall g \in \partial\delta(\lambda) : \mathrm{dist}(\lambda, S_D) \le L\|g\|. \quad (4.5)$$
PROOF. First stage: definition of L. We assume that $S_D \neq \emptyset$, since otherwise there is nothing to prove. Let $\mathcal{B}$ be a bounded set in $\mathbb{R}^m$, i.e., $\mathcal{B} \subset \beta\mathcal{B}$ for some $\beta > 0$. To make the proof rigorous, we now define $L > 0$, even though the motivation for its definition will not look quite clear at this point.

Let $\mathcal{K}$ be the collection of index sets $K \subset \{1, \dots, m\}$ such that $I_l \subset K$ and $I_u \subset K^c := \{1, \dots, m\} \setminus K$ (the index sets $I_l$ and $I_u$ are defined in (3.7)). With any index set $K \subset \{1, \dots, m\}$, we associate the orthant

$$ \mathcal{O}_K := \{\lambda \in \mathbb{R}^m : \lambda_K \ge 0, \lambda_{K^c} \le 0\}. $$

Define $\mathcal{O}$ and $\mathcal{A}$ by (3.8). For any index set $K \in \mathcal{K}$, $\mathcal{O}_K \cap \mathcal{A}$ is nonempty (since it contains $S_D = \mathcal{O} \cap \mathcal{A}$, see lemma 3.2). Therefore, with an index set $K \in \mathcal{K}$, corollary 3.6 associates two constants $\tau_K > 0$ and $\gamma_K > 0$ such that for any $\lambda \in \mathbb{R}^m$:

$$ \mathrm{dist}(\lambda, \mathcal{O}_K) \le \tau_K \mathrm{dist}(\lambda, \mathcal{O}_K \cap \mathcal{A}) \implies \mathrm{dist}(\lambda, \mathcal{A}) \ge \gamma_K \mathrm{dist}(\lambda, \mathcal{O}_K \cap \mathcal{A}). $$

Since $\mathcal{K}$ is finite, the constants

$$ \tau := \min_{K \in \mathcal{K}} \tau_K \quad \text{and} \quad \gamma := \min_{K \in \mathcal{K}} \gamma_K $$

are positive. Therefore, we have found two constants $\tau > 0$ and $\gamma > 0$ such that, for any $K \in \mathcal{K}$, there holds

$$ \begin{aligned} \lambda \in \mathcal{O}_K^\tau &:= \{\lambda' \in \mathbb{R}^m : \mathrm{dist}(\lambda', \mathcal{O}_K) \le \tau \mathrm{dist}(\lambda', \mathcal{O}_K \cap \mathcal{A})\} \\ &\implies \mathrm{dist}(\lambda, \mathcal{A}) \ge \gamma \mathrm{dist}(\lambda, \mathcal{O}_K \cap \mathcal{A}). \end{aligned} \quad (4.6) $$

We also introduce

$$ \sigma := \inf_{\substack{\mu \in \partial B \cap R(C) \\ C^\top \mu \in R(Q^\dagger)}} \mu^\top H \mu, \quad (4.7) $$

which is therefore $+\infty$ if $\{\mu \in R(C) : C^\top \mu \in R(Q^\dagger)\} = \{0\}$. Otherwise, it is a positive number (indeed, then $Q^\dagger \neq 0$ and $C \neq 0$; on the other hand, since $C^\top \mu \in R(Q^\dagger) = N(Q^\dagger)^\perp$, $\mu^\top H \mu = \mu^\top CQ^\dagger C^\top \mu \ge \zeta_{\min}(Q^\dagger)\|C^\top \mu\|^2$, where $\zeta_{\min}(Q^\dagger)$ is the smallest nonzero eigenvalue of $Q^\dagger$; finally, since $\mu \in R(C)$, $\|C^\top \mu\| \ge \sigma_{\min}(C)\|\mu\|$, where $\sigma_{\min}(C)$ is the smallest nonzero singular value of $C$). Note that $\sigma$ is the smallest nonzero eigenvalue of $H$ when $Q > 0$.

Recall the definition (4.4) of $\Delta > 0$. Then, the constant $L > 0$ is defined by (when $\sigma = +\infty$, $L$ is set to the second argument of the max below)

$$ L := \max \left( \frac{1}{\sigma\gamma^2}, \frac{\beta}{\tau\Delta} \right). \quad (4.8) $$

Second stage: proof of (4.5). Fix $\lambda \in S_D + \mathcal{B}$ and $g \in \partial\delta(\lambda)$ (necessarily, $\lambda \in \mathrm{dom}\,\delta$). Denote the projection of $\lambda$ onto $S_D$ by $\bar{\lambda} := P_{S_D}(\lambda)$. Observe that,

$$ \|\lambda - \bar{\lambda}\| \le \beta. \quad (4.9) $$
Let $\varepsilon \in ]0, \Delta[$ and define $L_\varepsilon$ by formula (4.8), but with $\Delta - \varepsilon$ in place of $\Delta$. Observe now that showing

$$g^\top(\lambda - \bar{\lambda}) \geq \frac{1}{L_\varepsilon} \|\lambda - \bar{\lambda}\|^2 \quad (4.10)$$

suffices to conclude the proof since then the inequality in (4.5) follows from the Cauchy-Schwarz inequality applied to the left hand side of (4.10) and the fact that $\varepsilon$ can be chosen arbitrarily close to zero.

From the form of the subdifferential $\partial\delta(\lambda)$ given by lemma 3.1, we have for some $y_\lambda \in Y_\lambda$, some $y_{\bar{\lambda}} \in Y_{\bar{\lambda}}$, and some $u, \bar{u} \in N(Q)$:

$$g = H\lambda - v - y_{\lambda} + Cu \quad \text{and} \quad 0 = H\bar{\lambda} - v - y_{\bar{\lambda}} + C\bar{u}.$$

According to lemma 3.2, $y_{\bar{\lambda}}$ can be chosen arbitrarily in $S_P^y$ and we take it such that

$$\min \left( \min_{i \in I_l \cup J} (u_i - \bar{y}_i), \min_{i \in J \cup I_u} (\bar{y}_i - l_i) \right) \geq \Delta - \varepsilon. \quad (4.11)$$

As in the proof of proposition 4.1, $(u - \bar{u})^\top C^\top(\lambda - \bar{\lambda}) = 0$, because $C^\top(\lambda - \bar{\lambda}) \in R(Q)$ (both $\lambda$ and $\bar{\lambda} \in \text{dom } \delta \subset \Lambda$) and $Q(u - \bar{u}) = 0$. Therefore, subtracting the identities above and taking the scalar product with $(\lambda - \bar{\lambda})$ yield

$$g^\top(\lambda - \bar{\lambda}) = (\lambda - \bar{\lambda})^\top H(\lambda - \bar{\lambda}) - (y_\lambda - y_{\bar{\lambda}})^\top(\lambda - \bar{\lambda}). \quad (4.12)$$

We will get (4.10) by finding a lower bound of the right hand side of (4.12). Note that the two terms are nonnegative (this is clear for the first one, since $H$ is positive semi-definite; for the second one, use the monotonicity property (3.4)).

Since $\bar{\lambda} = P_{A \cap O}(\lambda)$, by lemma 3.7, one can find an index set $K \subset \mathcal{K}$ such that $\bar{\lambda} = P_{A \cap O_K}(\lambda)$. We analyze successively two complementary cases, using the set $O_K^\tau$ defined in (4.6) and $\lambda^t := (1-t)\bar{\lambda} + t\lambda$ for $t \in \mathbb{R}$.

Case A: there exists a $t \in ]0, 1]$ such that $\lambda^t \in O_K^\tau$. In this case, we work on the first term in the right hand side of (4.12), discarding the second one. Because $P_{A \cap O_K}(\lambda^t) = \bar{\lambda}$, (4.6) gives

$$\gamma \|\lambda^t - \bar{\lambda}\| \leq \|\lambda^t - P_A(\lambda^t)\|. $$

Decompose $\lambda^t - \bar{\lambda} = \mu_0 + \mu_1$, where $\mu_0 \in N(C^\top)$ and $\mu_1 \in R(C)$, and observe that $C^\top\mu_1 = C^\top(\lambda^t - \bar{\lambda}) \in R(Q) = R(Q^\dagger)$ (since both $\lambda^t$ and $\bar{\lambda} \in \text{dom } \delta \subset \Lambda$). Then, using the definition (4.7) of $\sigma$, one finds

$$(\lambda^t - \bar{\lambda})^\top H(\lambda^t - \bar{\lambda}) = \mu_1^\top H \mu_1 \geq \sigma \|\mu_1\|^2.$$

Clearly $\mu_1 = \lambda^t - P_A(\lambda^t)$, so that

$$(\lambda^t - \bar{\lambda})^\top H(\lambda^t - \bar{\lambda}) \geq \sigma \gamma^2 \|\lambda^t - \bar{\lambda}\|^2.$$
Since $\lambda - \bar{\lambda} = (\lambda^t - \bar{\lambda})/t$, we also have

$$ (\lambda - \bar{\lambda})^{\top} H (\lambda - \bar{\lambda}) \geq \sigma \gamma^2 \|\lambda - \bar{\lambda}\|^2. $$

Discarding the second term in the right hand side of (4.12) (it is nonnegative) and using the definition of $L$ in (4.8), it follows that

$$ g^{\top}(\lambda - \bar{\lambda}) \geq (\lambda - \bar{\lambda})^{\top} H(\lambda - \bar{\lambda}) \geq \sigma \gamma^2 \|\lambda - \bar{\lambda}\|^2 \geq \frac{1}{L} \|\lambda - \bar{\lambda}\|^2, $$

which is the expected inequality (4.10) (of course $L_{\varepsilon} \geq L$).

Case B: for any $t \in ]0, 1]$, $\lambda^t \notin O_K^\tau$. In this case, we work on the second term in the right hand side of (4.12), discarding the first one. Let us start by choosing $t \in ]0, 1]$ sufficiently small such that $\lambda_i^t \bar{\lambda}_i > 0$ when $\bar{\lambda}_i \neq 0$; by assumption, this $\lambda^t \notin O_K^\tau$. Denote by $\lambda_K^t := P_{O_K}(\lambda^t)$ the projection of $\lambda^t$ onto $O_K$ and decompose

$$ -(y_{\lambda^t} - y_{\bar{\lambda}})^{\top}(\lambda^t - \bar{\lambda}) = -(y_{\lambda^t} - y_{\bar{\lambda}})^{\top}(\lambda^t - \lambda_K^t) - (y_{\lambda^t} - y_{\bar{\lambda}})^{\top}(\lambda_K^t - \bar{\lambda}). $$

The last term in the right hand side is nonnegative. Indeed, by the choice of $t$, if $\bar{\lambda}_i \neq 0$, one has $\lambda_i^t \bar{\lambda}_i > 0$ and therefore $(y_{\lambda^t} - y_{\bar{\lambda}})_i = 0$ (see the definition (3.3) of $Y_\lambda$). The only nonzero terms of the last scalar product are therefore of the form $(y_{\bar{\lambda}} - y_{\lambda^t})_i(\lambda_K^t)_i$. If $(\lambda_K^t)_i > 0$, one has $\lambda_i^t > 0$ and therefore $(y_{\lambda^t})_i = l_i$, so that the term can be written $((y_{\bar{\lambda}})_i - l_i)(\lambda_K^t)_i \geq 0$ (since $l_i \leq (y_{\bar{\lambda}})_i \leq u_i$). Similarly, the term is nonnegative when $(\lambda_K^t)_i < 0$. Therefore

$$ -(y_{\lambda^t} - y_{\bar{\lambda}})^{\top}(\lambda^t - \bar{\lambda}) \geq -(y_{\lambda^t} - y_{\bar{\lambda}})^{\top}(\lambda^t - \lambda_K^t) = \sum_{i \in I_{K, \lambda^t}} (y_{\bar{\lambda}} - y_{\lambda^t})_i (\lambda^t - \lambda_K^t)_i, $$

where we have introduced the index set

$$ I_{K, \lambda^t} := \{i : \lambda_i^t \neq (\lambda_K^t)_i\}. $$

Let us show that all the terms of the sum on the indices $i \in I_{K, \lambda^t}$ above are positive. Observe first that $(\lambda_K^t)_i = 0$ (since $\lambda_i^t \neq (\lambda_K^t)_i$ and $\lambda_K^t$ is the projection of $\lambda^t$ onto the orthant $O_K$). On the other hand, if $\lambda_i^t > 0$, $(y_{\lambda^t})_i = l_i$ and the fact that $\lambda_i^t > 0$ and $(\lambda_K^t)_i = 0$ implies that $i \in K^c$, so that $l_i < (y_{\bar{\lambda}})_i \leq u_i$ (since $I_l \subset K$); therefore $(y_{\bar{\lambda}} - y_{\lambda^t})_i = (y_{\bar{\lambda}})_i - l_i \geq \Delta - \varepsilon > 0$ (see (4.11)). Similarly, if $\lambda_i^t < 0$, then $(y_{\lambda^t})_i = u_i$, $i \in K$, and $(y_{\bar{\lambda}} - y_{\lambda^t})_i = (y_{\bar{\lambda}})_i - u_i \leq -(\Delta - \varepsilon) < 0$. In particular, all the terms of the sum are positive. Therefore, if we remind us that $\lambda^t \notin O_K^\tau$, we get ($\| \cdot \|_1$ denotes the $\ell_1$-norm)

$$ -(y_{\lambda^t} - y_{\bar{\lambda}})^{\top}(\lambda^t - \bar{\lambda}) \geq (\Delta - \varepsilon) \|\lambda^t - \lambda_K^t\|_1 \geq (\Delta - \varepsilon) \|\lambda^t - \lambda_K^t\| \geq \tau (\Delta - \varepsilon) \|\lambda^t - \bar{\lambda}\|. \quad (4.13) $$
We can now conclude. For $g^t \in \partial\delta(\lambda^t)$, there holds

$$
\begin{align*}
g^{\top}(\lambda - \bar{\lambda}) &\ge (g^t)^{\top}(\lambda - \bar{\lambda}) \quad \text{[monotonicity of the subdifferential]} \\
&= \frac{1}{t}(g^t)^{\top}(\lambda^t - \bar{\lambda}) \quad \text{[definition of } \lambda^t \text{]} \\
&\ge -\frac{1}{t}(y_{\lambda^t} - y_{\bar{\lambda}})^{\top}(\lambda^t - \bar{\lambda}) \quad \text{[(4.12)]} \\
&\ge \frac{\tau(\Delta - \varepsilon)}{t} \|\lambda^t - \bar{\lambda}\| \quad \text{[(4.13)]} \\
&\ge \tau(\Delta - \varepsilon) \|\lambda - \bar{\lambda}\| \quad \text{[definition of } \lambda^t \text{]} \\
&\ge \frac{\tau(\Delta - \varepsilon)}{\beta} \|\lambda - \bar{\lambda}\|^2 \quad \text{[(4.9)]} \\
&\ge \frac{1}{L_{\varepsilon}} \|\lambda - \bar{\lambda}\|^2 \quad \text{[definition of } L_{\varepsilon} \text{].}
\end{align*}
$$

This is the expected inequality (4.10). □

## 4.2 Global linear convergence

We can now state the global linear convergence of the constraint norm towards zero in the AL algorithm of section 2. Note that the rate of convergence $\min(L/r_k, 1)$ depends through $L$ on the distance from the initial iterate $\lambda_0$ to the dual solution set $S_D$.

**Theorem 4.4** Consider the AL algorithm of section 2. For any $\beta > 0$, there exists an $L > 0$, such that $\text{dist}(\lambda_0, S_D) \le \beta$ implies that

$$ \|y_{k+1} - Cx_{k+1}\| \le \min\left(\frac{L}{r_k}, 1\right) \|y_k - Cx_k\|, \quad \text{for all } k \ge 1. \quad (4.14) $$

In particular, if $r_k \ge \bar{r}$ for all $k \ge 1$ and some $\bar{r} > L$, the constraint norm tends to zero globally linearly.

**PROOF.** The proof gathers known techniques (see for example [27, 28]) with the result of proposition 4.3. We give the details for completeness.

Let us note $g_{k+1} := Cx_{k+1} - y_{k+1}$. Recall from (3.10) that $g_{k+1} \in \partial\delta(\lambda_{k+1})$. Subtracting two consecutive iteration identities (2.4) provides

$$ \frac{1}{r_{k+1}}(\lambda_{k+2} - \lambda_{k+1}) + (g_{k+2} - g_{k+1}) = \frac{1}{r_k}(\lambda_{k+1} - \lambda_k). $$

Taking norms, using the monotonicity of the subdifferential, and discarding $\|g_{k+2} - g_{k+1}\|^2 \ge 0$, we get $\|\lambda_{k+2} - \lambda_{k+1}\|^2/r_{k+1}^2 \le \|\lambda_{k+1} - \lambda_k\|^2/r_k^2$ or $\|g_{k+2}\| \le \|g_{k+1}\|$, which yields the second part of the min in (4.14).

Subtracting an arbitrary dual solution $\bar{\lambda} \in S_D$ from the iteration identity (2.4), taking norms, and using the monotonicity of the subdifferential lead to

$$ \|\lambda_{k+1} - \bar{\lambda}\|^2 + r_k^2 \|g_{k+1}\|^2 \le \|\lambda_k - \bar{\lambda}\|^2, \quad \text{for } k \ge 0. \quad (4.15) $$
This shows in particular that the sequence {$\lambda_k - \bar{\lambda}$}$_{k \ge 0}$ is nonincreasing. Since $\bar{\lambda}$ is arbitrary in $S_D$, {$\text{dist}(\lambda_k, S_D)$}$_{k \ge 0}$ is also nonincreasing, so that $\lambda_k \in S_D + \beta B$ for all $k \ge 0$. Now, let $L > 0$ be the constant that proposition 4.3 associates with $B := \beta B$. By this proposition, $\|\lambda_k - P_{S_D}(\lambda_k)\| \le L\|g_k\|$. Now discarding the first term in the left hand side of (4.15) and using $P_{S_D}(\lambda_k)$ for $\bar{\lambda}$, we get $\|g_{k+1}\| \le (L/r_k)\|g_k\|$. This yields the first part of the min in (4.14). □

# 5 Numerical experiments and discussion

The aim of this section is to illustrate by numerical experiments the global linear convergence property of the AL algorithm studied in this paper and to assess the quality of the bound given by theorem 4.4. The numerical experiments are taken from seismic reflection tomography applications. We conclude with a discussion on algorithmic implications.

## 5.1 A seismic reflection tomography problem

Seismic reflection tomography is a technique used to recover the geological structure of the subsoil from the measurements of the travel-times of seismic waves (see [10] for a description of the approach). From an optimization viewpoint, the problem consists in minimizing a nonlinear least-squares function subject to nonlinear constraints. In [5], a Gauss-Newton SQP method globalized by line-search is proposed and analyzed. At each iteration, a solution to a strictly convex quadratic model of the objective function subject to linear constraints is computed using our code QPAL [4, 6].

We have chosen here to present the results obtained with the problem KARINE, which is representative of those observed with our collection of 2D and 3D seismic reflection problems. These have a number of variables up to $15 \cdot 10^3$ and a number of constraints up to $10^4$. The features of the selected problem are summarized in table 1. It is a 2D model

<table><caption>Table 1: Description of the tomography problem KARINE</caption><thead><tr><th>n</th><th>m</th><th>m<sup>*</sup><sub>act</sub></th><th>&kappa;<sub>2</sub></th></tr></thead><tbody><tr><td>442</td><td>320</td><td>108</td><td>8.4 10<sup>5</sup></td></tr></tbody></table>

depending on $n = 442$ parameters and having $m = 320$ linear inequality constraints. The matrix $Q$ of the selected quadratic subproblem (1.1) is positive definite and has its $l_2$ condition number equal to $\kappa_2 = 8.4 \cdot 10^5$. Its constraint matrix $C$ has been balanced (the Euclidean norm of each of its rows is equal to 1). The number of active constraints at the solution is $m^*_{act} = 108$, which represents 33% of the number of constraints.

The results presented in section 5.2 have been obtained using the AL algorithm described in section 2, with a fixed augmentation parameter $r$. In order to study the dependence of the results on $r$, we have run the QP solver for 21 different values of $r$, ranging from $1$ to $10^5$. In each case, the AL algorithm is initialized with a null
Lagrangian multiplier ($\lambda_0 = 0$) and is stopped when the constraint norm is sufficiently small ($\|y_k - Cx_k\| \le 10^{-10}$).

## 5.2 Assessing the global linear convergence result

In this section, we illustrate the global linear convergence property of the AL algorithm established in theorem 4.4. The actual global rate of linear convergence is given by $\rho := \sup\{\|y_{k+1} - Cx_{k+1}\|/\|y_k - Cx_k\| : k \ge 1\}$ and can be estimated during a particular run by

$$ \rho_{\text{est}} := \max_{1 \le k \le n_{AL}} \frac{\|y_{k+1} - Cx_{k+1}\|}{\|y_k - Cx_k\|} \le \rho, \quad (5.1) $$

where $n_{AL}$ is the number of AL iterations actually performed to reach the required accuracy of $10^{-10}$ on the constraint norm.

Theorem 4.4 has shown that $\rho$ is bounded above by a function of $r$:

$$ \rho \le \min\left(\frac{L}{r}, 1\right) \quad \text{or} \quad \log \rho \le \min(\log L - \log r, 0). \quad (5.2) $$

A natural question is to know whether this bound is tight in practice. This is difficult to say, since the value of $L$ is generally unknown, but the appearance of $\rho_{\text{est}}$ as a function of $r$ in the considered problem may give a clue on this question.

<figure><img src="image_5.png" /><figcaption>Figure 1: Global linear convergence rate of the constraint norm as a function of $r$</figcaption></figure>

The plain line in figure 1 gives $\log \rho_{\text{est}}$ as a function of $\log r$ (double logarithmic scale). As predicted by the theory, we see that $\rho_{\text{est}} \le 1$ for all positive $r$. Furthermore, the larger is the augmentation parameter $r$, the faster is the convergence: $\rho \approx 1$ for $r \le 10$ (convergence is hardly detectable) and $\rho \approx 3.10^{-3}$ for $r = 10^5$ (convergence is obtained in very few AL iterations). We have represented by a dotted line the tangent
to the $\rho_{est}$ curve with a slope $-1$. This line crosses the top horizontal line of the graph at the horizontal coordinate $L_{inf} \approx 304$. According to (5.2), $L_{inf}$ provides a lower estimate of the value of $L$. As both curves (the plain and dotted ones) are quite close, it is likely that the dotted curve is close to the upper bound on $\rho$ given by (5.2). As a result, it is likely that the upper bound given by theorem 4.4 is tight. Note that the small discrepancy between both curves for large values of $r$ ($r \ge 7.5 10^4$) is due here to the fact that the AL algorithm reaches the required constraint norm accuracy in very few AL iterations ($n_{AL} \le 3$), so that the maximum in (5.1) is taken on that few number. In other cases, such a discrepancy can come from an inexact solve of the bound constraint problem (2.3), due to a large value of $r$.

## 5.3 Discussion

As shown in this paper, the global linear rate of convergence of the AL algorithm depends on the Lipschitz constant $L$ given by (4.8) and on the value of $\bar{r} := \inf r_k$, where $r_k$ is the current value of the augmentation parameter. More precisely, the decrease of the constraint norm at iteration $k$ is bounded above by $L/r_k$. It is usually impossible to compute $L$ in practice, since it depends on the constants $\gamma$, $\Delta$, and $\sigma$ (see lemma 3.4, (4.4), (4.7), and finally (4.8)), which are either unknown or too expensive to compute. As a Lipschitz constant, however, $L$ has easily computable lower estimates.

The estimate $L_{inf}$ of $L$ given in section 5.2 is not available at run time, since it requires to run the AL algorithm on a particular problem for various values of $r$. Nevertheless, the quantities

$$L_{inf,k} := \max_{1 \le i \le k} \left( r_i \frac{\|y_{i+1} - Cx_{i+1}\|}{\|y_i - Cx_i\|} \right)$$

satisfy $L_{inf,k} \le L$ and can therefore be used as a lower estimate of $L$, after iteration $k$ is completed. A given desired rate of convergence $\rho_{des} \in ]0, 1[$ is then likely to be obtained at iteration $k+1$ by taking

$$r_{k+1} \ge \frac{L_{inf,k}}{\rho_{des}}. \qquad (5.3)$$

It is the fact that the estimate (4.14) has a global validity that gives sense to an update of the value of $r_k$ in this way at each iteration. It should be clear at this point that the AL algorithm gains in efficiency by taking $r_k$ as large as possible, the only limitation being that problem (2.3) needs to be numerically solvable. Since it is sometimes difficult to tell what is a large value for a particular problem, the lower bound on $r_{k+1}$ in (5.3) may also be useful as a reference.

We conclude with a result providing an estimate of the number of iterations needed to reach a given tolerance on the constraint norm. Assume that a number $\rho_{des} \in ]0, 1[$ is given as a desired rate of convergence. Of course, since the Lipschitz constant $L$ is unknown, this rate of convergence cannot be ensured, but the algorithm can try to approach it by updating $r_k$ when it feels it is necessary. The next result gives an
estimate of the iterative complexity of the AL algorithm with an update rule based on (5.3). More precisely, defining

$$ \rho_k := \frac{\|y_{k+1} - Cx_{k+1}\|}{\|y_k - Cx_k\|}, $$

the AL algorithm is supposed to update the value of $r_k$, for $k \ge 1$, according to:

$$ \text{if } \rho_k \le \rho_{\text{des}}, \text{ then } r_{k+1} = r_k, \text{ else } r_{k+1} = \frac{\rho_k}{\rho_{\text{des}}} r_k. \quad (5.4) $$

There is nothing magic in the update rule of $r_k$ above. It could equally use $r_{k+1} = 10 \rho_k r_k / \rho_{\text{des}}$ or simply $r_{k+1} = 10 r_k$ when $r_k$ needs to be increased.

**Proposition 5.1** Suppose that the AL algorithm of section 2 uses the rule (5.4) to update the augmentation parameter $r_k$, for $k \ge 1$. Let $\varepsilon \in ]0, 1]$ and let $L$ be the positive constant given by theorem 4.4. Fix any $t \in ]\rho_{\text{des}}, 1[$. Then

$$ \|y_{k+1} - Cx_{k+1}\| \le \varepsilon \|y_1 - Cx_1\|, \quad (5.5) $$

as soon as

$$ k \ge \frac{\log \varepsilon}{\log t} + \max \left( 1 + \frac{\log(L/(tr_1))}{\log(t/\rho_{\text{des}})}, 0 \right). $$

**PROOF.** Let $t \in ]\rho_{\text{des}}, 1[$. Clearly, since $\rho_i \le 1$,

$$ \frac{\|y_{k+1} - Cx_{k+1}\|}{\|y_1 - Cx_1\|} = \prod_{1 \le i \le k} \rho_i \le \prod_{\substack{1 \le i \le k \\ \rho_i \le t}} \rho_i \le t^{k_t}, $$

where $k_t := |K_t|$ is the number of elements in $K_t := \{i \in \mathbb{N} : 1 \le i \le k, \rho_i \le t\}$. Taking logarithms, we see that (5.5) holds as soon as $k_t \ge (\log \varepsilon)/(\log t)$.

If $K_t^c := \{1, \dots, k\} \setminus K_t$ is empty, then $k = k_t$ and the result is proven.

Suppose now that $K_t^c \ne \emptyset$. Since $\rho_i \le L/r_i$ (by theorem 4.4), $i \in K_t$ as soon as $r_i \ge L/t$. Then the last index in $K_t^c$, namely the $(k-k_t)$th one, has a value of $r_i$ updated $k-k_t-1$ times from $r_1$ using a factor $\rho_i/\rho_{\text{des}} \ge t/\rho_{\text{des}}$ (see the update rule (5.4)). Hence we must have $(t/\rho_{\text{des}})^{k-k_t-1}r_1 \le L/t$. This gives an upper bound on the number of elements of $K_t^c$, namely

$$ k - k_t \le 1 + \frac{\log(L/(tr_1))}{\log(t/\rho_{\text{des}})}. $$

The total number of iterations to satisfy (5.5) is therefore at most this upper bound plus the lower bound on $k_t$ obtained above. □

Roughly expressed, the number of iterations needed to reach precision $\varepsilon > 0$ on the relative constraint norm is of order $O(\log \varepsilon) + O(\log L)$. As shown in the proof of proposition 5.1, the first term of order $O(\log \varepsilon)$ is due to the linear convergence of the constraint norm towards zero, which is triggered when the augmentation parameter is
large enough (a consequence of theorem 4.4). The second term of order $O(\log L)$ is
due to a possible too small value of $r_1$ and to the number of iterations that the rule
(5.4) needs to make $r_k$ large enough. This term can be made as small as desired by
choosing a large value for $r_1$ or by adopting an update rule of $r_k$ that increases these
values more rapidly than in (5.4). As a result, the computational complexity of the AL
algorithm of section 2 essentially rests on the one of the AL subproblems (2.3). When
strict complementarity holds, the finite identification of the active constraints in (2.3)
occurs and the computational complexity is then basically induced by the very first AL
subproblems. Our experience with the AL algorithm, limited to the seismic reflection
tomography problems described in section 5.1, supports that conclusion.

References

[1] K.J. Arrow, R.M. Solow (1958). Gradient methods for constrained maxima with weak-ened assumptions. In K.J. Arrow, L. Hurwicz, H. Uzawa (editors), *Studies in Linear and Nonlinear Programming*, pages 166–176. Stanford University Press, Stanford, Calif.

[2] J.F. Bonnans, J.Ch. Gilbert, C. Lemaréchal, C. Sagastizábal (2003). *Numerical Optimization – Theoretical and Practical Aspects*. Springer Verlag, Berlin

[3] J.V. Burke, M.C. Ferris (1991). Characterization of solution sets of convex programs. *Operations Research Letters*, 10, 57-60

[4] F. Delbos (2004). *Problèmes d'Optimisation de Grande Taille avec Contraintes en Tomographie de Réflexion 3D*. PhD Thesis, University Pierre et Marie Curie (Paris VI), Paris. (to appear)

[5] F. Delbos, J.Ch. Gilbert, R. Glowinski, D. Sinoquet (2003). Constrained optimization in reflection tomography: an augmented Lagrangian SQP approach. Technical report, Institut Français du Pétrole. (to appear)

[6] F. Delbos, J.Ch. Gilbert, D. Sinoquet (2003). QPAL: a solver of large-scale convex quadratic optimization problems using an augmented Lagrangian approach. Technical report, INRIA, BP 105, 78153 Le Chesnay, France. (to appear)

[7] D. den Hertog (1992). *Interior Point Approach to Linear, Quadratic and Convex Programming*. Mathematics and its Applications 277. Kluwer Academic Publishers, Dordrecht

[8] M. Fortin, R. Glowinski (1982). *Méthodes de Lagrangien Augmenté – Applications à la Résolution Numérique de Problèmes aux Limites*. Méthodes Mathématiques de l'Informatique 9. Dunod, Paris

[9] A. Friedlander, J.M. Martínez (1994). On the maximization of a concave quadratic function with box constraints. *SIAM Journal on Optimization*, 4, 177–192

[10] R. Glowinski, Q.-H. Tran (1993). Constrained optimization in reflexion tomography: the augmented Lagrangian method. *East-West J. Numer. Math.*, 1(3), 213–234

[11] M.R. Hestenes (1969). Multiplier and gradient methods. *Journal of Optimization Theory and Applications*, 4, 303–320

[12] J.-B. Hiriart-Urruty, C. Lemaréchal (1993). *Convex Analysis and Minimization Algorithms*. Grundlehren der mathematischen Wissenschaften 305-306. Springer-Verlag
[13] A.F. Izmailov, M.V. Solodov (2001). Error bounds for 2-regular mappings with Lipschitzian derivatives and their applications. *Mathematical Programming*, 89, 413–435

[14] B. Jansen (1997). *Interior Point Techniques in Optimization – Complementarity, Sensitivity and Algorithms*. Applied Optimization 6. Kluwer Academic Publishers, Dordrecht

[15] M. Kojima, N. Megiddo, T. Noma, A. Yoshise (1991). *A Unified Approach to Interior Point Algorithms for Linear Complementarity Problems*. Lecture Notes in Computer Science 538. Springer-Verlag, Berlin

[16] O.L. Mangasarian (1988). A simple characterization of solution sets of convex programs. *Operations Research Letters*, 7, 21–26

[17] J.J. Moré, G. Toraldo (1991). On the solution of large quadratic programming problems with bound constraints. *SIAM Journal on Optimization*, 1, 93–113

[18] Y.E. Nesterov, A.S. Nemirovskii (1994). *Interior-Point Polynomial Algorithms in Convex Programming*. SIAM Studies in Applied Mathematics 13. SIAM, Philadelphia

[19] J. Nocedal, S.J. Wright (1999). *Numerical Optimization*. Springer Series in Operations Research. Springer, New York

[20] J.S. Pang (1983). Methods for quadratic programming: a survey. *Computers and Chemical Engineering*, 7, 583–594

[21] J.S. Pang (1997). Error bounds in mathematical programming. *Mathematical Programming*, 79, 299–332

[22] M.J.D. Powell (1969). A method for nonlinear constraints in minimization problems. In R. Fletcher (editor), *Optimization*, pages 283–298. Academic Press, London

[23] R.T. Rockafellar (1970). *Convex Analysis*. Princeton Mathematics Ser. 28. Princeton University Press, Princeton, New Jersey

[24] R.T. Rockafellar (1971). New applications of duality in convex programming. In *Proceedings of the 4th Conference of Probability, Brasov, Romania*, pages 73–81

[25] R.T. Rockafellar (1973). A dual approach to solving nonlinear programming problems by unconstrained optimization. *Mathematical Programming*, 5, 354–373

[26] R.T. Rockafellar (1974). Augmented Lagrange multiplier functions and duality in nonconvex programming. *SIAM Journal on Control*, 12, 268–285

[27] R.T. Rockafellar (1976). Monotone operators and the proximal point algorithm. *SIAM Journal on Control and Optimization*, 14, 877–898

[28] R.T. Rockafellar (1976). Augmented Lagrangians and applications of the proximal point algorithm in convex programming. *Mathematics of Operations Research*, 1, 97–116

[29] C. Roos, T. Terlaky, J.-Ph. Vial (1997). *Theory and Algorithms for Linear Optimization – An Interior Point Approach*. John Wiley & Sons, Chichester

[30] T. Terlaky (editor) (1996). *Interior Point Methods of Mathematical Programming*. Kluwer Academic Press, Dordrecht

[31] S.J. Wright (1997). *Primal-Dual Interior-Point Methods*. SIAM Publication, Philadelphia
<figure><img src="image_6.png" /></figure>

Unité de recherche INRIA Rocquencourt
Domaine de Voluceau - Rocquencourt - BP 105 - 78153 Le Chesnay Cedex (France)

Unité de recherche INRIA Futurs : Parc Club Orsay Université - ZAC des Vignes
4, rue Jacques Monod - 91893 ORSAY Cedex (France)

Unité de recherche INRIA Lorraine : LORIA, Technopôle de Nancy-Brabois - Campus scientifique
615, rue du Jardin Botanique - BP 101 - 54602 Villers-lès-Nancy Cedex (France)

Unité de recherche INRIA Rennes : IRISA, Campus universitaire de Beaulieu - 35042 Rennes Cedex (France)

Unité de recherche INRIA Rhône-Alpes : 655, avenue de l'Europe - 38334 Montbonnot Saint-Ismier (France)

Unité de recherche INRIA Sophia Antipolis : 2004, route des Lucioles - BP 93 - 06902 Sophia Antipolis Cedex (France)

Éditeur

INRIA - Domaine de Voluceau - Rocquencourt, BP 105 - 78153 Le Chesnay Cedex (France)

http://www.inria.fr

ISSN 0249-6399