NPASA - Motivation and Global ConvergenceJ. Diffenderfer, W. W. Hager

# NPASA: An algorithm for nonlinear programming - Motivation and Global
Convergence††thanks: Submitted to the editors ADD DATE.The authors gratefully
acknowledge support by the National Science Foundation under Grant 1819002,
and by the Office of Naval Research under Grants N00014-15-1-2048 and
N00014-18-1-2100. This work was performed under the auspices of the U.S.
Department of Energy by Lawrence Livermore National Laboratory under Contract
DE-AC52-07NA27344, LLNL-JRNL-824568-DRAFT.

James Diffenderfer Lawrence Livermore National Laboratory, Livermore, CA ()
<EMAIL_ADDRESS>William W. Hager The University of Florida,
Gainesville, FL (, http://people.clas.ufl.edu/hager/<EMAIL_ADDRESS>
###### Abstract

In this paper, we present a two phase method for solving nonlinear programming
problems called Nonlinear Polyhedral Active Set Algorithm (NPASA) that has
global and local convergence guarantees under reasonable assumptions. The
first phase consists of an augmented Lagrangian method to ensure global
convergence while the second phase is designed to promote fast local
convergence by performing a balanced reduction of two error estimators for
nonlinear programs. After presenting error estimators for nonlinear programs
and our algorithm NPASA, we establish global convergence properties for NPASA.
Local quadratic convergence of NPASA is established in a companion paper [11].

###### keywords:

nonlinear programming, global convergence, local convergence

90C30, 65K05

##  1 Introduction

Some of the most successful approaches for solving nonlinear programs make use
of augmented Lagrangian, sequential quadratic programming (SQP), or interior-
point methods [28]. Early progress towards the implementation of a general
purpose algorithm for solving nonlinear programs was made for nonlinear
programs with linear constraints in the late 1970s [27] in the form of the
software MINOS. Some theoretical support for using the method of multipliers
to solve nonlinear programs was established in [30] and the first significant
general purpose algorithm for large-scale nonlinear programs with nonlinear
constraints was published in 1991 [7]. This algorithm, implemented under the
name LANCELOT [8], made use of augmented Lagrangian techniques and its success
encouraged the continued study [4, 9, 24, 25] and use [2] of augmented
Lagrangian techniques for solving nonlinear programs. Some notable
implementations of a SQP approach for solving problem nonlinear programs are
DNOPT [16], FILTERSQP [12], KNITRO [6], and SNOPT [14, 15]. As a note on the
usefulness of augmented Lagrangian methods, the SQP implementations SNOPT and
DNOPT make use of an augmented Lagrangian technique to ensure global
convergence. Originally introduced in the 1960s [37], early refinements to the
theory of SQP techniques in the 1970s [22, 29, 32] provided better performance
over augmented Lagrangian techniques by attaining superlinear convergence
rates under certain assumptions. Another approach that has been successful in
practice is interior-point methods. Originally referred to as barrier methods,
this approach was first introduced in the 1950s [13]. A recent implementation
of interior point techniques, called IPOPT [36], achieves global convergence
under few assumptions [35] and it also has a local superlinear convergence
rate under some assumptions after it was shown to not experience the Maratos
effect [34]. Benchmarking results comparing several of the listed
implementations can be found at [1].

In this paper, we present a method for solving a general nonlinearly
constrained optimization problem based on the recently developed [20]
polyhedral active set algorithm (PASA), an efficient algorithm for optimizing
a nonlinear function over a polyhedron. Let us consider the nonlinear
programming problem

(3)
$\displaystyle\begin{array}[]{cc}\displaystyle\min_{\bm{x}\in\mathbb{R}^{n}}&f(\bm{x})\\\
\text{s.t.}&\bm{h}(\bm{x})=\bm{0},\ \bm{r}(\bm{x})\leq\bm{0}\end{array}$

where $f:\mathbb{R}^{n}\to\mathbb{R}$ and
$\bm{h}:\mathbb{R}^{n}\to\mathbb{R}^{\ell}$ are nonlinear functions and
$\bm{r}:\mathbb{R}^{n}\to\mathbb{R}^{m}$ is a linear function defined by

(4) $\displaystyle\bm{r}(\bm{x}):=\bm{A}\bm{x}-\bm{b},$

where $\bm{A}\in\mathbb{R}^{m\times n}$ and $\bm{b}\in\mathbb{R}^{m}$. For
simplicity, we also denote the polyhedral constraint set for problem (3) by
$\Omega=\\{\bm{x}\in\mathbb{R}^{n}:\bm{r}(\bm{x})\leq\bm{0}\\}$. Any nonlinear
programming problem can be written in this form. It is assumed that $f$ and
$\bm{h}$ are at least continuously differentiable. The algorithm for solving
the nonlinearly constrained optimization problem (3) is based on the following
steps consisting of polyhedral constrained subproblems:

  1. 1.

An augmented Lagrangian step subject to the polyhedral constraints. This step
is used to ensure convergence to a stationary point even with a poor starting
guess.

  2. 2.

A Newton iteration applied to the nonlinear constraints. This step amounts to
a projection onto a polyhedral set.

  3. 3.

A quadratic program for estimating the multipliers.

  4. 4.

The optimization of an augmented Lagrangian over the tangent space associated
with the nonlinear constraints.

In this way, the general nonlinearly constrained optimization problem is
reduced to the solution of polyhedral constrained subproblems. Under suitable
assumptions, the solution of subproblems 2 – 4 is locally quadratically
convergent. If the starting guess is good enough, the algorithm simply loops
over the final three subproblems until the convergence tolerance is satisfied.
Hence, our algorithm for solving the nonlinearly constrained optimization
problems uses PASA to solve a series of polyhedral constrained optimization
problems and is denoted NPASA. This paper establishes a global convergence
result for NPASA, while the companion paper [11] establishes a local quadratic
convergence result.

The remainder of the paper is organized as follows. In Section 2, we discuss
the error estimators used by NPASA for solving problem (3). In Section 3, we
provide pseudocode for the NPASA algorithm together with a discussion
providing some motivation behind each step in the algorithm. Sections 4 and 5
focus on establishing convergence results for the various subproblems solved
within NPASA. We state and prove the main global convergence result for NPASA,
Theorem 6.1, in Section 6. We now specify some notation used throughout the
paper.

###  1.1 Notation

We will write $\mathbb{R}_{+}$ to denote the set $\\{x\in\mathbb{R}:x\geq
0\\}$. Scalar variables are denoted by italicized characters, such as $x$,
while vector variables are denoted by boldface characters, such as $\bm{x}$.
For a vector $\bm{x}$ we write $x_{j}$ to denote the $j$th component of
$\bm{x}$ and we write $\bm{x}^{\intercal}$ to denote the transpose of
$\bm{x}$. Outer iterations of NPASA will be enumerated by boldface characters
with subscript $k$, such as $\bm{x}_{k}$, and we will write $x_{kj}$ to denote
the $j$th component of the iterate $\bm{x}_{k}$. Inner iterations of NPASA
will be enumerated by boldface characters with subscript $i$. If $\mathcal{S}$
is a subset of indices of $\bm{x}$ then we write $\bm{x}_{\mathcal{S}}$ to
denote the subvector of $\bm{x}$ with indices $\mathcal{S}$. Additionally, for
a matrix $\bm{M}$ we write $\bm{M}_{\mathcal{S}}$ to denote the submatrix of
$\bm{M}$ with row indices $\mathcal{S}$. The ball with center $\bm{x}$ and
radius $r$ will be denoted by $\mathcal{B}_{r}(\bm{x})$. Real valued functions
are denoted by italicized characters, such as $f(\cdot)$, while vector valued
functions are denoted by boldface characters, such as $\bm{h}(\cdot)$. The
gradient of a real valued function $f(\bm{x})$ is denoted by $\nabla
f(\bm{x})$ and is a row vector. The Jacobian matrix of a vector valued
function $\bm{h}:\mathbb{R}^{n}\to\mathbb{R}^{\ell}$ is denoted by
$\nabla\bm{h}\in\mathbb{R}^{\ell\times n}$. Given a vector $\bm{x}$, we will
write $\mathcal{A}(\bm{x})$ denote the set of active constraints at $\bm{x}$,
that is
$\mathcal{A}(\bm{x}):=\\{i\in\mathbb{N}:\bm{a}_{i}^{\intercal}\bm{x}=0\\}$
where $\bm{a}_{i}^{\intercal}$ is the $i$th row of matrix $\bm{A}$. For an
integer $j$, we will write $f\in\mathcal{C}^{j}$ to denote that the function
$f$ is $j$ times continuously differentiable. We write $c$ to denote a generic
nonnegative constant that takes on different values in different inequalities.
Given an interval $[a,b]\subset\mathbb{R}$, we write $Proj_{[a,b]}(\bm{x})$ to
denote the euclidean projection of each component of the vector $\bm{x}$ onto
the interval $[a,b]$. That is, if $\bm{v}=Proj_{[a,b]}(\bm{x})$ then the $j$th
component of $\bm{v}$ is given by

(8) $\displaystyle v_{j}=\left\\{\begin{array}[]{lll}a&:&x_{j}\leq a\\\
x_{j}&:&a<x_{j}<b\\\ b&:&x_{j}\geq b\end{array}\right..$

Additionally, in order to simplify the statement of several results throughout
the paper, we provide abbreviations for assumptions that are used. Note that
each assumption will be clearly referenced when required.

  1. (LICQ)

Linear Independence Constraint Qualification: Given $\bm{x}$,
$\begin{bmatrix}\nabla\bm{h}(\bm{x})\\\
\bm{A}_{\mathcal{A}(\bm{x})}\end{bmatrix}$ has full row rank.

  2. (SCS)

Strict complementary slackness holds at a minimizer of problem (3). That is,
Definition A.3 holds at a minimizer of problem (3).

##  2 Error Estimators for Nonlinear Optimization

In [18], the authors presented two error estimators for optimization problems
of the form

(11)
$\displaystyle\begin{array}[]{cc}\displaystyle\min_{\bm{x}\in\mathbb{R}^{n}}&f(\bm{x})\\\
\text{s.t.}&\bm{h}(\bm{x})=\bm{0},\ \bm{x}\geq\bm{0}.\end{array}$

As a natural generalization to the work in [18], we provide minor updates to
these error estimators for use with problem (3). Both error estimators make
use of the Lagrangian function
$\mathcal{L}:\mathbb{R}^{n}\times\mathbb{R}^{\ell}\times\mathbb{R}^{m}\to\mathbb{R}$
for problem (3) defined by

(12)
$\displaystyle\mathcal{L}(\bm{x},\bm{\lambda},\bm{\mu})=f(\bm{x})+\bm{\lambda}^{\intercal}\bm{h}(\bm{x})+\bm{\mu}^{\intercal}\bm{r}(\bm{x}).$

The first error estimator is defined over the set
$\mathcal{D}_{0}:=\\{(\bm{x},\bm{\lambda},\bm{\mu}):\bm{x}\in\Omega,\bm{\lambda}\in\mathbb{R}^{\ell},\bm{\mu}\geq\bm{0}\\}$
and denoted $E_{0}:\mathcal{D}_{0}\to\mathbb{R}$ by

(13) $\displaystyle
E_{0}(\bm{x},\bm{\lambda},\bm{\mu})=\sqrt{\|\nabla_{x}\mathcal{L}(\bm{x},\bm{\lambda},\bm{\mu})\|^{2}+\|\bm{h}(\bm{x})\|^{2}-\bm{\mu}^{\intercal}\bm{r}(\bm{x})}.$

We note that $E_{0}$ here is a modified version of $E_{0}$ as it was
originally defined in [18]. In particular, if we were to generalize $E_{0}$
for problem (3) as defined in [18] then the term
$-\bm{\mu}^{\intercal}\bm{r}(\bm{x})$ would be squared in (13). This
modification to $E_{0}$ results in beneficial properties that can be exploited
in the global and local convergence analysis of NPASA. Additionally, in
Corollary 2.3 we will show that $E_{0}$ as defined in (13) satisfies an error
bound of the same form as $E_{1}$ regardless of whether or not strict
complementary slackness holds which was previously required for $E_{0}$ as
defined in [18].

A second error estimator is considered that removes the restricted domain
required by $E_{0}$. This estimator makes use of the componentwise minimum
function, denoted here by
$\bm{\Phi}:\mathbb{R}^{m}\times\mathbb{R}^{m}\to\mathbb{R}^{m}$, so that the
$i$th component is defined by

(14) $\displaystyle\Phi_{i}\left(\bm{x},\bm{y}\right)=\min\\{x_{i},y_{i}\\},$

for $1\leq i\leq m$. Then the second error estimator
$E_{1}:\mathbb{R}^{n}\times\mathbb{R}^{\ell}\times\mathbb{R}^{m}\to\mathbb{R}$
is defined by

(15) $\displaystyle
E_{1}(\bm{x},\bm{\lambda},\bm{\mu})=\sqrt{\|\nabla_{x}\mathcal{L}(\bm{x},\bm{\lambda},\bm{\mu})\|^{2}+\|\bm{h}(\bm{x})\|^{2}+\|\bm{\Phi}\left(-\bm{r}(\bm{x}),\bm{\mu}\right)\|^{2}}.$

The error bound property for $E_{1}$ is provided in Theorem 2.2.

In our discussion of NPASA, it will be useful to split the error estimators
$E_{0}$ and $E_{1}$ into two parts. To this end, we define the _multiplier_
portion of the error estimators $E_{0}$ and $E_{1}$ by

(16) $\displaystyle
E_{m,0}(\bm{x},\bm{\lambda},\bm{\mu}):=\|\nabla_{x}\mathcal{L}(\bm{x},\bm{\lambda},\bm{\mu})\|^{2}-\bm{\mu}^{\intercal}\bm{r}(\bm{x})$

and

(17) $\displaystyle
E_{m,1}(\bm{x},\bm{\lambda},\bm{\mu}):=\|\nabla_{x}\mathcal{L}(\bm{x},\bm{\lambda},\bm{\mu})\|^{2}+\|\bm{\Phi}\left(-\bm{r}(\bm{x}),\bm{\mu}\right)\|^{2},$

respectively, and the _constraint_ portion of the error estimators by

(18) $\displaystyle E_{c}(\bm{x}):=\|\bm{h}(\bm{x})\|^{2}.$

Then from the definitions of $E_{0}$ and $E_{1}$ in (13) and (15) it follows
that

(19) $\displaystyle
E_{j}(\bm{x},\bm{\lambda},\bm{\mu})^{2}=E_{m,j}(\bm{x},\bm{\lambda},\bm{\mu})+E_{c}(\bm{x}),$

for $j\in\\{0,1\\}$.

As will be evident throughout the analysis in the following sections, there
are practical benefits and drawbacks to each error estimator provided in this
section. We briefly highlight some of the key differences between the error
estimators to provide some insight on when each estimator is most useful. The
first difference illustrated is that $E_{1}$ provides a tighter bound than
$E_{0}$. This is established in the following Lemma. The proof is omitted as
it is straightforward.

######  Lemma 2.1.

Suppose
$(\bm{x},\bm{\lambda},\bm{\mu})\in\mathcal{D}_{0}:=\\{(\bm{x},\bm{\lambda},\bm{\mu}):\bm{x}\in\Omega,\bm{\lambda}\in\mathbb{R}^{\ell},\bm{\mu}\geq\bm{0}\\}$.
Then

(20) $\displaystyle E_{m,1}(\bm{x},\bm{\lambda},\bm{\mu})\leq
E_{m,0}(\bm{x},\bm{\lambda},\bm{\mu}).$

and

(21) $\displaystyle E_{1}(\bm{x},\bm{\lambda},\bm{\mu})\leq
E_{0}(\bm{x},\bm{\lambda},\bm{\mu}).$

We now provide the error bound properties for $E_{0}$ and $E_{1}$. We note
that the theorem and proof only require slight modifications to the statement
and proof of Theorem 3.1 in [18]. As such, we omit the proof. The error bound
result for the modified $E_{0}$ then immediately follows as a corollary to
Theorem 2.2 when Lemma 2.1 is applied.

######  Theorem 2.2.

Suppose that $\bm{x}^{*}$ is a local minimizer of problem (3) and that
$f,\bm{h}\in\mathcal{C}^{2}$ at $\bm{x}^{*}$. If there exists
$\bm{\lambda}^{*}$ and $\bm{\mu}^{*}\geq\bm{0}$ such that the KKT conditions
and the second-order sufficient optimality conditions hold satisfied at
$(\bm{x}^{*},\bm{\lambda}^{*},\bm{\mu}^{*})$ then there exists a neighborhood
$\mathcal{N}$ of $\bm{x}^{*}$ and a constant $c$ such that

(22)
$\displaystyle\|\bm{x}-\bm{x}^{*}\|+\|\bm{\lambda}-\hat{\bm{\lambda}}\|+\|\bm{\mu}-\hat{\bm{\mu}}\|\leq
cE_{1}(\bm{x},\bm{\lambda},\bm{\mu})$

for all $\bm{x}\in\mathcal{N}$ where $(\hat{\bm{\lambda}},\hat{\bm{\mu}})$
denotes the projection of $(\bm{\lambda},\bm{\mu})$ onto the set of KKT
multipliers at $\bm{x}^{*}$.

######  Corollary 2.3.

Suppose that the hypotheses of Theorem 2.2 are satisfied at
$(\bm{x}^{*},\bm{\lambda}^{*},\bm{\mu}^{*})$. Then there exists a neighborhood
$\mathcal{N}$ of $\bm{x}^{*}$ and a constant $c$ such that

(23)
$\displaystyle\|\bm{x}-\bm{x}^{*}\|+\|\bm{\lambda}-\hat{\bm{\lambda}}\|+\|\bm{\mu}-\hat{\bm{\mu}}\|\leq
cE_{0}(\bm{x},\bm{\lambda},\bm{\mu})$

for all $(\bm{x},\bm{\lambda},\bm{\mu})\in\mathcal{D}_{0}$ with
$\bm{x}\in\mathcal{N}$ where $(\hat{\bm{\lambda}},\hat{\bm{\mu}})$ denotes the
projection of $(\bm{\lambda},\bm{\mu})$ onto the set of KKT multipliers at
$\bm{x}^{*}$.

Based on the error bounds in (22) and (23) it follows that $E_{0}$ and $E_{1}$
are useful for determining when to terminate an iterative method for locating
a stationary point of problem (3). In particular, criterion for stopping an
iterative method for solving problem (3) could be when either
$E_{0}(\bm{x},\bm{\lambda},\bm{\mu})\leq\varepsilon$ or
$E_{1}(\bm{x},\bm{\lambda},\bm{\mu})\leq\varepsilon$ for some small constant
$\varepsilon$. A practical take-away from Lemma 2.1 is that $E_{1}$ is a
better choice than $E_{0}$ for stopping criterion in NPASA.

We now consider a setting in which $E_{0}$ is more beneficial than $E_{1}$.
Based on the bounds (22) and (23) established in Theorem 2.2 and Corollary
2.3, it is of interest of minimize either $E_{0}$ or $E_{1}$ subject to the
constraints $\bm{x}\in\Omega$ and $\bm{\mu}\geq\bm{0}$. In particular, for a
fixed vector $\bm{\bar{x}}\in\Omega$ near a stationary point it is of interest
to solve the problem

(26)
$\displaystyle\begin{array}[]{cc}\displaystyle\min_{\bm{\lambda},\bm{\mu}}&E_{j}(\bm{\bar{x}},\bm{\lambda},\bm{\mu})^{2}\\\
\text{s.t.}&\bm{\mu}\geq\bm{0},\end{array}$

for some $j\in\\{0,1\\}$. In such a setting, we claim that it may be easier to
choose $j=0$ instead of $j=1$. To support this claim, suppose that
$\mathcal{B}$ is a neighborhood of $\bm{\bar{x}}$ such that
$f,\bm{h}\in\mathcal{C}^{t}(\mathcal{B})$, for some positive integer $t$. By
definition of $E_{0}$, it follows immediately that
$E_{0}(\bm{x},\bm{\lambda},\bm{\mu})^{2}\in\mathcal{C}^{t}(\mathcal{B})$. As
such, differentiable methods can be used to solve (26) when $j=0$. However,
due to the nature of the componentwise minimum function in the definition of
$E_{1}$ in (15), it is not necessarily true that
$E_{1}(\bm{x},\bm{\lambda},\bm{\mu})^{2}\in\mathcal{C}^{t}(\mathcal{B})$.
Hence, alternative techniques, such as the active set algorithm in [18], would
be required to solve (26) when $j=1$. Additionally, differentiability of the
error estimator is useful in establishing desirable local convergence results
in the companion paper [11]. Ultimately, our method will leverage the
strengths of each error estimator for solving problem (3). We conclude this
comparison by providing conditions under which $E_{1}$ is differentiable in
the following result. The proof is omitted as it is straightforward.

######  Lemma 2.4.

Let $(\bm{x}^{*},\bm{\lambda}^{*},\bm{\mu}^{*})$ be a KKT point for problem
(3) that satisfies assumption (SCS) and suppose that
$f,\bm{h}\in\mathcal{C}^{2}$ in a neighborhood of $\bm{x}^{*}$. Then there
exists a neighborhood $\mathcal{N}$ about
$(\bm{x}^{*},\bm{\lambda}^{*},\bm{\mu}^{*})$ such that
$E_{1}(\bm{x},\bm{\lambda},\bm{\mu})\in\mathcal{C}^{2}$, for all
$(\bm{x},\bm{\lambda},\bm{\mu})\in\mathcal{N}$.

##  3 NPASA: Algorithm for Nonlinear Programming Problems

In this section, we present details of the Nonlinear Polyhedral Active Set
Algorithm (NPASA) designed to solve problem (3). Together with the pseudocode
for NPASA, we provide a high-level discussion that serves to justify and
motivate our approach.

###  3.1 NPASA Overview and Pseudocode

NPASA is organized into two phases: Phase one contains an algorithm designed
to ensure global convergence of NPASA, the Global Step (GS) algorithm, while
phase two contains an algorithm designed to achieve fast local convergence by
performing a balanced reduction of the constraint and multiplier error
estimators, the Local Step (LS) algorithm. NPASA improves upon a dual
algorithm established by Hager [17] through the use of enhanced error
estimators, a cleverly designed phase two that can provide quadratic local
convergence under reasonable assumptions, and by leveraging advances for
solving polyhedral constrained optimization problems [20, 21]. Note that we
will write $(\bm{x}_{k},\bm{\lambda}_{k},\bm{\mu}_{k})$ to denote the current
primal-dual iterate of NPASA.

1 Inputs: Initial guess $(\bm{x},\bm{\lambda},\bm{\mu})$ and scalar parameters
$\phi>1$, $\bar{\lambda}>0$, and $q$

2
$\displaystyle\bm{\bar{\lambda}}=Proj_{[-\bar{\lambda},\bar{\lambda}]}\left(\bm{\lambda}\right)$

3
$\displaystyle\bm{x}^{\prime}=\arg\min\left\\{\mathcal{L}_{q}\left(\bm{x},\bm{\bar{\lambda}}\right):\bm{x}\in\Omega\right\\}$

4 Set
$\bm{\lambda}^{\prime}\leftarrow\bm{\bar{\lambda}}+2q\bm{h}(\bm{x}^{\prime})$

5 Construct $\bm{\mu}(\bm{x}^{\prime},1)$ from PPROJ output111Using formula
(B.20) in companion paper [11] and set
$\bm{\mu}^{\prime}\leftarrow\bm{\mu}(\bm{x}^{\prime},1)$

Return $(\bm{x}^{\prime},\bm{\lambda}^{\prime},\bm{\mu}^{\prime})$

Algorithm 1 GS - Global Step Algorithm

Global Step (GS). The goal here is to compute a new iterate
$(\bm{x}_{k+1},\bm{\lambda}_{k+1},\bm{\mu}_{k+1})$ such that
$E_{m,0}(\bm{x}_{k+1},\bm{\lambda}_{k+1},\bm{\mu}_{k+1})\leq\theta
E_{c}(\bm{x}_{k+1})$. This is done by using a method of multipliers approach
together with PASA [20]. Specifically, we use $\bm{x}_{k}$ as a starting guess
in an iterative method applied to the polyhedral constrained augmented
Lagrangian problem

(29)
$\displaystyle\begin{array}[]{cc}\min&\mathcal{L}_{q}(\bm{x},\bm{\bar{\lambda}}_{k})\\\
\text{s.t.}&\bm{x}\in\Omega\end{array}$

where the objective function in problem (29) is the augmented Lagrangian
function given by

(30)
$\displaystyle\mathcal{L}_{q}(\bm{x},\bm{\lambda})=f(\bm{x})+\bm{\lambda}^{\intercal}\bm{h}(\bm{x})+q\|\bm{h}(\bm{x})\|^{2},$

$\bm{\bar{\lambda}}_{k}=Proj_{[-\bar{\lambda},\bar{\lambda}]}\left(\bm{\lambda}_{k}\right)$,
$\bar{\lambda}>0$ is a safeguarding parameter, and $q\in\mathbb{R}$ is a
penalty parameter. At a local minimizer of problem (29), there are easily
computable multipliers for which $E_{m,0}$ vanishes which we use to update
$\bm{\lambda}_{k+1}$ and $\bm{\mu}_{k+1}$. To encourage desired convergence
properties, we use a technique often referred to as modified augmented
Lagarangian with safeguarding. This approach has been useful in establishing
convergence in several implementations of augmented Lagrangian solvers [4, 26,
38].

1 Inputs: Initial guess $(\bm{x},\bm{\lambda},\bm{\mu})$ and scalar parameters
$\theta\in(0,1)$, $\alpha\in(0,1]$, $\beta\geq 1$, $\sigma\in(0,1)$,
$\tau\in(0,1)$, $p>>1$, $\delta\in(0,1)$, $\gamma\in(0,1)$

2 _Constraint Step_ : Set $\bm{w}_{0}\leftarrow\bm{x}$ and set $i\leftarrow
0$.

3 while  _$E_{c}(\bm{w}_{i}) >\theta E_{m,1}(\bm{x},\bm{\lambda},\bm{\mu})$_
do

4 Choose $p_{i}$ such that
$p_{i}\geq\max\left\\{\beta^{2},\|\bm{h}(\bm{w}_{i})\|^{-2}\right\\}$ and set
$s_{i}=1$

5
$\left[\bm{\overline{w}}_{i+1},\bm{y}_{i+1}\right]=\arg\min\left\\{\|\bm{w}-\bm{w}_{i}\|^{2}+p_{i}\|\bm{y}\|^{2}:\nabla\bm{h}(\bm{w}_{i})(\bm{w}-\bm{w}_{i})+\bm{y}=-\bm{h}(\bm{w}_{i}),\bm{w}\in\Omega\right\\}$

6 Set $\alpha_{i+1}\leftarrow 1-\|\bm{y}_{i+1}\|$

7 if  _$\alpha_{i+1} <\alpha$_ then

8 Return $(\bm{x},\bm{\lambda},\bm{\mu})$

9 else

10 while  _$\|\bm{h}(\bm{w}_{i}+s_{i}(\bm{\overline{w}}_{i+1}-\bm{w}_{i}))\|
>(1-\tau\alpha_{i+1}s_{i})\|\bm{h}(\bm{w}_{i})\|$_ do

11 Set $s_{i}\leftarrow\sigma s_{i}$

12 Set
$\bm{w}_{i+1}\leftarrow\bm{w}_{i}+s_{i}(\bm{\overline{w}}_{i+1}-\bm{w}_{i})$
and set $i\leftarrow i+1$

13

14Set $\bm{w}\leftarrow\bm{w}_{i}$

15 _Multiplier Step_ : Set $\bm{z}_{0}\leftarrow\bm{w}$, set $p_{0}\leftarrow
p$, and set $i\leftarrow 0$

16
$(\bm{\nu}_{0},\bm{\eta}_{0})\in\arg\min\left\\{E_{m,0}(\bm{z}_{0},\bm{\nu},\bm{\eta})+\gamma\|[\bm{\nu},\bm{\eta}]\|^{2}:\bm{\eta}\geq\bm{0}\right\\}$

17
$\bm{\eta}_{0}^{\prime}\in\arg\min\left\\{E_{m,1}(\bm{z}_{0},\bm{\nu}_{0},\bm{\eta}):\bm{\eta}\geq\bm{0}\right\\}$

18 while  _$E_{m,1}(\bm{z}_{i},\bm{\nu}_{i},\bm{\eta}_{i}^{\prime}) >\theta
E_{c}(\bm{w})$_ do

19 Increase $p_{i}$ if necessary

20
$\bm{z}_{i+1}=\arg\min\left\\{\mathcal{L}_{p_{i}}^{i}(\bm{z},\bm{\nu}_{i}):\nabla\bm{h}(\bm{z}_{i})(\bm{z}-\bm{z}_{i})=0,\bm{z}\in\Omega\right\\}$

21
$(\bm{\nu}_{i+1},\bm{\eta}_{i+1})\in\arg\min\left\\{E_{m,0}(\bm{z}_{i+1},\bm{\nu},\bm{\eta})+\gamma\|[\bm{\nu},\bm{\eta}]\|^{2}:\bm{\eta}\geq\bm{0}\right\\}$

22
$\bm{\eta}_{i+1}^{\prime}\in\arg\min\left\\{E_{m,1}(\bm{z}_{i+1},\bm{\nu}_{i+1},\bm{\eta}):\bm{\eta}\geq\bm{0}\right\\}$

23 if  _$E_{m,1}(\bm{z}_{i+1},\bm{\nu}_{i+1},\bm{\eta}_{i+1}^{\prime}) >\delta
E_{m,1}(\bm{z}_{i},\bm{\nu}_{i},\bm{\eta}_{i}^{\prime})$_ then

24 Return $(\bm{x},\bm{\lambda},\bm{\mu})$

25 else

26 Set $i\leftarrow i+1$

27

Return $(\bm{z}_{i},\bm{\nu}_{i},\bm{\eta}_{i}^{\prime})$

Algorithm 2 LS - Local Step Algorithm

Local Step (LS). LS, Algorithm 2, consists of two subproblems formulated to
achieve a balanced reduction of the constraint and multiplier error for
problem (3):

  1. 1.

Constraint Step. Use iterative method to generate $\bm{w}$ satisfying
$E_{c}(\bm{w})\leq\theta E_{m,1}(\bm{x}_{k},\bm{\lambda}_{k},\bm{\mu}_{k})$.

  2. 2.

Multiplier Step. Starting from $\bm{w}$, apply an iterative method to generate
$(\bm{z},\bm{\nu},\bm{\eta})$ which satisfies
$E_{m,1}(\bm{z},\bm{\nu},\bm{\eta})\leq\theta E_{c}(\bm{w})$.

To solve the constraint step we use a perturbed Newton’s method scheme. With a
unit step size, Newton’s method is given by the iterative scheme

(31)
$\displaystyle\bm{w}_{i+1}=\arg\min\left\\{\|\bm{w}-\bm{w}_{i}\|^{2}:\nabla\bm{h}(\bm{w}_{i})(\bm{w}-\bm{w}_{i})=-\bm{h}(\bm{w}_{i}),\bm{w}\geq\bm{0}\right\\},$

where we start with $\bm{w}_{0}=\bm{x}_{k}$. Problem (31) can be solved using
a polyhedral projection algorithm such as PPROJ [21]. Based on the analysis in
[17], the procedure in step 1 is locally quadratically convergent when the
second-order sufficient optimality conditions hold and the active constraint
gradients are linearly independent. However, as the constraint set in (31) may
be infeasible, we encourage robustness of this approach by introducing a
perturbation to the linear constraint. The modified problem becomes

(32)
$\displaystyle\left[\bm{w}_{i+1},\bm{y}_{i+1}\right]=\arg\min\left\\{\|\bm{w}-\bm{w}_{i}\|^{2}+p\|\bm{y}\|^{2}:\nabla\bm{h}(\bm{w}_{i})(\bm{w}-\bm{w}_{i})+\bm{y}=-\bm{h}(\bm{w}_{i}),\bm{w}\geq\bm{0}\right\\},$

where we start with $\bm{w}_{0}=\bm{x}_{k}$ and $p>0$ is a penalty parameter.
Under this formulation, this problem can still be solved using the algorithm
PPROJ [21] and we consider this method in detail in Subsection 5.1.

More details on the constraint step can be found in Sections 5.1 and 5.2. The
main problem can be solved using the PPROJ [21] algorithm followed by a line
search. The remaining lines of the constraint step are to ensure desired local
convergence properties are satisfied. In particular, the requirements on
$p_{i}$ and $\alpha_{i}$ are covered in detail during the convergence analysis
in Section 5. As a note to the reader, when $\alpha_{i+1}<\alpha$ then
Algorithm 2 results in phase two immediately branching to phase one in NPASA.
This branching is triggered when the perturbation in the constraint step
becomes too large thereby indicating that the constraint step may be
experiencing trouble minimizing the constraint error, $E_{c}$. In such a case,
NPASA works to find a better point by executing GS then branching back to LS.

Step 2, the multiplier step, can be decomposed into three parts:

  * 2a.

Dual Step for equality multiplier:
$\displaystyle(\bm{\nu}_{i},\bm{\eta}_{i})\in\arg\min_{\bm{\lambda},\bm{\eta}}\left\\{E_{m,0}(\bm{z}_{i},\bm{\nu},\bm{\eta})+\gamma\|[\bm{\nu},\bm{\eta}]\|^{2}:\bm{\eta}\geq\bm{0}\right\\}$

  * 2b.

Dual Step for inequality multiplier:
$\displaystyle\bm{\eta}_{i}^{\prime}\in\arg\min_{\bm{\eta}}\left\\{E_{m,1}(\bm{z}_{i},\bm{\nu}_{i},\bm{\eta}):\bm{\eta}\geq\bm{0}\right\\}$

  * 2c.

Primal Step:
$\displaystyle\bm{z}_{i+1}=\arg\min_{\bm{z}}\left\\{\mathcal{L}_{p_{i}}^{i}(\bm{z},\bm{\nu}_{i}):\nabla\bm{h}(\bm{z}_{i})(\bm{z}-\bm{z}_{i})=\bm{0},\bm{z}\in\Omega\right\\}$,

where Step 2c makes use of the function defined by

(33)
$\displaystyle\mathcal{L}_{p}^{i}(\bm{z},\bm{\nu}):=f(\bm{z})+\bm{\nu}^{\intercal}\bm{h}(\bm{z})+p\|\bm{h}(\bm{z})-\bm{h}(\bm{z}_{i})\|^{2}.$

For initialization, $\bm{z}_{0}$ is set equal to the final $\bm{w}$ generated
in the Constraint Step.

Noting that $E_{m,0}$ is twice continuously differentiable with respect to
$\bm{\lambda}$ and $\bm{\mu}$ and the hessian is positive semi-definite, there
are many techniques available to solve Step 2a. The regularization term is
added to ensure positive definiteness of the hessian which is used in
establishing fast local convergence [11]. We note that PASA [20] can be used
for Step 2a. The choice of $E_{m,0}$ in Step 2a over $E_{m,1}$ is required to
establish desirable local convergence properties of NPASA for nondegenerate
and degenerate minimizers of problem (3). More details can be found in the
local convergence analysis of the multiplier step in the companion paper [11].

The update to the inequality multiplier in Step 2b is required to achieve
local quadratic convergence of the multiplier error
$E_{m,1}(\bm{z}_{i},\bm{\nu}_{i},\bm{\eta}_{i}^{\prime})$. While it seems that
Step 2a should be able to minimize over $\bm{\nu}$ by fixing
$\bm{\eta}=\bm{\eta}_{i-1}^{\prime}$, there is difficulty in establishing a
stability result for the solution in Step 2b that is required for local
convergence analysis when assumption (SCS) fails to hold. Hence, the Step 2a
and Step 2b structure is designed to ensure that NPASA still satisfies fast
local convergence when approaching both degenerate and nondegenerate
minimizers of (3).

Step 2c is a polyhedral constrained optimization problem which can also be
solved using PASA [20]. It should be noted that in Step 2c the polyhedron will
always contain a feasible point since the point $\bm{z}=\bm{z}_{i}$ satisfies
all the constraints. Hence, no reformulation of the problem in Step 2c with
some perturbation of the equality constraint is required to ensure robustness.
A result on the global convergence properties of the multiplier step is found
in Section 5.3.

NPASA. As previously noted, NPASA, Algorithm 3, is split into two phases and
contains criterion for branching between the phases that encourage desirable
global and local convergence properties. Phase one contains the GS algorithm,
updates for the penalty parameter used in GS, and criterion for branching to
phase two. Phase two contains the LS algorithm and criterion for branching to
phase one. The goal of the branching criterion is to remain in phase two
whenever possible as fast local convergence is achieved there but to still
ensure global convergence by branching to phase one when sufficient reduction
of the error estimator is not occurring in phase two. More details on how
these branching conditions achieve this goal can be found in Section 6.

If the starting guess is good enough then NPASA can perform repeated
iterations of the LS algorithm without branching to phase one. We continue in
this manner until the error estimator $E_{1}$ is below the stopping tolerance,
$\varepsilon$, at which point NPASA terminates and returns the final iterate.
However, if the starting guess is not good or if the Newton iteration in the
constraint step does not converge quickly then we branch to phase one of
Algorithm 3 and run the GS algorithm. Note that upon entering phase one we
increase the penalty parameter $q_{k}$. As $q_{k}$ increases, $E_{c}$
typically decreases since the constraint is penalized in the objective
function while $E_{m}$ decreases as iterates approach a stationary point of
(29).

1 Inputs: Initial guess $(\bm{x}_{0},\bm{\lambda}_{0},\bm{\mu}_{0})$ and
scalar parameters $\varepsilon>0$, $\theta\in(0,1)$, $\phi>1$,
$\bar{\lambda}>0$, $q_{0}\geq 1$, $\alpha\in(0,1]$, $\beta\geq 1$,
$\sigma\in(0,1)$, $\tau\in(0,1)$, $p>>1$, $\delta\in(0,1)$, $\gamma>0$.

2Set $e_{0}=E_{1}(\bm{x}_{0},\bm{\lambda}_{0},\bm{\mu}_{0})$, $k=0$, and goto
phase one.

3 —— Phase one ——

4 Set $q_{k}\leftarrow\max\left\\{\phi,(e_{k-1})^{-1}\right\\}q_{k-1}$.

5 while  _$E_{m,1}(\bm{x}_{k},\bm{\lambda}_{k},\bm{\mu}(\bm{x}_{k},1))
>\varepsilon$_ do

6 _Global Step_ :
$\displaystyle\left(\bm{x}_{k+1},\bm{\lambda}_{k+1},\bm{\mu}_{k+1}\right)\leftarrow
GS\left(\bm{x}_{k},\bm{\lambda}_{k},\bm{\mu}_{k};\phi,\bar{\lambda},q_{k}\right)$

7 _Update parameters_ : $q_{k+1}\leftarrow q_{k}$,
$e_{k+1}\leftarrow\min\left\\{E_{1}(\bm{x}_{k+1},\bm{\lambda}_{k+1},\bm{\mu}_{k+1}),e_{k}\right\\}$,
$k\leftarrow k+1$.

8 _Check branching criterion_ :

9 if  _$E_{m,1}(\bm{x}_{k},\bm{\lambda}_{k},\bm{\mu}_{k})\leq\theta
E_{c}(\bm{x}_{k-1})$_ then

10 goto phase two

11

12—— Phase two ——

13 while  _$E_{m,1}(\bm{x}_{k},\bm{\lambda}_{k},\bm{\mu}_{k}) >\varepsilon$_
do

14 _Local Step_ :
$\displaystyle\left(\bm{z},\bm{\nu},\bm{\eta}\right)\leftarrow
LS\left(\bm{x}_{k},\bm{\lambda}_{k},\bm{\mu}_{k};\theta,\alpha,\beta,\sigma,\tau,p,\gamma,\delta\right)$

15 _Check branching criterion_ :

16 if  _$E_{1}(\bm{z},\bm{\nu},\bm{\eta}) >\theta
E_{1}(\bm{x}_{k},\bm{\lambda}_{k},\bm{\mu}_{k})$_ then

17 goto phase one

18 else

19 Set
$(\bm{x}_{k+1},\bm{\lambda}_{k+1},\bm{\mu}_{k+1})\leftarrow(\bm{z},\bm{\nu},\bm{\eta})$
and set
$e_{k+1}\leftarrow\min\left\\{E_{1}(\bm{x}_{k+1},\bm{\lambda}_{k+1},\bm{\mu}_{k+1}),e_{k}\right\\}$

20 Set $k\leftarrow k+1$.

21

Return $(\bm{x}_{k},\bm{\lambda}_{k},\bm{\mu}_{k})$

Algorithm 3 NPASA - Nonlinear Polyhedral Active Set Algorithm

Algorithm 3 contains several parameters which we briefly outline here. First,
we discuss $\varepsilon$ and $\theta$ as they are used in both phases.
$\varepsilon$ is the stopping tolerance used in NPASA and $\theta$ is a
scaling parameter used when checking if NPASA should branch to another phase.
The choice of $\theta\in(0,1)$ ensures that NPASA is decreasing the error
estimator a sufficient amount at each iteration and is used to establish
global convergence of NPASA. Phase one requires the parameters $\phi$,
$\bar{\lambda}$, and $q_{0}$. The parameter $\phi$ is used to increase the
penalty $q_{k}$ each time NPASA enters phase one and $q_{0}$ is the initial
value for the penalty parameter in the augmented Lagrangian. A typical choice
would be $\phi=10$ so that the penalty increases by a factor of 10 each time
NPASA branches to phase one. $\bar{\lambda}$ is used to safeguard the method
of multipliers approach used in the GS algorithm in phase one. The LS
algorithm in phase two requires the parameters $\beta$, $\alpha$, $\sigma$,
$\tau$, $p$, $\gamma$, and $\delta$. $\beta$ is used to ensure that the
penalty parameter in constraint step does not decrease at each iteration.
$\alpha$ is used to determine how much perturbation is permitted at a solution
of the constraint step. If $\alpha=1$ then no perturbation is allowed. As
$\alpha$ approaches zero, larger amounts of perturbation are permitted. More
details on the role of $\alpha$ can be found in Lemma 5.3. The parameter
$\tau$ is used during the line search in the constraint step. $\sigma$ is used
to decrease the step length during the line search in the constraint step. A
typical choice is $\sigma=1/2$. $p$ is used to reset the penalty parameter at
the start of the multiplier step in phase two. Updates to $p_{i}$ in the
constraint and multiplier steps will be handled determinstically based on
convergence analysis so that the user only needs to provide an initial value
for $p$. $\gamma$ is used to ensure that the hessian of the problem in step 2a
is positive definite. Lastly, $\delta$ is used to determine if the multiplier
error estimator in the multiplier step of LS has sufficiently decreased at
each iteration.

1 Inputs: Initial guess $\bm{x}_{k}$, scalars $\theta\in(0,1)$,
$\gamma\in(0,1)$, $\varepsilon>0$, $q_{k}>0$ and vector parameter
$\bm{\bar{\lambda}}_{k}$

2 Set $\bm{u}_{0}=\bm{x}_{k}$, $\bm{u}_{1}=\mathcal{P}_{\Omega}(\bm{u}_{0})$,
and $i=1$

3 —— Phase one ——

4 while
_$E_{m,0}(\bm{u}_{i},\bm{\bar{\lambda}}_{k}+2q_{k}\bm{h}(\bm{u}_{i}),\bm{\mu}(\bm{u}_{i},1))
>\varepsilon$_ do

5 Set $\bm{u}_{i+1}\leftarrow GPA(\bm{u}_{i})$ and $i\leftarrow i+1$

6 if  _$\mathcal{U}(\bm{u}_{i})=\emptyset$ and
$e_{\scaleto{PASA}{4pt}}(\bm{u}_{i})<\theta
E_{m,0}(\bm{u}_{i},\bm{\bar{\lambda}}_{k}+2q_{k}\bm{h}(\bm{u}_{i}),\bm{\mu}(\bm{u}_{i},1))$_
then

7 Set $\theta\leftarrow\gamma\theta$

8 if  _$e_{\scaleto{PASA}{4pt}}(\bm{u}_{i})\geq\theta
E_{m,0}(\bm{u}_{i},\bm{\bar{\lambda}}_{k}+2q_{k}\bm{h}(\bm{u}_{i}),\bm{\mu}(\bm{u}_{i},1))$_
then

9 goto phase two

10

11—— Phase two ——

12 while
_$E_{m,0}(\bm{u}_{i},\bm{\bar{\lambda}}_{k}+2q_{k}\bm{h}(\bm{u}_{i}),\bm{\mu}(\bm{u}_{i},1))
>\varepsilon$_ do

13 Set $\bm{u}_{i+1}\leftarrow LCO(\bm{u}_{i})$ and $i\leftarrow i+1$

14 if  _$e_{\scaleto{PASA}{4pt}}(\bm{u}_{i}) <\theta
E_{m,0}(\bm{u}_{i},\bm{\bar{\lambda}}_{k}+2q_{k}\bm{h}(\bm{u}_{i}),\bm{\mu}(\bm{u}_{i},1))$_
then

15 goto phase one

16

Return $\bm{u}_{i}$

Algorithm 4 Modified PASA Algorithm for GS Augmented Lagrangian Problem

###  3.2 Modified PASA Algorithm for Augmented Lagrangian Problem in Global
Step

Here, we provide Algorithm 4 which is PASA with modified stopping and
branching criterion designed specifically for the augmented Lagrangian problem
found in line 3 of Algorithm 1. Note that the algorithms GPA and LCO used in
Algorithm 4 refer to a _Gradient Projection Algorithm_ and a _Linearly
Constrained Optimizer_ , respectively. The original paper on PASA [20]
provides pseudocode for GPA and conditions that the LCO must satisfy in order
for PASA to satisfy their convergence results. The details of the GPA and LCO
used in Algorithm 4 are the same as specified in the original PASA paper. The
value $e_{\scaleto{PASA}{4pt}}(\bm{x})$ corresponds to the _local error
estimator for PASA_ as defined in [20].

Traditionally, the stopping criterion in PASA at the iterate $\bm{x}$ is given
by $E_{\scaleto{PASA}{4pt}}(\bm{x})\leq\varepsilon$, where $\varepsilon$ is a
user provided stopping tolerance value and $E_{\scaleto{PASA}{4pt}}(\bm{x})$
denotes the _global error estimator for PASA_ at the primal variable $\bm{x}$
originally defined in [20]. We have modified this criterion by replacing
$E_{\scaleto{PASA}{4pt}}$ with the multiplier error estimator for problem (3),
$E_{m,0}$, and intend to use the formula derived during the following
analysis, in (46), to compute $E_{m,0}$. The motivation for this modified
stopping and branching criterion is to reduce the number of assumptions
required to ensure convergence of Algorithm 1. While global convergence for
Algorithm 1 can be established when using PASA with the original stopping
criterion, the number of assumptions required for convergence are reduced when
using Algorithm 4. A more detailed analysis yielding the modifications can be
found in Section 4.2.

##  4 Convergence Analysis for Phase One of NPASA

We begin our convergence analysis by performing convergence analysis for phase
one of NPASA. Motivation for the global step in phase one is to ensure the
robustness of the NPASA algorithm. In particular, phase two of NPASA may
encounter problems when the linearization of $\bm{h}(\bm{x})$ at the current
iterate requires a large perturbation to be feasible resulting in an
unsuccessful minimization of $\|\bm{h}(\bm{x})\|$ in the Constraint Step. In
such an instance, phase one serves as a safeguard to move the primal iterate
to a new point before attempting phase two again. For phase one we have
decided to use augmented Lagrangian techniques. First introduced by Hestenes
[23] and Powell [31] and sometimes referred to as the method of multipliers,
augmented Lagrangian methods and their convergence properties have since been
well studied [3, 26, 28] and have been successful at solving nonlinear
programs in practice using software packages such as ALGENCAN [2], ALGLIB [5],
LANCELOT [7, 9], and MINOS [27]. Hence, the choice of an augmented Lagrangian
technique is a suitable safeguard to ensure NPASA will converge to a
stationary point in the case that problems arise during the constraint step in
phase two.

This section is organized as follows. First, we analyze the method of
multipliers scheme in Algorithm 1, GS algorithm. In particular, we establish
necessary conditions under which a subsequence of iterates of the GS algorithm
converge to a stationary point for problem (3). During this discussion, we
also develop a closed form expression for computing the multiplier error of
the original nonlinear program, problem (3), when solving the augmented
Lagrangian problem without requiring the explicit value of the inequality
multiplier. In Section 4.2, we provide theoretical analysis for Algorithm 4,
the modified PASA algorithm, which reduces the number of necessary conditions
required to establish convergence when solving the minimization problem in the
GS algorithm using Algorithm 4.

###  4.1 Preliminary Convergence Analysis for GS Algorithm

Consider the minimization problem in Algorithm 1 at the $k$th iteration of
NPASA given by

(36)
$\displaystyle\begin{array}[]{cc}\displaystyle\min_{\bm{x}\in\mathbb{R}^{n}}&\mathcal{L}_{q_{k}}(\bm{x},\bm{\bar{\lambda}}_{k})=f(\bm{x})+\bm{\bar{\lambda}}_{k}^{\intercal}\bm{h}(\bm{x})+q_{k}\|\bm{h}(\bm{x})\|^{2}\\\
\text{s.t.}&\bm{r}(\bm{x})\leq\bm{0}\end{array}$

where $\bm{r}(\bm{x})=\bm{A}\bm{x}-\bm{b}$ is the linear function from problem
(3) and
$\bm{\bar{\lambda}}_{k}=Proj_{[-\bar{\lambda},\bar{\lambda}]}(\bm{\lambda}_{k})$
is known. Now suppose that $\bm{x}_{k+1}$ is a local minimizer for problem
(36). Then by the KKT conditions there exists a $\bm{\mu}_{k+1}\geq\bm{0}$
such that

(37) $\displaystyle\bm{0}$
$\displaystyle=\nabla_{x}\mathcal{L}_{q_{k}}(\bm{x},\bm{\bar{\lambda}}_{k})+\bm{\mu}_{k+1}^{\intercal}\bm{A}=\nabla\mathcal{L}(\bm{x}_{k+1},\bm{\bar{\lambda}}_{k}+2q_{k}\bm{h}(\bm{x}_{k+1}),\bm{\mu}_{k+1})$

and $\bm{\mu}_{k+1}^{\intercal}\bm{r}(\bm{x}_{k+1})=0$. Using the formula
$\bm{\lambda}_{k+1}=\bm{\bar{\lambda}}_{k}+2q_{k}\bm{h}(\bm{x}_{k+1})$ from
line 4 of Algorithm 1 it now follows from the definition of $E_{m,0}$ that
$E_{m,0}(\bm{x}_{k+1},\bm{\lambda}_{k+1},\bm{\mu}_{k+1})=0$. Hence, we have
that the multiplier error is equal to zero at
$(\bm{x}_{k+1},\bm{\lambda}_{k+1},\bm{\mu}_{k+1})$, where $\bm{\mu}_{k+1}$ is
some vector that we do not explicitly know.

Now define $\bm{y}(\bm{x},\alpha)$ to be the minimizer for the problem

(40)
$\displaystyle\begin{array}[]{cc}\displaystyle\min_{\bm{y}}&\displaystyle\frac{1}{2}\|\bm{x}-\alpha\bm{g}(\bm{x})-\bm{y}\|^{2}\\\
\text{s.t.}&\bm{r}(\bm{y})\leq\bm{0}\end{array}$

for $\alpha>0$ and where $\bm{g}(\bm{x})$ is the transpose of the gradient of
the cost function in problem (36). Note that in our case we have

(41) $\displaystyle\bm{g}(\bm{x})^{\intercal}$
$\displaystyle=\nabla_{x}\mathcal{L}_{q_{k}}(\bm{x},\bm{\bar{\lambda}}_{k})=\nabla
f(\bm{x})+\bm{\bar{\lambda}}_{k}^{\intercal}\nabla\bm{h}(\bm{x})+2q_{k}\bm{h}(\bm{x})^{\intercal}\nabla\bm{h}(\bm{x}),$

however we will interchangeably write $\bm{g}(\bm{x})^{\intercal}$ to simplify
notation. Now fix $\bm{x}\in\Omega$ and consider problem (40). By the KKT
conditions, there exists a vector $\bm{\mu}(\bm{x},\alpha)$ such that

(42) $\displaystyle\text{Gradient of Lagrangian equals zero:}\
(\bm{y}(\bm{x},\alpha)-\bm{x})^{\intercal}+\alpha\bm{g}(\bm{x})^{\intercal}+\bm{\mu}(\bm{x},\alpha)^{\intercal}\bm{A}=\bm{0},$
(43) $\displaystyle\text{Satisfies inequality constraints:}\
\bm{r}(\bm{y}(\bm{x},\alpha))\leq\bm{0},$ (44)
$\displaystyle\text{Nonnegativity of inequality multipliers:}\
\bm{\mu}(\bm{x},\alpha)\geq\bm{0},$ (45) $\displaystyle\text{Complementary
slackness:}\ \mu_{i}(\bm{x},\alpha)r_{i}(\bm{y}(\bm{x},\alpha))=0,\
\text{for}\ 1\leq i\leq m.$

Our goal is to use the KKT conditions (42) – (45) for problem (40) to derive a
formula for the multiplier error of our original problem, $E_{m,0}$, in terms
of the current primal iterate, $\bm{x}$. This idea will be useful in
establishing a theoretical convergence result for phase one of NPASA and for
later establishing practical convergence results using PASA and modified PASA.
This is summarized in the following result.

######  Lemma 4.1.

Suppose $\bm{x}\in\Omega$. Then for each $\alpha>0$ there exists an inequality
multiplier $\bm{\mu}(\bm{x},\alpha)$ such that

(46) $\displaystyle
E_{m,0}(\bm{x},\bm{\bar{\lambda}}_{k}+2q_{k}\bm{h}(\bm{x}),\bm{\mu}(\bm{x},\alpha)/\alpha)$
$\displaystyle=-\bm{g}(\bm{x})^{\intercal}(\bm{y}(\bm{x},\alpha)-\bm{x})+\left(\frac{1}{\alpha^{2}}-\frac{1}{\alpha}\right)\|\bm{y}(\bm{x},\alpha)-\bm{x}\|^{2}.$

In particular, for $\alpha=1$ we have

(47) $\displaystyle\|\bm{y}(\bm{x},1)-\bm{x}\|^{2}$ $\displaystyle\leq
E_{m,0}(\bm{x},\bm{\bar{\lambda}}_{k}+2q_{k}\bm{h}(\bm{x}),\bm{\mu}(\bm{x},1))\leq\|\nabla_{x}\mathcal{L}_{q_{k}}(\bm{x},\bm{\bar{\lambda}}_{k})\|\|\bm{y}(\bm{x},1)-\bm{x}\|.$

######  Proof 4.2.

Let $\bm{x}\in\Omega$. By recalling that the multiplier error estimator
$E_{m,0}$ for problem (3) is given by
$E_{m,0}(\bm{x},\bm{\lambda},\bm{\mu})=\|\nabla
f(\bm{x})+\bm{\lambda}^{\intercal}\nabla\bm{h}(\bm{x})+\bm{\mu}^{\intercal}\bm{A}\|^{2}-\bm{\mu}^{\intercal}\bm{r}(\bm{x})$,
it follows that

(48) $\displaystyle
E_{m,0}(\bm{x},\bm{\bar{\lambda}}_{k}+2q_{k}\bm{h}(\bm{x}),\bm{\mu})$
$\displaystyle=\|\bm{g}(\bm{x})^{\intercal}+\bm{\mu}^{\intercal}\bm{A}\|^{2}-\bm{\mu}^{\intercal}\bm{r}(\bm{x}).$

Now note that from (42) we have
$\bm{g}(\bm{x})^{\intercal}+\frac{1}{\alpha}\bm{\mu}(\bm{x},\alpha)^{\intercal}\bm{A}=\frac{1}{\alpha}(\bm{x}-\bm{y}(\bm{x},\alpha))^{\intercal}$.
Taking the norm of both sides and squaring yields

(49)
$\displaystyle\left\|\bm{g}(\bm{x})^{\intercal}+\frac{1}{\alpha}\bm{\mu}(\bm{x},\alpha)^{\intercal}\bm{A}\right\|^{2}=\frac{1}{\alpha^{2}}\|\bm{y}(\bm{x},\alpha)-\bm{x}\|^{2}.$

Additionally, right multiplying both sides of (42) by
$\frac{1}{\alpha}(\bm{y}(\bm{x},\alpha)-\bm{x})$ and using (45) yields

(50)
$\displaystyle-\frac{1}{\alpha}\bm{\mu}(\bm{x},\alpha)^{\intercal}\bm{r}(\bm{x})=-\frac{1}{\alpha}\|\bm{y}(\bm{x},\alpha)-\bm{x}\|^{2}-\bm{g}(\bm{x})^{\intercal}(\bm{y}(\bm{x},\alpha)-\bm{x}).$

Combining (48), (49), and (50) yields (46).

Next, from P6 in Proposition 2.1 of [19] we have that
$-\bm{g}(\bm{x})^{\intercal}(\bm{y}(\bm{x},\alpha)-\bm{x})\geq\frac{1}{\alpha}\|\bm{y}(\bm{x},\alpha)-\bm{x}\|^{2}$,
for all $\bm{x}\in\Omega$. Combining this fact with (46) yields a lower bound
for the multiplier error

(51) $\displaystyle
E_{m,0}(\bm{x},\bm{\bar{\lambda}}_{k}+2q_{k}\bm{h}(\bm{x}),\bm{\mu}(\bm{x},\alpha)/\alpha)\geq\frac{1}{\alpha^{2}}\|\bm{y}(\bm{x},\alpha)-\bm{x}\|^{2}.$

In particular, by taking $\alpha=1$ in (46) and recalling (41) we have

(52) $\displaystyle
E_{m,0}(\bm{x},\bm{\bar{\lambda}}_{k}+2q_{k}\bm{h}(\bm{x}),\bm{\mu}(\bm{x},1))=-\nabla_{x}\mathcal{L}_{q_{k}}(\bm{x},\bm{\bar{\lambda}}_{k})(\bm{y}(\bm{x},1)-\bm{x}).$

Lastly, (51) and (52) together with Cauchy-Schwarz yield (47).

Hence, in (46) we have derived a closed form expression for the multiplier
error estimator for the general nonlinear program, problem (3), in terms of
the current primal iterate, $\bm{x}$, and a known KKT equality multiplier. The
inequality in (47) will be useful in establishing the following convergence
result.

######  Theorem 4.3.

Suppose $\Omega$ is compact and that $f,\bm{h}\in\mathcal{C}^{1}(\Omega)$.
Suppose the following are satisfied:

  1. 1.

$\\{\bm{x}_{k}\\}_{k=0}^{\infty}$ is a sequence with $\bm{x}_{k}\in\Omega$ and
$E_{m,0}(\bm{x}_{k},\bm{\lambda}_{k},\bm{\mu}(\bm{x}_{k},1))\leq\theta
E_{c}(\bm{x}_{k})$, for all $k\geq 0$.

  2. 2.

$\\{\bm{\lambda}_{k}\\}_{k=0}^{\infty}$ is updated using the formula for
$\bm{\lambda}^{\prime}$ in line 4 of Algorithm 1 (GS algorithm).

  3. 3.

$\\{q_{k}\\}_{k=0}^{\infty}$ is a sequence of real numbers such that
$q_{k}\to\infty$ as $k\to\infty$.

Then any subsequence of $\\{\bm{x}_{k}\\}_{k=0}^{\infty}$ with limit point
$\bm{x}^{*}$ such that assumption (LICQ) holds at $\bm{x}^{*}$ satisfies
$E_{1}(\bm{x}^{*},\bm{\lambda}^{*},\bm{\mu}(\bm{x}^{*},1))=E_{0}(\bm{x}^{*},\bm{\lambda}^{*},\bm{\mu}(\bm{x}^{*},1))=0$,
where $\bm{\lambda}^{*}$ is the limit of the corresponding subsequence for
$\\{\bm{\lambda}_{k}\\}$.

######  Proof 4.4.

Let $\\{\bm{x}_{k}\\}_{k=0}^{\infty}$ be a subsequence with limit point
$\bm{x}^{*}$ and let $\\{\bm{\lambda}_{k}\\}_{k=0}^{\infty}$ and
$\\{\bm{\mu}(\bm{x}_{k},1)\\}_{k=0}^{\infty}$ be corresponding subsequences of
multipliers. For simplicity, let $\mathcal{A}\>=\mathcal{A}(\bm{x}^{*})$.
First, we show that $\\{\bm{\lambda}_{k}\\}_{k=0}^{\infty}$ and
$\\{\bm{\mu}(\bm{x}_{k},1)\\}_{k=0}^{\infty}$ are uniformly bounded for $k$
sufficiently large. Rearranging the terms in equation (42) with $\alpha=1$
combined with the fact that
$\mathcal{A}(\bm{x}_{k})\subseteq\mathcal{A}(\bm{x}^{*})$ for sufficiently
large $k$, we have that

(53) $\displaystyle\begin{bmatrix}\nabla\bm{h}(\bm{x}_{k})^{\intercal}\
\Big{|}\
\bm{A}_{\mathcal{A}}^{\intercal}\end{bmatrix}\begin{bmatrix}\bm{\lambda}_{k}\\\
\bm{\mu}_{\mathcal{A}}(\bm{x}_{k},1)\end{bmatrix}=\bm{x}_{k}-\bm{y}(\bm{x}_{k},1)-\nabla
f(\bm{x}_{k})^{\intercal}$

By our hypothesis that (LICQ) holds at $\bm{x}^{*}$, there exists a $N$ such
that the matrix $\begin{bmatrix}\nabla\bm{h}(\bm{x}_{k})^{\intercal}\ |\
\bm{A}_{\mathcal{A}}^{\intercal}\end{bmatrix}$ is of full column rank, for all
$k\geq N$. Hence, there exists a constant $M$ such that

(54) $\displaystyle\left\|\begin{bmatrix}\bm{\lambda}_{k}\\\
\bm{\mu}_{\mathcal{A}}(\bm{x}_{k},1)\end{bmatrix}\right\|\leq
M\|\bm{x}_{k}-\bm{y}(\bm{x}_{k},1)-\nabla f(\bm{x}_{k})^{\intercal}\|,$

for $k$ sufficiently large. From property P6 in Proposition 2.1 of [19] we
have
$\|\bm{x}_{k}-\bm{y}(\bm{x}_{k},1)\|^{2}\leq\nabla_{x}\mathcal{L}_{q}(\bm{x}_{k},\bm{\bar{\lambda}}_{k-1})^{\intercal}(\bm{x}_{k}-\bm{y}(\bm{x}_{k},1))$.
Hence, using the hypothesis that
$E_{m,0}(\bm{x}_{k},\bm{\lambda}_{k},\bm{\mu}(\bm{x}_{k},1))\leq\theta
E_{c}(\bm{x}_{k})$ together with (47) it follows that
$\|\bm{x}_{k}-\bm{y}(\bm{x}_{k},1)\|^{2}\leq
E_{m,0}(\bm{x}_{k},\bm{\lambda}_{k},\bm{\mu}(\bm{x}_{k},1))\leq
E_{c}(\bm{x}_{k})=\|\bm{h}(\bm{x}_{k})\|^{2}$. Combining this with the
inequality in (54) yields

(55) $\displaystyle\left\|\begin{bmatrix}\bm{\lambda}_{k}\\\
\bm{\mu}_{\mathcal{A}}(\bm{x}_{k},1)\end{bmatrix}\right\|\leq
M\left(\|\bm{h}(\bm{x}_{k})\|+\|\nabla f(\bm{x}_{k})\|\right),$

for $k$ sufficiently large. As $f,\bm{h}\in\mathcal{C}^{1}$ and
$\mu_{j}(\bm{x}_{k},1)=0$ for $j\not\in\mathcal{A}$, it follows from (55) that
the sequences $\\{\bm{\lambda}_{k}\\}_{k=0}^{\infty}$ and
$\\{\bm{\mu}(\bm{x}_{k},1)\\}_{k=0}^{\infty}$ are uniformly bounded for $k$
sufficiently large. As the sequence
$\\{\bm{\bar{\lambda}}_{k}\\}_{k=0}^{\infty}$ is uniformly bounded by
construction, there exists a constant $\Lambda$ such that
$\|\bm{\lambda}_{k}-\bm{\bar{\lambda}}_{k-1}\|\leq\Lambda$ for all $k$
sufficiently large. Now using the update formula
$\bm{\lambda}_{k}=\bm{\bar{\lambda}}_{k-1}+2q_{k}\bm{h}(\bm{x}_{k})$, we have
that

(56)
$\displaystyle\|\bm{h}(\bm{x}_{k})\|=\frac{\|\bm{\lambda}_{k}-\bm{\bar{\lambda}}_{k-1}\|}{2q_{k}}\leq\frac{\Lambda}{2q_{k}},$

for all $k$ sufficiently large. So (56) together with the hypothesis that
$q_{k}\to\infty$ as $k\to\infty$, yields that

(57)
$\displaystyle\lim_{k\to\infty}E_{c}(\bm{x}_{k})=\lim_{k\to\infty}\|\bm{h}(\bm{x}_{k})\|^{2}\leq\lim_{k\to\infty}\frac{\Lambda^{2}}{4q_{k}^{2}}=0.$

As
$E_{0}(\bm{x}_{k},\bm{\lambda}_{k},\bm{\mu}(\bm{x}_{k},1))^{2}=E_{m,0}(\bm{x}_{k},\bm{\lambda}_{k},\bm{\mu}(\bm{x}_{k},1))+E_{c}(\bm{x}_{k})$
and $E_{m,0}(\bm{x}_{k},\bm{\lambda}_{k},\bm{\mu}(\bm{x}_{k},1))\leq\theta
E_{c}(\bm{x}_{k})$ for all $k$ by hypothesis, it now follows that
$E_{0}(\bm{x}^{*},\bm{\lambda}^{*},\bm{\mu}(\bm{x}^{*},1))=0$. Lastly, from
Lemma 2.1 it follows that
$E_{1}(\bm{x}^{*},\bm{\lambda}^{*},\bm{\mu}(\bm{x}^{*},1))=0$ which concludes
the proof.

###  4.2 Global Convergence of GS Algorithm using Modified PASA

We now focus on establishing a global convergence result for the augmented
Lagrangian problem when using Algorithm 4, PASA [21] with modified stopping
criterion. Algorithm 4 is obtained by replacing all instances of
$E_{\scaleto{PASA}{4pt}}(\bm{x})$ in the original PASA with
$E_{m,0}(\bm{x},\bm{\bar{\lambda}}_{k}+2q_{k}\bm{h}(\bm{x}),\bm{\mu}(\bm{x},1))=-\nabla\mathcal{L}_{q}(\bm{x},\bm{\bar{\lambda}}_{k})(\bm{y}(\bm{x},1)-\bm{x})$.
First, we provide a necessary lemma that considers the relationship of the
multiplier error estimator
$E_{m,0}(\bm{x},\bm{\bar{\lambda}}_{k}+2q_{k}\bm{h}(\bm{x}),\bm{\mu}(\bm{x},\alpha)/\alpha)$
to
$E_{m,0}(\bm{x},\bm{\bar{\lambda}}_{k}+2q_{k}\bm{h}(\bm{x}),\bm{\mu}(\bm{x},\beta)/\beta)$
for various values of $0<\alpha\leq\beta$.

######  Lemma 4.5.

For a given $\bm{x}$, let $\bm{\mu}(\bm{x},\alpha)$ satisfy the KKT conditions
in (42) – (45). Then

(58)
$\displaystyle-\nabla\mathcal{L}_{q_{k}}(\bm{x},\bm{\bar{\lambda}}_{k})(\bm{y}(\bm{x},\alpha)-\bm{x})\leq-\nabla_{x}\mathcal{L}_{q_{k}}(\bm{x},\bm{\bar{\lambda}}_{k})(\bm{y}(\bm{x},\beta)-\bm{x}),\
\text{for}\ 0<\alpha\leq\beta,$ (59) $\displaystyle
E_{m,0}(\bm{x},\bm{\bar{\lambda}}_{k}+2q_{k}\bm{h}(\bm{x}),\bm{\mu}(\bm{x},\beta)/\beta)\leq
E_{m,0}(\bm{x},\bm{\bar{\lambda}}_{k}+2q_{k}\bm{h}(\bm{x}),\bm{\mu}(\bm{x},\alpha)/\alpha),\
\text{for}\ 0<\alpha\leq\beta\leq 1.$

######  Proof 4.6.

Suppose $0<\alpha\leq\beta$. Let
$\bm{d}^{\alpha}(\bm{x})=\bm{y}(\bm{x},\alpha)-\bm{x}$. By P1 in [19] with
$\bm{x}\leftarrow\bm{x}-\alpha\bm{g}(\bm{x})$ and
$\bm{y}\leftarrow\bm{y}(\bm{x},\beta)$, we have that

(60) $\displaystyle 0$
$\displaystyle\leq(\bm{y}(\bm{x},\alpha)-(\bm{x}-\alpha\bm{g}(\bm{x})))^{\intercal}(\bm{y}(\bm{x},\beta)-\bm{y}(\bm{x},\alpha))$
(61)
$\displaystyle=(\bm{d}^{\alpha}(\bm{x})+\alpha\bm{g}(\bm{x}))^{\intercal}(\bm{d}^{\beta}(\bm{x})-\bm{d}^{\alpha}(\bm{x}))$
(62)
$\displaystyle=\bm{d}^{\alpha}(\bm{x})^{\intercal}\bm{d}^{\beta}(\bm{x})-\|\bm{d}^{\alpha}(\bm{x})\|^{2}+\alpha\bm{g}(\bm{x})^{\intercal}\bm{d}^{\beta}(\bm{x})-\alpha\bm{g}(\bm{x})^{\intercal}\bm{d}^{\alpha}(\bm{x}).$

Rearranging terms and dividing both sides by $\alpha$ yields

(63) $\displaystyle-\bm{g}(\bm{x})^{\intercal}\bm{d}^{\beta}(\bm{x})$
$\displaystyle\leq-\bm{g}(\bm{x})^{\intercal}\bm{d}^{\alpha}(\bm{x})+\frac{1}{\alpha}\bm{d}^{\alpha}(\bm{x})^{\intercal}\bm{d}^{\beta}(\bm{x})-\frac{1}{\alpha}\|\bm{d}^{\alpha}(\bm{x})\|^{2}.$

It now follows from Cauchy-Schwarz that

(64) $\displaystyle-\bm{g}(\bm{x})^{\intercal}\bm{d}^{\beta}(\bm{x})$
$\displaystyle\leq-\bm{g}(\bm{x})^{\intercal}\bm{d}^{\alpha}(\bm{x})+\frac{1}{\alpha}\|\bm{d}^{\alpha}(\bm{x})\|\|\bm{d}^{\beta}(\bm{x})\|-\frac{1}{\alpha}\|\bm{d}^{\alpha}(\bm{x})\|^{2}.$

By swapping $\alpha$ and $\beta$ in the construction of (64), we obtain the
inequality

(65) $\displaystyle-\bm{g}(\bm{x})^{\intercal}\bm{d}^{\alpha}(\bm{x})$
$\displaystyle\leq-\bm{g}(\bm{x})^{\intercal}\bm{d}^{\beta}(\bm{x})+\frac{1}{\beta}\|\bm{d}^{\beta}(\bm{x})\|\|\bm{d}^{\alpha}(\bm{x})\|-\frac{1}{\beta}\|\bm{d}^{\beta}(\bm{x})\|^{2}.$

As $\alpha\leq\beta$, by P4 in [19] we have that
$\|\bm{d}^{\alpha}(\bm{x})\|\leq\|\bm{d}^{\beta}(\bm{x})\|$. Substituting this
into (65) yields

(66) $\displaystyle-\bm{g}(\bm{x})^{\intercal}\bm{d}^{\alpha}(\bm{x})$
$\displaystyle\leq-\bm{g}(\bm{x})^{\intercal}\bm{d}^{\beta}(\bm{x})+\frac{1}{\beta}\|\bm{d}^{\beta}(\bm{x})\|^{2}-\frac{1}{\beta}\|\bm{d}^{\beta}(\bm{x})\|^{2}=-\bm{g}(\bm{x})^{\intercal}\bm{d}^{\beta}(\bm{x}).$

which is equivalent to (58).

Now suppose that $\alpha\leq\beta\leq 1$. By P5 in [19], we have that
$\|\bm{d}^{\beta}(\bm{x})\|\leq\frac{\beta}{\alpha}\|\bm{d}^{\alpha}(\bm{x})\|$.
Substituting this into (64) yields

(67) $\displaystyle-\bm{g}(\bm{x})^{\intercal}\bm{d}^{\beta}(\bm{x})$
$\displaystyle\leq-\bm{g}(\bm{x})^{\intercal}\bm{d}^{\alpha}(\bm{x})+\frac{\beta}{\alpha^{2}}\|\bm{d}^{\alpha}(\bm{x})\|^{2}-\frac{1}{\alpha}\|\bm{d}^{\alpha}(\bm{x})\|^{2}$
(68)
$\displaystyle=-\bm{g}(\bm{x})^{\intercal}\bm{d}^{\alpha}(\bm{x})+\left(\frac{1}{\alpha^{2}}-\frac{1}{\alpha}\right)\|\bm{d}^{\alpha}(\bm{x})\|^{2}-\frac{1-\beta}{\alpha^{2}}\|\bm{d}^{\alpha}(\bm{x})\|^{2}.$

Again, by P5 in [19] we have that
$-\frac{1}{\alpha^{2}}\|\bm{d}^{\alpha}(\bm{x})\|^{2}\leq-\frac{1}{\beta^{2}}\|\bm{d}^{\beta}(\bm{x})\|$
and since $1-\beta\geq 0$ it follows that

(69)
$\displaystyle-\frac{1-\beta}{\alpha^{2}}\|\bm{d}^{\alpha}(\bm{x})\|^{2}\leq-\frac{1-\beta}{\beta^{2}}\|\bm{d}^{\beta}(\bm{x})\|^{2}=-\left(\frac{1}{\beta^{2}}-\frac{1}{\beta}\right)\|\bm{d}^{\beta}(\bm{x})\|^{2}.$

Combining (68) and (69) and rearranging terms yields

(70)
$\displaystyle-\bm{g}(\bm{x})^{\intercal}\bm{d}^{\beta}(\bm{x})+\left(\frac{1}{\beta^{2}}-\frac{1}{\beta}\right)\|\bm{d}^{\beta}(\bm{x})\|^{2}$
$\displaystyle\leq-\bm{g}(\bm{x})^{\intercal}\bm{d}^{\alpha}(\bm{x})+\left(\frac{1}{\alpha^{2}}-\frac{1}{\alpha}\right)\|\bm{d}^{\alpha}(\bm{x})\|^{2}$

which by (46) is equivalent to (59).

Before providing a convergence result for this approach to phase one of NPASA,
we consider how convergence of Algorithm 4 for the minimization problem in the
GS algorithm relates to the multiplier error estimator of the original
nonlinear program, (3). In particular, since
$E_{m,1}(\bm{x}_{k},\bm{\lambda}_{k},\bm{\mu}_{k})\leq
E_{m,0}(\bm{x}_{k},\bm{\lambda}_{k},\bm{\mu}_{k})$ and our branching criterion
in phase one of NPASA is given by
$E_{m,1}(\bm{x}_{k},\bm{\lambda}_{k},\bm{\mu}_{k})\leq\theta
E_{c}(\bm{x}_{k})$, we would like some conditions under which Algorithm 4
converges to a point $\bm{x}_{k}$ at which
$E_{m,0}(\bm{x}_{k},\bm{\lambda}_{k},\bm{\mu}_{k})\approx 0$. Before providing
such a result, we need to set up the assumptions that are used in this result
as they are based on details of PASA [20].

Given $\bm{u}\in\mathbb{R}^{n}$, $\bm{\nu}\in\mathbb{R}^{\ell}$, and
$q\in\mathbb{R}$ we define the level set

(71)
$\displaystyle\mathcal{S}(\bm{u},\bm{\nu},q):=\\{\bm{x}\in\Omega:\mathcal{L}_{q}(\bm{x},\bm{\nu})\leq\mathcal{L}_{q}(\bm{u},\bm{\nu})\\}.$

Additionally, we let the set $\mathcal{D}(\bm{u},\bm{\nu},q)$ denote the
search directions generated by phase one of the algorithm PASA [20] when
solving the problem

(74)
$\displaystyle\begin{array}[]{cc}\displaystyle\min_{\bm{x}}&\displaystyle\mathcal{L}_{q}(\bm{x},\bm{\nu})\\\
\text{s.t.}&\bm{x}\in\Omega\end{array}$

starting at the initial guess $\bm{x}=\bm{u}$. Note that the definition of the
search directions generated in phase one of PASA can be found in Algorithm 1
in Section 2 of [20]. We are now ready to state the assumptions:

  1. (G1)

Given $\bm{u}\in\mathbb{R}^{n}$ and $\bm{\nu}\in\mathbb{R}^{\ell}$,
$\mathcal{L}_{q}(\bm{x},\bm{\nu})$ is bounded from below on the level set
$\mathcal{S}(\bm{u},\bm{\nu})$ and
$d_{\max}=\sup_{\bm{d}\in\mathcal{D}(\bm{u},\bm{\nu},q)}\|\bm{d}\|<\infty$.

  2. (G2)

Given $\bm{u}\in\mathbb{R}^{n}$ and $\bm{\nu}\in\mathbb{R}^{\ell}$, if
$\overline{\mathcal{S}}(\bm{u},\bm{\nu})$ is the collection of
$\bm{x}\in\Omega$ whose distance to $\mathcal{S}(\bm{u},\bm{\nu})$ is at most
$d_{\max}$, then $\nabla_{x}\mathcal{L}_{q}(\bm{x},\bm{\nu})$ is Lipschitz
continuous on $\overline{\mathcal{S}}(\bm{u},\bm{\nu})$.

With Lemma 4.5 established and (G1) and (G2) stated, we are now ready to
provide a convergence result for the nonmonotone gradient projection algorithm
(NGPA) in phase one of PASA. that is necessary to establish convergence when
using Algorithm 4 to solve problem (36).

######  Theorem 4.7.

Given $\bm{u}_{0}\in\mathbb{R}^{n}$, $\bm{\bar{\lambda}}_{k}$, and $q_{k}$
suppose that (G1) and (G2) hold at
$(\bm{u}_{0},\bm{\bar{\lambda}}_{k},q_{k})$. If the nonmonotone gradient
projection algorithm (NGPA) in [19] with $\varepsilon=0$ is used to solve
problem (36) then NGPA either terminates in a finite number of iterations at a
stationary point, or generates a sequence of iterates
$\\{\bm{u}_{i}\\}_{i=0}^{\infty}$ satisfying

(75)
$\displaystyle\liminf_{i\to\infty}E_{m,0}(\bm{u}_{i},\bm{\bar{\lambda}}_{k}+2q_{k}\bm{h}(\bm{u}_{i}),\bm{\mu}(\bm{u}_{i},1))=0.$

######  Proof 4.8.

Let $F(\bm{x})=\mathcal{L}_{q_{k}}(\bm{x},\bm{\bar{\lambda}}_{k})$. First,
note that from (3.6) in [21], the stepsize $s_{k}$ in the Armijo line search
is bounded below by
$s_{k}\geq\min\left\\{1,\frac{2\eta(1-\delta)}{\kappa\alpha_{\max}}\right\\}=:c$,
where $\kappa$ is the Lipschitz constant for
$\nabla_{x}\mathcal{L}_{q_{k}}(\bm{x},\bm{\bar{\lambda}}_{k})$. From equation
(2.16) in [19] we have that

(76) $\displaystyle-\delta
c\bm{g}(\bm{u}_{i})^{\intercal}(\bm{y}(\bm{u}_{i},\bar{\alpha}_{i})-\bm{u}_{i})\leq
F_{i}^{r}-F(\bm{u}_{i+1}),$

for some $\bar{\alpha}_{i}\in[\alpha_{\min},\alpha_{\max}]$. We now proceed by
way of contradiction. Accordingly, suppose that

(77)
$\displaystyle\liminf_{i\to\infty}E_{m,0}(\bm{u}_{i},\bm{\bar{\lambda}}_{k}+2q_{k}\bm{h}(\bm{u}_{i}),\bm{\mu}(\bm{u}_{i},1))>0.$

Then there exists a constant $\tau>0$ such that

(78) $\displaystyle\tau\leq
E_{m,0}(\bm{u}_{i},\bm{\bar{\lambda}}_{k}+2q_{k}\bm{h}(\bm{u}_{i}),\bm{\mu}(\bm{u}_{i},1))=-\bm{g}(\bm{u}_{i})^{\intercal}(\bm{y}(\bm{u}_{i},1)-\bm{u}_{i}),$

for all $i\geq 0$. Suppose $\bar{\alpha}_{i}\geq 1$. By (58) in Theorem 4.5
with $\alpha=1$ and $\beta=\bar{\alpha}_{i}$, we have

(79)
$\displaystyle\tau\leq-\bm{g}(\bm{u}_{i})^{\intercal}(\bm{y}(\bm{u}_{i},1)-\bm{u}_{i})\leq-\bm{g}(\bm{u}_{i})^{\intercal}(\bm{y}(\bm{u}_{i},\bar{\alpha}_{i})-\bm{u}_{i}).$

On the other hand, if $\bar{\alpha}_{i}\in(0,1)$ then by (59) in Theorem 4.5
with $\alpha=\bar{\alpha}_{i}$ and $\beta=1$ we have

(80)
$\displaystyle-\bm{g}(\bm{u}_{i})^{\intercal}(\bm{y}(\bm{u}_{i},1)-\bm{u}_{i})\leq-\bm{g}(\bm{u}_{i})^{\intercal}(\bm{y}(\bm{u}_{i},\bar{\alpha}_{i})-\bm{u}_{i})+\left(\frac{1}{\bar{\alpha}_{i}}-1\right)\frac{\|\bm{y}(\bm{u}_{i},\bar{\alpha}_{i})-\bm{u}_{i}\|^{2}}{\bar{\alpha}_{i}}.$

From P6 in [19] we have that

(81)
$\displaystyle\frac{\|\bm{y}(\bm{u}_{i},\bar{\alpha}_{i})-\bm{u}_{i}\|^{2}}{\bar{\alpha}_{i}}\leq-\bm{g}(\bm{u}_{i})^{\intercal}(\bm{y}(\bm{u}_{i},\bar{\alpha}_{i}).$

As $\bar{\alpha}_{i}\in(0,1)$, we have that
$\left(\frac{1}{\bar{\alpha}_{i}}-1\right)>0$ which combined with (81) yields

(82)
$\displaystyle\left(\frac{1}{\bar{\alpha}_{i}}-1\right)\frac{\|\bm{y}(\bm{u}_{i},\bar{\alpha}_{i})-\bm{u}_{i}\|^{2}}{\bar{\alpha}_{i}}$
$\displaystyle\leq\left(\frac{1}{\bar{\alpha}_{i}}-1\right)\left(-\bm{g}(\bm{u}_{i})^{\intercal}(\bm{y}(\bm{u}_{i},\bar{\alpha}_{i})\right)=\left(1-\frac{1}{\bar{\alpha}_{i}}\right)\bm{g}(\bm{u}_{i})^{\intercal}(\bm{y}(\bm{u}_{i},\bar{\alpha}_{i}).$

From (80) and (82) it now follows that

(83)
$\displaystyle-\bm{g}(\bm{u}_{i})^{\intercal}(\bm{y}(\bm{u}_{i},1)-\bm{u}_{i})\leq-\frac{1}{\bar{\alpha}_{i}}\bm{g}(\bm{u}_{i})^{\intercal}(\bm{y}(\bm{u}_{i},\bar{\alpha}_{i})-\bm{u}_{i}).$

As $\bar{\alpha}_{i}\geq\alpha_{\min}$, (78) and (83) now yield

(84)
$\displaystyle\alpha_{\min}\tau\leq-\bm{g}(\bm{u}_{i})^{\intercal}(\bm{y}(\bm{u}_{i},\bar{\alpha}_{i})-\bm{u}_{i}).$

Hence, by observing (79) and (84) we have established that for any choice of
$\bar{\alpha}_{i}\in[\alpha_{\min},\alpha_{\max}]$ there exists a constant
$\nu$ such that

(85)
$\displaystyle\nu\leq-\bm{g}(\bm{u}_{i})^{\intercal}(\bm{y}(\bm{u}_{i},\bar{\alpha}_{i})-\bm{u}_{i}).$

Combining (76) and (85) yields $F(\bm{u}_{i+1})\leq F_{i}^{r}-\delta c\nu$
which, by the same argument in the proof of Theorem 2.2 [19], results in a
contradiction. Thus,
$\liminf_{i\to\infty}E_{m,0}(\bm{u}_{i},\bm{\bar{\lambda}}_{k}+2q_{k}\bm{h}(\bm{u}_{i}),\bm{\mu}(\bm{u}_{i},1))=0$.

We now establish a global convergence result for Algorithm 4 when solving this
problem. We note that this is similar to Theorem 3.2 in [21].

######  Theorem 4.9.

Given $\bm{u}_{0}\in\mathbb{R}^{n}$, $\bm{\bar{\lambda}}_{k}$, and $q_{k}$
suppose that (G1) and (G2) hold at
$(\bm{u}_{0},\bm{\bar{\lambda}}_{k},q_{k})$. If Algorithm 4 with
$\varepsilon=0$ is used to solve problem (36) then Algorithm 4 either
terminates in a finite number of iterations at a stationary point, or
generates a sequence of iterates $\\{\bm{u}_{i}\\}_{i=0}^{\infty}$ satisfying

(86)
$\displaystyle\liminf_{i\to\infty}E_{m,0}(\bm{u}_{i},\bm{\bar{\lambda}}_{k}+2q_{k}\bm{h}(\bm{u}_{i}),\bm{\mu}(\bm{u}_{i},1))=0.$

######  Proof 4.10.

The proof is split into three cases.
Case 1: Phase two of Algorithm 4 is executed a finite number of times. In this
case, only phase one of Algorithm 4 is executed for $i$ sufficiently large. As
such, from Theorem 4.7 it follows that (86) holds.
Case 2: Phase one of Algorithm 4 is executed a finite number of times. Here,
only phase two is executed for $i$ sufficiently large. Since phase two
eventually does not branch back to phase one, by the branching criterion in
phase two we have that

(87) $\displaystyle\theta
E_{m,0}(\bm{u}_{i},\bm{\bar{\lambda}}_{k}+2q_{k}\bm{h}(\bm{u}_{i}),\bm{\mu}(\bm{u}_{i},1))\leq
e_{\scaleto{PASA}{4pt}}(\bm{u}_{i})$

for $i$ sufficiently large. Using an argument found in the proof of Theorem
3.2 from [21] we have that (86) holds.
Case 3: Phases one and two of Algorithm 4 are executed an infinite number of
times. By way of contradiction, assume that (86) does not hold. Then there
exists a constant $\tau>0$ such that

(88) $\displaystyle\tau\leq
E_{m,0}(\bm{u}_{i},\bm{\bar{\lambda}}_{k}+2q_{k}\bm{h}(\bm{u}_{i}),\bm{\mu}(\bm{u}_{i},1))=-\bm{g}(\bm{u}_{i})^{\intercal}(\bm{y}(\bm{u}_{i},1)-\bm{u}_{i}).$

Letting
$F(\bm{u}_{i+1})=\mathcal{L}_{q_{k}}(\bm{u}_{i+1},\bm{\bar{\lambda}}_{k})$ and
applying the same argument as in the proof of Theorem 4.7, we have that there
exists a constant $\nu$ such that

(89) $\displaystyle F(\bm{u}_{i+1})\leq F(\bm{u}_{i})-\delta c\nu$

for each iteration of GPA when $i$ is sufficiently large. Additionally, by
property F1 of the LCO in [21], we have that $F(\bm{u}_{i+1})\leq
F(\bm{u}_{i})$. Since there are an infinite number of iterations in phase one,
it follows that the cost function is unbounded from below which contradictions
assumption (G1). Thus, (86) holds.

Based on the analysis in this section, Algorithm 4 should be used to solve the
minimization problem in Algorithm 1, the GS algorithm. To conclude this
section, we note that the inequality multiplier $\bm{\mu}(\bm{x}_{k},1)$
required to check the convergence criterion in phase one of NPASA can be
explicitly constructed using values computed during the solution of the GS
algorithm. Details of this construction are provided in the companion paper
[11].

##  5 Convergence Analysis for Phase Two of NPASA

In this section, we provide conditions under which the constraint and
multiplier steps in the LS algorithm satisfy desirable convergence properties.
We first provide a result on the perturbed Newton step scheme used in the
constraint step in Section 5.1. Then we establish a convergence result for the
constraint step in Section 5.2 followed by a result for the multiplier step of
NPASA in Section 5.3.

###  5.1 Algorithm with Infeasibility Detection for Solving Constraint Step

One method for minimizing the violation of the equality constraints in problem
(3) is to generate a sequence of iterates $\\{\bm{w}_{i}\\}_{i=1}^{\infty}$
defined by $\bm{w}_{i+1}=\bm{w}_{i}+\bm{d}_{i}$, where $\bm{d}_{i}$ is the
solutioin to the minimization problem given by
$\min\left\\{\|\bm{w}-\bm{d}_{i-1}\|^{2}:\nabla\bm{h}(\bm{d}_{i-1})(\bm{w}-\bm{d}_{i-1})=-\bm{h}(\bm{d}_{i-1}),\bm{w}\in\Omega\right\\}$.
Given $\bar{\bm{x}}$, PPROJ [21] uses a dual method for solving the projection
problem

(92)
$\displaystyle\begin{array}[]{cc}\displaystyle\min_{\bm{x}}&\frac{1}{2}\|\bar{\bm{x}}-\bm{x}\|^{2}\\\
\text{s.t.}&\bm{l}\leq\bm{A}\bm{x}\leq\bm{u},\ \ \bm{x}\geq\bm{0}.\end{array}$

then reconstructs the primal solution. We would like to use this algorithm to
solve for the Newton direction, $\bm{d}_{i}$. However, since
$\\{\bm{w}\in\mathbb{R}^{n}:\bm{h}(\bm{w})=\bm{0},\bm{w}\in\Omega\\}\neq\emptyset$
does not guarantee

(93)
$\displaystyle\left\\{\bm{w}\in\mathbb{R}^{n}:\nabla\bm{h}(\bm{d}_{i})(\bm{w}-\bm{d}_{i})=-\bm{h}(\bm{d}_{i}),\bm{w}\in\Omega\right\\}\neq\emptyset,$

it may be necessary to perturb the constraint
$\nabla\bm{h}(\bm{d}_{i})(\bm{w}-\bm{d}_{i})=-\bm{h}(\bm{d}_{i})$ in order to
solve for $\bm{d}_{i+1}$.

In order to simplify the notation in our discussion of this problem, we
express the linearized constraint set as
$\left\\{\bm{w}\in\mathbb{R}^{n}:\bm{M}(\bm{w}-\bar{\bm{w}})=\bm{c},\bm{w}\geq\bm{0}\right\\}$
supposing that we are given $\bm{M}\in\mathbb{R}^{m\times n}$,
$\bm{c}\in\mathbb{R}^{m}$, and $\bar{\bm{w}}\in\mathbb{R}^{n}$. Ideally, we
would like to determine the solution to the problem

(94) $\displaystyle\min_{\bm{w}\in\mathcal{S}}\
\frac{1}{2}\|\bm{w}-\bar{\bm{w}}\|^{2},$

where
$\mathcal{S}=\arg\min_{\bm{w}\geq\bm{0}}\left\\{\|\bm{c}-\bm{M}(\bm{w}-\bar{\bm{w}})\|^{2}\right\\}$.
Since problem (94) involves minimizing a strongly convex objective function
over a convex set, there exists a unique solution which we will denote by
$\bm{w}_{\infty}$. As such, every search direction from $\bm{w}_{\infty}$
leading to a point in $\mathcal{S}$ must be an ascent direction. Hence, the
first order optimality conditions for $\bm{w}_{\infty}$ can be expressed as

(95)
$\displaystyle\langle\bm{w}_{\infty}-\bar{\bm{w}},\bm{w}-\bm{w}_{\infty}\rangle\geq
0,\ \text{for all}\ \bm{w}\in\mathcal{S}.$

Since determining elements of the constraint set for problem (94) requires
knowing all minimizers of $\|\bm{c}-\bm{M}(\bm{w}-\bar{\bm{w}})\|$ we consider
an alternative problem.

To simplify this problem, we define the sets
$\mathcal{F}=\left\\{[\bm{w},\bm{y}]^{\intercal}\in\mathbb{R}^{n+m}:\bm{M}(\bm{w}-\bar{\bm{w}})+\bm{y}=\bm{c},\bm{w}\geq\bm{0}\right\\}$.
Now, given a penalty parameter $p>0$, the problem is defined by

(96) $\displaystyle\min_{[\bm{w},\bm{y}]^{\intercal}\in\mathcal{F}}\ \
f_{p}(\bm{w},\bm{y}):=\frac{1}{2}\|\bm{w}-\bar{\bm{w}}\|^{2}+\frac{p}{2}\|\bm{y}\|^{2}.$

Since problem (96) involves minimizing a strongly convex objective function
over a convex set, it has a unique solution which we will denote by
$[\bm{w}_{p},\bm{y}_{p}]^{\intercal}$. Note that the first order optimality
conditions for $[\bm{w}_{p},\bm{y}_{p}]^{\intercal}$ can be expressed as

(97)
$\displaystyle\langle\bm{w}_{p}-\bar{\bm{w}},\bm{w}-\bm{w}_{p}\rangle+p\langle\bm{y}_{p},\bm{y}-\bm{y}_{p}\rangle\geq
0,\ \text{for all}\ [\bm{w},\bm{y}]^{\intercal}\in\mathcal{F}.$

While the number of primal variables in (96) is greater than that of problem
(94), the dimension of the dual variables is the same for problems (94) and
(96). Hence, by using the dual method PPROJ [21] we can solve problem (96)
without increasing the dimension of the optimization problem. The following
result provides valuable insight on the relationship between $\bm{w}_{p}$ and
$\bm{w}_{\infty}$. Note that we make use of the set
$\mathcal{F}_{\delta}=\left\\{[\bm{w},\bm{y}]^{\intercal}\in\mathbb{R}^{n+m}:\bm{M}(\bm{w}-\bar{\bm{w}})+\bm{y}=\bm{c},\bm{w}\geq\bm{0},\|\bm{y}\|^{2}\leq\delta\right\\}$.
in the proof of the following result.

######  Theorem 5.1.

Let $[\bm{w}_{p},\bm{y}_{p}]^{\intercal}$ be the unique solution to problem _(
96)_. Then $\|\bm{w}_{p}-\bm{x}_{\infty}\|=O\left(\frac{1}{p}\right)$, for all
$p>0$ sufficiently large.

######  Proof 5.2.

Letting
$\gamma=\min_{\bm{w}\geq\bm{0}}\|\bm{c}-\bm{M}(\bm{w}-\bar{\bm{w}})\|^{2}$ it
follows that
$\mathcal{S}=\\{\bm{w}\in\mathbb{R}^{n}:[\bm{w},\bm{y}]^{\intercal}\in\mathcal{F}_{\gamma}\\}$.
So $\mathcal{S}\subseteq\mathcal{F}$ and by defining
$\bm{y}_{\infty}=\bm{c}-\bm{M}(\bm{w}_{\infty}-\bar{\bm{w}})$ we have that
$[\bm{w}_{\infty},\bm{y}_{\infty}]^{\intercal}\in\mathcal{F}$. Since
$[\bm{w}_{p},\bm{y}_{p}]^{\intercal}$ minimizes the objective in problem (96),
we now have that

(98)
$\displaystyle\frac{1}{2}\|\bm{w}_{p}-\bar{\bm{w}}\|^{2}+\frac{p}{2}\|\bm{y}_{p}\|^{2}$
$\displaystyle\leq\frac{1}{2}\|\bm{w}_{\infty}-\bar{\bm{w}}\|^{2}+\frac{p}{2}\|\bm{y}_{\infty}\|^{2}=\frac{1}{2}\|\bm{w}_{\infty}-\bar{\bm{w}}\|^{2}+\frac{p}{2}\gamma.$

Using (98) we observe that
$\|\bm{w}_{p}-\bar{\bm{w}}\|^{2}\leq\|\bm{w}_{\infty}-\bar{\bm{w}}\|^{2}+p\left(\gamma-\|\bm{y}_{p}\|^{2}\right)\leq\|\bm{w}_{\infty}-\bar{\bm{w}}\|^{2}$,
where the final inequality holds since $\|\bm{y}_{p}\|^{2}\geq\gamma$, for all
$p\geq 1$. Hence,

(99)
$\displaystyle\|\bm{w}_{p}\|\leq\|\bm{w}_{\infty}-\bar{\bm{w}}\|+\|\bar{\bm{w}}\|.$

Also, from (98) it follows that
$\|\bm{y}_{p}\|^{2}\leq\frac{1}{p}\|\bm{w}_{\infty}-\bar{\bm{w}}\|^{2}+\gamma$
which yields that $\lim_{p\to\infty}\|\bm{y}_{p}\|^{2}=\gamma$. Since
$\|\bm{w}_{p}\|$ is bounded for all values of $p$ in (99), there exist
convergent subsequences of $\\{\bm{w}_{p}\\}_{p=1}^{\infty}$. Letting
$\\{\bm{w}_{p_{k}}\\}_{k=1}^{\infty}$ be a convergent subsequence of
$\\{\bm{w}_{p}\\}_{p=1}^{\infty}$ with limit $\tilde{\bm{w}}$ we observe that

(100) $\displaystyle\|\bm{c}-\bm{M}(\tilde{\bm{w}}-\bar{\bm{w}})\|^{2}$
$\displaystyle=\lim_{k\to\infty}\|\bm{c}-\bm{M}(\bm{w}_{p_{k}}-\bar{\bm{w}})\|^{2}=\lim_{k\to\infty}\|\bm{y}_{p_{k}}\|^{2}=\gamma.$

So, taking $\tilde{\bm{y}}=\bm{c}-\bm{M}(\tilde{\bm{w}}-\bar{\bm{w}})$, we
have
$\lim_{k\to\infty}[\bm{w}_{p_{k}},\bm{y}_{p_{k}}]^{\intercal}=[\tilde{\bm{w}},\tilde{\bm{y}}]^{\intercal}\in\mathcal{F}_{\gamma}$.
Hence, $\tilde{\bm{w}}\in\mathcal{S}$.

We now claim that $\tilde{\bm{w}}=\bm{w}_{\infty}$. Observing the optimality
conditions in (95) it follows that we need to show that
$\langle\tilde{\bm{w}}-\bar{\bm{w}},\bm{w}-\tilde{\bm{w}}\rangle\geq 0$, for
all $\bm{w}\in\mathcal{S}$. Let $p\geq 1$ be given. By the optimality
conditions in (97), we have that

(101) $\displaystyle\langle\bm{w}_{p}-\bar{\bm{w}},\bm{w}-\bm{w}_{p}\rangle$
$\displaystyle\geq p\langle\bm{y}_{p},\bm{y}_{p}-\bm{y}\rangle,$

for all $[\bm{w},\bm{y}]^{\intercal}\in\mathcal{F}$. As
$\mathcal{S}=\\{\bm{w}\in\mathbb{R}^{n}:[\bm{w},\bm{y}]^{\intercal}\in\mathcal{F}_{\gamma}\\}$
and $\mathcal{F}_{\gamma}\subseteq\mathcal{F}$, it suffices to show that
$\langle\bm{y}_{p},\bm{y}_{p}-\bm{y}\rangle\geq 0$, for all
$\bm{y}\in\mathbb{R}^{m}$ with $\|\bm{y}\|^{2}=\gamma$. Accordingly, let
$\bm{y}\in\mathbb{R}^{m}$ with $\|\bm{y}\|^{2}=\gamma$ be given and observe
that

(102) $\displaystyle\langle\bm{y}_{p},\bm{y}_{p}-\bm{y}\rangle$
$\displaystyle=\|\bm{y}_{p}\|^{2}-\frac{1}{2}\left(\|\bm{y}_{p}\|^{2}-\|\bm{y}_{p}-\bm{y}\|^{2}+\|\bm{y}\|^{2}\right)=\frac{1}{2}\left(\|\bm{y}_{p}\|^{2}+\|\bm{y}_{p}-\bm{y}\|^{2}-\gamma\right).$

By our choice of $\gamma$ we have that $\|\bm{y}_{p}\|^{2}\geq\gamma$. This
fact together with the inequality in (102) implies that
$\langle\bm{y}_{p},\bm{y}_{p}-\bm{y}\rangle\geq\frac{1}{2}\|\bm{y}_{p}-\bm{y}\|^{2}\geq
0$. Substituting this into (101) yields that
$\langle\bm{w}_{p}-\bar{\bm{w}},\bm{w}-\bm{w}_{p}\rangle\geq 0$, for all
$\bm{w}\in\mathcal{S}$. Thus,
$\langle\tilde{\bm{w}}-\bar{\bm{w}},\bm{w}-\tilde{\bm{w}}\rangle=\lim_{k\to\infty}\langle\bm{w}_{p_{k}}-\bar{\bm{w}},\bm{w}-\bm{w}_{p_{k}}\rangle\geq
0$, for all $\bm{w}\in\mathcal{S}$. So we conclude that
$\tilde{\bm{w}}=\bm{w}_{\infty}$.

Next we claim that the entire sequence $\\{\bm{w}_{p}\\}_{p=1}^{\infty}$
converges to $\bm{w}_{\infty}$. By way of contradiction, suppose that there
exists a subsequence $\\{\bm{w}_{p_{k}}\\}_{k=1}^{\infty}$ and a $\delta>0$
such that $\|\bm{w}_{p_{k}}-\bm{w}_{\infty}\|\geq\delta$ for all $k\geq 1$.
Since $\\{\bm{w}_{p_{k}}\\}_{k=1}^{\infty}$ is bounded it has a convergent
subsequence, say $\\{\bm{w}_{q_{k}}\\}_{k=1}^{\infty}$. As
$\\{\bm{w}_{q_{k}}\\}_{k=1}^{\infty}$ is a subsequence of
$\\{\bm{w}_{p}\\}_{p=1}^{\infty}$ and every convergent subsequence of
$\\{\bm{w}_{p}\\}_{p=1}^{\infty}$ converges to $\bm{w}_{\infty}$, there exists
a $M\in\mathbb{N}$ such that $\|\bm{w}_{q_{k}}-\bm{w}_{\infty}\|<\delta$ for
all $k\geq M$, a contradiction. Thus, every subsequence of
$\\{\bm{w}_{p}\\}_{p=1}^{\infty}$ converges to $\bm{w}_{\infty}$ implying that
$\bm{w}_{p}\to\bm{w}_{\infty}$ as $p\to\infty$.

Lastly, it remains to show that $\|\bm{w}_{p}-\bm{w}_{\infty}\|=O(1/p)$.
Accordingly, let $p\geq 1$ be given, let
$\mathcal{B}:=\\{i:\left(\bm{w}_{p}\right)_{i}=0\text{ and }1\leq i\leq n\\}$
denote the active set at point $\bm{w}_{p}$, and let
$\mathcal{I}:=\\{1,2,\ldots,n\\}\setminus\mathcal{B}$. We now consider the
problem

(103) $\displaystyle\min_{\bm{w}_{\mathcal{B}}=0}\
f(\bm{w}):=\frac{1}{2}\|\bm{w}-\bar{\bm{w}}\|^{2}+\frac{p}{2}\|\bm{c}+\bm{M}\bar{\bm{w}}-\bm{M}\bm{w}\|^{2}$

where $\bm{w}_{\mathcal{B}}$ is the vector with components $w_{i}$, for all
$i\in\mathcal{B}$. Letting $\bm{u}^{*}$ denote the solution to problem (103),
we claim that $\bm{u}^{*}=\bm{w}_{p}$. In order to justify this claim, first
note that since $\bm{w}_{p}$ satisfies the constraint for problem (103) we
have that

(104) $\displaystyle f(\bm{u}^{*})$ $\displaystyle\leq
f(\bm{w}_{p})=f_{p}(\bm{w}_{p},\bm{y}_{p}).$

Now suppose that $\bm{u}^{*}\neq\bm{w}_{p}$. By (104) and the strict convexity
of $f$ we have that $f(\bm{w})<f(\bm{w}_{p})$, for all
$\bm{w}\in[\bm{w}_{p},\bm{u}^{*}]$, that is for all $\bm{w}$ on the line
segment connecting $\bm{w}_{p}$ and $\bm{u}^{*}$. By our choice of
$\mathcal{B}$, $(\bm{w}_{p})_{i}>0$ for all $i\in\mathcal{I}$ and
$(\bm{w}_{p})_{i}=0=\bm{u}^{*}_{i}$ for all $i\in\mathcal{B}$. Hence,
$\bm{w}\geq 0$ for all $\bm{w}\in[\bm{w}_{p},\bm{u}^{*}]$ sufficiently close
to $\bm{w}_{p}$, a contradiction. Thus, $\bm{u}^{*}=\bm{w}_{p}$.

Now observe that we can reformulate problem (103) as the unconstrained
optimization problem

(105) $\displaystyle\min_{\bm{w}_{+}\in\mathbb{R}^{|\mathcal{I}|}}\ \ \ \
g(\bm{w}_{+}):=\frac{1}{2}\|\bm{w}_{+}-\bar{\bm{w}}_{+}\|^{2}+\frac{p}{2}\|\bm{c}+\bm{M}\bar{\bm{w}}-\bm{M}_{+}\bm{w}_{+}\|^{2}$

where $\bar{\bm{w}}_{+}\in\mathbb{R}^{|\mathcal{I}|}$ is the vector with
components $\bar{w}_{i}$, for all $i\in\mathcal{I}$, and $\bm{M}_{+}$ is the
matrix with columns $M_{i}$ of $\bm{M}$, for all $i\in\mathcal{I}$. Since
problem (105) is an unconstrained optimization problem with a strongly convex
objective function, the first order optimality conditions are given by

(106) $\displaystyle\bm{w}_{+}\ \text{minimizes}\ g(\bm{w})\ \ \text{if and
only if}\ \ \nabla
g(\bm{w}_{+})=\frac{1}{p}\left(\bm{w}_{+}-\bar{\bm{w}}_{+}\right)-\bm{M}_{+}^{\intercal}\left(\bm{c}+\bm{M}\bar{\bm{w}}-\bm{M}_{+}\bm{w}_{+}\right)=0.$

Letting $\bm{Q}\bm{\Sigma}\bm{P}^{\intercal}$ be a singular value
decomposition of $\bm{M}_{+}$ and performing the change of variables
$\bm{w}_{+}=\bm{P}\bm{z}$ we obtain $\nabla
g(\bm{w}_{+})=\bm{P}\left(\frac{1}{p}\bm{z}+\bm{\Sigma}^{\intercal}\bm{\Sigma}\bm{z}\right)-\bm{P}\bm{\Sigma}^{\intercal}\bm{Q}^{\intercal}\left(\bm{c}+\bm{M}\bar{\bm{w}}\right)-\frac{1}{p}\bar{\bm{w}}_{+}$.
Hence, multiplying both sides of equation (106) by $\bm{P}^{\intercal}$ and
simplifying yields

(107)
$\displaystyle\left(\frac{1}{p}\bm{I}+\bm{\Sigma}^{\intercal}\bm{\Sigma}\right)\bm{z}=\bm{\Sigma}^{\intercal}\bm{Q}^{\intercal}\left(\bm{c}+\bm{M}\bar{\bm{w}}\right)+\frac{1}{p}\bm{P}^{\intercal}\bar{\bm{w}}_{+}.$

Now let $r=rank(\bm{M}_{+})$ and let
$\sigma_{1}\geq\sigma_{2}\geq\ldots\geq\sigma_{r}$ denote the non-zero
singular values of $\bm{M}_{+}$. From equation (107), for $1\leq i\leq r$, we
have
$\left(\frac{1}{p}+\sigma_{i}^{2}\right)z_{i}=\sigma_{i}Q_{i}^{\intercal}(\bm{c}+\bm{M}\bar{\bm{w}})+\frac{1}{p}\left(\bm{P}^{\intercal}\bar{\bm{w}}_{+}\right)_{i}$,
or equivalently,

(108) $\displaystyle z_{i}$
$\displaystyle=\frac{\sigma_{i}}{\sigma_{i}^{2}+1/p}Q_{i}^{\intercal}(\bm{c}+\bm{M}\bar{\bm{w}})+\frac{1/p}{\sigma_{i}^{2}+1/p}\left(\bm{P}^{\intercal}\bar{\bm{w}}_{+}\right)_{i}=\frac{1}{\sigma_{i}}Q_{i}^{\intercal}(\bm{c}+\bm{M}\bar{\bm{w}})+O\left(\frac{1}{p}\right).$

For $i>r$ it also follows from equation (107) that
$z_{i}=\left(\bm{P}^{\intercal}\bar{\bm{w}}_{+}\right)_{i}$. Combining these
observations with (108) yields that there exists a vector
$\bm{u}^{\mathcal{B}}\in\mathbb{R}^{n}$ such that
$\bm{z}=\bm{u}^{\mathcal{B}}+O\left(1/p\right)$. As $\bm{w}_{+}=\bm{P}\bm{z}$,
we have that $\bm{w}_{+}=\bm{P}\bm{w}^{\mathcal{B}}+O\left(1/p\right)$. Since
$\bm{w}_{+}$ is comprised of the positive components of $\bm{w}_{p}$, we
conclude that

(109) $\displaystyle\bm{w}_{p}=\bm{v}^{\mathcal{B}}+O\left(1/p\right)$

where $\bm{v}^{\mathcal{B}}\in\mathbb{R}^{n}$ is the vector satisfying
$\left(\bm{v}^{\mathcal{B}}\right)_{\mathcal{B}}=0$ and
$\left(\bm{v}^{\mathcal{B}}\right)_{\mathcal{I}}=\bm{P}\bm{u}^{\mathcal{B}}$.
From equation (109) there exists a constant $c_{\mathcal{B}}\geq 0$ such that
$\left\|\bm{w}_{p}-\bm{v}^{\mathcal{B}}\right\|\leq\frac{c_{\mathcal{B}}}{p}$.
Moreover, since the number of subsets of $\\{1,2,\ldots,n\\}$ is finite,
taking $c:=\max_{\mathcal{D}\subseteq\\{1,\ldots,n\\}}c_{\mathcal{D}}$ yields
$\left\|\bm{w}_{p}-\bm{v}^{\mathcal{B}}\right\|\leq\frac{c}{p}$. As $p\geq 1$
was taken to be arbitrary, we conclude that this inequality holds for all
$p\geq 1$.

Finally, we claim that $\bm{v}^{\mathcal{B}}=\bm{w}_{\infty}$, for $p$
sufficiently large. By way of contradiction, suppose that there exists a
$\delta>0$ such that
$\left\|\bm{v}^{\mathcal{B}}-\bm{w}_{\infty}\right\|\geq\delta$, for all $p$
sufficiently large. Since $\bm{w}_{p}\to\bm{w}_{\infty}$ as $p\to\infty$ and
$\left\|\bm{w}_{p}-\bm{v}^{\mathcal{B}}\right\|\leq\frac{c}{p}$ there exists a
$N>0$ such that $\|\bm{w}_{p}-\bm{w}_{\infty}\|<\frac{1}{2}\delta$ and
$\left\|\bm{w}_{p}-\bm{v}^{\mathcal{B}}\right\|<\frac{1}{2}\delta$, for all
$p\geq N$. Hence, for all $p\geq N$, it follows that
$\left\|\bm{v}^{\mathcal{B}}-\bm{w}_{\infty}\right\|\leq\left\|\bm{v}^{\mathcal{B}}-\bm{w}_{p}\right\|+\|\bm{w}_{p}-\bm{w}_{\infty}\|<\delta$,
a contradiction. Since for subsets $\mathcal{B}\subseteq\\{1,2,\ldots,n\\}$
there are only finitely many such vectors $\bm{v}^{\mathcal{B}}$, we have that
$\bm{v}^{\mathcal{B}}=\bm{w}_{\infty}$ for all $p$ sufficiently large.
Therefore, $\|\bm{w}_{p}-\bm{w}_{\infty}\|\leq\frac{c}{p}$ for sufficiently
large $p$.

###  5.2 Convergence Analysis for Constraint Step

In this section, we provide a convergence result for the constraint step in
the LS algorithm. In particular, we establish a linear convergence rate for
the constraint step scheme. Recall the iterative scheme for determining
$\bm{w}\in\Omega$ such that $\bm{h}(\bm{w})=\bm{0}$ used in the constraint
step of the LS algorithm:

(110) $\displaystyle(\bm{w}_{i+1},\bm{y}_{i+1})$
$\displaystyle=\arg\min\left\\{\|\bm{w}-\bm{w}_{i}\|^{2}+\|\bm{y}\|^{2}:\nabla\bm{h}(\bm{w}_{i})(\bm{w}-\bm{w}_{i})+\frac{1}{\sqrt{p_{i}}}\bm{y}=-\bm{h}(\bm{w}_{i}),\bm{w}\in\Omega\right\\}$
(111) $\displaystyle\bm{w}_{i+1}$
$\displaystyle\leftarrow\bm{w}_{i}+s_{i}(\bm{w}_{i+1}-\bm{w}_{i})$

where $p_{i}\geq 1$ is a penalty parameter and where $s_{i}\leq 1$ is chosen
such that

(112)
$\displaystyle\|\bm{h}(\bm{w}_{i}+s_{i}(\bm{w}_{i+1}-\bm{w}_{i}))\|\leq(1-\tau\alpha_{i}s_{i})\|\bm{h}(\bm{w}_{i})\|.$

Here, $\alpha_{i}$ is defined by $\alpha_{i}:=1-\|\bm{y}_{i}\|$ and to
continue performing the scheme we require that $\alpha_{i}\in[\alpha,1]$ for
some fixed parameter $\alpha>0$. Note that this is the iterative scheme in the
constraint step of the LS algorithm, Algorithm 2, where we have made the
change of variables $\bm{y}\leftarrow\frac{1}{\sqrt{p_{i}}}\bm{y}$ to simplify
the analysis in Lemma 5.3. For problems in (110) where the constraint set is
nonempty for $\bm{y}=\bm{0}$, sufficiently large choices of $p$ will result in
$\|\bm{y}\|<1-\alpha$, or equivalently $\alpha_{i}\geq\alpha$. For problems in
(110) where the constraint set is empty for $\bm{y}=\bm{0}$ we may not be able
to satisfy the requirement $\alpha_{i}>0$, let alone $\alpha_{i}\geq\alpha$.
Recall that in such a case NPASA branches to phase one and performs a global
step. An approach similar to (110) – (111) was previously considered where
iterates are given by

(113) $\displaystyle\bm{w}_{i+1}$
$\displaystyle=\arg\min\left\\{\|\bm{w}-\bm{w}_{i}\|^{2}:\nabla\bm{h}(\bm{w}_{i})(\bm{w}-\bm{w}_{i})=-\bm{h}(\bm{w}_{i}),\bm{w}\in\Omega\right\\}$
(114) $\displaystyle\bm{w}_{i+1}$
$\displaystyle\leftarrow\bm{w}_{i}+s_{i}(\bm{w}_{i+1}-\bm{w}_{i})$

and $s_{i}\leq 1$ is chosen using (112). Motivation and discussion of the
scheme in (113) – (114) can be found in [10, 17, 33]. In particular, in
Section 3 of [10] Daniel considered the solvability of the constraint
$\nabla\bm{h}(\bm{w}_{i})(\bm{w}-\bm{w}_{i})=-\bm{h}(\bm{w}_{i})$ found in
(113). By introducing the slack variable, $\bm{y}$, into our update scheme, it
is no longer necessary to consider the solvability of this problem as the
constraint set is always nonempty. However, it is necessary to perform updated
convergence analysis for scheme (113) – (114) to account for the presence of
the slack variable $\bm{y}$ and penalty parameter $p_{i}$ in (110).

Before presenting this result, we consider some heuristics for the constraint
step that should encourage global convergence of NPASA. One of the updates in
this implementation of the primal-dual method from [17] is the use of the
Newton step with feasibility detection in the constraint step of the
algorithm. Previously, the convergence analysis required the assumption that
the linearized constraints in (113) had a solution at each iterate $i$. While
we no longer require this assumption, we do require some method for deciding
when a perturbation $\|\bm{y}_{i+1}\|$ in (110) is small enough for us to
continue with the constraint step iteration. From our analysis in Subection
5.1, in particular inequality (98), we have that if the linearized constraints
are feasible then $\|\bm{y}_{i+1}\|=O\left(\frac{1}{\sqrt{p_{i}}}\right)$.
With this in mind, suppose now that we fix two parameters $\alpha\in(0,1]$ and
$\beta\geq 1$ such that $\frac{1}{\beta}<1-\alpha$. Then we choose our penalty
parameter $p_{i}=\max\left\\{\beta^{2},\|\bm{h}(\bm{w}_{i})\|^{-2}\right\\}$
and consider the linearized constraints to be feasible when
$\|\bm{y}_{i+1}\|\leq 1-\alpha$. As $p_{i}\geq\beta^{2}$, if the linearized
constraint set is feasible without perturbation, then under the updated scheme
in (110) we have that
$\|\bm{y}_{i+1}\|=O\left(\frac{1}{\sqrt{p_{i}}}\right)=O\left(\frac{1}{\beta}\right)$.
By our choice of $\alpha$ and $\beta$, we now have a good chance of satisfying
the criterion $\|\bm{y}_{i+1}\|\leq 1-\alpha$ when the linearized constraint
is feasible without perturbation. If it is infeasible by a perturbation
$\gamma$ in norm, then from inequality (98) we have that
$\|\bm{y}_{i+1}\|^{2}\leq\frac{c}{p_{i}}+\gamma$. Hence, for small
perturbations we still may be able to accept iterates and move towards
satisfying the constraint step.

######  Lemma 5.3.

Suppose that there exists a constant $M$ such that
$\|\bm{w}_{i+1}-\bm{w}_{i}\|\leq M$ for every $i$ and there exists a constant
$r>0$ such that $\nabla\bm{h}$ is Lipschitz continuous with modulus $\kappa$
in the ball $B(\bm{w}_{i},r)$, for each $i$. Then for each $i$ where
$\bm{h}(\bm{w}_{i})\neq\bm{0}$ and $\|\bm{y}_{i}\|\leq 1-\alpha$, the iterate
$\bm{w}_{i+1}$ defined in (110) – (111) exists and satisfies

(115)
$\displaystyle\|\bm{h}(\bm{w}_{i+1})\|\leq(1-\gamma_{i}\alpha_{i}\tau)\|\bm{h}(\bm{w}_{i})\|,$

where $\alpha_{i}=1-\|\bm{y}_{i}\|$ and

(116)
$\displaystyle\gamma_{i}=\min\left\\{1,\frac{r\sigma}{M},\frac{2\alpha\sigma(1-\tau)\|\bm{h}(\bm{w}_{i})\|}{\kappa
M^{2}}\right\\}.$

Thus, either $\bm{h}(\bm{w}_{i})=\bm{0}$ after a finite number of iterations
or $\displaystyle\lim_{i\to\infty}\bm{h}(\bm{w}_{i})=\bm{0}$.

######  Proof 5.4.

First note that by the definition of $\alpha_{i}$ and our hypothesis that
$\|\bm{y}_{i}\|\leq 1-\alpha$ it follows that $\alpha\leq\alpha_{i}$. Now let
$\bm{w},\bm{d}\in\mathbb{R}^{n}$ and $s\in\mathbb{R}$ such that $\bm{h}$ is
$r$-Lipschitz continuous on the line segment $[\bm{w},\bm{w}+s\bm{d}]$. By the
fundamental theorem of calculus it follows that

(117) $\displaystyle\bm{h}(\bm{w}+s\bm{d})$
$\displaystyle=\bm{h}(\bm{w})+s\nabla\bm{h}(\bm{w})\bm{d}+\int_{0}^{s}\left(\nabla\bm{h}(\bm{w}+t\bm{d})-\nabla\bm{h}(\bm{w})\right)\bm{d}\
dt.$

By (110), we have that

(118)
$\displaystyle\nabla\bm{h}(\bm{w}_{i})(\bm{w}_{i+1}-\bm{w}_{i})=-\bm{h}(\bm{w}_{i})-\frac{1}{\sqrt{p}}\bm{y}_{i}.$

Hence, setting $\bm{w}=\bm{w}_{i}$, $\bm{y}=\bm{y}_{i}$, and
$\bm{d}=\bm{w}_{i+1}-\bm{w}_{i}$ in (117) and combining the resulting equation
with (118) we have that

(119) $\displaystyle\|\bm{h}(\bm{w}_{i}+s(\bm{w}_{i+1}-\bm{w}_{i}))\|$
$\displaystyle\leq\|\bm{h}(\bm{w}_{i})+s\nabla\bm{h}(\bm{w}_{i})(\bm{w}_{i+1}-\bm{w}_{i})\|+\frac{\kappa
s^{2}M^{2}}{2}$ (120)
$\displaystyle\leq\|\bm{h}(\bm{w}_{i})-s\bm{h}(\bm{w}_{i})\|+\frac{s}{\sqrt{p}}\|\bm{y}_{i}\|+\frac{\kappa
s^{2}M^{2}}{2}$ (121)
$\displaystyle\leq(1-s)\|\bm{h}(\bm{w}_{i})\|+s\|\bm{y}_{i}\|\|\bm{h}(\bm{w}_{i})\|+\frac{\kappa
s^{2}M^{2}}{2}$ (122)
$\displaystyle=\left(1-s\left(1-\|\bm{y}_{i}\|\right)\right)\|\bm{h}(\bm{w}_{i})\|+\frac{\kappa
s^{2}M^{2}}{2}$ (123)
$\displaystyle=\left(1-\alpha_{i}s\right)\|\bm{h}(\bm{w}_{i})\|+\frac{\kappa
s^{2}M^{2}}{2}.$

Now assume that $\bm{h}(\bm{w}_{i})\neq\bm{0}$. By our choice of $s_{i}$ in
NPASA, we have that either $s_{i}=1$ or

(124)
$\displaystyle(1-\tau\alpha_{i}\sigma^{-1}s_{i})\|\bm{h}(\bm{w}_{i})\|\leq\|\bm{h}(\bm{w}_{i}-\sigma^{-1}s_{i}(\bm{w}_{i+1}-\bm{w}_{i}))\|.$

If $s_{i}=1$ then there is nothing to be done. Hence, we suppose that (124) is
satisfied and consider two subcases.
_Case_ ($i$): $\sigma^{-1}s_{i}\|\bm{w}_{i+1}-\bm{w}_{i}\|>r$. By our
hypothesis that $\|\bm{w}_{i+1}-\bm{w}_{i}\|\leq M$ we have the lower bound

(125) $\displaystyle
s_{i}>\frac{r\sigma}{\|\bm{w}_{i+1}-\bm{w}_{i}\|}\geq\frac{r\sigma}{M}.$

_Case_ ($ii$): $\sigma^{-1}s_{i}\|\bm{w}_{i+1}-\bm{w}_{i}\|\leq r$. First,
note that by setting $s=\sigma^{-1}s_{i}$ in (123) we have that

(126)
$\displaystyle\|\bm{h}(\bm{w}_{i}+\sigma^{-1}s_{i}(\bm{w}_{i+1}-\bm{w}_{i}))\|$
$\displaystyle\leq\left(1-\alpha_{i}\sigma^{-1}s_{i}\right)\|\bm{h}(\bm{w}_{i})\|+\frac{\kappa
s_{i}^{2}M^{2}}{2\sigma^{2}}.$

Combining (124) with (126) yields
$(1-\tau\alpha_{i}\sigma^{-1}s_{i})\|\bm{h}(\bm{w}_{i})\|\leq\left(1-\alpha_{i}\sigma^{-1}s_{i}\right)\|\bm{h}(\bm{w}_{i})\|+\frac{\kappa
s_{i}^{2}M^{2}}{2\sigma^{2}}$, which, when simplified, gives us

(127)
$\displaystyle\frac{2\alpha_{i}\sigma(1-\tau)\|\bm{h}(\bm{w}_{i})\|}{\kappa
M^{2}}$ $\displaystyle\leq s_{i}.$

Lastly, note that by the definition of $\alpha_{i}$ and our hypothesis that
$\|\bm{y}_{i}\|\leq 1-\alpha$ it follows that $\alpha\leq\alpha_{i}$. Thus, by
(5.4), (125), and (127) we conclude that the choice of $\gamma_{i}$ stated in
the lemma satisfies the Armijo line search at iterate $i$.

###  5.3 Convergence Analysis for Multiplier Step

We now provide a convergence result for the multiplier step in the LS
algorithm. We note that this result is similar to Lemma 6.2 in [17] and we
omit the proof as it follows from the same argument provided in [17] when the
function $\bm{K}(\bm{\Lambda},\bm{x})$ is replaced with
$E_{m,1}(\bm{z},\bm{\nu},\bm{\eta})$.

######  Lemma 5.5.

Suppose that $f,\bm{h}\in\mathcal{C}^{1}(\Omega)$ and that $\Omega$ is
compact. Then there exists a sequence of natural numbers
$\\{n_{i}\\}_{i=0}^{\infty}$ such that $\\{\bm{z}_{n_{i}}\\}_{i=0}^{\infty}$
is a convergent subsequence with limit point $\bm{z}^{*}$. Furthermore, if
(LICQ) holds at $\bm{z}^{*}$ then there exists a vector $\bm{\nu}^{*}$ such
that $\min_{\bm{\eta}\geq\bm{0}}\
E_{m,1}(\bm{z}^{*},\bm{\nu}^{*},\bm{\eta})=0$.

We note that Lemma 5.5 could be used to modify the multiplier step of the LS
Algorithm, Algorithm 2, so that it is designed to identify a convergent
subsequence instead of a convergent sequence. However, as phase two of NPASA
is developed for fast local convergence we have intentionally designed the
multiplier step to continue only when the multiplier error, $E_{m,1}$, is
decreasing at a linear rate or higher. This concludes our convergence analysis
for phase two of NPASA.

##  6 Convergence Analysis of NPASA

At this point, we have established convergence results for each problem solved
in NPASA and we are ready to state and prove a global convergence result for
NPASA.

######  Theorem 6.1 (NPASA Global Convergence Theorem).

Suppose that $\Omega$ is compact and that
$f,\bm{h}\in\mathcal{C}^{1}(\Omega)$. Suppose that NPASA (Algorithm 3) with
$\varepsilon=0$ generates a sequence
$\\{(\bm{x}_{k},\bm{\lambda}_{k},\bm{\mu}_{k})\\}_{k=0}^{\infty}$ with
$\bm{x}_{k}\in\Omega$ and let $\mathcal{S}_{j}$ be the set of indices such
that if $k\in\mathcal{S}_{j}$ then $\bm{x}_{k}$ is generated in phase $j$ of
NPASA, for $j\in\\{1,2\\}$. Suppose that the following assumptions hold:

  1. _(H1)_

For every $k\in\mathcal{S}_{1}$, assumptions (G1) and (G2) in Section 4 hold
at $(\bm{x}_{k-1},\bm{\nu}_{k-1},q_{k})$, where
$\bm{\nu}_{k-1}=Proj_{[-\bar{\lambda},\bar{\lambda}]}(\bm{\lambda}_{k-1})$ and
$\bar{\lambda}>0$ is a scalar parameter.

  2. _(H2)_

If $\bm{x}^{*}$ is a subsequential limit point of
$\\{\bm{x}_{k}\\}_{k=0}^{\infty}$ then (LICQ) holds at $\bm{x}^{*}$.

Then

(128)
$\displaystyle\liminf_{k\to\infty}E_{1}(\bm{x}_{k},\bm{\lambda}_{k},\bm{\mu}_{k})=0.$

######  Proof 6.2.

To prove this result we consider three different scenarios. Before considering
each scenario, we claim that the hypotheses of Theorem 4.3 are satisfied
across all cases. To see this, from assumption (H1) it follows that Theorem
4.9 holds. By Theorem 4.9, Algorithm 1 can always generate a point satisfying
hypothesis (1) in Theorem 4.3. Additionally, (2) and (3) in Theorem 4.3 are
satisfied by the definition of Algorithm 3. The final assumption in Theorem
4.3 follows from (H2). Hence, all of the hypotheses of Theorem 4.3 are
satisfied. We now consider the three cases.
Case 1: $\mathcal{S}_{1}$ is finite. In this case, there exists an nonnegative
integer $N$ such that $k\geq N$ implies $k\in\mathcal{S}_{2}$. As the
branching criterion for phase two to phase one in line 42 of Algorithm 3 are
never satisfied for $k>N$, we have that
$E_{1}(\bm{x}_{k},\bm{\lambda}_{k},\bm{\mu}_{k})\leq\theta
E_{1}(\bm{x}_{k-1},\bm{\lambda}_{k-1},\bm{\mu}_{k-1})$, for all $k>N$. In
particular, we have
$E_{1}(\bm{x}_{N+k},\bm{\lambda}_{N+k},\bm{\mu}_{N+k})\leq\theta^{k}E_{1}(\bm{x}_{N},\bm{\lambda}_{N},\bm{\mu}_{N})$,
for all $k\geq 0$. As $\theta\in(0,1)$, it follows that
$\lim_{k\to\infty}E_{1}(\bm{x}_{k},\bm{\lambda}_{k},\bm{\mu}_{k})=0$.
Case 2: $\mathcal{S}_{2}$ is finite. In this case, there exists an nonnegative
integer $N$ such that $k\geq N$ implies $k\in\mathcal{S}_{1}$. Now suppose
phase one of NPASA is given the starting point
$(\bm{x}_{N},\bm{\lambda}_{N},\bm{\mu}(\bm{x}_{N},1))$. Then phase one of
NPASA generates all remaining iterates and, as the hypotheses of Theorem 4.3
are satisfied, it follows that
$E_{1}(\bm{x}^{*},\bm{\lambda}^{*},\bm{\mu}(\bm{x}^{*},1))=0$. As $\bm{x}^{*}$
is a limit point of a subsequence of $\\{\bm{x}_{k}\\}_{k=0}^{\infty}$, it
follows that
$\liminf_{k\to\infty}E_{1}(\bm{x}_{k},\bm{\lambda}_{k},\bm{\mu}(\bm{x}_{k},1))=E_{1}(\bm{x}^{*},\bm{\lambda}^{*},\bm{\mu}(\bm{x}^{*},1))=0$,
which proves (128).
Case 3: $\mathcal{S}_{1}$ and $\mathcal{S}_{2}$ are infinite. Suppose that
$k\in\mathcal{S}_{1}$. By (H1), phase one will always generate a point that
satisfies the criterion for branching to phase two, that is
$E_{m,0}(\bm{x}_{k},\bm{\lambda}_{k},\bm{\mu}_{k})\leq\theta
E_{c}(\bm{x}_{k})$. Combining this with
$E_{m,1}(\bm{x}_{k},\bm{\lambda}_{k},\bm{\mu}_{k})\leq
E_{m,0}(\bm{x}_{k},\bm{\lambda}_{k},\bm{\mu}_{k})$ from Lemma 2.1 yields

(129) $\displaystyle
E_{m,1}(\bm{x}_{k},\bm{\lambda}_{k},\bm{\mu}_{k})\leq\theta
E_{c}(\bm{x}_{k}).$

Now let $\\{n_{1},n_{2},\ldots\\}\subset\mathcal{S}_{1}$ be the indices such
that $\bm{x}_{n_{k}}$ is a term in the subsequence converging to $\bm{x}^{*}$.
From equation (56) in the proof of Theorem 4.3, there exists a constant
$\Lambda$ such that

(130) $\displaystyle\|\bm{h}(\bm{x}_{n_{k}})\|\leq\frac{\Lambda}{2q_{n_{k}}},$

for all $k$ sufficiently large. By the update formula for $q_{k}$ in line 4 of
Algorithm 3 and the fact that $\phi>1$, we have that $q_{k}\to\infty$ as
$k\to\infty$. Hence, for $k$ sufficiently large, $q_{n_{k}}\geq\Lambda$. In
particular, by the update formula for $q_{k}$ we have that

(131) $\displaystyle
q_{n_{k}}\geq\frac{q_{n_{k-1}}}{e_{n_{k}-1}}\geq\frac{\Lambda}{e_{n_{k}-1}},$

for all $k$ sufficiently large. As $E_{c}(\bm{x})=\|\bm{h}(\bm{x})\|^{2}$, it
now follows from (129), (130), and (131) that

(132) $\displaystyle
E_{1}(\bm{x}_{n_{k}},\bm{\lambda}_{n_{k}},\bm{\mu}_{n_{k}})\leq\sqrt{(1+\theta)E_{c}(\bm{x}_{n_{k}})}\leq\sqrt{1+\theta}\left(\frac{\Lambda}{2q_{n_{k}}}\right)\leq\left(\frac{\sqrt{1+\theta}}{2}\right)e_{n_{k}-1}.$

As $\displaystyle e_{n_{k}-1}=\min_{0\leq
i<n_{k}}E_{1}(\bm{x}_{i},\bm{\lambda}_{i},\bm{\mu}_{i})$, it follows from
(132) that

(133) $\displaystyle
E_{1}(\bm{x}_{n_{k}},\bm{\lambda}_{n_{k}},\bm{\mu}_{n_{k}})\leq\left(\frac{\sqrt{1+\theta}}{2}\right)E_{1}(\bm{x}_{n_{k-1}},\bm{\lambda}_{n_{k-1}},\bm{\mu}_{n_{k-1}}),$

for all $k$ sufficiently large.

On the other hand, suppose that $k\in\mathcal{S}_{2}$. By the branching
criterion for phase two, we have that

(134) $\displaystyle E_{1}(\bm{x}_{k},\bm{\lambda}_{k},\bm{\mu}_{k})$
$\displaystyle\leq\theta
E_{1}(\bm{x}_{k-1},\bm{\lambda}_{k-1},\bm{\mu}_{k-1}).$

As $\theta\in(0,1)$ by definition, we also have that
$(\sqrt{1+\theta})/2\in(0,1)$. Thus, combining (133) and (134) yields (128).

##  7 Conclusion

In this paper, we presented a method for solving nonlinear programs, NPASA,
and established global convergence properties for NPASA. In particular, under
a small set of assumptions we established global convergence properties for
each of the three subproblems solved during NPASA which, when combined, were
able to ensure global convergence of NPASA. As noted in the introduction, a
companion paper [11] focuses on establishing local quadratic convergence for
NPASA. In the future, we plan on implementing this approach and performing
benchmarking tests to compare the NPASA algorithm to leading methods for
solving nonlinear programs.

This work was performed, in part, under the auspices of the U.S. Department of
Energy by Lawrence Livermore National Laboratory under Contract DE-
AC52-07NA27344.

## References

  * [1] Decision tree for optimization software.  http://plato.asu.edu/bench.html.  Accessed: 2019-06-24.
  * [2] R. Andreani, E. Birgin, J. Martínez, and M. Schuverdt, On augmented lagrangian methods with general lower-level constraints, SIAM Journal on Optimization, 18 (2008), pp. 1286–1309, https://doi.org/10.1137/060654797, https://doi.org/10.1137/060654797, https://arxiv.org/abs/https://doi.org/10.1137/060654797.
  * [3] D. Bertsekas, Nonlinear Programming, Athena Scientific, 1995, https://books.google.com/books?id=QeweAQAAIAAJ.
  * [4] E. G. Birgin, R. A. Castillo, and J. M. MartÍnez, Numerical comparison of augmented lagrangian algorithms for nonconvex problems, Computational Optimization and Applications, 31 (2005), pp. 31–55, https://doi.org/10.1007/s10589-005-1066-7, https://doi.org/10.1007/s10589-005-1066-7.
  * [5] S. Bochkanov, ALGLIB.  http://www.alglib.net.  Accessed: 2019-05-07.
  * [6] R. H. Byrd, J. Nocedal, and R. A. Waltz, Knitro: An Integrated Package for Nonlinear Optimization, Springer US, Boston, MA, 2006, pp. 35–59, https://doi.org/10.1007/0-387-30065-1_4, https://doi.org/10.1007/0-387-30065-1_4.
  * [7] A. Conn, N. Gould, and P. Toint, A globally convergent augmented lagrangian algorithm for optimization with general constraints and simple bounds, SIAM Journal on Numerical Analysis, 28 (1991), pp. 545–572, https://doi.org/10.1137/0728030, https://doi.org/10.1137/0728030, https://arxiv.org/abs/https://doi.org/10.1137/0728030.
  * [8] A. R. Conn, N. Gould, and P. L. Toint, Numerical experiments with the lancelot package (release a) for large-scale nonlinear optimization, Mathematical Programming, 73 (1996), p. 73, https://doi.org/10.1007/BF02592099, https://doi.org/10.1007/BF02592099.
  * [9] A. R. Conn, N. I. M. Gould, and P. L. Toint, Lancelot: A Fortran Package for Large-Scale Nonlinear Optimization (Release A), Springer Publishing Company, Incorporated, 1st ed., 2010.
  * [10] J. W. Daniel, Newton’s method for nonlinear inequalities, Numerische Mathematik, 21 (1973), pp. 381–387, https://doi.org/10.1007/BF01436488, https://doi.org/10.1007/BF01436488.
  * [11] J. Diffenderfer and W. W. Hager, NPASA: An algorithm for nonlinear programming - local convergence, 2021.
  * [12] R. Fletcher and S. Leyffer, User manual for filtersqp, Technical Report NA/181, (1998).
  * [13] R. Frisch, The Logarithmic Potential Method of Convex Programming, with Particular Application to the Dynamics of Planning for National Development: Synopsis of a Communication Etc, Memoranda, Socialøkonomiske Institutt, 1955, https://books.google.com/books?id=Qq3kswEACAAJ.
  * [14] P. E. Gill, W. Murray, and M. A. Saunders, Snopt: An sqp algorithm for large-scale constrained optimization, SIAM Review, 47 (2005), pp. 99–131, https://doi.org/10.1137/S0036144504446096, http://dx.doi.org/10.1137/S0036144504446096, https://arxiv.org/abs/http://dx.doi.org/10.1137/S0036144504446096.
  * [15] P. E. Gill, W. Murray, M. A. Saunders, and E. Wong, User’s guide for SNOPT 7.7: Software for large-scale nonlinear programming, Center for Computational Mathematics Report CCoM 18-1, Department of Mathematics, University of California, San Diego, La Jolla, CA, 2018.
  * [16] P. E. Gill, M. A. Saunders, and E. Wong, User’s guide for DNOPT: Software for medium-scale nonlinear programming, (2017).
  * [17] W. W. Hager, Analysis and implementation of a dual algorithm for constrained optimization, Journal of Optimization Theory and Applications, 79 (1993), pp. 427–462, https://doi.org/10.1007/BF00940552, http://dx.doi.org/10.1007/BF00940552.
  * [18] W. W. Hager and D. Mico-Umutesi, Error estimation in nonlinear optimization, Journal of Global Optimization, 59 (2014), pp. 327–341, https://doi.org/10.1007/s10898-014-0186-y, http://dx.doi.org/10.1007/s10898-014-0186-y.
  * [19] W. W. Hager and H. Zhang, A new active set algorithm for box constrained optimization, SIAM J. on Optimization, 17 (2006), pp. 526–557, https://doi.org/10.1137/050635225, http://dx.doi.org/10.1137/050635225.
  * [20] W. W. Hager and H. Zhang, An active set algorithm for nonlinear optimization with polyhedral constraints, Sci. China Math., 59 (2016), pp. 1525–1542.
  * [21] W. W. Hager and H. Zhang, Projection onto a polyhedron that exploits sparsity, SIAM Journal on Optimization, 26 (2016), pp. 1773–1798, https://doi.org/10.1137/15M102825X, http://dx.doi.org/10.1137/15M102825X, https://arxiv.org/abs/http://dx.doi.org/10.1137/15M102825X.
  * [22] S.-P. Han, Superlinearly convergent variable metric algorithms for general nonlinear programming problems, Mathematical Programming, 11 (1976), pp. 263–282, https://doi.org/10.1007/BF01580395, https://doi.org/10.1007/BF01580395.
  * [23] M. R. Hestenes, Multiplier and gradient methods, Journal of Optimization Theory and Applications, 4 (1969), pp. 303–320, https://doi.org/10.1007/BF00927673, https://doi.org/10.1007/BF00927673.
  * [24] A. Izmailov, M. Solodov, and E. Uskov, Global convergence of augmented lagrangian methods applied to optimization problems with degenerate constraints, including problems with complementarity constraints, SIAM Journal on Optimization, 22 (2012), pp. 1579–1606, https://doi.org/10.1137/120868359, https://doi.org/10.1137/120868359, https://arxiv.org/abs/https://doi.org/10.1137/120868359.
  * [25] R. Lewis and V. Torczon, A globally convergent augmented lagrangian pattern search algorithm for optimization with general constraints and simple bounds, SIAM Journal on Optimization, 12 (2002), pp. 1075–1089, https://doi.org/10.1137/S1052623498339727, https://doi.org/10.1137/S1052623498339727, https://arxiv.org/abs/https://doi.org/10.1137/S1052623498339727.
  * [26] H. Luo, X. Sun, and D. Li, On the convergence of augmented lagrangian methods for constrained global optimization, SIAM Journal on Optimization, 18 (2008), pp. 1209–1230, https://doi.org/10.1137/060667086, https://doi.org/10.1137/060667086, https://arxiv.org/abs/https://doi.org/10.1137/060667086.
  * [27] B. A. Murtagh and M. A. Saunders, Large-scale linearly constrained optimization, Mathematical Programming, 14 (1978), pp. 41–72, https://doi.org/10.1007/BF01588950, https://doi.org/10.1007/BF01588950.
  * [28] J. Nocedal and S. J. Wright, Numerical optimization, Springer Series in Operations Research and Financial Engineering, Springer, Berlin, 2006, http://opac.inria.fr/record=b1120179.  NEOS guide http://www-fp.mcs.anl.gov/otc/Guide/.
  * [29] U. M. G. Palomares and O. L. Mangasarian, Superlinearly convergent quasi-newton algorithms for nonlinearly constrained optimization problems, Mathematical Programming, 11 (1976), pp. 1–13, https://doi.org/10.1007/BF01580366, https://doi.org/10.1007/BF01580366.
  * [30] E. Polak and A. L. Tits, A globally convergent, implementable multiplier method with automatic penalty limitation, Applied Mathematics and Optimization, 6 (1980), pp. 335–360, https://doi.org/10.1007/BF01442901, https://doi.org/10.1007/BF01442901.
  * [31] M. J. D. Powell, A method for nonlinear constraints in minimization problems, in Optimization, R. Fletcher, ed., Academic Press, New York, 1969, pp. 283–298.
  * [32] M. J. D. Powell, A fast algorithm for nonlinearly constrained optimization calculations, in Numerical Analysis, G. A. Watson, ed., Berlin, Heidelberg, 1978, Springer Berlin Heidelberg, pp. 144–157.
  * [33] S. M. Robinson, Extension of newton’s method to nonlinear functions with values in a cone, Numerische Mathematik, 19 (1972), pp. 341–347, https://doi.org/10.1007/BF01404880, https://doi.org/10.1007/BF01404880.
  * [34] A. Wächter and L. Biegler, Line search filter methods for nonlinear programming: Local convergence, SIAM Journal on Optimization, 16 (2005), pp. 32–48, https://doi.org/10.1137/S1052623403426544, https://doi.org/10.1137/S1052623403426544, https://arxiv.org/abs/https://doi.org/10.1137/S1052623403426544.
  * [35] A. Wächter and L. Biegler, Line search filter methods for nonlinear programming: Motivation and global convergence, SIAM Journal on Optimization, 16 (2005), pp. 1–31, https://doi.org/10.1137/S1052623403426556, https://doi.org/10.1137/S1052623403426556, https://arxiv.org/abs/https://doi.org/10.1137/S1052623403426556.
  * [36] A. Wächter and L. T. Biegler, On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming, Math. Program., 106 (2006), pp. 25–57, https://doi.org/10.1007/s10107-004-0559-y, https://doi.org/10.1007/s10107-004-0559-y.
  * [37] R. Wilson, A Simplicial Algorithm for Concave Programming, Graduate School of Business Administration, George F. Baker Foundation, Harvard University, 1963, https://books.google.com/books?id=Ec4oYAAACAAJ.
  * [38] H. Wu, H. Luo, and S. Li, The global convergence of augmented lagrangian methods based on ncp function in constrained nonconvex optimization, Applied Mathematics and Computation, 207 (2009), pp. 124 – 134, https://doi.org/https://doi.org/10.1016/j.amc.2008.10.015, http://www.sciencedirect.com/science/article/pii/S0096300308007716.  Includes Special issue on Emergent Applications of Fractals and Wavelets in Biology and Biomedicine.

##  Appendix A Constrained Optimization Definitions and Results

In this appendix, we highlight key definitions and results from the
constrained optimization theory that were used as assumptions in our analysis
of NPASA. Statements of these results are included here for completeness but
we note that more details on these theorems and proofs can be found by
referencing [28, 3]. As the hypotheses of many of these results are used as
assumptions in establishing global and local convergence properties for NPASA,
these results are only referenced in Section 1.1 where we provide simplified
abbreviations for the hypotheses of these theorems. In this appendix, given
$i\in\mathbb{N}$ note that we will write $[i]$ to denote the set of integers
$\\{1,2,\ldots,i\\}$.

######  Definition A.1.

A vector $\bm{x}$ is said to satisfy the linear independence constraint
qualification condition for problem (3) if the matrix
$\displaystyle\begin{bmatrix}\nabla\bm{h}(\bm{x})\\\
\bm{A}_{\mathcal{A}(\bm{x})}\end{bmatrix}$ is of full row rank.

As in Section 1.1, note that if $\bm{x}$ satisfies the linear independence
constraint qualification condition then we abbreviate this by writing _( LICQ)
holds at $\bm{x}$_. We now state the first order optimality or KKT conditions.

######  Theorem A.2 (Karush-Kuhn-Tucker Conditions [28]).

Suppose $\bm{x}^{*}$ is a local solution of problem (3), that
$f,\bm{h}\in\mathcal{C}^{1}$, and (LICQ) holds at $\bm{x}^{*}$. Then there
exist KKT multipliers $\bm{\lambda}^{*}$ and $\bm{\mu}^{*}$ such that the
following conditions hold at $(\bm{x}^{*},\bm{\lambda}^{*},\bm{\mu}^{*})$:

  1. (KKT.1)

Gradient of Lagrangian equals zero:
$\nabla_{x}\mathcal{L}(\bm{x}^{*},\bm{\lambda}^{*},\bm{\mu}^{*})^{\intercal}=\bm{0}$

  2. (KKT.2)

Satisfies equality constraints: $\bm{h}(\bm{x}^{*})=\bm{0}$

  3. (KKT.3)

Satisfies inequality constraints: $\bm{r}(\bm{x}^{*})\leq\bm{0}$

  4. (KKT.4)

Nonnegativity of inequality multipliers: $\bm{\mu}^{*}\geq\bm{0}$

  5. (KKT.5)

Complementary slackness: $r_{i}(\bm{x}^{*})\mu_{i}^{*}=0$ for $i\in[m]$.

######  Definition A.3.

A point $(\bm{x},\bm{\lambda},\bm{\mu})$ is said to satisfy strict
complementary slackness for problem (3) if it satisfies (KKT.1) – (KKT.5) and
exactly one of $r_{i}(\bm{x})$ and $\mu_{i}$ is zero for each $i\in[m]$.