# Constrained Optimization in the Presence of Noise

Figen Oztoprak Artelys Corporation Richard Byrd Computer Science Department,
University of Colorado, Boulder, USA  Jorge Nocedal Department of Industrial
Engineering and Management Sciences, Northwestern University, USA. This author
was supported by National Science Foundation grant DMS-2011494, AFOSR grant
FA95502110084, and ONR grant N00014-21-1-2675.

###### Abstract

The problem of interest is the minimization of a nonlinear function subject to
nonlinear equality constraints using a sequential quadratic programming (SQP)
method. The minimization must be performed while observing only noisy
evaluations of the objective and constraint functions. In order to obtain
stability, the classical SQP method is modified by relaxing the standard
Armijo line search based on the noise level in the functions, which is assumed
to be known. Convergence theory is presented giving conditions under which the
iterates converge to a neighborhood of the solution characterized by the noise
level and the problem conditioning. The analysis assumes that the SQP
algorithm does not require regularization or trust regions. Numerical
experiments indicate that the relaxed line search improves the practical
performance of the method on problems involving uniformly distributed noise.
One important application of this work is in the field of derivative-free
optimization, when finite differences are employed to estimate gradients.

##  1 Introduction

Let us consider the equality constrained nonlinear optimization problem

$\min_{x}f(x)\qquad\mbox{s.t.}\ \ c(x)=0,$ (1.1)

where $f:\mathbb{R}^{n}\rightarrow\mathbb{R}$ and
$c(x):\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}$ are smooth functions. We assume
that the minimization must be performed while observing approximate
evaluations $\tilde{f}(x),\tilde{c}(x)$ of the functions $f,c$ and their
derivatives.

We consider the application of a sequential quadratic programming (SQP)
algorithm that employs an $\ell_{1}$ merit function to control the stepsize.
The goal of the paper is to study the effect of noise on the behavior of the
SQP algorithm, particularly the achievable accuracy in the solution, and to
highlight the aspects of the algorithm that are most susceptible to errors (or
noise)—and redesign them. This work was motivated by applications in which the
derivatives of $f$ and $c$ are approximated by finite differences [15], and
thus contain errors, but the algorithm and analysis apply to the more general
setting when stochastic or deterministic noise are present in both the
function and derivative evaluations.

Let us define

$g(x)=\nabla f(x),\qquad J(x)=\nabla c(x)\in\mathbb{R}^{m\times n},\quad m<n,$
(1.2)

and let $\tilde{g}(x),\tilde{J}(x)$ be the corresponding noisy evaluations.
The iteration of the SQP algorithm is given by

$\displaystyle x_{k+1}=x_{k}+\alpha_{k}d_{k},$ (1.3)

where $d_{k}$ is the solution of the quadratic subproblem

$\displaystyle\min_{d\in\mathbb{R}^{n}}$ $\displaystyle\
\tfrac{1}{2}d^{T}H_{k}d+\tilde{g}_{k}^{T}d$ (1.4) $\displaystyle{\rm s.t.}$
$\displaystyle\ \tilde{c}_{k}+\tilde{J}_{k}d=0,$ (1.5)

and the steplength $\alpha_{k}>0$ is chosen so as to ensure sufficient
decrease in the merit function

$\tilde{\phi}(x)=\tilde{f}(x)+\pi\|\tilde{c}(x)\|_{1}$ (1.6)

when the iterates are far away from a solution. Here $\pi>0$ is a penalty
parameter that is adjusted during the course of the optimization. The
symmetric matrix $H_{k}$ is generally chosen as an approximation to the
Hessian of the Lagrangian. However, in this paper we assume that $H_{k}$ is a
multiple of the identity matrix,

$H_{k}=\beta_{k}I\qquad\beta_{k}>0,$ (1.7)

because allowing more general choices introduces more constants in the
analysis without contributing to the main goals of this investigation.

As in the noiseless case, the control of the penalty parameter in (1.6) is of
critical importance in the SQP algorithm. $\pi$ should be chosen so that the
SQP direction $d_{k}$ is a descent direction for $\tilde{\phi}$ at $x_{k}$,
and it should provide adequate control on the size of $\alpha_{k}$. The
proposed algorithm has the general form of a classical SQP method [7],
specialized to the case when $H_{k}$ is a multiple of the identity matrix, and
introduces a modification in the line search designed to handle noise.

We assume throughout that the noise in the function and gradient evaluations
is bounded by some constants $\epsilon_{f}$ and $\epsilon_{c}$. This is not
always the case in practice (e.g. when noise is Gaussian) but it covers many
important practical settings, including computational noise [11]. Furthermore,
we assume that $\epsilon_{f},\epsilon_{c}$ are known, or can be estimated, and
that the algorithm has access to them.

This study was motivated by some practical computations performed by the
authors using the knitro software package [5]. They selected a few challenging
nonlinear optimization problems involving equality and inequality constraints,
injected noise in the objective and constraints, and computed derivatives
using noise-aware finite difference formula; see e.g., Moré and Wild [11] and
Berahas et al. [1]. They observed that, for low levels of noise, knitro
returned acceptable answers, even though one might suspect the default
algorithm to be brittle in this setting. As the noise level was increased, the
quality of the solution deteriorated markedly, suggesting that classical
optimization methods should be redesigned to handle noise. To guide this
investigation, it is essential to develop a convergence theory. In this paper,
we focus on the case when noise cannot be diminished, and characterize the
accuracy of a noise tolerant optimization algorithm.

As a first step in this investigation, we find it convenient to consider
equality constrained optimization, and study the performance of a sequential
quadratic optimization method, which is a simple method in this setting and
must yet confront some important challenges raised by the presence of noise.

###  1.1 Contributions of this work

The main contribution of this paper is the development of a convergence theory
for a classical sequential quadratic programming (SQP) algorithm for equality
constrained optimization in the presence of noise. It is shown that, by
introducing a relaxation in the line search procedure while keeping all other
components of the SQP method unchanged, the iterates of the algorithm reach an
acceptable neighborhood $C_{1}$ of the solution defined by a stationarity
measure for the problem. Furthermore, once the iterates enter $C_{1}$ they
cannot escape a larger neighborhood $C_{2}$ and must revisit $C_{1}$ an
infinite number of times. The analysis gives a detailed characterization of
these neighborhoods in terms of the noise level and problem characteristics.
Numerical experiments show that the relaxed line search is, in fact,
beneficial in practice.

Our convergence results assume that errors in function and gradients are
bounded, and the analysis is deterministic, yielding somewhat pessimistic
bounds. We believe, however, that the results can be useful in the design of
robust constrained optimization methods. Specifically, our analysis suggests
that only slight modifications are needed so that a classical SQP method is
able to handle bounded noise.

###  1.2 Literature Review

Early work on constrained optimization in the presence of noise is reviewed by
Poljak (a.k.a. Polyak) [13]. His study includes penalty, Lagrange, or extended
Lagrange functions, and establishes probabilistic convergence theorems
provided the steplength is chosen small enough from the start. Hintermueller
[10] studies a penalty SQP method in which equality constraints are replaced
with upper and lower bounding surrogates. Assuming that the noise level in the
function is known, it is shown that in the limit the bounds contain a
solution. Schittkowski [14] uses a non-monotone line search to handle errors
due to approximate function and derivative evaluations. His algorithm was
implemented in the NLPQLP software, which is reported to be successful in
practice, but no convergence theory were presented.

The work that is most closely related to this study is [2, 3, 6]. In [3], an
SQP method for equality constrained optimization is presented to handle the
case when the objective function is stochastic and the constraints are
deterministic. The stepsize is obtained by adaptively estimating Lipschitz
constants in place of a line search. Conditions for convergence in expectation
are established. [2] considers the case when Jacobians can be rank deficient,
proposes a step decomposition approach, and presents compelling numerical
results. [6] studies an SQP algorithm with an inexact step computation for the
same problem setting. These three papers give careful attention to the
behavior of the penalty parameter. For example, in [3] the penalty parameter
is chosen in a way that provides sufficient descent in the quadratic model of
the merit function in the deterministic setting. In the stochastic setting,
they employ the stochastic gradient of the objective in the same formulae for
updating the penalty parameter, but they can no longer guarantee that the
resulting penalty parameter will be large enough and bounded. They prove their
convergence results assuming that the penalty parameter is well behaved. Then,
they discuss the probability of having small penalty values, and note that the
boundedness issue is resolved by making the same assumption as in this paper,
namely that noise is always bounded.

Notation. We let $\|\cdot\|$ denote the $\ell_{2}$ norm, unless otherwise
stated. As is the convention, $f_{k}$ stands for $f(x_{k})$ and similarly for
other functions. The terms error and noise in the functions is used
interchangeably. Since we assume absolute bounds on these quantities, the
distinction between them is not important in this study.

##  2 The Algorithm

Before presenting the algorithm, we introduce some notation. We model the
first-order change in the merit function $\phi$ at an iterate $x_{k}$ as

$\tilde{\ell}(x_{k};d_{k})=\tilde{g}_{k}^{T}d_{k}+\pi_{k}\|\tilde{c}_{k}+\tilde{J}_{k}d_{k}\|_{1}-\pi_{k}\|\tilde{c}_{k}\|_{1}.$
(2.1)

We also define

$\hat{\lambda}_{k}=(\tilde{J}_{k}\tilde{J}_{k}^{T})^{-1}\tilde{J}_{k}\tilde{g}_{k},$
(2.2)

which is the standard least squares multiplier estimate [12, eqn(18.21)],
accounting for noisy function evaluations. We assume that $\tilde{J}_{k}$ is
full rank for all $k$, hence $\hat{\lambda}_{k}$ is well defined.

The penalty parameter will be updated using the following classical formula
[12, eqn(18.32)]. Given a (fixed) parameter $\tau\in(0,1)$, we set at every
iteration

$\pi_{k}=\left\\{\begin{array}[]{rl}\pi_{k-1}&\mbox{if}\ \
{\pi_{k-1}}\geq\frac{1}{1-\tau}\|\hat{\lambda}_{k}\|_{\infty}\\\
\frac{2}{1-\tau}\|\hat{\lambda}_{k}\|_{\infty}&\mbox{otherwise.}\end{array}\right.$
(2.3)

The factor 2 in the second line of (2.3) is introduced so that when $\pi_{k}$
is increased, it is increased substantially. We will see that this rule
ensures that $\pi_{k}$ is eventually fixed. (In general, SQP methods do not
set $H_{k}=\beta_{k}I$. In that case, using the least squares multiplier
estimate in (2.3) will not lead to a convergent method.)

The algorithm for solving problem (1.1), when only noisy evaluations of the
functions $\tilde{f},\tilde{c},\tilde{g},\tilde{J}$ are available, is as
follows.

Algorithm 1 Noise Tolerant SQP Algorithm

Input: Initial iterate $x_{0}$, initial merit parameter $\pi_{-1}>0$, bounds
$\epsilon_{f},\epsilon_{c}$ on the noise (3.1), and constants
$\tau,\nu\in(0,1)$.
Set $k\leftarrow 0$
Repeat until a termination test is satisfied:

1:Compute $\beta_{k}>0$ and set $H_{k}=\beta_{k}I$ in (1.4)

2:Compute $d_{k}$ by solving (1.4)-(1.5)

3:Compute $\hat{\lambda}_{k}$ via (2.2)

4:Update penalty parameter $\pi_{k}$ by (2.3)

5:Compute $\tilde{\ell}(x_{k};d_{k})$ as in (2.1)

6:Set $\epsilon_{R}=2(\epsilon_{f}+\pi_{k}\epsilon_{c})$

7:Choose steplength $\alpha_{k}>0$ such that
$\tilde{\phi}(x_{k}+\alpha_{k}d_{k})\leq\tilde{\phi}(x_{k})+\nu\alpha_{k}\tilde{\ell}(x_{k};d_{k})+\epsilon_{R},$
(2.4)

8:Compute new iterate: $x_{k+1}=x_{k}+\alpha_{k}d_{k}$

9:Set $k\leftarrow k+1$

The steplength $\alpha_{k}$ is computed in Step 7 using a backtracking line
search. We refer to (2.4) as the _relaxed Armijo condition_. The term
$\epsilon_{R}$ introduces a margin that facilitates the convergence analysis
in the presence of noise, and as discussed in Section 4, is also useful in
practice. Note that the line search cannot fail since (2.4) is satisfied for
sufficiently small $\alpha_{k}$, by definition of $\epsilon_{R}$. In this
paper, we assume that the quadratic subproblem (1.4)-(1.5) has a unique
solution at every iteration—admittedly a strong assumption, but one that helps
us focus on the effect of noise without the complicating effects of
regularization parameters or trust regions. The study of a practical algorithm
that employs those globalization strategies will be the subject of future
work.

##  3 Global Convergence

In this section we show that the iterates generated by Algorithm 1 converge to
a neighborhood of the solution determined by the noise level and certain
characteristics of the problem. We also show that once the iterates reach this
neighborhood they cannot stray away from it (under normal circumstances). We
start by stating the assumptions upon which our analysis is built.

######  Assumptions 3.1.

The function $f$ has a Lipschitz continuous gradient with constant $L_{f}$.
The functions $\nabla c_{i}$ are Lipschitz continuous for $i=1,\ldots,m$ with
the corresponding constants held in the vector $L_{c}$.

We also assume that the error (or noise) in the evaluation of the functions is
bounded.

######  Assumptions 3.2.

There exist positive constants
$\epsilon_{f},\epsilon_{c},\epsilon_{g},\epsilon_{J}$ such that for all
$x\in\mathbb{R}^{n}$,

$\displaystyle|\tilde{f}(x)-f(x)|\leq\epsilon_{f},$
$\displaystyle\qquad\|\tilde{c}(x)-c(x)\|_{1}\leq\epsilon_{c},$ (3.1)
$\displaystyle\|\tilde{g}(x)-g(x)\|\leq\epsilon_{g},$
$\displaystyle\qquad\|\tilde{J}(x)-J(x)\|_{1,2}\leq\epsilon_{J}.$ (3.2)

Here, $\|\cdot\|$ denotes the Euclidean norm and $\|\cdot\|_{1,2}$ denotes the
matrix norm induced by the $\ell_{1}$ norm on $\mathbb{R}^{m}$ and the
Euclidean norm on $\mathbb{R}^{n}$.

As already mentioned, we assume that, for all $k$, the matrices
$\tilde{J}_{k}$ have full rank so that the quadratic problem (1.4)-(1.5) has a
unique solution. To state this precisely, we let $\sigma_{min}(A)$ denote the
smallest singular value of a matrix $A$.

######  Assumptions 3.3.

For all $k$, the scalar $\beta_{k}$ in (1.7) satisfies

$0<b_{l}\leq\beta_{k}\leq b_{u},$ (3.3)

for some constants $b_{l},b_{u}$, and there is a constant $\gamma>0$ such that

$\sigma_{min}(J_{k})\geq\gamma,\quad\mbox{with }\ \gamma>\epsilon_{J},\ \
\forall k.$ (3.4)

Furthermore, the sequences ${\\{f_{k}\\}},\\{\|c_{k}\|\\}$, $\\{\|g_{k}\|\\}$,
$\\{\|J_{k}\|\\}$ generated by the algorithm are bounded.

By the matrix inversion lemma [8] and (3.2), if $J_{k}$ has full rank and
$\gamma>\epsilon_{J}$, then $\tilde{J}_{k}$ is also full rank and

$\|\tilde{J}_{k}^{T}(\tilde{J}_{k}\tilde{J}_{k}^{T})^{-1}\|\leq\frac{1}{\gamma-\epsilon_{J}}\equiv\delta,\
\ \forall k.$ (3.5)

The assumption that the sequences ${\\{f_{k}\\}},\\{\|c_{k}\|\\}$,
$\\{\|g_{k}\|\\}$, $\\{\|J_{k}\|\\}$ generated by the algorithm are bounded is
fairly standard in the literature and is designed to avoid pathological
situations. For example, the merit function $\phi$ may be unbounded below away
from the solution if $\pi$ is not large enough. Although there are strategies
to avoid these situations (see e.g. [12, §18.5], we do not include them in our
algorithm, for simplicity.

Given these three sets of assumptions, we are ready to study the convergence
properties of Algorithm 1. Let us apply the well known descent lemma (see
e.g.[4]) to the true (noiseless) merit function

$\phi(x)=f(x)+\pi\|c(x)\|_{1}.$ (3.6)

We have that for any $(x,d)$

$\displaystyle\ \phi(x+\alpha d)\leq\phi(x)+\alpha
g(x)^{T}d+\pi\big{[}\|c(x)+\alpha
J(x)d\|_{1}-\|c(x)\|_{1}\big{]}+\tfrac{1}{2}\big{(}L_{f}+\pi\|L_{c}\|_{1}\big{)}\alpha^{2}\|d\|^{2}.$
(3.7)

Thus, we can write

$\displaystyle\ \phi(x+\alpha d)-\phi(x)\leq\ell(x;\alpha
d)+\tfrac{1}{2}\big{(}L_{f}+\pi\|L_{c}\|_{1}\big{)}\alpha^{2}\|d\|^{2},$ (3.8)

where

$\ell(x;s)=g(x)^{T}s+\pi\|c(x)+J(x)s\|_{1}-\pi\|c(x)\|_{1}.$ (3.9)

When function and derivatives are exact, it is easy to show that for $\pi$
sufficiently large and $\alpha$ sufficiently small we can guarantee a
reduction in $\phi$; see [12]. We must establish that this is also the case in
the noisy setting—before the iterates approach the region around the solution
where noise dominates. We begin by establishing bounds on the step $d_{k}$.

###  3.1 Preliminary results

The optimality conditions of the quadratic problem (1.4)-(1.5) are given by

$\displaystyle\begin{pmatrix}H_{k}&\tilde{J}_{k}^{T}\\\
\tilde{J}_{k}&0\end{pmatrix}\begin{pmatrix}d_{k}\\\
d_{y}\end{pmatrix}=-\begin{pmatrix}\tilde{g}_{k}+\tilde{J}_{k}^{T}y\\\
\tilde{c}_{k}\end{pmatrix},$ (3.10)

for some Lagrange multiplier $y\in\mathbb{R}^{m}$. The step $d_{k}$ can be
written as the sum of two orthogonal components,

$d_{k}=v_{k}+u_{k},$ (3.11)

where $v_{k}$ is in the range space of $\tilde{J}_{k}^{T}$ and $u_{k}$ is in
the null space of $J_{k}$. A simple computation from (3.10) shows that

$v_{k}=-\tilde{J}_{k}^{T}(\tilde{J}_{k}\tilde{J}_{k}^{T})^{-1}\tilde{c}_{k},\quad\mbox{}\quad
u_{k}=-\frac{1}{\beta_{k}}\tilde{P}_{k}\tilde{g}_{k},$ (3.12)

where

$\tilde{P}_{k}=I-\tilde{J}_{k}^{T}\left(\tilde{J}_{k}\tilde{J}_{k}^{T}\right)^{-1}\tilde{J}_{k}$
(3.13)

is an orthogonal projection matrix onto the tangent space of the constraints.
We now establish bounds on $u_{k},v_{k}$. In what follows, we let
$J^{\dagger}$ denote the Moore-Penrose generalized inverse of a matrix $J$,
and define $P_{k}=I-J_{k}^{T}\left(J_{k}J_{k}^{T}\right)^{-1}J_{k}$. Since
$\tilde{P}_{k}$ and $P_{k}$ are orthogonal projections, we have that
$\|\tilde{P}_{k}\|=\|P_{k}\|=1$.

######  Lemma 3.4.

Under Assumptions 3.1 and 3.2 we have both

$\displaystyle\|v_{k}\|_{1}$
$\displaystyle\leq\delta\|\tilde{c}_{k}\|_{1}\leq{\delta}(\|c_{k}\|_{1}+\epsilon_{c})$
(3.14) $\displaystyle\|u_{k}\|$
$\displaystyle\leq\frac{1}{\beta_{k}}\big{(}\|P_{k}g_{k}\|+\|g_{k}\|\eta\epsilon_{J}+\epsilon_{g}\big{)},$
(3.15)

where $\delta$ is defined in (3.5) and

$\eta=1/\gamma.$ (3.16)

Therefore,

$\displaystyle\|d_{k}\|\leq\delta(\|c_{k}\|_{1}+\epsilon_{c})+\frac{1}{\beta_{k}}\big{(}\|P_{k}g_{k}\|+\|g_{k}\|\eta\epsilon_{J}+\epsilon_{g}\big{)}.$
(3.17)

###### Proof.

The bounds (3.14) follow directly from (3.12), (3.5), and (3.1). By (3.2), we
can bound the norm of the tangential component as follows

$\displaystyle\|u_{k}\|$
$\displaystyle=\frac{1}{\beta_{k}}\|\tilde{P}_{k}\tilde{g}_{k}\|$
$\displaystyle\leq\frac{1}{\beta_{k}}\left(\|P_{k}g_{k}\|+\|(\tilde{P}_{k}-P_{k})g_{k}\|+\|\tilde{P}_{k}\|\|g_{k}-\tilde{g}_{k}\|\right)$
$\displaystyle\leq\frac{1}{\beta_{k}}\left(\|P_{k}g_{k}\|+\|(\tilde{P}_{k}-P_{k})\|\|g_{k}\|+\epsilon_{g}\right).$
(3.18)

Moreover, by the bounds on perturbed projection matrices [16, Theorems 2.3 and
2.4] we have that

$\|\tilde{P}_{k}-P_{k}\|\leq\frac{\epsilon_{J}}{\gamma}\equiv\eta\epsilon_{J}.$
(3.19)

This yields (3.15). ∎

###  3.2 Penalty Parameter and Model Decrease

We note from (3.8) that in order to obtain a decrease in the true merit
function $\phi$, we must ensure that $\ell(x_{k};\alpha_{k}d_{k})$ is
negative. We will see that this can be achieved for $\alpha_{k}=1$ by choosing
a sufficiently large penalty parameter $\pi$, and provided noise does not
dominate.

######  Lemma 3.5.

If at every iteration $k$ the penalty parameter satisfies

$\pi_{k}\geq\frac{1}{1-\tau}\|(\tilde{J}_{k}\tilde{J}_{k}^{T})^{-1}\tilde{J}_{k}\tilde{g}_{k}\|_{\infty},\quad\tau\in(0,1),$
(3.20)

then

$\displaystyle\ell(x_{k};d_{k})\leq$
$\displaystyle-\frac{1}{\beta_{k}}g_{k}^{T}P_{k}g_{k}+\frac{1}{\beta_{k}}(\|g_{k}\|^{2}\eta\epsilon_{J}+\epsilon_{g}\|g_{k}\|)-\tau\pi_{k}\|c_{k}\|_{1}+{\epsilon_{g}}{\delta}(\|c_{k}\|_{1}+\epsilon_{c})$
(3.21)
$\displaystyle+\pi_{k}\left((2-\tau)\epsilon_{c}+\epsilon_{J}\left(\delta(\|c_{k}\|_{1}+\epsilon_{c})+\frac{1}{\beta_{k}}(\|P_{k}g_{k}\|+\|g_{k}\|\eta\epsilon_{J}+\epsilon_{g})\right)\right).$

###### Proof.

Since
$d_{k}=-\frac{1}{\beta_{k}}\tilde{P}_{k}\tilde{g}_{k}-\tilde{J}_{k}^{T}(\tilde{J}_{k}\tilde{J}_{k}^{T})^{-1}\tilde{c}_{k},$
we have from (3.9), (1.5), (3.5), (3.1), (3.2), and the definition of the
$\|\cdot\|_{1,2}$ norm in (3.2), that

$\displaystyle\ell(x_{k};d_{k})=$ $\displaystyle\
g_{k}^{T}d_{k}+\pi_{k}\|c_{k}+J_{k}d_{k}\|_{1}-\pi_{k}\|c_{k}\|_{1}$ (3.22)
$\displaystyle\leq$
$\displaystyle-\frac{1}{\beta_{k}}g_{k}^{T}\tilde{P}_{k}\tilde{g}_{k}-g_{k}^{T}\tilde{J}_{k}^{T}(\tilde{J}_{k}\tilde{J}_{k}^{T})^{-1}\tilde{c}_{k}+\pi_{k}\|c_{k}+J_{k}d_{k}\|_{1}-\pi_{k}\|c_{k}\|_{1}$
$\displaystyle\leq$
$\displaystyle-\frac{1}{\beta_{k}}g_{k}^{T}\tilde{P}_{k}\tilde{g}_{k}-g_{k}^{T}\tilde{J}_{k}^{T}(\tilde{J}_{k}\tilde{J}_{k}^{T})^{-1}\tilde{c}_{k}+\pi_{k}\|(c_{k}-\tilde{c}_{k})+(J_{k}-\tilde{J}_{k})d_{k})\|_{1}-\pi_{k}\|c_{k}\|_{1}$
$\displaystyle\leq$
$\displaystyle-\frac{1}{\beta_{k}}g_{k}^{T}\tilde{P}_{k}\tilde{g}_{k}-\tilde{g}_{k}^{T}\tilde{J}_{k}^{T}(\tilde{J}_{k}\tilde{J}_{k}^{T})^{-1}\tilde{c}_{k}+{\epsilon_{g}}{\delta}\|\tilde{c}_{k}\|_{1}+\pi_{k}(\epsilon_{c}+\epsilon_{J}\|d_{k}\|)-\pi_{k}\|c_{k}\|_{1}$
$\displaystyle\leq$
$\displaystyle-\frac{1}{\beta_{k}}g_{k}^{T}\tilde{P}_{k}\tilde{g}_{k}-\tilde{g}_{k}^{T}\tilde{J}_{k}^{T}(\tilde{J}_{k}\tilde{J}_{k}^{T})^{-1}\tilde{c}_{k}+{\epsilon_{g}}{\delta}(\|c_{k}\|_{1}+\epsilon_{c})$
$\displaystyle+\pi_{k}\left[\epsilon_{c}+\epsilon_{J}\left(\delta(\|c_{k}\|_{1}+\epsilon_{c})+\frac{1}{\beta_{k}}(\|P_{k}g_{k}\|+\|g_{k}\|\eta\epsilon_{J}+\epsilon_{g})\right)\right]-\pi_{k}\|c_{k}\|_{1},$

the last line following by (3.17). Next, since $\|\tilde{P}_{k}\|=1$ and
recalling (3.19), we obtain

$\displaystyle-g_{k}^{T}\tilde{P}_{k}\tilde{g}_{k}$ $\displaystyle\leq-
g_{k}^{T}P_{k}g_{k}+\|g_{k}\|\|P_{k}g_{k}-\tilde{P}_{k}\tilde{g}_{k}\|$
$\displaystyle\leq-
g_{k}^{T}P_{k}g_{k}+\|g_{k}\|\|P_{k}g_{k}-\tilde{P}_{k}g_{k}+\tilde{P}_{k}g_{k}-\tilde{P}_{k}\tilde{g}_{k}\|$
$\displaystyle\leq-
g_{k}^{T}P_{k}g_{k}+\|g_{k}\|^{2}\|P_{k}-\tilde{P}_{k}\|+\|g_{k}\|\|g_{k}-\tilde{g}_{k}\|$
$\displaystyle\leq-
g_{k}^{T}P_{k}g_{k}+\|g_{k}\|^{2}\eta\epsilon_{J}+\|g_{k}\|\epsilon_{g}.$

Therefore,

$\displaystyle\ell(x_{k};d_{k})$
$\displaystyle\leq-\frac{1}{\beta_{k}}g_{k}^{T}P_{k}g_{k}+\frac{1}{\beta_{k}}(\|g_{k}\|^{2}\eta\epsilon_{J}+\epsilon_{g}\|g_{k}\|)-\tilde{g}_{k}^{T}\tilde{J}_{k}^{T}(\tilde{J}_{k}\tilde{J}_{k}^{T})^{-1}\tilde{c}_{k}+{\epsilon_{g}}{\delta}(\|c_{k}\|_{1}+\epsilon_{c})$
$\displaystyle+\pi_{k}\left[\epsilon_{c}+\epsilon_{J}\left(\delta(\|c_{k}\|_{1}+\epsilon_{c})+\frac{1}{\beta_{k}}(\|P_{k}g_{k}\|+\|g_{k}\|\eta\epsilon_{J}+\epsilon_{g})\right)\right]-\pi_{k}\|c_{k}\|_{1}.$

Now suppose that we choose the parameter $\pi_{k}$ so that (3.20) holds. Then

$-\tilde{g}_{k}^{T}\tilde{J}_{k}^{T}(\tilde{J}_{k}\tilde{J}_{k}^{T})^{-1}\tilde{c}_{k}\leq\|\tilde{g}_{k}^{T}\tilde{J}_{k}^{T}(\tilde{J}_{k}\tilde{J}_{k}^{T})^{-1}\|_{\infty}\|\tilde{c}_{k}\|_{1}\leq(1-\tau)\pi_{k}(\|c_{k}\|_{1}+\epsilon_{c}),$

and it follows that

$\displaystyle\ell(x_{k};d_{k})\leq$
$\displaystyle-\frac{1}{\beta_{k}}g_{k}^{T}P_{k}g_{k}+\frac{1}{\beta_{k}}(\|g_{k}\|^{2}\eta\epsilon_{J}+\epsilon_{g}\|g_{k}\|)-\tau\pi_{k}\|c_{k}\|_{1}+{\epsilon_{g}}{\delta}(\|c_{k}\|_{1}+\epsilon_{c})$
$\displaystyle+\pi_{k}\left[(2-\tau)\epsilon_{c}+\epsilon_{J}\left(\delta(\|c_{k}\|_{1}+\epsilon_{c})+\frac{1}{\beta_{k}}(\|P_{k}g_{k}\|+\|g_{k}\|\eta\epsilon_{J}+\epsilon_{g})\right)\right].$

∎

Lemma 3.5 implies that for any $x_{k}$ such that the right hand side of (3.21)
is negative, we have $\ell(x_{k};d_{k})<0.$ We now provide conditions under
which the decrease in $\ell$ is proportional to the optimality conditions of
the nonlinear problem (1.1). Specifically, since
$g_{k}^{T}P_{k}g_{k}=\|P_{k}g_{k}\|^{2}$ is the norm squared of the projected
gradient, a combination of $g_{k}^{T}P_{k}g_{k}$ and $\|c_{k}\|_{1}$ can be
regarded as a measure of stationarity of the constrained optimization problem.
The following result assumes that the optimality measure is not small compared
to the errors (or noise).

######  Corollary 3.6.

Choose any $\theta_{1}\in[0,1)$. For any $x_{k}$ sufficiently far from the
solution such that

$\displaystyle(1-\theta_{1})\Big{(}\frac{1}{\beta_{k}}g_{k}^{T}P_{k}g_{k}+$
$\displaystyle\tau\pi_{k}\|c_{k}\|_{1}\Big{)}\geq{E}(x_{k},\beta_{k},\pi_{k}),$
(3.23)

where

$\displaystyle{E}(x,\beta,\pi)=$
$\displaystyle\frac{1}{\beta}(\|g(x)\|^{2}\eta\epsilon_{J}+\epsilon_{g}\|g(x)\|)+{\epsilon_{g}}{\delta}{(\|c(x)\|_{1}+\epsilon_{c})}$
$\displaystyle+$
$\displaystyle\pi\left[{(2-\tau)\epsilon_{c}}+\epsilon_{J}\left(\delta(\|c(x)\|_{1}+\epsilon_{c})+\frac{1}{\beta}(\|P(x)g(x)\|+\|g(x)\|\eta\epsilon_{J}+\epsilon_{g})\right)\right],$
(3.24)

we have

$\ell(x_{k};d_{k})\leq-\theta_{1}\left(\frac{1}{\beta_{k}}g_{k}^{T}P_{k}g_{k}+\tau\pi_{k}\|c_{k}\|_{1}\right).$
(3.25)

###### Proof.

For any $\theta_{1}\in[0,1)$, we can rewrite (3.21) as

$\displaystyle\ell(x_{k};d_{k})\leq$
$\displaystyle-\theta_{1}(\frac{1}{\beta_{k}}g_{k}^{T}P_{k}g_{k}+\tau\pi_{k}\|c_{k}\|_{1})-(1-\theta_{1})(\frac{1}{\beta_{k}}g_{k}^{T}P_{k}g_{k}+\tau\pi_{k}\|c_{k}\|_{1})$
$\displaystyle+\frac{1}{\beta_{k}}(\|g_{k}\|^{2}\eta\epsilon_{J}+\epsilon_{g}\|g_{k}\|)+{\epsilon_{g}}{\delta}\|{(\|c_{k}\|_{1}+\epsilon_{c})}$
$\displaystyle+\pi_{k}\left({(2-\tau)\epsilon_{c}}+\epsilon_{J}\left(\delta(\|c_{k}\|_{1}+\epsilon_{c})+\frac{1}{\beta_{k}}(\|P_{k}g_{k}\|+\|g_{k}\|\eta\epsilon_{J}+\epsilon_{g})\right)\right),$

from which (3.25) follows by condition (3.23). ∎

In order to make this result, and similar results to be proved later, more
understandable and more convenient to use, we recall that
$g(x)^{T}P(x)g(x)=\|P(x)g(x)\|^{2}$, and define the function

$\psi_{\pi}(x)=\frac{1}{b_{u}}\|P(x)g(x)\|^{2}+\pi\tau\|c(x)\|_{1},$ (3.26)

where $b_{u}$ is given in (3.3). Clearly, $\psi_{\pi}$ may be viewed as a
measure of non-stationarity since $\psi_{\pi}(x^{\ast})=0$ when $x^{\ast}$ is
a stationary point of the problem (1.1). Given this notation we can restate a
slightly weaker version of Corollary 3.6.

######  Corollary 3.7.

Choose any $\theta_{1}\in[0,1)$. For any $x_{k}$ sufficiently far from the
solution such that

$\displaystyle\psi_{\pi_{k}}(x_{k})\geq
E(x_{k},\beta_{k},\pi_{k})/(1-\theta_{1}),$ (3.27)

we have

$\ell(x_{k};d_{k})\leq-\theta_{1}\left(\frac{1}{\beta_{k}}g_{k}^{T}P_{k}g_{k}+\tau\pi_{k}\|c_{k}\|_{1}\right)\leq-\theta_{1}\psi_{\pi_{k}}(x_{k}).$
(3.28)

###### Proof.

The result follows from the fact that

$\psi_{\pi_{k}}(x_{k})\leq\Big{(}\frac{1}{\beta_{k}}g_{k}^{T}P_{k}g_{k}+\tau\pi_{k}\|c_{k}\|_{1}\Big{)}.$

∎

###  3.3 Line search

Since $\pi_{k}$ is defined by (2.3) and (2.2), and by Assumptions 3.3, we have
that $\\{\|\hat{\lambda}_{k}\|\\}$ is bounded. Moreover, since $\\{\pi_{k}\\}$
is monotone and since $\pi_{k}-\pi_{k-1}$ is either zero or greater than
$\pi_{k-1}$, there exists values $k_{0}$ and $\bar{\pi}$ such that:

$\pi_{k}=\bar{\pi},\quad\forall k\geq k_{0},$ (3.29)

and (3.20) is satisfied. The rest of the analysis assumes that the penalty
parameter has attained that fixed value $\bar{\pi}$. Thus, for the rest of the
section

$\tilde{\phi}(x_{k})\equiv\tilde{f}(x_{k})+\bar{\pi}\|\tilde{c}(x_{k})\|_{1},\qquad\phi(x_{k})\equiv
f(x_{k})+\bar{\pi}\|c(x_{k})\|_{1},\ \forall k\geq k_{0}.$ (3.30)

Algorithm 1 sets $x_{k+1}=x_{k}+\alpha_{k}d_{k}$, where $\alpha_{k}$ is chosen
by repeated halving until the relaxed Armijo condition is satisfied:

$\tilde{\phi}(x_{k}+\alpha_{k}d_{k})\leq\tilde{\phi}(x_{k})+{\nu}\alpha_{k}\tilde{\ell}(x_{k};d_{k})+\epsilon_{R},$

for some constants ${\nu}\in(0,1)$ and
$\epsilon_{R}\,{\geq}\,2(\epsilon_{f}+\bar{\pi}\epsilon_{c}),$ where
$\tilde{\ell}(x_{k};d_{k})$ is defined in (2.1). In other words, we require
that the decrease in the noisy merit function be a fraction of the decrease of
the noisy first-order model $\tilde{\ell}$, plus a relaxation term.

To ensure that the line search yields significant progress toward a solution,
we need to show that $\alpha_{k}$ is bounded away from zero and that
$\tilde{\ell}(x_{k};d_{k})$ is sufficiently negative. To do so, we recall that
we have established in (3.28) that the noiseless first-order model
$\ell(x_{k};d)_{k}$ is sufficiently negative when condition (3.27) is
satisfied. To relate $\ell(x_{k};d_{k})$ to $\tilde{\ell}(x_{k};d_{k})$, we
recall (3.9) and (2.1), and measure the difference between these two
quantities. By (3.17)

$\displaystyle|\tilde{\ell}(x_{k};d_{k})-\ell(x_{k};d_{k})|\leq$
$\displaystyle\
\epsilon_{g}\|d_{k}\|+2\bar{\pi}\epsilon_{c}+\bar{\pi}\epsilon_{J}\|d_{k}\|$
$\displaystyle\leq$ $\displaystyle\
(\epsilon_{g}+\bar{\pi}\epsilon_{J})\left(\delta(\|c_{k}\|_{1}+\epsilon_{c})+\frac{1}{\beta_{k}}(\|P_{k}g_{k}\|+\|g_{k}\|\eta\epsilon_{J}+\epsilon_{g})\right)+2\bar{\pi}\epsilon_{c}$
(3.31) $\displaystyle\leq$ $\displaystyle\
{(\epsilon_{g}+\bar{\pi}\epsilon_{J})\left(\delta(C_{c}+\epsilon_{c})+\frac{1}{b_{l}}(C_{g}+C_{g}\eta\epsilon_{J}+\epsilon_{g})\right)+2\bar{\pi}\epsilon_{c}}$
$\displaystyle\equiv$ $\displaystyle\ {\epsilon_{\ell},}$ (3.32)

where $C_{g},C_{c}$ are constants such that

$\|g(x_{k})\|\leq C_{g},\qquad\|c(x_{k})\|_{1}\leq C_{c}\quad\forall k>k_{0}.$
(3.33)

We know that these constants exist because of Assumption 3.3. We now describe
conditions under which one can characterize the size of the steplength
$\alpha_{k}$. Let

$L=L_{f}+\bar{\pi}\|L_{c}\|_{1},$ (3.34)

where $L_{f},L_{c}$ are defined in Assumptions 3.1.

######  Theorem 3.8.

Let $\theta_{1}$ be defined as in Corollary 3.6, choose constants
$\theta_{2}<\theta_{1}$, $\nu\in(0,1)$ and

$\epsilon_{R}\,{\geq}\,2(\epsilon_{f}+\bar{\pi}\epsilon_{c})\equiv
2\epsilon_{\phi}.$ (3.35)

Then, for all iterates $x_{k}$ with $k\geq k_{0}$ that satisfy both (3.27) and

$(1-\nu)(\theta_{1}-\theta_{2})\left(\frac{1}{\beta_{k}}\|P_{k}g_{k}\|^{2}+{\bar{\pi}}\tau\|c_{k}\|_{1}\right)>2\nu\epsilon_{\ell},$
(3.36)

if the steplength satisfies

$\alpha_{k}<\frac{(1-\nu)\theta_{2}\left({\frac{1}{\beta_{k}}}\|P_{k}g_{k}\|^{2}+\bar{\pi}\tau\|c_{k}\|_{1}\right)}{\frac{L}{2}[\delta^{2}(\|c_{k}\|_{1}+\epsilon_{c})^{2}+\frac{1}{\beta_{k}^{2}}(\|P_{k}g_{k}\|+\|g_{k}\|\eta\epsilon_{J}+\epsilon_{g})^{2}]}\equiv\hat{\alpha}_{k},$
(3.37)

then

$\tilde{\phi}(x_{k}+\alpha_{k}d_{k})\leq\tilde{\phi}(x_{k})+\nu\alpha_{k}\tilde{\ell}(x_{k};d_{k})+\epsilon_{R}.$
(3.38)

###### Proof.

By the definition (3.35) of $\epsilon_{\phi}$, (3.8), (3.34), the convexity of
$\ell(x_{k};\cdot)$, (3.32), (3.28), the fact that $P_{k}^{2}=P_{k}$, and
(3.17), we get

$\displaystyle\tilde{\phi}(x_{k}+\alpha d_{k})-\tilde{\phi}(x_{k})\leq$
$\displaystyle\,\phi(x_{k}+\alpha d_{k})-\phi(x_{k})+2\epsilon_{\phi}$
$\displaystyle\leq$ $\displaystyle\,\ell(x_{k};\alpha
d_{k})+2\epsilon_{\phi}+\tfrac{L}{2}\alpha^{2}\|d_{k}\|^{2}$
$\displaystyle\leq$
$\displaystyle\,\alpha\ell(x_{k};d_{k})+2\epsilon_{\phi}+\tfrac{L}{2}\alpha^{2}\|d_{k}\|^{2}$
$\displaystyle=$
$\displaystyle\,{\nu\alpha\ell(x_{k};d_{k})+2\epsilon_{\phi}+(1-\nu)\alpha\ell(x_{k};d_{k})+\tfrac{L}{2}\alpha^{2}\|d_{k}\|^{2}}$
$\displaystyle\leq$
$\displaystyle\,{\nu\alpha\tilde{\ell}(x_{k};d_{k})+\nu\alpha\epsilon_{\ell}+2\epsilon_{\phi}+(1-\nu)\alpha\ell(x_{k};d_{k})+\tfrac{L}{2}\alpha^{2}\|d_{k}\|^{2}}$
$\displaystyle\leq$
$\displaystyle\,\nu\alpha\tilde{\ell}(x_{k};d_{k})+2\epsilon_{\phi}+2\nu\alpha\epsilon_{\ell}-(1-\nu)\theta_{1}\alpha\Big{(}\frac{1}{\beta_{k}}g_{k}^{T}P_{k}g_{k}+{\bar{\pi}}\tau\|c_{k}\|_{1}\Big{)}$
$\displaystyle+\tfrac{L}{2}\alpha^{2}\|d_{k}\|^{2}$ $\displaystyle\leq$
$\displaystyle\,\nu\alpha\tilde{\ell}(x_{k};d_{k})+2\epsilon_{\phi}+2\nu\alpha\epsilon_{\ell}-(1-\nu)\theta_{1}\alpha\Big{(}\frac{1}{\beta_{k}}\|P_{k}g_{k}\|^{2}+{\bar{\pi}}\tau\|c_{k}\|_{1}\Big{)}$
$\displaystyle+\tfrac{L}{2}\alpha^{2}\Big{[}\delta^{2}(\|c_{k}\|_{1}+\epsilon_{c})^{2}+\frac{1}{\beta_{k}^{2}}\big{(}\|P_{k}g_{k}\|+\|g_{k}\|\eta\epsilon_{J}+\epsilon_{g}\big{)}^{2}\Big{]},$

the last line following from the orthogonality of the components (3.11) of
$d_{k}$.

Now we choose a constant $\theta_{2}<\theta_{1}$, and consider iterates
$x_{k}$ such that (3.36) holds. For such iterates we have,

$\displaystyle\tilde{\phi}(x_{k}+\alpha d_{k}){-\tilde{\phi}(x_{k})}\leq$
$\displaystyle\nu\alpha\tilde{\ell}(x_{k};d_{k})+2\epsilon_{\phi}-(1-\nu)\theta_{2}\alpha\Big{(}{\frac{1}{\beta_{k}}}\|P_{k}g_{k}\|^{2}+\bar{\pi}\tau\|c_{k}\|_{1}\Big{)}$
$\displaystyle+\tfrac{L}{2}\alpha^{2}\left[\delta^{2}(\|c_{k}\|_{1}+\epsilon_{c})^{2}+\frac{1}{\beta_{k}^{2}}\Big{(}\|P_{k}g_{k}\|+\|g_{k}\|\eta\epsilon_{J}+\epsilon_{g}\Big{)}^{2}\right].$

Then, for any steplength satisfying (3.37) where $x_{k}$ satisfies the (3.23)
and (3.36), we have

$\tilde{\phi}(x_{k}+\alpha
d_{k})-\tilde{\phi}(x_{k})\leq\nu\alpha\tilde{l}(x_{k};d_{k})+2\epsilon_{\phi},$

and thus (3.38) holds since $\epsilon_{R}\,{\geq 2}\epsilon_{\phi}$. ∎

Note that condition (3.36) is implied by the slightly weaker inequality

$(1-\nu)(\theta_{1}-\theta_{2})\psi_{\bar{\pi}}(x_{k})>2\nu\epsilon_{\ell}.$
(3.39)

Since the numerator in (3.37) is bounded away from zero by (3.36), and the
denominator is bounded above given the assumed global upper bounds on
$c_{k},g_{k},$ and lower bound on $\beta_{k}$ stated in Assumptions 3.3, it
follows that there is a constant $\bar{\alpha}$ such that
$\hat{\alpha}_{k}>2\bar{\alpha}$ for all $k\geq k_{0}$. The algorithm employs
a backtracking line search that halves each trial step, hence we can conclude
that

$\bar{\alpha}\leq\alpha_{k},\qquad\mbox{for}\ k\geq k_{0}.$ (3.40)

This will allow us to show that, when the conditions in Theorem 3.8 are
satisfied, the algorithm will make non-negligible progress.

###  3.4 The Main Convergence Result

Now we show that Algorithm 1 will eventually generate iterates close to a
stationary point of the problem, as measured by the function
$\psi_{\bar{\pi}}(x)$ defined in (3.26). To do so, we note that condition
(3.27) implies that the linear model decrease $\ell$ is sufficiently negative
in the sense of (3.28), and we have established a bound in (3.32) for the
distance between $\ell$ and $\tilde{\ell}$. Furthermore, we have shown that
condition (3.39) ensures that the relaxed Armijo condition (3.38) is satisfied
for steplengths $\alpha_{k}$ that are bounded away from zero. Those two
conditions—(3.27), (3.39)—are necessary to ensure that the algorithm makes
significant progress, but they are not sufficient. To control the effect of
noise in testing (2.4) as well as the effect of the relaxation factor, we
impose one additional condition,

$\psi_{\bar{\pi}}(x_{k})\geq\frac{2\epsilon_{R}+4\epsilon_{\phi}}{\nu\bar{\alpha}\theta_{2}},$
(3.41)

to help define the region where Algorithm 1 progresses toward stationarity.

One more refinement is needed. The definition of the term
$E(x_{k},\beta_{k},\pi_{k})$ defined in (3.24) involves $c(x_{k})$ and
$g(x_{k})$, which makes the region defined by (3.27) difficult to interpret.
Therefore, we compute an upper bound for $E$. If we define

$\displaystyle{\cal E}=$
$\displaystyle\frac{1}{b_{l}}(C_{g}^{2}\eta\epsilon_{J}+\epsilon_{g}C_{g})+{\epsilon_{g}}{\delta}{(C_{c}+\epsilon_{c})}$
$\displaystyle+$
$\displaystyle\bar{\pi}\left[{(2-\tau)\epsilon_{c}}+\epsilon_{J}\left(\delta(C_{c}+\epsilon_{c})+\frac{1}{b_{l}}(C_{g}+C_{g}\eta\epsilon_{J}+\epsilon_{g})\right)\right],$
(3.42)

where $C_{g},C_{c}$ are given in (3.33), then we have that
${E}(x_{k},\beta_{k},\bar{\pi})\leq{\cal E}$ for all $k\geq k_{0}$. We can
thus state a condition that implies (3.27):

$\psi_{{\bar{\pi}}}(x_{k})\geq\frac{{\cal E}}{(1-\theta_{1})},\qquad\forall
k\geq k_{0}.$ (3.43)

In summary, the analysis presented above holds if conditions (3.43), (3.39)
are satisfied and we also impose condition (3.41). This allows us to
characterize a region, which we denote by $C_{1}$, where errors dominate and
improvement in the merit function $\phi$ cannot be guaranteed. In other words,
$C_{1}$ is the region where at least one of the three conditions—(3.43),
(3.39), (3.41)—is not satisfied.

######  Definition 3.9.

The critical region $C_{1}$ is defined as the set of $x\in\mathbb{R}^{n}$
satisfying

$\displaystyle\psi_{\bar{\pi}}(x)\leq$ $\displaystyle\max{\left\\{\frac{{\cal
E}}{(1-\theta_{1})},\frac{2{\nu}\epsilon_{\ell}}{({1-\nu})(\theta_{1}-\theta_{2})},\frac{2\epsilon_{R}+4\epsilon_{\phi}}{\nu\bar{\alpha}\theta_{2}}\right\\},}$
(3.44)

where $\cal E$ and $\epsilon_{\ell}$ are defined by (3.42) and (3.32),
respectively, and $\theta_{1},\theta_{2}$ are constants such that
$0<\theta_{2}<\theta_{1}<1$.

We also define the following set.

######  Definition 3.10.

Let $w=\sup\\{\phi(x):x\in C_{1}\\}$, and define the level set

$C_{2}=\\{x:\phi(x)\leq w+2\epsilon_{\phi}+\epsilon_{R}\\}.$

Note that by construction $C_{1}\subseteq C_{2}$. We are now ready to state
the main convergence result.

######  Theorem 3.11.

Suppose that Algorithm 1 generates a sequence $\\{x_{k}\\}$ from $x_{0}$
satisfying Assumptions 3.1-3.3. There is an iteration $k_{1}$ at which
$\\{x_{k}\\}$ enters the critical region $C_{1}$, and for all $k>k_{1}$ the
iterates remain in the critical level set $C_{2}$. The iterates may leave
$C_{1}$, but there must be infinitely many iterates in $C_{1}$.

###### Proof.

Recall that the index $k_{0}$ is defined in (3.29). If $k\not\in C_{1}$ and
$k\geq k_{0}$, then the assumptions of Theorem 3.8 are satisfied and (3.38)
holds. Therefore, by (3.40), (3.32), (3.25), (3.28), (3.26), (3.36)

$\displaystyle\phi(x_{k}+\alpha_{k}d_{k})-\phi(x_{k})\leq$
$\displaystyle\,{\tilde{\phi}(x_{k}+\alpha_{k}d_{k})-\tilde{\phi}(x_{k})+2\epsilon_{\phi}}$
(3.45) $\displaystyle\leq$
$\displaystyle\,\nu\bar{\alpha}\tilde{\ell}_{\phi}(x_{k};d_{k})+2\epsilon_{\phi}+\epsilon_{R}$
(3.46) $\displaystyle\leq$
$\displaystyle\,\nu\bar{\alpha}\ell_{\phi}(x_{k};d_{k})+\nu\bar{\alpha}\epsilon_{\ell}+2\epsilon_{\phi}+\epsilon_{R}$
$\displaystyle\leq$
$\displaystyle-\nu\bar{\alpha}\theta_{1}\big{(}\frac{1}{\beta_{k}}g_{k}^{T}P_{k}g_{k}+\tau\bar{\pi}\|c_{k}\|_{1}\big{)}+\nu\bar{\alpha}\epsilon_{\ell}+2\epsilon_{\phi}+\epsilon_{R}$
$\displaystyle\leq$
$\displaystyle-\nu\bar{\alpha}\theta_{1}\psi_{\bar{\pi}}(x_{k})+\nu\bar{\alpha}\epsilon_{\ell}+2\epsilon_{\phi}+\epsilon_{R}$
$\displaystyle=$
$\displaystyle-[\nu\bar{\alpha}\theta_{2}+\hat{\alpha}\nu(\theta_{1}-\theta_{2})]\psi_{\bar{\pi}}(x_{k})+\nu\bar{\alpha}\epsilon_{\ell}+2\epsilon_{\phi}+\epsilon_{R}$
$\displaystyle\leq$
$\displaystyle-\nu\bar{\alpha}\theta_{2}\psi_{\bar{\pi}}(x_{k})+2\epsilon_{\phi}+\epsilon_{R}.$
(3.47)

Combining this bound with (3.41), we have that if $x_{k}\notin C_{1}$ then

$\phi(x_{k+1})-\phi(x_{k})\leq-\frac{\nu\bar{\alpha}\theta_{2}}{2}\psi_{\bar{\pi}}(x_{k}).$
(3.48)

Since the sequence $\\{\phi(x_{k})\\}$ is bounded below by Assumptions 3.3,
$\psi_{\bar{\pi}}(x_{k})$ converges to zero and thus it follows that Algorithm
1 eventually generates an iterate in $C_{1}$.

Now if $x_{k}\in C_{1}$, then by Step 6 in Algorithm 1,
$\phi(x_{k+1})\leq\phi(x_{k})+2\epsilon_{\phi}+\epsilon_{R}\leq
w+2\epsilon_{\phi}+\epsilon_{R}$, so that $x_{k+1}\in C_{2}.$

On the other hand, if $x_{k}\in C_{2}\setminus C_{1}$, then by (3.48)

$\phi(x_{k+1})-\phi(x_{k})\,\leq\,0,$

which implies $x_{k+1}\in C_{2}$. Thus the rest of the sequence lies in
$C_{2}$, with infinitely many iterates in $C_{1}$. ∎

We should note that since we are not assuming that the objective function is
strongly convex or satisfies a quadratic growth condition, it is possible that
the supremum in Definition 3.10 is $w=\infty$. This is, however, an unlikely
scenario.

###  3.5 Discussion

Let us take a closer look the main result of this paper, Theorem 3.11, since
the critical region $C_{1}$ defined in (3.44) is complex.

By the definitions (3.42) and (3.32), we have that $\cal E$ and
$\epsilon_{\ell}$ are both of order
$O(\epsilon_{c},\epsilon_{g},\epsilon_{J}$), and so is the right hand side in
(3.44). This is as desired. The constants in these orders of magnitude matter,
so we must characterize them.

First note that the critical region $C_{1}$, the set $C_{2}$ and $\bar{\pi}$
depend on the starting point $x_{0}$. It is then possible that $\bar{\pi}$
could be very large in some cases, although in practice this does not seem to
be a major concern. The constants $C_{g},C_{c}$, which also enter in the
definition of $\cal E$ and $\epsilon_{\ell}$ could be quite large. One can,
however, give a tighter definition of $C_{1}$ by not introducing these
constants. In this case, we would define $\epsilon_{\ell}$ by (3.31) and
employ (3.27), rather than (3.43). This makes the main theorem more precise,
albeit more difficult to interpret.

Returning to the constants in (3.42) and (3.32), we have that

$\epsilon_{\ell},{\cal
E}\sim\left[\delta,\,\frac{1}{b_{l}},\,\frac{\eta}{b_{l}}\right],$

and from (3.4), (3.5), (3.16) we observe that

$\sigma_{min}(\tilde{J}_{k})\geq\gamma,\qquad\delta=\frac{1}{\gamma-\epsilon_{J}}\geq\frac{1}{\sigma_{min}(\tilde{J}_{k})-\epsilon_{J}},\qquad\mbox{and}\
\eta=\frac{1}{\gamma}=\frac{1}{\sigma_{min}(\tilde{J}_{k})}.$

The effect of a near rank-deficient Jacobian and Hessian approximations
$\beta_{k}I$ are now apparent.

It is interesting to compare $C_{1}$ with the region obtained by Berahas et
al. [1] for unconstrained strongly convex optimization. When constraints are
not present, i.e., $m=0$, conditions (3.23) and (3.36) defining $C_{1}$ reduce
to requirements of form $\|g_{k}\|\geq c_{1}\epsilon_{g}$ and
$\|g_{k}\|^{2}\geq c_{2}(\epsilon_{g}\|g_{k}\|+\epsilon_{g}^{2})$ for some
constants $c_{1}$ and $c_{2}$, respectively. That corresponds to _Case 1_ in
the analysis of [1], in which case $\epsilon_{g}$ is small as compared to
$\|g_{k}\|$ by some factor $\beta\in(0,1)$, so that the line search ensures an
improvement in the exact objective function – $f(x)$ in our notation. Similar
to the setting in this paper, [1] employs a relaxed line search which does not
fail even in the critical region; that is, when
$\|g_{k}\|\leq\beta\epsilon_{g}$. Their analysis then provides a level set
that the iterates cannot leave, which depends on the relaxation term
$\epsilon_{R}$ as well as $\epsilon_{\phi}$ (i.e. $\epsilon_{f}$ in the
unconstrained case) as in the definition of $C_{2}$ in our analysis. Since
strong convexity is assumed in [1], they can define this level set in terms of
a strong convexity parameter rather than a bound such as $w$ in Definition
3.10.

##  4 Numerical Experiments

We implemented Algorithm 1 in Python. We set $\nu=0.1$, $\tau=0.9$, and
$\beta_{k}=50$, for all $k$. The purpose of the numerical experiments is to
supplement the theoretical results, which are stated in terms of the merit
function $\phi$, by reporting the distance to the solution
$\|x_{k}-x^{\ast}\|$ as the iteration progresses. In order to gain an idea of
this behavior, it suffices to test only a few examples. We selected the
following three small-scale equality-constrained problems from the CUTEst set
[9].

problem | classification | objective | constraints
---|---|---|---
HS7 | OOR2-AN-2-1 | $\ln(1+x_{1}^{2})-x_{2}$ | $(1+x_{1}^{2})^{2}+x_{2}^{2}=4$
BT11 | OOR2-AN-4-3 | $-x_{1}x_{2}x_{3}x_{4}$ | $\begin{aligned} x_{1}^{3}+x_{2}^{2}&=1\\\ x_{1}^{2}x_{4}-x_{3}&=0\\\ x_{4}^{2}-x_{2}&=0\end{aligned}$
HS40 | OOR2-AY-5-3 | $\begin{aligned} (x_{1}-1)^{2}+(x_{1}-x_{2})^{2}+(x_{2}-x_{3})^{2}\\\ +(x_{3}-x_{4})^{4}+(x_{4}-x_{5})^{4}\end{aligned}$ | $\begin{aligned} x_{1}+x_{2}^{2}+x_{3}^{3}&=-2+\sqrt{18}\\\ x_{2}+x_{4}+x_{3}^{2}&=-2+\sqrt{8}\\\ x_{1}-x_{5}&=2\end{aligned}$

We add uniformly distributed random noise to the exact function values and to
each component of the exact gradients; i.e., for
$\xi_{i}\sim\mathcal{U}(-\epsilon_{1},\epsilon_{1})$, and
$\psi_{ij}\sim\mathcal{U}(-\epsilon_{2},\epsilon_{2})$ we set

$\displaystyle\tilde{f}(x)$
$\displaystyle=f(x)+\xi_{0},\qquad~{}~{}~{}\tilde{c}_{i}(x)=c_{i}(x)+\xi_{i}$
$\displaystyle\tilde{g}_{i}(x)$
$\displaystyle=g_{i}(x)+\psi_{0j},\qquad\tilde{J}_{ij}(x)=J_{ij}(x)+\psi_{ij}.$

In our tests, we vary $\epsilon_{1},\epsilon_{2}$, and report
$\|x_{k}-x^{\ast}\|$, where $x^{\ast}$ is a locally optimal solution obtained
by using exact gradients in the algorithm. For each of these problems,
$x^{\ast}$ is a a nondegenerate stationary point.

#### Asymptotic Behavior.

In Figure 1(c), we plot $\|x_{k}-x^{\ast}\|$ for 1000 iterations, for
$\epsilon_{1}=\epsilon_{2}=10^{-3}$ in the definitions of $\xi_{i}$, and
$\psi_{ij}$ We also display the values of
$\epsilon_{f},\epsilon_{c},\epsilon_{g},\epsilon_{J}$ defined in (3.1)-(3.2)
We should note that in each of the runs the penalty parameter $\pi_{k}$ became
fixed within the first 15 iterations. We observe that
$\\{\|x_{k}-x^{\ast}\|\\}$ is contained in a band whose upper bound is
frequently visited by the algorithm, whereas the lower bound is defined by
large irregular spikes. These results suggest that if one desires the highest
accuracy in the solution, the algorithm should continue beyond the point where
oscillations in the merit function occur, since there is little risk that the
iterates will stray away from the neighborhood of the solution, and there is a
chance that significantly higher accuracy is achieved at some iterates.

Figure 4.1: Distance to optimality ($\log_{2}(\|x_{k}-x^{\ast}\|)$) vs
iteration number for $\epsilon_{1}=\epsilon_{2}=10^{-3}$

(a) HS7. $\epsilon_{f}=10^{-3},\epsilon_{c}=10^{-3},\epsilon_{g}=1.41\times
10^{-3},\epsilon_{J}=1.41\times 10^{-3}$

(b) BT11. $\epsilon_{f}=10^{-3},\epsilon_{c}=3\times
10^{-3},\epsilon_{g}=2.24\times 10^{-3},\epsilon_{J}=6.71\times 10^{-3}$

(c) HS40. $\epsilon_{f}=10^{-3},\epsilon_{c}=3\times
10^{-3},\epsilon_{g}=2\times 10^{-3},\epsilon_{J}=6\times 10^{-3}$

#### Benefits of the relaxed line search.

The only unconventional part of Algorithm 1 is the relaxed line search (2.4).
To observe the effect of the relaxation, we solved the test problems with and
without it; the results are reported in Tables 4.1–4.3. We observe that when
the relaxation is disabled, the line search often fails in a neighborhood of
$x^{\ast}$ (we terminate the algorithm as soon as there is a line search
failure). When the relaxation is enabled, the line search is always
successful. In this case, we let the algorithm run for 100, 500, and 1000
iterations. It is apparent that the relaxed line search allows the algorithm
to continue iterating past the point where the traditional line search would
fail, yielding much better accuracy in the solution.

Table 4.1: $\min_{k}\\{\|x_{k}-x^{\ast}\|\\}$ when $\epsilon_{1}=\epsilon_{2}=10^{-5}$ | relaxation disabled | relaxation enabled
---|---|---
problem | iter. of failure | $\min_{k}\\{\|x_{k}-x^{\ast}\|\\}$ | $k_{\max}=100$ | $k_{\max}=500$ | $k_{\max}=1000$
HS7 | 77 | 7.8260E-3 | 1.0234E-3 | 4.9413E-8 | 4.9413E-8
BT11 | 64 | 4.8346E-2 | 3.9258E-3 | 1.9791E-6 | 1.4133E-6
HS40 | 26 | 3.4728E-2 | 2.1251E-3 | 1.09888E-6 | 1.0988E-6
Table 4.2: $\min_{k}\\{\|x_{k}-x^{\ast}\|\\}$ when $\epsilon_{1}=\epsilon_{2}=10^{-3}$ | relaxation disabled | relaxation enabled
---|---|---
problem | iter. of failure | $\min_{k}\\{\|x_{k}-x^{\ast}\|\\}$ | $k_{\max}=100$ | $k_{\max}=500$ | $k_{\max}=1000$
HS7 | 42 | 8.0390E-2 | 1.0401E-3 | 4.9328E-6 | 4.9328E-6
BT11 | 18 | 9.3324E-1 | 4.0003E-3 | 1.9804E-4 | 1.4060E-4
HS40 | 6 | 6.4144E-2 | 2.2293E-3 | 1.1183E-4 | 4.9328E-6
Table 4.3: $\min_{k}\\{\|x_{k}-x^{\ast}\|\\}$ when $\epsilon_{1}=\epsilon_{2}=10^{-1}$ | relaxation disabled | relaxation enabled
---|---|---
problem | iter. of failure | $\min_{k}\\{\|x_{k}-x^{\ast}\|\\}$ | $k_{\max}=100$ | $k_{\max}=500$ | $k_{\max}=1000$
HS7 | 10 | 3.7404E-1 | 1.3113E-3 | 4.5607E-4 | 2.5422E-4
BT11 | 8 | 1.7108 | 2.0598E-2 | 2.0598E-2 | 1.9451E-2
HS40 | 2 | 1.1817E-1 | 5.8202E-2 | 3.8673E-2 | 3.8673E-2

#### Effect of incorrect noise level estimations.

In Algorithm 1, estimations of $\epsilon_{f}$ and $\epsilon_{c}$ are needed to
set the relaxation bound $\epsilon_{R}$ in (2.4). It is clear that
underestimating the noise level can cause failure of the relaxed line search,
which never fails when the true level (or an overestimation) is provided. On
the other hand, overestimation can lead to large oscillations. The precise
behavior of the algorithm will depend on the stop test, and there is no
universally adopted stopping criterion in the noisy setting, to our knowledge.

Nevertheless, we performed the following experiments using a stop test that
that could be considered as a naive modification of termination tests in
standard packages. We simply terminate the algorithm when the observed (noisy)
feasibility and optimality errors are smaller than the (estimated) noise
provided for these quantities, i.e.,

$\|\tilde{c}(x_{k})\|_{1}\leq\epsilon_{c}^{est}\quad\mbox{and}\quad\|\tilde{g}(x_{k})+\tilde{J}(x_{k})^{T}\lambda_{k}\|\leq\epsilon_{g}^{est}+\|\lambda_{k}\|_{\infty}\epsilon_{J}^{est}.$
(4.1)

Figures 4.6–4.4 report the quantity $\min_{k}\\{\|x_{k}-x^{\ast}\|\\}$ when
the algorithm employs estimated noise levels $\epsilon_{1}^{est}$ and
$\epsilon_{2}^{est}$ that are 10, 100 and 1000 times larger or smaller than
the correct values. We perform this experiment for
$\epsilon_{i}=10^{-1},10^{-3},10^{-5}$. A termination due to the satisfaction
of the condition (4.1) is marked with _(opt)_ , and a line search failure is
marked with _(ls)_.

Table 4.4: $\min_{k}\\{\|x_{k}-x^{\ast}\|\\}$ when true $\epsilon_{i}=10^{-5};\ i=1,2$ | $\epsilon_{i}^{est}=\epsilon_{i}$ | $\epsilon_{i}^{est}=0.001\epsilon_{i}$ | $\epsilon_{i}^{est}=1000\epsilon_{i}$
---|---|---|---
problem | iter. | $\min_{k}\\{\|x_{k}-x^{\ast}\|\\}$ | iter. | $\min_{k}\\{\|x_{k}-x^{\ast}\|\\}$ | iter. | $\min_{k}\\{\|x_{k}-x^{\ast}\|\\}$
HS7 | 188 (opt) | 3.2017E-6 | 69 (ls) | 8.8000E-3 | 74 (opt) | 5.5704E-3
BT11 | 233 (opt) | 2.4010E-6 | 64 (ls) | 4.8346E-2 | 64 (opt) | 4.0112E-2
HS40 | 2703 (opt) | 8.0766E-7 | 26 (ls) | 3.4728E-2 | 27 (opt) | 2.8305E-2
Table 4.5: $\min_{k}\\{\|x_{k}-x^{\ast}\|\\}$ when true $\epsilon_{i}=10^{-3};\ i=1,2$ | $\epsilon_{i}^{est}=\epsilon_{i}$ | $\epsilon_{i}^{est}=0.01\epsilon_{i}$ | $\epsilon_{i}^{est}=100\epsilon_{i}$
---|---|---|---
problem | iter. | $\min_{k}\\{\|x_{k}-x^{\ast}\|\\}$ | iter. | $\min_{k}\\{\|x_{k}-x^{\ast}\|\\}$ | iter. | $\min_{k}\\{\|x_{k}-x^{\ast}\|\\}$
HS7 | 117 (opt) | 3.5750E-4 | 42 (ls) | 8.0390E-2 | 39 (opt) | 5.4414E-2
BT11 | 149 (opt) | 2.7466E-4 | 29 (ls) | 5.3925E-1 | 22 (opt) | 5.9597E-1
HS40 | 154 (opt) | 4.2653E-4 | 7 (ls) | 6.4142E-2 | 2 (opt) | 6.9002E-2
Table 4.6: $\min_{k}\\{\|x_{k}-x^{\ast}\|\\}$ when true $\epsilon_{i}=10^{-1};\ i=1,2$ | $\epsilon_{i}^{est}=\epsilon_{i}$ | $\epsilon_{i}^{est}=0.1\epsilon_{i}$ | $\epsilon_{i}^{est}=10\epsilon_{i}$
---|---|---|---
problem | iter. | $\min_{k}\\{\|x_{k}-x^{\ast}\|\\}$ | iter. | $\min_{k}\\{\|x_{k}-x^{\ast}\|\\}$ | iter. | $\min_{k}\\{\|x_{k}-x^{\ast}\|\\}$
HS7 | 51 (opt) | 2.6752E-2 | 556 (ls) | 2.5682E-4 | 5 (opt) | 4.2796E-1
BT11 | 20 (opt) | 6.6650E-1 | 3233 (ls) | 6.8738E-3 | 2 (opt) | 2.4045
HS40 | 210 (opt) | 5.82E-2 | 982 (ls) | 1.4785E-2 | 0 (opt) | 2.8877E-1

As expected, underestimations cause line search failures while overestimations
cause (4.1) to be triggered at earlier iterations. Another consequence of
underestimating $\epsilon_{2}$ is that the algorithm might never be able to
satisfy (4.1), even if a line search failure occurs sufficiently late in the
run; see for example the entry corresponding to
$\epsilon_{i}=10^{-1},\epsilon_{i}^{est}=0.1\epsilon_{i}$. In summary, over-
and underestimation of the noise levels can be harmful in ways that are
dependent on the implementation.

We must point out that an optimization algorithm may provide an indication
that the noise estimates must be re-computed. For example, the recovery
procedure described by Berahas et al. [1] uses information from the line
search to request a better estimate (e.g. through sampling or finite
difference tables), and can take precautions to avoid harmful iterations.
Robust implementations of methods for constrained optimization in the presence
of noise should include such features.

##  5 Final Remarks

Two questions guided this research. What is the best behavior one can expect
of a constrained optimization method when functions and constraints contain a
moderate amount of bounded noise that cannot be diminished at will? What are
the minimal modifications of a classical optimization algorithm that allow it
to tolerate noise, when the noise level can be estimated?

In this paper, we focused on a classical sequential quadratic programming
method applied to equality constrained problems. We showed that a modification
(relaxation) of the line search allows the iterates to approach a region
around the solution where noise dominates—and that the iterates remain in a
vicinity of this region, under normal circumstances. The analysis is presented
under benign assumptions, for example that the Jacobian of the constraints is
never close to singular, which facilitates the choice of the penalty
parameter. Nevertheless, we believe that the essence of the analysis captures
some of the main challenges to be confronted when functions and derivatives
contain noise. The accuracy bounds presented in this paper will be sharpened
in a forthcoming paper that studies the local behavior of the method near a
well behaved minimizer.

The thorny issue of how to design a proper stop test that reflects the desires
of the users has not been addressed in this paper and is worthy of research.
The treatment of singularity and the use of a nondiagonal Hessian also
requires attention, as well as the very important question of how to handle
noisy inequality constraints.

Acknowledgement. We thank Shigeng Sun for his careful reading of the paper and
useful suggestions.

## References

  * [1] Albert S Berahas, Richard H Byrd, and Jorge Nocedal.  Derivative-free optimization of noisy functions via quasi-newton methods.  SIAM Journal on Optimization, 29(2):965–993, 2019.
  * [2] Albert S Berahas, Frank E Curtis, Michael J O’Neill, and Daniel P Robinson.  A stochastic sequential quadratic optimization algorithm for nonlinear equality constrained optimization with rank-deficient jacobians.  arXiv preprint arXiv:2106.13015, 2021.
  * [3] Albert S Berahas, Frank E Curtis, Daniel Robinson, and Baoyu Zhou.  Sequential quadratic optimization for nonlinear equality constrained stochastic optimization.  SIAM Journal on Optimization, 31(2):1352–1379, 2021.
  * [4] Dimitri P Bertsekas.  Convex Optimization Algorithms.  Athena Scientific, 2015.
  * [5] R. H. Byrd, J. Nocedal, and R.A. Waltz.  KNITRO: An integrated package for nonlinear optimization.  In G. di Pillo and M. Roma, editors, Large-Scale Nonlinear Optimization, pages 35–59. Springer, 2006.
  * [6] Frank E Curtis, Daniel P Robinson, and Baoyu Zhou.  Inexact sequential quadratic optimization for minimizing a stochastic objective function subject to deterministic nonlinear equality constraints.  arXiv preprint arXiv:2107.03512, 2021.
  * [7] R. Fletcher.  Practical Methods of Optimization.  Wiley, second edition, 1987.
  * [8] Gene H. Golub and Charles F. Van Loan.  Matrix Computations.  Johns Hopkins University Press, Baltimore, second edition, 1989.
  * [9] Nicholas IM Gould, Dominique Orban, and Philippe L Toint.  CUTEst: a constrained and unconstrained testing environment with safe threads for mathematical optimization.  Computational Optimization and Applications, 60(3):545–557, 2015\.
  * [10] M Hintermüller.  Solving nonlinear programming problems with noisy function values and noisy gradients.  Journal of optimization theory and applications, 114(1):133–169, 2002.
  * [11] Jorge J Moré and Stefan M Wild.  Estimating derivatives of noisy simulations.  ACM Transactions on Mathematical Software (TOMS), 38(3):19, 2012\.
  * [12] Jorge Nocedal and Stephen Wright.  Numerical Optimization.  Springer New York, 2 edition, 1999.
  * [13] BT Poljak.  Nonlinear programming methods in the presence of noise.  Mathematical programming, 14(1):87–97, 1978.
  * [14] K Schittkowski.  Nlpqlp-nonlinear programming with non-monotone and distributed line search, 2014.
  * [15] Hao-Jun Michael Shi, Melody Qiming Xuan, Figen Oztoprak, and Jorge Nocedal.  On the numerical performance of derivative-free optimization methods based on finite-difference approximations.  arXiv preprint arXiv:2102.09762, 2021.
  * [16] Gilbert W Stewart.  On the perturbation of pseudo-inverses, projections and linear least squares problems.  SIAM review, 19(4):634–662, 1977.


# Evaluation of Embeddings of Laboratory Test Codes for Patients at a Cancer
Center

Lorenzo A. Rossi, Chad Shawber, Janet Munu  Center for Informatics
City of Hope National Medical CenterDuarte, CAUSA<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>and  Finly Zachariah  Department of Clinical
Supportive Care
City of Hope National Medical CenterDuarte, CAUSA<EMAIL_ADDRESS>
(2019)

###### Abstract.

Laboratory test results are an important and generally high dimensional
component of a patient’s Electronic Health Record (EHR). We train embedding
representations (via Word2Vec and GloVe) for LOINC codes of laboratory tests
from the EHRs of about 80,000 patients at a cancer center. To include
information about lab test outcomes, we also train embeddings on the
concatenation of a LOINC code with a symbol indicating normality or
abnormality of the result. We observe several clinically meaningful
similarities among LOINC embeddings trained over our data. For the embeddings
of the concatenation of LOINCs with abnormality codes, we evaluate the
performance for mortality prediction tasks and the ability to preserve
ordinality properties: i.e. a lab test with normal outcome should be more
similar to an abnormal one than to the a very abnormal one.

††copyright: none††conference: KDD Workshop on Applied Data Science for
Healthcare; August 2019; Anchorage, AK††journalyear: 2019

##  1\. Introduction

A patient’s Electronic Health Record (EHR) contains textual data (in form of
provider notes), numerical data (e.g. lab results and vitals) and sets of
codes for laboratory tests, diagnoses and medications. The codes for the
laboratory tests are expressed via the Logical Observation Identifiers Names
and Codes (LOINC) standard. There are 70,000 LOINC codes. Especially due to
the presence of free text and of different types clinical codes, EHR data
requires potentially very high dimensional representations of patient
information, with a challenge to design machine learning models that for many
institutions like ours can only be trained over a relatively limited number of
instances.

Word embeddings techniques such as Word2Vec (Mikolov et al., 2013) and GloVe
(Pennington et al., 2014) are unsupervised approaches to represent text in low
dimensional spaces. They are based on the principle that different words in
similar contexts may have similar meanings and therefore can be represented by
similar vectors. There is a growing interest in applying embeddings to the
healthcare domain to generate patient representations. In (Choi et al.,
2016a), the authors trained Word2Vec embeddings over LOINC, ICD-9 and NDC
codes associated to the insurance claims of 4 millions of patients and
evaluated them over a series of known relationships among medical concepts
(e.g disease A is treated by medication B). In (Choi et al., 2016b) embeddings
of the aforementioned types of codes are use to predict heart failure. To
date, the embeddings of medical concepts trained on the largest amount of data
(insurance claims and medical notes from 60 million patients) are those
presented in (Beam et al., 2018). Embeddings are trained from physician notes
of 100,000 and 250,000 general hospital patients and evaluated over predictive
tasks in (Dubois et al., 2017), (Wang et al., 2018) respectively. In contrast
with the settings of the aforementioned works, in our organization we have
access to data from fewer patients (in the order of tens of thousands), mostly
affected by cancer.

We focus on embedding representations for laboratory test data. In an analogy
with the application of embedding representations to natural language
processing, we can see each lab code as a word, a lab order as a sentence and
a visit as a document. In our data warehouse, a patient’s laboratory test
record includes a LOINC code to identify the type of lab ordered and an
abnormality code to indicate whether the lab result is normal or abnormal with
respect to a reference lab value. We train embeddings to represent LOINC codes
and the concatenations of LOINC and abnormality codes. We evaluate the lab
code embeddings qualitatively by looking at their groupings in bidimensional
t-SNE plots (Maaten and Hinton, 2008) (Ssec. 4.1). Moreover, we evaluate if
certain ordinality properties of the lab outcomes are preserved in the
embedding domain (Ssec. 4.2). We aggregate the lab code embeddings to form
feature representations of patients over certain time intervals and evaluate
those for mortality prediction (Ssec. 4.3).

To the best of our knowledge, this is the first work evaluating embeddings
trained on clinical codes and outcome information. Our qualitative and
quantitative evaluations suggest that:

  1. (1)

embeddings of codes associated to the same category of labs (e.g. respiratory
tests) have higher similarity to each other;

  2. (2)

GloVe is better than Word2Vec at preserving ordinality properties of lab
results;

  3. (3)

embeddings of the concatenation of LOINC and abnormality codes outperform bag
of words representations and embeddings of just LOINC codes in mortality
prediction tasks.

We share the Word2Vec embeddings of the LOINCs and the Python code used to
generate the t-SNE plot (Figure 3) on a GitHub
repository111github.com/elleros/DSHealth2019_loinc_embeddings. We believe that
this work may be useful to other healthcare institutions that are building
predictive models from the EHR without access to massive amounts of data.

##  2\. Cohorts

###  2.1. Embedding Training Cohort

We consider the laboratory tests ordered for a population of 79,081 inpatients
and outpatients, for a total of 802,238 visits and 8,280,820 lab orders, over
a period of 8 years between 2010 and 2017. This yields to a dictionary of 1098
LOINCs, after the removal of codes occurring less than 5 times. For every lab
test record, there is also available an “Abnormality Code” with 8 possible
values to indicate whether and how a result is abnormal with respect to a
range of reference values for that type of patient. Such ranges were
determined by the staff of the labs in our organization. Possible values of
the abnormality codes are: ‘N’, ‘A’, ‘AA’, ‘L’, ‘LL’, ‘H’, ‘HH’, ‘U’, standing
for normal, abnormal, extremely abnormal, low, very low, high, very high and
unknown, respectively. We consider the concatenation of a LOINC code with its
corresponding abnormality code as a whole word (e.g. 777-3_N would indicate a
normal platelet count). In this case, the dictionary size of the combined
LOINC and abnormality codes, after removing symbols appearing less than 5
times, is 2260 (so just about 2.1 times the number of distinct LOINC codes in
our data warehouse).

###  2.2. Predictive Evaluation Cohort

To evaluate the predictive performance of the embeddings, we consider a cohort
of about 21,000 inpatients and outpatients. For each patient’s EHR, a
prediction date was selected from her/his encounters. For deceased patients,
the prediction dates are set to achieve roughly a uniform distribution of the
temporal intervals to the date of death. This implies a slow decaying profile
of the survival curve. For instance, a cohort with a significant portion of
prediction dates within 30 days from the dates of death would not reflect a
realistic clinical situation. An alive patient was defined as not having a
recorded date of death and having at least one recorded encounter 1 year after
the chosen prediction date. The survival curves of deceased and alive patients
are shown in Figure 1. About 3% of the patients died within 90 days since
their own prediction date. The median ages of the alive and deceased patients
are 64 and 67 years, respectively. The sample includes patients affected by
solid tumors (69%) and hematology malignancies (28%). For each patient, we
consider an observation time of 30 days before the prediction time. LOINCs
associated to lab test recorded during the observation times are aggregated to
form feature representations.

Figure 1. Survival curves of deceased and alive patient populations.

##  3\. Methods

We trained Word2Vec (Mikolov et al., 2013) and GloVe (Pennington et al., 2014)
embeddings from LOINC codes grouped by visit and order. Codes within the same
order were shuffled. The shuffling is suggested in (Choi et al., 2016a) and
motivated by the fact that the ordering of LOINC codes within a given lab test
order is irrelevant, unlike the case of words within a sentence. For Word2Vec
embeddings we consider skip-gram and continuous bag of words (CBOW)
approaches. We produced set of embeddings with dimensions from 50 to 300.

Patients’ feature representations are obtained by averaging the laboratory
test embeddings associated to the labs for the visits during a certain period
of time. Of course, additional types of data such as demographics, diagnoses,
vitals would improve the prediction performance. However, our focus is the
evaluation of the lab test features alone. The test set consists in about
3,500 records of patients whose observation windows start from March
$22^{nd},2017$. The remaining 17,500 samples are used for cross-validation and
training. Classification is performed via Logistic Regression with randomnized
hyperparameter tuning. The embeddings for the prediction tasks have been
trained only with records prior the above date, so from a population of of
about 70,000 patients, to simulate making predictions and applying embeddings
to unseen data. The embeddings for qualitative evaluation and ordinality tests
are trained on the aforementioned population of 79,081 patients.

The implementations are written in Python, using open sources libraries such
as pandas, Gensim (Řehůřek and Sojka, 2010), seaborn and scikit-learn.

##  4\. Results

###  4.1. Qualitative Evaluation

For the qualitative evaluation we perform t-SNE dimensionality reduction and
look at groups of contiguous codes in a bidimensional plane. Figure 3 shows
the t-SNE plot for the 500 most frequent LOINC embeddings trained via
200-dimensional Word2Vec (Skip-gram) vectors with size of context window equal
to 5. The codes are represented in different colors according to classes of
lab tests (e.g. Urynanalysis, Pulmonary). The classes of lab tests are defined
in the ’LOINC Table’ CSV file downloaded from the LOINC.org website. We
display explicitly the 10 largest classes of labs in our dataset and group the
remaining ones into the class ’Others’. Several clinically meaningful
subclusters of embeddings can be noticed, e.g.: tests for respiratory
capacity, antibiotic susceptibility and infectious diseases. We also
highlighted the labs for Complete Blood Count panel. We highlight one of the
clusters in Figure 2, where chemical tests for arterial blood gas are next to
tests for pulmonary capacity. We provide the 200 dimensional Word2Vec
embeddings along with the code to generate the plots in the repository:
github.com/elleros/DSHealth2019_loinc_embeddings.

For comparison, we also trained embeddings of diagnoses expressed with ICD-9
codes (generally much less frequent than lab orders in a EHR) without
observing clinically meaningful similarities. This suggests that clinically
meaningful embedding representations depend on the quantity of training
samples.

Figure 2. t-SNE plot for cluster of respiratory lab test (enlarged from
cluster at the bottom of Fig. 3).

###  4.2. Ordinality Tests

We use similarity measures to check whether quantitative relationships of the
lab results are reflected also in the embedding models for the concatenation
of LOINC and abnormality codes. For example, given a lab test with normal,
abnormal and very abnormal results, we would like normal test embeddings to
have higher similarity to the abnormal test ones than to the very abnormal
ones, i.e. given a LOINC $code^{(k)}$ with possible test outcomes ‘N’, ‘A’ and
‘AA’:

$S(code^{(k)}\\_N,code^{(k)}\\_A)>S(code^{(k)}\\_N,code^{(k)}\\_AA),$

where $S(.,.)$ is a similarity measure (e.g. cosine similarity). We defined 68
binary tests similar to the one above on our LOINC-abnormality code
emebeddings. GloVe embeddings, with an error rate between 7% and 13%,
outperform Word2Vec embeddings, whose error rate is around 50%. Since this is
a binary problem, a 50% error rate means that no ordinality relationship is
preserved with Word2Vec embeddings, at least for our training data. This does
not exclude that with a larger training set, Word2Vec could have preserved the
ordinality properties.

###  4.3. Predictive Evaluation

For the quantitative evaluation, we consider the task of predicting 90 days
patient mortality. Mortality prediction is very relevant to a cancer center
population, because it enables better decisions to optimize the end of life
experience of the patients, e.g. to recommend clinical trials with the right
timing or to avoid overtreatment. We compare the results with the bag of word
(BOW) set of features and its truncated singular value decomposition
representation (which is equivalent to PCA, but more suitable for sparse
data). The resulting feature vectors are used to train and test logistic
regression classifiers. To compare performances, we use the area under the
receiver operating characteristic curve (ROC AUC).222If our goal were to have
a full assessment of the performance, we should also include the precision-
recall curve, since the labels are highly imbalanced.. Results are shown in
Figure 4, where the ROC AUCs for BOW, truncated SVD and embeddings are
compared. For the embeddings, we considered dimensions from 50 to 300 and
context windows equal to 5. Overall, embeddings perform better than the BOW
and truncated SVD baselines with an AUC ranging from 0.74 to 0.76. 300
dimensional embeddings generally outperformed lower dimensional embeddings.
GloVe and skip-gram perform slightly better than continuous bag of words. The
ordinality preservation property of the GloVe embeddings (Ssec 4.2) does not
seem to lead to a superior prediction performance. We also varied the duration
of the observation window (e.g. to 60 days) and tried embedding aggregations
such as median, min and/or max, but with inferior results.

Besides, we evaluated embeddings of LOINC codes without bnormality symbols,
obtaining AUC scores below 0.65. To the best of our knowledge, the only
publicly available set of embeddings encompassing LOINC codes is the one
presented in (Choi et al., 2016a).333github.com/clinicalml/embeddings The
embeddings are trained over sets of LOINC, ICD-9 and NDC codes from 4 million
insured subjects. The LOINCs in this embedding model are about 3,500, with 493
of them in common with our dataset. We compared the prediction performance
restricted to the aforementioned 493 embeddings. Overall the AUC score from
the public embeddings trained over insurance claims is slightly below the
performance of the embeddings trained over our patient records.

Figure 3. t-SNE plot of 500 LOINC embeddings trained via Word2Vec. The
dimension of the embeddings is 200. The colors represent classes of laboratory
tests defined on the LOINC.org website. Figure 4. Comparison of ROC AUC for
features derived from embeddings in a 90 day mortality prediction task. The
dimension of the representations is on the horizontal axis. Truncated Singular
Value Decomposition (SVD) is a dimensionality reduction technique applied to
the bag of words representation. Continuous Bag of Words (CBOW) is an approach
to train Word2Vec embeddings.

##  5\. Conclusion

We have trained embeddings for laboratory test data, in form of LOINC and
concatenations of LOINC and abnormality codes, for a population of 79,081
cancer patients. The embeddings have been trained via the skip-gram, CBOW and
GloVe approaches for different dimensions and sizes of the context window. The
embeddings of LOINC + abnormality codes yield to better prediction results
than the embeddings of just LOINC codes providing a potentially useful tool
for low dimensional representations of lab test outcome. Generally, embedding
based features outperform in prediction tasks features from baselines such as
bag of words and truncated SVD. Moreover, GloVe embeddings seem to preserve
ordinality properties of the lab test outcomes such as the pair (abnormal,
normal) having higher similarity than (very abnormal, normal).

## References

  * (1)
  * Beam et al. (2018) Andrew L Beam, Benjamin Kompa, Inbar Fried, Nathan P Palmer, Xu Shi, Tianxi Cai, and Isaac S Kohane. 2018.  Clinical Concept Embeddings Learned from Massive Sources of Medical Data. (2018).  arXiv preprint arXiv:1804.01486.
  * Choi et al. (2016b) Edward Choi, Andy Schuetz, Walter F. Stewart, and Jimeng Sun. 2016b.  Medical Concept Representation Learning from Electronic Health Records and its Application on Heart Failure Prediction.  CoRR abs/1602.03686 (2016).
  * Choi et al. (2016a) Youngduck Choi, Chill Yi-I Chiu, and David Sontag. 2016a.  Learning low-dimensional representations of medical concepts.  AMIA Summits on Translational Science Proceedings 2016 (2016), 41\.
  * Dubois et al. (2017) Sebastien Dubois, David C Kale, Nigam Shah, and Kenneth Jung. 2017\.  Learning Effective Representations from Clinical Notes. (2017).  arXiv preprint arXiv:1705.07025.
  * Maaten and Hinton (2008) Laurens van der Maaten and Geoffrey Hinton. 2008.  Visualizing data using t-SNE.  Journal of machine learning research 9, Nov (2008), 2579–2605.
  * Mikolov et al. (2013) Tomas Mikolov, Ilya Sutskever, Kai Chen, Gregory S. Corrado, and Jeffrey Dean. 2013.  Distributed Representations of Words and Phrases and their Compositionality. In NIPS. 3111–3119.
  * Pennington et al. (2014) Jeffrey Pennington, Richard Socher, and Christopher D. Manning. 2014.  Glove: Global Vectors for Word Representation. In EMNLP. ACLWEB, 1532–1543.
  * Řehůřek and Sojka (2010) Radim Řehůřek and Petr Sojka. 2010.  Software Framework for Topic Modelling with Large Corpora. In Proceedings of the LREC 2010 Workshop on New Challenges for NLP Frameworks. ELRA, Valletta, Malta, 45–50.  http://is.muni.cz/publication/884893/en.
  * Wang et al. (2018) Yanshan Wang, Sijia Liu, Naveed Afzal, Majid Rastegar-Mojarad, Liwei Wang, Feichen Shen, Paul Kingsbury, and Hongfang Liu. 2018\.  A comparison of word embeddings for the biomedical natural language processing.  Journal of biomedical informatics 87 (2018), 12–20.